The Age Distribution of Stellar Orbit Space Clumps

The orbit distribution of young stars in the Galactic disk is highly structured, from well-defined clusters to streams of stars that may be widely dispersed across the sky, but are compact in orbital action-angle space. The age distribution of such groups can constrain the timescales over which conatal groups of stars disperse into the “field.” Gaia data have proven powerful in identifying such groups in action-angle space, but the resulting member samples are often too small and have too narrow a color–magnitude diagram (CMD) coverage to allow robust age determinations. Here, we develop and illustrate a new approach that can estimate robust stellar population ages for such groups of stars. This first entails projecting the predetermined action-angle distribution into the 5D space of positions, parallaxes, and proper motions, where much larger samples of likely members can be identified over a much wider range of the CMD. It then entails isochrone fitting that accounts for: (a) widely varying distances and reddenings; (b) outliers and binaries; (c) sparsely populated main-sequence turnoffs, by incorporating the age information of the low-mass main sequence; and (d) the possible presence of an intrinsic age spread in the stellar population. When we apply this approach to 92 nearby stellar groups identified in 6D orbit space, we find that they are predominantly young (≲1 Gyr), mono-age populations. Many groups are established (known) localized clusters with possible tidal tails, while others tend to be widely dispersed and manifestly unbound. This new age-dating tool offers a stringent approach to understanding on which orbits stars form in the solar neighborhood and how quickly they disperse into the field.


INTRODUCTION
Stars are born in a highly clustered fashion, in position, velocity, and orbit space (e.g.Lada & Lada 2003;Lamb et al. 2010;Kounkel & Covey 2019a).It appears that the vast majority of them do not remain in bound clusters after the gas expulsion phase (e.g.Chandar et al. 2017).And even those that first remain in "open" clusters beyond this phase tend to eventually become unbound from these clusters, by internal twoor three-body interactions, or by external tides.Consequently, most stars with ages > 5 Gyrs have become "field stars".However, how clustered (bound or un-Corresponding author: Verena Fürnkranz fuernkranz@mpia.debound) the full population of stars was in their first, say, 10 Myrs, and how, or how rapidly, the transition from star clusters and associations to the Galactic field takes place has yet to be fully understood.Studying this transition process will open a new window on constraining whether all stars are born in clusters, give insights into how hierarchical star formation is, and allow one to tackle how rapidly the dissolution of clusters takes place.
The Gaia mission has provided data that are revolutionary for our understanding of the formation and structure of the stellar component of the Milky Way disk (Gaia Collaboration et al. 2016).The third Gaia data release (Gaia DR3) on June 13, 2022 has provided precise 6D position measurements (α, δ, ϖ), proper motion (µ α , µ δ ), and radial velocity measurements (v r ) for millions of stars in the extended solar neighborhood (Gaia Collaboration et al. 2022a).Recent work, e.g. by Trick et al. (2019) and Coronado et al. (2020), had shown that this 6D parameter space -or equivalently orbit spaceis highly clustered over a wide range of scales.Largescale structures might be caused by resonances of the Galactic bar and spiral arms, or by satellite infall into the Milky Way Galaxy (Gaia Collaboration et al. 2018).Small-scale clustering represents the concepts of open clusters, associations, and stellar (disk) streams.
In the immediate neighborhood, the space velocities can serve as orbit proxies.On larger scales, action-angle coordinates (J, Θ) stand out as a powerful coordinate system to search for substructure in the Galactic orbit and orbital phase space: the three actions denote the "orbit"; the three angles the "orbital phase".Recent work by Coronado et al. (2022) has shown that we can find numerous stellar groups by linking stars in this 6D parameter space: some are simply clusters, others extend over 100's of parsecs.In Coronado et al. (2022) we discovered that the same orbits that host one of these groups very frequently host numerous other distinct orbit-phase overdensities along these same orbits, which we call pearls 1 .This might open up a new way to probe the filamentary structure of the ISM.
The central overall goal of our present work is to associate ages to stellar action-angle groups in the Galactic disk.These are crucial data for understanding for how long ensembles of stars remember their initially near-identical orbits and orbital phases.Once sets of stars (clusters or associations) are no longer gravitationally bound but still tightly clustered in orbit space (J, Θ), they will continue to disperse into the field by simple phase-mixing and by further relative orbit diffusion, driven e.g. by the Galactic tidal field.To make predictions about the rate of this dispersal, age information about the ensembles of stars is needed.
For coeval ensembles of stars, fitting isochrones to their CMD positions is a long-established approach to age estimates.This suggests that one should simply fit the CMDs of the stellar ensembles identified (see Coronado et al. 2022) in 6D phase space.However, radial velocities from Gaia exist only for stars over a limited range of effective temperatures (3.100 to 14.500 K) and for stars brighter than G RV S ≈ 14 mag.Furthermore, the radial velocities for hot and cool very young stars are not very precise (δv r ≈ 5 km s −1 ), which is a serious impediment in orbit clustering studies.At any rate, this means that the CMDs of these 6D sample are incomplete or missing for the uppermost main sequence (in-1 These pearls are strung along one orbit, like pearls on a string. cluding potentially the turn-off for ages τ < 200 Myrs) and for the lowest-mass stars, which are also good age diagnostics as stars may take > 10 8 years to settle on the zero-age main sequence.
The practical goal of this paper is to overcome this limitation by projecting the tight 6D action-angle distribution of these stellar ensembles to 5D position-proper motion space, where we can then establish much larger membership samples with higher completeness at the extreme mass and temperature ranges that are so precious for age estimates.This is possible since Gaia is highly complete to G ≈ 19 mag across a very broad range of colors, if only positions, parallaxes, and proper motions are required.We aim to end up with CMD's that are much better populated in both total sample size and color range.This requires that the groups, high-density peaks in 6D (J, Θ), still have a distinct density contrast in their 5D projection (α, δ, ϖ, µ α , µ δ , but without v r ).
Conceptually, our approach to extending group membership is as follows: we start with a set of stellar groups identified in 6D space.Specifically, we get these by identifying 6D action-angle groups within 800 pc from the Sun in a similar fashion to Coronado et al. (2022).We use the DBSCAN clustering algorithm and choose a search setting that identifies 102 action-angle groups of stars, which is about twice the sample size as the 55 action-angle groups identified in Coronado et al. (2022).For our analysis, we pick all groups with at least 20 member stars.This reduced set of 94 groups largely overlaps with the open cluster catalog published in Cantat-Gaudin et al. (2020), as well as with the stellar groups detected in Hunt & Reffert (2023).Additionally, almost all of our stellar groups can be found in the catalog of Kounkel & Covey (2019b), which is known to be complete but highly contaminated (see Zucker et al. (2022)).We project these groups into the 5D positionproper motion space, where we have a factor ∼8 more stars in the Gaia catalog.For each action-angle group, we assign additional candidate members that are located in the same 5D region.Specifically, we run a kernel density estimate on the 5D distribution of each group and pick all stars in the Gaia catalog that show an excess at high densities, compared to stars in the Gaia mock catalog (Rybizki et al. 2020).
We then take 92 of these groups with their extended membership to develop, test and apply our new approach to age dating these groups of potentially conatal stars, accounting for their potentially sparse and widely dispersed nature.This approach entails likelihood isochrone fitting that takes into account a) differential reddening, b) the widely ranging member distances, c) the individual photometry and distance errors, as well as d) binaries and outliers.We also test whether groups can be described as effectively mono-age, rather than have evidence for a large age dispersion in their CMD.Accounting for all these aspects in age-dating is crucial in this regime, and has not been implemented before.
With this new approach, we determine the ages and possible intrinsic age spreads of the identified stellar action-angle groups, and find them to be predominately young (τ ≲ 2 τ orbit ) and mono-age.The analysis of the dispersal of these groups in orbit and orbital phase-space will be subject of upcoming work, as will be the age distribution of any "pearls" along the same orbit, as well as the smallest and most diffuse groups we can identify that still show mono-age CMD sequences.
The paper is organized as follows.In Section 2 we describe the data used in this study.Section 3 outlines our method for finding groups in 6D action-angle space and associating member stars by extending our search to the 5D position proper motion space.In Section 4 we present our likelihood isochrone fitting method which determines ages of the groups.Our results are discussed in Section 5. Finally, we summarize our findings and conclude in Section 6.

