The Gaia-ESO Survey: 3D dynamics of young groups and clusters from GES and Gaia EDR3

We present the first large-scale 3D kinematic study of ~2000 spectroscopically-confirmed young stars (<20 Myr) in 18 star clusters and OB associations (hereafter groups) from the combination of Gaia astrometry and Gaia-ESO Survey spectroscopy. We measure 3D velocity dispersions for all groups, which range from 0.61 to 7.4 km/s (1D velocity dispersions of 0.35 to 4.3 km/s). We find the majority of groups have anisotropic velocity dispersions, suggesting they are not dynamically relaxed. From the 3D velocity dispersions, measured radii and estimates of total mass we estimate the virial state and find that all systems are super-virial when only the stellar mass is considered, but that some systems are sub-virial when the mass of the molecular cloud is taken into account. We observe an approximately linear correlation between the 3D velocity dispersion and the group mass, which would imply that the virial state of groups scales as the square root of the group mass. However, we do not observe a strong correlation between virial state and group mass. In agreement with their virial state we find that nearly all of the groups studied are in the process of expanding and that the expansion is anisotropic, implying that groups were not spherical prior to expansion. One group, Rho Oph, is found to be contracting and in a sub-virial state (when the mass of the surrounding molecular cloud is considered). This work provides a glimpse of the potential of the combination of Gaia and data from the next generation of spectroscopic surveys.


INTRODUCTION
Stars form in turbulent molecular clouds with a hierarchical and highly substructured spatial distribution that follows the distribution of dense gas (Elmegreen 2002).Many recently-formed stars are observed in groups or 'clusters' of some sort, some showing substructure, while others are centrally concentrated with a smooth density distribution (Lada & Lada 2003;Gutermuth et al. 2008).The majority of groups quickly disperse, with very few young groups surviving as long-lived open clusters (Adams 2010).The stellar groups that disperse are often briefly observed as 'associ-ations', gravitationally unbound and low-density groups of young stars (Wright 2020;Wright et al. 2023).
Young star clusters may form monolithically within their parental molecular cloud (Banerjee & Kroupa 2015) or they may form hierarchically from mergers between subclusters (Bonnell et al. 2003;Arnold et al. 2022).Existing kinematic studies of young subgroups have found no evidence that they are in the process of merging (Kuhn et al. 2019) and therefore if this process takes place if must happen very early, most likely while the subgroups are still heavily embedded.Directly observing these mergers may therefore be very difficult without infrared astrometry.However, kinematic signatures of the merger process could be observable after the merger (e.g., Parker & Wright 2016;Arnold et al. 2022), including inverse energy equipartition, recently observed in the young cluster NGC 6530 (Wright & Parker 2019), and hinting at the possibility of observing such signatures in other groups.
The evolution of a young group depends on its mass, structure, dynamics and the gravitational potential from the surrounding molecular gas (e.g., Parker et al. 2014).Denser systems will usually evolve faster (Spitzer 1987), mixing, relaxing, and evolving towards energy equipartition.The initial structure and state of the system, the distribution of stellar masses, and the external gravitational potential will all affect the timescale on which a given system evolves.
The vast majority of young groups do not survive to maturity as long-lived open clusters.The process of residual gas expulsion was long thought to be responsible for unbinding young clusters by dispersing the residual gas and its associated gravitational potential (Tutukov 1978;Lada et al. 1984).However, recent studies have questioned the role that residual gas expulsion may play (Dale et al. 2015) or have considered other mechanisms for cluster dispersal such as tidal interactions (Kruijssen 2012).While some older studies of dispersing OB associations could not find clear evidence that they were dispersing from a more compact configuration (e.g., Wright et al. 2016;Dzib et al. 2017;Wright & Mamajek 2018;Ward et al. 2020), recent studies that identify members of the association kinematically have found the majority show strong expansion patterns (Cantat-Gaudin et al. 2019;Kuhn et al. 2019;Armstrong et al. 2020;Quintana & Wright 2021;Quintana et al. 2023).
The physical processes that drive the formation, evolution and dispersal of star clusters and associations are clearly poorly understood.This is due to the difficulty observing young stars that are often highly embedded in their parental molecular cloud, confirming the youth of stars in low-density OB associations, and diagnosing the highly-stochastic physical processes that are predominantly driven by gravitational attraction, a very difficult force to trace or directly observe.Kinematic information can allow the physical processes involved to be constrained, and this is quickly becoming more abundant thanks to Gaia and data from large-scale spectroscopic surveys such as the Gaia-ESO Survey (GES).
GES (Gilmore et al. 2022;Randich et al. 2022) is a large public survey programme that was carried out using the FLAMES (Fibre Large Array Multi Element Spectrograph, Pasquini et al. 2002) fibre-fed spectrograph at the VLT (Very Large Telescope).The aim of the survey was to understand the formation and evolution of all components of our Galaxy, achieved by measuring radial velocities (RVs) and chemical abundances.Over a 6 year period the survey observed approximately 10 5 stars in our galaxy, including in 62 young and open clusters, groups and associations (Randich et al. 2022).
In this paper we are targeting the 18 young (< 20 Myrs) clusters, associations and star-forming regions observed by the Gaia-ESO Survey, listed in Table 1.In Appendix A we provide a summary of the 18 groups studied with information from the literature.We combine GES spectroscopy and Gaia astrometry for thousands of stars towards these 18 young groups, confirming the youth of the stars using spec-troscopy and combining spectroscopic RVs with astrometric proper motions (PMs) to facilitate the largest 3D kinematic study of young stellar groups to date.We chose to focus on young groups as their dynamics provide a valuable probe of the star cluster formation process, relatively unspoilt by dynamical mixing and evolution.
Section 2 introduces the targeted groups and their properties.Section 3 describes the spectroscopic and astrometric data used, while Section 4 outlines our process for confirming the youth of the targeted stars.In Section 5 we measure 3D velocity dispersions for all the groups, exploring a variety of models to determine how best to measure the velocity dispersion of a system, and then calculating the virial state of each system.In Section 6 we search for evidence of expansion or contraction in our sample using different methods to quantify the presence of expansion and the rate at which each group is expanding.In Section 7 we compare the kinematic properties of these groups with their physical properties and in Section 8 we discuss our results and their implications for the formation and dispersal of star clusters.

OBSERVATIONAL DATA
Here we describe the observational data, GES spectroscopy and Gaia early data release 3 (EDR3) astrometry for stars towards our 18 targeted groups.In the work that follows we discuss the 18 groups as separate targets, although in this section groups that were part of the same observing block were processed together before being separated into their constituent groups.

Gaia-ESO Survey spectroscopy
The main source of data for this work is spectroscopy from the sixth and final internal data release (iDR6) of GES.GES targets are selected using homogeneous criteria based on the positions of stars in various colour magnitude diagrams (CMD), unbiased with respect to kinematics, though the exact strategy varies from group to group (Bragaglia et al. 2022).The observations of each group are not complete within any area or magnitude range, being limited by the field of view of the instrument and the difficulty placing fibres close to each other.This can introduce slight biases against targets in dense areas or on the periphery of groups, though thanks to the large number of observing blocks and the good spatial coverage of the blocks, these biases are not significant.
GES data were downloaded from the Edinburgh Wide Field Astronomy Unit 1 .GES data are reduced and analysed using common methods and software to produce a uniform set of spectra and extracted stellar parameters.These methods are described in Jeffries et al. (2014) and Sacco et al. (2014) for the data reduction, and in Lanzafame et al. (2015), Pancino et al. (2017) and Hourihane et al. (2023) for the data analysis, calibration and extraction of stellar parameters of stars in young groups.The final full GES data set is also now available through the ESO archive.
In total, spectroscopy for 11,247 targets in our 11 observational datasets was used.These were then subdivided into 18 groups (see Section A for more details) and trimmed of spatial non-members (e.g., in outlying fields not associated with the targeted group), which lead to 10,897 targets, as listed in Table 1.
Where possible, stellar effective temperatures derived directly from spectroscopy were used (available for 8871 or 79% of sources).Where these were not available we estimated temperatures using the τ spectral index introduced by Damiani et al. (2014) and calculated for most GES spectra (available for an extra 767 or 7% of sources).We calibrated a new relationship between T ef f and τ for young stars (< 50 Myr) using sources with both spectroscopic effective temperatures and τ values.Uncertainties were estimated based on the spread in temperature at a given τ value.For sources that lacked a spectroscopic temperature and a τ index we estimated effective temperatures from the available photometry, dereddened using the average group extinction (Table 1), extinction coefficients from Schlegel et al. (1998) and Danielski et al. (2018), and the tables of intrinsic colour from Pecaut & Mamajek (2016), prioritising the colour with the longest baseline.Uncertainties were estimated from the photometric uncertainty and the spread in temperature at a given colour.This was possible for a further 1569 or 14% of sources.For the remaining 38 (0.3%) of sources this was im-possible either due to a lack of photometry or the available photometry being inconsistent with group membership.
RVs are determined using a combination of a crosscorrelation method and a direct modelling approach (Gilmore et al. 2022).For the RV uncertainties we use the improved empirical precision provided by Jackson et al. (2015), which is calculated from a combination of the target signal-to-noise ratio, the projected equatorial velocity (v sin i) and empirically-derived constants (see Jackson et al. 2020, for the values and derivation of these).
Lithium equivalent widths and Hα full widths at half maximum are determined using direct profile integration, taking into account the star's RV, v sin i and S/N (see Franciosini et al. 2022a, for more details).Values of the gravitysensitive spectral index γ are defined and calibrated according to Damiani et al. (2014).All of these values are available in the GES data set at the ESO archive.
Table 1 lists the groups targeted in this study, including the number of targets that were observed with GES towards each group.In total we have spectroscopy for 11,247 objects, of which 10,897 are towards our 18 groups, and of these 10,859 (99.7%) have effective temperatures, 8668 (80%) have RVs, 7366 (68%) have lithium equivalent widths, and 9321 (86%) have γ indices.

