Young, wild, and free: The early expansion of star clusters

Early expansion plays a fundamental role in the dynamical evolution of young star clusters. However, until very recently most of our understanding of cluster expansion was based only on indirect evidence or on statistically limited samples of clusters. Here we present a comprehensive kinematic analysis of virtually all known young Galactic clusters ( t < 300Myr) based on the improved astrometric quality of the Gaia DR3 data. Such a large sample provides an unprecedented opportunity to robustly constrain the fraction of clusters and the timescale during which expansion has a prominent impact on the overall kinematics. We ﬁnd that a remarkable fraction (up to 80%) of clusters younger than ∼ 30Myr is currently experiencing signiﬁcant expansion, whereas older systems are mostly compatible with equilibrium conﬁgurations. We observe a trend in which the expansion speed increases with the cluster-centric distance, suggesting that clusters undergoing expansion will likely lose a fraction of their present-day mass. Also, most young expanding clusters are extended, possibly due to the expansion itself. A comparison with a set of N -body simulations of young star clusters shows that the observed expansion pattern is in general qualitative agreement with that found for systems undergoing violent relaxation and evolving toward a ﬁnal virial equilibrium state. However, we also note that additional processes likely associated with residual gas expulsion and mass loss due to stellar evolution likely also play a key role in driving