Excerpts of the Gaia DR3 catalog
We base our analysis exclusively on Gaia DR3 (Gaia Collaboration et al. 2022b).As our intention is to study the extended solar neighborhood, we limit the analysis to stars with parallax ≥ 1 mas, corresponding to a maximum distance of 1 kpc from the Sun.We clean this sample from sources with spurious astrometry, following Rybizki et al. (2022), and exclude sources with fidelity v2 < 0.9.We also limit the analysis to stars with parallax over error ≥ 10, to eliminate distance uncertainties as the main source of orbit uncertainty.This results in a database containing 41.138.028sources, and we will refer to this catalog as the Gaia sample throughout.
As we start our analysis by identifying groups in the 6D action-angle space, we also introduce a slightly more confined sample with ϖ ≥ 1.25 (within 800 pc from the Sun), comprising only the stars that also have a radial velocity (v r ) measurement.We adopt the same fidelity cut, but relax it to parallax over error ≥ 3 to get a larger 6D sample.This 6D database contains 9.080.884sources, and we will refer to it as Gaia RVS sample.
As we need to fit isochrones to dereddened colors and magnitudes, we query extinction values (A G , A BP , A RP ) from the astrophysical parameter catalog pro-duced by the Apsis processing chain developed by the Gaia DPAC (Fouesneau et al. 2022).Specifically, we query all stars with an extinction entry available, and within 1 kpc from the Sun.Additionally, we limit the catalog to zero-age main sequence stars and apply a 2 mag < M G < 8 mag cut (see Section 4.1), which reduces the catalog to 23.543.432sources and hereafter we will refer to it as Gaia extinction sample.For the most part, it encompasses the Gaia sample.

Gaia EDR3 mock catalog
Since we will identify overdensities in 5D positionproper motion space, we will need to know what densities we should expect from a smooth distribution.To this end, we make use of the Gaia EDR3 mock catalog by Rybizki et al. (2020), which presents an analogous dataset to the Gaia DR3 catalog, but with a smooth and fully phase-mixed orbit distribution.We apply the same distance cut by only selecting stars with ϖ ≥ 1, and add WHERE popid != 11 into the ADQL query, which removes clusters and moving groups.These are otherwise included in the catalog by default.The resulting dataset comprises 46.830.045stars, and will be named Gaia mock sample throughout the manuscript.The size of the catalog is comparable to the size of the Gaia sample catalog, and therefore it is well suited for comparison.

PARSEC isochrones
For age determinations via isochrone fitting, we need to choose one (or more) isochrone models.After considerable experimentation, we settled on PARSEC isochrones (Bressan et al. 2012) for our analysis, as they provide the most consistent age estimates between the turn-off and the low-mass main sequence.MIST isochrones (Dotter 2016;Choi et al. 2016) provide less consistent age estimates for our set of action-angle groups, as they are not optimized for the lower main sequence.
Specifically, we download the default isochrone tables for a set of 6 metallicities in the range of [Fe/H] = −0.4 to +0.1, and 107 ages from log τ = 5 to 10.3, resulting in an isochrone table comprising 224.951 entries.
Our analysis requires that each individual isochrone can return finely spaced magnitude and color values, which we ensure by linear interpolation.Specifically, for each individual isochrone we loop through all mass points and interpolate Mass, G mag , BP mag , and RP mag whenever ∆G mag > 0.03 mag between two mass points.After applying this method, we doubled the isochrone entries to a total of 458.181.We will refer to this setup as isochrone table in the text.

