The Kinematics of the Scorpius-Centaurus OB Association from Gaia DR1
Nicholas J. Wright Eric E. Mamajek
astro-ph.SR astro-ph.GA

*
    Joshua N. Winn
    December 30, 2023
=====================





We present a kinematic study of the Scorpius-Centaurus (Sco-Cen) OB association (Sco OB2) using Gaia DR1 parallaxes and proper motions. Our goal is to test the classical theory that OB associations are the expanded remnants of dense and compact star clusters disrupted by processes such as residual gas expulsion. Gaia astrometry is available for 258 out of 433 members of the association, with revised Hipparcos astrometry used for the remainder. We use this data to confirm that the three subgroups of Sco-Cen are gravitationally unbound and have non-isotropic velocity dispersions, suggesting they have not had time to dynamically relax. We also explore the internal kinematics of the subgroups to search for evidence of expansion. We test Blaauw's classical linear model of expansion, search for velocity trends along the Galactic axes, compare the expanding and non-expanding convergence points, perform traceback analysis assuming both linear trajectories and using an epicycle approximation, and assess the evidence for expansion in proper motions corrected for virtual expansion / contraction. None of these methods provide coherent evidence for expansion of the subgroups, with no evidence to suggest that the subgroups had a more compact configuration in the past. We find evidence for kinematic substructure within the subgroups that supports the view that they were not formed by the disruption of individual star clusters. We conclude that Sco-Cen was likely born highly substructured, with multiple small-scale star formation events contributing to the overall OB association, and not as single, monolithic bursts of clustered star formation.





stars: formation - stars: kinematics and dynamics - open clusters and associations: individual: Scorpius-Centaurus, Sco OB2, Upper Scorpius, Upper Centaurus-Lupus, Lower Centaurus-Crux.












§ INTRODUCTION




The environment where stars form and spend the first few million years of their lives has profound consequences for the rest of their lives and their potential to host habitable exoplanets. Ultraviolet (UV) radiation from nearby massive stars can lead to photoevaporation of protoplanetary disks <cit.>, while close encounters in dense clusters can affect the properties of binary and multiple systems <cit.>, protoplanetary disks <cit.>, and influence the formation of planetary systems <cit.>. Stars born in low-density environments <cit.> where dynamical interactions are rare and UV radiation fields are weaker may form planetary and binary systems with little or no external disruption.

The majority of young stars are observed in groups or clusters of some sort, but by an age of 10 Myr only ∼10% of stars are found in bound clusters <cit.>. The classical explanation for this is that star clusters form embedded within molecular clouds, held together by the gravitational potential of both the stars and the gas, but when feedback disperses the gas left over from star formation <cit.> then the cluster becomes supervirial and will expand and disperse <cit.>. During the expansion the system would briefly be visible as a low-density group of young stars known as an association <cit.>, before dispersing into the Galactic field.

This theoretical framework has existed for many decades, and while many observations support aspects of this model <cit.> it has been difficult to test kinematically due to the lack of high-precision proper motions necessary to assess evidence for expansion. <cit.> performed the first such study, using ground-based proper motions to study the kinematics of the Cygnus OB2 association. They found that the association lacked the coherent radial expansion pattern expected if it was an expanding star cluster. Its non-isotropic velocity dispersions and both physical <cit.> and kinematic substructure instead argued for an origin as a highly substructured and extended association of stars.

The recent release of data from the Gaia satellite <cit.> is set to change this picture, providing orders of magnitude improvement in the quality and quantity of proper motions. In this study we use Gaia astrometry to study the kinematics of the high- and intermediate-mass stars in the Sco-Cen OB association to assess its dynamical state, search for evidence of expansion, and constrain its initial conditions.

Sco-Cen (also known as Sco OB2) is the nearest OB association to the Sun and therefore the nearest region that has recently formed massive stars. Spanning distances of ∼100–150 pc <cit.> and being at least 100 pc across, the association spans almost 90 degrees on the sky (see Figure <ref>). The association consists of three subgroups <cit.> known as Upper Scorpius (US), Upper Centaurus-Lupus (UCL) and Lower Centaurus-Crux (LCC), with median ages of 11, 16, and 17 Myr respectively, though each has a considerable spread of ages <cit.>. Because of its age the region has likely already seen several supernovae <cit.>, one of which may have been responsible for the ejected runaway star ζ Oph <cit.>. 

The association has been the focus of numerous studies of its low-mass population <cit.>, with studies focussing on the binary properties <cit.>, circumstellar disks <cit.>, stellar ages <cit.>, eclipsing binaries <cit.>, and low-mass stars and brown dwarfs <cit.>. Its high-mass population has been studied by numerous authors, including comprehensive studies by <cit.> and later, using Hipparcos astrometry, by <cit.>, <cit.>, and <cit.>.

This paper is outlined as follows. In Section 2 we present the data used. In Section 3 we recalculate the distance to the subgroups and quantify their 3D structure. In Section 4 we use proper motions and radial velocities (RVs) to measure the 3D velocity dispersions and bulk motions of the subgroups. In Section 5 we assess evidence for the expansion of the association using various methods, and in Section 6 we study its kinematic substructure. Finally in Section 7 we discuss our results and their implications.



§ OBSERVATIONAL DATA