Gaia EDR3 astrometry and photometry
We use astrometry and photometry from Gaia (Gaia Collaboration et al. 2016) EDR3 (Gaia Collaboration et al. 2021), which contains parallaxes, ϖ, and PMs, µα and µ δ , calculated from the first 34 months of Gaia observations.These data achieve an extremely high level of astrometric precision for a sample of unprecedented size, though it is still calculated assuming single-star behaviour (Lindegren et al. 2021).EDR3 is effectively complete in the magnitude range G = 3 − 18.5 mag, with a detection threshold of G = 20.7 mag.Extremely bright (G < 3 mag) or high PM (> 0.6 arcsec/yr) stars, or objects in very crowded areas of the sky (predominantly in globular clusters), can suffer from incompleteness (Fabricius et al. 2021), though this should not significantly affect this work.The systematic uncertainties in parallax and PM that were present in DR2 have been greatly reduced, with less variation over the sky (Gaia Collaboration et al. 2021;Lindegren et al. 2021).Comparison with quasars and known binaries indicates that the parallax zero-point has been reduced compared to DR2 and is now ∼ -17 µas, a correction applied to all parallaxes when calculating distances.The effect of this correction is largest for the most distant group, Trumpler 16, where it results in a reduction in the inferred distance of ∼4%.The affect is smaller for all other groups and <2% for groups within 1 kpc.
GES sources were cross-matched with Gaia EDR3 using a radius of 1 ′′ , which resulted in a total of 11,179 matches (99.4%).Gaia astrometry and photometry was then filtered based on the criteria outlined in the Gaia data release papers and technical notes.For the astrometry we required the 're-normalised unit weighted error' (RUWE) to be less than or equal to 1.4 (Lindegren et al. 2021).This removes sources with spurious astrometry and helps filter contamination from double stars, astrometric effects from binary stars and other contamination problems2 .Applying this cut removed astrometry for 1407 sources (12.6%),leaving 9772 sources with Gaia EDR3 astrometry (sources without valid astrometry were not completely discarded as they might still have useful RVs from GES).Gaia photometry is believed to be well-behaved for sources with G < 18, as is the case for the vast majority of our sources, and so no filtering of the photometry was performed.
We do not use stellar effective temperatures from Gaia EDR3 since they are derived under the assumption of no reddening, which is invalid for our targets.We also don't use Gaia RVs for our sources as they are limited to the brightest sources and generally inferior to the GES data.
In the astrometric analysis that follows we always use the full (5 × 5) Gaia covariance matrix whenever propagating uncertainties, and consider all errors to follow a normal distribution (Arenou et al. 2018).

MEMBERSHIP SELECTION
In this section we outline how we have identified the young stars in each group studied.GES spectroscopy is not complete, either spatially or in any magnitude range, and therefore our approach to membership selection is not to attempt to be complete, but to derive a reliable list of members that is as free from contamination as possible.To avoid any kinematic biases we have not applied any sort of membership criteria based on RVs or PMs.
The main source of contamination for candidate young stars selected from a CMD (as GES targets were) are foreground main sequence stars and background giants.Distinguishing young stars from foreground main sequence stars requires an indicator of youth such as the presence of photospheric lithium or strong Hα emission.Identifying and removing background giants can be achieved using either a spectroscopic surface gravity indicator to discern giants from pre-main sequence stars or by using parallaxes.
GES spectroscopy provides three separate measures that can be used for this purpose: the equivalent width of the lithium 6708 Å line, the gravity index γ (Damiani et al. 2014), and the full width at zero intensity (FWZI) of the Hα line.These three parameters have been used by multiple GES studies to define samples of young stars free of kinematic bias (e.g., Jeffries et al. 2014;Rigliaco et al. 2016;Sacco et al. 2017;Bravi et al. 2018;Wright et al. 2019), and we follow the same approach here, complementing this with the use of Gaia parallaxes.

Separating foreground main-sequence stars
Main sequence stars of certain types can be separated from young, late-type stars based on the presence of lithium in their atmospheres.Late-type stars deplete their primordial lithium after approximately 30 (for M-type stars) to 300 Myr (for K-type stars) due to burning and subsequent mixing throughout the convection zone (Soderblom 2010).The presence of lithium in the atmosphere of K-M stars is therefore an effective indicator of youth.We follow the method of Wright et al. (2019) and require our targets to have EW(Li) values greater than observed in stars of the same temperature in the 30-50 Myr cluster IC 2602 (Randich et al. 1997).Out of our 10,897 targets, 7366 have measures of EW(Li) (excluding upper limits) and of these 2384 pass our membership criteria (32.4%).
We note that the EW(Li) can be underestimated in stars with high accretion rates if the continuum emission in excess produced by the accretion shock reduces the measured signal (e.g., Palla et al. 2005).From the 8513 stars that fail the EW(Li) membership test (or lack equivalent widths to perform such a test), 2659 have FWZI(Hα) measurements, and 1096 of these (41.2%) have FWZI(Hα) greater than 4 Å (Bonito et al. 2013;Prisinzano et al. 2019) and are therefore re-classified as members.This leaves 3480 stars as likely young stars based on their spectroscopic properties.
Figures B1 and B2 show both EW(Li) and FWZI(Hα) as a function of T ef f for all stars in the 18 groups studied.

Separating background giants
To separate young stars from background giants (which sometimes show lithium in their atmospheres) we use the surface gravity-sensitive index γ (Damiani et al. 2014), which is capable of distinguishing between main sequence stars, pre-main sequence stars and giants (cool giants with T ef f < 5600 K have γ > 1, Damiani et al. 2014).Of our 3480 current candidate members, 3326 have γ measurements, and of these 413 (12.4%) have γ values that suggest they could be background giants and are therefore removed from our sample.This leaves 2913 stars that we retain as likely young stars.We also retained as likely young stars the 154 stars that lack γ measurements, since the probability of these being giants is relatively low (∼12%) and will be removed by the next step.Figures B1 and B2 show γ as a function of T ef f for the stars in all 18 groups studied.
As a further check to remove non-members we implement a parallax cut using data from Gaia EDR3, which is available for 2532 of the 3067 likely young stars.For each group we fit the parallax distribution with a single Gaussian and derive central values and dispersions.The dispersions of the parallax distributions are corrected for the contribution of the nonuniform parallax errors using the method of Ivezić et al. (2014) to calculate parallax dispersions for each group.
We then exclude all sources with parallaxes more than two standard deviations from the central value, where the standard deviation is the combination of the intrinsic parallax dispersion for the group and the measurement uncertainty for each star.A cut at 2σ was chosen to balance the need to include the majority of members whilst rejecting a reasonable number of non-members.The Col. 197 group was found to be projected against another group of young stars in the background with parallaxes of ϖ = 0.5-0.9mas (see Figure B2) that we isolated by ensuring that the parallax distribution was fitted to the main peak at ϖ ∼ 1.1 mas that represents Col. 197.This cut excludes 384 (15.2%) of the 2532 likely young stars with parallaxes, leaving 2148 high-confidence young stars.We retained the 535 stars without parallaxes as the probability of these being non-members is low (∼14%) and many of the background giants were already removed by the previous step.Note that we do not use parallax to add in any new stars, only to remove background contaminants.Figures B1 and B2 show both γ as a function of T ef f and the parallax distribution for all candidate members.

Final sample of members
Our final sample contains 2683 spectroscopically-confirmed young stars distributed between the 18 groups, of which 1914 have 3D kinematics (71%), 503 only have RVs (18.7%), 234 only have PMs (8.7%), and 32 lack both PMs and RVs (1.2%), the latter of which are discarded to leave a sample of 2651 young stars with either RVs, PMs or both.We do not require the stars in our sample to have both RVs and PMs so as to maximise the numbers of stars with reliable kinematic data in at least one dimension that we can use to constrain the kinematic properties of the groups and clusters studied.Following the Bayesian approach used in this paper there should be no downsides to this, and we have tested this with additional simulations and can confirm that no systematic biases arise (see Section 4.2).
The median RV uncertainty for stars in our sample is 0.59 km s −1 , with some as low as 0.25 km s −1 .The median PM uncertainty is 0.062 mas yr −1 , which equates to between 0.04 and 0.80 km s −1 depending on the distance to the target group.The median parallax uncertainty is 0.074 mas.
Our full catalogue of members, including their photometry, astrometry, and spectroscopic parameters is included in an online table made available on VizieR.

3D VELOCITY DISTRIBUTIONS AND DISPERSIONS
Here we present and analyse the 3D kinematics of stars in the young groups studied.We start by considering the 3D velocity distributions for each group, calculating velocity dispersions in all three dimensions for each group using a forward-modelling Bayesian approach, and from this derive virial masses for each group that allows us to assess their gravitational boundedness.

Velocity distributions
The velocity distributions of stars in the 18 groups studied show a variety of morphologies.A small number of the velocity distributions are broadly Gaussian, albeit with extended wings that might constitute runaway stars, binaries or non-members of the group within our fields of view.However many of the velocity distributions deviate from normality, sometimes subtly and other times with very clear kinematic substructure that is evident even in 1D velocity distributions.An example is shown in Figure 1 for young stars in Barnard 30.The velocity distributions in all three dimensions show evidence for substructure, with clear bimodal distributions in µα and µ δ .We estimate that of the 18 groups, 8 show evidence for kinematic substructure in at least one dimension.We conducted Shapiro-Wilk tests of normality on the velocity distributions of stars in all 3 dimensions of the 18 groups.All of the groups except five show highly significant (>3σ) evidence for non-Gaussianity in all three dimensions (significance of rejecting the null hypothesis of a normal distribution).The groups whose velocity distributions are consistent with Gaussianity in at least one dimension are primarily the poorly-sampled groups (Cha I South, Rho Oph, Barnard 30 and Barnard 35) and therefore this may just be that there are insufficient stars to confidently rule out a normal distribution.The Gamma Vel cluster is the only well-sampled region (∼100 stars) whose velocity distribution is not significantly non-Gaussian in more than one dimension (µα and µ δ ), most likely because the cluster is relatively old (∼19.5 Myr) and therefore more dynamically processed.

