<figure><img src="image_1.png" /></figure>

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

Nicholas J. Wright<sup>1,★</sup> R. D. Jeffries<sup>1</sup> R. J. Jackson,<sup>1</sup> G. G. Sacco,<sup>2</sup> Becky Arnold<sup>1</sup>

E. Franciosini<sup>1,2</sup> G. Gilmore,<sup>3</sup> A. Gonneau,<sup>3,4</sup> L. Morbidelli,<sup>2</sup> L. Prisinzano,<sup>5</sup> S. Randich<sup>2</sup> and Clare C. Worley<sup>3,6</sup>

Accepted 2024 July 19. Received 2024 July 19; in original form 2023 November 2

## ABSTRACT

We present the first large-scale 3D kinematic study of ~2700 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⁻<fnref n="1" /> (1D velocity dispersions of 0.35–4.3 km s⁻<fnref n="1" />). 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.

**Key words:** stars: formation – stars: kinematics and dynamics.

## 1 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 ‘associations’, 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 sub-clusters (Bonnell, Bate & Vine 2003; Arnold, Wright & Parker 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, it 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 time-scale 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

<footnote n="1"><sup>1</sup>Astrophysics Group, Keele University, Keele ST5 5BG, UK</footnote>

<footnote n="2"><sup>2</sup>INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi, 5, F-50125 Firenze, Italy</footnote>

<footnote n="3"><sup>3</sup>Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="4"><sup>4</sup>CNRS, Observatoire astronomique de Strasbourg, Université de Strasbourg, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="5"><sup>5</sup>INAF – Osservatorio Astronomico di Palermo, Piazza del Parlamento 1, I-90134 Palermo, Italy</footnote>

<footnote n="6"><sup>6</sup>School of Physical and Chemical Sciences – Te Kura Matū, University of Canterbury, Private Bag 4800, Christchurch 8140, New Zealand</footnote>

<footnote>*E-mail: n.j.wright@keele.ac.uk</footnote>
<table><caption>**Table 1.** Young clusters and associations studied in this work in order of increasing distance. Ages and extinctions are taken from the literature, while distances are calculated as described in Section 3.2 (note that distance uncertainties combine the observational errors and fitting errors, but do not include any systematic errors, though the Gaia EDR3 parallax zero point of −17 µas has been corrected for). N<sub>obs</sub> is the number of stars with spectroscopy from GES, N<sub>data</sub> is the number of stars with either EW(Li) or FWZI(H α) necessary to determine membership, N<sub>yso</sub> is the number of high-confidence young stars towards each group, N<sub>1D+</sub> is the number of objects with either a PM or a RV, and N<sub>3D</sub> is the number of young stars with 3D kinematics.</caption><thead><tr><th>Group</th><th>Age<br>(Myr)</th><th>Ref.</th><th>Distance<br>(pc)</th><th>A<sub>V</sub><br>(mag.)</th><th>Ref.</th><th>N<sub>obs</sub></th><th>N<sub>data</sub></th><th>N<sub>yso</sub></th><th>N<sub>1D+</sub></th><th>N<sub>3D</sub></th></tr></thead><tbody><tr><td>Rho Ophiuchus</td><td>3.0</td><td>GR21</td><td>136<sup>+1</sup><sub>-1</sub></td><td>1.0</td><td>W08</td><td>310</td><td>281</td><td>38</td><td>38</td><td>30</td></tr><tr><td>Cha I (South)</td><td>1.5</td><td>GA21</td><td>187<sup>+1</sup><sub>-1</sub></td><td>3.0</td><td>L07</td><td>345</td><td>305</td><td>48</td><td>48</td><td>28</td></tr><tr><td>Cha I (North)</td><td>1.5</td><td>GA21</td><td>191<sup>+1</sup><sub>-1</sub></td><td>3.0</td><td>L07</td><td>362</td><td>334</td><td>45</td><td>45</td><td>26</td></tr><tr><td>Gamma Vel</td><td>19.5</td><td>J17</td><td>334<sup>+1</sup><sub>-1</sub></td><td>0.13</td><td>J09</td><td>124</td><td>119</td><td>97</td><td>97</td><td>95</td></tr><tr><td>Vela OB2</td><td>14</td><td>A22</td><td>367<sup>+2</sup><sub>-2</sub></td><td>0.13</td><td>J09</td><td>989</td><td>916</td><td>73</td><td>73</td><td>72</td></tr><tr><td>25 Ori</td><td>19</td><td>F22</td><td>339<sup>+1</sup><sub>-1</sub></td><td>0.3</td><td>Z19</td><td>245</td><td>227</td><td>149</td><td>149</td><td>115</td></tr><tr><td>Barnard 30</td><td>2.4</td><td>K18</td><td>384<sup>+2</sup><sub>-2</sub></td><td>0.4</td><td>Z19</td><td>227</td><td>196</td><td>57</td><td>57</td><td>36</td></tr><tr><td>λ Ori</td><td>10</td><td>B13</td><td>389<sup>+1</sup><sub>-1</sub></td><td>0.4</td><td>Z19</td><td>380</td><td>329</td><td>144</td><td>144</td><td>111</td></tr><tr><td>Barnard 35</td><td>2.6</td><td>K18</td><td>390<sup>+1</sup><sub>-1</sub></td><td>0.4</td><td>Z19</td><td>226</td><td>183</td><td>57</td><td>57</td><td>39</td></tr><tr><td>S Mon Cluster (NGC 2264)</td><td>2</td><td>V19</td><td>690<sup>+2</sup><sub>-2</sub></td><td>0.71</td><td>D08</td><td>863</td><td>780</td><td>246</td><td>244</td><td>192</td></tr><tr><td>Spokes Cluster (NGC 2264)</td><td>1</td><td>V19</td><td>691<sup>+2</sup><sub>-2</sub></td><td>0.71</td><td>D08</td><td>992</td><td>881</td><td>281</td><td>279</td><td>188</td></tr><tr><td>ASCC 50 / RCW33</td><td>5</td><td>This work</td><td>912<sup>+3</sup><sub>-3</sub></td><td>0.71</td><td>K13</td><td>1224</td><td>1131</td><td>193</td><td>190</td><td>155</td></tr><tr><td>Collinder 197</td><td>5</td><td>B10</td><td>925<sup>+7</sup><sub>-7</sub></td><td>1.05</td><td>B10</td><td>408</td><td>387</td><td>140</td><td>140</td><td>110</td></tr><tr><td>NGC 6530</td><td>1.5</td><td>B13</td><td>1242<sup>+33</sup><sub>-32</sub></td><td>1.1</td><td>S00</td><td>1982</td><td>1687</td><td>644</td><td>635</td><td>455</td></tr><tr><td>NGC 2244</td><td>2</td><td>B13</td><td>1368<sup>+24</sup><sub>-23</sub></td><td>1.5</td><td>M19</td><td>248</td><td>216</td><td>109</td><td>109</td><td>88</td></tr><tr><td>NGC 2237</td><td>2</td><td>W10</td><td>1494<sup>+44</sup><sub>-41</sub></td><td>1.5</td><td>M19</td><td>54</td><td>54</td><td>23</td><td>23</td><td>20</td></tr><tr><td>Trumpler 14</td><td>2.5</td><td>H12, D17</td><td>2211<sup>+16</sup><sub>-16</sub></td><td>2.0</td><td>H12</td><td>631</td><td>363</td><td>162</td><td>152</td><td>82</td></tr><tr><td>Trumpler 16</td><td>2</td><td>H12, D21</td><td>2260<sup>+23</sup><sub>-22</sub></td><td>2.0</td><td>H12</td><td>1269</td><td>616</td><td>208</td><td>202</td><td>100</td></tr><tr><td>Totals</td><td></td><td></td><td></td><td></td><td></td><td>10,879</td><td>9,005</td><td>2,683</td><td>2,651</td><td>1,914</td></tr></tbody></table>

Notes. References: A22 (Armstrong et al. 2022), B10 (Bonatto & Bica 2010), B13 (Bell et al. 2013), D08 (Dahm 2008), D17 (Damiani et al. 2017), D21 (Dias et al. 2021), F22 (Franciosini et al. 2022b), GA21 (Galli et al. 2021), GR21 (Grasser et al. 2021), H12 (Hur, Sung & Bessell 2012), J09 (Jeffries et al. 2009), J17 (Jeffries et al. 2017), K13 (Kharchenko et al. 2013), K18 (Kounkel et al. 2018), L07 (Luhman 2007), M19 (Mužić et al. 2019), S00 (Sung, Chun & Bessell 2000), V19 (Venuti, Damiani & Prisinzano 2019), W08 (Wilking, Gagné & Allen 2008), W10 (Wang et al. 2010), and Z19 (Zari, Brown & de Zeeuw 2019).

dispersing the residual gas and its associated gravitational potential (Tutukov 1978; Lada, Margulis & Dearborn 1984). However, recent studies have questioned the role that residual gas expulsion may play (Dale, Ercolano & Bonnell 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, Kruijssen & Rix 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, Wright & Jeffries 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 Fibre Large Array Multi Element Spectrograph (FLAMES; Pasquini et al. 2002) fibred spectrograph at the Very Large Telescope (VLT). 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⁵ 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 spectroscopy 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, which is 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.

## 2 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.

### 2.1 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.<fnref n="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 data-sets 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 per cent 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 per cent of sources). We calibrated a new relationship between T$_{eff}$ 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, Finkbeiner & Davis (1998) and Daniëlski et al. (2018), and the tables of intrinsic colour from Pecaut & Mamajek (2016), prioritizing 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 per cent of sources. For the remaining 38

(0.3 per cent) of sources this was impossible either due to a lack of photometry or the available photometry being inconsistent with group membership.

RVs are determined using a combination of a cross-correlation 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 \alpha$ 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 gravity-sensitive spectral index $\gamma$ 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 per cent) have effective temperatures, 8668 (80 per cent) have RVs, 7366 (68 per cent) have lithium equivalent widths, and 9321 (86 per cent) have $\gamma$ indices.

### 2.2 Gaia EDR3 astrometry and photometry

We use astrometry and photometry from Gaia (Gaia Collaboration 2016) and EDR3 (Gaia Collaboration 2021), which contains parallaxes, $\sigma$, and PMs, $\mu_\alpha$, and $\mu_\delta$, 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 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 per cent. The affect is smaller for all other groups and <2 per cent for groups within 1 kpc.

GES sources were cross-matched with Gaia EDR3 using a radius of 1 arcsec, which resulted in a total of 11 179 matches (99.4 per cent). Gaia astrometry and photometry were then filtered based on the criteria outlined in the Gaia data release papers and technical notes. For the astrometry, we required the 're-normalized 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 problems.<fnref n="2" /> Applying this cut removed

<footnote n="2">²Note that we do not discard the GES RVs for stars failing the RUWE cut, even though they might be binaries, because we do not know the fraction of binary stars that such a cut would discard and would prefer to model the effects of binarity on the velocity dispersion (see Section 4.2) rather than exclude them.</footnote>

