# The Effects of Starspots on Spectroscopic Mass Estimates of Low-mass Young
Stars

C. Flores M. S. Connelley B. Reipurth Institute for Astronomy, University of
Hawaii at Manoa, 640 N. Aohoku Place, Hilo, HI 96720, USA
G. Duchêne Astronomy Department, University of California, Berkeley, CA 94720,
USA  Université Grenoble Alpes, CNRS, Institut de Planétologie et
d’Astrophysique, IPAG, 38000 Grenoble, France

(Accepted Nov 5, 2021)

###### Abstract

Magnetic fields and mass accretion processes create dark and bright spots on
the surface of young stars. These spots manifest as surface thermal
inhomogeneities, which alter the global temperature measured on the stars. To
understand the effects and implications of these starspots, we conducted a
large iSHELL high-resolution infrared spectroscopic survey of T Tauri stars in
Taurus-Auriga and Ophiuchus star-forming regions. From the K band spectra, we
measured stellar temperatures and magnetic field strengths using a magnetic
radiative transfer code. We compared our infrared-derived parameters against
literature optical temperatures and found a) a systematic temperature
difference between optical and infrared observations, and b) a positive
correlation between the magnetic field strengths and the temperature
differences. The discrepant temperature measurements imply significant
differences in the inferred stellar masses from stellar evolutionary models.
To discern which temperature better predicts the mass of the star, we compared
our model-derived masses against dynamical masses measured from ALMA and PdBI
for a sub-sample of our sources. From this comparison we conclude that, in the
range of stellar masses from 0.3 to 1.3 $\rm M_{\odot}$, neither infrared nor
optical temperatures perfectly reproduce the stellar dynamical masses. But, on
average, infrared temperatures produce more precise and accurate stellar
masses than optical ones.

Starspots — Stellar magnetic fields — Pre-main sequence stars — High
resolution spectroscopy — Radiative transfer

††journal: ApJ††facilities: IRTF, ALMA, PdBI††software: astropy (Astropy
Collaboration et al., 2013) MoogStokes (Deen, 2013), emcee (Foreman-Mackey et
al., 2013), Spextools (Cushing et al., 2004), xtellcor (Vacca et al., 2003)

##  1 Introduction

Young stars have magnetic fields that directly and indirectly produce hot and
cold spots on their stellar surfaces. Hot spots produced through accretion
appear when circumstellar material channeled from the disk impacts the stellar
surface at almost free-fall speed (Gullbring et al., 2000; Hartmann et al.,
2016). Hot spots can also be produced in the form of stellar plages and
faculae due to an increase in surface magnetic activity, similar to what is
observed in the Sun (Strassmeier, 2009). Meanwhile cold stellar spots appear
when magnetic fields inhibit the transport of energy through convection,
producing localized dark regions (Berdyugina, 2005).

Although there has been substantial research on the photometric effects of
spots in stars, spectroscopic signatures have often been neglected, in part
because spectroscopic monitoring of stars is observationally expensive but
also because most evolved stars host weak magnetic fields. Indeed, in a study
of low-mass stars with ages from less than one Myr to over a Gyr old, Vidotto
et al. (2014) found that the average magnetic field strength derived from
polarimetry in young stars is thousands of times stronger than the magnetic
field value of main-sequence stars. This field difference means that the size
and/or temperature contrast of cold spots on young stars, could be
considerably larger than what is observed on their more evolved counterparts.
Hints of this phenomenon have emerged from individual sources, where
significant differences were reported for temperature measurements performed
at different wavelengths ($\rm\Delta\,T>300$ K) (e.g., Vacca & Sandell, 2011;
Gully-Santiago et al., 2017; Flores et al., 2019, 2020).

One of the main problems of these discrepant temperature measurements is that
the effective temperature of such stars can no longer be obtained from a
single small-range spectroscopic observation (Gully-Santiago et al., 2017).
For low-mass young stars, this additionally implies that their calculated
masses could be systematically offset, as their effective temperatures
strongly correlates with their stellar masses. This occurs because low-mass
stars contract almost isothermally during the first 5-20 Myrs of their
evolution as they descend the Hayashi track. Therefore, understanding and
correcting for this effect is of paramount importance as direct mass
measurements of young stars are often hard and expensive to obtain and, most
of the time, simply impossible (Guilloteau et al., 2014; Simon et al., 2017,
2019; Braun et al., 2021).

In this paper, we present iSHELL K band observations of 40 bright ($\rm K<11$)
young stellar sources from Ophiuchus and Taurus-Auriga star forming regions.
We mostly selected sources from Simon et al. (2019) and the ODISEA project
(Cieza et al., 2019) and did our best to try to avoid spectroscopic binaries.
The literature optical spectral types for the sources in our sample range from
K2 to M4.5. Optical spectral types for the sources in Taurus-Auriga were
collected from White & Hillenbrand (2004); Hartigan & Kenyon (2003); Luhman et
al. (2010); Herczeg & Hillenbrand (2014) while optical measurements for the
Ophiuchus sources were obtained from Bouvier & Appenzeller (1992); Wilking et
al. (2005); Torres et al. (2006); Erickson et al. (2011). Two stars do not
have optical spectral types (GSS 26 and GSS 39), and another source (WSB 82)
has only a broad classification of mid-K-type (Esplin et al., 2018). All but
one of the young stars in our sample are Class II sources, the exception is
Haro 6-13, which is a Class I source.

Our spectroscopic observations with iSHELL are presented in Section 2. We
derive stellar parameters for the stars using the radiative transfer code
MoogStokes (Deen, 2013) in Section 3. In Section 4, we combine optical and
infrared temperatures with the magnetic stellar evolutionary models of Feiden
(2016) to derive stellar masses. The model-derived masses are then compared to
independently measured stellar masses from interferometric observations. We
discuss our results in Section 5 and summarize our findings in Section 6.

##  2 Observations and Data Reduction

We observed 40 young stars between 2017 October 13 and 2020 October 23 using
the high-resolution near-infrared echelle spectrograph iSHELL on IRTF (Rayner
et al., 2016). The observations were performed in the $K2$ mode, thus from
2.09 to 2.38 $\mu$m; using the 0$\farcs$75 slit to achieve a spectral
resolution of R $\sim$ $50,000$. Due to a known fringing issue in the iSHELL
spectra, we decided to adopt the following observing strategy. First, we
guided on a young star and integrated on it long enough to obtain a spectrum
with $S/N>50$. We stopped guiding and immediately acquired a set of
calibration spectra consisting of $K$ band flats and arc-lamp spectra. Then,
we observed a nearby telluric A0 standard star within 0.1 airmasses of the
young source and integrated on it long enough to achieve a $S/N>100$ spectrum.
Afterward, we stopped the guider and obtained a second set of $K$ band
calibration files to correct the telluric standard star. The same process was
repeated for each star in the sample.

The SNR obtained for the young sources ranges between 56 and 175 with a median
value of 105, as reported by Spextool v5.0.2 (Cushing et al., 2004). The data
were reduced in the same manner as reported in Flores et al. (2020).

Figure 1: Left panel: temperature difference between optical and K band
measurements as a function of the optical temperature. Right panel:
temperature difference between optical and K band measurements as a function
the magnetic field strengths. Both relations are significantly correlated and
show a positive trend. Data points with thick black contours correspond to the
sample of 24 stars with measured dynamical masses presented in section 4.3.

##  3 Stellar Parameter Measurements

We derive stellar parameters for the 40 young stars using a synthetic stellar
spectrum technique that combines MARCS stellar atmospheric models (Gustafsson
et al., 2008), the magnetic radiative transfer code MoogStokes (Deen, 2013),
along with atomic and molecular lines from the VALD3 (Ryabchikova et al.,
2015; Pakhomov et al., 2019) and HITEMP databases (Rothman & Gordon, 2010).

More details of the modeling technique can be found in Flores et al. (2019,
2020) but in brief, we took a two-step approach to fit the observations.
First, we varied seven stellar parameters: temperature, gravity, magnetic
field strength, projected rotational velocity, K band veiling, micro-
turbulence, and CO abundance. We performed the fitting using eight wavelength
regions, one of which contains the CO overtone starting at 2.294 $\micron$
(see Figures in Appendix A). In the second step, we excluded the CO region to
avoid potential CO contamination from the circumstellar material (as seen, for
example, in Doppmann et al., 2005; Flores et al., 2020). We anchored the
projected rotational velocity to the best value found in the previous step and
re-calculated the temperature, gravity, magnetic field strength, K band
veiling, and micro-turbulence. In both steps, we fitted the spectra with a
Markov Chain Monte Carlo (MCMC) method using the emcee package (Foreman-Mackey
et al., 2013). We ran a total of 50 individual chain “walkers” advancing a
total of 500 steps. Examples of how sensitive the stellar spectra are to
parameters such as surface gravity and magnetic field strengths are depicted
in Appendix B and C.

##  4 Results

In Table 1, we report all of the derived stellar parameters for the 40 stars.
The measured K band temperatures range from about 3000 K to 4400 K, with a
median K band temperature value of 3600 K, magnetic field strengths range from
0.7 to 3.2 kG with a median value of 1.8 kG, and surface gravity values span
the range of 3.1 to 4.3 in log of $\rm cm\,s^{-2}$, with a median value of 3.8
in log of $\rm cm\,s^{-2}$. The parameter values are reported as the median of
the MCMC distributions while the uncertainty in the stellar parameters are
reported as 3$\sigma$ values (see example in Appendix D).

###  4.1 Optical and IR temperature differences

All the literature optical spectral types were transformed into temperatures
using the temperature scales of Herczeg & Hillenbrand (2014). Our measured K
band temperatures are different, and almost always lower, than the optical
temperatures. On the left panel of Figure 1, we see a significant temperature
difference between optical and IR measurements, which increases towards
earlier optical spectral types. From the 37 stars with determined optical
spectral types, only ROX 25 has an optical temperature measurement that is
significantly lower than our K band derived temperature ($\rm\Delta
T=-285\,K$). This intriguing case could be caused by a hidden multiplicity in
the system, which might also explain the unusually weak magnetic field
strength measured for this source ($B\sim 0.7$ kG). Follow-up studies would be
necessary to confirm this hypothesis. The right panel of Figure 1 shows a
broad yet increasing stellar temperature difference as a function of the
surface magnetic field strength. Put differently, stars with stronger magnetic
fields have, on average, a larger temperature contrast between both
measurements.

We assessed the statistical significance of the trends mentioned above using
the Pearson correlation coefficient. We obtained Pearson coefficients of 0.60
and 0.59 for the temperature difference as a function of the optical
temperature and as a function of the magnetic field strength, and P-values for
testing non-correlations of $8.3\times 10^{-5}$, and $1.3\times 10^{-4}$,
respectively. These tests shows that the variables are positively correlated
and both relations are statistically significant.

Additionally, we performed a linear fit to both relations. For the temperature
difference as a function of optical temperature, we obtained $\Delta
T=(0.36\pm 0.08)\,T_{opt}-1170\pm 297~{}K$. This means that for stars hotter
than 3200 K in the optical (M4 spectral type), there is, on average, a 36 K
increase in $\Delta T$ for every 100 K increase in optical temperature. For
the temperature difference vs. magnetic field strength in kilogauss, we found
$\Delta T=(206\pm 44)\,B-135\pm 80$ K. We note that there is a significant
scatter in both correlations. Thus, the linear relationships should be
considered as an ensemble correlation rather than an object-by-object
relationship.

Furthermore, these relationships are only valid for the temperatures
considered in our sample and are expected to break at hotter and cooler
temperatures. In Section 5, we discuss physical mechanisms that could produce
these effects in young stars and further speculate on the validity of these
first-order fits outside the temperature range investigated here.

###  4.2 Stellar masses from evolutionary models

Figure 2: Left panel: temperatures and gravities plotted in the HR diagram.
Blue circles represent optical temperature measurement, while the orange ones
correspond to K band derived temperatures. Thick black contours demarcate
stars with measured dynamical masses. Representative 3$\sigma$ uncertainties
for the K band temperature and gravity of 50 K and of 0.1 dex are shown on the
upper right. Right panel: Feiden’s magnetic mass tracks (red) and isochrones
(black) are shown in a $\rm T_{eff}$ vs. gravity space. The scale on both
panels is the same.

The effective temperature of a star is a mathematical construction that
satisfies $T_{eff}=(L_{\star}/4\pi\sigma R_{\star}^{2})^{1/4}$. The
computation of effective temperatures in this form is difficult because either
radii measurements or stellar luminosities over broad spectral ranges can be
hard to obtain. Thus, it is customary to instead measure the temperature of a
star from stellar spectra using a combination of atmospheric models and
radiative transfer codes (often called spectroscopic temperature $T_{spec}$).
The spectroscopic temperature is assumed to be a good representation of the
effective temperature of a star regardless of the observed wavelength.
Although there are technical differences for temperatures measured at
different wavelengths (such as the knowledge of the atomic and molecular line
lists, the characterizations of the instrument spectral profiles and
responses, among others) proper calibrations of the stellar atmospheric models
should provide consistent values.

Flores et al. (2019) demonstrated that iSHELL $K$ band spectroscopic
observations combined with the same radiative transfer code used here, provide
accurate effective temperatures for main and post-main sequence stars (see
Figure 6 of their paper). We thus argue that the observed temperature
differences displayed in Figure 1 cannot be solely accounted for by
uncertainties in the techniques used at the two different wavelengths, instead
we propose that they reveal a true physical peculiarity of the stellar
surfaces in low-mass young stars.