Introduction
It is commonly accepted that most Galactic star formation takes place in embedded proto-clusters as a result of turbulent, clumpy, and stochastic processes in molecular clouds (e.g.Lada & Lada 2003;Offner et al. 2009;Feigelson et al. 2013).The environments in which stars form determine to some extent a number of key properties characterizing the stellar clusters we observe today, such as their initial mass function, the stellar multiplicity, and the star-by-star chemical abundance differences routinely observed in massive stellar clusters (e.g.Bastian & Lardo 2018;Gratton et al. 2019;Milone & Marino 2022).
The past years have witnessed a renewed interest in the theoretical and observational study of the formation of stellar clusters.Particular attention has been paid to the very early evolutionary phases, when clusters are expected to undergo violent relaxation and to be evolving toward a virial equilibrium state.Studies of clusters' formation and very early evolution have shown that clusters can emerge from these evolutionary phases with internal rotation, radial anisotropy in the velocity distribution, and mass segregation.However, the detailed role of the various physical processes involved and the variety of different dynamical paths are still a matter of intense investigation (see e.g.McMillan et al. 2007;Allison et al. 2009;Moeckel & Bonnell 2009;Pfalzner & Kaczmarek 2013;Banerjee & Kroupa 2014;Vesperini et al. 2014;Fujii & Portegies Zwart 2016;Parker et al. 2016;Domínguez et al. 2017;Mapelli 2017;Daffern-Powell & Parker 2020;González-Samaniego & Vazquez-Semadeni 2020).
⋆ e-mail: alessandro.dellacroce@inaf.it The clusters' early evolutionary phases have been poorly constrained also by observations, mainly due to the difficulty of obtaining reliable kinematic information for statistically significant samples of cluster members.Data from the Gaia space mission and, in particular the DR3 (Gaia Collaboration et al. 2023), provided the deepest and most precise astrometric catalog ever obtained and opened the possibility of new studies about star cluster formation and dynamics.Indeed, Gaia's unprecedented kinematical mapping of the Galaxy and of its stellar components significantly enriched our knowledge about the number and the properties of star clusters in the Milky Way (see e.g.Cantat-Gaudin et al. 2020;Cantat-Gaudin & Anders 2020).Also, it enabled detailed studies of nearby star-forming regions and their use as ideal laboratories to shed new light on cluster formation and early evolution (e.g.Beccari et al. 2018;Kuhn et al. 2019;Meingast & Alves 2019;Dalessandro et al. 2021;Della Croce et al. 2023).
The evolution of a young cluster depends on many physical factors, such as the cluster's structural properties (e.g.mass and size), the strength of the external tidal fields, and their physical and dynamical ages (see e.g. Hills 1980;Parker et al. 2014;Sills et al. 2018).Early cluster expansion is expected to play a fundamental role in cluster evolution and to drive stellar systems toward complete dissolution (Elmegreen 1983;Mathieu 1983;Adams 2000;Kroupa 2001;Goodwin & Bastian 2006;Pelupessy & Portegies Zwart 2012;Dinnbier & Kroupa 2020;Dinnbier et al. 2022;Farias & Tan 2023).Theoretical studies show that cluster expansion due to processes like gas expulsion and the subsequent violent relaxation phase may last several Myr, depending on many cluster properties such as the cluster dynamical state at the expulsion, the specific star-formation efficiency, and cluster size, density, and mass (Kroupa et al. 2001;Baumgardt & Kroupa 2007;Pelupessy & Portegies Zwart 2012;Banerjee & Kroupa 2013;Pfalzner & Kaczmarek 2013;Brinkmann et al. 2017;Farias et al. 2017;Farias & Tan 2023;Li et al. 2019;Pang et al. 2021;Leveque et al. 2022).Eventually, the system can either settle into equilibrium, after losing some mass during expansion or dissolve completely.The subsequent long-term evolution of surviving clusters will be driven by two-body relaxation and the interplay between the effects of internal evolution and the external Galactic tidal field (see e.g. de La Fuente Marcos 1997; Hurley et al. 2005).
In particular, Kuhn et al. (2019) studied the kinematical properties of a sample of 28 young (1 − 5 Myr) clusters and associations by using Gaia (DR2) proper motions.According to those authors, observations are consistent with early cluster expansion driven by changes in the gravitational potential due to the dispersal of the molecular cloud.More recently, Guilherme-Garcia et al. ( 2023) investigated the kinematics of 1237 clusters using a technique that aims at reconstructing the underlying velocity field.They found 8 (and an additional 9 candidates) clusters displaying rotation patterns whereas 14 (15 more candidates) clusters show evidence of expansion in their velocity fields (Guilherme-Garcia et al. 2023).The vast majority of the expanding systems in their sample are younger than 100 Myr, with a peak around 10 Myr.While these studies provided important clues about the early evolution and survival of young clusters, the availability of more precise data from the Gaia DR3 and the possibility to assess cluster membership more robustly can allow us to disentangle critical aspects related to their evolution.
In this work, we present a comprehensive and systematic analysis of the kinematical properties of (virtually) all young (t < 300 Myr) stellar clusters identified so far in the Milky Way (using the cluster catalog by Cantat-Gaudin & Anders 2020; Cantat-Gaudin et al. 2020), with particular focus on the early cluster expansion phase.Such a large sample provides the unprecedented opportunity to robustly constrain the timescale during which expansion has a prominent impact on the overall cluster kinematics.Also, it allows us to trace how expansion affects or is linked to the cluster properties and formation mechanisms.
The paper is organized as follows.In Section 2 the cluster membership determinations and kinematic analyses are described.We present our results in Section 3 and we compare them with numerical simulations in Section 4. We assess the impact of different age estimates on our results in Section 5 and a detailed comparison with previous works is presented in Section 6.Finally, conclusions are drawn in Section 7.