Velocity dispersion fitting
To calculate the velocity dispersion for these groups we use Bayesian inference and model the distribution of observed velocities using a simple parameterised model that we compare with the observations in a probabilistic way (see e.g., Wright et al. 2019).The aim is to determine which of the various sets of parameters, θ, best explain the observations, d.In Bayes's theorem this is known as the posterior distribution, P (θ|d) = P (d|θ) P (θ)/P (d), where P (d|θ) is the likelihood model, P (θ) are the priors (which includes our a priori knowledge about the model parameters) and P (d) is a normalising constant.Bayesian inference allows the model predictions to be projected into observational space, where the measurement uncertainties are defined.This is particularly important when the observational uncertainties are both heteroskedastic and correlated, as is the case here.
We consider three different models for the velocity dispersion of each group, approximating their distributions as Gaussians3 .The first is a trivariate Gaussian aligned to the observed coordinate system (α, δ, and along the line of sight), hereafter known as the unrotated model.This model is used because it represents one of the most commonlyused fitting methods employed in the literature4 .The second model replaces the plane of the sky components with a ro-tated bivariate Gaussian, alongside the third, line of sight, component (the partially rotated model).The third model utilises a rotated trivariate Gaussian, allowed to fully rotate in all three dimensions (the fully rotated model).These models require 6, 7, and 9 model parameters, respectively (3 central velocities, 3 velocity dispersions, and 0, 1, and 3 rotation angles respectively).In addition we utilise an unrotated trivariate Gaussian to represent the non-cluster component (which may include runaway stars or nearby young stars not part of the group under study), requiring a further 6 model parameters for the Gaussians and a seventh parameter, f f ield to represent the fraction of stars that are not members of the group.For each of the three model cases, a population of N = 10 5 stars are modelled with 3D velocities sampled from either the cluster or non-cluster components and projected into the observational space as necessary, assuming all stars are at the same distance (a reasonable assumption given the line-of-sight distance dispersion is small compared to the group distance).For the first two models we use all 2651 stars in our sample for the fits, while for the third model we are limited to the 1867 stars with both by performing the same fit in the Galactic Cartesian coordinate system (U V W ), obtaining consistent results.RVs and PMs.We repeated the first to the first two models only using the 1867 stars with 3D kinematics and found fully consistent results, albeit with larger uncertainties, showing that this approach does not bias our results.Unresolved binary systems will broaden the observed RV distribution due to the contribution that binary orbital motion makes to the measured velocity.To simulate this process we follow Odenkirchen et al. (2002) and Cottaar et al. (2012) by assuming that a fraction of our sample are in binary systems5 .The fraction of binary stars in young stellar groups is poorly constrained and so we set the binary fraction to be 46%, appropriate for solar-type field stars (Raghavan et al. 2010).The primary star masses were sampled from a standard mass function (Maschberger 2013) in the mass range 0.5 -1.0 M⊙ (appropriate for the typical stars observed by GES), while the secondary masses were selected from a power-law companion mass ratio distribution with index γ = −0.5 over the range of mass ratios q = 0.1-1.0(Reggiani & Meyer 2011).The orbital periods were selected from a log-normal distribution with a mean period of 5.03 and a dispersion of 2.28 in log10 days (Raghavan et al. 2010).The eccentricities were selected from a flat distribution from e = 0 to a maximum that scales with the orbital period as proposed by Parker & Goodwin (2009).We then calculate instantaneous velocity offsets for the primary and secondary stars (relative to the centre of mass of the system) at a random phase within the binary's orbit, and then use the luminosity-weighted average of the two as the photo-centre velocity6 , which is then added to the modelled RV.Note that we do not correct for the effects of binarity on the measured PMs since this is estimated to be much smaller than the effect on RVs (e.g., Jackson et al. 2020).
Finally we add measurement uncertainties for the RVs and PMs for each star, randomly sampling these from the observed uncertainty distributions and include the correlated PM uncertainties from Gaia EDR3.This produces our fully forward-modelled velocity distribution models.
The three models used have 13, 14 and 16 free parameters for the unrotated, partially rotated and fully rotated models, respectively.Wide, uniform, and linear priors were used for each of these parameters covering 0 to 100 km s −1 for the velocity dispersions, −100 to +100 km s −1 for the central velocities, 0 to 0.2 for the field star component and 0 to 180 • for the rotation angles.
To sample the posterior distribution function we use the Markov-Chain Monte Carlo (MCMC) ensemble sampler (Goodman & Weare 2010) emcee (Foreman-Mackey et al. 2013) and compare the modelled probability density function to the observations using a 3D unbinned maximum likelihood test, which is made efficient by the smooth velocity distributions modelled.For the MCMC sampler we used 1000 walkers and 2000 iterations, discarding the first half as a burn-in.The parameters were found to have similar autocorrelation lengths, typically of the order of ∼200 iterations, resulting in ∼10 independent samples per walker.The posterior distribution functions were found to follow a normal distribution, and thus the median value was used as the best fit, with the 16 th and 84 th percentiles used for the 1σ uncertainties.

Velocity dispersion fitting results
Velocity dispersion fitting was performed for all 18 groups using all three models described above.Table 2 lists the fitted 1D velocity dispersions determined using the unrotated and partially-rotated velocity dispersion models.The 1D velocity dispersions fitted using these methods vary from 0.2 -6.1 km s −1 , with the majority around ∼1 km s −1 .Example forward-modelled fits are shown in Figure 2 for the unrotated model for the Gamma Velorum cluster and in Figure 3 for the partially rotated model for ASSC 50.
For the fully-rotated model, not all fits provided an improvement on the partially-rotated model fits.Bayesian Information Criterion (BIC, Schwarz 1978) was used to assess the improvement in the model fit achieved for the more complex models (those with partial or full rotation of the velocity ellipsoid) by applying a penalty to the likelihood of the more complex models compared to the unrotated model.If the BIC indicated that the more complex model did not provide a sufficiently better fit to the data than the simple model, then the former models were not considered.This typically happened when the model fitting process could not identify a rotation angle at which the velocity distributions were better fit than with an unrotated model.
Table 3 lists the five groups for which the fully-rotated model provided significantly improved fits to the data compared to the partially-rotated model.The groups for which this was the case are all nearby (within 500 pc) and therefore low-mass.This leads to a bias wherein the 1D velocity dispersions are lower than for the full sample, typically less than 1 km s −1 .Figure 4 shows the fully-rotated model fit for the cluster Cha I North.
The 3D velocity dispersions calculated using the three different models generally agree very well with each other.61% of the unrotated model fits agree with the partiallyrotated model fits within 1σ and 94% within 2σ, approximately as expected for a normal distribution.The agreement is poorer between the fully-rotated model fits and the partially-rotated model fits, with only 40% within 1σ and 60% within 2σ.The more complex models do give consistently smaller 3D velocity dispersions than the less complex models.Notably, the fully-rotated model fitted 3D velocity dispersions are, on average, 10% smaller than the unrotated 3D velocity dispersions.This implies that velocity dispersions previously calculated without considering the rotation of the velocity ellipsoid may have over-estimated both the 1D and 3D velocity dispersions by ∼10%.
Regardless of how the velocity dispersions are measured they are significantly anisotropic.The fraction of groups whose velocity dispersions are inconsistent with isotropy to 1σ is 94% for the unrotated models, 89% for the partiallyrotated models, and 100% for the fully-rotated models.Limiting this comparison to just the PM axis velocity dispersions (to remove any influence of binarity or distance uncertainty) then the fractions are slightly lower at 50% for the unrotated models and 67% for the partially-rotated models.Despite this it is clear that the vast majority of these groups do not have isotropic velocity dispersions.
Our results are generally in good agreement with previous studies when comparing like-for-like, albeit with a slight tendency to recover lower velocity dispersions.Comparing our results to the 1D RV dispersions calculated by Jeffries et al. (2014) and Rigliaco et al. (2016) for the Gamma Vel cluster, Vela OB2 and Rho Ophiuchus, we find excellent agreement (i.e., within 1 σ), as would be expected given that the RV data is the same (though it has been updated as the GES pipeline has improved).Our velocity dispersions for Cha I North and South are also consistent with the velocity dispersions of 0.95 ± 0.18 and 0.87 ± 0.24 km s −1 calculated by Sacco et al. (2017), respectively.Our results are in reasonable agreement with the 2D PM velocity dispersions calculated by Kuhn et al. (2019) for the S Mon cluster, NGC 2244, NGC 6530, Trumpler 14 and Trumpler 16, most agreeing within 1-2 σ, though our velocity dispersions are, on average, slightly smaller.This could be due to differences in the group samples or membership (Kuhn et al. 2019, did not have access to spectroscopic data to confirm the youth of their sample).Finally, our 3D velocity distribution is smaller than that calculated by Wright et al. (2016) for NGC 6530, despite using similar data.This can be attributed to different membership criteria since Wright et al. (2019) used a combination of spectroscopic and X-ray youth indicators, while we have only used spectroscopy.

Virial state
To measure the virial state of each group we use the virial equation, which in its three-dimensional form is given by where σ3D is the 3D velocity dispersion, G is the gravitational constant, M is the mass and rvir is the radius.We substitute the parameter η = 6rvir/r ef f , where r ef f is the effective (or half-light) radius in projection, to get The parameter η can be derived from the power-law index, γ, of an Elson et al. (1987, hereafter EFF) surface brightness profile (see e.g., Portegies Zwart et al. 2010), from which the effective radius, r ef f , can also be measured.Some studies assume a standard value of η = 10, a reasonable approximation for young clusters, but one that hides a significant level of uncertainty in the true value for a given group or cluster.We used the ancillary datasets mentioned in Section 2 to fit EFF profiles for all our groups, deriving r ef f and γ, with uncertainties calculated using a bootstrapping process.The power-law indexes, γ, and therefore also the η values, are not well constrained for these groups, with values ranging Table 2. Velocity dispersion fitting results for the unrotated and partially-rotated models.All PM velocity dispersions have been recalculated in physical units for ease of comparison with the RV dispersions, the uncertainties for which take into account the full distance uncertainties.µ 1 and µ 2 are the velocity dispersions along the semi-major and semi-minor axes of the partially-rotated velocity ellipsoid, the position angle of which is θ.