Furthermore, since young stars evolve almost isothermally during their first
few million years of evolution (Hayashi, 1961), the mass of a star can be
almost directly inferred from its effective temperature using ‘theoretical’ HR
diagram ($\log{g}$ vs. $T_{eff}$, see Figure 2). The problem of having two
different temperatures for the same star, then translate into having also two
discrepant stellar masses.

To measure how different these two masses could be, we computed stellar masses
using our gravities, the literature optical temperatures, the infrared
temperatures, and the magnetic stellar evolutionary models from Feiden (2016).
We chose these magnetic models because they provide more accurate stellar
masses than other evolutionary tracks for young sources in the mass range of
$\rm M_{\star}\lesssim 1.3\,M_{\odot}$ (Simon et al., 2019; Braun et al.,
2021). We use a ‘theoretical’ HR diagram ($\rm log{(g)}$ vs. $\rm T_{eff}$)
instead of the ‘observational’ HR diagram ($\rm L_{\star}$ vs. $\rm T_{eff}$)
because a) gravities are directly probed in the stellar spectra and b) because
it circumvents multiple sources of biases introduced by the measured
luminosity such as extinction, variability, and multiplicity.

Figure 3: Difference between masses calculated from Feiden (2016) magnetic
models using literature optical and our K band infrared temperatures. The
gravities were assumed to be the same for both measurements. Masses were only
calculated for the 35 sources that completely lie within the grid of stellar
evolutionary models.

In Figure 2, we plot our sample of stars with optical literature measurements
and infrared temperatures using our derived gravities in both cases. Since
optical temperatures are generally hotter than the IR-derived ones, stars with
optical temperatures lie on higher mass tracks than stars with K band
temperatures. From the figure, we see that only 35 of the stars in our sample
completely lie within the grid of stellar evolutionary models, and of those
stars, on average, masses obtained from optical temperatures are $28\pm 2\,\%$
higher than the stellar masses calculated from infrared temperatures.

Trends in the difference between stellar masses derived from optical
temperatures vs. infrared temperatures are depicted in Figure 3. Not
surprisingly, it shows that stellar masses derived from optical temperatures
are almost always higher than the ones obtained from infrared measurements and
this difference increases for stars with higher optically-derived masses.

###  4.3 Dynamical vs. spectroscopic masses

The natural question is then, does the optical or infrared temperature better
correspond with the mass of the star? To address this, we focus on a sub-
sample of our young stars that have dynamical masses derived from CO and CN
observations obtained with the Atacama Large Millimeter/submillimeter Array
(ALMA) and the Plateau de Bure Interferometer (PdBI) (Schaefer et al., 2009;
Guilloteau et al., 2014; Simon et al., 2019). In the following analysis, we
only consider sources that have dynamical mass measurements with a precision
better than $20\%$. After this selection criterion, 24 sources were left
within the mass range of 0.3 $\rm M_{\odot}$ to 1.3 $\rm M_{\odot}$.

Figure 4 shows the difference between the spectroscopic derived masses and the
dynamical masses calculated from interferometric observations as a percentage
of the dynamical masses. K band and optically derived masses reproduce
reasonably well the dynamical masses of sources with $\rm M_{dyn}\gtrsim
0.5\,M_{\odot}$. In fact, in this mass range, the mean percentage difference
for the IR-derived masses is $-8\pm 4\,\%$ with a standard deviation of
$18\,\%$. The optical measurements in the same range, produce a mean
percentage difference of $12\pm 6\,\%$, with a standard deviation of $24\,\%$.
For lower stellar masses $\rm M_{dyn}<0.5\,M_{\odot}$, both spectroscopic
temperatures produce masses in excess of the measured dynamical masses, with
the optically-derived masses performing significantly worse than the infrared
ones. The infrared derived masses, over-predict the stellar masses by $31\pm
8\,\%$, while the optical masses over-predict the stellar masses by $94\pm
16\,\%$.

Considering the full mass range ($\rm 0.3-1.3\,M_{\odot}$), the infrared-
derived masses have a mean percentage difference of $4\pm 5\,\%$ with a
standard deviation of $26\,\%$. The optical masses, in turn, have a mean
percentage difference of $36\pm 8\,\%$ with a standard deviation of $48\,\%$.
Although neither model-derived mass is a perfect representation of the
dynamical mass of the stars, this analysis implies that masses derived from K
band observations are not only more precise than masses derived from optical,
but they are also more accurate.

Figure 4: Both panels show the difference between model-derived stellar masses
and dynamical masses as a percentage of the dynamical masses for 24 young
stars. The upper panel shows the result for masses derived using our K band
temperatures. The bottom panel displays the results for masses calculated
using optical literature temperatures. The K band derived masses are not only
more precise but also more accurate than the optical masses.

##  5 Discussion

###  5.1 Starspots on young stars

Our interpretation of the temperature differences between optical and infrared
observations is that they are caused by thermal inhomogeneities on the surface
of the young stars, i.e., cold and hot spots. The presence of spots is
supported by the fact that all the observed young stars have strong ($\sim$kG)
magnetic fields, and the strength of the magnetic fields positively correlates
with the observed temperature difference. Although this alone cannot prove a
causation between magnetic fields and starspots, we note that temperature
inhomogeneities in main and post-main sequence stars are intrinsically linked
to magnetic effects. This is the case, for example, for the Sun, the well
studied RS CVn Stars (see reviews of Berdyugina, 2005; Strassmeier, 2009), the
BY Drac type stars (Bopp & Evans, 1973; Anderson et al., 1976), and the FK
Comae stars (Bopp & Stencel, 1981; Jetsu et al., 1993).

As shown in Figure 5, our measured magnetic field strengths positively
correlate with literature optical temperatures. This finding, however, is in
contradiction with results found in the previous study of Sokal et al. (2020)
who compiled magnetic field strengths of classical T Tauri stars (cTTS) and
weak-line T Tauri stars from different authors (e.g., Johns-Krull, 2007; Yang
& Johns-Krull, 2011; Lavail et al., 2017), and plotted them against their
optical temperatures. With a sample of 28 cTTS, 12 wTTS, and one Class I
source, they found an anti-correlation between the measured magnetic field
strengths and the optically measured temperatures, which was attributed to a
possible evolutionary change in the stars.

Figure 5: Magnetic field strengths as a function of literature optical
temperatures. The increasing magnetic field strengths as a function of optical
temperature are opposite to a trend found by Sokal et al. (2020).

To check for biases in the measured magnetic field strengths, we compared the
nine stars in common between the compilation of Sokal et al. (2020) and our
survey. These nine sources have magnetic field strengths ranging from 0.7 to
2.4 kG, and differ by $\sim 13\%$ (median difference) with our B-field
strengths, with the higher values presented in Sokal et al. (2020).

We propose two explanations for the observed differences; a) a real variations
in the B-field strength of the stars over decade-long periods. The
abovementioned sources were observed by Johns-Krull (2007) between 1996 and
2004, leading to a 13-24 years difference with our observations (see Section
2). Although no study has proven magnetic variability over such a long
timespan using a consistent technique, Lavail et al. (2019) measured magnetic
field strength variability of a sample of seven young stars over a period of
2-5 months. They found a typical variation of the order of 0.3 kG (or $\sim$
10%), which could suggest that magnetic field variability is possible over
longer timespans. b) another possible explanation for the discrepant magnetic
field strengths is the difference of the instruments, atmospheric lines,
radiative transfer codes, and underlying assumptions used in each study. The
abovementioned nine sources were observed with CSHELL which allowed to
observed three wavelength regions in the $K$ band: two Ti-line regions, and
one region that targets part of the CO overtone (Johns-Krull, 2007). The
wavelength range of the CSHELL observation was more than ten times smaller
than what we acquired with iSHELL and also at half the spectral resolution.
Furthermore, stellar parameters such as temperature, gravity, and rotational
velocity were not self-consistently derived but adopted from previous studies.

Since the two modeling techniques and data quality are significantly different
(due to the rapid technological change in detectors and computers), we lean
towards the latter explanation rather than the former. However, we do not
entirely discard real magnetic field strength variations as a possibility. In
any case, given that our results are opposite of the trend found in Sokal et
al. (2020) and because both studies contain a similar number of sources, only
a follow-up survey with a larger number of stars can confirm or support either
trend.

###  5.2 Temperature difference and stellar structure

The positive correlation between the optical temperature and the temperature
difference (shown in Figure 1) is not expected to extend indefinitely in
either temperature direction. Towards the cool end, at about an optical
temperature of $\sim 3200$ K, stars should manifest a negligible optical to
infrared temperature difference. One way to produce this is that the
temperature of the cool spots on the stellar surface approaches the
temperature of the warm photosphere. For example, they might follow a fixed
$T_{\rm spot}/T_{\rm phot}$ ratio, similar to what is seen in RS CVn stars
(e.g., Catalano et al., 2002; Frasca et al., 2008), which would reduce the
overall contrast between optical and IR observations. This alone, however,
would not completely explain the observed pattern towards the cool end, as it
would also be necessary that hot spots become less important.

For stars on the hot end, the temperature differences could keep increasing as
a function of their effective temperature. These differences, however, should
eventually decline as the convective envelope of the stars become shallower
and the surface magnetic fields of the stars weaken. Lavail et al. (2017)
measured the magnetic field strengths of a sample of six intermediate mass T
Tauri stars (IMTTS) and found evidence that their values were, on average,
lower than what has been measured for low-mass T Tauri stars.

The internal structure and thus magnetic field of the stars is expected to
strongly depend on the stellar temperature and age (see e.g., Gregory et al.,
2012; Villebrun et al., 2019). This means that observing earlier type young
stars ($T_{opt}>5000$ K) in both the infrared and optical could provide direct
observational evidence of convective layers in stars becoming shallower, thus
allowing us to understand the evolution of stellar internal structure better
and posing constraints on evolutionary models.

###  5.3 Evolutionary models and dynamical masses

In Section 4.3, we showed that the infrared-derived masses are more accurate
and precise than the ones derived from optical temperatures. However, the
underlying assumption in our analysis is that the magnetic models reported by
Feiden (2016) are correct. It could be argued that this analysis shows that
the magnetic models from Feiden provide incorrect masses for a given
temperature. Although we cannot prove otherwise, we can rely on previous
studies that have shown the reliability of Feiden (2016) magnetic models.
Simon et al. (2019) compared magnetic vs. non-magnetic models from different
authors (Baraffe et al., 2015; Feiden, 2016) and determined that magnetic
models produce better results for young stars in a similar mass range to what
is investigated here. Recently, Braun et al. (2021) also tested various
stellar evolutionary models against ALMA dynamical masses of young stars and
concluded that the Feiden (2016) magnetic models provide the best estimate of
the stellar masses in the mass range of 0.6 to 1.3 $\rm M_{\odot}$. For
sources less massive than 0.6 $\rm M_{\odot}$, Braun et al. (2021) found that
the Feiden (2016) models overestimate the stellar masses. Our study confirms
these findings but also suggests that this problem is more severe when using
optical temperatures (as was done in Braun et al., 2021), than when using
infrared-derived temperatures of the young stars.

There are two other aspects that must be kept in mind regarding our analysis.
First, we assumed that our measured gravities can be used in combination with
literature optical data to derive stellar masses. This assumption makes sense
from a physical perspective, as starspots should be in hydrostatic equilibrium
with the ambient photosphere, but in practice, the gravities and optical
temperatures were calculated independently from different datasets,
instruments, and techniques. Second, the dynamical mass values adopted in our
computations could be more uncertain than what was originally reported by
Simon et al. (2019). Braun et al. (2021) used archival ALMA data to measure
dynamical masses of young stars in Taurus. Their sample overlapped for 8
sources with the Taurus sources analyzed by Simon et al. (2019), and
discrepancies of dynamical mass measurements were as high as $46\%$. As
explained in Braun et al. (2021), these stellar-mass differences mainly arise
due to the difference in the adopted disk inclination angles (since stellar
masses scale as $M_{\star}\propto 1/\sin^{2}{i}$).

Future gas observations of these disks at higher resolution and $S/N$ could
help to decide which measurements we should favor.

##  6 summary and conclusion

We obtained iSHELL K band observations of 40 pre-main-sequence sources from
Taurus-Auriga and Ophiuchus star-forming regions. We used a magnetic radiative
transfer code to derive parameters of the stars such as temperature, gravity,
and magnetic field strength. We compared our derived values against optical
literature temperatures and used stellar evolutionary models to test whether
optical or infrared temperatures should be used to derive stellar masses of
young stars. A summary of our findings are as follows:

  * •

Our K band derived temperatures are almost always lower than optical
temperatures from the literature, and this difference increases as a function
of the optical temperatures. Additionally, the observed difference between
optical and infrared temperatures increases as a function of the magnetic
field strengths of the stars.

  * •

We interpret that the temperature differences between optical and infrared
observations are caused by hot and cold spots on the stellar surface. This
interpretation is supported by the fact that, on average, stars with stronger
magnetic fields have larger temperature differences and is also reinforced by
the fact that magnetic activity and thermal inhomogeneities are intrinsically
linked in a variety of main and post-main sequence stars.

  * •

Using the Feiden (2016) magnetic models, we computed the stellar masses for 35
young sources using both optical and infrared temperatures. The mean
spectroscopic mass calculated from optical temperatures was $28\pm 2\%$ higher
than from the K band temperatures. This implies that the choice of optical vs.
infrared temperatures when using stellar evolutionary models is important.

  * •

We used a sub-sample of 24 stars with masses calculated from ALMA and PdBI to
explore whether optical or infrared temperatures better reproduce the
dynamical stellar masses. In the mass range of 0.3 $\rm M_{\odot}$ to 1.3 $\rm
M_{\odot}$, we found that the IR temperatures produce more precise and
accurate stellar masses than optical temperatures.

