Starspots and Undetected Binary Stars Have Distinct Signatures in Young Stellar Associations

Young stars form in associations, meaning that young stellar associations provide an ideal environment to measure the age of a nominally coeval population. Isochrone fitting, which is the typical method for measuring the age of a coeval population, can be impacted by observational biases that obscure the physical properties of a population. One feature in isochrone fits of star-forming regions is an apparent mass-dependent age gradient, where lower-mass stars appear systematically younger than higher-mass stars. Starspots and stellar multiplicity are proposed mechanisms for producing the mass-dependent age gradient, but the relative importance of starspots versus multiplicity remains unclear. We performed a synthetic red-optical low-resolution spectroscopic survey of a simulated analog to a 10 Myr stellar association including mass-dependent multiplicity statistics and age-dependent starspot coverage fractions. We found that undetected starspots alone do not produce an apparent mass-dependent age gradient, but instead uniformly reduce the average measured age of the population. We also found that binaries continue to produce an apparent mass-dependent age gradient, and introduce more scatter in the age measurement than spots, but are easily removed from the population as long as there are good distance measurements to each target. We conclude that it is crucial to incorporate treatments of both starspots and undetected stellar multiplicity into isochrone fits of young stellar associations to attain reliable ages.


INTRODUCTION
Stars typically form in associations within giant molecular clouds (Lada & Lada 2003).Determining the ages of these nominally coeval young stellar associations is important for a variety of scientific goals.For example, ages are used to calibrate evolutionary models (e.g., MIST; Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015;;Choi et al. 2016;Dotter 2016), determine the duration and detailed history of star formation (e.g., Kerr et al. 2021), reassemble star-forming events after the fact (e.g., Zucker et al. 2022), and measure the protoplanetary disk dissipation timescale, which sets the timescale for exoplanet formation (e.g., Mamajek 2009).One of the most efficient techniques for measuring stellar ages is isochrone fitting, which compares sources on an HR diagram to evolutionary model tracks at various ages.Because young stars are still descending the Hayashi track (Hayashi 1961), Corresponding author: Kendall Sullivan ksulliv4@ucsc.edu it is possible to place them on isochrones to determine their ages with high relative precision.
One key site of local star formation is the Scorpius-Centaurus (Sco-Cen) OB association, which is the closest collection of massive stars and thus also the one of the largest nearby associations of roughly coeval low-mass stars.Sco-Cen is comprised of three main subgroups: Upper Scorpius (or Upper Sco), Lower Centaurus-Crux, and Upper Centaurus-Lupus.Upper Sco is the youngest region, although there is substantial substructure with age spreads in each of the regions (e.g., Kerr et al. 2021;Žerjal et al. 2023).Sco-Cen hosts ∼ 10 4 low-mass stars (Rizzuto et al. 2015), and falls at a distance of ∼ 140 pc (e.g., Rizzuto et al. 2015;Žerjal et al. 2023), making it an excellent location for studying young stellar populations, reconstructing the star formation history of stellar associations, calibrating evolutionary models, and measuring the circumstellar disk dissipation timescale.
However, one major barrier in using young stellar associations to study star and planet formation is difficulty in measuring their ages.In Upper Sco, early work (e.g., Preibisch & Zinnecker 1999;Preibisch et al. 2001Preibisch et al. , 2002) ) suggested a 5 Myr age using primarily low-mass stars, while more recent work has supported a 10 Myr age using higher-mass stars (e.g., Pecaut et al. 2012;Feiden 2016;Pecaut & Mamajek 2016;Sullivan & Kraus 2021).Some work has observed an apparent mass-dependent age gradient in Upper Sco, where low-mass stars almost always appear younger than high-mass stars even in a uniformly observed and characterized population (Pecaut et al. 2012;Rizzuto et al. 2015;Pecaut & Mamajek 2016), suggesting that disagreement between different analyses could be caused by a systematic feature of the population.If the observed mass-dependent age gradient is physical, it could indicate hierarchical star formation (with massive stars forming after low-mass stars; Elmegreen & Lada 1977;Preibisch & Zinnecker 1999;Krumholz et al. 2019).However, most work has focused on potential observational biases or astrophysical features of young stellar populations that could produce an apparent mass-dependent age gradient in an approximately coeval population.Some of the possible observational or astrophysical sources of this systematic feature in associations are starspots, undetected binary stars, uncertain distances, uncertain extinctions, and multimodal star formation.Each of these phenomena has an impact on the inferred HR diagram, and can alter age measurements.
Much of the focus on systematic effects that can alter inferred isochronal ages has fallen on magnetic activity, especially starspots (Somers & Pinsonneault 2015;Feiden 2016;Gully-Santiago et al. 2017;Somers et al. 2020;Cao et al. 2022;Pérez Paolino et al. 2023, 2024).Evolutionary models incorporating starspots can more accurately reproduce observed radii of young eclipsing binaries (e.g., Tofflemire et al. 2023), and young stars are known to be heavily spotted (e.g., Gully-Santiago et al. 2017;Morris 2020;Cao & Pinsonneault 2022;Pérez Paolino et al. 2023).The presence of starspots introduces a secondary temperature component into the stellar spectrum, which can cause the inferred stellar temperature to be cooler than the photospheric temperature (e.g.Gully-Santiago et al. 2017;Pérez Paolino et al. 2024).Because starspot coverage fraction is dependent on age and magnetic activity level (e.g., Morris 2020;Somers et al. 2020), undetected starspots could produce an apparent mass-dependent age gradient that may be similar to the gradient observed in Upper Sco.Magnetic activity can also impact isochrone fits by altering the interior structure of the star (Feiden 2016), and the Feiden (2016) analysis of Upper Sco that incorporated magnetic fields (but not starspots) found a consistent 10 Myr age for Upper Sco across all stellar masses.The magnitude of the effect of starspots on observations is not established, and it is unclear if starspots alone sufficiently explain the observed mass-dependent age gradient in Upper Sco.An alternative explanation for the observed massdependent age gradient of Upper Sco is undetected binary stars.Binaries can produce a similar observational effect to starspots, with the presence of a secondary star introducing a cooler component into the spectrum.Sullivan & Kraus (2021) simulated an Upper Sco analog without starspots but including undetected binary stars, and found that undetected binaries in a single-age 10 Myr population caused low-mass stars to appear several Myr younger than high-mass stars.They also found that undetected binaries introduce significant scatter into the HR diagram that increases toward lower masses, which is expected given mass-dependent binary statistics.
One key distinguishing feature between starspots and undetected binaries in an HR diagram is that starspots reduce the luminosity of a system while decreasing the apparent temperature, so systems with spots will appear fainter and cooler, while binaries increase the luminosity of a system while decreasing the apparent temperature, so systems with binaries will appear brighter and cooler.Stellar evolution leads to stars becoming less luminous at the same T eff as they age, meaning that a brighter, cooler star will appear younger than a fainter, cooler star, because the spotted star will more closely follow its original isochrone (see Figure 5 for an illustration of this effect).However, distinguishing between these scenarios requires an accurate distance to measure a precise luminosity, which was difficult for analyses before the Gaia mission (Gaia Collaboration et al. 2016, 2018, 2021, 2022), or for extincted, optically faint populations.Another complicating factor is the possible multimodal star-formation history of Upper Sco (e.g., Žerjal et al. 2023), but a multi-age population should only introduce additional scatter on top of that caused by undetected binaries and starspots, rather than introducing additional features like the mass-dependent age gradient.
To expand on the work of Sullivan & Kraus (2021) and explore the relative impact of starspots versus undetected binary stars on isochronal ages, we simulated a spectroscopic survey of an analog of the Upper Scorpius star-forming region including both undetected binary stars and starspots using the SPOTS evolutionary tracks (Somers et al. 2020).In this paper we present the results from simulated age and mass measurements of spotted single stars at an age of 10 Myr, as well as results from a similar simulation that incorporated both spots and unresolved binaries.

