Evaporation Ages: a New Dating Method for Young Star Clusters

The ages of young star clusters are fundamental clocks to constrain the formation and evolution of pre-main-sequence stars and their protoplanetary disks and exoplanets. However, dating methods for very young clusters often disagree, casting doubts on the accuracy of the derived ages. We propose a new method to derive the kinematic age of star clusters based on the evaporation ages of their stars. The method is validated and calibrated using hundreds of clusters identified in a supernova-driven simulation of the interstellar medium forming stars for approximately 40 Myr within a 250 pc region. We demonstrate that the clusters' evaporation-age uncertainty can be as small as about 10% for clusters with a large enough number of evaporated stars and small but realistic observational errors. We have obtained evaporation ages for a pilot sample of 10 clusters, finding a good agreement with their published isochronal ages. The evaporation ages will provide important constraints for modeling the pre-main-sequence evolution of low-mass stars, as well as to investigate the star-formation and gas-evaporation history of young clusters. These ages can be more accurate than isochronal ages for very young clusters, for which observations and models are more uncertain.


Introduction
The precise determination of stellar ages is a fundamental tool in cosmology and astrophysics.The ages of the oldest globular clusters provide a lower limit to the age of the universe (e.g., Chaboyer et al. 1996;Krauss & Chaboyer 2003;Jimenez et al. 2019), while at the other extreme, the ages of the youngest stars are fundamental clocks to time the evolution of protoplanetary disks and to guide our understanding of the origin of planets (e.g., Williams & Cieza 2011;Bell et al. 2013;Berger et al. 2023).
An alternative approach to determining the age of unbound star clusters, independent of stellar-evolution models, is the analysis of their kinematics.The first attempts to derive kinematic expansion ages, based on the correlation of the positions and motions, date from several decades ago (Blaauw 1964).A more sophisticated technique consists of tracing the positions of stars, members of a cluster, back in time using their current velocities and a Galactic potential.The dynamical traceback age is then the time corresponding to the highest stellar density.This technique requires the 3D positions and velocities of individual stars, which has been an observational limitation for years (Brown et al. 1997;Ortega et al. 2002;de la Reza et al. 2006;Ducourant et al. 2014;Donaldson et al. 2016;Miret-Roig et al. 2018).However, thanks to the precise astrometry of Gaia (Gaia Collaboration 2023) and extensive complementary radial velocity surveys like the Apache Point Observatory Galactic Evolution Experiment, APOGEE (Majewski et al. 2017), it is now possible to obtain significantly more precise kinematic ages (Crundall et al. 2019;Miret-Roig et al. 2020, 2022;Kerr et al. 2022a,b;Galli et al. 2023;Couture et al. 2023).
Although kinematic ages are too uncertain for old clusters, they have become complementary to the isochrone fitting method (hereafter CMD) for ages <50 Myr.This is the age range for the clusters' sample selected in this study.Very young clusters (<5 Myr) may still be partly embedded in gas and dust from their parent cloud producing differential extinction on the cluster members, and they often display photometric variability widening the isochrone and the uncertainties on the age.CMD ages for the youngest clusters are also uncertain because the evolution of young accreting stars in the CMD is difficult to capture with standard evolutionary tracks where accretion is not appropriately modeled (Froebrich et al. 2006; Baraffe et al. 2009;Hosokawa et al. 2011;Jensen & Haugbølle 2018).
A detailed comparison of modern dynamical traceback and CMD ages, including observational uncertainties, has shown that there is a systematic difference between the two, with the former shorter than the latter (Miret-Roig et al. 2024).That study suggests that the two methods have a different "time zero": a star cluster may be gravitationally bound before the dispersion of the parent gas cloud, so the time zero of the kinematic method that measures the expansion time would start a few million years after the time zero of the CMD method.If the timescale of cluster formation and gas dispersal were known, it would be possible to correct the dynamical traceback age to a similar time zero as the CMD age.Conversely, the age difference between the two methods may give us important clues about the cluster formation and gas dispersal processes.
This work proposes a new method of determining kinematic ages: instead of measuring the "expansion" age of the whole cluster, we evaluate the "evaporation" ages of its individual stars, and estimate the cluster age from that of the oldest evaporation ages.We subsequently show that this method results in ages on the order of the actual stellar ages, and we calibrated it using star clusters from a numerical simulation.We also applied the method to observational data, finding kinematic ages in reasonable agreement with the CMD ages.
In the following, the term "star cluster" is used to refer to both bound and unbound groups of stars, irrespective of their stellar density or total mass.In Sect.2, we briefly describe the star formation simulation, and in Sect. 3 we characterize the clusters identified from the numerical data.The method is outlined in Sect. 4 and then we explain how it was validated and calibrated in Sect. 5 using the simulated clusters.In Sect.6 we explain how we applied the method to observations of real star clusters, and we finally summarize our conclusions in Sect.7.