<footnote n="1">¹http://ges.roe.ac.uk</footnote>
astrometry for 1407 sources (12.6 per cent), 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 do not 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 \times 5)$ *Gaia* covariance matrix whenever propagating uncertainties, and consider all errors to follow a normal distribution (Arenou et al. 2018).

# 3 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 \alpha$ 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 $\gamma$ (Damiani et al. 2014), and the full width at zero intensity (FWZI) of the $H \alpha$ 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.

## 3.1 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 per cent).

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 \alpha$) measurements, and 1096 of these (41.2 per cent) have FWZI($H \alpha$) 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.

Figs B1 and B2 show both EW(Li) and FWZI($H \alpha$) as a function of $T_{\text{eff}}$ for all stars in the 18 groups studied.

## 3.2 Separating background giants

To separate young stars from background giants (which sometimes show lithium in their atmospheres) we use the surface gravity-sensitive index $\gamma$ (Damiani et al. 2014), which is capable of distinguishing between main-sequence stars, pre-main-sequence stars, and giants (cool giants with $T_{\text{eff}} < 5600$ K have $\gamma > 1$, Damiani et al. 2014). Of our 3480 current candidate members, 3326 have $\gamma$ measurements, and of these 413 (12.4 per cent) have $\gamma$ 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 $\gamma$ measurements, since the probability of these being giants is relatively low (~12 per cent) and will be removed by the next step. Figs B1 and B2 show $\gamma$ as a function of $T_{\text{eff}}$ 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 non-uniform 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\sigma$ 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 $\sigma = 0.5-0.9$ mas (see Fig. B2) that we isolated by ensuring that the parallax distribution was fitted to the main peak at $\sigma \sim 1.1$ mas that represents Col. 197.

This cut excludes 384 (15.2 per cent) 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 per cent) 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. Figs B1 and B2 show both $\gamma$ as a function of $T_{\text{eff}}$ and the parallax distribution for all candidate members.

## 3.3 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 per cent), 503 only have RVs (18.7 per cent), 234 only have PMs (8.7 per cent), and 32 lack both PMs and RVs (1.2 per cent), 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 maximize the numbers of stars with reliable kinematic data in at least one dimension that
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Velocity distributions for young stars in Barnard 30 (53 stars with RVs, 38 stars with PMs).</figcaption></figure>

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⁻¹, with some as low as 0.25 km s⁻¹. The median PM uncertainty is 0.062 mas yr⁻¹, which equates to between 0.04 and 0.80 km s⁻¹ 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.

# 4 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.

## 4.1 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 Fig. 1 for young stars in Barnard 30. The velocity distributions in all three dimensions show evidence for substructure, with clear bimodal distributions in $\mu_\alpha$ and $\mu_\delta$. 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 ($\mu_\alpha$ and $\mu_\delta$), most likely because the cluster is relatively old (~19.5 Myr) and therefore more dynamically processed.

## 4.2 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 parametrized 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, $\theta$, best explain the observations, $d$. In Bayes’s theorem this is known as the posterior distribution, $P(\theta|d) = P(d|\theta)P(\theta)/P(d)$, where $P(d|\theta)$ is the likelihood model, $P(\theta)$ are the priors (which includes our a priori knowledge about the model parameters) and $P(d)$ is a normalizing 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 Gaussians.<fnref n="3" /> The first is a trivariate Gaussian aligned to the observed coordinate system ($\alpha, \delta$, and along the line of sight), hereafter known as the *unrotated* model. This model is used because it represents one of the most commonly used fitting methods employed in the literature.<fnref n="4" /> The second model replaces the plane of the sky components with a rotated bivariate Gaussian, alongside the third, line of sight, component (the *partially rotated* model). The third model utilizes 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

<footnote n="3">³While many of the distributions appear non-Gaussian, a Gaussian distribution still represents a reasonable approximation to the underlying distribution.</footnote>

<footnote n="4">⁴We note that the results obtained are independent of the, relatively arbitrary, choice of coordinate system, which we verified by performing the same fit in the Galactic Cartesian coordinate system (UVW), obtaining consistent results.</footnote>
3 rotation angles, respectively). In addition, we utilize 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_{\text{field}}$ 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 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 is in binary systems.<fnref n="5" /> The fraction of binary stars in young stellar groups is poorly constrained and so we set the binary fraction to be 46 per cent, 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 $\gamma = -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 log$_{10}$ 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 photocentre velocity,<fnref n="6" /> 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–100 km s⁻¹ for the velocity dispersions, −100 to +100 km s⁻¹ for the central velocities, 0–0.2 for the field star component and 0–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 16th and 84th percentiles used for the 1σ uncertainties.

## 4.3 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 to 6.1 km s⁻¹, with the majority around ~1 km s⁻¹. Example forward-modelled fits are shown in Fig. 2 for the unrotated model for the Gamma Velorum cluster and in Fig. 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⁻¹. Fig. 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 per cent of the unrotated model fits agree with the partially rotated model fits within 1σ and 94 per cent 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 per cent within 1σ and 60 per cent 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 per cent smaller than the unrotated 3D velocity dispersions. This implies that velocity dispersions previously calculated without considering the rotation of the velocity ellipsoid may have overestimated both the 1D and 3D velocity dispersions by ~10 per cent.

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 per cent for the unrotated models, 89 per cent for the partially rotated models, and 100 per cent for the fully rotated models. Limiting this com-

<footnote n="5">⁵We do not consider triple systems because their properties are poorly constrained and are typically hierarchical, meaning that the third star is usually on a wide, long-period orbit that does not introduce a large velocity offset.</footnote>