SIMULATING A SPECTROSCOPIC SURVEY OF UPPER SCORPIUS
Our original simulated survey of an analog of Upper Scorpius including unresolved binaries is described in Sullivan & Kraus (2021).In the following subsections we provide a brief summary of our original methodology and detail the changes made to the simulation for this work.

Population Statistics
We assigned properties to our systems at the population level using an initial mass function (IMF), observed binary star properties, and estimated starspot properties.For the desired number of systems, we drew stellar masses using a Chabrier IMF (Chabrier 2003(Chabrier , 2005)), which is defined as a Salpeter power law (Salpeter 1955) at the high-mass end and a lognormal on the low-mass end: where we used the parameters from Chabrier (2005) and thus defined m c = 0.2M ⊙ , σ = 0.55M ⊙ , M cutoff = 1M ⊙ , and α = 1.3.We artificially reduced the number of low-mass stars at an arbitrary second cutoff mass to produce more high-mass stars with a smaller population size, but performed that reduction after creating the initial IMF.A typical second cutoff mass was 0.8M ⊙ , where stars below that mass were suppressed by a factor of 25.This reduction did not change our analysis or conclusions, and only served to increase our computational efficiency.
After drawing a primary star mass, we assigned an age to each system by randomly drawing a value from a normal distribution with µ = 10Myr and σ = 2Myr.These ages matched the age distribution in Sullivan & Kraus (2021) and were consistent with the age derived using F stars in Pecaut et al. (2012).We also assigned a distance for each system, which we drew from a normal distribution with µ = 140pc and σ = 20pc, based on results from Rizzuto et al. (2015).With the assigned distance, mass, and age of each primary star, we added starspots and stellar multiplicity to create a complete system.Figure 1 shows various distributions from the simulated primary star population.
To add stellar multiplicity, we assigned a system a secondary star based on its mass using a piecewise mass-dependent multiplicity fraction function of f = {0.25,0.35, 0.45, 0.7} with transitions between values at transition masses m = {0.15,0.4, 1.2}M ⊙ .These fractions were based on observations from the literature (e.g., Kraus & Hillenbrand 2008;Raghavan et al. 2010;Duchêne & Kraus 2013;De Rosa et al. 2014;Ward-Duong et al. 2015;Winters et al. 2019;Tokovinin & Briceño 2020).If a system was a multiple, we assigned it a separation and mass ratio that were again based in observations and mass-dependent, using functions derived in Sullivan & Kraus (2021).We assigned a resolution limit of 1", and discarded any secondary star that fell outside that separation limit, since our goal was to investigate unresolved and undetected binary stars.An arcsecond is approximately the spatial resolution of the Gaia mission (Gaia Collaboration et al. 2016) and is also a typical value for good seeing at most ground-based observatories, meaning that many binaries closer than 1" would not be detected without dedicated surveys to detect multiplicity.The properties of the generated secondary star population are shown in Figure 2.
To add starspots, we imposed a normally distributed spot coverage fraction with a mean value of 50% and a standard deviation of 10%.For drawn values below 0% or above 80% we resampled the distribution to avoid nonphysical values or starspot covering fractions that were outside the range of the SPOTS models (Somers et al. 2020).Only a small fraction of systems (< 1%) needed to be redrawn because they fell outside the limits.