STELLAR GROUPS IN ACTION-ANGLE SPACE
We now describe membership identification of the stellar groups whose ages are at the heart of our analysis: their initial identification as prominent overdensities in the 6D action-angle space and the subsequent extended membership identification after projection of each action-angle distribution ("patch") into the 5D position-proper motion space.

Action-angle computation
The calculation of actions (J) and angles (Θ) requires both the full 6D position-velocity information in Gaia's original coordinates (α, δ, ϖ, µ α , µ δ , v r ) and an assumed gravitational potential.For the potential, we adopt Galpy's axisymmetric MWPotential2014 (Bovy 2015), which has adopted a circular velocity of 220 km s −1 at the Solar radius of 8 kpc, and which comprises a bulge, a disk, and a dark-matter component (Bovy 2015).We then use the Stäckel Fudge algorithm by Binney (2012) to estimate the actions and angles for all stars in our Gaia RVS sample.

Identifying Stellar Overdensities in 6D
Action-Angle Space In Coronado et al. (2022) we showed that stellar overdensities -from clusters, to associations, to disk streams -can be identified as groups that are overdense in action-angle space.Here, we follow the same approach applied to our Gaia RVS sample, with modifications that we summarize here.First, we describe all orbital angles via their sine and cosine functions bypass complexities that arise from the periodicity of these variables.To reach a consistent metric across all six coordinates that serves our purpose, we first normalize each of the action-angle coordinates (J R , J ϕ , J Z , sin(Θ R,ϕ,Z ), cos(Θ R,ϕ,Z )) by removing the mean and scaling to unit variance.Additionally, we explore changes in the relative weight of "orbits" vs "orbital angles", and reduce the variance in the six angle coordinates (now sines and cosines of the angles) by multiplying the scaled data with 0.4.
To identify clusters of objects in this 6D space, Coronado et al. (2022) had used a version of the Friends-of-Friends algorithm to link objects.Here, we apply the density-based algorithm DBSCAN (Ester et al. 1996) to extract clustered ensembles of objects from the (J R , J ϕ , J Z , sin(Θ R,ϕ,Z ), cos(Θ R,ϕ,Z )).We choose this clustering algorithm instead of the Friends-of-Friends algorithm because it is much faster in handling large datasets.Running DBSCAN requires two choices: a minimum number of members (minPts), and a maximum distance scale ϵ.We have explored a wide range of ϵ and minPts parameters, and decided on a final setup of ϵ = 0.035 and minPts = 20.In combination with our orbital angle weighting, this parameter choice yields to 0.1% of the stars in the Gaia RVS sample to be clustered in about twice as many groups as we had identified in Coronado et al. (2022).With these choices we now identify 102 action-angle groups, and keep the 94 most prominent groups that have at least 20 member stars.This sample is large enough to develop, test, and apply our age estimates; and it is large enough to get a first census which fraction of groups are "young" and approximately "mono-age".
Fig. 1 illustrates these 94 action-angle groups: the top row shows their distribution (color-coded by group) in the three action plane projections, and the middle row shows their distribution in the three angle plane projections.By construction, all groups are very compact in this parameter space.The bottom row shows a top-down view on the 94 groups in Galactic Cartesian Coordinates (left), and their distribution on the sky in Galactic longitude and latitude (right).Although the groups are tightly clustered in action-angle space, most of them extend many degrees on the sky, and a substantive fraction over several 100 parsecs in physical extent.
For six representative action-angle groups Fig. 2 displays the CMDs for their members identified in 6D space.In all cases, these members form a narrow sequence in CMD space, as expected for simple stellar populations.However, the overlaid isochrones in Fig. 2 also illustrate that these stars only cover a portion of the expected isochrone locus.This is for two reasons: the limited magnitude range and the limited temperature range of the Gaia RVS sample, where only stars with G < 14 mag and effective temperatures in the range of 3.100 to 14.500 K are included.For the age determinations of these groups this is a serious limitation, as the (upper) main-sequence turn-off (for ages > 10 8 yrs) and the lower main sequence (for ages < 10 8 yrs) provide the best age constraints for stellar ensembles.This is apparent from the set of three isochrones overplotted in Fig. 2, spanning an age range of τ = 60 -630 Myrs, specifically log τ = 7.8, 8.4, 8.8.E.g. for Group 58 in the top row, the group members identified in 6D are located exclusively on the part of the main sequence where all these isochrones coincide.
Hence, determining ages by isochrone fitting requires identifying member stars at temperatures (or stellar masses) that are not included in the Gaia RVS sample.Such hotter and cooler stars that are missing here should, however, be in the overall Gaia sample, with valid 5D position -proper motion entries.