Clusters' membership
The present work makes use of the list of clusters identified by Cantat-Gaudin et al. (2020) using Gaia DR2 data.We focused on systems younger than ≤ 300 Myr (according to age estimates by Cantat-Gaudin et al. 2020) to investigate the evolution of star clusters in their very early stages.In this way, we selected 1179 clusters out of 2017 in the original catalog.
To take full advantage of the most recent and accurate Gaia DR3 data release for the definition of cluster member stars, we performed an independent membership analysis.For every cluster, we retrieved Gaia DR3 data for sources brighter than G = 18, and that have a 5-parameter solution (i.e. with sky position, proper motions, and parallax measurements).In particular, every query was centered on the cluster's centroid (as reported by Cantat-Gaudin et al. 2020) and the search radius was defined as R search, sky ≡ 2R 95, sky , where R 95, sky is the radius enclosing 95% of the member stars reported by Cantat-Gaudin et al. (2020).
For each cluster, we then selected stars according to their motion with respect to the cluster's bulk motion.Only sources within R search, PM ≡ 2 R 95, PM , with R 95, PM being the circle in proper motion space enclosing 95% of the members were retained in the subsequent analysis.We did not apply any preliminary parallax selection.These selections allowed us to include all the previously listed members in our starting Gaia DR3 catalog.
We performed the clustering analysis in the five-dimensional space of Galactic coordinates, proper motions, and parallax (i.e.ℓ, b, µ α * , µ δ , and ϖ).Since we are dealing with heterogeneous quantities, we preliminary scaled them all, so that the mean and standard deviation of their distributions are equal to zero and one respectively.We used the StandardScaler provided by the Python library sklearn.We then performed an unsupervised clustering analysis on the scaled coordinates by means of the HDBSCAN (McInnes et al. 2017;McInnes & Healy 2017) algorithm.
Based on the results of several tests, we set the algorithm's parameters as follows.min_cluster_size = 20 for N CG20 < 100 , and where N CG20 is the number of member stars reported by Cantat-Gaudin et al. (2020).These parameters appeared to be best suited for the unsupervised search for members in star clusters characterized by significantly different extensions (R search, sky ), velocity distributions (R search, PM ) and number of likely members (N CG20 ).As a test case, in Fig. 1 we show the spatial, parallax, proper motion, and color-magnitude diagram (CMD) distributions for the selected members of ASCC 113 As expected, selected member stars are clustered in the five-dimensional astrometric space (i.e.ℓ; b; ϖ; µ α * ; µ δ ).Also, they exhibit a well-defined, clusterlike sequence in the CMD, thereby confirming the ability of our clustering analysis to recover stellar cluster members among field stars.
While running the procedure for all the clusters in the sample, we found that the algorithm performed poorly on clusters with only a few member stars (N CG20 < 30).This is probably due Article number, page 2 of 10 80.0 82.5 85.0 N (Cantat-Gaudin+20) to the choice of min_cluster_size, which sets a lower limit to the number of member stars in our search for star clusters.However, since the kinematic analysis strongly benefits from the availability of relatively large samples of individual velocities, we decided to exclude these clusters from our analysis.In this way, we ended up with a final sample of 949 clusters for which the clustering analysis was performed.
Fig. 2 shows a direct comparison between the number of member stars obtained by using Gaia DR3 and that found by Cantat-Gaudin et al. (2020).While some differences are expected due to the different adopted clustering algorithms and Gaia data releases, Fig. 2 shows an overall agreement between the two compilations.