Data Generation
Using the statistics described above, we produced systems with the following known properties: mass, age, distance, and starspot coverage fraction, as well as mass ratio, separation, and a second starspot coverage fraction if the system was a binary.To produce simulated data using model spectra, we had to convert these physical properties to observables using evolutionary models.We used the SPOTS models (Somers et al. 2020), which incorporate starspots and provide both starspot and photospheric temperatures.The SPOTS models have a fixed range of possible spot coverage fractions, ranging from 0% to 85%, and a limited mass range of 0.1-1.3M⊙ .Therefore, for this analysis, we restricted the simulated population to primary star masses 0.16 ≤ M p ≤ 1.3 M ⊙ , and we redrew any secondary stars with masses that fell below 0.1 M ⊙ .Since the mass ratio distribution of lowmass binaries becomes increasingly biased toward unity (e.g., Winters et al. 2019), the number of secondaries that needed to be redrawn was small.
Using the SPOTS models, we linearly interpolated the models in spot coverage fraction, mass, and age, and calculated a photospheric temperature, spot temperature, luminosity, surface gravity, and VRIJHK magnitudes for each component of the system.We used ) for the primary and secondary.We created a spectrum for each component of each photosphere (i.e., model spectra at T phot and T spot for singles) by bilinearly interpolating each spectrum pixel-by-pixel in T eff and log(g) using the Husser et al. (2013) model spectra calculated with the PHOENIX code.We chose these model spectra for their wide range of T eff and surface gravity.We normalized each component to the appropriate surface area (e.g., ), where f spot,1 is the spot coverage fraction), then summed the two (for singles) or four (for binaries) spectra (e.g., Figure 3) to produce a single composite spectrum for each system.To match previous moderate-resolution studies of Upper Sco, we created the spectra at a resolution of R∼2100 by convolving with a Gaussian of the appropriate width, and covered the red-optical wavelength range from 6000-9000 Å (e.g.Preibisch et al. 2002;Pecaut et al. 2012;Rizzuto et al. 2015;Pecaut & Mamajek 2016).
We also analyzed SEDs for each system created by synthetic broadband photometry.The SPOTS models already incorporate the impact of spots into the model photometry, so we simply added the distance modu-lus and combined magnitudes for the two components of the binaries by summing the fluxes.We assumed a distance measurement with a typical parallax error of 0.05 mas (consistent with Gaia parallax measurements for relatively bright sources; Gaia Collaboration et al. 2021), so we inverted each assigned system distance, perturbed it with a Gaussian with σ = 0.05 mas, then inverted the perturbed parallax to achieve a perturbed distance "measurement".Finally, we extincted each component with a randomly assigned extinction ranging from 0 ≤ A V ≤ 2 drawn from a uniform distribution, and converted between wavelengths using the extinction law from Cardelli et al. (1989).At the conclusion of the data generation process we had a single spectrum and SED in VRIJHK for each system, incorporating starspots and undetected binaries with separations ρ < 1".