Associating additional member stars in the 5D space of positions and proper motions
In this Section we describe how to extend the stellar membership of the action-angle groups across a wider color or stellar-mass range, in order to overcome the just-described limitations in age dating these groups.
We do this by mapping their distribution from 6D (rescaled) action-angle space to 5D position-proper motion space (X, Y, Z, µ α , µ δ ), to search for the most likely members in this projected space, where we do not have the color limitations of the Gaia RVS sample.We start by constructing an approximate Probability Density Function (PDF) for each group in 6D space.We do this with the help of a kernel density estimate for each group in the rescaled 6D space of J R , J ϕ , J Z , sin(Θ R,ϕ,Z ), cos(Θ R,ϕ,Z ), using a Gaussian kernel with a unit bandwidth that we can then sample finely.Before doing that, we remove members in each action-angle group that are beyond 2σ from the median in any coordinate, to mitigate the impact of outliers; this leaves about 95% of the stars.We then sample this PDF estimate with 5000 random points that we can then map from action-angle space (J R , J ϕ , J Z , Θ R , Θ ϕ , Θ Z ) to positionproper motion space (α, δ, ϖ, µ α , µ δ , v r ).For this mapping we use the MWPotential2014 in the TORUSMAP-PER package (Binney & McMillan 2016), which is part of the AGAMA library (Vasiliev 2019), a software framework offering several tools applicable to galaxy modeling and stellar dynamics.We then remove some "nonphysical" sampled datapoints, e.g.datapoints with negative radial and vertical actions.
Finally, we take the sampled datapoints that have been mapped into 5D space and perform a kernel density estimate, where we choose (X, Y, Z, µ α , µ δ ) as coordinates, to have an intuitive metric for the spatial part.
In each dimension, we use a Gaussian kernel with a unit bandwidth to obtain a density estimate.This smooth 5D density now allows us to assign a local density value ρ to each star in the actual Gaia sample that is within the groups individual boundaries, e.g.minimum and maximum coordinates in (X, Y, Z, µ α , µ δ ).
Of course the 5D space just described can and will also contain unrelated field stars.We can only identify new members of this group with any confidence, where the 5D density of the actual Gaia sample is higher than the expectation from a smooth mock catalog.Therefore we apply the same procedure to the smooth Gaia mock sample described in Section 2.2, leaving us with a local density value assigned to each star in the mock sample.Now we can identify the domain of this 5D space where each group's 5D density is far higher than the density in the smooth Gaia mock sample.To this end, we create a histogram of the number of Gaia (n gaia ) and Gaia mock (n mock ) stars as a function of their assigned Kernel density ρ.We set a density threshold for each group with the following approach: The first density threshold T 1 is the lowest density bin where the local density of the Gaia sample is 25 times higher than expected from the Gaia mock sample.Given that the Gaia mock sample is an imperfect representation of the real data, we decide on the ratio number 25 after some experimentation.This selects all stars in a group as likely 5D members, if their local density is 25 times higher than expected from the Gaia mock sample.Second, we calculate the ratio between n gaia and n mock in each bin and set the second density threshold T 2 at the density with the highest contrast, with the additional condition to contain at least 100 stars.We pick a final density threshold T f inal = min(T 1 , T 2 ) for each group by picking the threshold at the lower density (minimum), e.g. to assign as many 5D candidate members as possible to each group.
We find that about half of the groups (41) get the density threshold T 1 assigned, whereas the remaining groups (53) get the density threshold T 2 assigned.From the groups with T 2 , 16 groups get the density threshold at 100 stars, and 37 groups get a density threshold at lower densities with a higher ratio.
For each group we then merge our original group member selection drawn from the Gaia RVS sample with the new group member selection obtained from the full Gaia sample.From this final group member selection, we manually remove two groups that show overpopulated CMD's and end up with a final set of 92 groups with extended memberships to test and apply the age dating method described in Section 4. The full "group member selection" is available at the CDS, including the Gaia DR3 Source ID, the Group membership, and the coordinates α, δ.
The above described procedure is illustrated in Fig. 3.It shows the density histograms for all Gaia stars (turquoise) and Gaia mock stars (violet) for six representative action-angle groups, selected within the individual group's volume and proper motion range.Most Gaia and Gaia mock stars have small density values, with a tail extending to higher densities.As expected, the actual Gaia stars -in this 6D group selected patch -far exceed the number of Gaia mock stars in the high density regime.This illustrates immediately that we find these groups also as prominent overdensities in the 5D Gaia sample.The dotted blue line marks the location of the density cut T f inal , which varies considerably between the groups, depending on the distribution of the Gaia mock stars.Some groups exhibit a more pronounced excess of Gaia stars in the high density regime than other groups.This could be for two reasons: either some groups might be more compact, or some groups may occur far away from the disk plane, which both provides naturally a higher density contrast.

AGE DETERMINATION
We now want to determine stellar ages of the 92 action-angle groups.In the end this estimate comes from a forward model of the photometry, ⃗ data, of the group members, p τ |{ ⃗ data} , where ⃗ data incorporates the magnitudes and parallaxes of the stars, as well as their errors.The magnitude errors are not provided in the Gaia catalog, therefore we obtain them from the G, BP and RB band fluxes and flux errors.In this modelling, we want to account for the following effects: The stars are widely spread across the sky (see Fig. 1), they may have widely ranging member distances (we obtain distances by inverting the parallax), and they have varying photometric and distance uncertainties.In addition, we need to allow for (unresolved) binaries, possible outliers, and we need to take into account differential reddening.