Group Unrotated velocity dispersions Partially rotated velocity dispersions
Rho Oph 1.087 Table 3. Velocity dispersion fit results for the fully-rotated model where the fit provides an improvement on the partially-rotated model fit (Table 2), assessed using BIC.v 1 , v 2 and v 2 are the velocity dispersions along the three axes, with rotation angles of θ, ψ and ϕ.
There are no direct measurements of η or γ in the literature for these clusters (or almost any cluster for that matter), so this is difficult to check and represents one of the main uncertainties in virial calculations.
The fitted values of r ef f are also listed in Table 4.These are generally in good agreement with measurements from the literature, though different studies often define or measure the radius in different ways.For example, for Cha I North and South, Luhman (2007) estimate radii of ∼0.1 and ∼0.2 deg.(∼0.3 and ∼0.6 pc) respectively, while we measure radii of 0.42 and 0.40 pc, in broad agreement.For the Spokes Cluster, Maíz Apellániz (2019) estimate a full diameter of ∼1.3 pc, which is consistent with our effective radius of 0.39 +0.03 −0.03 pc, while for Trumpler 14 Figer (2008) measure an effective radius of ∼0.5 pc, consistent with our effective radius of 0.53 +0.05 −0.07 pc.The effective radii are combined with the stellar and gas masses listed in Section A and Equation 2 to give virial velocity dispersions.Table 4 compares these values with our measured velocity dispersions from Tables 2 and 3.For ease of comparison we also provide in Table 4 values of the virial ratio, α = σ3D/σ3D,vir, for which α < 1 indicates a subvirial system, α = 1 indicates a system in virial equilibrium, and α > 1 indicates a super-virial system.Where relevant we have calculated virial velocity dispersions and ratios using just the stellar mass and the sum of stellar and gas masses.
When considering only the contribution of the stellar mass to the gravitational potential, nearly all of the groups are super-virial (α > 1), with only λ Ori and the S Mon cluster in NGC 2264 being in virial equilibrium or sub-virial (α < 1).A small number of systems are close to being in virial equilibrium with α < 2, such as Cha I South, Gamma Vel, 25 Ori, and Trumpler 14.For these clusters if their stellar mass has been under-estimated or their radii or η values over-estimated, then this may bring them into  Table 4. Velocity dispersions, stellar and gas masses, effective radii and virial states for all groups.The 3D velocity dispersion, σ 3D , is the partially-rotated model from Table 2 unless an improved fit was achieved using the fully-rotated model from Table 3.The stellar and molecular gas masses were gathered from the literature, as described in Section A. The effective radii were fitted as described in the text.The virial velocity dispersions were calculated using Equation 2 using either the stellar mass or the sum of the stellar and gas masses.The virial ratio is given by α = σ 3D /σ 3D,vir , again using either the stellar mass or the sum of the stellar and gas masses.virial equilibrium7 .We note that 25 Ori and the Gamma Vel cluster are highly clustered and sufficiently old (19 and 19.5 Myr, respectively) that we would expect them to be gravitationally bound (if they were not then they would have expanded and dispersed), which suggests that certainly these clusters have their masses under-estimated or r ef f or η overestimated.
A different picture emerges if we consider the gravitational potential resulting from both the stellar and gas parts of the local system.Molecular gas masses were gathered from the literature for any group still associated with or embedded within a molecular cloud and for which such data were available (as described in Section A).This was the case for six groups, Rho Oph, Cha I North and South, the two clusters in NGC 2264 and NGC 6530.We note that certain other groups are associated with molecular gas, but estimates of the total gas mass were not available from the literature (Barnard 30 and 35, ASSC 50, NGC 2244, Trumpler 14 and 16), while the remaining groups are not associated with any molecular gas (Gamma Vel, Vela OB2, 25 Ori, λ Ori, Collinder 197 and NGC 2237).When the molecular gas mass is taken into account, all of the groups with such information available are found to be sub-virial, with typical virial ratios of ∼0.5.Some of these groups may not be fully embedded within their molecular cloud, or the molecular gas may be more spatially extended than the stellar part of the system, both of which would mean that these gas masses would be over-estimates in this context.Nonetheless, it is notable that when the molecular gas masses in these regions are taken into account, all such systems are found to be sub-virial and some would still be in virial equilibrium if the gravitational potentials were half that estimated here.

GROUP EXPANSION
Many young stellar systems, particularly OB associations, have recently been found to be expanding (Cantat-Gaudin et al. 2019;Armstrong et al. 2020;Wright 2020).Younger star forming regions have presented mixed results, with some young groups expanding, while others are not (e.g., Kuhn et al. 2019).Inspection of velocity vector maps of these groups show many tend to exhibit a preference for coherent outward motion, particularly NGC 6530, Vela OB2, and ASCC 50.In this section we quantify the level of expansion in the groups in our sample using a variety of measures, correcting all transverse velocities for radial streaming motions (virtual expansion) using their central radial velocities and equation A3 in Brown et al. (1997), before performing this analysis.

Fitting 1D expansion gradients
The traditional method of measuring the expansion of a group of stars is to search for velocity gradients, i.e., correlations between position and velocity in the same axis that would suggest a ballistic expansion.We follow the method 99.8 99.9 100.0 100. of Wright & Mamajek (2018) and fit linear relationships of the form v = Ax+B between the velocity, v, and spatial position, x, in each dimension (α, δ, and along the line of sight, ϖ, if the group is resolved).The gradient, A, and zero-point, B, were fitted by maximising the likelihood function, using the MCMC ensemble sampler emcee to sample the posterior distribution.A third parameter (f ) was introduced to represent the scatter in the relationship (see Hogg et al. 2010) and which was marginalised over to calculate the fit and uncertainties on the fitted velocity gradients.
The results of the velocity gradient fits are listed in Table 5.The fitted expansion gradients vary from −0.99 +0.33   −0.29 (indicating contraction) up to 0.90 +0.45  −0.33 km s −1 pc −1 (indicating expansion), for Trumpler 14 and NGC 2237 respectively.These values are broadly consistent with, but typically larger than, estimates of expansion made for other groups such as OB associations, which typically extend up to 0.1-0.2km s −1 pc −1 (Wright & Mamajek 2018;Quintana & Wright 2021;Armstrong et al. 2022).We explored the effects of removing 2 or 3σ outliers in position or velocity from the samples before performing the fits, but found that it did not have a significant effect on the results.Figure 5 shows an example of a fitted 1D expansion gradients for the S Mon cluster of NGC 2264, with fitted gradients of 0.27 +0.16  −0.16 mas yr −1 deg −1 = 0.07 +0.04 −0.04 km s −1 pc −1 in RA and 0.74 +0.11  −0.13 mas yr −1 deg −1 = 0.20 +0.03 −0.04 km s −1 pc −1 in Dec.These ∼1.5 and ∼5.5 σ measurements of expansion are in strong contrast to the lack of expansion measured along the line of sight (due to the S Mon cluster not being fully resolved with the available Gaia data).
Table 5. Expansion gradients and indicators of expansion for the groups studied in this work.1D expansion gradients were calculated in RA (α), Dec (δ) and along the line of sight (ϖ).Rotated expansion gradients were calculated in the plane of the sky and rotated in steps of 5 • until the largest single expansion gradient (positive or negative) was found.The gradient in this direction and the gradient perpendicular to this direction are given, as well as the angle of maximum expansion gradient.Weighted-median expansion velocities are calculated from the radial component of the 2D (plane of the sky) velocities of all stars in each group.All transverse velocities were corrected for the effects of radial streaming motions.Along the line of sight the results are limited.We do not fit or report expansion gradients for the most distant or compact groups studied, where the precision of the parallax measurements do not allow us to probe the line of sight distribution of sources.Of the nearest 6 groups, Gamma Vel is not resolved, and only 3 of the other 5 groups show evidence of expansion, all at only ∼1σ significance (Cha I North, Vela OB2 and 25 Ori).

Group
In the RA and Dec directions we find more significant evidence for expansion, with 11 out of 18 groups showing evidence for expansion in RA and 11 out of 18 groups showing evidence for expansion in Dec. Notable examples of highsignificance indications of expansion are B35 (3σ evidence for expansion in both RA and Dec), NGC 6530 (6σ evidence for expansion in Dec, but no evidence for expansion in RA, see also Wright et al. 2019), the S Mon cluster in NGC 2264 (6σ evidence for expansion in Dec, but only 1σ in RA), and ASSC 50 (9σ evidence for expansion in both RA and Dec).Notably, we also find that Trumpler 14 is contracting at 3σ significance in Dec and Rho Oph is contracting at 2σ significance in Dec.We will return to these cases later.
In the groups with evidence for expansion in one or both of the plane of the sky directions, approximately half (9 out of 15) have expansion rates that differ significantly (>2σ) between axes (i.e., they have anisotropic expansion), while the remaining 6 are consistent with having isotropic expansion or exhibit anisotropy at only the 1σ level.This is in slight contrast with the strong evidence for anisotropic expansion found in OB associations (Wright 2020).

Fitting rotated expansion gradients
Given that more than half of the groups studied exhibit evidence for anisotropic expansion, and that there is no reason to assume that the strongest expansion would occur along one of our arbitrarily-defined observational axes, we also explored the evidence for expansion along arbitrary axes in the plane of the sky.To achieve this we rotated the 2D plane of the sky PM axes in steps of 5 • , reprojecting the PMs, and repeating the expansion gradient fits as described in Section 5.1.
Table 5 reports the angle at which the strongest evidence for expansion was found, the expansion gradient fit along that axis, as well as along the perpendicular axis.We fit expansion gradients that vary from −0.53 +0.24  −0.21 up to 0.93 +0.46  −0.35 km s −1 pc −1 , for Rho Oph and NGC 2237, respec-  tively.Again, these values are typically larger than previous estimates for OB associations, but in some cases our results are less significant due to the difficulty measuring expansion in a compact cluster.Again we explored the effects of removing 2 or 3σ outliers in position or velocity from the samples, but found that it did not have a significant effect on the results.Figure 6 shows an example of fitted rotated expansion gradients for the S Mon cluster of NGC 2264.These fits were performed along the axis with the strongest measured expansion gradient (with a position angle of 120 • ) and the perpendicular axis.The fitted gradients are 1.30 +0.13 −0.14 mas yr −1 deg −1 = 0.35 +0.04 −0.04 km s −1 pc −1 along the axis of strongest expansion and 0.39 +0.11  −0.13 mas yr −1 deg −1 = 0.11 +0.03 −0.04 km s −1 pc −1 along the perpendicular axis.These ∼8.5 and ∼3 σ measurements of expansion are notably both more significant and stronger (by approximately 50%) than when expansion is measure along the equatorial axes.
We find that the evidence for expansion is significantly stronger when the axis of expansion is allowed to vary, as expected.In Section 5.1 we found evidence for expansion in 11 out of 18 groups along each dimension (albeit with many showing evidence of expansion at only the 1 σ level), while when the axes are allowed to rotate we find evidence for expansion in 17 out of 18 groups (with only 2 at the 1 σ level).Notable examples are λ Ori, Barnard 35 and the S Mon cluster (all at 6σ), NGC 6530 (8σ) and ASSC 50 (13σ).The degree of expansion in the primary expansion axis is strongly correlated with the degree of expansion in the perpendicular axis, with Kendall's rank correlation test giving a correlation coefficient of 0.490 (p-value = 0.0039).
The only group not found to be expanding is Rho Oph, which we have found to be contracting along one or more axes, regardless of the orientation of the axes.
In the groups with evidence for expansion, the level of anisotropy is broadly the same as when expansion was explored along the equatorial axes.10 out of 17 groups show evidence for significantly (>2σ) anisotropic expansion, while the remaining 7 are consistent with either isotropic expansion or mildly significant expansion (at the 1σ level).

Outward motion
Another method of quantifying the presence and significance of expansion is to separate the velocities of stars into their radial and azimuthal components (relative to the centre of each group) in the plane of the sky (e.g., Wright et al. 2016;Kuhn et al. 2019).To do this we must estimate the centre of each group, which we do using the results from the EFF profile fitting (Section 4.4).
We follow Kuhn et al. (2019) and calculate the weighted-median outward velocity, vout, for each group, calculating uncertainties using a Monte Carlo process that takes into account all observational uncertainties as well as the inherent uncertainties on the calculation of a median.The results can be found in Table 5.We measure median outward velocities of up to 1.37 +0.09 −0.28 km s −1 (for Trumpler 16) with most values around 0.1-0.5 km s −1 .These values are consistent with those found by Kuhn et al. (2019), albeit with smaller uncertainties, and for the five clusters in both samples (S Mon cluster, NGC 6530, NGC 2244, Trumpler 14 and Trumpler 16), four of them agree within 1σ.
The median outward velocity method gives broadly similar results to our previous methods, with 15 out of 18 groups showing evidence for expansion by this metric.Notable examples that are consistent with the method of fitting linear expansion gradients include ASSC 50 (20σ) and NGC 6530 (36σ), while notably different results are obtained, for example, for 25 Ori, which is found to be expanding with 16σ significance by the median outward velocity method, but only 4σ using linear expansion gradients.
As with the previous methods, Rho Oph is again found to be contracting using the median outward velocity method, at a significance of 8σ, providing further evidence of interesting kinematic behaviour in this young group.