System Analysis
To analyze each composite spectrum, we used a Markov Chain Monte Carlo (MCMC) method implemented as a modified Gibbs sampler to retrieve a temperature for each system, then used the VRIJHK SED, best-fit T eff and A V and assigned distance to measure the luminosity.Using the best-fit luminosity and temperature, we inferred a mass and age. .Composite and component spectra for a typical binary system.The composite spectrum (black), primary photosphere (blue), primary spot (orange), secondary photosphere (yellow) and secondary spot (red) spectra are all shown.The primary flux dominates at the short wavelengths, but the cooler components progressively contribute more flux toward longer wavelengths.
To fit each system, we initialized a walker with a guess T eff , surface gravity, and A V defined by the initial (primary star) values.To avoid beginning in an exactly correct location for the single stars, we perturbed the initial T eff and A V using a normal distribution with σ T eff = 100 K and σ A V = 0.1 mag.We also fit for a normalization factor, which we initially perturbed by σ norm = 1%, to account for the extra flux introduced by the secondary star.The normalization factor is required because we flux-normalized the "observed" spectra, rather than flux-calibrating the model spectra.We found that at the relatively low spectral resolution of our data our fit was not very sensitive to surface gravity changes, so we did not perturb the initial log g guess, and we did not fit for a new surface gravity value.Our fit results were not sensitive to the choice of initialization values, as discussed in Sullivan & Kraus (2021).To perform the modified Gibbs sampling, we took the perturbed guess values, produced a single-star spectrum using those values, and compared the model to the data using a χ 2 value.If the χ 2 was lower than the previous best value, we always accepted the guess, and if the χ 2 was higher we always rejected the guess.In contrast, a traditional Gibbs sampler occasionally accepts worse guesses than the current recorded best fit.
We fit in two stages, each of which had a coarsely sampled optimization fit followed by a fit that was more finely sampled to achieve the best-fit parameter values.In the initial coarse fit we drew random guesses from relatively broad normal distributions (σ T eff = 100 K, σ A V = 0.05 mag, σ norm = 1%) until the fitter ran for half the requested number of steps without finding a better fit.We then switched to finer sampling, using narrow normal distributions (σ T eff = 5 K, σ A V = 0.005 mag, σ norm = 0.05%) until the fitter ran for the other half of the requested number of steps without finding a better fit.We implemented an initial "burn-in" fit of 40 steps (i.e., a fit that reached 40 steps with one guess without achieving a better fit), then a full fit that ran for 200 steps.The typical reduced-χ 2 values after the second fitting stage were χ 2 << 1 for single stars without spots, which was expected given that we did not introduce additional artificial random noise, meaning that in the ideal case (i.e., an infinite run time) the reduced-χ 2 for a single star trend toward zero.We tested adding artificial noise to more closely match real observations, but found that it did not change our results and substantially slowed down our fits, so we chose to run our analysis in the optimal zero-noise scenario.
Our fitting method retrieved a best-fit composite T eff , A V , and normalization value for each system.Validation of this technique on single stars and tests of the extinction fitting are described in Sullivan & Kraus (2021).In general, for single stars, we were able to retrieve temperatures within a few K, and extinctions within 0.1 mag, although the fits became worse at higher temperatures (above ∼ 4500 K) because of a reduction in the number of spectral features present in our wavelength range at low spectral resolution.However, because our young stars were almost all relatively low-mass, almost all of our stars were expected to have precise fits in the spot-free single star scenario.
We also fit the VRIJHK SED to measure a composite luminosity for each system using the MESA Isochrones and Stellar Tracks (MIST) stellar evolutionary models (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015;;Choi et al. 2016;Dotter 2016).Using the best-fit temperature, we interpolated the luminosity at that temperature as a function of each point in the SED, then took the median of the interpolated luminosities as the best-fit value.Using the luminosity and temperature, we also measured an age and a mass for each composite system.Because the MIST evolutionary models do not go below 3000 K, we removed any systems that had a best-fit T eff cooler than 3000 K from the final results.This step typically eliminated fewer than 100 systems from the total sample of ∼ 10 4 systems, a decrease of 1% of our total sample size at most.