<footnote n="6">⁶This compensates for the fact that for high mass-ratio binaries some of the light from the secondary contributes to the spectral features used to measure the RV, and thus the measured RV will be intermediate between that of the two stars. From simulations, we found that this reduces the broadening of the RV distribution due to binaries by ~25 per cent.</footnote>
<table><caption>**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. $\mu_1$ and $\mu_2$ are the velocity dispersions along the semimajor and semiminor axes of the partially rotated velocity ellipsoid, the position angle of which is $\theta$.</caption><thead><tr><th rowspan="2">Group</th><th colspan="4">Unrotated velocity dispersions</th><th colspan="4">Partially rotated velocity dispersions</th></tr><tr><th>RV [km s<sup>-1</sup>]</th><th>μ<sub>α</sub> [km s<sup>-1</sup>]</th><th>μ<sub>δ</sub> [km s<sup>-1</sup>]</th><th>σ<sub>3D</sub> [km s<sup>-1</sup>]</th><th>RV [km s<sup>-1</sup>]</th><th>μ<sub>1</sub> [km s<sup>-1</sup>]</th><th>μ<sub>2</sub> [km s<sup>-1</sup>]</th><th>θ [deg]</th><th>σ<sub>3D</sub> [km s<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Rho Oph</td><td>1.087<sup>+0.080</sup><sub>-0.111</sub></td><td>0.67<sup>+0.07</sup><sub>-0.07</sub></td><td>0.96<sup>+0.09</sup><sub>-0.08</sub></td><td>1.60<sup>+0.08</sup><sub>-0.09</sub></td><td>1.147<sup>+0.103</sup><sub>-0.100</sub></td><td>0.96<sup>+0.09</sup><sub>-0.08</sub></td><td>0.74<sup>+0.06</sup><sub>-0.07</sub></td><td>91<sup>+24</sup><sub>-26</sub></td><td>1.67<sup>+0.10</sup><sub>-0.08</sub></td></tr><tr><td>Cha I (south)</td><td>0.859<sup>+0.147</sup><sub>-0.165</sub></td><td>0.40<sup>+0.08</sup><sub>-0.06</sub></td><td>0.52<sup>+0.04</sup><sub>-0.06</sub></td><td>1.08<sup>+0.13</sup><sub>-0.13</sub></td><td>0.844<sup>+0.105</sup><sub>-0.099</sub></td><td>0.54<sup>+0.05</sup><sub>-0.06</sub></td><td>0.52<sup>+0.06</sup><sub>-0.05</sub></td><td>92<sup>+30</sup><sub>-30</sub></td><td>1.13<sup>+0.09</sup><sub>-0.08</sub></td></tr><tr><td>Cha I (north)</td><td>0.649<sup>+0.037</sup><sub>-0.056</sub></td><td>0.54<sup>+0.07</sup><sub>-0.05</sub></td><td>0.76<sup>+0.07</sup><sub>-0.08</sub></td><td>1.13<sup>+0.06</sup><sub>-0.06</sub></td><td>0.624<sup>+0.052</sup><sub>-0.051</sub></td><td>0.68<sup>+0.03</sup><sub>-0.05</sub></td><td>0.48<sup>+0.02</sup><sub>-0.03</sub></td><td>31<sup>+6</sup><sub>-5</sub></td><td>1.05<sup>+0.04</sup><sub>-0.05</sub></td></tr><tr><td>Gamma Vel</td><td>0.344<sup>+0.060</sup><sub>-0.065</sub></td><td>0.36<sup>+0.03</sup><sub>-0.03</sub></td><td>0.37<sup>+0.04</sup><sub>-0.03</sub></td><td>0.62<sup>+0.05</sup><sub>-0.04</sub></td><td>0.353<sup>+0.024</sup><sub>-0.024</sub></td><td>0.36<sup>+0.02</sup><sub>-0.02</sub></td><td>0.35<sup>+0.02</sup><sub>-0.02</sub></td><td>176<sup>+6</sup><sub>-7</sub></td><td>0.61<sup>+0.03</sup><sub>-0.02</sub></td></tr><tr><td>Vela OB2</td><td>1.490<sup>+0.109</sup><sub>-0.131</sub></td><td>0.54<sup>+0.05</sup><sub>-0.04</sub></td><td>1.17<sup>+0.07</sup><sub>-0.06</sub></td><td>1.97<sup>+0.10</sup><sub>-0.10</sub></td><td>1.411<sup>+0.106</sup><sub>-0.119</sub></td><td>1.09<sup>+0.08</sup><sub>-0.08</sub></td><td>0.51<sup>+0.07</sup><sub>-0.07</sub></td><td>35<sup>+6</sup><sub>-6</sub></td><td>1.86<sup>+0.10</sup><sub>-0.10</sub></td></tr><tr><td>25 Ori</td><td>0.501<sup>+0.045</sup><sub>-0.048</sub></td><td>0.40<sup>+0.05</sup><sub>-0.04</sub></td><td>0.39<sup>+0.04</sup><sub>-0.04</sub></td><td>0.75<sup>+0.05</sup><sub>-0.04</sub></td><td>0.501<sup>+0.058</sup><sub>-0.061</sub></td><td>0.39<sup>+0.02</sup><sub>-0.02</sub></td><td>0.32<sup>+0.04</sup><sub>-0.03</sub></td><td>193<sup>+11</sup><sub>-11</sub></td><td>0.71<sup>+0.05</sup><sub>-0.04</sub></td></tr><tr><td>Barnard 30</td><td>0.811<sup>+0.258</sup><sub>-0.206</sub></td><td>1.40<sup>+0.34</sup><sub>-0.22</sub></td><td>1.67<sup>+0.34</sup><sub>-0.26</sub></td><td>2.33<sup>+0.38</sup><sub>-0.20</sub></td><td>0.831<sup>+0.121</sup><sub>-0.150</sub></td><td>1.94<sup>+0.57</sup><sub>-0.52</sub></td><td>0.85<sup>+0.25</sup><sub>-0.19</sub></td><td>47<sup>+14</sup><sub>-12</sub></td><td>2.27<sup>+0.53</sup><sub>-0.41</sub></td></tr><tr><td>λ Ori</td><td>0.717<sup>+0.052</sup><sub>-0.048</sub></td><td>0.52<sup>+0.06</sup><sub>-0.04</sub></td><td>0.33<sup>+0.04</sup><sub>-0.03</sub></td><td>0.95<sup>+0.06</sup><sub>-0.04</sub></td><td>0.696<sup>+0.038</sup><sub>-0.036</sub></td><td>0.45<sup>+0.04</sup><sub>-0.04</sub></td><td>0.30<sup>+0.04</sup><sub>-0.03</sub></td><td>136<sup>+7</sup><sub>-7</sub></td><td>0.88<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>Barnard 35</td><td>1.053<sup>+0.120</sup><sub>-0.105</sub></td><td>0.85<sup>+0.05</sup><sub>-0.06</sub></td><td>0.46<sup>+0.06</sup><sub>-0.06</sub></td><td>1.43<sup>+0.10</sup><sub>-0.08</sub></td><td>1.058<sup>+0.122</sup><sub>-0.109</sub></td><td>0.83<sup>+0.06</sup><sub>-0.05</sub></td><td>0.42<sup>+0.12</sup><sub>-0.04</sub></td><td>97<sup>+29</sup><sub>-5</sub></td><td>1.41<sup>+0.12</sup><sub>-0.08</sub></td></tr><tr><td>S Mon Cluster</td><td>1.993<sup>+0.146</sup><sub>-0.117</sub></td><td>0.59<sup>+0.06</sup><sub>-0.06</sub></td><td>0.58<sup>+0.05</sup><sub>-0.04</sub></td><td>2.16<sup>+0.14</sup><sub>-0.11</sub></td><td>1.903<sup>+0.083</sup><sub>-0.080</sub></td><td>0.88<sup>+0.24</sup><sub>-0.25</sub></td><td>0.64<sup>+0.12</sup><sub>-0.12</sub></td><td>172<sup>+20</sup><sub>-19</sub></td><td>2.19<sup>+0.14</sup><sub>-0.11</sub></td></tr><tr><td>Spokes Cluster</td><td>2.106<sup>+0.128</sup><sub>-0.131</sub></td><td>1.48<sup>+0.07</sup><sub>-0.07</sub></td><td>0.92<sup>+0.07</sup><sub>-0.08</sub></td><td>2.73<sup>+0.11</sup><sub>-0.11</sub></td><td>2.011<sup>+0.078</sup><sub>-0.084</sub></td><td>1.70<sup>+0.09</sup><sub>-0.08</sub></td><td>0.56<sup>+0.13</sup><sub>-0.12</sub></td><td>93<sup>+1</sup><sub>-1</sub></td><td>2.69<sup>+0.09</sup><sub>-0.08</sub></td></tr><tr><td>ASSC 50</td><td>0.573<sup>+0.037</sup><sub>-0.030</sub></td><td>1.16<sup>+0.09</sup><sub>-0.07</sub></td><td>1.05<sup>+0.08</sup><sub>-0.06</sub></td><td>1.67<sup>+0.08</sup><sub>-0.06</sub></td><td>0.554<sup>+0.044</sup><sub>-0.042</sub></td><td>1.69<sup>+0.09</sup><sub>-0.10</sub></td><td>0.16<sup>+0.03</sup><sub>-0.03</sub></td><td>56<sup>+3</sup><sub>-3</sub></td><td>1.79<sup>+0.09</sup><sub>-0.09</sub></td></tr><tr><td>Col. 197</td><td>0.655<sup>+0.035</sup><sub>-0.040</sub></td><td>1.15<sup>+0.10</sup><sub>-0.09</sub></td><td>1.19<sup>+0.12</sup><sub>-0.09</sub></td><td>1.78<sup>+0.11</sup><sub>-0.08</sub></td><td>0.657<sup>+0.032</sup><sub>-0.030</sub></td><td>1.12<sup>+0.11</sup><sub>-0.12</sub></td><td>0.77<sup>+0.10</sup><sub>-0.07</sub></td><td>47<sup>+7</sup><sub>-8</sub></td><td>1.51<sup>+0.11</sup><sub>-0.09</sub></td></tr><tr><td>NGC 6530</td><td>2.242<sup>+0.218</sup><sub>-0.222</sub></td><td>2.51<sup>+0.14</sup><sub>-0.13</sub></td><td>1.80<sup>+0.09</sup><sub>-0.09</sub></td><td>3.82<sup>+0.17</sup><sub>-0.15</sub></td><td>2.320<sup>+0.251</sup><sub>-0.264</sub></td><td>2.63<sup>+0.12</sup><sub>-0.13</sub></td><td>2.63<sup>+0.14</sup><sub>-0.13</sub></td><td>100<sup>+20</sup><sub>-19</sub></td><td>4.38<sup>+0.18</sup><sub>-0.17</sub></td></tr><tr><td>NGC 2244</td><td>2.658<sup>+0.018</sup><sub>-0.015</sub></td><td>1.69<sup>+0.27</sup><sub>-0.36</sub></td><td>1.78<sup>+0.23</sup><sub>-0.28</sub></td><td>3.62<sup>+0.18</sup><sub>-0.20</sub></td><td>2.698<sup>+0.063</sup><sub>-0.062</sub></td><td>1.26<sup>+0.13</sup><sub>-0.13</sub></td><td>1.23<sup>+0.13</sup><sub>-0.14</sub></td><td>102<sup>+11</sup><sub>-11</sub></td><td>3.22<sup>+0.09</sup><sub>-0.09</sub></td></tr><tr><td>NGC 2237</td><td>0.492<sup>+0.071</sup><sub>-0.063</sub></td><td>1.00<sup>+0.41</sup><sub>-0.21</sub></td><td>1.10<sup>+0.37</sup><sub>-0.26</sub></td><td>1.56<sup>+0.43</sup><sub>-0.20</sub></td><td>0.501<sup>+0.068</sup><sub>-0.068</sub></td><td>1.14<sup>+0.20</sup><sub>-0.27</sub></td><td>1.12<sup>+0.21</sup><sub>-0.26</sub></td><td>46<sup>+32</sup><sub>-32</sub></td><td>1.67<sup>+0.21</sup><sub>-0.23</sub></td></tr><tr><td>Trumpler 14</td><td>6.084<sup>+0.234</sup><sub>-0.192</sub></td><td>3.07<sup>+0.15</sup><sub>-0.16</sub></td><td>3.13<sup>+0.19</sup><sub>-0.18</sub></td><td>7.50<sup>+0.22</sup><sub>-0.18</sub></td><td>6.025<sup>+0.193</sup><sub>-0.161</sub></td><td>3.12<sup>+0.18</sup><sub>-0.16</sub></td><td>2.85<sup>+0.13</sup><sub>-0.13</sub></td><td>81<sup>+8</sup><sub>-8</sub></td><td>7.36<sup>+0.19</sup><sub>-0.15</sub></td></tr><tr><td>Trumpler 16</td><td>2.879<sup>+0.146</sup><sub>-0.139</sub></td><td>3.26<sup>+0.30</sup><sub>-0.27</sub></td><td>2.29<sup>+0.21</sup><sub>-0.19</sub></td><td>4.91<sup>+0.25</sup><sub>-0.20</sub></td><td>2.724<sup>+0.181</sup><sub>-0.166</sub></td><td>2.97<sup>+0.23</sup><sub>-0.28</sub></td><td>1.80<sup>+0.13</sup><sub>-0.14</sub></td><td>163<sup>+5</sup><sub>-5</sub></td><td>4.41<sup>+0.20</sup><sub>-0.21</sub></td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Velocity distributions (blue histograms) for young stars in the Gamma Velorum cluster (97 stars with RVs, 95 stars with PMs) with our forward modelled, non-rotated velocity dispersion fits superimposed (black lines).</figcaption></figure>

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 per cent for the unrotated models and 67 per cent 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\sigma$), 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 \pm 0.18$ and $0.87 \pm 0.24$ km s⁻¹ 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\sigma$, 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.
<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