In conclusion, we have shown that magnetically-induced starspots can
significantly alter the temperature measured on young stars. Although it is
common to use optical temperatures to infer stellar masses, we have
demonstrated that infrared temperatures should be preferred in the mass range
of 0.3 to 1.3 $M_{\odot}$. Future studies on a broader range of stellar masses
will be necessary to test whether our findings also apply to higher mass young
stars.

We thank the referee for an insightful review which helped to improve this
paper. CF and MSC acknowledge support from the NASA Infrared Telescope
Facility, which is operated by the University of Hawaii under contract
80HQTR19D0030 with the National Aeronautics and Space Administration. GD
acknowledges support from NASA grants NNX15AC89G and NNX15AD95G/NExSS as well
as 80NSSC18K0442. The technical support and advanced computing resources from
University of Hawaii Information Technology Services – Cyberinfrastructure,
funded in part by the National Science Foundation MRI award $\\#1920304$, are
gratefully acknowledged. We are particularly grateful to the iSHELL team for
all the help provided during the observation and data reduction process. This
research has made use of NASA’s Astrophysics Data System Bibliographic
Services. This research has made use of the SIMBAD database, operated at CDS,
Strasbourg, France. This work has made use of the VALD database, operated at
Uppsala University, the Institute of Astronomy RAS in Moscow, and the
University of Vienna.

Table 1: Stellar parameters of young stars derived from K band.

Source | $\rm T_{\rm K-band}$ | log(g) | $\rm r_{K}$ | $v_{micro}$ | $\langle\rm B\rangle$ | $v\sin(i)$
---|---|---|---|---|---|---
name | (K) |  |  | (km s-1) | (kG) | (km s-1)
AA Tau | $3593^{+61}_{-76}$ | $3.82^{+0.08}_{-0.11}$ | $0.99^{+0.05}_{-0.05}$ | $0.28^{+0.57}_{-0.18}$ | $1.88^{+0.18}_{-0.14}$ | $14.44^{+0.85}_{-0.68}$
BP Tau | $3703^{+48}_{-35}$ | $4.26^{+0.07}_{-0.07}$ | $1.36^{+0.07}_{-0.06}$ | $2.03^{+0.22}_{-0.03}$ | $2.38^{+0.15}_{-0.14}$ | $9.94^{+2.61}_{-0.99}$
CI Tau | $3888^{+57}_{-73}$ | $3.74^{+0.10}_{-0.14}$ | $1.87^{+0.08}_{-0.08}$ | $1.19^{+0.69}_{-0.69}$ | $1.78^{+0.11}_{-0.12}$ | $12.71^{+0.40}_{-0.48}$
CX Tau | $3498^{+47}_{-53}$ | $3.64^{+0.10}_{-0.12}$ | $0.26^{+0.03}_{-0.02}$ | $2.32^{+0.43}_{-0.44}$ | $1.10^{+0.16}_{-0.16}$ | $21.41^{+0.32}_{-0.34}$
CY Tau | $3403^{+43}_{-21}$ | $3.82^{+0.06}_{-0.06}$ | $0.58^{+0.04}_{-0.04}$ | $1.84^{+0.56}_{-0.44}$ | $1.39^{+0.08}_{-0.09}$ | $10.55^{+0.54}_{-0.33}$
DE Tau | $3459^{+40}_{-52}$ | $3.67^{+0.07}_{-0.08}$ | $1.15^{+0.04}_{-0.04}$ | $2.49^{+0.34}_{-0.33}$ | $0.75^{+0.10}_{-0.10}$ | $10.33^{+0.33}_{-0.40}$
DK Tau A | $3501^{+43}_{-29}$ | $3.71^{+0.08}_{-0.09}$ | $1.90^{+0.08}_{-0.08}$ | $0.98^{+0.53}_{-0.56}$ | $2.33^{+0.16}_{-0.16}$ | $17.40^{+0.88}_{-0.72}$
DK Tau B | $3411^{+67}_{-18}$ | $3.75^{+0.06}_{-0.06}$ | $1.10^{+0.06}_{-0.04}$ | $1.58^{+0.37}_{-0.40}$ | $0.81^{+0.12}_{-0.12}$ | $9.47^{+0.34}_{-0.68}$
DL Tau | $3915^{+39}_{-41}$ | $3.81^{+0.07}_{-0.07}$ | $2.97^{+0.07}_{-0.08}$ | $0.40^{+0.59}_{-0.30}$ | $2.05^{+0.07}_{-0.09}$ | $9.77^{+0.66}_{-0.37}$
DM Tau | $3302^{+23}_{-14}$ | $3.90^{+0.02}_{-0.03}$ | $0.03^{+0.02}_{-0.01}$ | $0.14^{+0.20}_{-0.04}$ | $1.71^{+0.06}_{-0.08}$ | $5.15^{+0.07}_{-1.28}$
DN Tau | $3520^{+44}_{-27}$ | $3.65^{+0.07}_{-0.06}$ | $0.62^{+0.02}_{-0.02}$ | $0.90^{+0.34}_{-0.38}$ | $1.59^{+0.05}_{-0.06}$ | $11.01^{+0.21}_{-0.21}$
DoAr 25 | $3696^{+58}_{-56}$ | $3.52^{+0.13}_{-0.10}$ | $0.50^{+0.03}_{-0.03}$ | $0.38^{+0.56}_{-0.28}$ | $2.09^{+0.13}_{-0.12}$ | $15.90^{+0.60}_{-0.34}$
DoAr 33 | $3697^{+22}_{-35}$ | $3.86^{+0.04}_{-0.05}$ | $0.70^{+0.01}_{-0.02}$ | $0.16^{+0.25}_{-0.06}$ | $2.16^{+0.07}_{-0.06}$ | $12.68^{+0.35}_{-0.33}$
DoAr 43 SW | $4221^{+67}_{-56}$ | $3.65^{+0.13}_{-0.20}$ | $1.96^{+0.14}_{-0.17}$ | $0.19^{+0.43}_{-0.09}$ | $2.04^{+0.39}_{-0.47}$ | $34.57^{+1.13}_{-2.06}$
DS Tau | $3744^{+66}_{-67}$ | $3.81^{+0.10}_{-0.11}$ | $1.61^{+0.05}_{-0.05}$ | $0.90^{+0.47}_{-0.67}$ | $2.23^{+0.13}_{-0.11}$ | $14.07^{+0.43}_{-0.42}$
FM Tau | $3378^{+68}_{-88}$ | $3.88^{+0.06}_{-0.08}$ | $1.86^{+0.11}_{-0.11}$ | $0.28^{+0.73}_{-0.18}$ | $2.86^{+0.18}_{-0.18}$ | $10.25^{+0.74}_{-1.20}$
FP Tau | $3395^{+82}_{-63}$ | $3.40^{+0.20}_{-0.12}$ | $0.26^{+0.03}_{-0.02}$ | $1.81^{+0.71}_{-0.51}$ | $1.14^{+0.25}_{-0.21}$ | $33.47^{+0.65}_{-0.62}$
FX Tau A | $3552^{+51}_{-53}$ | $3.81^{+0.07}_{-0.08}$ | $0.69^{+0.04}_{-0.03}$ | $1.78^{+0.45}_{-0.43}$ | $0.99^{+0.07}_{-0.07}$ | $5.62^{+0.34}_{-0.19}$
GK Tau A | $3655^{+57}_{-59}$ | $3.61^{+0.12}_{-0.11}$ | $1.74^{+0.05}_{-0.05}$ | $0.90^{+0.41}_{-0.42}$ | $2.23^{+0.13}_{-0.13}$ | $20.35^{+0.55}_{-0.54}$
GM Aur | $3843^{+43}_{-57}$ | $3.91^{+0.06}_{-0.08}$ | $0.59^{+0.02}_{-0.03}$ | $0.17^{+0.31}_{-0.07}$ | $1.99^{+0.10}_{-0.12}$ | $14.06^{+0.34}_{-0.52}$
GO Tau | $3429^{+48}_{-31}$ | $3.86^{+0.05}_{-0.05}$ | $0.28^{+0.03}_{-0.02}$ | $0.93^{+0.34}_{-0.49}$ | $1.48^{+0.09}_{-0.11}$ | $13.44^{+0.48}_{-0.28}$
GSS 26 | $3498^{+87}_{-83}$ | $3.44^{+0.24}_{-0.10}$ | $0.53^{+0.03}_{-0.04}$ | $0.78^{+0.56}_{-0.67}$ | $2.57^{+0.16}_{-0.17}$ | $17.94^{+0.53}_{-0.71}$
GSS 37 W | $3590^{+31}_{-47}$ | $3.45^{+0.12}_{-0.06}$ | $0.47^{+0.03}_{-0.02}$ | $0.91^{+0.44}_{-0.44}$ | $1.87^{+0.10}_{-0.10}$ | $16.73^{+0.43}_{-0.45}$
GSS 39 | $3396^{+63}_{-83}$ | $3.20^{+0.10}_{-0.16}$ | $1.82^{+0.07}_{-0.07}$ | $1.57^{+0.47}_{-0.43}$ | $1.85^{+0.18}_{-0.19}$ | $26.28^{+0.75}_{-0.63}$
Haro 1-16 | $3953^{+109}_{-118}$ | $3.93^{+0.09}_{-0.14}$ | $2.33^{+0.15}_{-0.17}$ | $0.31^{+0.77}_{-0.21}$ | $2.47^{+0.26}_{-0.27}$ | $16.17^{+1.14}_{-1.01}$
Haro 6-13 | $3700^{+96}_{-109}$ | $3.54^{+0.21}_{-0.25}$ | $1.55^{+0.12}_{-0.13}$ | $1.96^{+1.02}_{-1.02}$ | $1.15^{+0.28}_{-0.28}$ | $22.20^{+0.83}_{-0.80}$
Haro 6-33 | $3622^{+62}_{-51}$ | $3.69^{+0.12}_{-0.14}$ | $0.55^{+0.03}_{-0.03}$ | $1.52^{+0.67}_{-0.56}$ | $1.35^{+0.19}_{-0.24}$ | $23.70^{+0.63}_{-0.36}$
HK Tau A | $3509^{+48}_{-46}$ | $3.62^{+0.09}_{-0.10}$ | $0.74^{+0.03}_{-0.03}$ | $2.48^{+0.47}_{-0.41}$ | $1.15^{+0.14}_{-0.14}$ | $22.39^{+0.47}_{-0.26}$
HO Tau | $3427^{+50}_{-45}$ | $4.34^{+0.07}_{-0.09}$ | $1.08^{+0.05}_{-0.05}$ | $2.54^{+0.90}_{-0.53}$ | $0.93^{+0.30}_{-0.29}$ | $22.66^{+0.72}_{-0.61}$
HP Tau | $4119^{+39}_{-39}$ | $3.81^{+0.09}_{-0.11}$ | $1.80^{+0.09}_{-0.08}$ | $0.46^{+0.68}_{-0.36}$ | $2.13^{+0.12}_{-0.14}$ | $17.27^{+0.70}_{-0.51}$
IP Tau | $3659^{+37}_{-41}$ | $4.34^{+0.05}_{-0.05}$ | $1.24^{+0.05}_{-0.04}$ | $2.03^{+0.17}_{-0.03}$ | $2.64^{+0.11}_{-0.11}$ | $12.67^{+0.39}_{-1.23}$
IQ Tau | $3549^{+60}_{-54}$ | $3.60^{+0.12}_{-0.13}$ | $0.91^{+0.04}_{-0.04}$ | $1.72^{+0.53}_{-0.47}$ | $1.85^{+0.11}_{-0.11}$ | $15.59^{+0.33}_{-0.39}$
LkCa 15 | $4093^{+41}_{-45}$ | $4.16^{+0.08}_{-0.19}$ | $1.07^{+0.05}_{-0.06}$ | $2.04^{+0.36}_{-0.49}$ | $1.97^{+0.13}_{-0.13}$ | $15.15^{+0.49}_{-0.89}$
ROX 25 | $4400^{+46}_{-51}$ | $3.95^{+0.22}_{-0.13}$ | $1.05^{+0.19}_{-0.15}$ | $2.56^{+1.10}_{-1.50}$ | $0.71^{+0.57}_{-0.61}$ | $28.81^{+0.60}_{-0.57}$
ROX 27 | $3760^{+51}_{-58}$ | $3.60^{+0.11}_{-0.13}$ | $1.28^{+0.04}_{-0.05}$ | $0.94^{+0.49}_{-0.55}$ | $1.83^{+0.11}_{-0.13}$ | $18.21^{+0.44}_{-0.41}$
UY Aur NE | $3603^{+45}_{-64}$ | $3.42^{+0.08}_{-0.10}$ | $1.23^{+0.05}_{-0.05}$ | $1.06^{+0.58}_{-0.62}$ | $1.75^{+0.12}_{-0.13}$ | $21.03^{+0.40}_{-0.48}$
V710 Tau N | $3484^{+20}_{-35}$ | $3.76^{+0.04}_{-0.05}$ | $0.27^{+0.02}_{-0.02}$ | $1.99^{+0.45}_{-0.41}$ | $1.83^{+0.09}_{-0.09}$ | $17.96^{+0.34}_{-0.44}$
V710 Tau S | $3400^{+10}_{-7}$ | $3.78^{+0.03}_{-0.03}$ | $0.22^{+0.01}_{-0.02}$ | $1.67^{+0.27}_{-0.31}$ | $0.88^{+0.17}_{-0.14}$ | $25.12^{+0.30}_{-0.22}$
WSB 82 | $3968^{+142}_{-169}$ | $3.99^{+0.22}_{-0.44}$ | $4.15^{+0.44}_{-0.70}$ | $3.64^{+0.36}_{-2.08}$ | $3.24^{+0.63}_{-0.75}$ | $32.03^{+0.83}_{-0.96}$
YLW 58 | $3004^{+110}_{-37}$ | $3.10^{+0.37}_{-0.10}$ | $1.68^{+0.11}_{-0.15}$ | $1.06^{+1.00}_{-0.70}$ | $1.18^{+0.15}_{-0.20}$ | $12.14^{+1.00}_{-1.20}$