SIMULATION RESULTS
We simulated two populations, with one run composed of only single stars with starspots and the other run composed of single and undetected binary stars, all with starspots.Both runs can be compared against Sullivan & Kraus (2021), with simulations for single stars and binaries without spots.The single-star-only simulation provided a benchmark for the expected outcome if spots were the dominant factor in producing scatter or apparent mass-dependent age gradients in a single-age population, while the simulation including single and binary stars allowed a comparison between the impact of undetected multiplicity and the impact of starspots.Each run consisted of 5000 simulated systems with primary star masses ranging from 0.16 M ⊙ to 1.3M ⊙ .

Single Star Simulation
To demonstrate the impact of starspots without the presence of binary stars, we simulated a population of single stars with spots.Figure 4 shows the T eff and A V residuals for single stars with spots after fitting the single-only, spotted stars using a single-T eff spectrum.The left panel includes the T eff residual plotted against the input photospheric temperature.Nearly all of the systems at all temperatures have a derived T eff that is cooler than the photospheric T eff , demonstrating the well-established impact of starspots on inferred stellar T eff (e.g., Gully-Santiago et al. 2017;Cao et al. 2022).At lower T eff (T eff < 4300 K), there is a systematic residual; the average residual is 154 (64) K.At slightly higher T eff (T eff > 4300 K), there is an even larger systematic effect of spots, with an average residual of 326 (rms = 152) K.The scatter at low T eff is driven by the T eff bias produced by large starspot fractional coverage, while the scatter at high T eff is a result of spots, in addition to fewer T effsensitive spectral features being present at higher T eff at our chosen spectral resolution and wavelength range, as discussed in Sullivan & Kraus (2021).
In contrast to the T eff measurement, there is no systematic residual with temperature in A V , with residuals of -0.01 (rms = 0.1) mag for T eff < 4300 K and -0.06 (rms = 0.1) mag for T eff >4300 K.This excellent A V recovery remains consistent with the results from Sullivan & Kraus (2021), which found that A V was recovered accurately for both single stars and undetected binary stars in systems without spots.The A V residual is correlated with the spot coverage fraction, with a higher spot coverage fraction being correlated with an underestimated observed A V value.The reason for this underestimation is unclear, but is likely related to the presence of starspots, given that the same feature does not occur in the spot-free single star simulations in Sullivan & Kraus (2021).The apparent gaps in the residual are caused by small gaps in the mass function produced by the artificial reduction in the number of low-mass stars, and are not a physical feature of the A V residual structure.To demonstrate the observational impact of starspots on derived ages, Figure 5 shows an HR diagram of our input and output simulated populations with MIST isochrones underlaid.The input population (in orange) falls along the correct isochrone with no apparent gradient, while the output population appears younger and shows more scatter than the single population.