<table><caption>**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_3$ are the velocity dispersions along the three axes, with rotation angles of $\theta$, $\psi$, and $\phi$.</caption><thead><tr><th rowspan="3">Group</th><th colspan="7">Fully rotated velocity dispersions</th></tr><tr><th>v<sub>1</sub></th><th>v<sub>2</sub></th><th>v<sub>3</sub></th><th>θ</th><th>ψ</th><th>φ</th><th>σ<sub>3D</sub></th></tr><tr><th>[km s<sup>-1</sup>]</th><th>[km s<sup>-1</sup>]</th><th>[km s<sup>-1</sup>]</th><th>[deg]</th><th>[deg]</th><th>[deg]</th><th>[km s<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Rho Oph</td><td>0.85<sup>+0.08</sup><sub>-0.07</sub></td><td>0.35<sup>+0.07</sup><sub>-0.05</sub></td><td>0.91<sup>+0.09</sup><sub>-0.07</sub></td><td>39<sup>+14</sup><sub>-18</sub></td><td>40<sup>+15</sup><sub>-15</sub></td><td>3<sup>+1</sup><sub>-1</sub></td><td>1.30<sup>+0.09</sup><sub>-0.06</sub></td></tr><tr><td>Cha I (south)</td><td>1.01<sup>+0.24</sup><sub>-0.17</sub></td><td>0.39<sup>+0.04</sup><sub>-0.02</sub></td><td>0.06<sup>+0.02</sup><sub>-0.01</sub></td><td>33<sup>+19</sup><sub>-15</sub></td><td>54<sup>+22</sup><sub>-19</sub></td><td>53<sup>+18</sup><sub>-16</sub></td><td>1.09<sup>+0.23</sup><sub>-0.15</sub></td></tr><tr><td>Cha I (north)</td><td>1.32<sup>+0.11</sup><sub>-0.18</sub></td><td>0.25<sup>+0.09</sup><sub>-0.04</sub></td><td>0.42<sup>+0.08</sup><sub>-0.06</sub></td><td>45<sup>+20</sup><sub>-15</sub></td><td>54<sup>+14</sup><sub>-10</sub></td><td>83<sup>+18</sup><sub>-20</sub></td><td>1.41<sup>+0.12</sup><sub>-0.16</sub></td></tr><tr><td>Vela OB2</td><td>1.40<sup>+0.10</sup><sub>-0.10</sub></td><td>0.60<sup>+0.07</sup><sub>-0.07</sub></td><td>0.80<sup>+0.11</sup><sub>-0.13</sub></td><td>38<sup>+8</sup><sub>-8</sub></td><td>40<sup>+6</sup><sub>-6</sub></td><td>43<sup>+10</sup><sub>-10</sub></td><td>1.72<sup>+0.10</sup><sub>-0.10</sub></td></tr><tr><td>Barnard 30</td><td>0.45<sup>+0.06</sup><sub>-0.05</sub></td><td>0.66<sup>+0.06</sup><sub>-0.06</sub></td><td>1.54<sup>+0.20</sup><sub>-0.28</sub></td><td>19<sup>+29</sup><sub>-13</sub></td><td>24<sup>+31</sup><sub>-17</sub></td><td>46<sup>+11</sup><sub>-10</sub></td><td>1.73<sup>+0.18</sup><sub>-0.24</sub></td></tr></tbody></table>

## 4.4 Virial state

To measure the virial state of each group, we use the virial equation, which in its 3D form is given by

$$ \sigma_{3D, vir}^2 = \frac{G M}{2 r_{vir}}, \qquad (1) $$

where $\sigma_{3D}$ is the 3D velocity dispersion, $G$ is the gravitational constant, $M$ is the mass, and $r_{vir}$ is the radius. We substitute the parameter $\eta = 6r_{vir}/r_{eff}$, where $r_{eff}$ is the effective (or half-light) radius in projection, to get

$$ \sigma_{3D, vir}^2 = \frac{3 G M_{vir}}{\eta r_{eff}}. \qquad (2) $$

The parameter $\eta$ can be derived from the power-law index, $\gamma$, of an Elson, Fall & Freeman (1987, hereafter EFF) surface brightness profile (see e.g. Portegies Zwart, McMillan & Gieles 2010), from which the effective radius, $r_{eff}$, can also be measured. Some studies assume a standard value of $\eta = 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 data-sets mentioned in Section 2 to fit EFF profiles for all our groups, deriving $r_{eff}$ and $\gamma$, with uncertainties calculated using a bootstrapping process. The power-law indexes, $\gamma$, and therefore also the $\eta$ values, are not well constrained for these groups, with values ranging from $\eta = 6-11$, with varying levels of precision (see Table 4). There are no direct measurements of $\eta$ or $\gamma$ 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_{eff}$ 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
<figure><img src="image_5.png" /><figcaption>**Figure 4.** 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.</figcaption></figure>

<table><caption>**Table 4.** Velocity dispersions, stellar and gas masses, effective radii and virial states for all groups. The 3D velocity dispersion, $\sigma_{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 $\alpha = \sigma_{3D}/\sigma_{3D, vir}$, again using either the stellar mass or the sum of the stellar and gas masses.</caption><thead><tr><th rowspan="2">Group</th><th rowspan="2">σ<sub>3D</sub><br>[km s<sup>-1</sup>]</th><th colspan="2">M [M<sub>☉</sub>]</th><th rowspan="2">r<sub>eff</sub><br>[pc]</th><th rowspan="2">η</th><th colspan="2">σ<sub>3D, vir</sub> [km s<sup>-1</sup>]</th><th colspan="2">α</th></tr><tr><th>(stars)</th><th>(gas)</th><th>(stars)</th><th>(stars + gas)</th><th>(stars)</th><th>(stars + gas)</th></tr></thead><tbody><tr><td>Rho Oph</td><td>1.30<sup>+0.09</sup><sub>-0.06</sub></td><td>106</td><td>1750</td><td>0.69<sup>+0.08</sup><sub>-0.08</sub></td><td>10.1<sup>+0.3</sup><sub>-0.3</sub></td><td>0.44<sup>+0.01</sup><sub>-0.04</sub></td><td>1.86<sup>+0.06</sup><sub>-0.16</sub></td><td>2.93<sup>+0.38</sup><sub>-0.13</sub></td><td>0.70<sup>+0.09</sup><sub>-0.03</sub></td></tr><tr><td>Cha I (south)</td><td>1.09<sup>+0.23</sup><sub>-0.15</sub></td><td>59</td><td>1000</td><td>0.4<sup>+0.1</sup><sub>-0.09</sub></td><td>6.0<sup>+0.1</sup><sub>-3.5</sub></td><td>0.56<sup>+0.05</sup><sub>-0.16</sub></td><td>2.39<sup>+0.19</sup><sub>-0.68</sub></td><td>1.93<sup>+1.00</sup><sub>-0.24</sub></td><td>0.46<sup>+0.24</sup><sub>-0.06</sub></td></tr><tr><td>Cha I (north)</td><td>1.41<sup>+0.12</sup><sub>-0.16</sub></td><td>54</td><td>1000</td><td>0.42<sup>+0.15</sup><sub>-0.19</sub></td><td>6.0<sup>+0.5</sup><sub>-3.5</sub></td><td>0.53<sup>+0.01</sup><sub>-0.20</sub></td><td>2.33<sup>+0.01</sup><sub>-0.88</sub></td><td>2.67<sup>+1.69</sup><sub>-0.01</sub></td><td>0.60<sup>+0.38</sup><sub>-0.01</sub></td></tr><tr><td>Gamma Vel</td><td>0.61<sup>+0.03</sup><sub>-0.02</sub></td><td>152</td><td>—</td><td>1.88<sup>+0.27</sup><sub>-0.23</sub></td><td>10.6<sup>+1.0</sup><sub>-3.6</sub></td><td>0.31<sup>+0.10</sup><sub>-0.01</sub></td><td>—</td><td>1.96<sup>+0.08</sup><sub>-0.47</sub></td><td>—</td></tr><tr><td>Vela OB2</td><td>1.72<sup>+0.10</sup><sub>-0.10</sub></td><td>1285</td><td>—</td><td>23.5<sup>+1.1</sup><sub>-2.0</sub></td><td>9.5<sup>+0.1</sup><sub>-0.1</sub></td><td>0.27<sup>+0.02</sup><sub>-0.01</sub></td><td>—</td><td>6.32<sup>+0.25</sup><sub>-0.58</sub></td><td>—</td></tr><tr><td>25 Ori</td><td>0.71<sup>+0.05</sup><sub>-0.04</sub></td><td>400</td><td>—</td><td>2.55<sup>+0.3</sup><sub>-0.18</sub></td><td>9.8<sup>+0.5</sup><sub>-1.9</sub></td><td>0.46<sup>+0.02</sup><sub>-0.04</sub></td><td>—</td><td>1.56<sup>+0.18</sup><sub>-0.11</sub></td><td>—</td></tr><tr><td>Barnard 30</td><td>1.73<sup>+0.18</sup><sub>-0.24</sub></td><td>60</td><td>—</td><td>1.75<sup>+0.8</sup><sub>-1.71</sub></td><td>8.7<sup>+0.1</sup><sub>-3.5</sub></td><td>0.23<sup>+0.01</sup><sub>-0.09</sub></td><td>—</td><td>7.68<sup>+5.74</sup><sub>-0.10</sub></td><td>—</td></tr><tr><td>λ Ori</td><td>0.88<sup>+0.04</sup><sub>-0.04</sub></td><td>650</td><td>—</td><td>0.23<sup>+0.29</sup><sub>-0.96</sub></td><td>8.7<sup>+1.1</sup><sub>-3.9</sub></td><td>2.06<sup>+0.01</sup><sub>-1.39</sub></td><td>—</td><td>0.43<sup>+0.92</sup><sub>-0.01</sub></td><td>—</td></tr><tr><td>Barnard 35</td><td>1.41<sup>+0.12</sup><sub>-0.08</sub></td><td>73</td><td>—</td><td>4.66<sup>+0.8</sup><sub>-0.93</sub></td><td>9.5<sup>+0.1</sup><sub>-3.5</sub></td><td>0.15<sup>+0.07</sup><sub>-0.01</sub></td><td>—</td><td>9.64<sup>+0.42</sup><sub>-3.04</sub></td><td>—</td></tr><tr><td>S Mon Cluster</td><td>2.19<sup>+0.14</sup><sub>-0.11</sub></td><td>425</td><td>3000</td><td>0.17<sup>+0.05</sup><sub>-0.09</sub></td><td>6.0<sup>+1.2</sup><sub>-3.5</sub></td><td>2.33<sup>+0.01</sup><sub>-0.96</sub></td><td>6.62<sup>+0.01</sup><sub>-2.74</sub></td><td>0.94<sup>+0.68</sup><sub>-0.01</sub></td><td>0.33<sup>+0.24</sup><sub>-0.01</sub></td></tr><tr><td>Spokes Cluster</td><td>2.69<sup>+0.09</sup><sub>-0.08</sub></td><td>425</td><td>3000</td><td>0.39<sup>+0.03</sup><sub>-0.03</sub></td><td>10.6<sup>+2.0</sup><sub>-2.7</sub></td><td>1.15<sup>+0.45</sup><sub>-0.01</sub></td><td>3.27<sup>+1.26</sup><sub>-0.03</sub></td><td>2.34<sup>+0.03</sup><sub>-0.66</sub></td><td>0.82<sup>+0.01</sup><sub>-0.23</sub></td></tr><tr><td>ASSC 50</td><td>1.79<sup>+0.09</sup><sub>-0.09</sub></td><td>200</td><td>—</td><td>4.24<sup>+0.37</sup><sub>-0.28</sub></td><td>6.0<sup>+3.6</sup><sub>-0.1</sub></td><td>0.32<sup>+0.01</sup><sub>-0.06</sub></td><td>—</td><td>5.62<sup>+1.35</sup><sub>-0.26</sub></td><td>—</td></tr><tr><td>Col. 197</td><td>1.51<sup>+0.11</sup><sub>-0.09</sub></td><td>79</td><td>—</td><td>2.64<sup>+0.55</sup><sub>-0.44</sub></td><td>9.6<sup>+0.1</sup><sub>-0.1</sub></td><td>0.20<sup>+0.01</sup><sub>-0.02</sub></td><td>—</td><td>7.52<sup>+1.23</sup><sub>-0.44</sub></td><td>—</td></tr><tr><td>NGC 6530</td><td>4.38<sup>+0.18</sup><sub>-0.17</sub></td><td>3125</td><td>40000</td><td>1.49<sup>+0.17</sup><sub>-0.19</sub></td><td>6.0<sup>+3.5</sup><sub>-0.1</sub></td><td>2.13<sup>+0.06</sup><sub>-0.58</sub></td><td>7.90<sup>+0.21</sup><sub>-2.15</sub></td><td>2.06<sup>+0.78</sup><sub>-0.07</sub></td><td>0.55<sup>+0.21</sup><sub>-0.02</sub></td></tr><tr><td>NGC 2244</td><td>3.22<sup>+0.09</sup><sub>-0.09</sub></td><td>1300</td><td>—</td><td>1.19<sup>+0.1</sup><sub>-0.24</sub></td><td>6.0<sup>+4.7</sup><sub>-0.1</sub></td><td>1.54<sup>+0.07</sup><sub>-0.40</sub></td><td>—</td><td>2.10<sup>+0.75</sup><sub>-0.11</sub></td><td>—</td></tr><tr><td>NGC 2237</td><td>1.67<sup>+0.21</sup><sub>-0.23</sub></td><td>250</td><td>—</td><td>2.07<sup>+0.54</sup><sub>-0.52</sub></td><td>6.0<sup>+1.2</sup><sub>-3.5</sub></td><td>0.51<sup>+0.09</sup><sub>-0.14</sub></td><td>—</td><td>3.29<sup>+1.42</sup><sub>-0.65</sub></td><td>—</td></tr><tr><td>Trumpler 14</td><td>7.36<sup>+0.19</sup><sub>-0.15</sub></td><td>5400</td><td>—</td><td>0.53<sup>+0.05</sup><sub>-0.07</sub></td><td>6.0<sup>+0.3</sup><sub>-4.7</sub></td><td>4.70<sup>+0.19</sup><sub>-1.30</sub></td><td>—</td><td>1.57<sup>+0.61</sup><sub>-0.06</sub></td><td>—</td></tr><tr><td>Trumpler 16</td><td>4.41<sup>+0.20</sup><sub>-0.21</sub></td><td>3250</td><td>—</td><td>1.72<sup>+0.12</sup><sub>-0.11</sub></td><td>9.5<sup>+0.2</sup><sub>-0.1</sub></td><td>1.60<sup>+0.08</sup><sub>-0.06</sub></td><td>—</td><td>2.75<sup>+0.17</sup><sub>-0.18</sub></td><td>—</td></tr></tbody></table>
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.07}^{+0.05}$ 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, $\alpha = \sigma_{3D}/\sigma_{3D,vir}$, for which $\alpha < 1$ indicates a sub-virial system, $\alpha = 1$ indicates a system in virial equilibrium, and $\alpha > 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 ($\alpha > 1$), with only $\lambda$ Ori and the S Mon cluster in NGC 2264 being in virial equilibrium or sub-virial ($\alpha < 1$). A small number of systems are close to being in virial equilibrium with $\alpha < 2$, such as Cha I South, Gamma Vel, 25 Ori, and Trumpler 14. For these clusters if their stellar mass has been underestimated or their radii or $\eta$ values overestimated, then this may bring them into virial equilibrium.<fnref n="7" /> 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 underestimated or $r_{eff}$ or $\eta$ 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, $\lambda$ 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 overestimates in this context. None the less, 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.

## 5 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 ASSC 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, Dekker & de Zeeuw (1997), before performing this analysis.

### 5.1 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 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 ($\alpha$, $\delta$, and along the line of sight, $\tau$, if the group is resolved). The gradient, $A$, and zero point, $B$, were fitted by maximizing 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, Bovy & Lang 2010) and which was marginalized 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.29}^{+0.33}$ (indicating contraction) up to $0.90_{-0.33}^{+0.45}$ 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.2$ km s$^{-1}$ pc$^{-1}$ (Wright & Mamajek 2018; Quintana & Wright 2021; Armstrong et al. 2022). We explored the effects of removing $2\sigma$ or $3\sigma$ 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.

Fig. 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.13}^{+0.11}$ mas yr$^{-1}$ deg$^{-1} = 0.20_{-0.04}^{+0.03}$ km s$^{-1}$ pc$^{-1}$ in Dec. These ~1.5$\sigma$ and ~5.5 $\sigma$ 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).

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$\sigma$ significance (Cha I North, Vela OB2, and 25 Ori).

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 high-significance indications of expansion are B35 (3$\sigma$ evidence for expansion in both RA and Dec.), NGC 6530 (6$\sigma$ 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$\sigma$ evidence for expansion in Dec., but only 1$\sigma$ in RA), and ASSC 50 (9$\sigma$ evidence for expansion in both RA and Dec.). Notably, we also find that Trumpler 14 is contracting at 3$\sigma$ significance in Dec. and Rho Oph is contracting at 2$\sigma$ 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$\sigma$) between axes