Cluster kinematics analysis
First, we accounted for perspective effects induced by the cluster bulk motion following the equations reported by van Leeuwen (2009).To this aim, systemic line-of-sight (LOS) velocities from Tarricq et al. (2021) were used.Only 509 out of 949 clusters (see Section 2.1) had reported LOS velocity and were thus retained in the subsequent analysis.In Figure 3  tributions within 300 Myr of the original catalog and the final sample of clusters.The distributions populate the investigated age range in a similar fashion and our final sample of 509 clusters is representative of the initial age distribution.
We then used stars with a membership probability larger than 70% based on the clustering analysis performed in this work (Section 2.1).We also selected cluster members having ruwe ≤ 1.4, astrometric_gof_al ≤ 1, and astrometric_excess_noise ≤ 1 mas (if Article number, page 3 of 10 astrometric_excess_noise_sig > 2), thus excluding stars for which the standard five-parameter solution did not provide a reliable fit of the observed data (Lindegren et al. 2021).
We inferred the mean radial velocity, ⟨v R ⟩, and the radial velocity dispersion, σ R , in a fully Bayesian framework properly accounting for errors on individual velocities.We explored the parameters space employing a Markov Chain Monte Carlo (MCMC) technique.In particular, we used the Python package emcee (Foreman-Mackey et al. 2013).For each system, we assumed the likelihood (Pryor & Meylan 1993) with v R,k and δv R,k being the radial velocity and the respective error of the k-th member.Eq. 3 assumes that the intrinsic distribution along the radial component of the velocity is a Gaussian with mean ⟨v R ⟩ and velocity dispersion σ R .We used uniform priors within [−10; +10] mas yr −1 and [0.001; 15] mas yr −1 for ⟨v R ⟩ and σ R respectively.For each cluster, we initialized 50 walkers and ran the algorithm for 500 steps, which was found to be sufficient for both ensuring convergence (for which the first half of samples was discarded) and accounting for correlations between samples (typically of the order of 20).Median values and 16% and 84% quantiles (corresponding to the 1σ interval if the distributions were Gaussian) were then computed for each quantity directly from posterior samples.The kinematical analysis was performed by adopting the clusters' geometric center, defined as the median of the positions of member stars.