Note. — The reported uncertainties correspond to 3$\sigma$ deviations from the
median value obtained from the MCMC distributions.

Table 2: Optical parameters from literature and derived masses.

Source | Adopted | $\rm T_{optical}$ | Reference | $\rm M_{\rm K-band}$ | $\rm M_{\rm optical}$ | $\rm M_{\rm dyn}$
---|---|---|---|---|---|---
name | Sp Type | (K) |  | ($\rm M_{\odot}$) | ($\rm M_{\odot}$) | ($\rm M_{\odot}$)
AA Tau | M0.6 | $3792\pm 36$ | (1) | $0.66_{-0.07}^{+0.07}$ | $0.90_{-0.05}^{+0.06}$ | $0.84_{-0.04}^{+0.04}$
BP Tau | M0.5 | $3810\pm 36$ | (1) | $0.73_{-0.05}^{+0.05}$ | $0.83_{-0.05}^{+0.04}$ | $1.10_{-0.04}^{+0.04}$
CI Tau | K5.5 | $4020\pm 48$ | (1) | $1.03_{-0.08}^{+0.07}$ | $1.30_{-0.05}^{+0.04}$ | $0.90_{-0.02}^{+0.02}$
CX Tau | M2.5 | $3485\pm 30$ | (1) | $0.59_{-0.07}^{+0.07}$ | $0.57_{-0.05}^{+0.06}$ | $0.38_{-0.02}^{+0.02}$
CY Tau | M2.3 | $3515\pm 30$ | (1) | $0.45_{-0.05}^{+0.05}$ | $0.57_{-0.04}^{+0.04}$ | $0.30_{-0.02}^{+0.02}$
DE Tau | M2.3 | $3515\pm 30$ | (1) | $0.55_{-0.07}^{+0.06}$ | $0.60_{-0.04}^{+0.05}$ | $0.41_{-0.03}^{+0.03}$
DK Tau A | K8.5 | $3930\pm 30$ | (1) | $0.58_{-0.06}^{+0.07}$ | $1.08_{-0.04}^{+0.04}$ | $0.55_{-0.13}^{+0.13}$
DK Tau B | M1.7 | $3608\pm 32$ | (1) | $0.47_{-0.05}^{+0.05}$ | $0.68_{-0.04}^{+0.05}$ | $1.19_{-0.27}^{+0.27}$
DL Tau | K5.5 | $4163\pm 48$ | (1) | $1.05_{-0.06}^{+0.05}$ | $1.29_{-0.05}^{+0.05}$ | $1.04_{-0.02}^{+0.02}$
DM Tau | M3 | $3410\pm 44$ | (1) | $0.35_{-0.04}^{+0.05}$ | $0.45_{-0.04}^{+0.04}$ | $0.55_{-0.02}^{+0.02}$
DN Tau | M0.3 | $3846\pm 36$ | (1) | $0.63_{-0.08}^{+0.06}$ | $1.01_{-0.05}^{+0.05}$ | $0.87_{-0.15}^{+0.15}$
DoAr 25 | K5 | $4210\pm 262$ | (3) | $0.89_{-0.09}^{+0.09}$ | $1.38_{-0.23}^{+0.13}$ |
DoAr 33 | K5.5 | $4163\pm 130$ | (3) | $0.76_{-0.06}^{+0.07}$ | $1.28_{-0.13}^{+0.10}$ |
DoAr 43 SW | K2 | $4710\pm 166$ | (6) | $1.39_{-0.10}^{+0.02}$ | $1.61_{-0.08}^{>0.09}$ |
DS Tau | M0.4 | $3828\pm 36$ | (1) | $0.83_{-0.08}^{+0.10}$ | $0.95_{-0.06}^{+0.05}$ | $0.83_{-0.02}^{+0.02}$
FM Tau | M0 | $3900\pm 72$ | (2) | $0.42_{-0.07}^{+0.08}$ | $1.01_{-0.09}^{+0.08}$ | $0.36_{-0.02}^{+0.02}$
FP Tau | M2.6 | $3470\pm 30$ | (1) | $0.56_{-0.11}^{+0.09}$ | $0.65_{-0.16}^{+0.04}$ | $0.36_{-0.02}^{+0.02}$
FX Tau A | M2.2 | $3530\pm 30$ | (1) | $0.61_{-0.06}^{+0.07}$ | $0.59_{-0.04}^{+0.04}$ |
GK Tau A | K6.5 | $4068\pm 48$ | (1) | $0.79_{-0.09}^{+0.09}$ | $1.24_{-0.06}^{+0.05}$ | $0.73_{-0.06}^{+0.06}$
GM Aur | K6 | $4115\pm 48$ | (1) | $0.94_{-0.07}^{+0.06}$ | $1.22_{-0.05}^{+0.05}$ | $1.14_{-0.02}^{+0.02}$
GO Tau | M2.3 | $3900\pm 72$ | (2) | $0.47_{-0.05}^{+0.05}$ | $1.02_{-0.09}^{+0.08}$ | $0.45_{-0.01}^{+0.01}$
GSS 26 |  |  |  | $0.68_{-0.18}^{+0.11}$ |  |
GSS 37 W | M1 | $3720\pm 170$ | (4) | $0.79_{-0.10}^{+0.06}$ | $0.93_{-0.24}^{+0.16}$ |
GSS 39 |  |  |  | $0.40_{-0.01}^{+0.18}$ |  | $0.47_{-0.04}^{+0.04}$
Haro 1-16 | K3 | $4543\pm 290$ | (7) | $1.07_{-0.16}^{+0.13}$ | $1.49_{-0.21}^{>0.21}$ |
Haro 6-13 | M0.5 | $3810\pm 90$ | (5) | $0.88_{-0.14}^{+0.12}$ | $1.00_{-0.13}^{+0.10}$ | $0.93_{-0.14}^{+0.14}$
Haro 6-33 | M1.6 | $3620\pm 36$ | (1) | $0.72_{-0.07}^{+0.09}$ | $0.72_{-0.05}^{+0.07}$ | $0.60_{-0.10}^{+0.10}$
HK Tau A | M1.5 | $3640\pm 32$ | (1) | $0.61_{-0.06}^{+0.08}$ | $0.77_{-0.06}^{+0.06}$ | $0.53_{-0.03}^{+0.03}$
HO Tau | M3.2 | $3810\pm 72$ | (2) | $0.43_{-0.05}^{+0.06}$ | $0.80_{-0.07}^{+0.07}$ | $0.43_{-0.03}^{+0.03}$
HP Tau | K4 | $4377\pm 60$ | (1) | $1.25_{-0.05}^{+0.05}$ | $1.43_{-0.04}^{+0.03}$ |
IP Tau | M0.6 | $3792\pm 36$ | (1) | $0.68_{-0.05}^{+0.05}$ | $0.79_{-0.06}^{+0.06}$ | $0.94_{-0.05}^{+0.05}$
IQ Tau | M1.1 | $3704\pm 32$ | (1) | $0.67_{-0.09}^{+0.09}$ | $0.86_{-0.06}^{+0.06}$ | $0.74_{-0.02}^{+0.02}$
LkCa 15 | K5.5 | $4163\pm 48$ | (1) | $1.09_{-0.12}^{+0.08}$ | $1.13_{-0.14}^{+0.09}$ | $1.14_{-0.03}^{+0.03}$
ROX 25 | K6 | $4115\pm 95$ | (4) | $1.41_{-0.11}^{+0.05}$ | $1.21_{-0.10}^{+0.10}$ | $1.27_{-0.07}^{+0.07}$
ROX 27 | K5 | $4210\pm 262$ | (3) | $0.93_{-0.08}^{+0.08}$ | $1.36_{-0.24}^{+0.14}$ |
UY Aur NE | M0 | $3900\pm 90$ | (8) | $0.80_{-0.12}^{+0.06}$ |  |
V710 Tau N | M1.7* | $3608\pm 44$ | (1) | $0.55_{-0.06}^{+0.05}$ | $0.68_{-0.05}^{+0.07}$ | $0.67_{-0.06}^{+0.06}$
V710 Tau S | M3.3* | $3344\pm 32$ | (1) | $0.45_{-0.05}^{+0.05}$ | $0.40_{-0.03}^{+0.03}$ |
WSB 82 |  |  |  | $1.06_{-0.20}^{+0.16}$ |  |
YLW 58 | M4.5 | $3085\pm 210$ | (3) |  |  | $0.10_{-0.01}^{+0.01}$

, (3) is Wilking et al. (2005), (4) is Erickson et al. (2011), (5) is White &
Hillenbrand (2004), (6) is Torres et al. (2006), (7) is Bouvier & Appenzeller
(1992), and (8) is Hartigan & Kenyon (2003).

Note. — (1) is Herczeg & Hillenbrand (2014), (2) is Luhman et al. (2010)

Note. — As discussed in Manara et al. (2019), Herczeg & Hillenbrand (2014)
defined V710 Tau A as the southern companion and V710 Tau B as the Northern
companion.

##  Appendix A Individual Stellar Spectra and Best Fit Models

We present portions of the $K$ band spectrum of 3 sources and their best fit
models in Figures 6, 7, and 8. The figures for the remaining 37 stars can be
accessed in the online journal as Figure sets. The plots display seven panels,
each covering wavelength regions of about 150 Å, that contain the strongest
absorption lines in the spectra of the stars. Green faded boxes delimit the
regions that were used to derive the stellar parameters of the stars (see
Section 3). For visual clarity the spectrum of the stars was smoothed with a
5-pixel wide Savitzky-Golay function of order 2. This smoothing was performed
only for figure display and not during the data fitting process.

The three sources selected for the printed journal, were chosen to display a
range in $K$ band temperature and projected rotational velocities. CY Tau
(shown in Figure 6) has a $K$ band temperature of
$T_{K-band}=3403^{+43}_{-21}$ K, and displays a “forest” of molecular water
absorption lines (not included in our models) which affect mostly, but are not
limited to, the wavelength regions shortward of 2.2 $\micron$. At these cool
temperatures (compare to the rest of the sample), the spectra of the stars are
also characterized by a large line-depth difference between the Si and the Al
lines at 2.179 $\micron$ (second panel from top to bottom), and weaker Fe
lines at 2.2265 $\micron$ compared to the Ti lines at 2.222, 2.224, 2.228, and
2.232 $\micron$ (fourth panel). Hotter stars such as DoAr 43SW (shown in
Figure 7) with $T_{K-band}=4223^{+64}_{-57}$ K, display no molecular
absorption bands from the water molecule, a weaker ratio between Ti and Si at
2.179 $\micron$, and comparatively stronger Fe to Ti lines in the region of
2.221 to 2.232 $\micron$.

The differences in the projected rotational velocity of the sources can be
assessed from almost any absorption line in the spectrum of the stars, but the
unequivocal determination come from the width of the CO lines starting at
$\sim 2.3\,\micron$, and the slope of the CO bandhead at 2.2937 $\micron$
(last panel). The CO absorption lines in the $K$ band are excellent indicators
of the projected rotational velocity of a star as their width are not affected
by changes in the magnetic field strength of the star, and very little
affected by thermal and pressure broadening (e.g., Johns-Krull & Valenti,
2001; Doppmann et al., 2005). In Figure 8, we show the spectrum of FX Tau A,
which display well resolved CO spectral lines, an almost vertical CO overtone
at 2.2937 $\micron$, and has a projected rotational velocity of
$v\sin{i}=5.6^{+0.34}_{-0.19}$ $\rm km\,s^{-1}$. In comparison, the spectrum
of DoAr 43SW (Figure 7) shows broad absorption lines throughout the full $K$
band spectrum, a first CO overtone with a slope that deviates from a vertical
line, and clearly rotationally broaden CO lines beyond 2.296 $\micron$.
Concordantly, the projected rotational velocity measured for this star is of
$v\sin{i}=34.6^{+1.13}_{-2.06}$ $\rm km\,s^{-1}$, which is largest velocity
value measured in our sample.

Figure 6: Comparison between the spectrum of CY Tau (in gray) and its best fit
model (in orange). The seven panels show regions with the strongest
photospheric absorption in lines. Green faded boxes denote the regions used in
the determination of stellar parameters. In each panel, some of the absorption
lines are labeled.  Figure 7: Same as before but for DoAr 43SW

Fig. Set1. Young Stellar Spectra

Figure 8: Same as before but for FX Tau A. The complete figure set of 40
images is available in the online journal.

##  Appendix B Gravity measurement

Unlike other stellar parameters such as temperature, projected rotational
velocity, and magnetic field strength, the change in surface gravity has a
more subtle effect on the spectrum of the stars. However, given the high
spectral resolution of our observations (R$\sim$50000), and the high quality
of our data (S/N$>50$), we can make accurate and precise measurements of this
stellar parameter.