<footnote n="7">⁷Note that these calculations do not take into account the uncertainty on the stellar mass, which is difficult to quantify and could be as high as 20–30 per cent for some of these systems.</footnote>
<table><caption>**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.</caption><thead><tr><th rowspan="2">Group</th><th colspan="3">1D expansion gradients</th><th colspan="3">Rotated expansion gradients</th><th rowspan="2">ṽ<sub>out</sub><br>[km s<sup>-1</sup>]</th></tr><tr><th>ω<br>[km s<sup>-1</sup> pc<sup>-1</sup>]</th><th>α<br>[km s<sup>-1</sup> pc<sup>-1</sup>]</th><th>δ<br>[km s<sup>-1</sup> pc<sup>-1</sup>]</th><th>x<sub>1</sub><br>[km s<sup>-1</sup> pc<sup>-1</sup>]</th><th>x<sub>2</sub><br>[km s<sup>-1</sup> pc<sup>-1</sup>]</th><th>θ<br>[°]</th></tr></thead><tbody><tr><td>Rho Oph</td><td>-0.061<sup>+0.138</sup><sub>-0.309</sub></td><td>-0.14<sup>+0.21</sup><sub>-0.20</sub></td><td>-0.54<sup>+0.24</sup><sub>-0.22</sub></td><td>-0.18<sup>+0.22</sup><sub>-0.20</sub></td><td>-0.53<sup>+0.24</sup><sub>-0.21</sub></td><td>50</td><td>-0.34<sup>+0.01</sup><sub>-0.12</sub></td></tr><tr><td>Cha I (south)</td><td>-0.016<sup>+0.054</sup><sub>-0.234</sub></td><td>0.28<sup>+0.20</sup><sub>-0.26</sub></td><td>0.07<sup>+0.14</sup><sub>-0.12</sub></td><td>0.46<sup>+0.18</sup><sub>-0.20</sub></td><td>0.19<sup>+0.18</sup><sub>-0.18</sub></td><td>170</td><td>0.10<sup>+0.13</sup><sub>-0.08</sub></td></tr><tr><td>Cha I (north)</td><td>0.21<sup>+0.17</sup><sub>-0.19</sub></td><td>0.28<sup>+0.12</sup><sub>-0.14</sub></td><td>0.40<sup>+0.24</sup><sub>-0.15</sub></td><td>0.56<sup>+0.12</sup><sub>-0.15</sub></td><td>0.29<sup>+0.11</sup><sub>-0.11</sub></td><td>90</td><td>0.29<sup>+0.16</sup><sub>-0.10</sub></td></tr><tr><td>Gamma Vel</td><td>-</td><td>0.00<sup>+0.04</sup><sub>-0.04</sub></td><td>0.01<sup>+0.04</sup><sub>-0.04</sub></td><td>0.05<sup>+0.04</sup><sub>-0.04</sub></td><td>0.01<sup>+0.03</sup><sub>-0.03</sub></td><td>15</td><td>-0.06<sup>+0.06</sup><sub>-0.01</sub></td></tr><tr><td>Vela OB2</td><td>0.046<sup>+0.020</sup><sub>-0.048</sub></td><td>-0.14<sup>+0.09</sup><sub>-0.09</sub></td><td>0.43<sup>+0.10</sup><sub>-0.12</sub></td><td>0.43<sup>+0.11</sup><sub>-0.12</sub></td><td>-0.03<sup>+0.06</sup><sub>-0.07</sub></td><td>75</td><td>0.48<sup>+0.01</sup><sub>-0.08</sub></td></tr><tr><td>25 Ori</td><td>0.006<sup>+0.006</sup><sub>-0.008</sub></td><td>0.11<sup>+0.01</sup><sub>-0.03</sub></td><td>0.01<sup>+0.01</sup><sub>-0.01</sub></td><td>0.12<sup>+0.03</sup><sub>-0.03</sub></td><td>0.07<sup>+0.04</sup><sub>-0.04</sub></td><td>50</td><td>0.16<sup>+0.03</sup><sub>-0.01</sub></td></tr><tr><td>Barnard 30</td><td>-</td><td>0.38<sup>+0.13</sup><sub>-0.14</sub></td><td>-0.13<sup>+0.10</sup><sub>-0.13</sub></td><td>0.39<sup>+0.11</sup><sub>-0.12</sub></td><td>-0.12<sup>+0.09</sup><sub>-0.09</sub></td><td>0</td><td>0.13<sup>+0.16</sup><sub>-0.11</sub></td></tr><tr><td>λ Ori</td><td>-</td><td>0.05<sup>+0.04</sup><sub>-0.04</sub></td><td>0.15<sup>+0.03</sup><sub>-0.04</sub></td><td>0.20<sup>+0.03</sup><sub>-0.04</sub></td><td>-0.01<sup>+0.03</sup><sub>-0.02</sub></td><td>75</td><td>0.24<sup>+0.06</sup><sub>-0.01</sub></td></tr><tr><td>Barnard 35</td><td>-</td><td>0.50<sup>+0.11</sup><sub>-0.14</sub></td><td>0.43<sup>+0.10</sup><sub>-0.12</sub></td><td>0.79<sup>+0.08</sup><sub>-0.14</sub></td><td>0.57<sup>+0.11</sup><sub>-0.13</sub></td><td>120</td><td>0.61<sup>+0.02</sup><sub>-0.05</sub></td></tr><tr><td>S Mon Cluster</td><td>-</td><td>0.07<sup>+0.04</sup><sub>-0.04</sub></td><td>0.20<sup>+0.03</sup><sub>-0.04</sub></td><td>0.35<sup>+0.04</sup><sub>-0.04</sub></td><td>0.11<sup>+0.03</sup><sub>-0.04</sub></td><td>120</td><td>0.40<sup>+0.01</sup><sub>-0.06</sub></td></tr><tr><td>Spokes Cluster</td><td>-</td><td>-0.02<sup>+0.06</sup><sub>-0.06</sub></td><td>0.19<sup>+0.04</sup><sub>-0.05</sub></td><td>0.24<sup>+0.04</sup><sub>-0.05</sub></td><td>-0.02<sup>+0.06</sup><sub>-0.06</sub></td><td>85</td><td>0.23<sup>+0.03</sup><sub>-0.07</sub></td></tr><tr><td>ASSC 50</td><td>-</td><td>0.21<sup>+0.03</sup><sub>-0.03</sub></td><td>0.23<sup>+0.02</sup><sub>-0.03</sub></td><td>0.27<sup>+0.02</sup><sub>-0.02</sub></td><td>0.15<sup>+0.03</sup><sub>-0.03</sub></td><td>65</td><td>1.16<sup>+0.08</sup><sub>-0.06</sub></td></tr><tr><td>Col. 197</td><td>-</td><td>0.27<sup>+0.10</sup><sub>-0.12</sub></td><td>0.60<sup>+0.12</sup><sub>-0.13</sub></td><td>0.59<sup>+0.12</sup><sub>-0.13</sub></td><td>0.29<sup>+0.10</sup><sub>-0.11</sub></td><td>85</td><td>0.59<sup>+0.08</sup><sub>-0.12</sub></td></tr><tr><td>NGC 6530</td><td>-</td><td>-0.05<sup>+0.06</sup><sub>-0.06</sub></td><td>0.40<sup>+0.06</sup><sub>-0.07</sub></td><td>0.49<sup>+0.06</sup><sub>-0.06</sub></td><td>-0.03<sup>+0.04</sup><sub>-0.04</sub></td><td>80</td><td>0.76<sup>+0.10</sup><sub>-0.02</sub></td></tr><tr><td>NGC 2244</td><td>-</td><td>0.37<sup>+0.11</sup><sub>-0.13</sub></td><td>-0.01<sup>+0.06</sup><sub>-0.07</sub></td><td>0.62<sup>+0.14</sup><sub>-0.14</sub></td><td>0.16<sup>+0.11</sup><sub>-0.12</sub></td><td>140</td><td>0.88<sup>+0.07</sup><sub>-0.23</sub></td></tr><tr><td>NGC 2237</td><td>-</td><td>0.05<sup>+0.18</sup><sub>-0.18</sub></td><td>0.90<sup>+0.45</sup><sub>-0.33</sub></td><td>0.93<sup>+0.46</sup><sub>-0.35</sub></td><td>0.05<sup>+0.17</sup><sub>-0.17</sub></td><td>90</td><td>0.16<sup>+0.39</sup><sub>-0.13</sub></td></tr><tr><td>Trumpler 14</td><td>-</td><td>-0.13<sup>+0.18</sup><sub>-0.17</sub></td><td>-0.99<sup>+0.33</sup><sub>-0.29</sub></td><td>0.12<sup>+0.10</sup><sub>-0.10</sub></td><td>-0.13<sup>+0.21</sup><sub>-0.21</sub></td><td>85</td><td>0.07<sup>+0.01</sup><sub>-0.45</sub></td></tr><tr><td>Trumpler 16</td><td>-</td><td>0.20<sup>+0.19</sup><sub>-0.19</sub></td><td>0.27<sup>+0.14</sup><sub>-0.15</sub></td><td>0.49<sup>+0.19</sup><sub>-0.21</sub></td><td>0.23<sup>+0.13</sup><sub>-0.14</sub></td><td>30</td><td>1.37<sup>+0.09</sup><sub>-0.28</sub></td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**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 Dec., 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 is listed in each panel.</figcaption></figure>