Simulation
The kinematic-age method presented in this work is tested with star clusters from the same supernova (SN) driven magnetohydrodynamic (MHD) simulation as in Padoan et al. (2017Padoan et al. ( , 2020)).The reader is referred to those papers for details of the numerical methods.The 3D MHD equations are solved with the Ramses adaptive-mesh-refinement (AMR) code (Teyssier 2002(Teyssier , 2007;;Fromang et al. 2006) within a cubic region of size L box = 250 pc, total mass M box = 1.9 × 10 6 M , and periodic boundary conditions.The initial conditions are taken from a SNdriven simulation that was integrated for 45 Myr without selfgravity (Padoan et al. 2016b) with a mean density n H,0 = 5 cm −3 and a mean magnetic field B 0 = 4.6 µG.The root mean square (rms) magnetic field generated by the turbulence has a value of 7.2 µG and an average of |B| of 6.0 µG, consistent with the value of 6.0 ± 1.8 µG derived from the "Millennium Arecibo 21-cm Absorption-Line Survey" by Heiles & Troland (2005).
The only driving force is from SN feedback, with SNe determined by the position and age of the massive sink particles formed when self-gravity is included (the SNe are randomly generated prior to star formation).When gravity is introduced, starting at t = 55.5 Myr from the initial conditions, the minimum cell size is dx = 0.0076 pc, obtained through a uniform root-grid of 512 3 cells and six AMR levels.With this setup, the simulation is run for a period of approximately 40 Myr.To follow the collapse of prestellar cores, sink particles are created in cells where the gas density is larger than 10 6 cm −3 , according to several criteria designed to avoid creating spurious sink particles in regions where the gas is not collapsing (see Haugbølle et al. 2018, for details).When a sink particle of mass larger than 7.5 M has an age equal to the corresponding stellar lifetime for that mass (Schaller et al. 1992), a sphere of 10 51 erg of thermal energy is injected at the location of the sink particle to simulate the SN explosion Padoan et al. (2016b).
With a total mass of 1.9×10 6 M , the mean column density of the simulation is 30 M pc −2 , typical for spiral arms in the outer Galaxy (Heyer & Terebey 1998).In fact, a lower-resolution version of this simulation (Padoan et al. 2016b,a;Pan et al. 2016) has been shown to produce dense clouds with properties consistent with those of real molecular clouds in the 12 CO FCRAO Outer Galaxy Survey (Heyer et al. 1998(Heyer et al. , 2001)).These starforming clouds are formed ab initio in the simulation, as a result of the large-scale dynamics driven by SNe.Although their disruption is also caused by SNe alone (HII regions and stellar winds are not included in the simulation yet), the clouds in the simulation have realistic lifetimes (Lu et al. 2020), comparable to estimates from observations of nearby galaxies (e.g., Chevance et al. 2022;Lee et al. 2023).
The simulation has generated ∼7000 stars with mass >1.0 M and ∼600 stars with mass >8 M .The star formation is distributed over many different clouds with realistic values of the SFR, and the global SFR corresponds to a mean gas depletion time in the computational volume of almost 1 billion years, also realistic for a 250-pc scale (Padoan et al. 2017).Young stars are found inside the densest filaments, while older ones have already left their parent clouds.Most of the stars in the simulation are formed in clusters, some of which have cleared their surrounding gas thanks to SN explosions of their most massive members.

Identification of star clusters
For this work, we select six snapshots at time intervals of approximately 4.4 Myr, with the first one at 17.2 Myr after gravity is introduced in the simulation, and the final one at 39.2 Myr.This selection yields a total sample of 26 842 stars with mass >1.0 M .The stars in each snapshot are assigned into clusters using their 6D phase space information, using a Gaussian Mixing Model (GMM) from the Scikit-learn software (Pedregosa et al. 2011).We explored models with different numbers of Gaussians (between 1 and 70) and chose the one where the Bayesian information criterion (BIC) stabilized to a minimum value.This procedure results in a total of 339 star clusters (35, 42, 63, 67, 65, and 67 in each of the snapshots in chronological order).
The distributions of the overall properties of our sample of star clusters are shown in the histograms in Fig. 1, where the cluster age is defined as the median age of the stars and the radius as the median distance of the stars to the center of mass of the cluster.The age distribution has a median value of 13 Myr, is relatively flat in the range between approximately 1 and 20 Myr, and decreases between 20 and 30 Myr.The mass and radius distributions have median values of 217 M and 15.1 pc, and cover the approximate ranges from 20 to 2000 M and from 1 to 100 pc, respectively.The median rms velocity of the stars in the clusters is 2.9 km s −1 , ranging from 1 to 10 km s −1 .We compared our numbers to the Gaia DR3 all-sky cluster catalog made by Hunt & Reffert (2023), making a cut to high-reliability, young (<50 Myr), nearby (<500 pc) clusters.We found that the median age was comparable to ours at 17 Myr, while the radius was significantly smaller at 5 pc, which is not surprising as they have a large fraction of bound clusters.On the other hand, the clusters in our observational sample presented in Sect.6 have a median radius of 7.3 pc but reach as high as 24 pc, comparable to that of the simulated clusters.Castro-Ginard et al. (2022) found a median velocity dispersion of 2.3 km s −1 for their class A open clusters, which is comparable to our median velocity dispersion of 2.9 km s −1 .
The scatter plots in the panels of Fig. 1 show the square root of the ratio of kinetic and gravitational energies versus the age, mass, size, and rms velocity of the clusters.As can be seen from the E k /E g ratio and the sizes, the great majority of our systems are unbound clusters that have expanded to rather large sizes.However, we also have 17 bound clusters (5% of the total sample), and a few rather concentrated clusters with M ∼ 10 3 M , R ∼ 1 pc and E k /E g ∼ 1.As mentioned above, in this work we refer to both bound and unbound star groups as clusters, and our evaporation-age method works for both types of clusters.
Because of the limited spatial resolution, the stellar initial mass function (IMF) in the simulation is incomplete below a few solar masses, and only stars >1.0 M are used in this work.The lack of low-mass stars, as well as the incomplete binary statistics, may affect the early dynamical evolution of some of the densest clusters in the simulation, while its effect on most of the diffuse clusters is probably not significant.However, because our sample contains a large number of star clusters formed and evolved ab initio from a self-consistent simulation of a very large ISM region, it represents a significant improvement over previous ad hoc models of single star clusters from pure N-body simulations without gas-dynamics, or star formation simulations of very small regions of just a few pc that were used to test kinematicage methods (e.g., Crundall et al. 2019).Although this should be tested with future simulations including low-mass stars, we do not expect the calibration of the method to be sensitive to details of the evolution of the clusters.If anything, the inclusion of low-mass stars would make the method even more successful than concluded here, due to improved statistics.A significant fraction of the low-mass stars would undergo early evaporation, so their inclusion would increase the number of reliable evaporation ages, improving the accuracy of the method, particularly in the case of bound clusters.

Evaporation-age method
The traditional way to determine the dynamical age of a star cluster assumes that the stars are coeval and that the cluster was never gravitationally bound, or at least started to expand very soon after its formation.The ages can then be derived by adopting a realistic Galactic potential and tracing the stellar orbits back in time to find the time of maximum stellar density (de la Reza et al. 2006;Miret-Roig et al. 2020), or using forward modeling of the stellar orbits to avoid propagation of observational errors (Crundall et al. 2019).Besides their uncertainties due to the observational errors, these methodologies have important limitations related to the above assumptions: (i) if the cluster is initially bound for a significant time relative to its age, the kinematic age largely underestimates its real age; (ii) if the cluster is still mostly gravitationally bound, the methods do not even apply (the cluster's potential is dominant over the Galactic one); (iii) if the stars have a significant age dispersion relative to the age of the cluster (in the case of very young clusters) or if they become unbound at different times, the time zero of the methods is poorly defined.
To avoid these limitations, we propose a new method of determining dynamical ages: instead of measuring the expansion age of the whole cluster, we evaluate the evaporation ages of its individual stars, and identify the cluster's age with that of the oldest evaporation ages.The evaporation age is the time since a star escaped the gravitational potential of the cluster.We define this time ignoring the cluster potential, as this requires careful modeling of the time evolution of both the gas and stars in the specific cluster, which is generally unknown, and the effect of the cluster potential is accounted for statistically when the method is calibrated with the simulated clusters.Younger evaporation ages may also correspond to the actual stellar ages if those stars were formed at a later time, or may represent a lower limit to the stellar ages if the stars had remained gravitationally bound to the system for a longer time and evaporated only more recently.For these reasons, we only use the oldest evaporation ages to establish the cluster age.However, the complete age and spatial distributions of the evaporated stars contain valuable information about the formation and dispersion process of a star cluster.
This method has a time zero comparable to that of the CMD method if the first evaporated stars are gravitationally bound to the cluster for a short time relative to their age, and if the observations can identify them, despite their large position and velocity dispersions.The time zero of the evaporation ages is independent of the duration of the bound phase of the cluster (even bound clusters undergo evaporation) and of its age dispersion.Therefore, this method can be applied to bound A165, page 3 of 13 Example of the application of the evaporation-age method to a simulated cluster.a) Projected xz positions of the cluster members in the present (t 0 ).Green triangles are initial cluster members rejected because they do not trace back into the cluster core (black circle).Black dots are cluster members in the cluster core in the present.Blue diamonds and red stars are evaporated cluster members.They are currently outside the cluster core but were inside in the past.Vectors show the 2D projected velocities back in time.b) The diagram illustrating the evaporation age determination.The x-axis is the current distance from the cluster center, d, with the gray vertical line showing the current core radius of the cluster, and the y-axis is the evaporation age, t eva , for each star.Red stars are evaporated stars with t eva higher than the 70th percentile of the t eva distribution, and the red line is the median t eva for these stars.Blue diamonds are evaporated stars with t eva lower than the 70th percentile and are not used for the evaporation age determination.The gray horizontal line is the simulation age of the cluster, t sim .c) Projected XZ positions of evaporated stars at half t sim .d) Positions of evaporated stars at the beginning of the cluster.Only 50% of all the cluster stars had been born by this time.
clusters as long as enough evaporated stars can be identified observationally.The method becomes increasingly uncertain for clusters of increasing age or increasing escape velocity, because of the increasing size of the spatial and velocity windows that must be explored to identify the oldest evaporated stars, and because the past trajectories become more uncertain.
The specific implementation of our evaporation-age method can be summarized by the following steps that are applied to each of the individual clusters identified by the GMM in the simulation snapshots.In the case of the clusters from the simulation, we assume constant velocities equal to the current values because the simulation does not include a Galactic potential.In general, when applied to real clusters, the orbits and the cluster center should be traced back in time accounting for the Galactic potential (see Sect. 6).
1. Rejection of divergent stars: We determine the position and velocity of a cluster's center by taking the mean position and the mean velocity of all the stars in the cluster.The core radius, R 50 , is defined as the median distance of the stars from the center.We then trace each star back in time, until it reaches its closest 3D point to the center.As we do not know apriori the age of a real cluster, we go as far back in time as necessary.If both the star and this closest point are outside of the radius, we reject the star for the rest of the analysis.
2. Subclustering: It is not easy to find a single, optimal number of components for the GMM.The GMM sometimes fails to separate all the clusters, in particular younger clusters that are close together.In order to break these "joined clusters" apart to their actual individual clusters, we perform additional subclustering.We identify dense concentrations of stars, with five stars or more and a maximum mutual distance of 1 pc, as potential subcluster centers (see Appendix A for more details).The cluster's center from the previous step is always included as the first subcluster center, too.This way, if there are no additional subclusters, the original cluster is treated as a single subcluster.We assign each star in the original cluster to the new potential subcluster it passes closest to.For each star member of the new subclusters, we calculate the closest distance to the subcluster center, the star evaporation age, defined as the time elapsed since the star was closest to the subcluster center, and the current distance to the center.With this step, a single GMM cluster may break into several subclusters, if a sufficiently high number of evaporated stars are found for more than one subcluster (see below).
3. Evaporation age: For each subcluster, we select all the stars of the subcluster that trace back inside its core radius and are currently located outside; we refer to the number of these evaporated stars as n eva .We compute the 70th percentile of the star evaporation ages, t 70 , and remove any evaporated stars with evaporation age larger than 1.4 × t 70 + 2 Myr, to reject potential outliers.We iterate until no more stars are removed and then select the stars with evaporation ages larger than the 70th percentile.We compute the evaporation age of the subcluster, t eva , as the median of the evaporation ages of this selection of stars, and record the number of stars used to obtain this value, n eva,70 .We only keep the subclusters with n eva,70 > 1.
The steps of the evaporation-age method are illustrated in Fig. 2, with one of the simulated clusters as an example.The topleft panel shows the current positions of the stars, where we have identified a few divergent stars (green triangles) that do not trace back (green vector) to the core radius.The other stars (blue diamonds and red stars, n eva ) trace back to the core radius (blue and red vectors) or are already inside (black points).In this particular example, we did not find subclusters, so we can directly plot the evaporation ages of individual stars as a function of the current distance, d, to the center, on the top-right panel.We used the stars outside the core radius to obtain the complete evaporationage distribution, but only the stars with an evaporation age equal to or larger than the 70th percentile (red stars, n eva,70 ) are used to obtain the final evaporation-age of the cluster, t eva .In this example, we did not find any outliers with the t 70 criterion mentioned above, and the evaporation age that we determine (red line) is slightly underestimated compared to the simulation age (gray line).The two bottom panels show the positions of the stars that are currently outside the cluster's core radius when traced back in time to half of t sim (bottom-left) and then by a full t sim (bottomright).At t sim , only half of all the stars were born and they may have been gravitationally bound.Thus, some stars overshoot the radius, when traced back to a time before their evaporation time.The dashdot red line corresponds to φ eva,0 t sim , where φ eva,0 is the median value of the age ratio t eva /t sim (Eq.( 1)) in each bin.The scatter, σ eva,0 , is the standard deviation of the ratio of the corrected t * eva (see Eq. ( 2)) of the individual clusters with respect to their t sim , as explained in Eq. ( 3).The black dashed line is a least-squares fit, f (x), showing that the slope is very close to unity and thus φ eva,0 does not depend on the age of the cluster.