Binary Star Simulation
To explore the relative impact of undetected binaries and undetected starspots, we simulated a population that was identical to the single-star-only simulation except for the inclusion of a population of undetected binary stars with separations ρ < 1 ′′ .
Figure 6 shows the T eff and A V residuals plotted against the input T eff and A V values.The T eff residuals for single stars match those in the single-star-only simulation, which is expected because we did not change any of the population parameters for single stars.The binary stars systematically have larger T eff residuals than the single stars, with moderate-contrast binaries having the largest T eff residuals.This is because moderate-T effcontrast systems have secondaries with sufficient additional flux and spectral morphology changes to bias the T eff measurement.Similar to the single stars, the binary bias decreases toward higher masses as the spot coverage fraction decreases, but then increases again as the scatter grows at earlier spectral types.Regardless of the input primary star photospheric T eff , the binaries consistently have larger T eff residuals than the single stars.
In contrast to the T eff measurements, the A V residuals shown in Figure 6 are uniformly distributed, although they do show a slight bias that is dependent on the T eff contrast between the components of the binary, unlike the results from Sullivan & Kraus (2021).The impact on the A V measurement is not large (a maximum of ∼0.2 mag deviation from the true A V value), so we do  4, but for the simulation including unresolved binary stars (cross symbols, color-coded by T eff difference between primary and secondary star).Left: The single-star T eff residual is identical to the single-star-only run, as expected, while the systems hosting unresolved binary stars typically have T eff measurements that are biased toward even cooler temperatures, with typical residuals of ∼150 K for the cooler stars.Right: The AV measurement is not biased by the presence of binaries, demonstrated by the residuals that match those of the single-star-only case.However, for binaries, the residual is dependent on the T eff contrast between the two stellar components, likely because of the extra apparent reddening of the SED caused by the secondary star.Figure 7.The same as Figure 5, but for the simulation including binary stars (shown as cross markers, with color coding indicating the T eff contrast between the primary and secondary star).The binary stars systematically fall at younger ages and higher luminosities than either the single star input or output populations, even with spots included.The binaries have a large amount of scatter and another apparent mass-dependent age gradient, both produced by the massdependent population statistics of binary stars.not expect it to significantly impact our analysis or conclusions.
To compare to a typical HR diagram of a star-forming region, we plotted the simulated population on an HR diagram along with MIST isochrones, shown in Figure 7.The binary stars systematically fall on younger isochrones and with larger luminosities than the single star input or output populations.The binaries have a large amount of scatter and show a mass-dependent age gradient matching that seen in Sullivan & Kraus (2021), where systems with hotter primary stars appear older than systems with cooler primaries.With the improved distance precision from Gaia a binary sequence is apparent, with most binaries falling above the sequence of single stars.This is the result of spots causing single stars to appear fainter as well as cooler, while the binaries appear brighter and cooler, creating a more distinct separation than for systems without spots.