Results
Fig. 4 shows the ratio between the mean radial velocity and the radial velocity dispersion (hereafter ⟨v R ⟩/σ R ) as a function of clusters' ages from Cantat-Gaudin et al. (2020).This quantity provides an indication of the amplitude of the ordered to the disordered motion of stars along the radial component, thus directly tracing ongoing expansion or contraction.A positive value of ⟨v R ⟩/σ R implies expansion.
The first key result highlighted by Fig. 4 is that about 80% of clusters younger than ∼ 30 Myr show positive ⟨v R ⟩/σ R values.More in general, the total fraction of young (< 30 Myr) systems having positive ⟨v R ⟩/σ R at the 3σ level is 43% (see Table 2 for the fraction of expanding systems in different age bins).The distribution attains maximum values ⟨v R ⟩/σ R = 1.5 − 2 (such clusters have typically ⟨v R ⟩ ∼ 2 km s −1 ) for clusters with age ∼ 10 Myr.Then it progressively decreases.For clusters older than ∼ 30 Myr, the distribution of ⟨v R ⟩/σ R flattens around 0 and it shows an intrinsic standard deviation of about 0.2 (corresponding to ⟨v R ⟩ ≲ 0.1 km s −1 ).Such a clear trend allows us to identify for the first time the timescale (of about 30 Myr) during which expansion plays a significant role in the overall cluster kinematics.We verified that these results do not change significantly if different age estimates were adopted (see Section 5).
In the left panel of Fig. 4 clusters are color-coded according to the number of member stars.We note that the amplitude of the scatter around zero observed in systems older than 30 Myr strongly depends on the number of members identified.Results of numerical experiments of equilibrium stellar clusters (see Section 4.2) suggest that the observed spread around zero can be largely explained in terms of statistical fluctuations due to a low number of stars (grey shaded area in Fig. 4).This in turn confirms that the distribution observed for older clusters is consistent with what is expected for systems in equilibrium.
Most of the young (< 30 Myr) clusters with clear ongoing expansion are characterized by some of the largest values (up to ∼ 10 pc) of R 50 (Fig. 4, middle panel).We note that among the expanding systems, those with smaller R 50 values preferentially attain smaller ⟨v R ⟩/σ R by up to a factor of two.Finally, old (> 30 Myr) expanding clusters have mostly fewer members and larger extensions suggesting that they might have been expanding for several tens of Myr.
Furthermore, we performed a linear fit to the distribution of v R as a function of the cluster-centric distance for each cluster and we derived the angular coefficient m and its relative error σ m .Interestingly, we found that the majority of the expanding clusters show a positive and significant (m/σ m ≥ 3) ranking in their v R distributions as a function of their positions in the cluster (see the right panel in Fig. 4).These patterns could suggest that young expanding clusters are likely losing a fraction of their original mass as stars in the outskirts become unbound.However, we note in passing that this does not necessarily imply that all the expanding clusters will become unbound.In fact, stars in the external regions moving away from the cluster can produce a positive slope even if the inner parts are not expanding.
Finally, we present the kinematic properties of a few prototypical systems in Figure 5.We selected two young (< 30 Myr) clusters, namely NGC 6193, and NGC 4103, and two older systems, LP 2219, and NGC 3114.Within each group, we picked one system showing significant evidence of expansion (NGC 6193, and LP 2219), whereas the other one is compatible with equilibrium (NGC 4103, and NGC 3114).Relevent cluster properties are reported in Table 1 for reference.Notes.From left to right, cluster name, the ratio between the mean radial velocity and the radial velocity dispersion, number of members, cluster age (according to Cantat-Gaudin et al. 2020), and radius enclosing half of the members.Errors on ⟨v R ⟩/σ R are directly obtained from the MCMC sampling of the posterior distribution, whereas errors on R 50 are obtained by bootstrap resampling the radial distribution of stellar members.The values reported correspond to the 16th and 84th percentiles of the distributions (i.e.corresponding to the 1σ value if the distributions were Gaussian).
In particular, the left panels of Fig. 5 show the spatial distribution of members used to compute ⟨v R ⟩/σ R (see Fig. 4), as well as the velocity vectors on the plane of the sky.Arrow lengths are proportional to their total speeds whereas the color coding traces the amplitude of the radial velocity component (Fig. 5).Right panels, on the other hand, show the distribution of individual velocities as a function of the cluster-centric distance.The linear regression used to obtain the m/σ m parameter (see Fig. 4) is also shown and the value of the slope is reported.To ease the comparison, spatial coordinates are scaled to R 50 for each cluster, and their velocities are shown in the same velocity scale.Expanding clusters (NGC 6193, and LP 2219) show a preferential alignment of velocity vectors along the positive (i.e.pointing outward) radial direction, as highlighted by the color coding (Fig. 5).This feature is particularly evident in the external regions.A clear trend is also observed in v R vs R distribution, and it is reflected in the positive value of the slope from the linear regression (as already pointed out in Fig. 4).
Non-expanding clusters (NGC 4103, and NGC 3114), on the other hand, present radial velocities that are scattered around zero, without any preferential alignment along the radial direction or significant radial trend.All these features are consistent with the systems being in equilibrium.
Lastly, we looked for any dependence of the expansion properties on the cluster masses.In particular, we crossmatched our cluster catalog with the recent compilation of cluster masses provided by Almeida et al. (2023).We found 227 clusters in common.For those clusters, we show ⟨v R ⟩/σ R as a function of cluster age, color-coded according to the mass reported by Almeida et al. (2023, see Fig. 6).Qualitatively, expansion affects clusters of any mass, from about 10 2 M ⊙ to a few 10 3 M ⊙ .To elaborate more on this point, we split the population of clusters younger than 30 Myr (45 clusters) into two sub-samples: expanding (27 clusters, about 60% of the sample), and non-expanding clusters (18 out of 45 systems, i.e. the 40%).For the purpose of this analysis, we classified a cluster to be expanding if the expansion signal is significant at the 1σ level.The only purpose here is to obtain two almost equally populated samples of clusters.We then compared the mass distribution of the two sub-populations with each other and with the full sample of young clusters.A Kolmogorov-Smirnov test shows that there is no significant statistical difference between the three populations.Although based on a few dozen clusters, this result suggests that the physical processes driving the expansion of star clusters in their early stages of formation and evolution are effective irrespective of the cluster mass.