Validation and calibration of the method
different ages, we need to calibrate the derived t eva for the specific percentile value we have chosen.We calibrate and evaluate the accuracy of t eva by comparison with the true age of the clusters in the simulation, t sim , which we define as the median age of all the stars in the cluster.We first perform this analysis by applying the method without adding position or velocity errors to the stars in the simulation.To ensure the quality of the evaporation ages, we consider only the clusters with t eva > 2 Myr, as our tests showed that the accuracy decreases significantly with younger clusters.
The values of t eva and t sim for all the clusters in the six snapshots of the simulation are plotted in Fig. 3, with the three panels showing clusters within three n eva,70 bins: 5-9, 10-19, and 20 or more, from left to right.The figure shows a strong correlation between t eva and t sim , increasing with increasing values of n eva,70 .In addition, the least-squares fit (black dashed lines) has a slope very close to unity, so we can calibrate the method by dividing the derived t eva by a scaling factor, φ eva , independent of cluster age, which we compute as the median of all the age ratios, shown by the red dashed-dotted lines in Fig. 3.We will derive the full dependence of φ eva on both n eva,70 and the size of the position and velocity errors in Sect.5.2.For the three n eva,70 bins shown in Fig. 3, the calibration factor without observational noise, φ eva,0 , varies from 0.81 to 0.93.With the established calibration, the corrected evaporation age of a cluster, t * eva , is defined as and the statistical uncertainty, σ eva , is estimated as the standard deviation of the ratio of the individual corrected ages divided by their simulation ages, t * eva /t sim , a measure of the scatter in Fig. 3, As seen in Fig. 3, this uncertainty is actually symmetrical in logarithmic space rather than in linear space.Thus, it should be interpreted as a factor: if σ eva = 0.3, the upper 1σ limit should be 1.3 × t * eva , while the lower limit should be t * eva /1.3.As the parameter φ eva , σ eva is also approximately independent of cluster age for t eva > 2 Myr, and it only depends on the number of evaporated stars (n eva,70 ) and on the size of the position and velocity errors (see Sect. 5.2).For the three n eva,70 bins shown in Fig. 3, σ eva is quite sensitive to n eva,70 , decreasing from 35% to 9%.Thus, with sufficiently small observational errors, the evaporation-age method could in principle yield ages with a statistical error approaching 9%.Because n eva,70 is defined by the 70th percentile, achieving a 9% error or lower requires the identification of at least 67 stars outside of the cluster's core radius with measured evaporation ages (stars tracing back to the cluster's core), so a minimum of 134 stars in total (assuming no star is rejected).Figure 3 shows that only a few of the clusters in our sample have n eva,70 values large enough to achieve this high accuracy.However, if we had all the lower-mass stars in the simulation, our total number of stars would increase by ∼10, yielding a much larger number of clusters in the right panel of Fig. 3.
Figure 4 shows σ eva,0 (σ eva without observational errors) as a function of n eva,70 , where the clusters are binned according to their n eva,70 so that each bin has ten or more clusters, with the exception of the last bin that has 8 clusters.We do a least-squares fit in log-log space, which results in the fitting formula: σ eva,0 (n eva,70 ) = 1.5 n −0.8 eva,70 , in linear space.This formula expresses the statistical error intrinsic to the method, that is in the absence of observational errors.