In Figure 9, we show that gravity changes of 0.5 dex in the spectrum of stars
can be distinguished from the depth and strength of particular absorption
lines. As example, we compared the spectrum of GSS 37W, DM Tau, and LkCa 15
with gravities of $\log{g}=3.45^{+0.13}_{-0.06}$, $3.90^{+0.02}_{-0.03}$, and
$4.17^{+0.09}_{-0.17}$, respectively, against MoogStokes models with gravities
ranging from 3.0 to 4.5 in $\log_{10}$ of $\rm cm\,s^{-2}$. For each source,
all the other stellar parameters were fixed to their best fit value obtained
in Section 4. We chose these three sources because they comprise a substantial
range in gravities and $K$ band temperatures, and are relatively slow rotators
$v\sin{i}<18$ $\rm km\,s^{-1}$. Three small spectral regions were chosen to
exemplify the how changes in gravity result in better or worse fit to the
data. A small subplot at the bottom of each wavelength region shows the
squared residuals, but with the sign preserved, i.e, $sign(data-
model)\times(data-model)^{2}$. The residuals were squared as they better
approximate what chi-square calculates.

In the top panel, we see that the purple line ($\log{g}=3.5$) better matches
the spectrum of GSS37 W (black line) in each wavelength region, and the second
closest match is the pink line ($\log{g}=3.0$). These matches are consistent
with the derived gravity for GSS37 W of $\log{g}=3.45$. The same situation
happens in the middle panel, where the model that closest match the data for
DM Tau is the green line with a gravity of $\log{g}=4.0$. For the bottom
panel, the green and blue lines (of $\log{g}=4.0$ and $\log{g}=4.5$,
respectively) closely reproduce the data for LKCa 15, which has a measured
$\log{g}=4.17$. In these panels, regions that allow to better discriminate
between different gravities are the depth and wings of the Na lines at 2.206
$\micron$, and the strength of the Ca line at 2.614 $\micron$ and Ca doublet
at 2.2660 $\micron$.

For additional information on the precision and accuracy our derived
gravities, we refer the reader to Figure 5 of Flores et al. (2019), where we
compared our measurements against literature gravities of main and post-main
sequence stars. It is worth emphasizing, that the modeling approach described
in Section 3 does not only use the lines shown in Figure 9, but all the
information of the green boxes displayed in Figures 6 to 8.

Figure 9: Comparison between MoogStokes models with different gravities vs.
the spectrum of (a) GSS37W, (b) DM Tau, and (c) LkCa 15, in the Na 2.206
$\micron$ (left), the Ca 2.261 $\micron$ (middle), and the Ca 2.265 $\micron$
regions (right). Residuals are squared as they better trace the chi-square
statistic, but the signs are preserved for graphical purposes.

##  Appendix C Magnetic field measurement

One of the fundamental stellar parameters presented in this work is the
magnetic field strengths of the stars. It is key for our interpretation of the
temperature differences as caused by magnetically induced starspots. We
therefore show in Figure 10 examples of how the magnetic field strength of the
stars changes their spectra. We selected four stars with magnetic field
strengths of $\langle B\rangle$ = 0.75, 1.15, 1.60, and 2.16 kG because they
are a representative sample of the values found for the rest of the sources.
The $K$ band temperature range of the selected sources varies between 3450 K
and 3700 K, and the projected rotational velocities change between 10 $\rm
km\,s^{-1}$ and 22 $\rm km\,s^{-1}$. Although the limited sample size makes it
hard to select stars with similar stellar parameters, we emphasize that the
plots shown in the appendix section are just meant to be a demonstration of
our method. In the following paragraphs, we focused on atomic lines with
strong magnetic sensitivity to illustrate that it is possible to detect
magnetic fields as weak as 0.7 kG, and how line shapes can be significantly
altered under the presence of a 2.0 kG field strength.

In Figure 10, we display four sub-figures that correspond to small portions of
the spectrum of DE Tau, HK Tau A, DN Tau, and DoAr 33. The sources were sorted
from weakest to strongest magnetic field strengths, the stellar spectra is
displayed in solid black, the best fit model is shown in orange, and a
comparative spectrum with a null magnetic field strength is displayed as
dotted and thinner black line.

For DE Tau (top row), we see almost no difference between the null magnetic
field strength model (dotted line) and the best fit model with $\langle
B\rangle=0.75$ kG in the first panel, and only a small difference in the last
three panels. For this source, the main observed difference is the width of
the Ti lines for panels three and four, and a slight line depth difference in
the second panel. The fact that we only see little differences between the
models and the observations for this source, reflects that our method is
approaching the detection limit. Indeed in Flores et al. (2019), we used
observations of non-magnetic giant stars to infer that the minimum magnetic
field strength we can measure is $\langle B_{limit}\rangle=0.3$ kG

The case of HK Tau A (second row) is similar to DE Tau, where the null
magnetic field model is not too different from the best fit model for this
source. The largest difference for this star is seen for the Ti lines at
2.2010 and 2.2217 $\micron$ (second and third panel). For last two stars DN
Tau and DoAr 33, which have the strongest magnetic fields, the models with
zero magnetic fields and the best fit models are substantially different in
all four panels. In particular, the Ti line at 2.2317 $\micron$ (fourth panel)
shows a much broader profile, which is reproduced by the model through
magnetic splitting of the line (Zeeman splitting). Although the S/N of the
observations (or perhaps the simplistic one magnetic component assumption in
our models) does not allow us to see the line splitting, the null-field model
with same stellar parameters (temperature, gravity, rotational velocity, etc.)
fails to reproduce both the width and depth of the Ti at 2.2317 $\micron$. In
the case of DoAr 33, even the Ti line at 2.2217 $\micron$ splits (third
panel), which indicates that the magnetic field of the star is significantly
stronger than the other sources displayed here.

Figure 10: Demonstration of the effects of magnetic fields on the spectra of
DE Tau, HK Tau A, DN Tau, and DoAr 33. Four panels that include the highly
magnetically sensitive Ti lines are shown in the figure. The spectra of the
stars are displayed in solid black, the best fit synthetic models are shown in
orange, and comparative synthetic models with null field strengths ($\langle
B\rangle=0$ kG) are displayed as dotted black lines. In all cases, even for
the weakest magnetic field strength detected, there are visible differences
between the magnetic and null field models.

##  Appendix D Stellar parameter uncertainty and degeneracy

In Figure 11, we present the corner plot for CY Tau made from its second MCMC
run. The corner plot for the rest of the sources can be accessed in the online
journal as Figure sets. As explained in Section 3, the second MCMC step does
no include the CO region (starting at $\sim 2.293\,\micron$) neither it fits
the $v\sin{i}$ parameter. For every source, the corner plots only include the
last 50% of the $\sim$50,000 models. We visually assessed the convergence of
the MCMC runs, which typically happens within the first 20$\%$ of each run.
There are small correlation between some stellar parameters, for example, the
$v_{micro}$ and the veiling parameters are weakly correlated for this source.
This makes sense, as both of these parameters are known to affect the depth of
the lines. There are also some minor asymmetries in the marginalized two
dimensional distribution for $v_{micro}$ vs. $\log{g}$, and for $\langle
B\rangle$ vs. $\log{g}$. For some of the stars the $T_{K-band}$ and the
$\log{g}$ can also weakly correlate. The values quoted on top of the 1-D
histograms correspond to the median of the distribution (central black dashed
line), the side black dashed lines represent the 3$\sigma$ uncertainties. Red
lines point to the specific best fit model of the MCMC chain. This best-fit
model is well within the 3$\sigma$ uncertainty for each stellar parameter in
every case.

Figure 11: We present the results for the second MCMC run for CY Tau
(rotational velocity is not included). The corner plot shows only 50% of the
models, i.e., about 12500 models. The central black lines correspond to the
median of the distributions, while the black lines on the side are the
3$\sigma$ uncertainties. The specific best fit model is traced by a red line.

Fig. Set2. Cornerplots of MCMC runs

## References

  * Anderson et al. (1976) Anderson, C. M., Hartmann, L. W., & Bopp, B. W. 1976, ApJ, 204, L51, doi: 10.1086/182052
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42, doi: 10.1051/0004-6361/201425481
  * Berdyugina (2005) Berdyugina, S. V. 2005, Living Reviews in Solar Physics, 2, 8, doi: 10.12942/lrsp-2005-8
  * Bopp & Evans (1973) Bopp, B. W., & Evans, D. S. 1973, MNRAS, 164, 343, doi: 10.1093/mnras/164.4.343
  * Bopp & Stencel (1981) Bopp, B. W., & Stencel, R. E. 1981, ApJ, 247, L131, doi: 10.1086/183606
  * Bouvier & Appenzeller (1992) Bouvier, J., & Appenzeller, I. 1992, A&AS, 92, 481
  * Braun et al. (2021) Braun, T. A. M., Yen, H.-W., Koch, P. M., et al. 2021, ApJ, 908, 46, doi: 10.3847/1538-4357/abd24f
  * Catalano et al. (2002) Catalano, S., Biazzo, K., Frasca, A., et al. 2002, Astronomische Nachrichten, 323, 260, doi: 10.1002/1521-3994(200208)323:3/4<260::AID-ASNA260>3.0.CO;2-9
  * Cieza et al. (2019) Cieza, L. A., Ruíz-Rodríguez, D., Hales, A., et al. 2019, MNRAS, 482, 698, doi: 10.1093/mnras/sty2653
  * Cushing et al. (2004) Cushing, M. C., Vacca, W. D., & Rayner, J. T. 2004, PASP, 116, 362, doi: 10.1086/382907
  * Deen (2013) Deen, C. P. 2013, AJ, 146, 51, doi: 10.1088/0004-6256/146/3/51
  * Doppmann et al. (2005) Doppmann, G. W., Greene, T. P., Covey, K. R., & Lada, C. J. 2005, AJ, 130, 1145, doi: 10.1086/431954
  * Erickson et al. (2011) Erickson, K. L., Wilking, B. A., Meyer, M. R., Robinson, J. G., & Stephenson, L. N. 2011, AJ, 142, 140, doi: 10.1088/0004-6256/142/4/140
  * Esplin et al. (2018) Esplin, T. L., Luhman, K. L., Miller, E. B., & Mamajek, E. E. 2018, AJ, 156, 75, doi: 10.3847/1538-3881/aacce0
  * Feiden (2016) Feiden, G. A. 2016, A&A, 593, A99, doi: 10.1051/0004-6361/201527613
  * Flores et al. (2019) Flores, C., Connelley, M. S., Reipurth, B., & Boogert, A. 2019, ApJ, 882, A75, doi: 10.3847/1538-4357/ab35d4
  * Flores et al. (2020) Flores, C., Reipurth, B., & Connelley, M. S. 2020, arXiv e-prints, arXiv:2006.10139.  https://arxiv.org/abs/2006.10139
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Frasca et al. (2008) Frasca, A., Biazzo, K., Tas, G., Evren, S., & Lanzafame, A. C. 2008, A&A, 479, 557, doi: 10.1051/0004-6361:20077915
  * Gregory et al. (2012) Gregory, S. G., Donati, J. F., Morin, J., et al. 2012, ApJ, 755, 97, doi: 10.1088/0004-637X/755/2/97
  * Guilloteau et al. (2014) Guilloteau, S., Simon, M., Piétu, V., et al. 2014, A&A, 567, A117, doi: 10.1051/0004-6361/201423765
  * Gullbring et al. (2000) Gullbring, E., Calvet, N., Muzerolle, J., & Hartmann, L. 2000, ApJ, 544, 927, doi: 10.1086/317253
  * Gully-Santiago et al. (2017) Gully-Santiago, M. A., Herczeg, G. J., Czekala, I., et al. 2017, ApJ, 836, A200, doi: 10.3847/1538-4357/836/2/200
  * Gustafsson et al. (2008) Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951, doi: 10.1051/0004-6361:200809724
  * Hartigan & Kenyon (2003) Hartigan, P., & Kenyon, S. J. 2003, ApJ, 583, 334, doi: 10.1086/345293
  * Hartmann et al. (2016) Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135, doi: 10.1146/annurev-astro-081915-023347
  * Hayashi (1961) Hayashi, C. 1961, PASJ, 13, 450
  * Herczeg & Hillenbrand (2014) Herczeg, G. J., & Hillenbrand, L. A. 2014, ApJ, 786, A97, doi: 10.1088/0004-637X/786/2/97
  * Jetsu et al. (1993) Jetsu, L., Pelt, J., & Tuominen, I. 1993, A&A, 278, 449
  * Johns-Krull (2007) Johns-Krull, C. M. 2007, ApJ, 664, 975, doi: 10.1086/519017
  * Johns-Krull & Valenti (2001) Johns-Krull, C. M., & Valenti, J. A. 2001, ApJ, 561, 1060, doi: 10.1086/323257
  * Lavail et al. (2019) Lavail, A., Kochukhov, O., & Hussain, G. A. J. 2019, A&A, 630, A99, doi: 10.1051/0004-6361/201935695
  * Lavail et al. (2017) Lavail, A., Kochukhov, O., Hussain, G. A. J., et al. 2017, A&A, 608, A77, doi: 10.1051/0004-6361/201731889
  * Luhman et al. (2010) Luhman, K. L., Allen, P. R., Espaillat, C., Hartmann, L., & Calvet, N. 2010, ApJS, 186, 111, doi: 10.1088/0067-0049/186/1/111
  * Manara et al. (2019) Manara, C. F., Tazzari, M., Long, F., et al. 2019, A&A, 628, A95, doi: 10.1051/0004-6361/201935964
  * Pakhomov et al. (2019) Pakhomov, Y. V., Ryabchikova, T. A., & Piskunov, N. E. 2019, Astronomy Reports, 63, 1010, doi: 10.1134/S1063772919120047
  * Rayner et al. (2016) Rayner, J., Tokunaga, A., Jaffe, D., et al. 2016, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, iSHELL: a construction, assembly and testing, 990884, doi: 10.1117/12.2232064
  * Rothman & Gordon (2010) Rothman, L., & Gordon, I. 2010, in EGU General Assembly Conference Abstracts, EGU General Assembly Conference Abstracts, 5561
  * Ryabchikova et al. (2015) Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015, Phys. Scr, 90, 054005, doi: 10.1088/0031-8949/90/5/054005
  * Schaefer et al. (2009) Schaefer, G. H., Dutrey, A., Guilloteau, S., Simon, M., & White, R. J. 2009, ApJ, 701, 698, doi: 10.1088/0004-637X/701/1/698
  * Simon et al. (2017) Simon, M., Guilloteau, S., Di Folco, E., et al. 2017, ApJ, 844, 158, doi: 10.3847/1538-4357/aa78f1
  * Simon et al. (2019) Simon, M., Guilloteau, S., Beck, T. L., et al. 2019, ApJ, 884, 42, doi: 10.3847/1538-4357/ab3e3b
  * Sokal et al. (2020) Sokal, K. R., Johns-Krull, C. M., Mace, G. N., et al. 2020, ApJ, 888, 116, doi: 10.3847/1538-4357/ab59d8
  * Strassmeier (2009) Strassmeier, K. G. 2009, A&A Rev., 17, 251, doi: 10.1007/s00159-009-0020-6
  * Torres et al. (2006) Torres, C. A. O., Quast, G. R., da Silva, L., et al. 2006, A&A, 460, 695, doi: 10.1051/0004-6361:20065602
  * Vacca et al. (2003) Vacca, W. D., Cushing, M. C., & Rayner, J. T. 2003, PASP, 115, 389, doi: 10.1086/346193
  * Vacca & Sandell (2011) Vacca, W. D., & Sandell, G. 2011, ApJ, 732, 8, doi: 10.1088/0004-637X/732/1/8
  * Vidotto et al. (2014) Vidotto, A. A., Gregory, S. G., Jardine, M., et al. 2014, MNRAS, 441, 2361, doi: 10.1093/mnras/stu728
  * Villebrun et al. (2019) Villebrun, F., Alecian, E., Hussain, G., et al. 2019, A&A, 622, A72, doi: 10.1051/0004-6361/201833545
  * White & Hillenbrand (2004) White, R. J., & Hillenbrand, L. A. 2004, ApJ, 616, 998, doi: 10.1086/425115
  * Wilking et al. (2005) Wilking, B. A., Meyer, M. R., Robinson, J. G., & Greene, T. P. 2005, AJ, 130, 1733, doi: 10.1086/432758
  * Yang & Johns-Krull (2011) Yang, H., & Johns-Krull, C. M. 2011, ApJ, 729, 83, doi: 10.1088/0004-637X/729/2/83