(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).

## 5.2 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<sup>+0.24</sup><sub>-0.21</sub> up to 0.93<sup>+0.46</sup><sub>-0.35</sub> km s<sup>-1</sup> pc<sup>-1</sup>, for Rho Oph and NGC 2237, respectively. 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.

Fig. 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<sup>+0.13</sup><sub>-0.14</sub> mas yr<sup>-1</sup> deg<sup>-1</sup> = 0.35<sup>+0.04</sup><sub>-0.04</sub> km s<sup>-1</sup> pc<sup>-1</sup> along the axis
<figure><img src="image_7.png" /><figcaption>**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° (top 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 is listed in each panel.</figcaption></figure>

of strongest expansion and $0.39_{-0.13}^{+0.11}$ mas yr$^{-1}$ deg$^{-1}$ = $0.11_{-0.04}^{+0.03}$ 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 per cent) 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\sigma$ 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\sigma$ level). Notable examples are $\lambda$ Ori, Barnard 35, and the S Mon cluster (all at $6\sigma$), NGC 6530 (8$\sigma$), and ASSC 50 (13$\sigma$). 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. Ten out of 17 groups show evidence for significantly (> 2$\sigma$) anisotropic expansion, while the remaining 7 are consistent with either isotropic expansion or mildly significant expansion (at the $1\sigma$ level).

## 5.3 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, $v_{\text{out}}$, 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.28}^{+0.09}$ 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\sigma$.

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$\sigma$) and NGC 6530 (36$\sigma$), while notably different results are obtained, for example, for 25 Ori, which is found to be expanding with 16$\sigma$ significance by the median outward velocity method, but only 4$\sigma$ 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$\sigma$, providing further evidence of interesting kinematic behaviour in this young group.

## 5.4 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 per cent) showing evidence for expansion when the axis of expansion is allowed to rotate (Section 5.2) and 15 out of 18 groups (83 per cent) 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$\sigma$ 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 per cent), 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\sigma$, 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_{-6}^{+4}$ per cent of systems are expanding, while using the median outward velocity method that $99_{-1}^{+1}$ per cent of systems are expanding.

# 6 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
<figure><img src="image_8.png" /><figcaption>**Figure 7.** 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 $\tau = 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 $\sigma_{3D} = 1.18^{+0.24}_{-0.32} + 1.16^{+0.61}_{-0.50} \times (M_{\text{cluster}}/1000 M_{\odot})^{0.89^{+0.41}_{-0.30}} \text{ km s}^{-1}$, with grey lines showing 100 randomly sampled fits from the posterior distribution of fitted gradients.</figcaption></figure>

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 per cent and cluster masses by 20–30 per cent, which we have included in the fits performed here.<fnref n="8" /> We also calculate dynamical and relaxation time-scales 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 time-scale.

Fig. 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 $\sigma_{3D} = A + B \times M_{\text{cluster}}^C$ between these two quantities using Bayesian inference and an MCMC ensemble sampler to derive a fit of

$$ \sigma_{3D} = 1.18^{+0.24}_{-0.32} + 1.16^{+0.61}_{-0.50} \left( \frac{M_{\text{cluster}}}{1000 M_{\odot}} \right)^{0.89^{+0.41}_{-0.30}} \text{ km s}^{-1}, \quad (3) $$

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Relationship between the 3D velocity dispersion of all systems in our sample and their age, with uncertainties shown for the former (uncertainties for the latter are likely to be highly correlated so are not shown). Kendall's rank correlation test reveals a correlation of $\tau = -0.385$ and a $p$-value of 0.0293 indicating a strong inverse correlation.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Relationship between the virial ratio, $\alpha_{vir} = \sigma_{3D}/\sigma_{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 $\tau = 0.516$ and a $p$-value of 0.0022 indicating a strong positive correlation.</figcaption></figure>

which is shown in Fig. 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 $\le 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$ (Fig. 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 overdensities 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 Fig. 9. The outlier in Fig. 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 $\alpha_{vir} > 1$) will expand to larger radii, with the more unbound systems expanding