Summary of expansion results
Our results show that the vast majority of groups are expanding in at least one dimension, with 17 out of 18 groups (94%) showing evidence for expansion when the axis of expansion is allowed to rotate (Section 5.2) and 15 out of 18 groups (83%) having positive values of the median outward velocity (Section 5.3).Only one group shows clear and consistent evidence for contraction and that is Rho Oph, which is contracting according to all our methods (with a significance of 1-8 σ depending on the method).
More than half of all groups show evidence for expansion in at least two dimensions (11 out of 18 or 61%), particularly when allowing the axes of expansion to rotate.The majority of groups that are expanding are doing so asymmetrically, with only 2 out of 17 of the expanding groups being consistent with symmetric expansion within 1σ, those being Cha I South (for which the uncertainties on the expansion gradients are very large) and Gamma 2 Vel (which has very low levels of expansion).
To estimate the uncertainties on the fraction of systems that are expanding (since measurement errors play a large role in the expansion gradients for some systems), we perform a Monte Carlo experiment to determine the underlying fraction of systems that are expanding.We find that the effect of measurement uncertainty is to reduce the fraction of systems observed to be expanding, particularly for the median outward velocity method.We find that, using the rotated expansion gradient fitting method, that 95 +4 −6 % of systems are expanding, while using the median outward velocity method that 99 +1 −1 % of systems are expanding.

COMPARISON OF GROUP AND KINEMATIC PROPERTIES
In this section we compare the physical properties of our sample of groups (their mass, radius and age) with their kinematic properties (velocity dispersion, virial state and expansion rates) to search for correlations that might expose the physical processes at work.The values of group age and mass are included in Tables 1 and 4, respectively, compiled from the literature, many of which do not report uncertainties.In our experience, cluster ages may be inaccurate by up to 50% and cluster masses by 20-30%, which we have included in the fits performed here8 .We also calculate dynamical and relaxation timescales according to the equations in Portegies Zwart et al. (2010).A number of obvious correlations are observed that we do not discuss in detail here, such as strong correlations between the 3D velocity dispersion and the 3D virial velocity dispersion, as well as between the virial ratio and the relaxation timescale.
Figure 7 shows the relationship between the 3D velocity dispersion and the stellar mass of the system.We find a strong correlation between these two quantities using Kendall's rank correlation test with a p-value of 0.0035.A similar correlation between mass and velocity dispersion was observed by Kuhn et al. (2019), who also found a strong correlation between velocity dispersion and group radius, but found that since mass and radius were related (a correlation that we do not find), argued that this correlation was driven by this interdependency.A correlation between these two quantities would be expected based on the assumption of virial equilibrium, though very few of these systems were found to be in virial equilibrium.We fit a relationship of the form σ3D = A + B × M C cluster between these two quantities using Bayesian inference and an MCMC ensemble sampler to derive a fit of 0 1000 2000 3000 4000 5000 6000 7000  (3) which is shown in Figure 7.The posterior distribution on the power-law index, C, derived from this fit gives a probability of 0.09 that the power-law index is ⩽0.5, the dependence that would be expected according to the virial relationship.
We find a strong inverse correlation between velocity dispersion and age with p = 0.0293 (Figure 8).This could be due to both dynamical evolution (the most rapidly moving stars will escape a group over time, causing the measured velocity dispersion of stars within the group to reduce) and an evolutionary bias (systems with smaller velocity dispersions should survive for longer times as visible over-densities of young stars that would be selected and observed).
We observe a strong positive correlation between the virial ratio of a system and its radius (p-value = 0.003), as shown in Figure 9.The outlier in Figure 9 with a large radius is the OB association Vela OB2 and if it is removed the correlation is stronger and suggests a linear correlation between the two quantities.There are two possible reasons for this correlation.The first is an evolutionary effect whereby gravitationally unbound systems (those with αvir > 1) will expand to larger radii, with the more unbound systems expanding faster.This should lead to a correlation between αvir and radius, though larger systems are generally less dense and harder to detect, which might introduce a counter bias.The second reason is that αvir is linearly dependent on the radius, and therefore if all other parameters that αvir depends on (σ3D and M specifically) either remained constant or had dependencies that cancelled each other out, one would expect to observe a linear correlation between αvir and radius.Since this does not appear to be the case (see e.g., Figure 10), the former evolutionary reason may be the dominant cause of this correlation.
We observe a weak inverse correlation between the virial ratio of a system and its stellar mass (Figure 10) with a pvalue of 0.081.However, this relationship is likely affected by an observational bias whereby more massive, unbound groups (such as massive OB associations, which would occupy the top-right of this diagram) would be larger than less massive unbound things, and therefore harder to observe by GES due to the relatively small field-of-view of the FLAMES instrument.It will be necessary to expand our sample of clusters and associations to determine the validity of this correlation.
We observe a strong positive correlation between the virial ratio and the primary rotated expansion gradient, with a p-value of 0.068 (no correlation is found with any of the other measures of expansion), which improves to 0.027 when the only non-expanding system, Rho Oph, is removed (see Figure 11).This suggests that the more super-virial a system is (the higher its virial ratio), the higher its expansion rate will be (at least when determined by fitting rotated expansion gradients), which is consistent with a picture wherein the expansion of a group of stars is dictated by how far out of virial equilibrium they are.Notably we do not observe any correlation between the velocity dispersion and any measure of the expansion of a system (either the expansion gradients or the median outward velocity).We do however observe strong correlations between the primary and secondary expansion gradients and vout, suggesting they are all measuring similar properties of a system.We also observe a very strong positive correlation between the virial ratio of a system and its dynamical timescale (p-value = 7.0 × 10 −5 ).This is most likely a product of the previous correlation, in which the more super-virial a system is, the more it expands, the larger it becomes and the more its density is reduced, and therefore the longer its dynamical timescale becomes.
Finally we also observe a strong positive correlation (pvalue = 0.027) between radius and age, with older systems being larger.Correlations such as this have been known about, in many guises, for over half a century (e.g., Blaauw 1964;Pfalzner 2009;Kuhn et al. 2015;Getman et al. 2018) and are generally interpreted as due to the relaxation or expansion of systems as they age.Related to this is a strong positive correlation between age and the dynamical timescale (p-value = 0.014), which is likely to be a product of this since the dynamical timescale is a strong function of the group density (and therefore radius).We also observe a strong inverse correlation (p-value = 0.0013) between gas mass associated with the system and radius, such that the larger the system the less gas it is associated with.As larger groups are generally older, they are less likely to be associ- Kendall's = 0.397 (P-value = 0.027) Figure 11.Relationship between the virial ratio, α vir = σ 3D /σ vir , of all systems in our sample and their expansion gradient along their primary axis of expansion (excluding the only non-expanding system, Rho Oph), with uncertainties for both.Kendall's rank correlation test reveals a correlation of τ = 0.397 and a p-value of 0.027 indicating a strong correlation.
ated with molecular gas (as the parental molecular cloud is either consumed or dispersed).

DISCUSSION
We have performed the first large-scale 3D kinematic study of multiple stellar groups to determine and compare their structural, kinematic and evolutionary properties.We have measured 3D velocity dispersions, anisotropy levels, virial states, expansion or contraction rates and directions, as well as half-mass radii and expansion timescales.We have compared these properties with each other, as well as with literature ages and group masses, and dynamical and relaxation timescales, to identify possible correlations.In this section we discuss the meaning and implications of these results on our understanding of the dynamics, formation, evolution and dispersal of young stellar groups.

Dynamical and virial state of young groups
The young stellar groups studied in this work typically have 3D velocity dispersions of 1-2 km s −1 , ranging from 0.61 +0.03 −0.02 km s −1 for Gamma Vel to 7.36 +0.19 −0.15 km s −1 for Trumpler 14.These results are in good agreement with previous estimates, but are typically slightly lower.This difference is likely due to our finding that simple 1D or 2D (unrotated) velocity dispersion models are likely to over-estimate the velocity dispersion by 10% or more compared to more advanced 2D or 3D (rotated) velocity dispersion models.
We find that nearly all groups exhibit anisotropic velocity dispersions.This is generally taken to mean that the group is not sufficiently dynamically mixed to have developed an isotropic velocity dispersion.Anisotropic velocity dispersions have been observed in many regions, particularly in OB associations (Wright et al. 2016;Wright & Mamajek 2018) where it can be exacerbated by kinematic substructure within the association.
We find that the majority of groups are super-virial (when considering the gravitational potential due to the stellar mass), with only two groups in virial equilibrium (λ Ori and the S Mon cluster in NGC 2264).When the mass of both stars and the surrounding molecular cloud are taken into account, all six groups with estimated molecular cloud masses in the literature are sub-virial.However, these molecular clouds are significantly more extended than their associated groups and therefore it may not always be appropriate to consider their full mass when estimating their gravitational potential.
In Section 6 we observed a positive correlation between the velocity dispersion of stellar groups and their stellar mass.Such a relation would be expected if all groups were in virial equilibrium since the virial velocity dispersion scales as M 0.5 (Equation 1).However, we observe and fit an approximately linear relationship between group mass and velocity dispersion with a power-law index of 0.94 +0.42  −0.31 and were able to rule out a scaling of M 0.5 with a confidence of 93%.Given this, the ratio of the velocity dispersion to the virial velocity dispersion would be expected to scale as The mass dependence of cluster radii has been studied by various authors with dependencies varying from approximately R ∝ M 0.25 (Brown & Gnedin 2021;Dobbs et al. 2022) to R ∝ M 0.5 (Adams et al. 2006;Pfalzner 2011).This leads to a mass dependency of α ∝ M 0.625 to M 0.75 , implying that less massive groups are more likely to be born in a virial or sub-virial state than more massive groups.When comparing the virial ratio with the group mass, we were unable to identify a correspondingly strong correlation due possibly to selection biases.