N-body simulations of cluster formation
In this section, we present a brief analysis of a few N-body simulations of young star clusters undergoing the violent relaxation phase and evolving toward their final virial equilibrium state.
The goal here is to illustrate the time evolution of their global expansion pattern and establish a general connection with the observational results presented in the previous sections.The simulations considered here are part of the suite discussed in detail in Livernois et al. (2021) who explored the early evolution of systems starting with the homogeneous or fractal spatial distribution and investigated the role of initial rotation on the cluster's early evolutionary phases.Here we focus on four models: two models with homogeneous initial spatial distributions, one without and one with rotation (hereafter referred to as H0 and H075 respectively; see Livernois et al. 2021, for further details about the initial conditions) and two models with an initial fractal spatial distribution without and with the initial rotation (hereafter F0 and F025).Fig. 7 shows the time evolution of ⟨v R ⟩/σ R for the selected models as obtained by using all stars in the system within the tidal radius.All models show an initial contracting phase followed by significant expansion.The initially homogeneous models display more rapid and extreme collapse and expansion phases than the initially fractal models as the clumps within the latter models merge and interact with other clumps before arriving at the center of the system.Although these specific models do not reach the extreme values of ⟨v R ⟩/σ R found in our study (Fig. 4), the range of expansion values found in these simulations spans those attained by most of the observed clusters, thus suggesting that violent relaxation can play a key role in triggering and driving early cluster expansion.We emphasize that the idealized models presented here are included just to illustrate the general kinematic behavior during these early evolutionary phases and they are not meant to provide a detailed fit to the observational data.Different initial conditions and more realistic simulations including additional processes such as gas expulsion and mass loss due to stellar evolution might be necessary to reach the most extreme values found in the observational sample and to constrain the range of timescales of the early ex-pansion, of the settling to the final equilibrium as well as the timescale associated with the possible cluster's dissolution.We note also that the theoretical lines shown in Fig. 7 represent the evolutionary path of the ⟨v R ⟩/σ R ratio of clusters surviving the early evolutionary phases.We point out that it is likely that not all the clusters with age < 30 Myr in the observational sample will follow this path as some of them will continue expanding and will eventually dissolve.
Fig. 7. Time evolution of the ratio between the mean radial velocity and the radial velocity dispersions for a selection of models from Livernois et al. ( 2021) for all stars within the tidal radius.

Distribution of ⟨v R ⟩/σ R for star clusters in equilibrium
Results presented in Fig. 4 suggest that for t > 30 Myr the distribution of the ⟨v R ⟩/σ R is compatible with what is expected for clusters in equilibrium and that the broadening of the distribution might be mostly driven by statistical fluctuations.
To check whether this is indeed the case, we have created 100 random realizations of a population of clusters.Each population contains the same number of clusters as in our observed sample and in each realization the clusters have the same number of stars as in the observed sample.Positions and velocity of stars in each cluster follow those of a King model with central dimensionless potential W 0 = 5 (but the results do not have any significant dependence on the particular model adopted).For each realization, we have then calculated the dispersion in the distribution of the values of ⟨v R ⟩/σ R of the clusters in the sample.The average value of the dispersion found in the 100 realizations is equal to about 0.11 and is shown as a grey-shaded area in Fig. 4. The comparison demonstrates that 1) the spread observed for older systems can be largely accounted for by statistical fluctuations and 2) equilibrium models cannot account for large positive ⟨v R ⟩/σ R (i.e.expansion) observed in clusters with t < 30 Myr and that those systems are therefore out of equilibrium.
We note, however, that the observed spread is larger (about a factor 2) than that derived from the simulations.While these residuals might suggest that some of the clusters older than ∼ 30 Myr are still oscillating around equilibrium configurations, it is important to emphasize that in the analysis of the ⟨v R ⟩/σ R from the numerical realizations of cluster populations, not all the possible effects that can determine the spread of this quantity are in-Article number, page 7 of 10 A&A proofs: manuscript no.main_article cluded.They do not account, for example, for a number of possible uncertainties associated with the observational data such as the uncertainties induced by wrong LOS systemic velocities (mainly due to low-number statistics).Further investigation is needed to properly study this issue.