<footnote n="8">⁸These uncertainties are likely to be correlated between clusters as they depend on stellar evolutionary models (for ages) or the initial mass function and binary properties (for masses), each of which have standard values that are commonly used between studies, and therefore the impact of such uncertainties is greatly reduced for the rank correlation tests used in this section.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Relationship between the virial ratio, $\alpha_{vir} = \sigma_{3D}/\sigma_{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 $\tau = -0.302$ and a *p*-value of 0.081 indicating a weak inverse correlation.</figcaption></figure>

faster. This should lead to a correlation between $\alpha_{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 $\alpha_{vir}$ is linearly dependent on the radius, and therefore if all other parameters that $\alpha_{vir}$ depends on ($\sigma_{3D}$ and *M* specifically) either remained constant or had dependencies that cancelled each other out, one would expect to observe a linear correlation between $\alpha_{vir}$ and radius. Since this does not appear to be the case (see e.g. Fig. 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 (Fig. 10) with a *p*-value 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 Fig. 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 $v_{out}$, 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 time-scale (*p*-value $= 7.0 \times 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 time-scale becomes.

Finally, we also observe a strong positive correlation (*p*-value $= 0.027$) between radius and age, with older systems being larger.

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Relationship between the virial ratio, $\alpha_{vir} = \sigma_{3D}/\sigma_{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 $\tau = 0.397$ and a *p*-value of 0.027 indicating a strong correlation.</figcaption></figure>

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 time-scale (*p*-value $= 0.014$), which is likely to be a product of this since the dynamical time-scale 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 associated with molecular gas (as the parental molecular cloud is either consumed or dispersed).

# 7 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 time-scales. We have compared these properties with each other, as well as with literature ages and group masses, and dynamical and relaxation time-scales, 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.

## 7.1 Dynamical and virial state of young groups

The young stellar groups studied in this work typically have 3D velocity dispersions of $1–2 \text{ km s}^{-1}$, ranging from $0.61_{-0.02}^{+0.03} \text{ km s}^{-1}$ for Gamma Vel to $7.36_{-0.15}^{+0.19} \text{ 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 overestimate the velocity dispersion by 10 per cent 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 suffi-
ciently 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 per cent. Given this, the ratio of the velocity dispersion to the virial velocity dispersion would be expected to scale as

$$ \alpha = \frac{\sigma_{3D}}{\sigma_{3D, vir}} \propto \frac{M^{0.94}}{(M/R)^{0.5}} \propto (MR)^{0.5}. \quad (4) $$

The mass dependence of cluster radii has been studied by various authors with dependencies varying from approximately $R \propto M^{0.25}$ (Brown & Gnedin 2021; Dobbs et al. 2022) to $R \propto M^{0.5}$ (Adams et al. 2006; Pfalzner 2011). This leads to a mass dependency of $\alpha \propto 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.

## 7.2 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 per cent using the median outward velocity, 89 per cent using 1D expansion gradients, and 94 per cent using 2D rotated expansion gradient fits. The significance of these results vary, but 12 of the 17 groups with evidence for expansion (71 per cent) do so at the $3\sigma$ 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 per cent of their groups had positive median outward velocities (compared to 83 per cent here), but that only 57 per cent of their groups had positive median outward velocities to a confidence of > 1σ (compared to 78 per cent 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, 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 ($\alpha_{vir} = 0.94^{+0.68}_{-0.01}$) and λ Ori ($\alpha_{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).

## 7.3 Expansion time-scales and ages

Expansion time-scales can be calculated from the expansion gradients by simple inversion (with a correction factor to account for units). We calculate expansion time-scales for all our groups using the largest expansion rate (from the rotated expansion gradient fit) as this will give the smallest expansion time-scale for the two expansion axes. We find that the expansion time-scales 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 time-scales 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 time-scale 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 time-scales that are significantly smaller than their evolutionary ages; Vela OB2 (2.3 versus 14 Myr), λ Ori (4.9 versus 10 Myr), and Col. 197 (1.6 versus 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 time-scale 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 time-scale 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 time-scales that are significantly larger than their evolutionary ages; Trumpler 14 (8 versus 2.5 Myr) and the Spokes cluster within NGC 2264 (4.1 versus 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.

## 7.4 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 $\alpha_{vir} = 2.93^{+0.38}_{-0.13}$, suggesting that the system is gravitationally unbound. However, when one takes into account the approximately 1750 M$_{\odot}$ (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 per cent 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 on to 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.

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, González-Samaniego & Colín 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

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Relationship between the dynamical time-scale, $t_{dyn} = (GM/r_{vir}^3)^{-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 time-scale of ~120 Myr) and Barnard 35 (age 2.6 Myr and dynamical time-scale ~45 Myr) are not shown. A 1:1 relationship is shown as a grey dashed line.</figcaption></figure>

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.

## 7.5 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 time-scale set by their dynamical time-scale. 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 time-scale in Fig. 12. This figure shows that the majority of groups (12 out of 18) are older than their dynamical time-scale. The groups that are younger than their dynamical time-scale (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 per cent.

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 time-scale) 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.

## 7.6 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 per cent of young groups are in the process of expanding is consistent with the observation by Lada & Lada (2003) that 90 per cent 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, virialized 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 per cent (only two groups are consistent with being in virial equilibrium). The fraction of groups that are inconsistent with having isotropic velocity distributions is >90 per cent (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 under-estimated or their radii overestimated 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.

## 8 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⁻¹ (1D velocity dispersions range from 0.36 to 4.3 km s⁻¹). 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-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 virial 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 expansion. 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 on to a central cluster.

We conclude that, since the majority of clusters are not dynamically evolved, despite being older than their dynamical time-scales, that these clusters did not form as we observe them now, but originally were larger and had lower densities (and thus had longer dynamical time-scales). 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.

## ACKNOWLEDGEMENTS

We are grateful to Núria Miret-Roig, Alexis Quintana, and Steve Stahler, for reading this manuscript and sending valuable comments. NJW acknowledges an STFC Ernest Rutherford Fellowship (grant number ST/M005569/1) and a Leverhulme Trust Research Project Grant (RPG-2019-379). This research has made use of NASA's Astrophysics Data System and the Simbad and VizieR databases, operated at CDS, Strasbourg.

This work is based on data products from observations made with ESO Telescopes on the La Silla Paranal Observatory under program ID 188.B-3002. These data products have been processed by the Cambridge Astronomy Survey Unit (CASU) at the Institute of Astronomy, University of Cambridge, and by the FLAMES-UVES reduction team at INAF - Osservatorio Astrofisico di Arcetri. These data have been obtained from the Gaia-ESO Survey Data Archive, prepared and hosted by the Wide Field Astronomy 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 MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia.

## DATA AVAILABILITY

The GES and Gaia data used in this article are available from the ESO and Gaia archives, respectively. The sample used in this article and the physical quantities calculated and derived in this article will be uploaded to Vizier.

## REFERENCES

* Adams F. C., 2010, ARA&A, 48, 47

* Adams F. C., Proszkow E. M., Fatuzzo M., Myers P. C., 2006, ApJ, 641, 504

* Arenou F. et al., 2018, A&A, 616, A17

* Armstrong J. J., Wright N. J., Jeffries R. D., 2018, MNRAS, 480, L121

* Armstrong J. J., Wright N. J., Jeffries R. D., Jackson R. J., 2020, MNRAS, 494, 4794

* Armstrong J. J., Wright N. J., Jeffries R. D., Jackson R. J., Cantat-Gaudin T., 2022, MNRAS, 517, 5704

* Arnold B., Wright N. J., 2024, MNRAS, 531, 1191

* Arnold B., Wright N. J., Parker R. J., 2022, MNRAS, 515, 2266

* Ascenso J., Alves J., Beletsky Y., Lago M. T. V. T., 2007, A&A, 466, 137

* Banerjee S., Kroupa P., 2015, MNRAS, 447, 728

* Bell C. P. M., Naylor T., Mayne N. J., Jeffries R. D., Littlefair S. P., 2013, MNRAS, 434, 806

* Blaauw A., 1964, ARA&A, 2, 213

* Blitz L., Thaddeus P., 1980, ApJ, 241, 676

* Bonatto C., Bica E., 2010, A&A, 516, A81

* Bonito R., Prisinzano L., Guarcello M. G., Micela G., 2013, A&A, 556, A108

* Bonnell I. A., Bate M. R., Vine S. G., 2003, MNRAS, 343, 413

* Bragaglia A. et al., 2022, A&A, 659, A200

* Bravi L. et al., 2018, A&A, 615, A37

* Brown A. G. A., Dekker G., de Zeeuw P. T., 1997, MNRAS, 285, 479

* Brown G., Gnedin O. Y., 2021, MNRAS, 508, 5935

* Cantat-Gaudin T. et al., 2019, A&A, 626, A17

* Cantat-Gaudin T. et al., 2020, A&A, 640, A1

* Chevance M. et al., 2022, MNRAS, 509, 272

* Cottaar M., Meyer M. R., Andersen M., Espinoza P., 2012, A&A, 539, A5

* Dahm S. E., 2008, in Reipurth B.ed., Handbook of Star Forming Regions, Vol. I., University of Hawaii, 390 Ashton Avenue San Francisco, CA, USA, p. 966

* Dale J. E., Ercolano B., Bonnell I. A., 2015, MNRAS, 451, 987

* Damiani F. et al., 2017, A&A, 604, A135

* Damiani F. et al., 2014, A&A, 566, A50

* Daniëlski C., Babusiaux C., Ruiz-Dern L., Sartoretti P., Arenou F., 2018, A&A, 614, A19

* Dias W. S., Monteiro H., Moitinho A., Lépine J. R. D., Carraro G., Paunzen E., Alessi B., Villela L., 2021, MNRAS, 504, 356

* Dobbs C. L., Bending T. J. R., Pettitt A. R., Buckner A. S. M., Bate M. R., 2022, MNRAS, 517, 675

* Dolan C. J., Mathieu R. D., 2002, AJ, 123, 387

* Duchêne G., Kraus A., 2013, ARA&A, 51, 269

* Dzib S. A. et al., 2017, ApJ, 834, 139

* Ekström S. et al., 2012, A&A, 537, A146

* Elmegreen B. G., 2002, ApJ, 564, 773

* Elson R. A. W., Fall S. M., Freeman K. C., 1987, ApJ, 323, 54 (EFF)

* Erickson K. L., Wilking B. A., Meyer M. R., Robinson J. G., Stephenson L. N., 2011, AJ, 142, 140

* Fabricius C. et al., 2021, A&A, 649, A5

* Figer D. F., 2008, in Bresolin F., Crowther P. A., Puls J., eds., Massive Stars as Cosmic Engines. Cambridge University Press, p. 247

* Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

* Franciosini E., Sacco G. G., Jeffries R. D., Damiani F., Roccatagliata V., Fedele D., Randich S., 2018, A&A, 616, L12

* Franciosini E. et al., 2022a, A&A, 668, A49

* Franciosini E. et al., 2022b, A&A, 659, A85

* Gaia Collaboration, 2021, A&A, 649, A1

* Gaia Collaboration, 2016, A&A, 595, A1

* Galli P. A. B. et al., 2021, A&A, 646, A46

* Getman K. V., Kuhn M. A., Feigelson E. D., Broos P. S., Bate M. R., Garmire G. P., 2018, MNRAS, 477, 298

* Gilmore G. et al., 2022, A&A, 666, A120

* Goodman J., Weare J., 2010, Comm. App. Math. Comp. Sci., 5, 65

* Grasser N. et al., 2021, A&A, 652, A2

* Gutermuth R. A. et al., 2008, ApJ, 674, 336

* Hills J. G., 1980, ApJ, 235, 986

* Hogg D. W., Bovy J., Lang D., 2010, preprint (arXiv:1008.4686)

* Hourihane A. et al., 2023, A&A, 676, A129

* Hur H., Sung H., Bessell M. S., 2012, AJ, 143, 41

* Ivezić Ž., Connelly A. J., VanderPlas J. T., Gray A., 2014, Statistics, Data Mining, and Machine Learning in Astronomy. Princeton Univ. Press, Princeton

* Jackson R. J. et al., 2015, A&A, 580, A75

* Jackson R. J. et al., 2020, MNRAS, 496, 4701

* Jeffries R. D., Naylor T., Walter F. M., Pozzo M. P., Devey C. R., 2009, MNRAS, 393, 538

* Jeffries R. D. et al., 2014, A&A, 563, A94

* Jeffries R. D. et al., 2017, MNRAS, 464, 1456

* Kharchenko N. V., Piskunov A. E., Röser S., Schilbach E., Scholz R. D., 2005, A&A, 438, 1163

* Kharchenko N. V., Piskunov A. E., Schilbach E., Röser S., Scholz R. D., 2013, A&A, 558, A53

* Kounkel M. et al., 2018, AJ, 156, 84

* Kruijssen J. M. D., 2012, MNRAS, 426, 3008

* Kruijssen J. M. D., Maschberger T., Moeckel N., Clarke C. J., Bastian N., Bonnell I. A., 2012, MNRAS, 419, 841

* Krumholz M. R., McKee C. F., 2020, MNRAS, 494, 624

* Kuhn M. A. et al., 2014, ApJ, 787, 107

* Kuhn M. A., Feigelson E. D., Getman K. V., Sills A., Bate M. R., Borissova J., 2015, ApJ, 812, 131

* Kuhn M. A., Hillenbrand L. A., Sills A., Feigelson E. D., Getman K. V., 2019, ApJ, 870, 32
Lada C. J., Lada E. A., 2003, ARA&A, 41, 57

Lada C. J., Margulis M., Dearborn D., 1984, ApJ, 285, 141

Lanzafame A. C. et al., 2015, A&A, 576, A80

Lindegren L. et al., 2021, A&A, 649, A2

Loren R. B., 1989, ApJ, 338, 902

Luhman K. L., 2007, ApJS, 173, 104

Maíz Apellániz J., 2019, A&A, 630, A119

Maschberger T., 2013, MNRAS, 429, 1725

Mizuno A., Yamaguchi R., Tachihara K., Toyoda S., Aoyama H., Yamamoto H., Onishi T., Fukui Y., 2001, PASJ, 53, 1071

Mužić K., Scholz A., Peña Ramírez K., Jayawardhana R., Schödel R., Geers V. C., Cieza L. A., Bayo A., 2019, ApJ, 881, 79

Odenkirchen M., Grebel E. K., Dehnen W., Rix H.-W., Cudworth K. M., 2002, AJ, 124, 1497

Oliver R. J., Masheder M. R. W., Thaddeus P., 1996, A&A, 315, 578

Palla F., Randich S., Flaccomio E., Pallavicini R., 2005, ApJ, 626, L49

Pancino E. et al., 2017, A&A, 598, A5

Parker R. J., Goodwin S. P., 2009, MNRAS, 397, 1041

Parker R. J., Wright N. J., 2016, MNRAS, 457, 3430

Parker R. J., Wright N. J., Goodwin S. P., Meyer M. R., 2014, MNRAS, 438, 620

Pasquini L. et al., 2002, The Messenger, 110, 1

Pecaut M. J., Mamajek E. E., 2016, MNRAS, 461, 794

Pettersson B., Reipurth B., 1994, A&AS, 104, 233

Pfalzner S., 2009, A&A, 498, L37

Pfalzner S., 2011, A&A, 536, A90

Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A, 48, 431

Prisinzano L., Damiani F., Guarcello M. G., Micela G., Sciortino S., Tognelli E., Venuti L., 2018, A&A, 617, A63

Prisinzano L. et al., 2019, A&A, 623, A159

Proszkow E.-M., Adams F. C., Hartmann L. W., Tobin J. J., 2009, ApJ, 697, 1020

Quintana A. L., Wright N. J., 2021, MNRAS, 508, 2370

Quintana A. L., Wright N. J., Jeffries R. D., 2023, MNRAS, 522, 3124

Raghavan D. et al., 2010, ApJS, 190, 1

Randich S., Aharpour N., Pallavicini R., Prosser C. F., Stauffer J. R., 1997, A&A, 323, 86

Randich S. et al., 2022, A&A, 666, A121

Reggiani M. M., Meyer M. R., 2011, ApJ, 738, 60

Reiter M., Parker R. J., 2019, MNRAS, 486, 4354

Rigliaco E. et al., 2016, A&A, 588, A123

Roccatagliata V., Sacco G. G., Franciosini E., Randich S., 2018, A&A, 617, L4

Sacco G. G. et al., 2014, A&A, 565, A113

Sacco G. G. et al., 2017, A&A, 601, A97

Sana H., Momany Y., Gieles M., Carraro G., Beletsky Y., Ivanov V. D., de Silva G., James G., 2010, A&A, 515, A26

Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525

Schwarz G., 1978, Ann. Stat., 6, 461

Sills A., Rieder S., Scora J., McCloskey J., Jaffa S., 2018, MNRAS, 477, 1903

Soderblom D. R., 2010, ARA&A, 48, 581

Spitzer L., 1987, Dynamical Evolution of Globular Clusters. Princeton Univ. Press, Princeton

Suárez G., Downes J. J., Román-Zúñiga C., Cerviño M., Briceño C., Petr-Gotzens M. G., Vivas K., 2019, MNRAS, 486, 1718

Sung H., Chun M.-Y., Bessell M. S., 2000, AJ, 120, 333

Sung H., Stauffer J. R., Bessell M. S., 2009, AJ, 138, 1116

Takeuchi T. et al., 2010, PASJ, 62, 557

Teixeira P. S. et al., 2006, ApJ, 636, L45

Tutukov A. V., 1978, A&A, 70, 57

Vázquez-Semadeni E., González-Samaniego A., Colín P., 2017, MNRAS, 467, 1313

Venuti L., Damiani F., Prisinzano L., 2019, A&A, 621, A14

Wang J., Townsley L. K., Feigelson E. D., Broos P. S., Getman K. V., Román-Zúñiga C. G., Lada E., 2008, ApJ, 675, 464

Wang J., Feigelson E. D., Townsley L. K., Broos P. S., Román-Zúñiga C. G., Lada E., Garmire G., 2010, ApJ, 716, 474

Ward J. L., Kruijssen J. M. D., Rix H.-W., 2020, MNRAS, 495, 663

Weidner C., Kroupa P., 2006, MNRAS, 365, 1333

Wilking B. A., Gagné M., Allen L. E., 2008, in Reipurth B.ed., Handbook of Star Forming Regions, Volume II. University of Hawaii, 390 Ashton Avenue San Francisco, CA, USA, p. 351

Wolk S. J. et al., 2011, ApJS, 194, 12

Wright N. J., 2020, New Astron. Rev., 90, 101549

Wright N. J., Mamajek E. E., 2018, MNRAS, 476, 381

Wright N. J., Parker R. J., 2019, MNRAS, 2232, 2694

Wright N. J., Bouy H., Drew J. E., Sarro L. M., Bertin E., Cuillandre J.-C., Barrado D., 2016, MNRAS, 460, 2593

Wright N. J. et al., 2019, MNRAS, 486, 2477

Wright N. J., Kounkel M., Zari E., Goodwin S., Jeffries R. D., 2023, in Inutsuka S., Aikawa Y., Muto T., Tomida K., Tamura M., eds., ASP Conf. Ser. Vol. 534, Protostars and Planets VII. Astron. Soc. Pac., San Francisco, p. 129

Zamora-Avilés M., Ballesteros-Paredes J., Hernández J., Román-Zúñiga C., Lora V., Kounkel M., 2019, MNRAS, 488, 3406

Zari E., Brown A. G. A., de Zeeuw P. T., 2019, A&A, 628, 123

## APPENDIX A: THE OBSERVED GROUPS

Here, we summarize 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 summarize how the targeted stars in the GES observations were separated into groups.

(i) **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 deg on the sky (14 × 7 pc), though the nine FLAMES pointings obtained by GES cover an area of only 1.4 × 1.4 deg (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_\odot$ (Erickson et al. 2011; Rigliaco et al. 2016), which rises to 106 $M_\odot$ to account for unresolved binaries<fnref n="9" /> with an additional ~1750 $M_\odot$ 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).

(ii) **Chamaeleon I (North and South)**, our most reddened targets with $A_V \sim 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 deg 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_\odot$ and accounting for unresolved binaries equates to ~113 $M_\odot$. Cha I North has ~10 per cent fewer stars than Cha I South (Galli et al. 2021) and so we estimate their total stellar masses as 54 and 59 $M_\odot$, respectively. The mass of the molecular cloud in which Cha I is embedded is estimated to be 1000 $M_\odot$ (Mizuno et al. 2001). We divide the data for Cha I into the North ($\delta > -77°$) and South ($\delta < -77°$) clusters (Roccatagliata et al. 2018). The GES data for Cha I were previously studied by Sacco et al. (2017).

<footnote n="9">⁹Assuming a 50 per cent binary fraction and a mean mass ratio of 0.5, a reasonable approximation across multiple primary star masses (e.g. Duchêne & Kraus 2013).</footnote>
(iii) **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_\odot$ for Gamma Vel (Jeffries et al. 2014) and 1285 ± 110 $M_\odot$ for Vela OB2 (Armstrong, Wright & Jeffries 2018). We follow Armstrong et al. (2020) and select stars having $(\mu_\alpha + 6.53)^2 + (\mu_\delta - 9.8)^2 \le 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).

(iv) **25 Ori** is a young, non-embedded group in the Orion OB1 association, which Franciosini et al. (2022b) recently determined an age of $19_{-7.0}^{+1.5}$ 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 \pm 25 M_\odot$ (Suárez et al. 2019), and when accounting for unresolved binaries this equates to $400 \pm 30 M_\odot$. We exclude stars outside of the central group as non-members (those with $\delta < 1.2°$). The GES data for 25 Ori has not previously been studied dynamically.

(v) **λ 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_\odot$. When accounting for unresolved binaries and revising for the nearer distance from Gaia, we estimate the total stellar mass to be $650 \pm 120 M_\odot$. The group radius is approximately 3° or ~20 pc. The GES data for λ Ori has not previously been studied dynamically.

(vi) **Barnard 30 and 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 \pm 1.3$ and $2.6 \pm 1.3$ Myr (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 per cent. Combining this with a mean stellar mass of $0.4 M_\odot$ and accounting for binaries gives estimated stellar masses of 60 and $73 M_\odot$, 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.

(vii) **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 an ~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, Stauffer & Bessell 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_\odot$, when accounting for unresolved binaries. Oliver, Masheder & Thaddeus (1996) estimated the molecular cloud that NGC 2264 is embedded within has a total mass of ~28000 $M_\odot$ (when scaled to our Gaia distance), but from molecular maps we estimate that only ~10 per cent of the molecular cloud’s mass (~3000 $M_\odot$) is centred around each cluster and therefore relevant for its dynamics. We select members of the Spokes cluster as those with $\delta < 9.72°$ and members of the S Mon cluster with $\delta > 9.72°$. The GES data for NGC 2264 has not previously been studied dynamically.

(viii) **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) in the RCW 33 H II region and has a Gaia distance of 912 ± 3 pc. The age of the group has been estimated from ~5 (Prisinzano et al. 2018) to ~11.5 Myr (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_\odot$ and if on the main sequence an age $\lesssim 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_\odot$ (Weidner & Kroupa 2006). Prisinzano et al. (2018) estimate a total mass of 50–86 $M_\odot$, which when scaled to account for binarity gives a mass of $85 \pm 23 M_\odot$. We compromise and use a mass of 200 ± 100 $M_\odot$. The group has a radius of ~0.25 deg, or 4 pc at a distance of 912 pc. The GES data for ASCC 50 have not previously been studied dynamically.

(ix) **Collinder 197** 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 by Bonatto & Bica (2010) to be $660_{-59}^{+102} M_\odot$, though Prisinzano et al. (2018) estimate a mass of 45–81 $M_\odot$, or $79 \pm 23 M_\odot$ 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_\odot$ 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.

(x) **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_\odot$ (Wright et al. 2019), while the molecular cloud it is associated with has an estimated mass of 40 000 $M_\odot$ (Takeuchi et al. 2010). The GES data for NGC 6530 was previously studied by Wright et al. (2019) and Wright & Parker (2019).

(xi) **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_\odot$ based on a mass of 40–50 $M_\odot$ 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_\odot$ (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.

(xii) **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_\odot$ (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 per cent closer than its Gaia EDR3 parallax implies). Its radius is ~0.2 deg or ~5 pc at 1.49 kpc. The GES data for NGC 2237 have not previously been studied dynamically.

(xiii) **Trumpler 14 and 16**, our most distant targets ($d \sim 2.2$ kpc), are both young (2.5 and 2 Myr, respectively; Hur et al. 2012), compact (radii ~ 1 pc; Ascenso et al. 2007), but not centrally concentrated (Reiter & Parker 2019). Trumpler 14 has an approximate total mass of $5400_{-1900}^{+4100} M_\odot$ (Sana et al. 2010, once unresolved binaries are accounted for). Trumpler 16 is comparable, but slightly less massive. Wolk et al. (2011) estimate a total stellar population of $6500 \pm 650$ stars, which suggests a total mass, including binaries, of $3250 \pm 325 M_\odot$. We separate the two clusters spatially, with Trumpler 14 members having $\delta > -59.76^\circ + 0.35 (\alpha - 160.8^\circ)$. The GES data for both Trumpler 14 and Trumpler 16 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 Figs 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 \alpha$ full width at zero intensity (FWZI) measurements greater than 4 Å (Bonito et al. 2013).

(ii) The star must not have a gravity-sensitive index $\gamma > 1$ and $T_{\text{eff}} < 5600$ K (which together indicate the star is a cool giant, Damiani et al. 2014).

(iii) The star must have a parallax within $2\sigma$ 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 $\le 1.4$ quality cut have their astrometry discarded but they themselves, and their RVs, are not discarded. This decision was made to maximize 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 Fig. C1 the proper motion 'vector point' diagrams for all groups and clusters, showing the members of each group with valid proper motions.
<figure><img src="image_14.png" /><figcaption>**Figure B1.** Quantities derived from GES spectroscopy used for membership selection and parallax distributions for all groups. Left-hand panels show the gravity index $\gamma$, centre-left panels show the EW of the lithium 6708 Å line, and centre-right-hand panels show the FWZI of $H\alpha$, 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 $\gamma-T_{\text{eff}}$ plot and to identify young stars above the thresholds in the EW(Li)-$T_{\text{eff}}$ and FWZI($H\alpha$)-$T_{\text{eff}}$ plots. A typical error bar is shown in the top-left corner of each panel illustrating the typical uncertainties of 100 K in $T_{\text{eff}}$, 13 mÅ for EW(Li), 0.011 for $\gamma$ and 1.1 Å for FWZI($H\alpha$). 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.</figcaption></figure>
NGC 2264, Stars = 1857, Members = 524

<figure><img src="image_15.png" /></figure>

ASSC 50, Stars = 1225, Members = 193

<figure><img src="image_16.png" /></figure>

Col. 197, Stars = 409, Members = 112

<figure><img src="image_17.png" /></figure>

NGC 6530, Stars = 1983, Members = 644

<figure><img src="image_18.png" /></figure>

NGC 2244, Stars = 452, Members = 161

<figure><img src="image_19.png" /></figure>

Trumpler 14, Stars = 1902, Members = 372

<figure><img src="image_20.png" /><figcaption>Figure B2. As per Fig. B1.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**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⁻¹, which in the vast majority of cases is smaller than the symbol size used.</figcaption></figure>

This paper has been typeset from a TeX/LA TeX file prepared by the author.