Expanding star clusters and groups
We find that the vast majority of groups show evidence that they are expanding.The exact fraction depends on the method used to measure expansion, varying from 83% using the median outward velocity, 89% using 1D expansion gradients, and 94% using 2D rotated expansion gradient fits.The significance of these results vary, but 12 of the 17 groups with evidence for expansion (71%) do so at the 3σ level.The strongest expansion is found for the groups NGC 2237 and Barnard 35, which are expanding at rates of 0.93 +0.46 −0.35 and 0.79 +0.08 −0.14 km s −1 pc −1 (at significances of 2.7 and 5.7σ), as determined using 2D rotated expansion gradient fits.
We find that a larger fraction of groups are expanding compared to previous studies.Kuhn et al. (2019) found that 75% of their groups had positive median outward velocities (compared to 83% here), but that only 57% of their groups had positive median outward velocities to a confidence of >1σ (compared to 78% here).In studies of larger OB associations, many historical studies have struggled to identify large-scale expansion patterns (e.g., Wright et al. 2016;Wright & Mamajek 2018), but more recent studies that dissected OB associations into subgroups using kinematic data have almost universally been able to identify expansion (e.g., Cantat-Gaudin et al. 2019;Armstrong et al. 2022;Quintana et al. 2023).Our results support the view that the subgroups of OB associations represent the expanded remnants of compact groups similar to the clusters studied in this work and that OB associations are therefore composed of multiple expanding star clusters (e.g., Wright et al. 2023).
The majority of groups that are expanding are doing so anisotropically, i.e., with different gradients along different axes.This is the case whether the expansion gradients are fitted in 1D or 2D.This implies that, if no other forces are acting on the stars during their expansion, that these groups either started their expansion from non-spherical initial conditions or had anisotropic velocity dispersions prior to expansion.This is consistent with previous studies that have found that the majority of young groups are non-spherical (e.g., Kuhn et al. 2014) and our result that the majority of groups have anisotropic velocity dispersions.Gravitational tidal forces from the surrounding molecular cloud (Kruijssen et al. 2012) or the residual molecular gas that has since been dispersed by feedback (Zamora-Avilés et al. 2019) could both lead to asymmetric expansion of the group if the gas is non-spherically distributed.
There are two groups in our sample that are in virial equilibrium but observed to be expanding; the S Mon cluster of NGC 2264 (αvir = 0.94 +0.68  −0.01 ) and λ Ori (αvir = 0.43 +0.92  −0.01 ).If these systems are in virial equilibrium then these results would appear to be counter-intuitive.However, we note that both of these groups have (primary) rotated expansion gradients at the lower end of those measured (respectively 0.35 and 0.20 km s −1 pc −1 ).This is part of a strong correlation between the virial ratio and the (primary) rotated expansion gradient that has a p-value of 0.034.If these systems are in virial equilibrium then their observed expansion may simply be due to the system settling down into an equilibrium configuration following formation (Parker & Wright 2016;Sills et al. 2018).

Expansion timescales and ages
Expansion timescales can be calculated from the expansion gradients by simple inversion (with a correction factor to account for units).We calculate expansion timescales for all our groups using the largest expansion rate (from the rotated expansion gradient fit) as this will give the smallest expansion timescale for the two expansion axes.We find that the expansion timescales vary from ∼1 Myr for NGC 2237 and Barnard 35, to ∼8 Myr for 25 Ori and Trumpler 14, and ∼20 Myr for Gamma Vel (though the expansion rates for Trumpler 14 and Gamma Vel are of marginal significance).
The majority of groups have expansion timescales that are very similar to their evolutionary ages, though there are some outliers that mean the measured correlation is very weak with a p-value of only 0.25.Multiple factors can affect the agreement between expansion timescale and age of the system, including delayed expansion (the group not expanding immediately after star formation), non-compact initial conditions (the group not expanding from an initially compact configuration), or additional forces acting on the stars (which may act to either accelerate or decelerate the expansion, e.g., Zamora-Avilés et al. 2019).These factors can each cause the expansion age to be over-or under-estimated relative to the true age of the system.
There are three systems with expansion timescales that are significantly smaller than their evolutionary ages; Vela OB2 (2.3 vs 14 Myr), λ Ori (4.9 vs 10 Myr) and Col. 197 (1.6 vs 5 Myr).Vela OB2 is an OB association with considerable kinematic substructure (e.g., Cantat-Gaudin et al. 2019;Armstrong et al. 2022) that would not have expanded from initially compact initial conditions, which goes a long way towards explaining this disagreement.λ Ori has an expansion timescale of ∼5 Myr, significantly smaller than its evolutionary age of ∼10 Myr from Bell et al. (2013), but there is some disagreement in the age of this system as Kounkel et al. (2018) estimate an evolutionary age of 4-5 Myr, which would be in better agreement with the expansion age.Finally, the expansion timescale of ∼1.6 Myr for Col. 197 is in sharp contrast to its evolutionary age of ∼5 Myr.Since the group is embedded within an H ii region (Gum 15) this might indicate that the group is younger than previously thought, or alternatively that the expansion of the system did not begin immediately after formation.
There are two systems with expansion timescales that are significantly larger than their evolutionary ages; Trumpler 14 (8 vs 2.5 Myr) and the Spokes cluster within NGC 2264 (4.1 vs 1 Myr).The former has large uncertainties on its evolutionary age, so this is not significant, but the latter is more significant.The most likely explanation for this is that this cluster did not begin its expansion from highly compact initial conditions, but rather did so from conditions more similar to how it is currently observed and has been slowly expanding since formation.

Contracting systems: Rho Ophiuchus
There is one system in our sample that is contracting rather than expanding: Rho Ophiuchus.Some systems exhibit negative expansion gradients along one axis, negative 1D expansion gradients or negative median outward velocities, but for nearly all of these systems these measurements are either insignificant or give different results depending on how the expansion is measured.Rho Oph however appears to be contracting according to all three methods used to measure expansion and contraction, with approximately 2.5σ significance using the 1D or 2D expansion gradient fits, and with ∼4σ significance using the median outward velocity method.This implies that Rho Oph is almost definitely contracting.
The virial ratio of Rho Oph, when considering only stars, is αvir = 2.93 +0.38  −0.13 , suggesting that the system is gravitationally unbound.However, when one takes into account the approximately 1750 M⊙ (Loren 1989) of gas in the molecular cloud the virial ratio drops to 0.70 +0.09 −0.03 suggesting that the group is gravitationally bound, and in particular is sub-virial.Rigliaco et al. (2016) came to a similar conclusion regarding the virial state of Rho Oph comparing their velocity dispersion to the mass of the surrounding molecular cloud and estimated an 80% probability that the group is gravitationally bound.The young stars observed by GES are part of the optically visible 3 Myr (Grasser et al. 2021) population that surrounds (on the plane of the sky) the L1688 molecular cloud and so presumably would feel the gravitational pull of its mass.The molecular cloud could therefore be responsible for the contraction of the surrounding group of young stars.This may lead to the accretion of these young stars onto the group of forming stars within L1688.If so, this would provide some of the first direct kinematic evidence of a group accreting other, already formed, young stars.Age [Myr] Figure 12.Relationship between the dynamical timescale, t dyn = (GM/r 3 vir ) −1/2 , of all systems in our sample and their literature age, with uncertainties for the former.Vela OB2 (which has an age of 14 Myr and a dynamical timescale of ∼120 Myr) and Barnard 35 (age 2.6 Myr and dynamical timescale ∼45 Myr) are not shown.A 1:1 relationship is shown as a grey dashed line.
Simulations of star cluster formation show that many systems undergo an initial collapse during the first crossing time of the system (Proszkow et al. 2009;Parker & Wright 2016), that could be observed as contraction.This collapse potentially leads to mergers between sub-groups (Bonnell et al. 2003;Vázquez-Semadeni et al. 2017).Observational support for the model of star cluster formation through mergers has been mixed however.Parker & Wright (2016) and Arnold et al. (2022) identified possible kinematic signatures of past subgroup mergers that could be used to determine if mergers have taken place after the fact, but so far only NGC 6530 has been observed to exhibit any such signature (Wright & Parker 2019).In their study of the dynamics of young groups, Kuhn et al. (2019) identified the group M17 as a system potentially undergoing collapse and mergers due to its negative median outward velocity and highly clumpy structure.However, they found no evidence for the converging motion of subgroups towards other subgroups in their sample, suggesting that if mergers take place they do so predominantly at early ages.It is therefore even more notable that Rho Ophiuchus is in the process of contracting given its estimated age of ∼3 Myr.

The formation of star clusters
How young and compact star clusters form is still an open question, with competing theories suggesting they either form monolithically in a dense, compact distribution (e.g., Banerjee & Kroupa 2015) or that they form over a wider area and collapse down to form a compact cluster (e.g., Bonnell et al. 2003;Arnold & Wright 2024).Studies suggest that the question of how bound star clusters form is related to both the star formation efficiency (Kruijssen 2012) and the collapse of the giant molecular cloud (Chevance et al. 2022).
If star clusters form monolithically then they should have existed in more-or-less their current configuration since birth, and if they are sufficiently dense then they will have had time to dynamically evolve on a timescale set by their dynamical timescale.However, the majority of groups we have studied are not dynamically evolved (lacking for example, isotropic velocity distributions), suggesting that they haven't existed in a grouped configuration for very long.
We compare the age of all the groups in this study with their dynamical timescale in Figure 12.This figure shows that the majority of groups (12 out of 18) are older than their dynamical timescale.The groups that are younger than their dynamical timescale (those to the right of the 1:1 line) are predominantly OB associations (Vela OB2 and ASSC 50) or very low-mass groups (Barnard 30 and 35), though it also includes NGC 2237 and Col. 197.If we exclude the known OB associations from this comparison and focus only on the known clusters or cluster candidates then the fraction increases to 12 out of 16, or 75%.This suggests that these systems have not existed in their current configuration since birth, otherwise they would have had time to dynamically evolve and establish some degree of isotropy.This argues against these clusters having formed monolithically and instead favours a picture whereby these clusters formed at a lower density (with a longer dynamical timescale) and have since collapsed down to form more compact clusters.
Further evidence in favour of the star cluster formation model of collapse and mergers is that all of the systems that are associated with a molecular cloud for which masses are available in the literature are in a subvirial state once the mass of the surrounding molecular cloud is taken into consideration.This shows that in these cases there is sufficient mass in the molecular cloud to introduce a state of contraction in these groups, which could then lead to the formation of a compact star cluster.This picture is further bolstered by the observation that Rho Oph is currently in a state of contraction.