Reddening correction
We start with the reddening correction for each individual star.It is tempting to simply take each group member's reddening estimate directly from the Gaia extinction sample.These estimates are based on the idea to deredden each star until it gets to its most plausible place in the CMD.However, for stars that lie above the solar-metallicity, zero-age main sequence -whether they are young stars or binaries -this often leads to overesti-mates of their reddening.Instead, we adopt the ensemble median of other stars in the same direction and at the same distance.To obtain robust estimates, we limit these reddening estimates to stars in the Gaia extinction sample where we can assume that they are firmly on the zero-age main sequence, requiring 2 mag < M G < 8 mag.We select all stars that are located within a sphere of radius 25 pc around each group member star under consideration for dereddening, and pick the 50 th percentiles in A G , A BP , and A RP to subtract from the member star's uncorrected Gaia magnitudes G, BP, and RP.Depending on the location of each group and the location of each stars within the group, the spheres contain between 266 and 2059 stars, typically around 800 stars.The effect of this extinction correction is illustrated for Group 68 in the left panel of Fig. 4. The colored data points represent the CMD before applying the reddening correction, and the black data points show the reddening corrected sequence.E.g. for Group 68, the average shift in color is 0.1 mag.Fig. 4 illustrates that reddening correction is important for precise isochrone fitting, even for stars located within the extended solar neighborhood of 800 pc from the Sun.

Basic Fitting Methodology
We now perform isochrone fitting by calculating the likelihood (here χ 2 ) of the reddening-corrected photometric data and parallax for each group, { ⃗ data} i=1,N⋆ .Here, N ⋆ is the number of likely 5D members in the group and each star's data is given by its apparent magnitudes and parallaxes, as well as their errors: Each isochrone, specified by the two parameters log τ , [Fe/H], provides model predictions where band stands for G, BP and RP, and M is the initial stellar mass.
A priori, we neither know any star's mass nor its exact distance, and we should treat these quantities as parameters, even if for the current context they are "nuisance parameters".This means we have to optimize both the global parameters logτ, [Fe/H] and the individual masses M and distance moduli (DM), which are mainly constrained by the parallaxes, ϖ.We perform this optimization by a multi-dimensional grid search.To speed up the computation we only consider a very limited range of DM-M combinations for each star, constrained by their apparent magnitude and parallax.
We only consider distance moduli that are 2σ consistent with the parallax measurement: The plot illustrates that very young as well as very old ages are not likely to fit the group's distribution, as they lead to a high χ 2 value.The best fit isochrone for this group is the one with the lowest χ 2 value.The stars in magenta mark group members that are labelled as 4-σ outliers by our code setup, and are not considered in the likelihood isochrone fitting.We show the isochrones not as solid lines, but as datapoints to illustrate our finely spaced magnitude and color setup, which we obtain by linear interpolation (see Sec. 2.3).
and we only consider (pre-computed) M G (M) in the range (2) Finally, to treat the parallax measurements as data we need to have a parallax prediction, which is simply (3) On this bases, we can calculate for each assumed parameter set and each star χ 2 DM, M; logτ, [Fe/H] : (4) For any set of logτ, [Fe/H] (the outer loop in the optimization) we optimize DM, M for each star, and then sum all stars' χ2 to obtain χ 2 tot logτ, [Fe/H] .This procedure then results in a mapping of the whole data set for any given group as a function of logτ, [Fe/H], including their best fit values.
The right panel in Fig. 4 illustrates the outcome of this process for Group 68.In color are shown isochrones for 107 different log τ , at a fixed metallicity of [Fe/H] = 0.1.The color coding represents the χ 2 value obtained for the different age setups.The plot illustrates that the ages are well constrained both the upper and lower main sequence, as both regimes lead to a high χ 2 value.
This procedure is straightforward and accounts for the varying distances, reddenings and data uncertainties (including the parallax).We now turn to the generalization of this procedure to incorporate outliers and to consider the possibility of an intrinsic age spread in the stellar groups.