With position and velocity errors
In this section, we investigate the dependence of φ eva and σ eva on the observational errors in the positions and velocities.Assuming Gaussian error distributions, we added random errors to the current positions and velocities of the stars in the simulation, with six different 1σ values in both position and velocity, hence a 6 × 6 grid of 1σ error pairs (in the ranges 0.15-4.80pc in positions and 0.15-4.80km s −1 in velocities).These values represent the present errors achievable with Gaia only (astrometry and radial velocities) for nearby stars (<500 pc) and also the Gaia astrometry complemented with ground-based highresolution spectroscopy.For example, in Upper Scorpius, the tangential errors are around 0.2 km s −1 , and with ground-based observations plus strict quality criteria (like the ones applied in our samples), it is possible to get radial velocity uncertainties of 1 km s −1 .At the same time, the uncertainties of Gaia's Fig. 4. Scatter σ eva,0 as a function of n eva,70 .The clusters are binned according to their n eva,70 so that each bin has ten or more clusters, with the exception of the last bin that has 8 clusters.The center of each bin is the mean of n eva,70 of the clusters in that bin.φ eva and σ eva are then derived for each bin using Eqs.( 1) and ( 3).The least-squares fit to the [bin center, σ eva ] pairs (blue points) is done in log-log space (see the inset panel), weighing by the square root of the number of clusters in each bin, which results in a fit f (x) = 1.5 x −0.8 in linear space (red line).
radial velocities are a few kilometers per second.We assumed that these error values, σ p and σ v , were the same in the three directions, σ p = σ p,x = σ p,y = σ p,z and σ v = σ v,x = σ p,y = σ p,z .In general, the observational errors in the line-of-sight direction are different (larger for nearby clusters) than in the 2D tangential plane.We have explored a case where the 1D error along the line-of-sight was three times the 1D error in another axis, and found differences in φ eva and σ eva of only 1-10% for the same 1D equivalent total errors, defined as: and Thus, although the following calculations were carried out using equal errors in the three directions, the results hold valid for different error distributions with the same 1D equivalent total errors defined above.We performed 1000 Monte-Carlo runs for each of the 36 error pairs, but we computed the optimal number of clusters (the GMM procedure described in Sect.3) only once for each error pair to reduce the computational time.
Figure 5 shows the impact of the observational errors on the fraction of clusters that we can identify, f GMM , compared to the case with no errors), the fraction of clusters for which we can compute an evaporation age, f ages , the calibration factor, φ eva , and statistical error, σ eva , of the evaporation ages.All these values are the medians of the 1000 individual realizations of each observational error pair.The fraction of identified clusters decreases significantly with increasing σ p and σ v values, going down to nearly 50% with the largest error pair.The fraction of retained clusters, f ages , decreases with increasing σ v , but stays constant or increases with increasing σ p .The increase is particularly noticeable in the cases with n eva,70 ≥ 10 and 20.The reason is two-fold.Firstly, the larger positional errors mean that GMM sorts all the stars into fewer clusters, as seen in the top row of Fig. 5, hence the individual clusters are bigger and retain more stars, so more clusters meet the minimum n eva,70 condition.Secondly, increasing random position errors cause a larger scatter of dense clumps of stars, reducing the chance that large GMM clusters get subclustered into smaller ones.The calibration factor, φ eva , is not very sensitive to the value of n eva,70 or σ p , but decreases significantly with increasing σ v , for σ v > 1 km s −1 .This is due to the shedding of the farthest and slowest evaporated stars: if a star is far away or has a low velocity, even a small deviation of its velocity may make it miss the cluster core and get rejected; if it is close to the cluster core or it has a high velocity the star is more likely retained.This results in a stronger reduction for the cluster's evaporation age than for the cluster's simulation age, as the latter is less sensitive to the loss of the oldest evaporated stars (they may not even be the oldest stars in the cluster).The value of the statistical age uncertainty, σ eva , increases with increasing position and velocity errors, as expected, with some exceptions in the case of n eva,70 > 20, due to the poor statistics (too few clusters).As in the case without errors, σ eva decreases with increasing n eva,70 , approaching 10% for small enough observational errors and large enough number A165, page 6 of 13  of stars.Histograms of φ eva and σ eva for all 36 error pairs are given in Appendix B, as well as the values of the medians of those histograms for five bins of n eva,70 .