The dispersal of young star clusters and the survival of long-lived open clusters
While the clustering of very young stars is an almost ubiquitous phenomenon, most stars do not find themselves in bound star clusters by an age of 10 Myr or older (Lada & Lada 2003).The explanation for this is believed to be a combination of the unbinding of gravitationally-bound (embedded) clusters by residual gas expulsion(e.g., Hills 1980) and the fact that most groups of clustered stars in star forming regions are actually gravitationally unbound.Our study has provided the strongest evidence to date that the majority of young groups are in the process of expanding and do so from an age as young as ∼1-2 Myr.Our finding that >90% of young groups are in the process of expanding is consistent with the observation by Lada & Lada (2003) that 90% of young 'clusters' disperse within 10 Myr.Given this, what determines which clusters will go to become long-lived open clusters?Krumholz & McKee (2020) suggest that the clusters that do survive their early evolution and go on to become gravitationally bound are distinct kinematically, with isothermal, spherically symmetric density distributions, virialised velocity distributions and are neither expanding or contracting.Our study has very few groups that meet these criteria.The fraction of groups that are out of virial equilibrium (when only considering the stellar part of the system) is 94% (only two groups are consistent with being in virial equilibirum).The fraction of groups that are inconsistent with having isotropic velocity distributions is >90% (depending on the method used to measure the velocity dispersion), with only the Gamma Vel cluster and NGC 6530 being consistent with isotropy.And finally, the vast majority of groups show significant evidence for expansion or contraction, with only Gamma Vel and Trumpler 14 being consistent with not being in a state of expansion or contraction.
We conclude that the vast majority of the systems studied will expand and disperse within the next 10-20 Myr.Even the six groups that are in virial equilibrium when considering the surrounding molecular cloud are likely to disperse once they spatially decouple from the gas (including Rho Oph, despite its currently contracting state).λ Ori may survive as a long-lived open cluster since it appears to be in virial equilibrium, though it is in a state of weak expansion (this may represent the system settling down into a stable state following formation).Both Gamma Vel and Trumpler 14 are close to being in virial equilibrium and do not exhibit strong expansion.If their stellar masses have been underestimated or their radii over-estimated then they may be in virial equilibrium.They appear to be the two best candidates for becoming long-lived open clusters from the sample of groups studied here.

CONCLUSIONS
We have presented the first large-scale 3D kinematic study of multiple young star clusters and OB associations.The combination of Gaia astrometry with Gaia-ESO Survey spectroscopy provides 3D positions and kinematics, as well as kinematically-unbiased indicators of youth, for ∼2700 young stars in 18 clusters or associations.
We measure 3D velocity dispersions for all 18 groups that range from 0.62 to 7.5 km s −1 (1D velocity dispersions range from 0.36 to 4.3 km s −1 ).We find that the majority of groups have anisotropic velocity dispersions, suggesting they are not dynamically relaxed.From the 3D velocity dispersions, measured radii and estimates of total mass in the literature we determine the virial state of all groups and find that all but two systems are super-virial when only the stellar mass is considered, but that some systems are in virial equilibrium when the mass of the surrounding molecular cloud is taken into account.We observe an approximately linear correlation between the 3D velocity dispersion and the group mass, implying that the virial state of groups should scale M 0.625 -M 0.75 .However, we do not observe a strong correlation between virial state and group mass.
In agreement with their virial state we find that nearly all of the groups are in the process of expanding, as indicated by both linear expansion models and the median outward velocity of stars in the group.Given their viral state and expansion, these systems are expected to continue to expand and form OB associations (or subgroups within OB associations).In the majority of cases the expansion is anisotropic, implying that either groups are not spherical or have anisotropic velocity dispersions prior to expansion, or that additional forces act on the group during their expan-sion.Given that most groups are not dynamically relaxed and that other observations find many young groups to have substantial levels of ellipticity, the former explanation is argued to play at least a contributing role.
One group, Rho Oph, is found to be contracting using all measures of expansion/contraction.The group is in a super-virial state when only the stellar mass is considered, but is sub-virial when the mass of its molecular cloud is considered.Whether or not the group is currently gravitationally bound, it is currently contracting, which may lead to mergers between subgroups or the accretion of stars onto a central cluster.
We conclude that, since the majority of clusters are not dynamically evolved, despite being older than their dynamical timescales, that these clusters did not form as we observe them now, but originally were larger and had lower densities (and thus had longer dynamical timescales).Combined with other evidence, we conclude that most clusters form through the collapse of an extended distribution of stars, with mergers between subgroups, and did not form monolithically.
We also conclude that the majority of the groups studied here will not survive as long-lived open clusters, being super-virial, having non-isotropic velocity dispersions and showing significant evidence for expansion.We conclude that they will expand and disperse into the field population of our galaxy.The best candidates for survival as long-lived open clusters are Gamma Vel, Trumpler 14 and λ Ori, each of which are either in or close to virial equilibrium.
The data and sample presented here provide a powerful illustration of the scientific potential that will arise from the combination of Gaia and data from the next generation of multi-object spectroscopic surveys such as WEAVE, 4MOST and SDSS-V.Spectroscopy from these surveys for tens to hundreds of thousands of young stars will allow this work to be extended, better sampling individual groups while also targeting a larger number of groups with a greater range of ages, masses and densities.This will help address a range of outstanding issues that this work has probed including how star clusters form, the dynamical processes at work within them, and the physical processes that drive their disruption and dissolution.
omy Unit, Institute for Astronomy, University of Edinburgh, which is funded by the UK Science and Technology Facilities Council.This work also made use of results from the European Space Agency (ESA) space mission Gaia.Gaia data are being processed by the Gaia DPAC.Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia Mul-tiLateral Agreement (MLA).The Gaia mission website is https://www.cosmos.esa.int/gaia.The Gaia archive website is https://archives.esac.esa.int/gaia.This paper has been typeset from a T E X/ L A T E X file prepared by the author.

APPENDIX A: THE OBSERVED GROUPS
Here we summarise the information, from the literature, on the 18 young clusters, associations and star-forming regions observed by the Gaia-ESO Survey that we have studied in this work.We also briefly summarise how the targeted stars in the GES observations were separated into groups.
• Rho Ophiuchus is our closest target (d = 136 pc), and one of the closest star-forming regions.The entire Ophiuchus cloud complex spans ∼ 6 × 3 degrees on the sky (14 × 7 pc), though the nine FLAMES pointings obtained by GES cover an area of only 1.4 × 1.4 degrees (3.4 × 3.4 pc) centred on the dense molecular cloud L1688, but focussed on the optically-visible non-embedded stars (Rigliaco et al. 2016).The total optically-visible mass of stars towards L1688 is ∼85 M⊙ (Erickson et al. 2011;Rigliaco et al. 2016), which rises to 106 M⊙ to account for unresolved binaries9 with an additional ∼1750 M⊙ in gas in L1688 (Rigliaco et al. 2016).The young stars studied here are part of Population 1 from Grasser et al. (2021) with an average age of ∼3 Myr.The GES data for Rho Oph were previously studied by Rigliaco et al. (2016).
• Chamaeleon I (North and South), our most reddened targets with AV ∼ 3 (Luhman 2007) and at distances of 191 and 187 pc, respectively, are also the youngest targets with an age of ∼1.5 Myr (Galli et al. 2021).Most of the known members are confined to an area of ∼ 0.7 × 1.7 degrees in the form of two groups, Cha I North and South (Galli et al. 2021), each with a radius of 0.35 • (1.1 pc).The total number of primary stars in Cha I has been estimated to be 226 (Luhman 2007), which for a mean stellar mass of 0.4 M⊙ and accounting for unresolved binaries equates to ∼113 M⊙.Cha I North has ∼10% fewer stars than Cha I South (Galli et al. 2021) and so we estimate their total stellar masses as 54 and 59 M⊙, respectively.The mass of the molecular cloud in which Cha I is embedded is estimated to be 1000 M⊙ (Mizuno et al. 2001).We divide the data for Cha I into the North (δ > −77 • ) and South (δ < −77 • ) clusters (Roccatagliata et al. 2018).The GES data for Cha I were previously studied by Sacco et al. (2017).
• Gamma Vel and Vela OB2, aged ∼ 19.5 and 14 Myr respectively (Jeffries et al. 2017;Armstrong et al. 2022), are our oldest targets (respectively Populations A and B from Jeffries et al. 2014).Gamma Vel is a small, non-embedded cluster towards the Vela OB2 association (Jeffries et al. 2014;Armstrong et al. 2020), lying at distances of 334 and 367 pc, respectively.The cluster has a radius of approximately 1.37 pc (Jeffries et al. 2014), while the association spans 50-100 pc (Armstrong et al. 2022).The total mass of stars is estimated to be ∼152 M⊙ for Gamma Vel (Jeffries et al. 2014) and 1285 ± 110 M⊙ for Vela OB2 (Armstrong et al. 2018).We follow Armstrong et al. (2020) and select stars having (µα + 6.53) 2 + (µ δ − 9.8) 2 ⩽ 0.7 2 mas/yr as being members of Gamma Vel, while all other stars are considered part of Vela OB2.The GES data for Gamma Vel and Vela OB2 were previously studied by Jeffries et al. (2014) and Franciosini et al. (2018).
• 25 Ori is a young, non-embedded group in the Orion OB1 association, which Franciosini et al. (2022b) recently determined an age of 19 +1.5  −7.0 Myr for.It lies at a distance of ∼339 pc and has a radius of approximately 0.62 • or 3.7 pc (Kharchenko et al. 2005).The mass of primary stars within a 1 • radius is estimated to be 324±25 M⊙ (Suárez et al. 2019), and when accounting for unresolved binaries this equates to 400 ± 30 M⊙.We exclude stars outside of the central group as non-members (those with δ < 1.2 • ).The GES data for 25 Ori has not previously been studied dynamically.
• λ Ori is a young (∼10 Myr, Dolan & Mathieu 2002), non-embedded group in the Orion OB1 association at a distance of ∼389 pc.The total mass of primary stars was estimated by Dolan & Mathieu (2002) to be 450-650 M⊙.
When accounting for unresolved binaries and revising for the nearer distance from Gaia we estimate the total stellar mass to be 650±120 M⊙.The group radius is approximately 3 • or ∼20 pc.The GES data for λ Ori has not previously been studied dynamically.
• Barnard 30 & Barnard 35 are two dark clouds in the vicinity of λ Ori (north-west and south-east of λ Ori respectively) at similar distances (384 and 390 pc respectively), and both containing young stellar groups.They were observed as part of the λ Ori observations.Their ages are estimated to be 2.4±1.3 and 2.6±1.3Myr (Kounkel et al. 2018), notably younger than the nearby λ Ori group.Their total stellar or gas masses are not precisely known, but Kounkel et al. (2018) identify 96 and 117 optically-exposed members, respectively, with an estimated completeness of 80%.Combining this with a mean stellar mass of 0.4 M⊙ and accounting for binaries gives estimated stellar masses of 60 and 73 M⊙, respectively, for Barnard 30 and 35.No estimates of the gas or dust mass in either dark cloud could be found in the literature.The GES data for Barnard 30 and 35 have not previously been studied dynamically.
• The S Mon Cluster and the Spokes Cluster (NGC 2264), are two young clusters in the NGC 2264 H ii region.NGC 2264 is a highly substructured region elongated along a ∼10 pc NW-SE orientation.The Spokes cluster (sometimes called NGC 2264-C) is centred around the Class I protostar IRS2 in the south (Teixeira et al. 2006), while the S Mon cluster (sometimes referred to as the Christmas Tree cluster) is associated with the O7 binary S Mon (Sung et al. 2009) in the north.Both clusters are young, with ages of ∼1 Myr for the Spokes cluster and ∼2 Myr for the S Mon cluster (Venuti et al. 2019), and partially embedded.They each have approximate sizes of ∼1 pc in diameter.The total stellar population in NGC 2264 is estimated to be ∼1700 members, approximately equally divided between the two clusters (Venuti et al. 2019).This equates to a total mass of each cluster of 425 M⊙, when accounting for unresolved binaries.Oliver et al. (1996) estimated the molecular cloud that NGC 2264 is embedded within has a total mass of ∼28000 M⊙ (when scaled to our Gaia distance), but from molecular maps we estimate that only ∼10% of the molecular cloud's mass (∼3000 M⊙) is centred around each cluster and therefore relevant for its dynamics.We select members of the Spokes cluster as those with δ < 9.72 • and members of the S Mon cluster with δ > 9.72 • .The GES data for NGC 2264 has not previously been studied dynamically.
• ASCC 50 (Alessi 43) is a young group in the Vela T2 association (Pettersson & Reipurth 1994).The group was first detected by Kharchenko et al. (2005) 2018) estimate a mass of 45-81 M⊙, or 79 ± 23 M⊙ when accounting for binaries.The latter seems more consistent with the group's most massive member being the B3/5II star HD 74804 (that we estimate has a mass of 5.5 M⊙ from fitting its photometry to evolutionary models, Ekström et al. 2012), and therefore we use that mass here.The radius of the group is approximately 12 pc (Bonatto & Bica 2010).The GES data for Collinder 197 has not previously been studied dynamically.
• NGC 6530 (the Lagoon Nebula), is the group with the most GES spectra, with ∼650 members identified.The group is ∼1.5 Myr old (Bell et al. 2013), with a radius of ∼2 pc (Wright et al. 2019) at a distance of 1.24 kpc.The total stellar mass of NGC 6530, including binaries, has been estimated to be 3125 M⊙ (Wright et al. 2019), while the molecular cloud it is associated with has an estimated mass of 40,000 M⊙ (Takeuchi et al. 2010).The GES data for NGC 6530 was previously studied by Wright et al. (2019) and Wright & Parker (2019).
• NGC 2244 (the Rosette Nebula) is a 2 Myr (Bell et al. 2013) rich cluster containing >70 OB stars (Wang et al. 2008).The total stellar mass has not been well measured in the literature, but is approximately 1300 M⊙ based on a mass of 40-50 M⊙ for its most massive member, the O4V star HD 46223.Such a mass would be consistent with suggestions that the cluster contains ∼2000 young stars (Wang et al. 2008).The cluster has a radius of approximately 4 pc (Wang et al. 2008) at a distance of ∼1.37 kpc.The entire Rosette Molecular Cloud has been estimated to have a total mass of ∼10 5 M⊙ (Blitz & Thaddeus 1980), though NGC 2244 is not embedded within the cloud, which actually surrounds the cluster.The GES data for NGC 2244 have not previously been studied dynamically.
• NGC 2237 is a young star cluster projected against the periphery of the Rosette Nebula, but ∼130 pc behind it (at a distance of 1.49 kpc).It was observed as part of the GES observations of NGC 2244.It is estimated to be ∼2 Myr old (Wang et al. 2010) and contain 400-600 stars (Wang et al. 2010), suggesting a total mass, including binaries, of ∼250 M⊙ (though note that both age and mass estimates were derived on the assumption that NGC 2237 was at the same distance as NGC 2244, 10% closer than its Gaia EDR3 parallax implies).Its radius is ∼0.2 degrees or ∼5 pc at 1.49 kpc.The GES data for NGC 2237 have not previously been studied dynamically.