The sample of Sco-Cen members used for our kinematic study is based on the Bayesian membership analysis of <cit.>, who used the updated Hipparcos catalogue <cit.> to revise the kinematic selection of <cit.>. Their linear model exploits positions, proper motions and parallaxes, as well as RVs from the 2nd version of the Catalogue of Radial Velocities with Astrometric Data <cit.>[<cit.> note that some of the RVs used by <cit.> are non-spectroscopic `astrometric' RVs that are inappropriate for calculating space motions or membership. We investigated the phase-space distribution of these sources and found that their proper motions and parallaxes are still consistent with being members of Sco-Cen and so have retained them.] to produce a list of 436 high- and intermediate-mass members of Sco-Cen (the majority of which are of A & B spectral types). Many of the sources discarded by <cit.> were rejected based on their parallax or RVs, two parameters that were not used by <cit.> in their original membership selection. 

We refined this sample based on recent studies that suggest some of these stars are unlikely to be members of Sco-Cen. <cit.> identify HIP 62428 as an F0III star whose position in the Hertzsprung-Russell diagram is inconsistent with being a member of Sco-Cen, while HIP 70833 and HIP 75824 exhibit low lithium equivalent widths that imply they are older than typical Sco-Cen stars (in the case of HIP 70833 the lithium is measured in its K-type binary companion). The vast majority of stars identified as non-members by <cit.>, <cit.> and <cit.> that were included in the original sample of <cit.> were not selected by <cit.> as members, reinforcing the validity of their membership selection. <cit.> recommend excluding a number of stars based on their kinematics, one of which is retained in the <cit.> catalogue, HIP 56227. We choose to retain this star for the time being and minimise the impact kinematic selection may have on our sample.

This leaves us with a membership list of 433 stars (after removing 3 stars), 107 in US, 179 in UCL, and 147 in LCC, as shown in Figure <ref>. For the handful of stars that cannot be assigned to a subgroup based on the division of <cit.> we put them in the nearest subgroup on the plane of the sky. While <cit.> concluded that the strict boundaries between the three subgroups are somewhat arbitrary we continue to use them here for the ease of analysing such an extended OB association and to allow comparison with previous studies. Due to the current lack of high-precision astrometry for most lower-mass members of the association and to avoid being spatially biased by the various pointings used for the spectroscopic observations of the low-mass stars we limit our sample in this study to the high-mass members observed by Hipparcos and Gaia.



 §.§ Gaia and Hipparcos astrometry




Astrometric data for our sample comes from Gaia <cit.> data release 1 <cit.>, which contains results from the Tycho-Gaia Astrometric Solution <cit.> that combines astrometry from Gaia, Hipparcos <cit.> and Tycho-2 <cit.>. Of particular relevance for our sample is the subset of 93,635 stars in TGAS where Hipparcos positions at epoch J1991.25 were combined with Gaia observations from the first 14 months of the mission (2014 to 2015). The parallaxes and proper motions calculated from this data are of particularly high precision, with uncertainties of 0.3 mas and 0.06 mas/yr, respectively, though additional systematic uncertainties are still present <cit.>. This data has some limitations since the astrometric solutions were calculated assuming single star behaviour (i.e., no consideration was made for the sources being binary systems) and perspective acceleration was not accounted for <cit.>, but despite this the data should be more than sufficient for our needs, providing over an order of magnitude improvement in the internal kinematics of the association compared to previous works.

Of the 433 stars in our sample Gaia data are available for 258 (60%). Figure <ref> shows colour-magnitude diagrams for the three Sco-Cen subgroups with the sources lacking Gaia astrometry indicated. These sources are predominantly bright (V < 6 mag) and are likely missing from DR1 due to saturation <cit.>. Additionally, sources with extremely red or blue colours, or those with high proper-motions (>3.5 arcsec/yr) are not included in DR1 <cit.>, though this won't include any of our targets. Where Gaia data are unavailable we take parallaxes and proper motions from the revised Hipparcos catalogue <cit.>.



The precision of the astrometric data used is shown in Figure <ref>. The median proper motion precision is 0.06–0.11 mas/yr (≃0.04–0.08 km/s at a representative distance of 150 pc) in α and 0.05-0.08 mas/yr (≃0.04–0.06 km/s) in δ for all the data for the three subgroups. When just Gaia data is considered the median proper motion precision is 0.04–0.06 mas/yr (0.03–0.06 km/s) in α and ∼0.04 mas/yr (0.03 km/s) in δ. The median parallax uncertainty is ∼0.4 mas for all subsets of the data[Gaia DR1 provides considerably improved proper motions compared to Hipparcos because of the increased baseline of observations used, but the parallax precision is not significantly improved because the additional Gaia observations do not cover sufficient time to resolve the parallax motion of stars.]. Note that systematics in the astrometric solution used for DR1 result in an additional parallax systematic error of 0.3 mas that cannot be reduced by averaging parallaxes for groups of stars <cit.>. Correlations between astrometric parameters for a single source can reach significant levels over large areas of the sky so in this work we make use of the full covariance matrix when calculating uncertainties on all astrometric parameters.



 §.§ Radial velocities


Results from Gaia's RV spectrometer are not included in DR1 and so we gathered RVs from the literature. We followed <cit.> in the prioritisation of RVs from different samples, but chose not to exclude RVs from known binaries[Our knowledge of binarity in Sco-Cen is likely to be incomplete and so excluding only the known binaries will not eliminate the impact of binaries on our measurements, only reduce it by an unknown amount. Our approach instead is to model the impact binaries have (such as on the velocity dispersion) and exclude sources known to be binaries where this approach is not feasible.]. Our primary source of RVs was the Magellan/MIKE study of <cit.>, from which 89 matches were found with our catalogue. Our next sources were the Pulkovo Compilation of Radial Velocities <cit.>, from which 118 matches were found, and CRVAD-2 <cit.>, from which 53 matches were found. For the latter catalogue we made sure to remove all sources included from the pre-publication PCRV list that are non-spectroscopic `astrometric' binaries not suitable for kinematic analysis. Finally we took 14 RVs from the Keck and Magellan spectra presented by <cit.>. In total this provided us with RVs for 274 out of 433 stars in our sample. The median RV uncertainty in each of the three subgroups is 2.3 (US), 2.0 (UCL), and 1.5 km/s (LCC).






§ DISTANCE AND 3-DIMENSIONAL STRUCTURE


In this section we calculate the distances to the centres of each subgroup, necessary for calculating the 3-dimensional bulk motions (Section <ref>) and for assessing the evidence for linear expansion (Section <ref>). We also combine the parallaxes with spatial positions to fit 3-dimensional radial profiles to each subgroup, which are necessary for calculating their virial state (Section <ref>).



 §.§ Individual distances



Distances to individual sources and their uncertainties were calculated using the Bayesian inference method of <cit.>, which overcomes the biases and limitations of traditional distance estimates derived from parallaxes. We do this for each star in our sample individually, employing a prior that converges asymptotically towards zero as the distance goes to infinity, P(r) = r^2 exp (-r / L), where L = 1.35 kpc is the scale length parameter, as recommended by <cit.> and <cit.>[The use of such a prior is not recommended when determining the distance to a cluster of stars and therefore we do not use it in Section <ref> when modelling the distances to each subgroup. The effect on our sample of stars of using such a prior over not using a prior is to increase individual distances by an average of ∼0.1 pc.]. To fully sample the posterior distribution we use the Markov Chain Monte Carlo (MCMC) ensemble sampler emcee <cit.>, and use the mode of the corresponding posterior as our distance estimate and the 68% confidence interval as the 1σ uncertainty.

Due to the proximity of Sco-Cen and its large extent on the sky it is often necessary to consider the physical structure and spatial distribution in Galactic cartesian coordinates, XYZ, rather than in α, δ and distance. We therefore calculate XYZ coordinates for all sources using the same method as for the distance and accounting for the correlated uncertainties between α, δ, and ϖ. Note that the distance and XYZ coordinate uncertainties are purely derived from the quoted random errors and do not include the systematic uncertainty.



Figure <ref> shows the 3-dimensional spatial structure of Sco-Cen with the internal structure resolved in all three subgroups for the first time ( was unable to resolve the internal structure of any of the subgroups, while  was only able to resolve the internal structures of UCL and LCC with their improved secular parallaxes). The physical extents of the three subgroups overlap to varying degrees in the three dimensions, being most clearly separated in the Y-Z plane, which is the closest equivalent to the observational l-b plane. A number of spatial outliers from the three subgroups can be seen, especially in the X-Y plane, though their uncertainties are sufficiently high that their membership of Sco-Cen cannot be rejected based solely on this.



 §.§ Modelling the subgroup distances and radial profiles using Bayesian inference



We use Bayesian inference to estimate the distance, size and radial profile of each subgroup. It is necessary to fit these quantities at the same time because the physical extent of each subgroup is a not-insignificant fraction of their distance and thus fits for the two will be correlated. 

The idea of Bayesian inference is to create a parameterised model that can reproduce the observations, and then compare that model (for different sets of parameters) to the observations in a probabilistic way. The aim of this process is to determine which of the various sets of parameters, θ, best explain the observations, d. In Bayes's theorem this is known as the posterior distribution, P(θ | d), and is given by


    P(θ | d) =  P(d | θ)×P(θ)/P(d)


where P(d | θ) is the likelihood model, P(θ) are the priors (which includes our a priori knowledge about the model parameters) and P(d) is a normalising constant.

We use Bayesian inference because it is better to project the model predictions into observational space, where the measurement uncertainties are defined, rather than deriving physical quantities and then trying to calculate appropriate physical uncertainties to use in comparison. This is particularly the case when the observations have highly correlated uncertainties, as in the case of DR1, or when parallaxes are involved <cit.>.



  §.§.§ Method


Our likelihood model constructs the 3D distribution of sources in each subgroup in Galactic cartesian coordinates (XYZ) and then reprojects the positions of individual stars into equatorial coordinates and distances, calculating parallaxes from the latter. We then add measurement uncertainties (including correlated uncertainties) to the parallax, α and δ, randomly sampling these from the observed uncertainty distributions for our sources (a mixture of Gaia and Hipparcos uncertainties).

The structure of each subgroup is modelled using a 3D <cit.> radial profile, where the stellar surface density, Σ(r), at a radius, r, is given by


    Σ (r) = Σ_0 ( 1 + r^2/a^2)^-γ / 2


wherein Σ_0 is the central surface density, a is the scale length, and γ is the power-law index. The scale length is related to the half-mass, or effective, radius according to r_hl = a √(4^1/γ - 1), while the power-law index can be used to calculate the quantity η in the virial equation <cit.>.

We tried two variants of this model, the first with an isotropic radial profile (i.e., single values of a and γ) and the second where the radial profiles in the three dimensions were allowed to vary. Our two models therefore have 5 (distance, central α and δ, a, and γ) and 9 parameters (where a and γ become a_X, a_Y, a_Z, γ_X, γ_Y, and γ_Z).

We adopt liberal priors where possible, to minimise their impact on the fit. For the EFF parameters we use linear priors of 1–50 pc for the scale factors and 2–20 for the power-law indices. For the central subgroup distances we apply linear priors of 50–250 pc, while the priors on α and δ are left unconstrained.

To fully sample the posterior distribution function we use the MCMC ensemble sampler emcee, computing the likelihood model for each point in the parameter space and comparing the model to the observations using an unbinned maximum likelihood test. We found that the model fits would typically converge within 10,000 iterations, with a similar number of iterations necessary to fully sample the posterior distribution function. The posterior distribution function was typically found to follow a normal distribution, and thus the median value was used as the best fit, with the 68% confidence interval used as the 1σ uncertainty.







  §.§.§ Results


The results of the model fits are provided in Table <ref> and illustrated in Figure <ref>. Both the 1D and 3D EFF models show reasonable fits to the data for US and LCC, though the fits for UCL are not as good, most likely because this subgroup is clearly less centrally concentrated, particularly in α, and therefore an EFF model is probably not appropriate. However, our purpose here is to obtain a reliable estimate of the distance to each subgroup and their physical sizes and so these results should be sufficient.

The central coordinates of each subgroup in all three dimensions show good agreement between the 1D and 3D model fit results, with all values within 1σ of their corresponding value. Combining the systematic 0.3 mas parallax uncertainty with the uncertainties already calculated leads to distances for the three subgroups of 143 ± 6, 136 ± 5, and 115 ± 4 pc (though if the systematic uncertainties have been over-estimated then these uncertainties will also be over-estimated). Our values are consistent with the mean distances found by <cit.> of 145 ± 2, 140 ± 2, and 118 ± 2, with our distances slightly smaller than theirs, most likely because of the changes to the membership of Sco-Cen performed by <cit.>, and with slightly larger uncertainties <cit.>. Note that because the physical extent of each subgroup is larger than the uncertainty on the central distances calculated here these distances should only be considered as representative distances for each subgroup and not as distances for the individual stars.

The power-law indices of the 1D EFF profiles are found to be 14.5, 17.6 and 15.2 for US, UCL and LCC respectively, all of which agree with their 3D EFF indices and with each other to within ∼1σ. The 1D EFF scale lengths are 35.8, 69.8 and 50.1 pc for US, UCL and LCC respectively. The scale lengths show good agreement between the 1D and 3D models. Both US and LCC appear close to spherical, with similar scale lengths in all three dimensions that differ by less than 10% from each other and from the 1D scale length.

UCL is the least spherical of the subgroups, elongated in the X direction relative to the Y and Z directions with scale lengths of 74, 57, and 48 pc, respectively. If the subgroups of Sco-Cen do not represent distinct episodes of star formation as some recent studies have suggested <cit.> but are actually made up of smaller substructures, this would explain the non-spherical structure of UCL. The true 3D structure is likely to be much more complex, and also not aligned with the XYZ coordinate system. Since the differences between the 1D and 3D scale lengths is not dramatic, and because we only require a simple estimate of the size of our subgroup, we will use the 1D scale lengths throughout the rest of this work.



§ BULK MOTIONS AND VELOCITY DISPERSIONS



Here we use proper motions and RVs to measure the basic kinematic quantities of each subgroup, their bulk motions, velocity dispersions, and convergence points. These are necessary for a full understanding of the dynamics of the association, to study their relative motions, assess their virial states, and to test whether the subgroups are expanding.



 §.§ Forward modelling the velocity dispersion



To measure the velocity dispersion and bulk motion of each subgroup we use Bayesian inference, as in Section <ref>. We model these quantities in the Galactic cartesian system XYZ, with velocities UVW, and then convert them to the observational frame. We adopt the classical definition of this coordinate system with X directed towards the Galactic centre and increasing in that direction, Y positive in the direction of Galactic rotation, and Z positive towards the north Galactic pole.



  §.§.§ Method


Our forward model begins by modelling a population of stars with 3D coordinates, XYZ, sampled from the observed distribution of these coordinates for each subgroup, as calculated in Section <ref>. The 3D velocities of each star are randomly sampled from a trivariate Gaussian velocity distribution and then reprojected into the observational frame to obtain modelled proper motions and RVs for each star.

Because we are modelling RVs it is necessary to model the contribution that binary orbital motions make to the observed velocity distribution. To do this we simulate a population of randomly-aligned binaries in a manner similar to that of <cit.> and <cit.>, randomly selecting primary and secondary masses, semi-major axes, and eccentricities for each binary. The primary masses were selected from a standard α = 2.3 initial mass function <cit.> in the mass range of 1.4–17.5 M_⊙ (appropriate for the A and B-type stars that constitute our sample) using the equations of <cit.>.



We use the orbital properties of intermediate-mass stars in Sco-Cen determined by <cit.>, who performed an extensive study of the visual, spectroscopic, and astrometric binaries in the association. They inferred a binary fraction of 100%, a mass ratio distribution f_q(q) ∝ q^-0.4, and a semi-major axis distribution of the form f_a(a) ∝ a^-1 (in the range 5 ≥ a / R_⊙≥ 5 × 10^6). The authors were unable to sufficiently constrain the eccentricity distribution from the available observations and so we adopt a flat distribution in the range of 0–1, which <cit.> found to be consistent with their observations. We then calculate instantaneous velocity offsets for the primary star relative to the centre of mass of the system for a random phase within the binary's orbit. Following this method we measure the velocity offsets for 10^6 random binaries, providing us with a distribution from which we can sample. We don't consider triple systems because their properties are poorly constrained <cit.> and are typically hierarchical, meaning that the third star is usually on a wide, long-period orbit that does not have a large velocity.

Stars that are modelled as being in a binary system have velocity offsets added to their RV, randomly sampled from our velocity offset distribution. Finally we add measurement errors for the proper motions and RVs for each star, randomly sampling these from the empirical uncertainty distribution for the stars in each subgroup.

This model has 6 free parameters, the central velocity (U_0, V_0, W_0) and velocity dispersion (σ_U, σ_V, σ_W) in each dimension, for which we adopt flat and wide priors of -100 ≤ v_0 [km/s]≤ +100 for the central velocities and 0 ≤σ[km/s]≤ 100 for the velocity dispersions. As in Section <ref> we sampled the posterior distribution function using an MCMC ensemble sampler, computed the likelihood model for each point in parameter space, and compared the model to the observations using an unbinned maximum likelihood test. The posterior distribution functions followed a broadly normal distribution, albeit with wide tails, and so the median value was used as the best fit and the 68% confidence interval used for the 1σ uncertainty.



  §.§.§ Results




The results are provided in Table <ref> and illustrated in Figure <ref>, with the central velocities, or bulk motions, also shown in Figure <ref>. The bulk motions are different from those presented by <cit.>, particularly in U and V, implying that the difference is probably due to the increased availability and quality of RVs for these stars <cit.>. Our bulk motions are in better agreement with those presented by <cit.>, with most agreeing within 1–2σ, and the remaining differences probably due to the change in membership.

The best-fitting velocity dispersions vary from 0.5 – 2.5 km/s between the three subgroups of Sco-Cen and along the three axes. None of the subgroups have isotropic velocity dispersions, each having one dimension with a significantly larger dispersion than the other two dimensions. For US the velocity dispersion is largest in W, while for UCL and LCC it is largest in U. This is unlikely to be due to the influence of incorrectly-modelled binary systems since RVs contribute most to U for US and UCL, and V for LCC.

Our velocity dispersions are mostly consistent with previous estimates in the literature. <cit.> found that the 1D velocity dispersions for all the subgroups are ≤ 1.0-1.5 km/s, while <cit.> measured velocity dispersions of 1.1–1.3 km/s for the subgroups. Our weighted mean 1D velocity dispersions for the three subgroups are 1.86 ± 0.21 (US), 1.38 ± 0.21 (UCL) and 1.21 ± 0.28 km/s (LCC), for which only the US velocity dispersion is inconsistent with previous findings (to approximately 1.9 and 2.8σ for the two studies cited previously). The disagreement for US can be attributed to a combination of the revised membership of Sco-Cen[We repeated these fits using the <cit.> membership list and find that the change in sample to <cit.> has caused the velocity dispersion to increase by ∼10%.] and our inclusion of RVs to reveal the full 3D velocity dispersions (because of the anisotropy this causes a large change in the mean 1D velocity dispersion).

<cit.> find that the observations of Sco-Cen are best fit with a binary fraction of 100%, though slightly smaller fractions are possible. We explored using binary fractions of 70–100% and found that lower fractions lead to larger velocity dispersions in each dimension, with the largest difference seen along the axis contributing most to the RVs (U for US and UCL, V for LCC). The velocity dispersion increases by ∼10% along that axis for each 10% by which the binary fraction is reduced. Figure <ref> shows the model fit results for a 70% binary fraction, with little difference in the fit quality compared to that for 100% binaries.



 §.§ Convergent point motion



Nearby groups of stars with a common space motion have proper motions that appear to `converge' towards a single point on the sky. This `convergent point' can be used to establish membership of these groups, measure their mean distances, and test models for their expansion (Section <ref>). In this section we measure the convergent points of the three subgroups using both the refurbished <cit.> method <cit.> and directly from the 3D bulk motions.



  §.§.§ Method


Jones's method iterates over a hemispherical grid of points to find the convergent point coordinates that minimises the proper motions perpendicular to the great circle joining each star to the convergent point. The best fitting convergent point is that which minimises the quantity ∑ t^2_⊥, where t_⊥ = μ_⊥ / σ_⊥ is the ratio of the perpendicular proper motion to its uncertainty, and the sum is over all stars in the group. This method was refined by <cit.> who changed this quantity to t_⊥ = μ_⊥ / √(σ^2_⊥ + σ^2_v), recognising that moving groups and associations have an intrinsic, one-dimensional velocity dispersion, σ_v, that prevents all their proper motions from being directed exactly towards the convergent point and should therefore be accounted for. We adopt this approach, using an MCMC method to search for the convergent point, instead of a grid-based technique.

Since the t_⊥^2 value can be treated as the classic χ^2 statistic it can be used to assess the reliability of a fit using the probability, ϵ, that t_⊥^2 will exceed the observed value of t_⊥^2 by chance even for a correct model. Therefore, if the value of t_⊥^2 is too high the probability ϵ will be below the threshold value, which <cit.> recommend setting to a value of 0.954 (giving a ∼5% probability of falsely rejecting the null hypothesis). If ϵ is below this value then the star with the highest value of | t_⊥ | is removed from the sample and the fitting process restarted. This process continues until a set of stars returns a convergent point with a sufficiently low t_⊥^2 value.



  §.§.§ Results


Convergent points were found for all three subgroups of Sco-Cen using the Jones method and the mean 1D velocity dispersions determined in Section <ref>[Since the velocity dispersions are not isotropic one should ideally use the component of the velocity ellipsoid perpendicular to the great circle connecting the subgroup with the convergent point. However this calculation would be complex to perform for every convergent point considered, and since the exact velocity dispersion used does not significantly affect the final convergent point calculated we have instead used the mean 1D velocity dispersion.]. For US a convergent point for the entire subgroup was found with a χ^2 probability of 53%, rising to 95.4% after removing the 20 stars with the highest | t_⊥ | values. The two convergent points vary by ∼0.05 deg, significantly less than their uncertainties. Higher velocity dispersions of 2.0 or 2.5 km/s resulted in convergence after rejecting only 6 or 0 stars, but these do not show significant variations from those already determined, so we choose the convergent point calculated using the mean 1D velocity dispersion of 1.86 km/s. Figure <ref> shows this convergent point compared to those determined previously by <cit.> and <cit.> using just proper motions, which agree well despite the change in membership of the subgroup, and the convergent point determined from the mean space motion, which agrees to within 2σ.

For UCL a convergent point for the entire subgroup was found with a χ^2 probability of 98% without discarding any stars, and is therefore much better constrained than that of US. For LCC a convergent point for the entire subgroup was found with a χ^2 probability of 20%, rising to 95.4% after removing the 10 stars with the highest | t_⊥ | values. A higher velocity dispersion of 1.5 km/s would result in an immediate fit with a χ^2 probability of >99% for all the stars in the subgroup. However, these three convergent points all agree within ∼0.05 deg, significantly less than the uncertainty, so the fit with the mean velocity dispersion of 1.21 km/s was used. For both UCL and LCC our convergent points agree with previous estimates in the literature to within 1–2σ, with any differences most likely due to changes in the subgroup membership. The agreement between the convergent points calculated using only proper motions and those calculated using the mean space motions are not as good however, disagreeing by more than 3σ. These disagreements suggest that some degree of expansion (or contraction) may be present.






§ EXPANSION OF THE ASSOCIATION



In this section we explore whether the subgroups are in the process of expanding, and if so whether they are expanding from a single point. Due to the proximity of Sco-Cen we cannot use proper motions alone to study the expansion of the association on the plane of the sky since a radial motion of the association towards (or away from) the observer will cause an apparent expansion (contraction) of the association, even when none is actually present <cit.>. We must therefore either combine the proper motions with RVs, which we can do using many different methods, or use the bulk radial motion to correct this virtual expansion / contraction.



 §.§ Linear expansion





If an association is expanding then the observed RVs of individual stars will vary from those predicted from the moving group method for parallel motion by an amount that varies with distance. Stars on the near-side (far-side) of the association will have smaller (larger) RVs than would be predicted by the moving group method as these stars are moving away from the centre of the association and towards (away from) us.

In the <cit.> linear expansion model individual RVs, v_rad, are predicted to follow the relation


    v_rad = v^' cos λ^' + κ d + K


where λ^' is the angular separation between star and convergent point (calculated using just the proper motions, see Section <ref>), d is the distance to the star in pc, κ is an expansion term, and K is a zero point (that can represent gravitational redshift or convective blueshift). v^' is the barycentric velocity of a star in the association, given by


    v^' = μ_v A /ϖ sin λ^'


where μ_v is the proper motion in the direction of the convergent point (in mas/yr), ϖ is the parallax (in mas), and A = 4.74 km yr s^-1 is the astronomical unit in km yr s^-1.

An expanding association will have κ > 0, from which an expansion age, τ = (γκ)^-1, (in Myrs) can be derived (γ is a conversion factor of 1.0227 pc Myr^-1 km^-1 s).

Figure <ref> shows the difference between the observed RVs and those predicted from the moving group method plotted against the distance to each star. The solid lines show the best-fitting linear relationships between these quantities with 3σ outliers on the ordinate excluded. Fits were obtained using the MCMC code emcee, accounting for the distance uncertainties by randomly varying the distances according to their uncertainties. Given the possibility that the uncertainties on the ordinate have been underestimated, these fits were obtained by introducing a factor, f, by which the uncertainties have been underestimated <cit.> and then marginalising over this parameter to obtain the uncertainties on the resulting fit parameters.

The best-fitting slopes for the three subgroups are -0.023^+0.019_-0.024 (US), -0.006^+0.016_-0.017 (UCL), and -0.034^+0.030_-0.031 km s^-1 pc^-1 (LCC). All three subgroups exhibit marginally negative slopes that are consistent with a slope of zero within a little over 1σ. The results for US are consistent with the fit of κ = -0.01 ± 0.04 km s^-1 pc^-1 found by <cit.>. Based on these slopes and their confidence intervals we can rule out expansion ages equal to or smaller than the median ages of 10, 16 and 15 Myr found by <cit.> with confidences of approximately 7, 4, and 3σ for the three subgroups respectively. The fits were repeated with known binaries excluded (10, 1, and 2 sources in the three subgroups), with no significant difference.



 §.§ 3D linear expansion





Using the three-dimensional positions and velocities (for the stars with RVs) we can search for evidence of expansion in each of the three axes of the Galactic Cartesian coordinate system by plotting velocity versus position and searching for evidence of a positive slope between the two quantities that would provide evidence for expansion along that axis.

Figure <ref> shows the velocities UVW plotted against the positional coordinates XYZ for each of the three subgroups. We fit linear relationships to these quantities using MCMC, as above (with 3σ outliers on the ordinate excluded), to estimate the slopes, κ_X, κ_Y, κ_Z, in the three dimensions. For US the slopes are (κ_X, κ_Y, κ_Z) = (-0.01 ± 0.02, 0.07 ± 0.03, -0.03 ± 0.03) km s^-1 pc^-1, for UCL they are (0.01 ± 0.01, 0.05 ± 0.02, 0.01 ± 0.01) km s^-1 pc^-1 and for LCC they are (0.02 ± 0.02, 0.05 ± 0.03, 0.01 ± 0.01) km s^-1 pc^-1.

In the X and Z dimensions all three subgroups show either marginally negative or positive slopes that are consistent with a slope of zero within 1σ, implying no evidence for either expansion or contraction. Interestingly however all three subgroups show significant evidence for expansion along the Y axis, with significances of approximately 3σ for US and UCL and 2σ for LCC. The fits were repeated with known binaries removed, as above, but there was no significant difference in the results.



 §.§ Expanding versus non-expanding convergent points





Our calculation of the convergent points that used only proper motions were based on the assumption that the associations were not expanding. However, a group of stars that has a small linear expansion will appear to converge to a point further away (i.e. with higher λ) than for a group of stars without expansion. This could explain the disagreement between the convergent points calculated using just proper motions and those calculated from the centroid space motion. To test whether the subgroups are expanding we add linear expansion to the centroid space motions (adding a term κ d to the RVs, as per equation <ref>) and compare these with the proper motion convergent points.

Figure <ref> shows the centroid space motion convergent points for various expansion ages in the range 0–100 Myr. For US adding expansion actually worsens the agreement between the two convergent points, which argues against this subgroup exhibiting significant levels of expansion. For UCL small levels of expansion lead to an improved agreement between the convergent points for expansion ages ≳20 Myr, though this is larger than its estimated age of 15 ± 3 Myr <cit.>. For LCC an expansion age of 18–25 Myr leads to a good agreement between the two convergent points, which is in reasonable agreement with the estimated age of 16 ± 2 Myr from <cit.>. Based on this analysis both US and UCL appear to be inconsistent with being in a state of expansion, whilst LCC is consistent with being an expanding association.



 §.§ Tracing back individual stellar motions



An alternative method to assess the expansion of a group of stars is to trace back their individual motions in 3D and quantify the spatial extent at each time-step. For OB associations a simple traceback is generally valid because the densities are not high enough for there to be a significant number of close encounters and scatterings that would invalidate this approach. We do this considering both linear trajectories, i.e., X(t) = X_0 + γ U t for the X/U dimension (where the time, t, is in Myr and γ = 1.022 s pc km^-1 Myr^-1), and using an epicycle approximation.

Figure <ref> shows the 1σ dispersions (16th to 84th percentiles) in X, Y, Z, and in 3D using linear trajectories. For all three subgroups the dispersion in X steadily increases as one goes further back into the past. The dispersion in Y is similar for US, but for UCL and LCC it is approximately constant over the past ∼5 and ∼10 Myr, respectively, before showing signs of increase. The dispersion in Z varies the most between subgroups. For LCC it is steadily increasing, for UCL it is either unchanged or slightly decreasing over the last ∼5 Myr and then increases, while for US it shows evidence for a slight decrease over the last ∼3 Myr before increasing. The 3D dispersions steadily increase as one goes back into the past, with only UCL showing evidence for a minimum not at the present day.

For a more accurate traceback of the stellar motions we use the epicycle approximation, employing the orbital equations from <cit.>. We use the Oort A and B constants from <cit.>, the local disc density from <cit.>, the local standard of rest velocity from <cit.> and a solar Z distance above the Galactic plane of 17 pc <cit.>. Figure <ref> shows the 1σ spatial dispersions using the epicycle approximation. The results do not significantly differ from those calculated using linear trajectories, with all subgroups and all dimensions showing a broad increase in dispersion as we look further back in time. The largest increases are universally in the X dispersion, with the smallest increases in the Y dispersion (and in fact for UCL and LCC the dispersions in Y are fairly steady over the lifetimes of the subgroups), and significant increases in the 3D dispersions.

Uncertainties in velocity can cause an expanding subgroup to appear larger in the past as the measured velocities differ from a true expansion pattern. To test the significance of this effect we performed a Monte Carlo simulation by randomly varying the measured velocities in accordance with their uncertainties and then recalculating the traceback motions. We repeated this 10,000 times and in Figure <ref> we show the 90% confidence interval around the 3D dispersions. The uncertainties in the size of each subgroup in the past are notably smaller than their actual size, suggesting that this effect is not due to measurement error.

To conclude, using either linear trajectories or an epicycle approximation and defining the size of the subgroups from their 1σ spatial dispersions there is no evidence that the three subgroups are expanding or that they have ever had a more compact configuration in the past. It is worth noting that <cit.> explored the validity of similar methods to this by producing synthetic proper motion observations of expanding OB associations and found that the smallest size of the OB association predicted from a traceback analysis using proper motions alone was always overestimated. The improvements in proper motion precision from Gaia since this work, as well as the increased availability and precision of RVs that allow this analysis to be performed in 3D rather than in the plane of the sky (and thus overcome the effects of radial streaming), do resolve many of the issues raised by <cit.> however.





 §.§ Corrected proper motion vector maps



If a group of stars has a non-zero radial motion this will cause an apparent expansion or contraction of the group in their proper motions as their spatial extent on the sky changes <cit.>. If the bulk RV of the stars is known one can perform a correction for this virtual expansion / contraction, effectively calculating the contribution to the proper motions due to the radial motion of the association and then subtracting these from the observed proper motions to provide corrected proper motions. To do this we use the equations in Appendix A1 of <cit.>, using the median RV of stars in each subgroup (to minimise the impact of binarity).

The corrected proper motion vector maps for the three subgroups are shown in Figure <ref>. None of the three subgroups show evidence for a coherent expansion pattern. UCL does appear to show some coherent patterns consistent with expansion, with stars in the north-east and south-west appearing to move away from the centre of the subgroup, though stars elsewhere in the group do not show this behaviour. We will argue, in Section <ref>, that these patterns are due to the kinematic substructure within the subgroups.

To quantify the expansion of each subgroup we follow the method of <cit.> and calculate, for each star, the amount of kinetic energy in the radial and transverse directions (using the subgroup centres calculated in Section <ref>), dividing the former between expansion and contraction, and then sum these for each subgroup to provide global energy fractions (for simplicity we assume all stars have equal masses). If a subgroup was expanding radially from a single point we would expect the majority of kinetic energy to be in the radial direction and for the stars to be moving outwards.

We find that the kinetic energy is divided almost equally between the radial and transverse directions, with 50% (US), 66% (UCL), and 45% (LCC) of energy in the radial direction. Only UCL exhibits a preference for kinetic energy in the radial direction. In the radial direction there is a preference for the motions of stars to be directed away from the centres of each subgroup, with 59% (US), 67% (UCL) and 90% (LCC) of the kinetic energy in the direction of expansion. This suggests that, while there isn't evidence for coherent expansion patterns in the correct proper motions of stars in each subgroup, they do appear to be expanding.






§ KINEMATIC SUBSTRUCTURE IN THE ASSOCIATION





Kinematic substructure is the observed tendency for stars in the same area of a star forming region or association to have motions that are more similar to their neighbours than to other stars. It may manifest as sizeable substructures within the association <cit.> or as smaller, more numerous substructures <cit.>. In the case of OB associations, kinematic substructure is believed to be a remnant of the primordial phase-space structure that was present when the association formed <cit.>. In this Section we explore the evidence for kinematic substructure in Sco-Cen.

Figure <ref> shows a corrected proper motion vector map for Sco-Cen with the vectors colour-coded by the position angle of their motion on the celestial sphere <cit.>. In this figure the three subgroups stand out relatively clearly due to their different kinematics, with some overlap between the groups that suggests their membership could benefit from some adjustments.

In addition to the three subgroups there is also a group of stars visible in this diagram with above average velocities moving in the direction of decreasing latitude and longitude (they appear a creamy-orange colour in this diagram). These stars have different proper motions from the stars they are projected against because, despite having similar kinematics, they are in the foreground at distances of 50–100 pc. Four of these stars are significantly in the foreground (d < 70 pc) and so are unlikely to be members of Sco-Cen; HIP 68413, HIP 69995, HIP 70931 and HIP 76063. We can find no evidence in the literature that these stars are particularly young <cit.> and therefore suggest that these are field stars that should be removed from the membership of Sco-Cen.

Figure <ref> shows corrected proper motion vector maps for the three subgroups of Sco-Cen. From this it is evident that the proper motions of stars in each subgroup are not random but are correlated with position. Each subgroup shows evidence for multiple groups of stars each in their own area of the sky and each with stars that are all moving in a similar direction. UCL and LCC are particularly noteworthy as they each exhibit a small number of groups (∼4–5) each with 10–30 members and in a distinct part of the association. In UCL and LCC there are hints of expansion in the α axis in that stars in the east of the subgroups are moving eastwards and stars in the west are moving westwards. Since this axis is close to the Galactic Y axis (in UCL and LCC) this is probably the same feature we observed in Section <ref>. In US the kinematic substructures are not as distinct and the subgroup appears to be made up of a larger number of smaller substructures, similar to the kinematic substructure observed in Cygnus OB2.



 §.§ Quantifying kinematic substructure




  §.§.§ Method


In an attempt to quantify this substructure we follow <cit.> and use Moran's spatial correlation index I <cit.>, which quantifies the overall degree of similarity between spatially close regions with respect to a numeric variable (in this case the proper motion in either dimension). The expectation value for no spatial correlation is I_0 = -1 / (n-1), with I > I_0 indicating positive spatial correlation and I < I_0 indicating negative spatial correlation. We use this instead of Geary's C index <cit.>, since Moran's I statistic is a global measure that can identify and quantify large-scale spatial correlation (as appears to be evident in Figure <ref>), while Geary's C is a local measure of spatial correlation. For brevity we don't provide quantification of the index, but instead refer the reader to these papers and <cit.>.

We calculate the index for the proper motions in both directions for all three subgroups of Sco-Cen, using all the stars in our sample. We did this in proper motion space rather than in cartesian UVW space as this would have required the use of individual RVs that are distorted by binary motions leading to an artificial `smearing' of the phase-space structure in this direction. We take into account the effect that the projection of motions onto a curved surface have on the measured kinematic substructure by creating synthetic populations of each subgroup (consisting of stars at the observed positions and distances of each star but with proper motions calculated from a uniform UVW space motion with an added velocity dispersion equal to those measured in Section <ref>). We then calculated Moran's I index for the simulated dataset and repeated this 10,000 times for each subgroup to determine the typical level of kinematic substructure that would be measured solely due to this.



  §.§.§ Results


For US we calculate values of I_α = 0.135 ± 0.048 and I_δ = -0.013 ± 0.048, indicating large and significant (3.0σ) positive spatial correlation (nearby regions tend to exhibit similar velocities) for μ_α, and a very weak measure of negative spatial correlation (nearby regions exhibit dissimilar velocities) for μ_δ that is consistent with no spatial correlation. The large difference between the two dimensions is probably because the compact kinematic substructures evident in US (such as the groups at α, δ = 243, -20 or 237, -25) have most of their motion in α leading to significant spatial correlation in α but not δ. Our Monte Carlo simulations indicate that the projection of a bulk motion into proper motions would give typical values of I_α = 0.002 (0.24σ) and I_δ = 0.006 (0.31σ), which at least for the α dimension is significantly smaller than the measured value.

For UCL we calculate values of I_α = 0.000 ± 0.022 and I_δ = -0.034 ± 0.022, indicating weak spatial correlation in both dimensions that is  consistent with no spatial correlation in μ_α or has a low significance (1.3σ) of negative spatial correlation in μ_δ. The latter is likely due to the presence of an extended group of rapidly moving stars in pink in Figure <ref> that are moving in a southerly direction. The presence of these stars in between stars moving in different directions will cause negative spatial correlation, particularly in μ_δ. Our Monte Carlo simulations predict that with no kinematic substructure we should measure substructure with significances of 0.21 and 0.11σ, which in α is consistent with the level measured.

For LCC we calculate I_α = 0.024 ± 0.019 and I_δ = 0.029 ± 0.019, indicating positive spatial correlation in both dimensions at 1.6 and 1.9σ, respectively. Our Monte Carlo simulations predict that with no kinematic substructure we should measure substructure with significances of 0.40 and 0.31σ, which slightly reduces the impact of the kinematic substructure we have measured.

Taken together these results suggest that while there is kinematic substructure in each of the three subgroups of Sco-Cen, the level of kinematic substructure varies between subgroups and is lower than that measured in Cyg OB2 <cit.>.





§ DISCUSSION




 §.§ The virial state of Sco-Cen's subgroups



To determine the virial state of each subgroup of Sco-Cen we use the virial equation, which in its 3D form is given by


    σ_3D^2 = G M_vir/2 r_vir


where σ_3D is the 3D velocity dispersion, G is the gravitational constant, M_vir is the virial mass, and r_vir is the virial radius. We substitute the parameter η = 6 r_vir / r_eff <cit.>, where r_eff is the effective (or half-light) radius, giving


    M_vir = η σ_3D^2   r_eff/3G


where the parameters η and r_eff can be calculated from the <cit.> parameters γ and a (see Section <ref>). For the large γ values calculated for the subgroups we find η values of 9.1 (US), 9.0 (UCL), and 9.1 (LCC)[The dependence of η on γ is very weak at large γ values and therefore the uncertainties on η are almost zero.], while the effective radii, r_eff are 11.3^+0.5_-0.3 (US), 20.0^+0.7_-0.7 (UCL), and 15.5^+0.9_-0.8 pc (LCC).

Combining the effective radii, 3D velocity dispersions (from Table <ref>), and η values, we calculate virial masses of 8.2^+1.5_-1.2× 10^4 (US), 8.3^+1.8_-1.5× 10^4 (UCL), and 5.1^+2.9_-1.3× 10^4 M_⊙ (LCC). These values can be compared to the total mass of stars in each subgroup to estimate the virial state of each group. <cit.> estimate US to have a total mass of ∼2000 M_⊙, while <cit.> estimate UCL and LCC to consist of approximately 2200 and 1200 stars, respectively, which equates to total stellar masses of about 1250 and 700 M_⊙ <cit.>. All of these estimates are significantly smaller than the calculated virial masses, by at least an order of magnitude, implying that the three subgroups are in a supervirial state, as expected given their low stellar density. We would therefore expect the subgroups to be in a state of expansion.

It is interesting to note that the virial mass estimated for US is very close to the mass of atomic H i in the shell surrounding the subgroup, which <cit.> estimate to have a mass of 8 × 10^4 M_⊙. They argue that the shell consists of swept-up material from the primordial GMC from which US formed, which suggests that the subgroup could have been born close to virial equilibrium within the GMC (if the stars would have felt the gravitational potential of the entire GMC). The masses of the shells surrounding UCL and LCC, 3 × 10^5 and 1 × 10^5 M_⊙ respectively, are more massive than their virial masses, implying that the subgroups might have initially been in a subvirial state.



 §.§ Expansion of the Sco-Cen subgroups


Ever since <cit.> we have known that OB associations are unbound and will therefore expand. This has lead many authors to postulate that they were more compact in the past and have been expanding for a while <cit.>. These ideas lead to suggestions that OB associations may be the expanded remnants of compact star clusters <cit.>, disrupted by processes such as residual gas expulsion <cit.>.

In Section <ref> we used five different methods to determine whether the subgroups of Sco-Cen are expanding, and if so whether they are expanding from a previously compact configuration such as a dense star cluster. Many of these methods assume that the cluster would expand radially (such as the expansion models explored in Sections <ref> and <ref>), while others better account for the Galactic potential (Section <ref>) or provide an overall view of the stellar motions (Section <ref>). Some of these methods use individual RVs (where available), while others use the bulk RV (and are thus less affected by binary motions).

None of the methods provide evidence that the subgroups are expanding from a more compact configuration. The linear expansion model considered in Section <ref> finds negative slopes that imply contraction rather than expansion, though all three are consistent with no expansion or contraction to ∼1σ. The 3D expansion model (Section <ref>) again shows slopes consistent with no expansion or contraction, except along the Y axis where all three subgroups exhibit evidence for expansion to a confidence of 2–3σ. Interestingly, <cit.> found for the β Pictoris moving group that the most significant evidence for expansion also came in the Y direction. This could be caused by galactic shear since the surface densities of OB associations like Sco-Cen are low enough for shear to be effective, though the timescales involved are probably too short <cit.>. Alternatively this shear pattern may have been imprinted on the molecular gas in the primordial GMC and then inherited by the OB association that formed. Since the OB association does not appear to be dynamically evolved such a kinematic pattern could have survived the early evolution of the system.

By tracing back the stellar motions (Section <ref>) there is evidence that the subgroups were actually larger in the past than they are today. This result is supported by the corrected proper motion vector maps in Figure <ref> that show no preference for motions radially away from the subgroup centres, though amongst the radial motions there is a preference for expanding motions over contracting motions. Similar results have been found for other associations and moving groups, such as the β Pictoris moving group <cit.> and Tuc-Hor <cit.>, both of which do not appear to have been significantly smaller in the past.

In summary, while the three subgroups are all gravitationally unbound and will therefore expand in the future, they all lack evidence for having been in a more compact configuration in the past. This goes against the view that associations expand as they age from an almost universal compact configuration <cit.>. This is particularly clear in the case of the three subgroups of Sco-Cen, which <cit.> suggest had an original size of 1–3 pc, approximately an order of magnitude smaller than their current size. The kinematics of stars in the subgroups shows that this is definitely not the case. Furthermore there is no evidence for a coherent radial expansion pattern amongst the members of each subgroup and the kinematic substructure evident in Figure <ref> suggest a more complex expansion pattern consisting of multiple expanding substructures within each subgroup.



 §.§ Kinematic substructure and the formation of Sco-Cen


Molecular clouds are known to be highly substructured, both spatially <cit.> and kinematically <cit.>, with this hierarchical or fractal structure passed on to the formed stars <cit.> and often evident in their initial spatial distribution <cit.>. Dynamical interactions between groups of stars can erase this substructure <cit.>, a process that acts on the dynamical timescale of a region. If the dynamical timescale is longer than the age of the region then the primordial substructure can be preserved <cit.>, though dynamical timescales can change as a region expands or collapses (so an instantaneous measure of the dynamical timescale isn't always so revealing). The amount of substructure present in a region can constrain the past dynamical evolution of a group of stars, since it requires that a region cannot have had periods with a short dynamical time otherwise dynamical mixing would have erased its substructure. Compact star clusters have significantly shorter dynamical timescales than extended OB associations, so the presence of physical or dynamical substructure in a region means it can never have been in a highly compact phase in the past <cit.>.

The three subgroups of Sco-Cen each display evidence for kinematic substructure, as has been observed in other OB associations <cit.> and across large star forming complexes <cit.>. US exhibits kinematic substructure similar to that seen in Cyg OB2, with many small subgroups containing 5–10 stars from our sample, and some stars that don't appear to be part of any visible substructures, while UCL and LCC exhibit much larger substructures with 10s of stars in each group, but only 2–4 clear groups in each subgroup. These substructures may be responsible for the complex star formation history revealed by <cit.>, with the observed age spreads within each subgroup due to the combination of individual substructures within it, each of which may not have a significant age spread. We therefore agree with the conclusions of those authors that the current division of Sco-Cen into three subgroups is too simplistic and at least UCL and LCC should be broken down into smaller subgroups based on their ages or kinematics.

A number of different formation scenarios have been considered for Sco-Cen, ranging from the sequential star formation process originally proposed by <cit.> and extended by <cit.>, the impact of a high-velocity cloud on the disk <cit.>, to its formation in a large ring-like structure containing other young groups of stars known as the Gould Belt <cit.>. The high-velocity cloud scenario and the Gould Belt model both predict stellar motions that are opposite to those observed and so have been ruled out <cit.>. The sequential star formation model suggests that star formation began in UCL and propagated to LCC, US, and then the ρ Ophiuchus star forming region through triggering by supernova shockwaves. The presence of kinematic substructure and age spreads within the subgroups that suggests they would be better subdivided into smaller groups with a more complex star formation history complicates this model, but doesn't necessarily disprove it.



§ CONCLUSIONS


Our kinematic study of the Sco-Cen OB association using Gaia DR1 parallaxes and proper motions has lead to the following findings.



  * We use Bayesian inference and forward modelling to calculate 3D velocity dispersions for the three subgroups. All three have non-isotropic velocity dispersions, suggesting they are not dynamically relaxed and likely have not been in the past. The 3D velocity dispersions are 3.20^+0.22_-0.20 (US), 2.45^+0.20_-0.20 (UCL), and 2.15^+0.47_-0.24 km/s (LCC). These imply virial masses that are over an order of magnitude larger than the known stellar mass, confirming the subgroups are gravitationally unbound.

  * We have explored multiple methods for testing if and how the subgroups are expanding, including the <cit.> linear expansion model, 3D linear expansion, a comparison of expanding and non-expanding convergence points, tracing back the individual stellar motions to identify the smallest size of each subgroup in the past, and studying corrected proper motion vector maps. The kinematic data are inconsistent with the subgroups being the expanded remnants of individual star clusters, with no coherent expansion pattern evident and no evidence that the subgroups had a more compact configuration in the past.

  * The subgroups all show evidence for kinematic substructure, which we quantify using a spatial correlation test. US appears to consist of multiple small substructures, similar to that observed in Cyg OB2 <cit.>, thought not as strong, while UCL and LCC appear to be made up of a smaller number of larger substructures with a lower significance. The presence of these substructures places constraints on the past structure and dynamical evolution of the subgroups.


To conclude, the three subgroups of Sco-Cen are not the expanded remnants of individual star clusters and instead are likely to have formed with considerable physical and kinematic sub- structure, such as from numerous smaller clusters. Much of this kinematic substructure remains today, most likely because the subgroups have not undergone a densely clustered phase during which the substructure would have been erased. US retains considerable substructure that may be very similar to its primordial structure, while UCL and LCC appear to consist of a few larger substructures that may indicate either that they were born as such groups or that these substructures have undergone some mergers already. Combined with the recent evidence from <cit.> showing that Sco-Cen has considerable substructure in its age distribution (which is likely to be closely related to its kinematic substructure) our results suggest that the subgroups did not form as individual bursts of star formation but are instead composed of multiple smaller subgroups, each of which probably formed in a different place and at a different time.

Combined with recent studies <cit.> our results suggest that the classical view of OB associations (and moving groups) being the expanded remnants of star clusters is incorrect. Instead these groups appear to have formed with more extended and substructured spatial and kinematic distributions. This implies that the assumption that clusters and associations expand uniformly as they age is not true <cit.> and that star cluster disruption mechanisms, such as residual gas expulsion <cit.>, are not necessary for dispersing the majority of young stars into the Galactic field. This has significant implications for the formation of the Galactic field and the birth environment of stars and proto-planetary disks.



§ ACKNOWLEDGMENTS


NJW acknowledges an STFC Ernest Rutherford Fellowship (grant number ST/M005569/1). Part of this work was completed at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. The authors would like to thank Mark Pecaut and Rob Jeffries for discussions and comments on this paper, as well as the anonymous referee for comments and suggestions that improved the quality of this paper. This work has made use of data from the ESA space mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of the Simbad and Vizier catalogue access tools (provided by CDS, Strasbourg, France), Astropy <cit.> and TOPCAT <cit.>. 

mn2e