Data and method
In this section, we apply the evaporation-age method to nine young (<50 Myr), nearby (<150 pc) clusters, namely, β Pictoris (Miret-Roig et al. 2020), Tucana-Horologium (Galli et al. 2023), and seven clusters in Upper Scorpius and Ophiucus detailed in Miret-Roig et al. (2022).We selected clusters for which we have a recent dynamical traceback age, t DT , obtained from Gaia astrometry plus a precise selection of radial velocities (errors 0.5 km s −1 ) from devoted ground-based observations, APOGEE, and Gaia.To obtain the 3D Cartesian heliocentric positions and velocities for the stars in our sample, we used the Gaia DR3 astrometry and the same radial velocities used for the dynamical traceback ages.This selection was carefully made to include only precise radial velocities ( 1 km s −1 ) and exclude potential binary stars which hinder the traceback.
A CMD age, t CMD , was also computed for these clusters, using the same stars as for the traceback analysis, and the isochronefitting methodology presented in Ratzenböck et al. (2023a), using the PARSEC v1.2S models (Marigo et al. 2017).For most clusters, the values of t CMD were presented in Miret-Roig et al. (2024), while for α Sco, π Sco and σ Sco they are reported here for the first time.The membership and dynamical traceback age for these three clusters are less accurate since α Sco contains two subclusters identified in Ratzenböck et al. (2023b), π Sco is incomplete, and σ Sco is more contaminated (Miret-Roig et al. 2022).In this work, we used the GMM code to separate α Sco into two subclusters, α Sco 1 and 2, by imposing a two-component solution on the original membership list, and apply the evaporationage method to the two subclusters separately.The properties of all the clusters are summarized in Table 1.
To apply the method described in Sect.4, the orbits of the stars and the cluster center are traced back in time using a 3D Galactic potential.We tested three axisymmetric potentials, namely MWPotential14 (Bovy 2015), McMillan17 (McMillan 2017), and Irrgang13 (Irrgang et al. 2013).The differences in the evaporation ages between the Galactic potential models are very small, generally of <1%, and only up to a few percent at most.The results that we report are obtained with the MWPotential14, following the orbits back by a maximum of 50 Myr (100 Myr in the case of Tuc-Hor), which is enough to find the evaporation ages.The initial number of stars used to determine the evaporation ages is the same as that used for the DT and the CMD ages (n stars in Table 1), except for the subclusters of α Sco, where the cluster was split in two.This sample is limited to stars with precise radial velocities, essential for the traceback analysis, but is only a lower limit to the total number of cluster members.Figure 6 shows the number of core, evaporated, and divergent stars for each cluster.Thanks to the careful selection of members  (and 6D phase-space data) in these clusters, the number of stars rejected (divergent) in the evaporation-age method is relatively small (see Table 1).The exceptions are the α Sco and σ Sco clusters, which are likely more contaminated (see Miret-Roig et al. 2022), and ρ Oph and ν Sco, which are the youngest groups in the sample and might still be gravitationally bound.By contrast, the more isolated clusters β Pic and Tuc-Hor only have two divergent stars (in Tuc-Hor); this shows the impact of the proximity of other clusters to assigning the correct membership to the stars.We derived the calibration factor, φ eva , the corrected age, t * eva , and its uncertainty, σ eva , using the 1D equivalent total observational errors and the values of n eva,70 as described in Sect. 5.All of these values are reported in Table 1.As explained in Sect.5, the age uncertainty is symmetrical in logarithmic space, hence used as a factor to get the error bars in linear space.