Testing the impact of different age estimates on the expansion pattern
The reference age compilation adopted in this analysis is the one by Cantat-Gaudin et al. (2020).Here we test the robustness of our results against different compilations of clusters ages from literature.In particular, we used cluster ages obtained by Kharchenko et al. (2013), Bossini et al. (2019), andDias et al. (2021) who determined ages for several star clusters in our Galaxy through isochrone fitting.We note that Kharchenko et al. (2013) exploited pre-main-sequence stars as a further age indicator to obtain more reliable ages in the young end.Hunt & Reffert (2023), performed an all-sky search for open clusters and they provided ages for every cluster in their sample obtained by means of a convolutional neural network.We crossmatched the sample of 509 clusters adopted in this study with these catalogs.
Figure 8 shows the distribution of ⟨v R ⟩/σ R as a function of age (as in Figure 4) for the clusters in common with the four compilations.Also, in Table 2 we report the fractions (as well as the total numbers) of clusters exhibiting evidence (at the 3σ level) for expansion in different age bins.Although different catalogs span different age ranges for the same clusters, the trend of expanding clusters for ages below ≃ 30 Myr is clearly visible in all the catalogs.This test demonstrates that the results presented in the manuscript are robust against different age estimates, also in terms of time scale during which expansion has an important role in cluster kinematics.Kuhn et al. (2019) investigated the expansion properties of twenty-eight stellar clusters and associations, concluding that at least 75% of clusters in their sample are expanding.Among the clusters in common, we find a significant expansion for NGC 1893, NGC 2244, and NGC 2362 and a mild expansion for IC 348, and NGC 6231, consistent with the results by Kuhn et al. (2019) for these systems.
Besides dedicated studies (such as Kuhn et al. 2019;Guilherme-Garcia et al. 2023), evidence of expansion was found in several works.For instance, Bravi et al. (2018) investigated the kinematical properties of four young (≳ 30 Myr) open clusters, namely IC 2602, IC 2391, IC 4665, and NGC 2547.They found that all the clusters but IC 4665 (for which they only put an upper limit) are super virial, concluding that this is consistent with the residual gas expulsion scenario.Clusters indeed expand eventually returning to an equilibrium state after unbound stars are dispersed.This process might take several tens of system crossing times (Baumgardt & Kroupa 2007).Consistently, we found that all these four clusters show slow expansion speeds (typically ⟨v R ⟩/σ R ≲ 0.3 ± 0.1), suggesting they might be close to equilibrium.Lim et al. (2020) studied the star-forming region W4, where the cluster IC 1805 is located.They found that the cluster is composed of an isotropic core and an external region showing clear evidence of expansion.These features suggest that the cluster is experiencing expansion after an early, initial collapsing phase (Lim et al. 2020).We also found a strong indication of expansion in IC 1805, although we did not find evidence for a central isotropic core in the cluster.
The cluster NGC 2244 lies at the center of the Rosette Nebula.The region shows a complex interplay between stellar and gas kinematics, and feedback-driven star formation in substructured environments.Both rotation and expansion were found in NGC 2244 (Lim et al. 2021).In our study we also found NGC 2244 to be significantly expanding (⟨v R ⟩/σ R ≃ +0.98 +0.14  −0.13 ).Pang et al. (2021) investigated the connection between the cluster's internal kinematics and their morphology.They found younger clusters to exhibit filament-like substructures while older ones show tidal-tail features.The majority of the systems are inferred gravitationally unbound and expanding (Pang et al. 2021).In the present study, mild expansion has been directly detected in NGC 2422, NGC 2451a, NGC 2451b, NGC 2232.On the other hand, neither expansion nor contraction has been found in NGC 2516, suggested to be in a super virial state (Pang et al. 2021).Investigating internal kinematic patterns, or dependencies on cluster morphology for each cluster is beyond the scope of this work.We defer this topic to a follow-up study.
In summary, the excellent agreement with previous, detailed studies that focussed on a handful of systems suggests that our results are solid and that we are effectively probing the expansion of young stellar clusters.