COMPARING AGE MEASUREMENTS BETWEEN SIMULATED POPULATIONS
For each simulation, we inferred an age for each star using the derived T eff and luminosity to compare to evolutionary models.In this subsection we present the inferred ages for a series of simulations.Figure 8 shows summary histograms for the input and output age distributions of our simulations.We present results from four different simulated scenarios.
Simulation A (top panel of Figure 8) presents a singlestar population with spots, but with the bilinearlyinterpolated f spot perturbed by a Gaussian with a fractional width of 20%.The average K and M ages are both 7.1 (rms ∼ 2) Myr.Therefore, the spot-only scenario leads to ages that are uniformly younger by ∼ 30%, but with age spreads that match the input 2 Myr spread.Thus, starspots alone do not produce an apparent massdependent age gradient.The ∼ 30% reduction in age is smaller than the ∼ 50% reduction in age found by Pérez Paolino et al. ( 2024 (2024) would thus be more heavily spotted than our simulated Upper Sco population.A larger spot coverage fraction would likely lead to an even larger effect on apparent age for the younger stars.
Simulation B (lower left panel of Figure 8) shows the population of single stars and binaries, with perturbed f spot .The binaries appear much younger than the rest of the population, producing a second peak of very young stars.Because the impact of multiplicity is more distinctive for lower-mass stars, the M star population contributes almost all the systems that comprise the young peak, while the K star population appears more spread out with a smaller young peak.The mean ages and spreads reflect the presence of the second peak that is caused by binaries, with the average M star age τ M = 5.7 (rms = 3.3) Myr, and the average K star age τ K = 6.4 (rms = 3.2) Myr.
It is likely that in a real set of observations, the secondary peak caused by the presence of binaries would be identified as a "binary main sequence" and removed from the analysis.Therefore, we also wished to examine the population with binaries included but removing the easily-detectable binary systems (i.e., the young peak of the retrieved-age histogram).The lower right panel of Figure 8 shows Simulation B, but with the age distributions truncated at the local minimum for the age distribution of the M stars at ∼3 Myr (Simulation C hereafter).In this scenario, mimicking removal of the obvious binaries, the average M star age is τ M = 7.2 (rms = 2.6) Myr, and the average K star age τ K = 7.1 (rms = 2.8) Myr.
Because adding scatter to the starspot coverage fraction does not introduce additional scatter into the age measurements (i.e., the rms scatter of Simulation A does not change from the initial input scatter of 2 Myr), any additional scatter beyond the intrinsic 2 Myr scatter is caused by the presence of additional undetected binaries present even in the truncated sample.Although the rms scatter is increased in Simulation C relative to Simulation A, the mean ages for both K and M stars are identical (within error) to Simulation A, which did not include binaries.This suggests that if the majority of binaries can be removed, age measurements using isochrone fitting should retrieve mean or median ages that are predominantly biased by the presence of undetected starspots, rather than binaries, but the age spread will be impacted by undetected multiplicity.
In our previous paper exploring the impact of binary stars on age measurements (Sullivan & Kraus 2021), we did not truncate the age distribution at the location of any apparent bimodality in the age distribution, although we did observe a double-peaked age distribution for the M stars.Sullivan & Kraus (2021) found that in the exact-distance scenario (which we essentially replicated in this work by mimicking Gaia observations), the average M star age was τM,2021 = 7.0 (rms = 4.3) Myr.The simulation in this paper that most closely matches the analysis of Sullivan & Kraus (2021) is Simulation C, where the average age for M stars was τM,2023 = 7.2 (rms = 2.6) Myr.In Sullivan & Kraus (2021) the rms is increased relative to Simulation C. In this work, the starspots drive the single-star age distribution to younger ages, decreasing the total rms scatter and also decreasing the average age of the population in conjunction with the young-star peak caused by the binaries.
The starspots in Simulation B and C are also the cause of the more defined binary sequence in this work relative to Sullivan & Kraus (2021).Multiplicity decreases measured T eff but increases system luminosity, while starspots decrease both measured T eff and luminosity.The combination of these effects causes more separation between the single and binary stars in the simulations of this work than in Sullivan & Kraus (2021).This suggests that it may be easier to remove binaries from low-mass stellar populations dominated by M stars than to remove them from even slightly higher-mass samples such as those dominated by K stars.Also, high-mass stellar populations are more likely than low-mass stellar populations to be hosting binaries that are challenging to detect via isochrone fitting, meaning that it will be more difficult to remove binaries from G or K star samples than from M star samples.
Additional challenges to age measurements will include uncertain distances, which will blur out the bimodality that enables removal of binaries for low-mass stars.Gaia astrometric solutions become worse for young stars because of disks and variability (e.g., Fitton et al. 2022), extinction, binarity itself (Wood et al. 2021), and more distant stellar populations, meaning that typical parallax errors for young stellar associations may be larger than what we assume in this work.

CONCLUSIONS
To explore the relative impact of starspots and undetected binary stars on age measurements of young stellar associations, we simulated two spectroscopic surveys of an analog of the Upper Sco association.One simulated survey consisted of only single stars with starspots that were not corrected for, while the other survey was identical but also included undetected binary stars.We found that starspots do not produce an apparent massdependent age gradient, but a uniform spot coverage fraction of 50% leads to an age for both M and K stars that is ∼ 30% younger than the true age, with a slightly larger age spread, assuming the obvious binaries are removed.If the binary sequence is not removed (as in situations with less certain distances), the lower-mass stars appear systematically younger than the highermass stars, and both populations appear younger than their true ages.Thus, we found that starspots produce a statistically significant age bias, in keeping with other studies, but found that they cannot replicate the apparent mass-dependent age gradient observed in starforming regions.Ideally, any census of a young stellar population would identify binaries and either analyze them separately or fully exclude them, but if that is not possible a forward-modeling routine that includes the impacts of both binaries and spots will be important for accurately retrieving mass functions and star formation histories.

Figure 1 .
Figure 1.Histograms for the simulated input population of single stars.Top row: stellar mass, age, and extinction.Bottom row: stellar T eff , distance, and spot coverage fraction.The apparent bimodality in mass, T eff , and spot coverage fraction is a result of the suppression of low-mass stars below 0.8M⊙ to improve computational efficiency.thedefault SPOTS spot temperature contrast value of T spot /T phot = 0.8.Using the luminosity and temperature, we inferred a stellar radius (R ⋆ =

Figure 2 .Primary
Figure 2. Histograms for the simulated input population of secondary stars.Top row: secondary star mass, separation, and eccentricity.Bottom row: period, secondary star T eff , and starspot coverage fraction.

Figure 4 .Figure 5 .
Figure4.Fundamental retrieved parameters for single stars with spots.Left: T eff residual plotted against the input (true) primary star T eff .Cooler stars (M and late K) have larger systematic T eff -dependent offsets because of spots than hotter stars (early K and G).The average T eff residual is 154 (rms = 64) K (T eff < 4300 K) and 326 (rms = 152) K (T eff > 4300 K).The increased scatter in hotter stars results from a reduction in the number of T eff -sensitive spectral features at hotter temperatures with our spectral resolution and wavelength range.The measured T eff for stars with spots are typically systematically cooler than the true photospheric T eff .Right: AV residual versus input AV for all single stars with spots.The presence of spots does not bias the AV retrieval, with residuals of -0.01 (rms = 0.1) mag (T eff < 4300 K) and -0.06 (rms = 0.1) mag (T eff > 4300 K).The AV residual is correlated with fspot because of the reddening introduced by the presence of uncorrected spots.The gap in the residual around 0.03 is caused by the artificial reduction of low-mass stars, and is not reflective of any bias in the AV recovery.

Figure 6 .
Figure6.The same as Figure4, but for the simulation including unresolved binary stars (cross symbols, color-coded by T eff difference between primary and secondary star).Left: The single-star T eff residual is identical to the single-star-only run, as expected, while the systems hosting unresolved binary stars typically have T eff measurements that are biased toward even cooler temperatures, with typical residuals of ∼150 K for the cooler stars.Right: The AV measurement is not biased by the presence of binaries, demonstrated by the residuals that match those of the single-star-only case.However, for binaries, the residual is dependent on the T eff contrast between the two stellar components, likely because of the extra apparent reddening of the SED caused by the secondary star.

Figure 8 .
Figure8.Histograms of input (dark blue) and output (light blue) ages for different simulations.Subsamples of K stars and M stars are plotted as the red and black histograms, respectively.The black (red) solid and dashed lines show the M (K) star mean and median age.(Top:) Single stars with 50% mean fspot with 20% intrinsic Gaussian scatter.The mean age and rms scatter for K and M stars are both 7.1 Myr (rms ∼ 2) Myr.(Bottom left:) A simulation with binary stars and spots, with 10% Gaussian scatter added to fspot.The average K star age is 6.4 (rms = 3.2) Myr, and the average M star age is 5.8 (rms = 3.3) Myr.The young peak in the output histogram is caused by undetected binary stars and produces a large age spread.(Bottom right:) The same as (b) but with the output age distribution truncated at the local minimum (∼ 3 Myr) between the young peak and the rest of the population.This population with apparent binaries removed has a mean age of 7.1 (rms = 2.8) Myr for K stars, and 7.2 (rms = 2.6) Myr for M stars.The presence of binaries increases the apparent age spread.This is likely because the spot coverage fraction is agedependent, and the Taurus targets of PérezPaolino et al. (2024) would thus be more heavily spotted than our simulated Upper Sco population.A larger spot coverage fraction would likely lead to an even larger effect on apparent age for the younger stars.Simulation B (lower left panel of Figure8) shows the population of single stars and binaries, with perturbed f spot .The binaries appear much younger than the rest of