Evaporation ages
Figure 7 shows the distributions of the evaporation ages of individual stars (red histograms), obtained from the 1000 error realizations of each cluster.The corresponding cumulative distributions are also shown (gray curves), as well as the estimated age, t * eva (black vertical lines).Because the method estimates the cluster age from the oldest evaporation ages of individual stars, the peak of the distribution is expected to be at younger ages than t * eva (to the left of the black vertical line), as shown in Fig. 7.The difference between the time of the peak and t * eva , typically a few Myr, is a rough estimate of the characteristic time during which the evaporated stars were bound to the cluster, which should increase over time, as more recently evaporated stars are gradually released from the cluster's potential.In addition, the width of the distribution sets an upper limit to the actual age spread of the stars (it would be equal to the age spread only if the stars were never gravitationally bound to the cluster).Future studies should attempt a detailed analysis of the distributions of the evaporation ages of individual stars, including a comparison of the evaporation-age spread with the age spread inferred from the clusters' CMD.
Figure 8 shows the comparison of the evaporation ages, t * eva , with the CMD ages, t CMD (left panel), and the dynamical trace-back ages, t DT (right panel).Nine out of ten clusters have an evaporation age compatible with the CMD age within the 1-σ uncertainty.Only for π Sco the evaporation age is not compatible with the CMD age.However, the census of π Sco is incomplete, limited to the stars inside the field of view analyzed in Miret-Roig et al. (2022).The truncation of the field of view leads to a shift of its center toward one side, reducing the stellar evaporation ages for the stars on that side and hence reducing the cluster's evaporation age.This example stresses the importance of obtaining clusters' membership lists as reliable and complete as possible.
In the left panel of Fig. 8, we have performed a least-squares fit to the data in the logarithmic scale, discarding π Sco due to the previously stated reasons.We find a slope of 1.2, although the 1-to-1 line is still within the 1-sigma confidence of the fit.We caution that this slope is based on mainly the two lowest and the two highest age clusters.For the two young clusters, it is possible that we miss some of the evaporated cluster members due to the crowded nature of Upper Scorpius, or perhaps for these clusters there was a longer bound period before any stars were evaporated, either of which would easily explain the shift of about 1 Myr.Also, in all four cases, the clusters' evaporation ages are relying on just 2-4 evaporated stars above the 70th percentile.Thus, the statistics are quite uncertain.
As shown in the right panel of Fig. 8, the evaporation ages are always larger than the dynamical traceback ages.This is expected since our method is designed to find the time when the first stars evaporated.In contrast, the dynamical traceback age is designed to find the overall expansion time of the cluster, which may start well after the evaporation of the first stars.The overall expansion is better represented by the evaporation of the overall stellar population, whose peak evaporation time (the peak of the red histograms in Fig. 7) is a few million years lower than the cluster's age, t * eva , as mentioned above.This result confirms the recent findings from the comparison between CMD and dynamical traceback ages by Miret-Roig et al. (2024).
We have found that the evaporation-age method results in ages compatible with those from isochrone fitting in the CMD.With accurate enough observations (e.g., position errors <0.6 pc and velocity errors <0.6 km s −1 ) and a large enough number The dotted black line is one-to-one line, showing that the evaporation ages are in general agreement with the CMD ages, but the DT ages are clearly offset.On the left panel, the blue dashed line is a least-squares fit to the data (except for π Sco) in logarithmic space, y = 1.23 +0.09 −0.09 x − 0.24 +0.09 −0.09 , with the shaded region showing the 1-sigma confidence of the fit.Apart from the youngest ages at <3 Myr, the 1-to-1 line is still within the shaded region.The values used and references can be found in Table 1.The t CMD and t DT ages of α Sco 1 and 2 have been shifted right and left by a factor of 1.02 to prevent the error bars from overlapping; the age reported in the table is in between.On the right panel, the symbol for ρ Oph is out of the plot to the left, as it has t DT = 0. of stars (n eva,70 ≥ 20), the evaporation-age uncertainty could approach a value as low as 10%, comparable to the typical statistical uncertainty in t cmd from isochrone-fitting.The real uncertainty in t cmd , particularly for t cmd < 10 Myr, is significantly larger than 10%, due to uncertainties from pre-main sequence modeling, such as the effect of starspots and magnetic fields (e.g., Simon et al. 2019;Somers et al. 2020;Cao et al. 2022), or the effect of extended and variable accretion histories (e.g., Froebrich et al. 2006;Baraffe et al. 2009;Hosokawa et al. 2011;Jensen & Haugbølle 2018).However, the good correlation between our evaporation ages and the CMD ages may suggest that systematic errors in the latter may not be as large as the full range of age values allowed by the largest differences between evolutionary models.As follow-up observations continue to increase the samples of accurate radial velocities in young clusters, isochrone-fitting ages should be carefully compared with evaporation ages, making sure that the same stars are used in the two methods.This comparison may shed light on the pre-main-sequence evolution (e.g., the duration of the accretion phase), on the age spread within a cluster, and on the gasdispersal mechanism.