APPENDIX B: SPECTROSCOPIC MEMBERSHIP INDICATORS
Here we provide a summary of our young star membership selection process, including figures illustrating the process in Figures B1 and B2.For a star to be included as a young star in our catalogue it must pass the following tests: (i) The star must have either a lithium equivalent width greater than observed in stars of the same temperature in the 30-50 Myr cluster IC 2602 (Randich et al. 1997) or they must have Hα full width at zero intensity (FWZI) measurements greater than 4 Å (Bonito et al. 2013).
(ii) The star must not have a gravity-sensitive index γ > 1 and T ef f < 5600 K (which together indicate the star is a cool giant, Damiani et al. 2014).
(iii) The star must have a parallax within 2σ of the central value determined for the group from Gaussian fitting to the parallax dispersion.
Stars are not required to have valid RVs or PMs to be included in our overall sample of 2683 young stars, and therefore stars whose astrometry does not pass the Gaia RUWE ⩽ 1.4 quality cut have their astrometry discarded but they themselves, and their RVs, are not discarded.This decision was made to maximise the number of stars with reliable kinematic data in at least one dimension that we can use to constrain the kinematic properties of the groups and clusters studied.Our kinematic analysis is of course limited to stars with at least an RV or a PM.

APPENDIX C: PROPER MOTION VECTOR POINT DIAGRAMS
Here we provide in Figure C1 the proper motion "vector point" diagrams for all groups and clusters, showing the members of each group with valid proper motions.ASSC 50 Figure C1.PM "vector point" diagrams for the members of all groups with available proper motions.A reminder that the median PM uncertainty is 0.062 mas yr −1 , which in the vast majority of cases is smaller than the symbol size used.

Figure 3 .
Figure 3. Velocity distributions for young stars in ASSC 50 (190 stars with RVs, 155 stars with PMs) with our forward-modelled, partially-rotated velocity dispersion fits superimposed.The top panels show the two PM velocity distributions (blue histograms) with projected velocity dispersion fits (black lines), the bottom left panel shows the RV velocity distribution (blue histogram) with the velocity dispersion fit (black line), and the bottom right panel shows a 2D projection of the PM velocity distribution (black points) with the projected velocity dispersion fit as a blue 2D histogram.

Figure 4 .
Figure4.Velocity distributions for young stars in Cha I North (45 stars with RVs, 26 stars with PMs) with our forward-modelled, fully-rotated velocity dispersion fits superimposed.The top panels show the three velocity distributions (blue histograms) with projected velocity dispersion fits (black lines), while the bottom panels show the three 2D projections of the velocity distributions (black points) with the projected velocity dispersion fits as a blue 2D histogram.

Figure 6 .
Figure 6.Expansion gradients in the S Mon cluster of NGC 2264 for 200 stars with PMs and rotated with a position angle of 120 • panel) and the perpendicular axis (bottom panel).The red lines show the best-fitting expansion gradients, with 100 additional fits sampled from the posterior distribution shown in light red.The best-fitting gradient, with uncertainties are listed in each panel.

Figure 7 .Figure 8 .
Figure7.Relationship between the 3D velocity dispersion of all systems in our sample and their total stellar mass, with uncertainties shown for both.Kendall's rank correlation test reveals a correlation of τ = 0.507 and a p-value of 0.0035 indicating a strong positive correlation.A linear fit between the two quantities is shown in black, of the form σ 3D = 1.18+0.24−0.32 + 1.16 +0.61 −0.50 × (M cluster /1000 M ⊙ ) 0.89 +0.41 −0.30 km s −1 , with grey lines showing 100 randomly sampled fits from the posterior distribution of fitted gradients.

Figure 9 .
Figure9.Relationship between the virial ratio, α vir = σ 3D /σ vir , of all systems in our sample and their radius, with uncertainties shown for both quantities.Kendall's rank correlation test reveals a correlation of τ = 0.516 and a p-value of 0.0022 indicating a strong positive correlation.

Figure 10 .
Figure10.Relationship between the virial ratio, α vir = σ 3D /σ vir , of all systems in our sample and their stellar mass, with uncertainties shown for both.Kendall's rank correlation test reveals a correlation of τ = −0.302and a p-value of 0.081 indicating a weak inverse correlation.

Figure B1 .
Figure B1.Quantities derived from GES spectroscopy used for membership selection and parallax distributions for all groups.Left-hand panels show the gravity index γ, centre-left panels show the EW of the lithium 6708 Å line, and centre-right-hand panels show the FWZI of Hα, all plotted against effective temperature.In all three panels red circles show sources that pass our spectroscopic membership criteria and blue circles show sources that do not.The dashed lines show the thresholds used to identify giants in the upper-right corner of the γ-T ef f plot and to identify young stars above the thresholds in the EW(Li)-T ef f and FWZI(Hα)-T ef f plots.A typical error bar is shown in the top-left corner of each panel illustrating the typical uncertainties of 100 K in T ef f , 13 m Å for EW(Li), 0.011 for γ and 1.1 Åfor FWZI(Hα).The right-hand panel shows the parallax distribution for the likely members towards each group.The black histogram shows the distribution of spectroscopically-identified members, the blue line shows a Gaussian fit to this distribution, and the red histogram shows the distribution of the sources that fall within 2 standard deviations (also accounting for errors) of the median parallax and therefore that constitute our final sample.
Figure 5. Expansion gradients in the S Mon cluster of NGC 2264 for 200 stars with PMs and 196 stars with RVs.From top to bottom are, respectively, PM in RA plotted against RA, PM in declination plotted against Declination, and RV plotted against the distance.The red line shows the best-fitting 1D expansion gradient, with 100 additional fits sampled from the posterior distribution shown in light red.The best-fitting gradient, with uncertainties are listed in each panel.
Bonatto & Bica (2010)88)n and has a Gaia distance of 912 ± 3 pc.The age of the group has been estimated from ∼5 Myrs(Prisinzano et al. 2018) to ∼11.5 Myrs (Cantat-Gaudin et al.  2020).The most massive member is the O9V+B0V binary HD 75759, which would have a mass of 17-20 M⊙ and if on the main sequence an age ≲ 5−6 Myr(Ekström et al. 2012), therefore we adopt an age estimate of ∼5 Myr.If this were the most massive star in the group it would imply a group mass of 300-400 M⊙(Weidner & Kroupa 2006).Prisinzano et al. (2018)estimate a total mass of 50-86 M⊙, which when scaled to account for binarity gives a mass of 85 ± 23 M⊙.We compromise and use a mass of 200 ± 100 M⊙.The group has a radius of ∼0.25 degrees, or 4 pc at a distance of 912 pc.The GES data for ASCC 50 have not previously been studied dynamically.•Collinder197 is a group of young (∼6 Myr,Prisinzano et al. 2018) stars at a distance of 925 pc in Vela.The total mass of the group was estimated byBonatto & Bica (2010)to be 660 +102 −59 M⊙, thoughPrisinzano et al. (