Discussion and Conclusion
We performed a comprehensive analysis of the internal kinematics of young star clusters (t < 300 Myr) in the Milky Way with the aim of reconstructing the key properties and possible physical mechanisms shaping the early cluster expansion.We emphasize that this analysis is based on a sample 20 times larger and spanning a cluster age range 60 times larger than previous analyses (see e.g.Kuhn et al. 2019), thus enabling for the first time the possibility of constraining the timescale during which expansion has a dominant impact on cluster kinematics and the fraction of stellar systems significantly affected by the expansion.
Our analysis reveals a clear trend in which the fraction of expanding clusters increases for younger clusters (see Table 2): a significant fraction of clusters younger than ∼30 Myr are characterized by a significant expansion attaining values as large as ison of the observed kinematic patterns and the timescales associated with the early evolution.Finally, we note that extremely young clusters (with ages smaller than a few Myr) not included in our sample would be necessary to probe the kinematic patterns associated with the systems' very early dynamics.The lack of these systems is likely a selection effect, as these very young systems are probably embedded clusters and could hardly be observed with Gaia.In this respect, future data from infra-red surveys would provide relevant insights into the embedded cluster population (see for example the VISIONS survey, Meingast et al. 2023) and allow us to build a more complete dynamical picture of the evolution of these systems.

Fig. 1 .
Fig. 1.Astrometric and photometric properties of ASCC 113 members (shown in blue), whereas in gray we show the properties of the initial sample of Gaia sources (see Section 2.1).In the left panel are the distributions in Galactic coordinates and parallax (top-right corner), while the middle panel shows the PM distribution.Finally, the color-magnitude diagram is shown in the right panel.The parallax distributions were scaled for visualization purposes only.

Fig. 2 .
Fig. 2. Comparison of member stars between this work and Cantat-Gaudin et al. (2020, N GC+20 ), for all the clusters included in the subsequent analysis.Diagonal lines show the one-to-one relation and the relations corresponding to half and twice the members (as also reported in the plot).

Fig. 3 .
Fig. 3. Age distributions for clusters in the starting catalog (in red, Cantat-Gaudin et al. 2020) and those retained in the kinematic analysis (in blue).Histograms were normalized such that their areas sum to unity.

Fig. 4 .
Fig. 4. The ratio between the mean radial velocity and the radial velocity dispersion for clusters younger than 300 Myr (according toCantat- Gaudin et al. 2020).Errors on the y-axes were obtained directly from the MCMC samples, while ages are from Cantat-Gaudin et al. (2020).In the left panel, clusters are color-coded according to the number of members.Also, the standard deviation (and twice the value) of the ⟨v R ⟩/σ R ratio obtained from numerical realizations of equilibrium star clusters (see Section 4.2) are shown as gray shaded areas.Values obtained from numerical realizations were convolved with the median observational error to allow for a direct comparison with the underlying data.In the middle panel colors depict R 50 , whereas in the right panel, the color coding represents the ratio between the slope and the corresponding error from the linear regression of radial velocities as a function of cluster-centric distance.

Fig. 5 .Fig. 6 .
Fig. 5. Kinematic properties of NGC 6139, NGC 4103, LP 2219, and NGC 3114 (from top to bottom).The left panels show the spatial distribution of members in Cartesian coordinates normalized to the radius enclosing half of the members, with arrows showing the velocity vectors on the plane of the sky.The arrow lengths are proportional to the speed on the plane of the sky (velocity scale reported in the top-right corners), while their colors map the radial component (v R ) of the velocity.Positive values point outward.The right panels show the distribution of members in the v R − R plane.The linear regressions of velocities as a function of the cluster-centric distances are shown (orange lines), as well as the posterior values on the slopes at the bottom.Article number, page 6 of 10

Table 1 .
The main properties of the four clusters whose kinematic features are shown in Figure5.