Conclusions
1. Using the clusters from the simulation, without adding observational errors, we find a strong correlation between the estimated age of a cluster and the median simulation age of its stars.We derive the ratio between estimated and real age, φ eva,0 , which we use to calibrate the method's age.2. Using the simulation we also estimate the intrinsic (not due to observational errors) statistical uncertainty of the method, σ eva,0 .We find it has no significant age dependence, and provide a fitting function for its dependence on the number of evaporated stars.We show that σ eva,0 can be as small as 9% for clusters with a large number of evaporated stars.3.By adding observational errors to the simulated clusters, we derive both φ eva and σ eva as a function of both the number of evaporated stars and the observational errors in position and velocity, with σ eva being the total uncertainty of the derived age, including both the intrinsic and observational random errors.4. We find that, assuming small observational errors (σ p < 1 pc and σ v < 1 km s −1 ) similar to the ones achievable with Gaia astrometry plus complementary high-resolution spectroscopy, cluster ages have a 1-σ uncertainty <30% for clusters with a minimum of ∼70 stars, and approaching 10% for clusters with ∼150 stars or more, if only a small fraction of the stars are rejected.5.A pilot application to ten nearby clusters results in evaporation ages consistent with CMD ages, and larger than previously determined dynamical traceback ages.The uncertainties are currently significantly larger than 10%, due to the limited number of confirmed evaporated stars in the observational samples.Given the large number of current and upcoming radial-velocity surveys complementing Gaia's database, evaporation ages with accuracy close to 10% may become the norm in a few years.However, it will be necessary to explore rather wide space and velocity windows around each cluster, to capture a significant fraction of the evaporated stars, particularly the earliest ones.In practice, this will limit the method to relatively young clusters, A165, page 9 of 13 probably with ages <50-100 Myr, for which we can determine membership and precise stellar orbits back in time of their earliest evaporated stars.Because of theoretical uncertainties (e.g., magnetic fields and starspots) and stochasticity (e.g., extended and episodic accretion) of the pre-main sequence evolution of stars, CMD ages for very young clusters (<10-20 Myr) are very uncertain, despite the often-quoted statistical error of ∼10% from isochrone fitting.As it is completely independent of stellar evolutionary models, and particularly suited for very young clusters, the evaporation-age method will provide important constraints for the modeling of pre-main sequence evolution.A165, page 13 of 13 Fig.2.Example of the application of the evaporation-age method to a simulated cluster.a) Projected xz positions of the cluster members in the present (t 0 ).Green triangles are initial cluster members rejected because they do not trace back into the cluster core (black circle).Black dots are cluster members in the cluster core in the present.Blue diamonds and red stars are evaporated cluster members.They are currently outside the cluster core but were inside in the past.Vectors show the 2D projected velocities back in time.b) The diagram illustrating the evaporation age determination.The x-axis is the current distance from the cluster center, d, with the gray vertical line showing the current core radius of the cluster, and the y-axis is the evaporation age, t eva , for each star.Red stars are evaporated stars with t eva higher than the 70th percentile of the t eva distribution, and the red line is the median t eva for these stars.Blue diamonds are evaporated stars with t eva lower than the 70th percentile and are not used for the evaporation age determination.The gray horizontal line is the simulation age of the cluster, t sim .c) Projected XZ positions of evaporated stars at half t sim .d) Positions of evaporated stars at the beginning of the cluster.Only 50% of all the cluster stars had been born by this time.