Accounting for Binaries and Outliers
Isochrones reflect the color-magnitude locus of single stars, while observed CMDs have "outliers" for a num- The magenta isochrone at the lower main-sequence represents the age of the Galactic disk.The three groups in the top row are very young (< 25 Myrs), and their lower main sequences are clearly offset from the age of the Galactic disk.This again illustrates the power of the lower main sequence to determine the ages of very young populations.The three groups in the bottom row are older, and for them the main-sequence turn-off is the main age determinant.ber of reasons.Such outliers may be spurious data or stellar systems whose physics is not captured by the isochrone: (unresolved) binary stars, evolved stars such as white dwarfs, or actual interlopers of different ages and metallicities.In addition, the color prediction of the lowest mass stars may be less robust for isochrone models, which may explain some of the deviations.
In our code, we implement two different approaches to deal with binaries on the one hand, and different outliers on the other.For binaries, we follow Coronado et al. (2018), where we generalize the isochrone prediction of a star's absolute magnitude M iso to be either that of a single star within a magnitude width σ M , or -with probability ϵ -to be that of a near-equal mass binary that is 0.7 mag brighter, or to be a binary of any kind between M iso and M iso + 0.7 mag.Note that some groups illustrated in Fig. 5 exhibit well populated binary sequences, illustrating the need for this methodological extension.
This results in a probability ( where G is a Gaussian and I is the "in-between" term defined as 1/0.7 between M iso and M iso + 0.7. The logarithm of p(M) in Eq. 5 then simply gets added for each object to the χ 2 in Eq. 4.
We deal with the other outliers by limiting all of them to be "n-σ outliers" at worst: any individual χ 2 -term cannot be larger than 2 × n 2 .We adopt n = 4, or χ 2 max = 32.This is illustrated in Fig. 4 for example Group 68, where stars in magenta mark group members that are labeled as 4-σ outliers by our code setup, and are not considered in the likelihood isochrone fitting.

Incorporating Age dispersion
The approach described so far determines ages, assuming a priori that we have mono-age populations.However, we also want to consider the case that groups may have some finite age dispersion: p logτ | logτ , σ logτ ≡ G logτ − logτ , σ logτ , where the mean (log) age is logτ and its dispersion is σ logτ .We implement this by an extension of Eq. 4: We consider the "linearized" version of χ 2 logτ, [Fe/H] , i.e. exp χ 2 (logτ, [Fe/H]), and then integrate over all logτ within ±3σ logτ : In this case we optimize globally over three parameters, the mean age logτ , the metallicity [Fe/H], and the intrinsic age spread σ logτ .
We then estimate the uncertainty δ log τ of the isochrone fit for each group and compare it to the in- dividual intrinsic age dispersions σ log τ , in order to determine what value of σ log τ is significant ("effectively mono-age").We use the χ 2 values to determine the data probabilities and then look at the marginalised distributions to get the uncertainties in the individual parameters.Rather than quoting them individually, we find that the fit uncertainty δ log τ is in 90% (83) of the cases smaller than 0.1, which is less than our log τ grid size of 0.1.In 9 cases the fit uncertainties are larger than 0.1, whereas in only 1 case δ log τ is larger than 0.15.Therefore, we choose a value of σ log τ = 0.2 to declare "effectively mono-age" as significant.

RESULTS
With these isochrone fitting tools in hand, we can now proceed to presenting the results from the CMD's for 92 groups described in Sect.3.

Illustrating the Age Constraints from Isochrone Fits
We start by illustrating the results derived from the above procedure for the same six groups from Fig. 2, now shown in Fig. 5, with the full set of 5D proper motion members: the isochrones are now populated over a wider age range, providing better age constraints.All stars have been individually reddening corrected.
On top of each group the best fitting single-age isochrones are shown in turquoise (log τ , [Fe/H]), derived following Section 4.
Across all groups, the isochrones fit both the upper and lower main sequences very well.In particular, this Figure illustrates that neither the binary sequence nor any outliers (e.g.below the main sequence) have noticeably biased the fit.The magenta isochrone at the lower main-sequence represents the mean age of the Galactic disk (log τ Disk = 9.7).
The three groups in the top row are very young (< 25 Myrs), and low-mass stars have clearly not yet reached the true zero age main sequence.This illustrates well the power of low-mass stars for constraining very young population ages, for which the main sequence turn-off is sparsely sampled.The three groups in the bottom row are older.For them low-mass stars provide no two-sided age constraint.But for them the main-sequence turn-off is at lower masses, hence wellpopulation and providing the main age constraint.
While these fits are pleasing, they need to be validated against external results, which we do by comparing them to the ages from Cantat-Gaudin et al. (2020).We crossmatch our groups with known clusters in the Cantat-Gaudin et al. ( 2020) catalog.We keep all groups for which we obtain at least 20 crossmatches with one literature cluster.Additionally, for groups that have matching sources with more than one literature cluster, we keep the literature cluster with the highest number of -but at least 20 -crossmatches.This applies to 62 of our groups.Figure 6 compares our age determinations to those of Cantat-Gaudin et al. (2020), showing overall very good agreement, with no systematic bias and very little scatter.The log τ uncertainties for individual clusters in the Cantat-Gaudin et al. (2020) catalog range from 0.1 -0.25, depending on the age itself and the number of stars (Cantat-Gaudin et al. 2020).As described in more detail in Section 4.4, the log τ uncertainties of our isochrone fitting approach are typically smaller than 0.1, which is less than our log τ grid of 0.1.

Ages, Age Dispersions, [M/H] and Spatial
Structure of the Action-Angle Groups These fits, listed in Table 1 for the 92 action-angle groups, now put us in a position to address the basic questions we have about the properties of these actionangle groups.Table 1 includes the group number of each group, the number of member stars N * , the median action and angle coordinates, the best age estimate log τ , the intrinsic age spread estimate σ log τ , and the metallicity [Fe/H] of the likelihood isochrone fitting.Additionally, for the 62 groups we identified in the cluster catalog by Cantat-Gaudin et al. ( 2020) as described in Section 5.1, we list the name of the identified literature cluster.Among the listed groups we find established clusters such as the Pleiades (Melotte 22), the Hyades (Melotte 25), Praesepe (NGC 2632), and Coma Berenices (Melotte 111).We also identify the Meingast 1 (Pisces-Eridanus) stellar stream by crossmatching our sources with Meingast et al. (2019).However, a direct comparison and assignment of our groups to literature cluster is difficult, because each literature work uses a different approach to identify and select member stars of clusters, which leads to a lot of variation in the member selections (e.g.some of our action-angle groups are overlapping with the same literature cluster identified in Cantat-Gaudin et al. ( 2020)).A more detailed look at these different member selections additionally shows that our method selects more extended and loose structures (associations, tidal tails, other mechanisms of dispersal) as compared to Cantat-Gaudin et al. ( 2020) and Hunt & Reffert (2023), whose member selections focus on cluster cores mainly.
Table 1.This table provides an overview of the main parameters of the 92 action-angle groups.We provide the group number, the number of member stars N * , and the median action and angle coordinates of all groups.Additionally, we state the age log τ , the intrinsic age spread σ log τ , and the metallicity [F e/H] as determined by our likelihood isochrone fitting for each group.For the subset of 62 groups which we identified in the cluster catalog by Cantat-Gaudin et al. ( 2020), we list the literature name in the last column, as well as for the Meingast 1 (Pisces-Eridanus) stellar stream that we identified from Meingast et al. (2019).A slightly modified version of this "group center table" is available at the CDS. Figure 7 shows the distribution of ages (log τ ) among the 92 groups, which range from 20 Myrs to 2 Gyrs: all are much younger than the mean age of Galactic disk stars (∼ 3 − 5 Gyrs).For our purpose, the definitions of "young" vs. "old" revolve around dynamical timescales: one central timescale, the orbital period at the Solar radius is indicated by the blue dotted vertical line, for which we adopt logτ = 8.33.We find that half (43) of the groups are young in this sense, younger than one dynamical period.The 49 older groups have ages τ = 1 − 7 t orbit ≪ t Hubble .

Group
Our fitting also shows that almost 80% (72) of the groups seem to be "effectively mono-age", with σ log τ ≤ 0.2; only 20 groups show a modest age dispersion of σ log τ > 0.2.The value of 0.2 to distinguish these two regimes was chosen by comparing σ log τ to the uncertainty δ log τ of the isochrone fit, which is typically around 0.1 (see Section 4.4).The intrinsic age spread is indicated in the color coding of the bars in Figure 7, where groups with significant intrinsic age spread are shown in magenta, and mono-age groups are shown in turquoise.
This result immediately shows that none of these groups are just stars herded together by orbital resonances, as a wide intrinsic age spread would be expected in this case.
The best fitting metallicities fall almost exclusively on [F e/H] = 0.1, with only one group showing a metallicity of [F e/H] = 0.0.They are essentially all Solar [Fe/H], as expected for young stars near the Solar radius.Whether the slightly super-Solar values, surprising in light of the [Fe/H] of the present-day ISM at the Solar radius, reflect slight [M/H]-systematics in the isochrone fitting, remains open.Our metallicities are purely photometric, and draw in many cases heavily on the colors of lower main-sequence objects, whose isochrone-predicted colors are known to be uncertain.Therefore, our metallicities may well be subject to systematic errors at the 0.2 dex level.We consider this to be not relevant for our analysis, as the covariance between ages and metallicities are not a source of serious systematics in the age determinations.

Spatial Structure of the Groups
The spatial extent of the 92 groups (color-coded) is shown in Figure 8 in a top-down view onto the Milky Way.For a better overview, the left panel shows the 43 groups younger than one Galactic revolution (τ < 250 Mrys), and the right panel shows the 49 old groups with ages > 250 Myrs.This figure reveals a variety of spatial distributions and the internal spatial extents of the groups: some heterogeneous morphologies with little center-to-edge structure, and often an extent of ≳ 100 pc; some show a tight clustered "core" with or without tidal (or uncertainty-driven) tails.For the most part, these tails extend along the direction of the Galactic rotation as extended for tidal tails, as found and discussed previously in various papers (Meingast & Alves 2019;Röser et al. 2019;Röser & Schilbach 2019;Meingast et al. 2021;Fürnkranz et al. 2019;Tang et al. 2019).Overall, most groups extend across several hundred parsecs on the sky.

Overlap with Existing Cluster Catalogs
As was already apparent from the age comparison with Cantat-Gaudin et al. (2020), many of our groups  Kounkel & Covey (2019b), known to be complete by also contaminated (Zucker et al. 2022).
Among our 92 groups, 22 have no stars in common with any clusters in Cantat-Gaudin et al. (2020); the majority of our groups (51) have stars that overlap with exactly one cluster in Cantat-Gaudin et al. (2020).Remarkably, for 19 of our groups the stars within one group overlap with several (up to ten) clusters in the Cantat-Gaudin et al. ( 2020) catalog.This shows that our algorithm -designed to also identify loose orbit-space groups -can and does encompass several clusters if they are on near-identical orbits at very similar orbital phases.This is very much in line with the extended structure of our groups, already discussed in Section 5.2.2.

SUMMARY AND DISCUSSION
We have set out to devise and implement a new stringent two-step method to determine isochrone ages of action-angle groups in the nearby Galactic disk, even if they are sparse and widely dispersed.First, we overcome the limitation that the initial sets of stars, found to be clustered in 6D action-angle space, have a limited T ef f -range that hinders age dating with isochrones.We do this by projecting the 6D action-angle patch into the 5D space of position and proper motion, where the Gaia catalog is far more extensive.We can identify typically an order of magnitude more likely members.The second step is the actual isochrone fitting, as a function of both age and metallicity.This fitting accounts for a number of factors important here: the group members may be widely distributed in sky position, distance and reddening; they may have sizeable parallax uncertainties; and we must account for binaries and other "outliers", such as stars that are not group members.We do not only determine the mean age of each group, but we also test whether any group's CMD provides evidence for a significant intrinsic age spread.We find that our approach provides well-constrained ages, age dispersions and metallicities.
We have found that all of the identified 92 action-angle groups are "young" in a global Galaxy formation sense: τ ≪ τ disk ≈ 8 − 11 Gyrs, as they all have τ ≤ 2 Gyrs.About half of them are even younger than a dynamical period (< t orbit ≈ 250 Myrs).The large majority of them is consistent with a mono-age population, σ log τ ≤ 0.2 mag.
Not accounting for differential selection effects, we find that the distribution p log τ is broadly constant between log τ = 7 and 9, which would naively imply that p(τ ) ∝ τ −1 .
We also find that some groups are spatially extended and most likely entirely unbound, while others are spatially compact.
These findings fit together in a picture in which these action-angle groups indeed reflect ensembles of stars born at nearly the same time in nearly the same place, and are now dispersing to become field stars.How rapidly these groups of stars disperse, and how compact they once were, would require quantitative modelling of their actions and angles, which we will carry out in a forthcoming paper.
We have also found that we are not severely limited in the age determination by the group member sample size with our approach, as our algorithm provides almost the exact same age estimates for all groups if we randomly sample only half of their stars.This bodes well for applying the same approach in a forthcoming paper to the smallest and most diffuse groups we can identify, such as the pearls on a string found by Coronado et al. (2022) for a set of action angle groups: these are -often numerous -further clumps in angle (or orbital phase) space on nearly the same orbit as the initially identified action-angle group.

Figure 1 .
Figure 1.The member stars of the 94 most prominent action-angle groups identified with DBSCAN.The various colors represent the different groups.The black background corresponds to all sources in the Gaia RVS sample.Top/Middle row: The distribution of the groups in the three projections of the action/angle space.By construction, all groups are tightly clustered in this parameter space.Bottom left: Top-down view on the 94 groups in Galactic Cartesian coordinates, where the Sun is located at (X,Y) = (0,0) pc.Bottom right: The sky distribution of the action-angle groups in Galactic longitude and latitude, illustrated in a Hammer projection.Although all groups are very tightly clustered in action-angle space, most of them are extended across several hundred degrees on the sky.

Figure 2 .
Figure2.Observational CMDs for six representative action-angle groups.The member stars of each group are arranged in narrow, but short, sequences, where stars are clearly missing in the upper and lower main-sequence parts.This is caused by the limited magnitude and temperature range of the Gaia RVS sample, where only stars with G < 14 mag and with effective temperatures in the range of 3.100 to 14.500 K are included.In color, we show a set of three isochrones(log τ = 7.8, 8.4, 8.8), illustrating that stars located on the main sequence only have a small age-diagnosing value.

Figure 3 .
Figure 3. Histograms of the allocated density values for all Gaia (turquoise) and Gaia mock stars (violet) within the individual (X, Y, Z, µα, µ δ ) ranges for six representative action-angle groups.The dotted blue line marks the location of the final density cut T f inal as described in Subsection 3.3.Most Gaia and Gaia mock stars have a density value close to zero, with a tail extending to higher densities, where typically the amount of Gaia stars exceeds the amount of Gaia mock stars.This is an indication that we find stars clumped in position-proper motion space in the Gaia sample.

Figure 4 .
Figure4.Left: The effect of the extinction correction method for stars in Group 68.The correction reveals a significant shift of the uncorrected (colored) stars to the reddening corrected (black) stars, emphazising the importance of extinction correction for precice age dating via isochrone fitting of CMD's, even for stars located in the solar neighborhood.The average shift in color for Group 68 is 0.1 mag.Right: Illustration of the likelihood isochrone fitting for Group 68.All isochrones are colored by age (107 different logarithmic τ ) for a fixed metallicity of [Fe/H] = 0.1.The color coding represents the χ 2 value obtained for the different age setups.The plot illustrates that very young as well as very old ages are not likely to fit the group's distribution, as they lead to a high χ 2 value.The best fit isochrone for this group is the one with the lowest χ 2 value.The stars in magenta mark group members that are labelled as 4-σ outliers by our code setup, and are not considered in the likelihood isochrone fitting.We show the isochrones not as solid lines, but as datapoints to illustrate our finely spaced magnitude and color setup, which we obtain by linear interpolation (see Sec. 2.3).

Figure 5 .
Figure 5. Observational CMD's for six representative groups, with additional member stars associated from the 5D positionproper motion space.All stars are individually extinction corrected.The turquoise curve on top of each sequence shows the best-fitting isochrone (log τ , [Fe/H]).Across all groups, the isochrones fit both the upper and lower main sequences very well.The magenta isochrone at the lower main-sequence represents the age of the Galactic disk.The three groups in the top row are very young (< 25 Myrs), and their lower main sequences are clearly offset from the age of the Galactic disk.This again illustrates the power of the lower main sequence to determine the ages of very young populations.The three groups in the bottom row are older, and for them the main-sequence turn-off is the main age determinant.

Figure 6 .
Figure 6.Comparison between the age determination of this work with the age determination of Cantat-Gaudin et al. (2020) for 62 established clusters that we could identify in the cluster catalog by Cantat-Gaudin et al. (2020).The plot shows that the age estimates from this work are in very good agreement with the literature values, with only very little scatter The log τ uncertainties of the age fit range from 0.1 -0.25 for Cantat-Gaudin et al. (2020), and are typically smaller than 0.1 for the method presented in this work.

Figure 7 .
Figure 7. Histograms of the distribution of age (log τ ) across the set of 92 groups.About half of them (43) are younger than an orbital period, and 49 of them are older.The central timescale, the orbital period at the Solar radius is indicated by the blue dotted vertical line, for which we adopt logτ = 8.33.The intrinsic age spread is indicated in the colorcoding of the bars, where groups with an intrinsic age spread are shown in magenta, and mono-age groups are illustrated in turquoise.

Figure 8 .
Figure 8. Top-down views on the 92 action-angle groups in Galactic Cartesian coordinates, where the Sun is located at (X,Y) = (0,0).Left: 43 young groups with ages < 250 Myrs.Right: 49 old groups with ages > 250 Myrs.Their morphologies are either heterogeneous, or show a tight clustered "core" with or without tidal (or uncertainty-driven) tails.are in published catalogs.Specifically, 33% of stars that are in one of our groups have a crossmatch with the Cantat-Gaudin et al. (2020) cluster catalog, 52% with the Hunt & Reffert (2023) catalog, and 49% with the catalog of Kounkel & Covey (2019b), known to be complete by also contaminated (Zucker et al. 2022).Among our 92 groups, 22 have no stars in common with any clusters in Cantat-Gaudin et al. (2020); the majority of our groups (51) have stars that overlap with exactly one cluster in Cantat-Gaudin et al. (2020).Remarkably, for 19 of our groups the stars within one group overlap with several (up to ten) clusters in the Cantat-Gaudin et al. (2020) catalog.This shows that our algorithm -designed to also identify loose orbit-space groups -can and does encompass several clusters if they are on near-identical orbits at very similar orbital phases.This is very much in line with the extended structure of our groups, already discussed in Section 5.2.2.

Table 1 continued
Table 1 (continued) How old are the groups?Are they mono-age?