# Blind Asynchronous Over-the-Air
Federated Edge Learning

Saeed Razavikia†, Jaume Anguera Peris†, José Mairton B. da Silva Jr∗,†, and
Carlo Fischione† †School of Electrical Engineering and Computer Science, KTH
Royal Institute of Technology, Stockholm, Sweden
∗Department of Electrical and Computer Engineering, Princeton University, New
Jersey, USA
Email: †{sraz, jaumeap<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Federated Edge Learning (FEEL) is a distributed machine learning technique
where each device contributes to training a global inference model by
independently performing local computations with their data. More recently,
FEEL has been merged with over-the-air computation (OAC), where the global
model is calculated over the air by leveraging the superposition of analog
signals. However, when implementing FEEL with OAC, there is the challenge on
how to precode the analog signals to overcome any time misalignment at the
receiver.

In this work, we propose a novel synchronization-free method to recover the
parameters of the global model over the air without requiring any prior
information about the time misalignments. For that, we construct a convex
optimization based on the norm minimization problem to directly recover the
global model by solving a convex semi-definite program. The performance of the
proposed method is evaluated in terms of accuracy and convergence via
numerical experiments. We show that our proposed algorithm is close to the
ideal synchronized scenario by $10\%$, and performs $4\times$ better than the
simple case where no recovering method is used.

###### Index Terms:

Asynchronous, federated edge learning, over-the-air computation, time
misalignment

##  I Introduction

The popularity of mobile devices and the evolution of sensing and machine
learning have led to rapid access to mobile data for training artificial
intelligent models. Yet, the variability and the distributed nature of the
data do not guarantee that the inference model obtained locally from one
device is able to perform well on data from other devices. To allow for a more
collaborative approach, the devices can offload and combine their local data
in the cloud, but this solution has been regarded as impractical, mainly due
to the communication costs, data privacy concerns, and hardware variability of
the devices [1]. A more prominent solution is federated learning, in which a
centralized node periodically receives and combines the local models from
multiple devices to obtain a more accurate global model.

Indeed, the above solutions have gathered interest for federated edge learning
(FEEL), which combines federated learning with edge intelligence. Here, the
devices periodically offload their models to a server located at the edge of
the network, such that these are combined in a federated-learning fashion to
obtain a global model that performs well across multiple devices. Specially,
FEEL leverages the computational resources at the edge server while minimizing
the communication costs at the network [2].

Yet, the wireless resources located at the edge may become limited as the
number of mobile devices participating in FEEL increases or as the dimension
of the local gradients increases [3]. To alleviate the communication costs and
accelerate the training process, a new paradigm in communication and
computation, called over-the-air computation (OAC), concurrently allocates all
users over the same frequency resources and leverages the inherent waveform
superposition of wireless signals. Specially, by appropriate precoding the
analog transmitted signals, it is possible to calculate over the air the class
of functions known as nomographic functions [4], such as the arithmetic mean
or the weighted sum. This joint communication-and-computation scheme is
particularly useful for FEEL, as the edge server is only interested in
constructing a global model from the weighted average of the local models from
all devices, so calculating this aggregation over the air using the
superposition of analog signals results in much more efficient use of the
communication resources and a faster federated-learning system than
transmitting the digital signals, decoding them at the receiver, and
calculating the aggregation later [5, 6].

However, there are many challenges upon implementing an OAC system. Recent
works have made significant progress on this topic regarding the signal
processing [5], the design of the transceivers [7], or the acquisition of
channel state information [8]. Nevertheless, there are challenges that are
still largely unsolved, such as how to accurately precode the transmitted
signals and how to synchronize the clock of all the devices at the level of
nanoseconds to maximize the number of participating devices while minimizing
the distortions introduced by the channel or the devices themselves. These
challenges are particularly important because, in practice, signal precoding
is imperfect due to inaccurate channel estimation and non-ideal hardware, and
the synchronization across different devices is costly because it requires all
devices to time their transmissions accurately based on the propagation time.
Therefore, in a realistic scenario, there appears some channel-gain
mismatches, time asynchronies, or both [9].

To address these issues, the literature offers different solutions to the
misalignment of OAC both for digital and analog systems. Early works [10]
propose different techniques to achieve clock synchronization using the IEEE
802.11 protocol. These techniques are robust to propagation delays on the
channel or within the devices, but they are susceptible to unreliable network
transmissions. The work in [11] implements an analog-modulated system based on
coarse block-synchronization to be robust against synchronization errors.
Specially, their model is optimal over medium-access strategies that separate
the transmissions of nodes in time or code space but not in frequency. The
work in [12] utilizes digital modulation and channel coding to combat any
channel or time misalignments at the receiver, but it still requires decoding
the transmitted signals, thus not fully leveraging the communication benefits
of analog signals with OAC. The most related work up to date is [9], which
tackles the problem of misalignment of analog OAC by constructing matched
filters at the receiver, but it assumes that the sampling process knows the
exact time misalignments from the different users.

$\bm{\Delta w}_{1}^{(m)}=\nabla
F_{1}(\bm{w}^{(m)})$$\bm{x}_{1}=\bm{r}(t)\bm{\Delta
w}_{1}^{(m)}$$\bm{w}^{(m)}=\bm{w}^{(m-1)}-\eta\bm{\Delta w}^{(m-1)}$Local
Dataset$\bm{\Delta w}_{K}^{(m)}=\nabla
F_{K}(w^{(m)})$$\bm{x}_{K}=\bm{r}(t)\bm{\Delta
w}_{K}^{(m)}$$\bm{w}^{(m)}=\bm{w}^{(m-1)}-\eta{\bm{\Delta w}}^{(m-1)}$Local
Dataset$\vdots$$\bm{\Delta w}^{(m)}=\sum_{k}\bm{\Delta
w}_{k}^{(m)}/K$$\bm{w}^{(m+1)}=\bm{w}^{(m)}-\eta\bm{\Delta w}^{(m)}$

Figure 1: Diagram of federated edge learning with over-the-air computation

###  I-A Our contribution

In this work, we focus on the time misalignment of analog OAC and present a
novel approach to recovering the arithmetic sum of the transmitted symbols. In
particular,

  * •

We do not consider any synchronization overhead, such that the receiver has no
prior information about the time misalignments of any participating user.

  * •

We consider an analog OAC scheme where all devices transmit each parameter of
their local model in an uncoded fashion and over a fading, multiple-access
wireless channel, such that the receiver sees the superposition of
asynchronous signals.

  * •

We reformulate the measurement at the receiver in terms of atoms and propose a
convex optimization problem to recover the delays and the elements of the
gradient, such that we obtain the summation to update the global gradient.

  * •

We evaluate our proposed algorithm against the MNIST dataset and show how the
proposed system achieves competitive results in this task, closing the
performance gap down to $10\%$ with respect to the ideal synchronized
scenario, and being $4\times$ better than the simple case where no recovery
method is implemented.

###  I-B Document organization

The remainder of this paper is organized as follows. Section II presents the
learning and communication models for the proposed OAC system. In Section III,
we present an optimization problem to recover the arithmetic sum of the
transmitted symbols. Then, we provide the numerical results in Section IV,
followed by the concluding remarks in Section V.

###  I-C Notation

Throughout this paper, scalars are denoted by lower-case letters $x$, and
vectors and matrices are denoted by lower-case $\bm{x}$ and upper-case
boldface letters $\bm{X}$, respectively. The transpose and Hermitian of a
matrix $\bm{X}$ are represented by $\bm{X}^{\mathsf{T}}$ and
$\bm{X}^{\mathsf{H}}$, respectively. We further use
$f(t)~{}\circledast~{}g(t)$ to show the convolution between two continuous
signals $f(t)$ and $g(t)$. For a set $S$, its cardinally is represented by
$|S|$. For an integer $N$, $[N]$ stands for $\\{1,2,\dots,N\\}$.
$\bm{X}\succeq\bm{0}$ means that $\bm{X}$ is a positive semidefinite matrix.
Finally, the pseudo-inverse of matrix $\bm{X}$ is denoted by
$\bm{X}^{\dagger}$, and the element-wise inverse is denoted by
$(\cdot)^{\circ-1}$,

##  II System Model

This section is divided into two parts. First, we present the federated
learning model, and then we describe the communication model for transmitting
the gradients of the inference model in both the uplink and the downlink.

###  II-A Learning model

Following the same structure as in [13, 2, 14, 15, 16], consider a federated
learning scenario with $K$ devices. For each device $k$, let $\mathcal{D}_{k}$
denote its local dataset, and let the function $F_{k}(\bm{w})$ represent the
average empirical loss at device $k$ with respect to the model parameters
$\bm{w}\in\mathbb{R}^{N}$, with $N$ denoting the number of parameters, such
that

$\displaystyle
F_{k}(\bm{w})=\frac{1}{|\mathcal{D}_{k}|}\sum_{j\in\mathcal{D}_{k}}f_{j}(\bm{w}),$
(1)

where $f_{j}(\bm{w})$ represents the empirical loss function at the data
sample $j$ of the $k$-th local dataset $\mathcal{D}_{k}$. Upon calculating
$F_{k}(\bm{w})$, each device offloads their local model to the edge server for
it to construct the global loss function of the model vector $\bm{w}$ from
each $\mathcal{D}_{k}$ as

$\displaystyle
F(\bm{w})=\frac{1}{|\mathcal{D}|}\sum_{k}|\mathcal{D}_{k}|F_{k}(\bm{w})$ (2)

where $\mathcal{D}=\sum_{k}|\mathcal{D}_{k}|$ represents the total number of
samples from all devices. With that information, the server trains the global
model by minimizing the following empirical cost function in a distributed
manner

$\displaystyle\bm{w}^{*}=\underset{\bm{w}}{\rm argmin}~{}F(\bm{w}).$ (3)

With the idea of preserving privacy, consider further a FEEL framework where
each device uses its local dataset $\mathcal{D}_{k}$ to perform stochastic
gradient descent (SGD) to minimize their local loss function $f_{j}(\bm{w})$.
More specifically, let $\bm{w}^{(m)}\in\mathbb{R}^{N}$ and $\bm{\Delta
w}_{k}^{(m)}\in\mathbb{R}^{N}$ represent the parameters of the gradient and
the estimate of the gradient of device $k$ at the $m$-th communication round,
respectively. Then,

$\displaystyle\bm{\Delta
w}_{k}^{(m)}=\frac{1}{|\mathcal{D}_{k}|}\sum_{j\in\mathcal{D}_{k}}\nabla
f_{j}(\bm{w}^{(m)})$ (4)

where $\nabla$ denotes the gradient operator. After this computation, the
$k$-th device offloads its local model update $\bm{\Delta w}_{k}^{(m)}$ to the
server. Upon receiving the information, the server calculates the global model
of the gradient $\bm{\Delta w}^{(m)}$ of the loss function $F(\bm{w})$ from
all $K$ users as

$\displaystyle\bm{\Delta w}^{(m)}=\frac{1}{K}\sum_{k=1}^{K}\bm{\Delta
w}_{k}^{(m)}.$ (5)

Finally, the server updates the current global model following the gradient
descent as

$\displaystyle\bm{w}^{(m+1)}=\bm{w}^{(m)}-\eta\bm{\Delta w}^{(m)},$ (6)

where $\eta$ is the learning rate. The global model is then broadcast back to
the devices, and the same procedure is repeated until the model converges to a
local minimum. Based on (5), we observe that the edge server only needs the
aggregation of the local estimations $\bm{\Delta w}_{k}^{(m)}$, but it does
not need the individual gradients from each users. Therefore, this aggregation
can be calculated over-the-air following an OAC scheme. Finally, and for
simplicity, since the communication procedure is the same for every round, we
omit the index $m$ from $\bm{\Delta w}^{(m)}$ and show the gradient at each
iteration as $\bm{\Delta w}$.

###  II-B Communication model

In each communication round, consider all devices to transmit a vector
$\bm{x}_{k}\in\mathbb{R}^{N}$ simultaneously to the edge server over a
broadband multi-access channel (MAC). In a typical communication system, since
the communications are not ideal because of the characteristics of the
receiver and the fact that all devices are not synchronized with each other,
let the channel introduce noise, and let a time misalignment appear between
the different received signals. Altogether, we express the received signal at
the edge server at time $t$ as

$\displaystyle\bm{y}(t)=\sum_{k=1}^{K}h_{k}\bm{x}_{k}(t-\tau_{k})+\bm{z}(t).$
(7)

The coefficients $h_{k}$ originate from transmitting the symbols from the
$k$-th device to the edge server over a fading channel, which are assumed to
be Rayleigh distributed and known at the transmitter. The parameters
$\tau_{k}$ represent the unknown delays of $k$-th device, and are considered
to take arbitrary continuous values in $[0,T]$. Finally, the additive term
$\bm{z}(t)$ represents the white Gaussian noise (AWGN), which is considered to
be distributed according to $\mathcal{N}(\bm{0},\sigma_{z}^{2}\bm{I})$.

Considering this uplink transmission model, we are interested in designing a
transmitter and a receiver that reverts the effect of the channel and allows
all $K$ devices to transmit all $N$ parameters from $\bm{\Delta w}_{k}$ at
once, such that the receiver recovers (5) from the received signals. For the
resource allocation, consider $\bar{B}$ to be the total available bandwidth in
our system, and let $B=\bar{B}/N$ be the bandwidth allocated to simultaneously
transmit one of the elements of the gradient $\bm{\Delta w}_{k}$ for all
devices $k\in[K]$111For large $N$, the whole uplink transmission can be
separated in different batches, such that all $N$ parameters are transmitted
to the server over different time periods.. Then, consider that we construct
$N$ orthogonal band-limited signals $r_{i}(t)$222The waveform $r(t)$ can be
constructed over the band $[-B/2,B/2]$ using the sinc kernel ${\rm sinc}(Bt)$.
Specially, since the sinc kernel decays relatively fast, $r(t)$ would be
approximately supported on the interval $[-3T/2,3T/2]$. for $i\in[N]$ to
simultaneously transmit the $i$-th element of the gradient ${\Delta
w}_{i,k}:=[\bm{\Delta w}]_{i,k}$ of all devices $k\in[K]$ over the $i$-th
frequency band, with

$R_{i}(f)=\mathcal{F}\\{r_{i}(t)\\}=0,\quad\forall|f|>B/2,$

where $\mathcal{F}\\{\cdot\\}$ denotes the Fourier transform.

At the transmitter side, let each device $k$ construct a positive
representation of the gradient ${\Delta w}_{i,k}$ as

${\Delta w}_{i,k}^{+}:={\Delta w}_{i,k}+\gamma\geq 0.$

where $\gamma$ is selected sufficiently large to satisfy the inequality
above333The parameter $\gamma$ can be selected before initiating the FEEL
process and be known at the receiver, so it suffices to make $\gamma$ large
enough to satisfy the inequality, with the only condition that it needs to be
larger than the possible minimum of all gradients, i.e.,
$\gamma\geq\min_{i,k}{\Delta w}_{i,k}$. This minimum can also be shared by the
edge server, in which case the overhead of transmitting one scalar per device
is negligible. With that, each device encodes the information of the
gradient444This encoding step also includes quantization, so the encoded
information in the transmitted signal loses accuracy. Here we do not consider
this loss in accuracy, but one could refer to [17] for more information about
the effect of the quantization error on the inference problem. into the signal
$\bm{x}_{k}(t)=\bm{\Delta
w}_{k}^{+}\frac{h^{*}_{k}}{|h_{k}|}\sqrt{p_{k}}\bm{r}(t)$ from the known
Rayleigh fading coefficients $h_{k}$ and the $N\times 1$ vector signal
$\bm{r}(t)$ whose elements $r_{i}(t)$ satisfy $r_{i}(t)\circledast r_{j}(t)=0$
for any $i\neq j$. Here we assume that all devices compensate for the effect
of the fading $h_{k}$ by setting $p_{k}=|h_{k}|^{-2}$; otherwise it would
result in weaker signal strength.

$\Delta
w_{k,1}$$\bm{Y}_{1}$$1/B$QuantizationandEncodingModulation$\bm{r}(t)$Carrier$\Delta
w_{k,N}$CarrierMatched-filter$\bm{r}^{*}(-t)$$\bm{Y}_{N}$$\vdots$$\vdots$

Figure 2: The block diagram of the communication system in the over-the-air
computation problem. The left side shows the architecture of the transmitter
of $k$-th node, and right side architecture of the transceiver is depicted

The task of the receiver is then to estimate the global gradient $\bm{\Delta
w}^{+}$ from the received signal (7). From the description above, the $i$-th
element of the received signal ${y}_{i}(t):=[\bm{y}]_{i}(t)$ can be expressed
as

$\displaystyle{y}_{i}(t)=\sum_{k=1}^{K}{\Delta
w}_{i,k}^{+}\,{r}_{i}(t)+z_{i}(t),\quad\forall i\in[N].$ (8)

After the signal goes through a matched filter, we obtain

$\displaystyle\tilde{y}_{i}(t)$ $\displaystyle=y_{i}(t)\circledast
r_{i}^{*}(-t)$ $\displaystyle=\Big{(}\sum_{k=1}^{K}{\Delta
w}_{i,k}^{+}\,{r}_{i}(t-\tau_{k})+z_{i}(t)\Big{)}\circledast r_{i}^{*}(-t)$
$\displaystyle=\sum_{k=1}^{K}{\Delta
w}_{i,k}^{+}\,g_{i}(t-\tau_{k})+\tilde{z}_{i}(t),$ (9)

where $g(t)$ represents the convoluted waveform at time $t$. Since
$\tilde{y}_{i}(t)$ is band-limited and approximately time-limited, the
receiver samples $\tilde{y}_{i}(t)$ at rate $1/B$ in the interval $[-T/2,T/2]$
to collect all its degree of freedom $L=BT=:2M+1$ samples. If we
mathematically manipulate (9) by applying the discrete Fourier transform (DFT)
and the inverse DFT, the $\ell$-th sample of the signal $\tilde{y}_{i}(t)$,
denoted by $Y_{i}(\ell)$, in the interval $\frac{\ell}{B}\in[-T/2,T/2]$ can be
expressed by [18]

$\displaystyle Y_{i}(\ell)=\sum_{k=1}^{K}{\Delta
w}_{i,k}^{+}\;\theta_{i,k}(\ell)+Z_{i}(\ell),$ (10)

where $Z_{i}(\ell)=\tilde{z}_{i}(\tfrac{\ell}{B})$,

$\theta_{i,k}(\ell)=\sum_{q=-M}^{M}[\bm{a}(\tau_{k})]_{q}\;g_{i}{(\tfrac{\ell}{B})}\,e^{\frac{j2\pi\ell
q}{L}},$ (11)

and

$[\bm{a}(\tau_{k})]_{q}=\frac{1}{L}\sum_{r=-M}^{M}e^{\tfrac{j2\pi(q-L\tau_{k})r}{L}},$
(12)

with $\bm{a}(\tau)\in\mathbb{C}^{L}$ being the atom or building block [19].
Without loss of generality, we assume all $g_{i}(t)$ to be periodic with
duration $T$555 We assume periodicity to obtain (10). Clearly, this assumption
is not practical because the waveform signal is not time-limited. However, we
can consider the quasi-periodic waveform satisfying $g(\frac{\ell}{B})=0$ for
$\ell\notin[-\lfloor L/2\rfloor,\lfloor L/2\rfloor]$., $T_{\rm max}=1$, and
$\tau_{k}\in[0,1)$ for $k\in[K]$. We further map the frequency indices $\ell$
from $-\lfloor B/2\rfloor,\ldots,0,\ldots,\lfloor B/2\rfloor$ to
$0,\dots,L-1$, where $L=2\lfloor B/2\rfloor+1$ is the total number of
frequency samples. Furthermore, to alleviate the notation, we rewrite (10) as

$\displaystyle\bm{Y}_{i}=\bm{G}_{i}\bm{X}_{i}+\bm{Z}_{i},$ (13)

where $\bm{X}_{i}=\sum_{k=1}^{K}{\Delta
w}_{i,k}^{+}\,\bm{a}(\tau_{k})\in\mathbb{C}^{L}$,
$\bm{G}_{i}\in\mathbb{C}^{L\times L}$ is the waveform matrix whose
$(\ell,q)$-th element is defined as
$[\bm{G}_{i}]_{(\ell,q)}=g_{i}(\tfrac{\ell}{B})e^{\frac{j2\pi\ell q}{L}}$,
$\bm{Z}_{i}=[Z_{i}(1),\ldots,Z_{i}(L)]^{\mathsf{T}}$, and
$\bm{Y}_{i}=[Y_{i}(1),\ldots,Y_{i}(L)]^{\mathsf{T}}$. Altogether, the sampling
scheme is depicted in Fig. 2. Since the waveform $r_{i}(t)$ is known, we can
remove the effect of $g_{i}(t)$ in (13) by multiplying the measurements to the
Fourier transform inverse $(\bm{G}_{i})^{-1}$ (assuming $\bm{G}_{i}$ is
invertible for any $i\in[N]$). Then, we observe

$\displaystyle\bm{V}_{i}=\bm{X}_{i}+\hat{\bm{Z}}_{i},\quad\forall i~{}\in[N],$
(14)

where $\bm{V}_{i}=(\bm{G}_{i})^{-1}\bm{Y}_{i}$, and $\hat{\bm{Z}}_{i}$ is the
AWGN additive noise with bounded variance $\sigma_{z}^{2}$. Considering this,
the next section presents a method to perfectly recover $\sum_{k=1}^{K}{\Delta
w}_{i,k}$ from $\bm{V}_{i}$.

##  III Recovery Method

The problem presented in Section II and summarized in (14) is closely related
to the classical problem of line spectrum estimation or super resolution
problem [20], in which we seek to estimate all the delays and the amplitudes
$\\{\tau_{k},{\Delta w}_{i,k}^{+}\\}_{k=1}^{K}$ from a mixture of sinusoids.
The main difference is that the edge server is interested in the mean
$\tfrac{1}{K}\sum_{k=1}^{K}{\Delta w}_{i,k}^{+}$ instead of each individual
component ${\Delta w}_{i,k}^{+}$. Our proposed solution is to estimate
directly the mean of gradients $s_{i}:=\tfrac{1}{K}\sum_{k=1}^{K}{\Delta
w}_{i,k}^{+}$ from $\eqref{eq:Measurment}$ using the notion of atomic norm of
the signal $\bm{X}_{i}$.

To begin, define the atomic set as

$\mathcal{A}:=\\{\bm{a}(\tau):\tau\in[0,1)\\},$ (15)

and the associated Minkowski functional over the set $\mathcal{A}$ as [19]

$\displaystyle\|\bm{X}\|_{\mathcal{A}}:=\inf_{\begin{subarray}{c}\beta_{j}\in\mathbb{R}^{+}\\\
\tau_{j}\in[0,1)\end{subarray}}\Big{\\{}\sum_{j}\beta_{j}~{}:~{}\bm{X}=\sum_{j}\beta_{j}\bm{a}(\tau_{j})\Big{\\}}.$
(16)

Finding the optimal parameters in (16) is not an easy task because it involves
an infinite-dimensional variable optimization due to the continuity of the set
$\mathcal{A}$. Alternatively, we can rewrite (16) using the Carathéodory-
Fejér-Pisarenko decomposition [21] as the following semi-definite program

$\displaystyle\|\bm{X}\|_{\mathcal{A}}=\min_{\bm{u}\in\mathbb{C}^{L},t>0}{\;\frac{1}{2}}(\mathsf{Re}(u_{1})+t)$
$\displaystyle{\rm s.t.}~{}~{}\begin{bmatrix}{\rm
Toep}(\bm{u})&\bm{F}^{-\mathsf{H}}\bm{X}\\\
\bm{X}^{\mathsf{H}}\bm{F}^{-1}&t\end{bmatrix}\succeq\bm{0},$ (17)

where $\bm{F}$ denotes the DFT matrix whose elements are given by
$[\bm{F}]_{(k,l)}:=\frac{1}{L}{\rm e}^{-{\rm j}2\pi lk/L}$, and the Hermitian
Toeplitz matrix ${\rm Toep}(\bm{u})$ for a vector $\bm{u}\in\mathbb{C}^{L}$,
is defined as

$\displaystyle{\rm
Toep}(\bm{u}):=\left[\begin{matrix}u_{1}&u_{2}&\dots&u_{L}\\\
u_{2}^{*}&u_{1}&\dots&u_{L-1}\\\ \vdots&\vdots&\ddots&\vdots\\\
u_{L}^{*}&u_{L-1}^{*}&\dots&u_{1}\end{matrix}\right].$ (18)

Hence, the optimization problem to jointly recover the delays $\tau_{k}$ and
the gradients ${\Delta w}_{i,k}^{+}$ can be obtained by searching for signals
$\bm{X}_{i}$ that are both sparse in the continuous atom set $\mathcal{A}$
(small atom norms) and close to the observation $\bm{V}_{i}$, i.e.,

$\displaystyle\min_{\bm{X}_{i}}\|\bm{X}_{i}\|_{\mathcal{A}}+\lambda\|\bm{X}_{i}-\bm{V}_{i}\|_{2}^{2}.$
(19)

Here, the regularization parameter $\lambda>0$ needs to be chosen
appropriately. By employing the atomic decomposition, (19) can be reformulated
as follows

$\displaystyle\hat{\bm{X}}_{i}=\min_{\begin{subarray}{c}\bm{x},t>0\\\
\bm{u}\in\mathbb{C}^{L}\end{subarray}}\;\frac{1}{2}(\mathsf{Re}(u_{1})+t)+\lambda\|\bm{X}_{i}-\bm{V}_{i}\|_{2}^{2}$
$\displaystyle{\rm s.t.}~{}~{}\begin{bmatrix}{\rm
Toep}(\bm{u})&\bm{F}^{-\mathsf{H}}\bm{X}_{i}\\\
\bm{X}_{i}^{H}\bm{F}^{-1}&t\end{bmatrix}\succeq\bm{0}.$ (20)

This convex optimization problem can be solved efficiently using standard
tools from convex optimization [22].

After obtaining the Toeplitz matrix ${\rm Teop}(\bm{u})$, we employ the
Vandermonde decomposition via solving a generalized eigenvalue problem [23] as

${\rm Toep}(\bm{u})=\sum_{k=1}^{K^{\prime}}{\Delta
w}_{i,k}^{+^{\prime}}\bm{a}(\tau_{k}^{{}^{\prime}})\bm{a}(\tau_{k}^{{}^{\prime}})^{\mathsf{H}}$

to identify the support $\\{\tau_{k}^{{}^{\prime}}\\}_{k=1}^{K^{\prime}}$ as
well as the atomic norm
$\|\hat{\bm{X}}_{i}\|_{\mathcal{A}}=\sum_{k=1}^{K}{\Delta
w^{\prime}}_{i,k}^{+}$. With that, the amplitudes, and subsequently the mean
value
$\hat{s}_{i}=\tfrac{1}{K}\|\hat{\bm{X}}_{i}\|_{\mathcal{A}}=\tfrac{1}{K}\sum_{k=1}^{K}{\Delta
w}_{i,k}^{+}-\gamma$, are obtained. Finally, the average gradient is estimated
from $\widehat{\bm{\Delta w}}=[\hat{s}_{1},\ldots,\hat{s}_{N}]^{\mathsf{T}}$.

######  Remark 1.

It is established that problem (III) can recover the positive mean value
$s_{i}$ without requiring any separation between the delays [24], as long as
$K\leq\lfloor\frac{L-1}{2}\rfloor$, and the estimated parameters $\hat{s}_{i}$
satisfy

$\displaystyle|s_{i}-\hat{s}_{i}|=\mathcal{O}\bigg{(}\sigma_{z}\sqrt{\frac{\log{L}}{L}}\bigg{)},\quad\forall
i\in[N]$ (21)

where $\sigma_{z}$ denotes the variance of noise and $L$ is the number of
samples.

##  IV Numerical Results

This section examines the performance of our proposed algorithm in terms of
the normalized mean squared error (NMSE) between the true sum of models and
the estimated one for different numbers of users $K$ and Fourier samples $L$.
We also evaluate our method against the MNIST dataset and show how the
proposed system can achieve competitive results in this task while using a
simple recovery method. In all the experiments, the channel coefficients
$\\{h_{k}\\}_{k=1}^{K}$ are i.i.d generated uniformly on the unit sphere, and
the delays $\\{\tau_{k}\\}_{k=1}^{K}$ are uniformly distributed at random
between $0$ and $1$, i.e., $\tau_{k}\sim\mathcal{U}[0,1)$. Moreover, all
waveforms $r_{i}(t)$ are randomly generated using
$\sum_{\ell=-L-M}^{L+M}a_{i}{\rm sinc(tB-\ell)}$ with
$a_{i}\sim\mathcal{N}(0,1/L)$ for all $i\in[N]$. The optimization problem in
(III) is implemented using the SDPT3 package of CVX in MATLAB.

$4$$6$$8$$10$$12$$14$$16$$18$$20$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$SNR(dB)NMSE$L=8$$L=16$$L=32$$L=64$$L=128$
(a) $K=5$ users

$4$$6$$8$$10$$12$$14$$16$$18$$20$$0$$0.1$$0.2$$0.3$SNR(dB)NMSE$K=50$$K=40$$K=30$$K=20$$K=10$
(b) $L=128$ samples

Figure 3: NMSE as a function of the SNR (dB) for (a) different number of
Fourier samples $L$, and (b) different number of users $K$ by averaging over
$100$ Monte Carlo trials. The results show the effect of the number of users
on the performance of the proposed algorithm for estimating the parameters of
the global inference model using OAC with unknown communication delays.

For the first experiment, we check the performance of the proposed algorithm
for estimating the summation of each individual element of the gradient vector
${\Delta w}_{i}$ for different signal to noise ratio (SNR), defined as ${\rm
SNR}=20\log{(\tfrac{\|\bm{G}\bm{X}\|_{2}}{\|\bm{Z}\|_{2}})}$. Figure 3 (a)
shows the NMSE for $K=5$ users over different SNRs within $[4,20]$ dB and over
different number of Fourier samples $L=\\{8,16,32,64,128\\}$. As we can
observe, the error decreases for increasing SNR and for increasing number of
samples. We further repeat the simulation to check the effect of the number of
users on the recovery problem. The results are shown in Figure 3 (b) for
$L=128$ samples. As the number of users increases, the NMSE error also
increases, which shows the need for more samples. Specially, the error
aggravates for a large number of users, e.g., $K\geq 40$. This comes from the
fact that having more users results in two or more parameters being close to
each other, thus making it more difficult to satisfy the conditions from
Remark 1. Consequently, the optimization in (III) cannot perfectly separate
the amplitude and delays of each individual device $k$, resulting in a worse
NMSE.

For the last experiment, we evaluate the performance of the estimation for the
FEEL problem for $K=10$ devices when the task is a multi-label classification
on the MNIST dataset [25]. We used a simple neural network that consists of an
input layer of $784$ nodes, a hidden layer of $100$ nodes, and an output layer
of $10$ nodes. We compare the performance of our proposed method with the
standard method when the communication system uses ideal versus imperfect
synchronization with $\text{SNR}=5$ dB. Figure 4 depicts the accuracy for
these cases with different number of samples $L=128$ and $L=256$.  Note that
the accuracy of our proposed algorithm increases for increasing number of
samples, closing the gap with respect to the ideal synchronized scenario down
to $10\%$, and being up to $4$ times better than just recovering the global
model directly from the received data without doing any processing. Overall,
the results suggest that our proposed method can successfully learn the global
model even in a blind imperfect synchronization scenario.

$20$$40$$60$$80$$100$$0.2$$0.4$$0.6$$0.8$$1$Communication RoundsAccuracyIdeal
synchImperfect synch with $L=256$Imperfect synch with $L=128$Imperfect synch
Figure 4: Accuracy of the MNIST task as a function of the communication rounds
for $K=10$ users and SNR $5$ dB. The results show the performance of our
proposed recovery method for imperfect and perfect synchronization.

##  V Conclusions

In this paper, we consider the FEEL problem where each device contributes to
training a global inference model by independently performing local
computations with their data. We explained that typical over-the-air
computation methods usually make the idealistic assumption that there is
perfect synchronization between the devices and the receiver. Other existing
methods solve this misalignment problem, but they assume the delays to be
known. Here, we instead focus on the OAC problem considering that there is no
prior information on channel delays over the AWGN channel. We proposed a novel
synchronization-free method to recover the global inference model without
requiring any prior information about the devices’ delays. For that, we
developed an atomic norm minimization problem in order to recover the
summation of the gradient by solving a convex semi-definite program.  Finally,
we evaluated the performance of the recovery method in terms of accuracy and
convergence via numerical experiments, and showed that our proposed system is
close to the ideal synchronized scenario by $10\%$, and performs $4\times$
better than using no recovery methods.

## References

  * [1] T. Li et al., “Federated learning: Challenges, methods, and future directions,” IEEE Sig. Proc. Mag., vol. 37, no. 3, pp. 50–60, 2020.
  * [2] B. McMahan et al., “Communication-efficient learning of deep networks from decentralized data,” in Artificial Intelligence and Statistics, pp. 1273–1282, PMLR, 2017.
  * [3] A. Tak and S. Cherkaoui, “Federated edge learning: Design issues and challenges,” IEEE Network, vol. 35, no. 2, pp. 252–258, 2020.
  * [4] M. Goldenbaum, H. Boche, and S. Stańczak, “Nomographic functions: Efficient computation in clustered gaussian sensor networks,” IEEE Trans. Wireless Commun., vol. 14, no. 4, pp. 2093–2105, 2014.
  * [5] M. Goldenbaum, H. Boche, and S. Stańczak, “Harnessing interference for analog function computation in wireless sensor networks,” IEEE Trans. Sig. Proc., vol. 61, no. 20, pp. 4893–4906, 2013.
  * [6] H. Hellström et al., “Wireless for machine learning: A survey,” Foundations and Trends® in Sig. Proc., vol. 15, no. 4, pp. 290–399, 2022.
  * [7] L. Chen et al., “Over-the-air computation for IoT networks: Computing multiple functions with antenna arrays,” IEEE Internet of Things J., vol. 5, no. 6, pp. 5296–5306, 2018.
  * [8] F. Ang et al., “Robust design for massive CSI acquisition in analog function computation networks,” IEEE Trans. Veh. Tech., vol. 68, no. 3, pp. 2361–2373, 2019.
  * [9] Y. Shao, D. Gunduz, and S. C. Liew, “Bayesian over-the-air computation,” arXiv preprint arXiv:2109.03780, 2021.
  * [10] A. Mahmood et al., “Clock synchronization over IEEE 802.11—a survey of methodologies and protocols,” IEEE Trans. on Industrial Informatics, vol. 13, no. 2, pp. 907–922, 2016.
  * [11] M. Goldenbaum and S. Stanczak, “Robust analog function computation via wireless multiple-access channels,” IEEE Trans. on Commun., vol. 61, no. 9, pp. 3863–3877, 2013.
  * [12] X. Zhao et al., “Broadband digital over-the-air computation for asynchronous federated edge learning,” in Proc. IEEE International Conf. Commun., pp. 5359–5364, 2022.
  * [13] M. I. Jordan et al., “Communication-efficient distributed statistical inference,” Journal of the American Statistical Association, 2018.
  * [14] J. Konečnỳ et al., “Federated learning: Strategies for improving communication efficiency,” arXiv preprint arXiv:1610.05492, 2016\.
  * [15] R. Tandon, Q. Lei, A. G. Dimakis, and N. Karampatziakis, “Gradient coding: Avoiding stragglers in distributed learning,” in Int. Conf. Mach. Learn., pp. 3368–3376, PMLR, 2017.
  * [16] T. Chen et al., “LAG: Lazily aggregated gradient for communication-efficient distributed learning,” Advances in Neural Information Processing Systems, vol. 31, 2018.
  * [17] G. Zhu et al., “One-bit over-the-air aggregation for communication-efficient federated edge learning: Design and convergence analysis,” IEEE Trans. Wireless Commun., vol. 20, no. 3, pp. 2120–2135, 2020.
  * [18] S. Bayat and S. Daei, “Separating radar signals from impulsive noise using atomic norm minimization,” IEEE Trans. Circuits Syst., II, Exp. Briefs, vol. 68, no. 6, pp. 2212–2216, 2020.
  * [19] V. Chandrasekaran et al., “The convex geometry of linear inverse problems,” Foundations of Computational mathematics, vol. 12, no. 6, pp. 805–849, 2012.
  * [20] E. J. Candès and C. Fernandez-Granda, “Towards a mathematical theory of super-resolution,” Communications on Pure and Applied Mathematics, vol. 67, no. 6, pp. 906–956, 2014.
  * [21] T. T. Georgiou, “The Carathéodory–Fejér–Pisarenko decomposition and its multivariable counterpart,” IEEE Trans. Automatic Control, vol. 52, no. 2, pp. 212–228, 2007.
  * [22] S. Boyd, S. P. Boyd, and L. Vandenberghe, Convex optimization.  Cambridge university press, 2004.
  * [23] Y. Hua and T. K. Sarkar, “Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise,” IEEE Trans. Acoust., Speech, Signal Process, vol. 38, no. 5, pp. 814–824, 1990.
  * [24] V. I. Morgenshtern and E. J. Candes, “Super-resolution of positive sources: The discrete setup,” SIAM Journal on Imaging Sciences, vol. 9, no. 1, pp. 412–444, 2016.
  * [25] Y. LeCun et al., “Gradient-based learning applied to document recognition,” Proceedings of the IEEE, vol. 86, no. 11, pp. 2278–2324, 1998.