Fig. 3 .
Fig.3.Evaporation age, t eva , versus simulation age, t sim , for the clusters in the simulation within three different bins of the number of stars, n eva,70 , from which the evaporation age was calculated.The left panel has 86 clusters, the middle one has 45 clusters, and the right one has 15 clusters.The dashdot red line corresponds to φ eva,0 t sim , where φ eva,0 is the median value of the age ratio t eva /t sim (Eq.(1)) in each bin.The scatter, σ eva,0 , is the standard deviation of the ratio of the corrected t * eva (see Eq. (2)) of the individual clusters with respect to their t sim , as explained in Eq. (3).The black dashed line is a least-squares fit, f (x), showing that the slope is very close to unity and thus φ eva,0 does not depend on the age of the cluster.

Fig. 5 .
Fig. 5. Impact of the observational errors.Top row: fraction of GMM clusters, f GMM , derived from the number of clusters found by GMM in each error pair, divided by the number of GMM clusters in the no errors case (339).Second row: fraction of clusters for which we derived ages, f ages , normalized to 339.Red, blue, and black lines are for n eva,70 bins of 5-9, 10-19, and 20+ stars, respectively.Third row: φ eva .Bottom row: σ eva .Columns correspond to different 3D position errors, σ p , and the x-axis is the 3D velocity error, σ v .

Fig. 6 .
Fig.6.Number of stars in each of the ten observed clusters found within the core radius, n core (blue), found outside of the core radius and tracing back to it, n eva (orange), or not tracing back to it, n div (green).

Fig. 7 .
Fig.7.Histograms (peaks scaled to 1.0) of the evaporation ages of individual evaporated stars (red histograms) and the corresponding cumulative distributions (gray curves), from the 1000 Monte-Carlo error realizations and using stellar orbits computed with MWPotential14.The derived t * eva Fig. B.3.Maps of the medians of φ eva (left column) and σ eva (right column) histograms in the previous figures as functions of σ v and σ p , for five bins of n eva,70 (rows).
Fig.1.Histograms (in red) of the ages, masses, radius, and rms velocities and scatter plots of the square root (to use the same scale as the histograms) of the ratio of kinetic and gravitational energies versus the same quantities (blue dots) of the 339 clusters identified in the six simulation snapshots used in this work.The black dashed line marks the value of E k = E g , showing that 5% of the clusters (17 of them) are gravitationally bound.

Table 1 .
Properties of the clusters considered in this study.R 50 n stars n div n eva n eva,70 φ eva σ eva Comparison of the corrected evaporation ages, t * eva , with the color-magnitude diagram ages, t CMD , and the dynamical traceback ages, t DT .