# The kinematics of the Scorpius-Centaurus OB association from Gaia DR1

The kinematics of the Scorpius-Centaurus OB association from Gaia DR1 Abstract 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 these data to confirm that the three subgroups of Sco-Cen are gravitationally unbound and have non-isotropic velocity dispersions, suggesting that 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 to have been born highly substructured, with multiple small-scale star formation events contributing to the overall OB association, and not as single, monolithic burst 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 1 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 discs (e.g. O'dell & Wen 1994; Wright et al. 2012; Guarcello et al. 2016), while close encounters in dense clusters can affect the properties of binary and multiple systems (Kroupa, Petr & McCaughrean 1999; Marks & Kroupa 2012; Parker & Goodwin 2012), protoplanetary discs (Scally & Clarke 2001; Olczak, Pfalzner & Eckart 2008; Rosotti et al. 2014), and influence the formation of planetary systems (Adams et al. 2006; Parker & Quanz 2012). Stars born in low-density environments (e.g. Bressert et al. 2010; Wright et al. 2014) 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 per cent of stars are found in bound clusters (Lada & Lada 2003). 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 (which can account for well over half the mass of the system; Elmegreen et al. 2000) then the cluster becomes supervirial and will expand and disperse (e.g. Tutukov 1978; Hills 1980; Lada, Margulis & Dearborn 1984; Baumgardt & Kroupa 2007). During the expansion the system would briefly be visible as a low-density group of young stars known as an association (e.g. Ambartsumian 1947; Blaauw 1964a; Brown, Dekker & de Zeeuw 1997; Kroupa, Aarseth & Hurley 2001), before dispersing into the Galactic field. This theoretical framework has existed for many decades, and while many observations support aspects of this model (e.g. Lada & Lada 2003, and references therein) it has been difficult to test kinematically due to the lack of high-precision proper motions necessary to assess evidence for expansion. Wright et al. (2016) 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 (Wright et al. 2014) 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 (Gaia Collaboration 2016b) 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 (de Zeeuw et al. 1999) and being at least 100 pc across, the association spans almost 90 deg on the sky (see Fig. 1). The association consists of three subgroups (first defined by Blaauw 1946) 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 (Pecaut & Mamajek 2016). Because of its age, the region has likely already seen several supernovae (e.g. Breitschwerdt et al. 2016), one of which may have been responsible for the ejected runaway star ζ Oph (de Geus 1992; Hoogerwerf, de Bruijne & de Zeeuw 2000). The association has been the focus of numerous studies of its low-mass population (e.g. Preibisch & Zinnecker 1999; Mamajek, Meyer & Liebert 2002), with studies focusing on the binary properties (e.g. Kouwenhoven et al. 2007; Janson et al. 2013), circumstellar discs (e.g. Carpenter et al. 2006; Chen et al. 2011), stellar ages (e.g. Pecaut, Mamajek & Bubar 2012; Pecaut & Mamajek 2016), eclipsing binaries (e.g. Kraus et al. 2015), and low-mass stars and brown dwarfs (e.g. Lodieu 2013). Its high-mass population has been studied by numerous authors, including comprehensive studies by Blaauw (1946, 1964b) and later, using Hipparcos astrometry, by de Zeeuw et al. (1999), de Bruijne (1999b), and Rizzuto et al. (2011). 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. 2 OBSERVATIONAL DATA The sample of Sco-Cen members used for our kinematic study is based on the Bayesian membership analysis of Rizzuto et al. (2011), who used the updated Hipparcos catalogue (van Leeuwen 2007) to revise the kinematic selection of de Zeeuw et al. (1999). 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 (CRVAD-2; Kharchenko et al. 2007)1 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 Rizzuto et al. (2011) were rejected based on their parallax or RVs, two parameters that were not used by de Zeeuw et al. (1999) 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. Chen et al. (2011) 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 Mamajek et al. (2002), Chen et al. (2011) and Pecaut et al. (2012) that were included in the original sample of de Zeeuw et al. (1999) were not selected by Rizzuto et al. (2011) as members, reinforcing the validity of their membership selection. Pecaut et al. (2012) recommend excluding a number of stars based on their kinematics, one of which is retained in the Rizzuto et al. (2011) catalogue, HIP 56227. We choose to retain this star for the time being and minimize 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 Fig. 1. For the handful of stars that cannot be assigned to a subgroup based on the division of de Zeeuw et al. (1999), we put them in the nearest subgroup on the plane of the sky. While Rizzuto et al. (2011) 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. Figure 1. View largeDownload slide Spatial distribution of Sco-Cen members (Rizzuto, Ireland & Robertson 2011) projected on to an inverted H α image from Finkbeiner (2003). Large circles represent sources with Gaia DR1 proper motions and parallaxes, while small circles are sources not detected by Gaia for which Hipparcos proper motions and parallaxes were used. The dashed lines show the division of the association into the three subgroups according to de Zeeuw et al. (1999), with the symbols coloured according to which of the three subgroups they nominally belong: red for US, green for UCL, and blue for LCC. Sources outside of these boundaries have been assigned to the nearest subgroup. Figure 1. View largeDownload slide Spatial distribution of Sco-Cen members (Rizzuto, Ireland & Robertson 2011) projected on to an inverted H α image from Finkbeiner (2003). Large circles represent sources with Gaia DR1 proper motions and parallaxes, while small circles are sources not detected by Gaia for which Hipparcos proper motions and parallaxes were used. The dashed lines show the division of the association into the three subgroups according to de Zeeuw et al. (1999), with the symbols coloured according to which of the three subgroups they nominally belong: red for US, green for UCL, and blue for LCC. Sources outside of these boundaries have been assigned to the nearest subgroup. 2.1 Gaia and Hipparcos astrometry Astrometric data for our sample comes from Gaia (Gaia Collaboration 2016b) data release 1 (DR1; Gaia Collaboration 2016a), which contains results from the Tycho-Gaia Astrometric Solution (TGAS; Michalik, Lindegren & Hobbs 2015; Lindegren et al. 2016) that combines astrometry from Gaia, Hipparcos (ESA 1997) and Tycho-2 (Høg et al. 2000). 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–2015). The parallaxes and proper motions calculated from these data are of particularly high precision, with uncertainties of 0.3 mas and 0.06 mas yr−1, respectively, though additional systematic uncertainties are still present (Gaia Collaboration 2016a; Lindegren et al. 2016). These data have 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 (Lindegren et al. 2016), 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 per cent). Fig. 2 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 (the necessary processing required to accurately measure astrometry for such sources was not established in DR1; Lindegren et al. 2016). Additionally, sources with extremely red or blue colours, or those with high proper-motions (>3.5 arcsec yr−1) are not included in DR1 (Gaia Collaboration 2016a; Arenou et al. 2017), though this will not include any of our targets. Where Gaia data are unavailable, we take parallaxes and proper motions from the revised Hipparcos catalogue (van Leeuwen 2007). Figure 2. View largeDownload slide Colour-magnitude diagrams for our Sco-Cen samples illustrating the completeness of the Gaia data (photometry from Anderson & Francis 2012). Sources with astrometry from Gaia are shown in red and those for which Hipparcos data was used are shown in blue. A 1 magnitude reddening vector is shown for illustrative purposes. The mean extinctions of members in each subgroup are AV = 0.73, 0.17, and 0.23 mag for US, UCL and LCC, respectively (de Geus, de Zeeuw & Lub 1989). Figure 2. View largeDownload slide Colour-magnitude diagrams for our Sco-Cen samples illustrating the completeness of the Gaia data (photometry from Anderson & Francis 2012). Sources with astrometry from Gaia are shown in red and those for which Hipparcos data was used are shown in blue. A 1 magnitude reddening vector is shown for illustrative purposes. The mean extinctions of members in each subgroup are AV = 0.73, 0.17, and 0.23 mag for US, UCL and LCC, respectively (de Geus, de Zeeuw & Lub 1989). The precision of the astrometric data used is shown in Fig. 3. The median proper motion precision is 0.06–0.11 mas yr−1 (≃0.04–0.08 km s−1 at a representative distance of 150 pc) in α and 0.05-0.08 mas yr−1 (≃0.04–0.06 km s−1) in δ for all the data for the three subgroups. When just Gaia data are considered, the median proper motion precision is 0.04–0.06 mas yr−1 (0.03–0.06 km s−1) in α and ∼0.04 mas yr−1 (0.03 km s−1) in δ. The median parallax uncertainty is ∼0.4 mas for all subsets of the data.2 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 (Lindegren et al. 2016). 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. Figure 3. View largeDownload slide Cumulative uncertainty distributions for the samples of stars in US (red), UCL (green), and LCC (blue) for parallax, proper motions, and RVs. Figure 3. View largeDownload slide Cumulative uncertainty distributions for the samples of stars in US (red), UCL (green), and LCC (blue) for parallax, proper motions, and RVs. 2.2 Radial velocities Results from Gaia's RV spectrometer are not included in DR1 and so we gathered RVs from the literature. We followed Murphy et al. (2015) in the prioritization of RVs from different samples, but chose not to exclude RVs from known binaries.3 Our primary source of RVs was the Magellan/MIKE study of Chen et al. (2011), from which 89 matches were found with our catalogue. Our next sources were the Pulkovo Compilation of Radial Velocities (PCRV, Gontcharov 2006), from which 118 matches were found, and CRVAD-2 (Kharchenko et al. 2007), 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 Dahm, Slesnick & White (2012). 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−1 (LCC). 3 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 4) and for assessing the evidence for linear expansion (Section 5). 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 7.1). 3.1 Individual distances Distances to individual sources and their uncertainties were calculated using the Bayesian inference method of Bailer-Jones (2015), 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) = r2 exp(−r/L), where L = 1.35 kpc is the scale length parameter, as recommended by Bailer-Jones (2015) and Astraatmadja & Bailer-Jones (2016).4 To fully sample the posterior distribution, we use the Markov Chain Monte Carlo (MCMC) ensemble sampler emcee (Foreman-Mackey et al. 2013), and use the mode of the corresponding posterior as our distance estimate and the 68 per cent 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. Fig. 4 shows the 3-dimensional spatial structure of Sco-Cen with the internal structure resolved in all three subgroups for the first time (de Zeeuw et al. 1999 were unable to resolve the internal structure of any of the subgroups, while de Bruijne 1999b 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. Figure 4. View largeDownload slide 3-dimensional spatial structure of Sco-Cen shown in the X − Y (top), X − Z (middle), and Y − Z planes (bottom panel). The sources are coloured according to their subgroup members: red for US, green for UCL, and blue for LCC. 1σ uncertainties on the positions of all sources are shown in grey. Also shown are the centres of each subgroup (large circles; see Section 3.2), the bulk motion of each subgroup over 2 Myr (arrows with black outline; see Section 4.1), and the relative motion of each subgroup over 10 Myr with the overall Sco-Cen bulk motion (U = −7.0 ± 0.2, V = −19.6 ± 0.2, W = −6.1 ± 0.1 km s−1) subtracted (arrow with grey outline). Figure 4. View largeDownload slide 3-dimensional spatial structure of Sco-Cen shown in the X − Y (top), X − Z (middle), and Y − Z planes (bottom panel). The sources are coloured according to their subgroup members: red for US, green for UCL, and blue for LCC. 1σ uncertainties on the positions of all sources are shown in grey. Also shown are the centres of each subgroup (large circles; see Section 3.2), the bulk motion of each subgroup over 2 Myr (arrows with black outline; see Section 4.1), and the relative motion of each subgroup over 10 Myr with the overall Sco-Cen bulk motion (U = −7.0 ± 0.2, V = −19.6 ± 0.2, W = −6.1 ± 0.1 km s−1) subtracted (arrow with grey outline). 3.2 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 parametrized 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, $$\boldsymbol{\theta }$$, best explain the observations, $$\boldsymbol{d}$$. In Bayes's theorem, this is known as the posterior distribution, $$\boldsymbol{P(\theta | d)}$$, and is given by   $$\boldsymbol{P(\theta | d)} = \frac{ \boldsymbol{P(d | \theta )} \times \boldsymbol{P(\theta )} }{ \boldsymbol{P(d)} }$$ (1)where $$\boldsymbol{P(d | \theta )}$$ is the likelihood model, $$\boldsymbol{P(\theta )}$$ are the priors (which include our a priori knowledge about the model parameters) and $$\boldsymbol{P(d)}$$ is a normalizing 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 (see discussion in Bailer-Jones 2015). 3.2.1 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 Elson, Fall & Freeman (EFF, 1987) radial profile, where the stellar surface density, Σ(r), at a radius, r, is given by   $$\Sigma (r) = \Sigma _0 \left( 1 + \frac{r^2}{a^2} \right)^{-\gamma / 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 \sqrt{4^{1/\gamma } - 1}$$, while the power-law index can be used to calculate the quantity η in the virial equation (e.g. Portegies Zwart, McMillan & Gieles 2010). 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 five (distance, central α and δ, a, and γ) and nine parameters (where a and γ become aX, aY, aZ, γX, γY, and γZ). We adopt liberal priors where possible, to minimize 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 per cent confidence interval used as the 1σ uncertainty. 3.2.2 Results The results of the model fits are provided in Table 1 and illustrated in Fig. 5. 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. Figure 5. View largeDownload slide Source distributions (black histogram) in parallax (top row), α and δ (bottom row) for the three subgroups of Sco-Cen compared to the EFF model fit results for 1D (red) and 3D (green) models. The model fit distributions were calculated by sampling the best-fitting forward models 1 000 000 times and plotting the resulting probability distribution functions. Figure 5. View largeDownload slide Source distributions (black histogram) in parallax (top row), α and δ (bottom row) for the three subgroups of Sco-Cen compared to the EFF model fit results for 1D (red) and 3D (green) models. The model fit distributions were calculated by sampling the best-fitting forward models 1 000 000 times and plotting the resulting probability distribution functions. Table 1. Structural and kinematic properties of the three Sco-Cen subgroups determined in this work. The structure parameters are all determined from the 1D EFF model, but these are generally in close agreement with the 3D EFF model. The distance uncertainties do not take into account the 0.3 mas systematic parallax uncertainty present in DR1. The velocity dispersion fitting results are taken from the model using a binary fraction of 100 per cent. Subgroup  US  UCL  LCC  Distance (pc)  $$143.0^{+0.3}_{-0.4}$$  $$135.9^{+0.5}_{-0.4}$$  $$115.2^{+0.3}_{-0.3}$$  α0 (deg)  $$243.02^{+0.15}_{-0.10}$$  $$228.36^{+0.27}_{-0.24}$$  $$188.57^{+0.23}_{-0.19}$$  δ0 (deg)  $$-24.19^{+0.08}_{-0.18}$$  $$-40.50^{+0.13}_{-0.13}$$  $$-54.53^{+0.20}_{-0.06}$$  X0 (pc)  $$133.1^{+0.4}_{-0.3}$$  $$114.0^{+0.5}_{-0.5}$$  $$57.7^{+0.3}_{-0.3}$$  Y0 (pc)  $$-21.3^{+0.2}_{-0.4}$$  $$-65.5^{+0.4}_{-0.5}$$  $$-98.3^{+0.3}_{-0.3}$$  Z0 (pc)  $$47.6^{+0.3}_{-0.3}$$  $$34.6^{+0.4}_{-0.3}$$  $$16.6^{+0.3}_{-0.2}$$  a (pc)  $$35.8^{+1.9}_{-1.3}$$  $$69.8^{+3.3}_{-3.2}$$  $$50.1^{+4.0}_{-3.2}$$  γ  $$14.5^{+1.4}_{-1.0}$$  $$17.6^{+1.3}_{-1.3}$$  $$15.2^{+1.9}_{-1.5}$$  rhl (pc)  $$11.4^{+0.2}_{-0.2}$$  $$20.0^{+0.3}_{-0.2}$$  $$15.5^{+0.3}_{-0.2}$$  U0 (km s−1)  $$-6.16_{-0.13}^{+0.14}$$  $$-5.90_{-0.12}^{+0.17}$$  $$-8.96_{-0.24}^{+0.14}$$  V0 (km s−1)  $$-16.89_{-0.10}^{+0.08}$$  $$-20.00_{-0.10}^{+0.12}$$  $$-20.55_{-0.13}^{+0.16}$$  W0 (km s−1)  $$-7.05_{-0.08}^{+0.09}$$  $$-5.80_{-0.09}^{+0.09}$$  $$-6.29_{-0.11}^{+0.12}$$  σU (km s−1)  $$1.63_{-0.20}^{+0.20}$$  $$1.96_{-0.12}^{+0.13}$$  $$1.89_{-0.13}^{+0.21}$$  σV (km s−1)  $$1.14_{-0.14}^{+0.13}$$  $$0.73_{-0.19}^{+0.15}$$  $$0.90_{-0.30}^{+0.50}$$  σW (km s−1)  $$2.51_{-0.09}^{+0.11}$$  $$1.27_{-0.14}^{+0.11}$$  $$0.51_{-0.18}^{+0.25}$$  σ3D (km s−1)  $$3.20_{-0.20}^{0.22}$$  $$2.45_{-0.20}^{+0.20}$$  $$2.15_{-0.24}^{+0.47}$$  αCP (deg)  $$116.22_{-9.46}^{+10.70}$$  $$100.88_{-1.52}^{+1.69}$$  $$95.32_{-1.31}^{+1.29}$$  δCP (deg)  $$-55.29_{-4.56}^{+5.37}$$  $$-35.88_{-1.61}^{+1.63}$$  $$-28.68_{-1.48}^{+1.53}$$  Subgroup  US  UCL  LCC  Distance (pc)  $$143.0^{+0.3}_{-0.4}$$  $$135.9^{+0.5}_{-0.4}$$  $$115.2^{+0.3}_{-0.3}$$  α0 (deg)  $$243.02^{+0.15}_{-0.10}$$  $$228.36^{+0.27}_{-0.24}$$  $$188.57^{+0.23}_{-0.19}$$  δ0 (deg)  $$-24.19^{+0.08}_{-0.18}$$  $$-40.50^{+0.13}_{-0.13}$$  $$-54.53^{+0.20}_{-0.06}$$  X0 (pc)  $$133.1^{+0.4}_{-0.3}$$  $$114.0^{+0.5}_{-0.5}$$  $$57.7^{+0.3}_{-0.3}$$  Y0 (pc)  $$-21.3^{+0.2}_{-0.4}$$  $$-65.5^{+0.4}_{-0.5}$$  $$-98.3^{+0.3}_{-0.3}$$  Z0 (pc)  $$47.6^{+0.3}_{-0.3}$$  $$34.6^{+0.4}_{-0.3}$$  $$16.6^{+0.3}_{-0.2}$$  a (pc)  $$35.8^{+1.9}_{-1.3}$$  $$69.8^{+3.3}_{-3.2}$$  $$50.1^{+4.0}_{-3.2}$$  γ  $$14.5^{+1.4}_{-1.0}$$  $$17.6^{+1.3}_{-1.3}$$  $$15.2^{+1.9}_{-1.5}$$  rhl (pc)  $$11.4^{+0.2}_{-0.2}$$  $$20.0^{+0.3}_{-0.2}$$  $$15.5^{+0.3}_{-0.2}$$  U0 (km s−1)  $$-6.16_{-0.13}^{+0.14}$$  $$-5.90_{-0.12}^{+0.17}$$  $$-8.96_{-0.24}^{+0.14}$$  V0 (km s−1)  $$-16.89_{-0.10}^{+0.08}$$  $$-20.00_{-0.10}^{+0.12}$$  $$-20.55_{-0.13}^{+0.16}$$  W0 (km s−1)  $$-7.05_{-0.08}^{+0.09}$$  $$-5.80_{-0.09}^{+0.09}$$  $$-6.29_{-0.11}^{+0.12}$$  σU (km s−1)  $$1.63_{-0.20}^{+0.20}$$  $$1.96_{-0.12}^{+0.13}$$  $$1.89_{-0.13}^{+0.21}$$  σV (km s−1)  $$1.14_{-0.14}^{+0.13}$$  $$0.73_{-0.19}^{+0.15}$$  $$0.90_{-0.30}^{+0.50}$$  σW (km s−1)  $$2.51_{-0.09}^{+0.11}$$  $$1.27_{-0.14}^{+0.11}$$  $$0.51_{-0.18}^{+0.25}$$  σ3D (km s−1)  $$3.20_{-0.20}^{0.22}$$  $$2.45_{-0.20}^{+0.20}$$  $$2.15_{-0.24}^{+0.47}$$  αCP (deg)  $$116.22_{-9.46}^{+10.70}$$  $$100.88_{-1.52}^{+1.69}$$  $$95.32_{-1.31}^{+1.29}$$  δCP (deg)  $$-55.29_{-4.56}^{+5.37}$$  $$-35.88_{-1.61}^{+1.63}$$  $$-28.68_{-1.48}^{+1.53}$$  View Large 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 overestimated then these uncertainties will also be overestimated). Our values are consistent with the mean distances found by de Zeeuw et al. (1999) 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 Rizzuto et al. (2011), and with slightly larger uncertainties (because of the systematic and spatially-correlated 0.3 mas parallax uncertainty in DR1 that is not reduced by averaging, Gaia Collaboration 2016a). 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 per cent 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 (e.g. Rizzuto et al. 2011; Pecaut & Mamajek 2016) 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 are 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. 4 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. 4.1 Forward modelling the velocity dispersion To measure the velocity dispersion and bulk motion of each subgroup, we use Bayesian inference, as in Section 3.2. 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. 4.1.1 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 3.1. 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 Odenkirchen et al. (2002) and Cottaar et al. (2012), 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 (Kroupa et al. 2001) 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 Maschberger (2013). We use the orbital properties of intermediate-mass stars in Sco-Cen determined by Kouwenhoven et al. (2007), who performed an extensive study of the visual, spectroscopic, and astrometric binaries in the association. They inferred a binary fraction of 100 per cent, a mass ratio distribution fq(q) ∝ q−0.4, and a semi-major axis distribution of the form fa(a) ∝ a−1 (in the range 5 ≥ a/R⊙ ≥ 5 × 106). 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 Kouwenhoven et al. (2007) 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 106 random binaries, providing us with a distribution from which we can sample. We don't consider triple systems because their properties are poorly constrained (Duchêne & Kraus 2013) 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 six free parameters, the central velocity (U0, V0, W0) and velocity dispersion (σU, σV, σW) in each dimension, for which we adopt flat and wide priors of − 100 ≤ v0[km s− 1] ≤ +100 for the central velocities and 0 ≤ σ[km s−1] ≤ 100 for the velocity dispersions. As in Section 3.2, 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 per cent confidence interval used for the 1σ uncertainty. 4.1.2 Results The results are provided in Table 1 and illustrated in Fig. 6, with the central velocities, or bulk motions, also shown in Fig. 4. The bulk motions are different from those presented by de Zeeuw et al. (1999), particularly in U and V, implying that the difference is probably due to the increased availability and quality of RVs for these stars (though the change in Sco-Cen membership by Rizzuto et al. 2011, will have had a small effect). Our bulk motions are in better agreement with those presented by Chen et al. (2011), with most agreeing within 1σ–2σ, and the remaining differences probably due to the change in membership. Figure 6. View largeDownload slide Source distributions (black histograms) for proper motions (left and centre) and RVs (right) for the three subgroups, compared to the best-fitting velocity dispersion model fit results for binary fractions of 100 per cent (green) and 70 per cent (red). The model fit distributions were calculated by sampling the best-fitting forward models 1 000 000 times and plotting the resulting probability distribution functions. Figure 6. View largeDownload slide Source distributions (black histograms) for proper motions (left and centre) and RVs (right) for the three subgroups, compared to the best-fitting velocity dispersion model fit results for binary fractions of 100 per cent (green) and 70 per cent (red). The model fit distributions were calculated by sampling the best-fitting forward models 1 000 000 times and plotting the resulting probability distribution functions. The best-fitting velocity dispersions vary from 0.5 to 2.5 km s−1 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. de Bruijne (1999b) found that the 1D velocity dispersions for all the subgroups are ≤1.0–1.5 km s−1, while Madsen et al. (2002) measured velocity dispersions of 1.1–1.3 km s−1 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−1 (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-Cen5 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). Kouwenhoven et al. (2007) find that the observations of Sco-Cen are best fit with a binary fraction of 100 per cent, though slightly smaller fractions are possible. We explored using binary fractions of 70–100 per cent 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 per cent along that axis for each 10 per cent by which the binary fraction is reduced. Fig. 6 shows the model fit results for a 70 per cent binary fraction, with little difference in the fit quality compared to that for 100 per cent binaries. 4.2 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 5.3). In this section we measure the convergent points of the three subgroups using both the refurbished Jones (1971) method (de Bruijne 1999a) and directly from the 3D bulk motions. 4.2.1 Method Jones's method iterates over a hemispherical grid of points to find the convergent point coordinates that minimizes the proper motions perpendicular to the great circle joining each star to the convergent point. The best-fitting convergent point is that which minimizes the quantity $$\sum t^2_\perp$$, 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 de Bruijne (1999a), who changed this quantity to $$t_\perp = \mu _\perp / \sqrt{\sigma ^2_\perp + \sigma ^2_v}$$, recognizing 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_\perp ^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_\perp ^2$$ will exceed the observed value of $$t_\perp ^2$$ by chance even for a correct model. Therefore, if the value of $$t_\perp ^2$$ is too high, the probability ε will be below the threshold value, which de Bruijne (1999a) recommends setting to a value of 0.954 (giving an ∼5 per cent 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_\perp ^2$$ value. 4.2.2 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 4.1.6 For US a convergent point for the entire subgroup was found with a χ2 probability of 53 per cent, rising to 95.4 per cent 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−1 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−1. Fig. 7 shows this convergent point compared to those determined previously by de Zeeuw et al. (1999) and Madsen et al. (2002) 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σ. Figure 7. View largeDownload slide Convergent points for the three subgroups of Sco-Cen determined using the Jones (1971) method (red circle and contours showing 1σ, 2σ, and 3σ uncertainties) and from the bulk motions (green circles with 1σ error bars typically equal to or smaller than the size of the symbols). Yellow circles show previously published convergent points from de Zeeuw et al. (1999, DZ99) and Madsen, Dravins & Lindegren (2002, M02), with 1σ error bars shown for the latter (uncertainties for the former are expected to be several degrees along the great circle connecting the convergent point with the subgroup). Dashed lines with black markers show the predicted convergent points derived from the mean space motions at a range of expansion ages. Figure 7. View largeDownload slide Convergent points for the three subgroups of Sco-Cen determined using the Jones (1971) method (red circle and contours showing 1σ, 2σ, and 3σ uncertainties) and from the bulk motions (green circles with 1σ error bars typically equal to or smaller than the size of the symbols). Yellow circles show previously published convergent points from de Zeeuw et al. (1999, DZ99) and Madsen, Dravins & Lindegren (2002, M02), with 1σ error bars shown for the latter (uncertainties for the former are expected to be several degrees along the great circle connecting the convergent point with the subgroup). Dashed lines with black markers show the predicted convergent points derived from the mean space motions at a range of expansion ages. For UCL a convergent point for the entire subgroup was found with a χ2 probability of 98 per cent 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 per cent, rising to 95.4 per cent after removing the 10 stars with the highest |t⊥| values. A higher velocity dispersion of 1.5 km s−1 would result in an immediate fit with a χ2 probability of >99 per cent 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−1 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. 5 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 (Blaauw 1964b). 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. 5.1 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 Blaauw (1964b) linear expansion model individual RVs, vrad, are predicted to follow the relation   $$v_\mathrm{rad} = v^\prime \, \mathrm{cos} \, \lambda ^\prime + \kappa d + K,$$ (3)where λ΄ is the angular separation between star and convergent point (calculated using just the proper motions; see Section 4.2), 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^\prime = \frac{ \mu _v A }{ \varpi \, \mathrm{sin} \, \lambda ^\prime },$$ (4)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). Fig. 8 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 (see Hogg, Bovy & Lang 2010) and then marginalizing over this parameter to obtain the uncertainties on the resulting fit parameters. Figure 8. View largeDownload slide Distance to sources in each subgroup of Sco-Cen plotted against the difference between the observed RV and that predicted from the moving group method (v΄ cos λ). The error bars for individual sources show the 16–84 per cent confidence intervals, calculated from Monte Carlo simulations. The solid lines show the best-fitting linear relationship between the plotted quantities, with 1σ uncertainties shown with dashed lines. The best-fitting slopes, κ, and uncertainties are noted in each panel. If the subgroups are not expanding, we would expect a slope of zero, while if they expanding, we would expect a positive slope. 3σ outliers on the ordinate, shown with non-coloured symbols, were excluded from the fit. Figure 8. View largeDownload slide Distance to sources in each subgroup of Sco-Cen plotted against the difference between the observed RV and that predicted from the moving group method (v΄ cos λ). The error bars for individual sources show the 16–84 per cent confidence intervals, calculated from Monte Carlo simulations. The solid lines show the best-fitting linear relationship between the plotted quantities, with 1σ uncertainties shown with dashed lines. The best-fitting slopes, κ, and uncertainties are noted in each panel. If the subgroups are not expanding, we would expect a slope of zero, while if they expanding, we would expect a positive slope. 3σ outliers on the ordinate, shown with non-coloured symbols, were excluded from the fit. 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 Pecaut et al. (2012). 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 Pecaut & Mamajek (2016) 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. 5.2 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. Fig. 9 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. Figure 9. View largeDownload slide Positions versus velocities along the three Galactic Cartesian axes XYZ for each of the three subgroups of Sco-Cen (US on the left, UCL in the centre, and LCC on the right). 1σ error bars are shown for all sources. The solid lines show the best-fit linear relationships between the plotted quantities, with 1σ uncertainties shown with dashed lines. The best-fitting slopes, κ, and uncertainties are noted in each panel. If the subgroups are not expanding we would expect slopes of zero, while if they are expanding we would expect a positive slope (a negative slope implies contraction along that axis). 3σ outliers on the ordinate, shown with non-coloured symbols, were excluded from the fits. Figure 9. View largeDownload slide Positions versus velocities along the three Galactic Cartesian axes XYZ for each of the three subgroups of Sco-Cen (US on the left, UCL in the centre, and LCC on the right). 1σ error bars are shown for all sources. The solid lines show the best-fit linear relationships between the plotted quantities, with 1σ uncertainties shown with dashed lines. The best-fitting slopes, κ, and uncertainties are noted in each panel. If the subgroups are not expanding we would expect slopes of zero, while if they are expanding we would expect a positive slope (a negative slope implies contraction along that axis). 3σ outliers on the ordinate, shown with non-coloured symbols, were excluded from the fits. 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 significance 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. 5.3 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 3) and compare these with the proper motion convergent points. Fig. 7 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 (Pecaut & Mamajek 2016). 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 Pecaut & Mamajek (2016). 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. 5.4 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) = X0 + γUt, 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. Fig. 10 shows the 1σ dispersions (16th–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. Figure 10. View largeDownload slide 1σ dispersions (16th–84th percentiles used to minimize the impact of outliers) of the size of each subgroup in each axis of the XYZ coordinate system (σX shown with a dashed line, σY with a dotted line, and σZ with a dash-dotted line) as well as the quadrature sum of all three dimensions (solid line). Black lines show the results of traceback analysis performed with linear trajectories and the red lines show the analysis performed using the epicycle orbit approximation. Shaded areas around the 3D sum lines show the 90 per cent confidence interval in the quadrature sum dispersion of each subgroup determined from Monte Carlo simulations exploring the impact of uncertainties in the velocities. Figure 10. View largeDownload slide 1σ dispersions (16th–84th percentiles used to minimize the impact of outliers) of the size of each subgroup in each axis of the XYZ coordinate system (σX shown with a dashed line, σY with a dotted line, and σZ with a dash-dotted line) as well as the quadrature sum of all three dimensions (solid line). Black lines show the results of traceback analysis performed with linear trajectories and the red lines show the analysis performed using the epicycle orbit approximation. Shaded areas around the 3D sum lines show the 90 per cent confidence interval in the quadrature sum dispersion of each subgroup determined from Monte Carlo simulations exploring the impact of uncertainties in the velocities. For a more accurate traceback of the stellar motions, we use the epicycle approximation, employing the orbital equations from Fuchs et al. (2006). We use the Oort A and B constants from Feast & Whitelock (1997), the local disc density from Holmberg & Flynn (2004), the local standard of rest velocity from Schönrich, Binney & Dehnen (2010) and a solar Z distance above the Galactic plane of 17 pc (Karim & Mamajek 2017). Fig. 10 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 Fig. 10 we show the 90 per cent 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 Brown et al. (1997) 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 Brown et al. (1997) however. 5.5 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 (referred to as virtual expansion /contraction; Brown et al. 1997). 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 Brown et al. (1997), using the median RV of stars in each subgroup (to minimize the impact of binarity). The corrected proper motion vector maps for the three subgroups are shown in Fig. 11. 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 6, that these patterns are due to the kinematic substructure within the subgroups. Figure 11. View largeDownload slide Proper motion vector maps for the three subgroups of Sco-Cen projected on to inverted H α images from Finkbeiner (2003). The proper motions have been corrected for radial streaming motion according to the equations in Brown et al. (1997) and using bulk (median) RVs of −6.2, 2.8, and 13.0 km s−1 for the three subgroups. Points show the current positions of stars and vectors show the proper motions over 0.5 Myr with the bulk motion of each subgroup subtracted to show the motion in the reference frame of each subgroup. The vectors have been colour-coded based on their direction of motion (see colour wheel in lower-right corner) to highlight kinematic substructure. Note that because the kinematic reference frame shown for each subgroup is different, the colours of the proper motions in each subgroup do not represent the same absolute proper motion as stars in other subgroups. Figure 11. View largeDownload slide Proper motion vector maps for the three subgroups of Sco-Cen projected on to inverted H α images from Finkbeiner (2003). The proper motions have been corrected for radial streaming motion according to the equations in Brown et al. (1997) and using bulk (median) RVs of −6.2, 2.8, and 13.0 km s−1 for the three subgroups. Points show the current positions of stars and vectors show the proper motions over 0.5 Myr with the bulk motion of each subgroup subtracted to show the motion in the reference frame of each subgroup. The vectors have been colour-coded based on their direction of motion (see colour wheel in lower-right corner) to highlight kinematic substructure. Note that because the kinematic reference frame shown for each subgroup is different, the colours of the proper motions in each subgroup do not represent the same absolute proper motion as stars in other subgroups. To quantify the expansion of each subgroup, we follow the method of Wright et al. (2016) and calculate, for each star, the amount of kinetic energy in the radial and transverse directions (using the subgroup centres calculated in Section 3.2), 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 per cent (US), 66 per cent (UCL), and 45 per cent (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 per cent (US), 67 per cent (UCL), and 90 per cent (LCC) of the kinetic energy in the direction of expansion. This suggests that, while there is not evidence for coherent expansion patterns in the correct proper motions of stars in each subgroup, they do appear to be expanding. 6 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 (e.g. the Gamma Velorum cluster shows two distinct kinematic components; Jeffries et al. 2014) or as smaller, more numerous substructures (as observed in Cygnus OB2, Wright et al. 2016). 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 (e.g. Larson 1981; Wright et al. 2016). In this section we explore the evidence for kinematic substructure in Sco-Cen. Fig. 12 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 (as per Wright et al. 2016). 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. Figure 12. View largeDownload slide Proper motion vector map for the entire Sco-Cen association projected on to an inverted H α image from Finkbeiner (2003). The proper motions were corrected for radial streaming motion according to the equations in Brown et al. (1997) and using bulk (median) RVs of −6.2, 2.8, and 13.0 km s−1 for the three subgroups. Points show the current positions of stars and vectors show the proper motions over 0.5 Myr with the bulk (median) proper motion of the entire association subtracted to show the relative motion. The vectors have been colour-coded based on their direction of motion (see colour wheel in lower left) to highlight kinematic substructure. Figure 12. View largeDownload slide Proper motion vector map for the entire Sco-Cen association projected on to an inverted H α image from Finkbeiner (2003). The proper motions were corrected for radial streaming motion according to the equations in Brown et al. (1997) and using bulk (median) RVs of −6.2, 2.8, and 13.0 km s−1 for the three subgroups. Points show the current positions of stars and vectors show the proper motions over 0.5 Myr with the bulk (median) proper motion of the entire association subtracted to show the relative motion. The vectors have been colour-coded based on their direction of motion (see colour wheel in lower left) to highlight kinematic substructure. 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 (e.g. David & Hillenbrand 2015) and therefore suggest that these are field stars that should be removed from the membership of Sco-Cen. Fig. 11 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 5.2. 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. 6.1 Quantifying kinematic substructure 6.1.1 Method In an attempt to quantify this substructure we follow Wright et al. (2016) and use Moran's spatial correlation index I (Moran 1950), 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 I0 = −1/(n − 1), with I > I0, indicating positive spatial correlation and I < I0, indicating negative spatial correlation. We use this instead of Geary's C index (Geary 1954), since Moran's I statistic is a global measure that can identify and quantify large-scale spatial correlation (as appears to be evident in Fig. 11), 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 Wright et al. (2016). 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 on to 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 4.1). We then calculated Moran's I index for the simulated data set 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. 6.1.2 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 Fig. 11 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 (Wright et al. 2016). 7 DISCUSSION 7.1 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   $$\sigma _{3{\rm D}}^2 = \frac{G M_{{\rm vir}}}{2 r_{{\rm vir}}}$$ (5)where σ3D is the 3D velocity dispersion, G is the gravitational constant, Mvir is the virial mass, and rvir is the virial radius. We substitute the parameter η = 6rvir/reff (e.g. Portegies Zwart et al. 2010), where reff is the effective (or half-light) radius, giving   $$M_{{\rm vir}} = \eta \, \frac{\sigma _{3D}^2 \, r_{{\rm eff}}}{3G}$$ (6)where the parameters η and reff can be calculated from the Elson et al. (1987) parameters γ and a (see Section 3.2). For the large γ values calculated for the subgroups we find η values of 9.1 (US), 9.0 (UCL), and 9.1 (LCC),7 while the effective radii, reff 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 1), and η values, we calculate virial masses of $$8.2^{+1.5}_{-1.2} \times 10^4$$ (US), $$8.3^{+1.8}_{-1.5} \times 10^4$$ (UCL), and $$5.1^{+2.9}_{-1.3} \times 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. Preibisch & Mamajek (2008) estimate US to have a total mass of ∼2000 M⊙, while Mamajek et al. (2002) 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⊙ (using the same initial mass function and binary fraction as Preibisch & Mamajek 2008). 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 de Geus (1992) estimates to have a mass of 8 × 104 M⊙. The author argues 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 × 105 and 1 × 105 M⊙, respectively, are more massive than their virial masses, implying that the subgroups might have initially been in a subvirial state. 7.2 Expansion of the Sco-Cen subgroups Ever since Ambartsumian (1947) 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 (e.g. Blaauw 1964a, 1991; Brown et al. 1999). These ideas lead to suggestions that OB associations may be the expanded remnants of compact star clusters (Lada & Lada 2003), disrupted by processes such as residual gas expulsion (e.g. Hills 1980; Lada et al. 1984). In Section 5 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 5.1 and 5.2), while others better account for the Galactic potential (Section 5.4) or provide an overall view of the stellar motions (Section 5.5). 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 5.1 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 5.2) 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, Mamajek & Bell (2014) 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 time-scales involved are probably too short (Dobbs & Pringle 2013, estimate a shear time-scale of ∼70 Myr in the vicinity of the Sun). 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 5.4), 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 Fig. 11 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 (Mamajek & Bell 2014) and Tuc-Hor (Makarov 2007), 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 (e.g. Pfalzner 2009). This is particularly clear in the case of the three subgroups of Sco-Cen, which Pfalzner & Kaczmarek (2013) 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 Fig. 11 suggest a more complex expansion pattern consisting of multiple expanding substructures within each subgroup. 7.3 Kinematic substructure and the formation of Sco-Cen Molecular clouds are known to be highly substructured, both spatially (Falgarone, Phillips & Walker 1991) and kinematically (Hacar et al. 2013), with this hierarchical or fractal structure passed on to the formed stars (Elmegreen & Elmegreen 2001) and often evident in their initial spatial distribution (Gutermuth et al. 2008). Dynamical interactions between groups of stars can erase this substructure (e.g. Scally & Clarke 2002; Parker et al. 2014), a process that acts on the dynamical time-scale of a region. If the dynamical time-scale is longer than the age of the region then the primordial substructure can be preserved (e.g. Goodwin & Whitworth 2004), though dynamical time-scales can change as a region expands or collapses (so an instantaneous measure of the dynamical time-scale is not 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 time-scales 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 (e.g. Parker et al. 2014). The three subgroups of Sco-Cen each display evidence for kinematic substructure, as has been observed in other OB associations (Wright et al. 2014) and across large star-forming complexes (e.g. Fűrész et al. 2006, 2008; Tobin et al. 2015). 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 do not 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 Pecaut & Mamajek (2016), 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 Blaauw (1964a) and extended by Preibisch & Zinnecker (1999), the impact of a high-velocity cloud on the disc (Lepine & Duvert 1994), to its formation in a large ring-like structure containing other young groups of stars known as the Gould Belt (Lindblad et al. 1973). 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 (e.g. Sartori, Lépine & Dias 2003; Preibisch & Mamajek 2008). 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 does not necessarily disprove it. 8 CONCLUSIONS Our kinematic study of the Sco-Cen OB association using Gaia DR1 parallaxes and proper motions has led 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 that 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−1 (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 Blaauw (1964b) 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 (Wright et al. 2016), 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 Pecaut & Mamajek (2016) 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 (e.g. Makarov 2007; Mamajek & Bell 2014; Wright et al. 2016) 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 (Pfalzner 2009) and that star cluster disruption mechanisms, such as residual gas expulsion (Kroupa et al. 2001), 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 discs. Acknowledgements NJW acknowledges an Science and Technology Facilities Council 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 European Space Agency 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 Centre de Données astronomiques de Strasbourg, Strasbourg, France), Astropy (Astropy Collaboration 2013) and TOPCAT (Tool for OPerations on Catalogues And Tables, Taylor 2005). NOTE ADDED IN PROOF At the time this paper was accepted two papers came to our attention that, respectively, find evidence for expansion of the OB associations Per OB1 and Car OB1 (Mel'nik & Dambis 2017) and no evidence for OB association expansion (Ward & Kruijssen 2018), both based on Gaia proper motions. Footnotes 1 Murphy, Lawson & Bento (2015) note that some of the RVs used by Rizzuto et al. (2011) 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. 2 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. 3 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. 4 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 3.2 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. 5 We repeated these fits, using the de Zeeuw et al. (1999) membership list, and find that the change in sample to Rizzuto et al. (2011) has caused the velocity dispersion to increase by ∼10 per cent. 6 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. 7 The dependence of η on γ is very weak at large γ values and therefore the uncertainties on η are almost zero. REFERENCES Adams F. C., Proszkow E. M., Fatuzzo M., Myers P. C., 2006, ApJ , 641, 504 https://doi.org/10.1086/500393 CrossRef Search ADS   Ambartsumian V. A., 1947, Stellar Evolution and Astrophysics . Armenian Acad. Sci. Anderson E., Francis C., 2012, Astron. Lett. , 38, 331 https://doi.org/10.1134/S1063773712050015 CrossRef Search ADS   Arenou F. et al.  , 2017, A&A , 599, A50 CrossRef Search ADS   Astraatmadja T. L., Bailer-Jones C. A. L., 2016, ApJ , 832, 137 https://doi.org/10.3847/0004-637X/832/2/137 CrossRef Search ADS   Astropy Collaboration, 2013, A&A , 558, A33 CrossRef Search ADS   Bailer-Jones C. A. L., 2015, PASP , 127, 994 https://doi.org/10.1086/683116 CrossRef Search ADS   Baumgardt H., Kroupa P., 2007, MNRAS , 380, 1589 https://doi.org/10.1111/j.1365-2966.2007.12209.x CrossRef Search ADS   Blaauw A., 1946, Publi. Kapteyn Astron. Lab. Groningen , 52, 1 Blaauw A., 1964a, ARA&A , 2, 213 CrossRef Search ADS   Blaauw A., 1964b, in Kerr F. J., ed., IAU Symp., Vol. 20, The Galaxy and the Magellanic Clouds . Australian Academy of Science, Canberra, p. 50 Blaauw A., 1991, in Lada C. J. Kylafis N. D., eds., NATO ASIC Proc. 342: The Physics of Star Formation and Early Stellar Evolution . Kluwer, Dordrecht, p. 125 Google Scholar CrossRef Search ADS   Breitschwerdt D., Feige J., Schulreich M. M., Avillez M. A. D., Dettbarn C., Fuchs B., 2016, Nature , 532, 73 https://doi.org/10.1038/nature17424 CrossRef Search ADS PubMed  Bressert E. et al.  , 2010, MNRAS , 409, L54 https://doi.org/10.1111/j.1745-3933.2010.00946.x CrossRef Search ADS   Brown A. G. A., Dekker G., de Zeeuw P. T., 1997, MNRAS , 285, 479 https://doi.org/10.1093/mnras/285.3.479 CrossRef Search ADS   Brown A. G. A. Blaauw A. Hoogerwerf R. de Bruijne J. H. J. de Zeeuw P. T., 1999, in Series C, Lada C. J., Kylafis N. D., eds., NATO Advanced Science Institutes (ASI) Vol. 540, p. 411 Carpenter J. M., Mamajek E. E., Hillenbrand L. A., Meyer M. R., 2006, ApJ , 651, L49 https://doi.org/10.1086/509121 CrossRef Search ADS   Chen C. H., Mamajek E. E., Bitner M. A., Pecaut M., Su K. Y. L., Weinberger A. J., 2011, ApJ , 738, 122 https://doi.org/10.1088/0004-637X/738/2/122 CrossRef Search ADS   Cottaar M., Meyer M. R., Andersen M., Espinoza P., 2012, A&A , 539, A5 CrossRef Search ADS   Dahm S. E., Slesnick C. L., White R. J., 2012, ApJ , 745, 56 https://doi.org/10.1088/0004-637X/745/1/56 CrossRef Search ADS   David T. J., Hillenbrand L. A., 2015, ApJ , 804, 146 https://doi.org/10.1088/0004-637X/804/2/146 CrossRef Search ADS   de Bruijne J. H. J., 1999a, MNRAS , 306, 381 https://doi.org/10.1046/j.1365-8711.1999.02643.x CrossRef Search ADS   de Bruijne J. H. J., 1999b, MNRAS , 310, 585 https://doi.org/10.1046/j.1365-8711.1999.02953.x CrossRef Search ADS   de Geus E. J., 1992, A&A , 262, 258 de Geus E. J., de Zeeuw P. T., Lub J., 1989, A&A , 216, 44 de Zeeuw P. T., Hoogerwerf R., de Bruijne J. H. J., Brown A. G. A., Blaauw A., 1999, AJ , 117, 354 https://doi.org/10.1086/300682 CrossRef Search ADS   Dobbs C. L., Pringle J. E., 2013, MNRAS , 432, 653 https://doi.org/10.1093/mnras/stt508 CrossRef Search ADS   Duchêne G., Kraus A., 2013, ARA&A , 51, 269 CrossRef Search ADS   Elmegreen B. G. Efremov Y. Pudritz R. E. Zinnecker H., 2000, Protostars Planets IV , 179 Elmegreen B. G., Elmegreen D. M., 2001, AJ , 121, 1507 https://doi.org/10.1086/319416 CrossRef Search ADS   Elson R. A. W., Fall S. M., Freeman K. C., 1987, ApJ , 323, 54 https://doi.org/10.1086/165807 CrossRef Search ADS   ESA ed., 1997, ESA Special Publication, Vol. 1200, The HIPPARCOS and TYCHO Catalogues. Astrometric and Photometric Star Catalogues Derived from the ESA HIPPARCOS Space Astrometry Mission . ESA. Available at: https://www.cosmos.esa.int/documents/532822/552851/vol1_all.pdf/99adf6e3-6893-4824-8fc2-8d3c9cbba2b5 Falgarone E., Phillips T. G., Walker C. K., 1991, ApJ , 378, 186 https://doi.org/10.1086/170419 CrossRef Search ADS   Feast M., Whitelock P., 1997, MNRAS , 291, 683 https://doi.org/10.1093/mnras/291.4.683 CrossRef Search ADS   Finkbeiner D. P., 2003, ApJS , 146, 407 https://doi.org/10.1086/374411 CrossRef Search ADS   Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP , 125, 306 https://doi.org/10.1086/670067 CrossRef Search ADS   Fuchs B., Breitschwerdt D., de Avillez M. A., Dettbarn C., Flynn C., 2006, MNRAS , 373, 993 https://doi.org/10.1111/j.1365-2966.2006.11044.x CrossRef Search ADS   Fűrész G. et al.  , 2006, ApJ , 648, 1090, https://doi.org/10.1086/506140 CrossRef Search ADS   Fűrész G., Hartmann L. W., Megeath S. T., Szentgyorgyi A. H., Hamden E. T., 2008, ApJ , 676, 1109 https://doi.org/10.1086/525844 CrossRef Search ADS   Gaia Collaboration, 2016a, A&A , 595, A2 CrossRef Search ADS   Gaia Collaboration, 2016b, A&A , 595, A1 CrossRef Search ADS   Geary R. C., 1954, Incorp. Stat. , 5, 115 https://doi.org/10.2307/2986645 Gontcharov G. A., 2006, Astron. Lett. , 32, 759 https://doi.org/10.1134/S1063773706110065 CrossRef Search ADS   Goodwin S. P., Whitworth A. P., 2004, A&A , 413, 929 CrossRef Search ADS   Guarcello M. G. et al.  , 2016, preprint (arXiv:1605.01773) Gutermuth R. A. et al.  , 2008, ApJ , 674, 336 https://doi.org/10.1086/524722 CrossRef Search ADS   Hacar A., Tafalla M., Kauffmann J., Kovács A., 2013, A&A , 554, A55 CrossRef Search ADS   Hills J. G., 1980, ApJ , 235, 986 https://doi.org/10.1086/157703 CrossRef Search ADS   Hogg D. W. Bovy J. Lang D., 2010, preprint (arXiv:1008.4686) Holmberg J., Flynn C., 2004, MNRAS , 352, 440 https://doi.org/10.1111/j.1365-2966.2004.07931.x CrossRef Search ADS   Hoogerwerf R., de Bruijne J. H. J., de Zeeuw P. T., 2000, ApJ , 544, L133 https://doi.org/10.1086/317315 CrossRef Search ADS   Høg E. et al.  , 2000, A&A , 355, L27 Janson M., Lafrenière D., Jayawardhana R., Bonavita M., Girard J. H., Brandeker A., Gizis J. E., 2013, ApJ , 773, 170 https://doi.org/10.1088/0004-637X/773/2/170 CrossRef Search ADS   Jeffries R. D. et al.  , 2014, A&A , 563, A94 CrossRef Search ADS   Jones D. H. P., 1971, MNRAS , 152, 231 https://doi.org/10.1093/mnras/152.2.231 CrossRef Search ADS   Karim M. T., Mamajek E. E., 2017, MNRAS , 465, 472 https://doi.org/10.1093/mnras/stw2772 CrossRef Search ADS   Kharchenko N. V., Scholz R.-D., Piskunov A. E., Röser S., Schilbach E., 2007, Astron. Nachr. , 328, 889 https://doi.org/10.1002/asna.200710776 CrossRef Search ADS   Kouwenhoven M. B. N., Brown A. G. A., Portegies Zwart S. F., Kaper L., 2007, A&A , 474, 77 CrossRef Search ADS   Kraus A. L., Cody A. M., Covey K. R., Rizzuto A. C., Mann A. W., Ireland M. J., 2015, ApJ , 807, 3 https://doi.org/10.1088/0004-637X/807/1/3 CrossRef Search ADS   Kroupa P., Petr M. G., McCaughrean M. J., 1999, New Astron. , 4, 495 https://doi.org/10.1016/S1384-1076(99)00038-X CrossRef Search ADS   Kroupa P., Aarseth S., Hurley J., 2001, MNRAS , 321, 699 https://doi.org/10.1046/j.1365-8711.2001.04050.x CrossRef Search ADS   Lada C. J., Lada E. A., 2003, ARA&A , 41, 57 CrossRef Search ADS   Lada C. J., Margulis M., Dearborn D., 1984, ApJ , 285, 141 https://doi.org/10.1086/162485 CrossRef Search ADS   Larson R. B., 1981, MNRAS , 194, 809 https://doi.org/10.1093/mnras/194.4.809 CrossRef Search ADS   Lepine J. R. D., Duvert G., 1994, A&A , 286, 60 Lindblad P. O., Grape K., Sandqvist A., Schober J., 1973, A&A , 24, 309 Lindegren L. et al.  , 2016, A&A , 595, A4 CrossRef Search ADS   Lodieu N., 2013, MNRAS , 431, 3222 https://doi.org/10.1093/mnras/stt402 CrossRef Search ADS   Madsen S., Dravins D., Lindegren L., 2002, A&A , 381, 446 CrossRef Search ADS   Makarov V. V., 2007, ApJS , 169, 105 https://doi.org/10.1086/509887 CrossRef Search ADS   Mamajek E. E., Bell C. P. M., 2014, MNRAS , 445, 2169 https://doi.org/10.1093/mnras/stu1894 CrossRef Search ADS   Mamajek E. E., Meyer M. R., Liebert J., 2002, AJ , 124, 1670 https://doi.org/10.1086/341952 CrossRef Search ADS   Marks M., Kroupa P., 2012, A&A , 543, A8 CrossRef Search ADS   Maschberger T., 2013, MNRAS , 429, 1725 https://doi.org/10.1093/mnras/sts479 CrossRef Search ADS   Mel'nik A. M., Dambis A. K., 2017, MNRAS , 472, 3887 CrossRef Search ADS   Michalik D., Lindegren L., Hobbs D., 2015, A&A , 574, A115 CrossRef Search ADS   Moran P. A. P., 1950, Biometrika , 37, 17 https://doi.org/10.1093/biomet/37.1-2.17 CrossRef Search ADS PubMed  Murphy S. J., Lawson W. A., Bento J., 2015, MNRAS , 453, 2220 O'dell C. R., Wen Z., 1994, ApJ , 436, 194 https://doi.org/10.1086/174892 CrossRef Search ADS   Odenkirchen M., Grebel E. K., Dehnen W., Rix H.-W., Cudworth K. M., 2002, AJ , 124, 1497 https://doi.org/10.1086/342287 CrossRef Search ADS   Olczak C., Pfalzner S., Eckart A., 2008, A&A , 488, 191 CrossRef Search ADS   Parker R. J., Goodwin S. P., 2012, MNRAS , 424, 272 https://doi.org/10.1111/j.1365-2966.2012.21190.x CrossRef Search ADS   Parker R. J., Quanz S. P., 2012, MNRAS , 419, 2448 https://doi.org/10.1111/j.1365-2966.2011.19911.x CrossRef Search ADS   Parker R. J., Wright N. J., Goodwin S. P., Meyer M. R., 2014, MNRAS , 438, 620 https://doi.org/10.1093/mnras/stt2231 CrossRef Search ADS   Pecaut M. J., Mamajek E. E., 2016, MNRAS , 461, 794 https://doi.org/10.1093/mnras/stw1300 CrossRef Search ADS   Pecaut M. J., Mamajek E. E., Bubar E. J., 2012, ApJ , 746, 154 https://doi.org/10.1088/0004-637X/746/2/154 CrossRef Search ADS   Pfalzner S., 2009, A&A , 498, L37 CrossRef Search ADS   Pfalzner S., Kaczmarek T., 2013, A&A , 559, A38 CrossRef Search ADS   Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A , 48, 431 CrossRef Search ADS   Preibisch T. Mamajek E., 2008, The Nearest OB Association: Scorpius-Centaurus (Sco OB2), Handbook of Star Forming Regions , p. 235 Preibisch T., Zinnecker H., 1999, AJ , 117, 2381 https://doi.org/10.1086/300842 CrossRef Search ADS   Rizzuto A. C., Ireland M. J., Robertson J. G., 2011, MNRAS , 416, 3108 https://doi.org/10.1111/j.1365-2966.2011.19256.x CrossRef Search ADS   Rosotti G. P., Dale J. E., de Juan Ovelar M., Hubber D. A., Kruijssen J. M. D., Ercolano B., Walch S., 2014, MNRAS , 441, 2094 https://doi.org/10.1093/mnras/stu679 CrossRef Search ADS   Sartori M. J., Lépine J. R. D., Dias W. S., 2003, A&A , 404, 913 CrossRef Search ADS   Scally A., Clarke C., 2001, MNRAS , 325, 449 https://doi.org/10.1046/j.1365-8711.2001.04274.x CrossRef Search ADS   Scally A., Clarke C., 2002, MNRAS , 334, 156 https://doi.org/10.1046/j.1365-8711.2002.05503.x CrossRef Search ADS   Schönrich R., Binney J., Dehnen W., 2010, MNRAS , 403, 1829 https://doi.org/10.1111/j.1365-2966.2010.16253.x CrossRef Search ADS   Taylor M. B., 2005, in Shopbell P. Britton M. Ebert R., eds., ASP Conf. Ser., Vol. 347, Astronomical Data Analysis Software and Systems XIV , p. 29 Tobin J. J., Hartmann L., Fűrész G., Hsu W.-H., Mateo M., 2015, AJ , 149, 119 https://doi.org/10.1088/0004-6256/149/4/119 CrossRef Search ADS   Tutukov A. V., 1978, A&A , 70, 57 van Leeuwen F., 2007, A&A , 474, 653 CrossRef Search ADS   Ward J. L., Kruijssen J. M. D., 2018, MNRAS , preprint (arXiv:1801.03938) Wright N. J., Drake J. J., Drew J. E., Guarcello M. G., Gutermuth R. A., Hora J. L., Kraemer K. E., 2012, ApJ , 746, L21 https://doi.org/10.1088/2041-8205/746/2/L21 CrossRef Search ADS   Wright N. J., Parker R. J., Goodwin S. P., Drake J. J., 2014, MNRAS , 438, 639 https://doi.org/10.1093/mnras/stt2232 CrossRef Search ADS   Wright N. J., Bouy H., Drew J. E., Sarro L. M., Bertin E., Cuillandre J.-C., Barrado D., 2016, MNRAS , 460, 2593 https://doi.org/10.1093/mnras/stw1148 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# The kinematics of the Scorpius-Centaurus OB association from Gaia DR1

, Volume 476 (1) – May 1, 2018
18 pages

/lp/ou_press/the-kinematics-of-the-scorpius-centaurus-ob-association-from-gaia-dr1-0q2O60nMDw
Publisher
The Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty207
Publisher site
See Article on Publisher Site

### Abstract

Abstract 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 these data to confirm that the three subgroups of Sco-Cen are gravitationally unbound and have non-isotropic velocity dispersions, suggesting that 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 to have been born highly substructured, with multiple small-scale star formation events contributing to the overall OB association, and not as single, monolithic burst 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 1 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 discs (e.g. O'dell & Wen 1994; Wright et al. 2012; Guarcello et al. 2016), while close encounters in dense clusters can affect the properties of binary and multiple systems (Kroupa, Petr & McCaughrean 1999; Marks & Kroupa 2012; Parker & Goodwin 2012), protoplanetary discs (Scally & Clarke 2001; Olczak, Pfalzner & Eckart 2008; Rosotti et al. 2014), and influence the formation of planetary systems (Adams et al. 2006; Parker & Quanz 2012). Stars born in low-density environments (e.g. Bressert et al. 2010; Wright et al. 2014) 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 per cent of stars are found in bound clusters (Lada & Lada 2003). 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 (which can account for well over half the mass of the system; Elmegreen et al. 2000) then the cluster becomes supervirial and will expand and disperse (e.g. Tutukov 1978; Hills 1980; Lada, Margulis & Dearborn 1984; Baumgardt & Kroupa 2007). During the expansion the system would briefly be visible as a low-density group of young stars known as an association (e.g. Ambartsumian 1947; Blaauw 1964a; Brown, Dekker & de Zeeuw 1997; Kroupa, Aarseth & Hurley 2001), before dispersing into the Galactic field. This theoretical framework has existed for many decades, and while many observations support aspects of this model (e.g. Lada & Lada 2003, and references therein) it has been difficult to test kinematically due to the lack of high-precision proper motions necessary to assess evidence for expansion. Wright et al. (2016) 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 (Wright et al. 2014) 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 (Gaia Collaboration 2016b) 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 (de Zeeuw et al. 1999) and being at least 100 pc across, the association spans almost 90 deg on the sky (see Fig. 1). The association consists of three subgroups (first defined by Blaauw 1946) 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 (Pecaut & Mamajek 2016). Because of its age, the region has likely already seen several supernovae (e.g. Breitschwerdt et al. 2016), one of which may have been responsible for the ejected runaway star ζ Oph (de Geus 1992; Hoogerwerf, de Bruijne & de Zeeuw 2000). The association has been the focus of numerous studies of its low-mass population (e.g. Preibisch & Zinnecker 1999; Mamajek, Meyer & Liebert 2002), with studies focusing on the binary properties (e.g. Kouwenhoven et al. 2007; Janson et al. 2013), circumstellar discs (e.g. Carpenter et al. 2006; Chen et al. 2011), stellar ages (e.g. Pecaut, Mamajek & Bubar 2012; Pecaut & Mamajek 2016), eclipsing binaries (e.g. Kraus et al. 2015), and low-mass stars and brown dwarfs (e.g. Lodieu 2013). Its high-mass population has been studied by numerous authors, including comprehensive studies by Blaauw (1946, 1964b) and later, using Hipparcos astrometry, by de Zeeuw et al. (1999), de Bruijne (1999b), and Rizzuto et al. (2011). 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. 2 OBSERVATIONAL DATA The sample of Sco-Cen members used for our kinematic study is based on the Bayesian membership analysis of Rizzuto et al. (2011), who used the updated Hipparcos catalogue (van Leeuwen 2007) to revise the kinematic selection of de Zeeuw et al. (1999). 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 (CRVAD-2; Kharchenko et al. 2007)1 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 Rizzuto et al. (2011) were rejected based on their parallax or RVs, two parameters that were not used by de Zeeuw et al. (1999) 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. Chen et al. (2011) 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 Mamajek et al. (2002), Chen et al. (2011) and Pecaut et al. (2012) that were included in the original sample of de Zeeuw et al. (1999) were not selected by Rizzuto et al. (2011) as members, reinforcing the validity of their membership selection. Pecaut et al. (2012) recommend excluding a number of stars based on their kinematics, one of which is retained in the Rizzuto et al. (2011) catalogue, HIP 56227. We choose to retain this star for the time being and minimize 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 Fig. 1. For the handful of stars that cannot be assigned to a subgroup based on the division of de Zeeuw et al. (1999), we put them in the nearest subgroup on the plane of the sky. While Rizzuto et al. (2011) 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. Figure 1. View largeDownload slide Spatial distribution of Sco-Cen members (Rizzuto, Ireland & Robertson 2011) projected on to an inverted H α image from Finkbeiner (2003). Large circles represent sources with Gaia DR1 proper motions and parallaxes, while small circles are sources not detected by Gaia for which Hipparcos proper motions and parallaxes were used. The dashed lines show the division of the association into the three subgroups according to de Zeeuw et al. (1999), with the symbols coloured according to which of the three subgroups they nominally belong: red for US, green for UCL, and blue for LCC. Sources outside of these boundaries have been assigned to the nearest subgroup. Figure 1. View largeDownload slide Spatial distribution of Sco-Cen members (Rizzuto, Ireland & Robertson 2011) projected on to an inverted H α image from Finkbeiner (2003). Large circles represent sources with Gaia DR1 proper motions and parallaxes, while small circles are sources not detected by Gaia for which Hipparcos proper motions and parallaxes were used. The dashed lines show the division of the association into the three subgroups according to de Zeeuw et al. (1999), with the symbols coloured according to which of the three subgroups they nominally belong: red for US, green for UCL, and blue for LCC. Sources outside of these boundaries have been assigned to the nearest subgroup. 2.1 Gaia and Hipparcos astrometry Astrometric data for our sample comes from Gaia (Gaia Collaboration 2016b) data release 1 (DR1; Gaia Collaboration 2016a), which contains results from the Tycho-Gaia Astrometric Solution (TGAS; Michalik, Lindegren & Hobbs 2015; Lindegren et al. 2016) that combines astrometry from Gaia, Hipparcos (ESA 1997) and Tycho-2 (Høg et al. 2000). 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–2015). The parallaxes and proper motions calculated from these data are of particularly high precision, with uncertainties of 0.3 mas and 0.06 mas yr−1, respectively, though additional systematic uncertainties are still present (Gaia Collaboration 2016a; Lindegren et al. 2016). These data have 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 (Lindegren et al. 2016), 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 per cent). Fig. 2 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 (the necessary processing required to accurately measure astrometry for such sources was not established in DR1; Lindegren et al. 2016). Additionally, sources with extremely red or blue colours, or those with high proper-motions (>3.5 arcsec yr−1) are not included in DR1 (Gaia Collaboration 2016a; Arenou et al. 2017), though this will not include any of our targets. Where Gaia data are unavailable, we take parallaxes and proper motions from the revised Hipparcos catalogue (van Leeuwen 2007). Figure 2. View largeDownload slide Colour-magnitude diagrams for our Sco-Cen samples illustrating the completeness of the Gaia data (photometry from Anderson & Francis 2012). Sources with astrometry from Gaia are shown in red and those for which Hipparcos data was used are shown in blue. A 1 magnitude reddening vector is shown for illustrative purposes. The mean extinctions of members in each subgroup are AV = 0.73, 0.17, and 0.23 mag for US, UCL and LCC, respectively (de Geus, de Zeeuw & Lub 1989). Figure 2. View largeDownload slide Colour-magnitude diagrams for our Sco-Cen samples illustrating the completeness of the Gaia data (photometry from Anderson & Francis 2012). Sources with astrometry from Gaia are shown in red and those for which Hipparcos data was used are shown in blue. A 1 magnitude reddening vector is shown for illustrative purposes. The mean extinctions of members in each subgroup are AV = 0.73, 0.17, and 0.23 mag for US, UCL and LCC, respectively (de Geus, de Zeeuw & Lub 1989). The precision of the astrometric data used is shown in Fig. 3. The median proper motion precision is 0.06–0.11 mas yr−1 (≃0.04–0.08 km s−1 at a representative distance of 150 pc) in α and 0.05-0.08 mas yr−1 (≃0.04–0.06 km s−1) in δ for all the data for the three subgroups. When just Gaia data are considered, the median proper motion precision is 0.04–0.06 mas yr−1 (0.03–0.06 km s−1) in α and ∼0.04 mas yr−1 (0.03 km s−1) in δ. The median parallax uncertainty is ∼0.4 mas for all subsets of the data.2 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 (Lindegren et al. 2016). 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. Figure 3. View largeDownload slide Cumulative uncertainty distributions for the samples of stars in US (red), UCL (green), and LCC (blue) for parallax, proper motions, and RVs. Figure 3. View largeDownload slide Cumulative uncertainty distributions for the samples of stars in US (red), UCL (green), and LCC (blue) for parallax, proper motions, and RVs. 2.2 Radial velocities Results from Gaia's RV spectrometer are not included in DR1 and so we gathered RVs from the literature. We followed Murphy et al. (2015) in the prioritization of RVs from different samples, but chose not to exclude RVs from known binaries.3 Our primary source of RVs was the Magellan/MIKE study of Chen et al. (2011), from which 89 matches were found with our catalogue. Our next sources were the Pulkovo Compilation of Radial Velocities (PCRV, Gontcharov 2006), from which 118 matches were found, and CRVAD-2 (Kharchenko et al. 2007), 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 Dahm, Slesnick & White (2012). 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−1 (LCC). 3 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 4) and for assessing the evidence for linear expansion (Section 5). 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 7.1). 3.1 Individual distances Distances to individual sources and their uncertainties were calculated using the Bayesian inference method of Bailer-Jones (2015), 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) = r2 exp(−r/L), where L = 1.35 kpc is the scale length parameter, as recommended by Bailer-Jones (2015) and Astraatmadja & Bailer-Jones (2016).4 To fully sample the posterior distribution, we use the Markov Chain Monte Carlo (MCMC) ensemble sampler emcee (Foreman-Mackey et al. 2013), and use the mode of the corresponding posterior as our distance estimate and the 68 per cent 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. Fig. 4 shows the 3-dimensional spatial structure of Sco-Cen with the internal structure resolved in all three subgroups for the first time (de Zeeuw et al. 1999 were unable to resolve the internal structure of any of the subgroups, while de Bruijne 1999b 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. Figure 4. View largeDownload slide 3-dimensional spatial structure of Sco-Cen shown in the X − Y (top), X − Z (middle), and Y − Z planes (bottom panel). The sources are coloured according to their subgroup members: red for US, green for UCL, and blue for LCC. 1σ uncertainties on the positions of all sources are shown in grey. Also shown are the centres of each subgroup (large circles; see Section 3.2), the bulk motion of each subgroup over 2 Myr (arrows with black outline; see Section 4.1), and the relative motion of each subgroup over 10 Myr with the overall Sco-Cen bulk motion (U = −7.0 ± 0.2, V = −19.6 ± 0.2, W = −6.1 ± 0.1 km s−1) subtracted (arrow with grey outline). Figure 4. View largeDownload slide 3-dimensional spatial structure of Sco-Cen shown in the X − Y (top), X − Z (middle), and Y − Z planes (bottom panel). The sources are coloured according to their subgroup members: red for US, green for UCL, and blue for LCC. 1σ uncertainties on the positions of all sources are shown in grey. Also shown are the centres of each subgroup (large circles; see Section 3.2), the bulk motion of each subgroup over 2 Myr (arrows with black outline; see Section 4.1), and the relative motion of each subgroup over 10 Myr with the overall Sco-Cen bulk motion (U = −7.0 ± 0.2, V = −19.6 ± 0.2, W = −6.1 ± 0.1 km s−1) subtracted (arrow with grey outline). 3.2 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 parametrized 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, $$\boldsymbol{\theta }$$, best explain the observations, $$\boldsymbol{d}$$. In Bayes's theorem, this is known as the posterior distribution, $$\boldsymbol{P(\theta | d)}$$, and is given by   $$\boldsymbol{P(\theta | d)} = \frac{ \boldsymbol{P(d | \theta )} \times \boldsymbol{P(\theta )} }{ \boldsymbol{P(d)} }$$ (1)where $$\boldsymbol{P(d | \theta )}$$ is the likelihood model, $$\boldsymbol{P(\theta )}$$ are the priors (which include our a priori knowledge about the model parameters) and $$\boldsymbol{P(d)}$$ is a normalizing 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 (see discussion in Bailer-Jones 2015). 3.2.1 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 Elson, Fall & Freeman (EFF, 1987) radial profile, where the stellar surface density, Σ(r), at a radius, r, is given by   $$\Sigma (r) = \Sigma _0 \left( 1 + \frac{r^2}{a^2} \right)^{-\gamma / 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 \sqrt{4^{1/\gamma } - 1}$$, while the power-law index can be used to calculate the quantity η in the virial equation (e.g. Portegies Zwart, McMillan & Gieles 2010). 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 five (distance, central α and δ, a, and γ) and nine parameters (where a and γ become aX, aY, aZ, γX, γY, and γZ). We adopt liberal priors where possible, to minimize 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 per cent confidence interval used as the 1σ uncertainty. 3.2.2 Results The results of the model fits are provided in Table 1 and illustrated in Fig. 5. 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. Figure 5. View largeDownload slide Source distributions (black histogram) in parallax (top row), α and δ (bottom row) for the three subgroups of Sco-Cen compared to the EFF model fit results for 1D (red) and 3D (green) models. The model fit distributions were calculated by sampling the best-fitting forward models 1 000 000 times and plotting the resulting probability distribution functions. Figure 5. View largeDownload slide Source distributions (black histogram) in parallax (top row), α and δ (bottom row) for the three subgroups of Sco-Cen compared to the EFF model fit results for 1D (red) and 3D (green) models. The model fit distributions were calculated by sampling the best-fitting forward models 1 000 000 times and plotting the resulting probability distribution functions. Table 1. Structural and kinematic properties of the three Sco-Cen subgroups determined in this work. The structure parameters are all determined from the 1D EFF model, but these are generally in close agreement with the 3D EFF model. The distance uncertainties do not take into account the 0.3 mas systematic parallax uncertainty present in DR1. The velocity dispersion fitting results are taken from the model using a binary fraction of 100 per cent. Subgroup  US  UCL  LCC  Distance (pc)  $$143.0^{+0.3}_{-0.4}$$  $$135.9^{+0.5}_{-0.4}$$  $$115.2^{+0.3}_{-0.3}$$  α0 (deg)  $$243.02^{+0.15}_{-0.10}$$  $$228.36^{+0.27}_{-0.24}$$  $$188.57^{+0.23}_{-0.19}$$  δ0 (deg)  $$-24.19^{+0.08}_{-0.18}$$  $$-40.50^{+0.13}_{-0.13}$$  $$-54.53^{+0.20}_{-0.06}$$  X0 (pc)  $$133.1^{+0.4}_{-0.3}$$  $$114.0^{+0.5}_{-0.5}$$  $$57.7^{+0.3}_{-0.3}$$  Y0 (pc)  $$-21.3^{+0.2}_{-0.4}$$  $$-65.5^{+0.4}_{-0.5}$$  $$-98.3^{+0.3}_{-0.3}$$  Z0 (pc)  $$47.6^{+0.3}_{-0.3}$$  $$34.6^{+0.4}_{-0.3}$$  $$16.6^{+0.3}_{-0.2}$$  a (pc)  $$35.8^{+1.9}_{-1.3}$$  $$69.8^{+3.3}_{-3.2}$$  $$50.1^{+4.0}_{-3.2}$$  γ  $$14.5^{+1.4}_{-1.0}$$  $$17.6^{+1.3}_{-1.3}$$  $$15.2^{+1.9}_{-1.5}$$  rhl (pc)  $$11.4^{+0.2}_{-0.2}$$  $$20.0^{+0.3}_{-0.2}$$  $$15.5^{+0.3}_{-0.2}$$  U0 (km s−1)  $$-6.16_{-0.13}^{+0.14}$$  $$-5.90_{-0.12}^{+0.17}$$  $$-8.96_{-0.24}^{+0.14}$$  V0 (km s−1)  $$-16.89_{-0.10}^{+0.08}$$  $$-20.00_{-0.10}^{+0.12}$$  $$-20.55_{-0.13}^{+0.16}$$  W0 (km s−1)  $$-7.05_{-0.08}^{+0.09}$$  $$-5.80_{-0.09}^{+0.09}$$  $$-6.29_{-0.11}^{+0.12}$$  σU (km s−1)  $$1.63_{-0.20}^{+0.20}$$  $$1.96_{-0.12}^{+0.13}$$  $$1.89_{-0.13}^{+0.21}$$  σV (km s−1)  $$1.14_{-0.14}^{+0.13}$$  $$0.73_{-0.19}^{+0.15}$$  $$0.90_{-0.30}^{+0.50}$$  σW (km s−1)  $$2.51_{-0.09}^{+0.11}$$  $$1.27_{-0.14}^{+0.11}$$  $$0.51_{-0.18}^{+0.25}$$  σ3D (km s−1)  $$3.20_{-0.20}^{0.22}$$  $$2.45_{-0.20}^{+0.20}$$  $$2.15_{-0.24}^{+0.47}$$  αCP (deg)  $$116.22_{-9.46}^{+10.70}$$  $$100.88_{-1.52}^{+1.69}$$  $$95.32_{-1.31}^{+1.29}$$  δCP (deg)  $$-55.29_{-4.56}^{+5.37}$$  $$-35.88_{-1.61}^{+1.63}$$  $$-28.68_{-1.48}^{+1.53}$$  Subgroup  US  UCL  LCC  Distance (pc)  $$143.0^{+0.3}_{-0.4}$$  $$135.9^{+0.5}_{-0.4}$$  $$115.2^{+0.3}_{-0.3}$$  α0 (deg)  $$243.02^{+0.15}_{-0.10}$$  $$228.36^{+0.27}_{-0.24}$$  $$188.57^{+0.23}_{-0.19}$$  δ0 (deg)  $$-24.19^{+0.08}_{-0.18}$$  $$-40.50^{+0.13}_{-0.13}$$  $$-54.53^{+0.20}_{-0.06}$$  X0 (pc)  $$133.1^{+0.4}_{-0.3}$$  $$114.0^{+0.5}_{-0.5}$$  $$57.7^{+0.3}_{-0.3}$$  Y0 (pc)  $$-21.3^{+0.2}_{-0.4}$$  $$-65.5^{+0.4}_{-0.5}$$  $$-98.3^{+0.3}_{-0.3}$$  Z0 (pc)  $$47.6^{+0.3}_{-0.3}$$  $$34.6^{+0.4}_{-0.3}$$  $$16.6^{+0.3}_{-0.2}$$  a (pc)  $$35.8^{+1.9}_{-1.3}$$  $$69.8^{+3.3}_{-3.2}$$  $$50.1^{+4.0}_{-3.2}$$  γ  $$14.5^{+1.4}_{-1.0}$$  $$17.6^{+1.3}_{-1.3}$$  $$15.2^{+1.9}_{-1.5}$$  rhl (pc)  $$11.4^{+0.2}_{-0.2}$$  $$20.0^{+0.3}_{-0.2}$$  $$15.5^{+0.3}_{-0.2}$$  U0 (km s−1)  $$-6.16_{-0.13}^{+0.14}$$  $$-5.90_{-0.12}^{+0.17}$$  $$-8.96_{-0.24}^{+0.14}$$  V0 (km s−1)  $$-16.89_{-0.10}^{+0.08}$$  $$-20.00_{-0.10}^{+0.12}$$  $$-20.55_{-0.13}^{+0.16}$$  W0 (km s−1)  $$-7.05_{-0.08}^{+0.09}$$  $$-5.80_{-0.09}^{+0.09}$$  $$-6.29_{-0.11}^{+0.12}$$  σU (km s−1)  $$1.63_{-0.20}^{+0.20}$$  $$1.96_{-0.12}^{+0.13}$$  $$1.89_{-0.13}^{+0.21}$$  σV (km s−1)  $$1.14_{-0.14}^{+0.13}$$  $$0.73_{-0.19}^{+0.15}$$  $$0.90_{-0.30}^{+0.50}$$  σW (km s−1)  $$2.51_{-0.09}^{+0.11}$$  $$1.27_{-0.14}^{+0.11}$$  $$0.51_{-0.18}^{+0.25}$$  σ3D (km s−1)  $$3.20_{-0.20}^{0.22}$$  $$2.45_{-0.20}^{+0.20}$$  $$2.15_{-0.24}^{+0.47}$$  αCP (deg)  $$116.22_{-9.46}^{+10.70}$$  $$100.88_{-1.52}^{+1.69}$$  $$95.32_{-1.31}^{+1.29}$$  δCP (deg)  $$-55.29_{-4.56}^{+5.37}$$  $$-35.88_{-1.61}^{+1.63}$$  $$-28.68_{-1.48}^{+1.53}$$  View Large 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 overestimated then these uncertainties will also be overestimated). Our values are consistent with the mean distances found by de Zeeuw et al. (1999) 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 Rizzuto et al. (2011), and with slightly larger uncertainties (because of the systematic and spatially-correlated 0.3 mas parallax uncertainty in DR1 that is not reduced by averaging, Gaia Collaboration 2016a). 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 per cent 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 (e.g. Rizzuto et al. 2011; Pecaut & Mamajek 2016) 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 are 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. 4 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. 4.1 Forward modelling the velocity dispersion To measure the velocity dispersion and bulk motion of each subgroup, we use Bayesian inference, as in Section 3.2. 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. 4.1.1 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 3.1. 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 Odenkirchen et al. (2002) and Cottaar et al. (2012), 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 (Kroupa et al. 2001) 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 Maschberger (2013). We use the orbital properties of intermediate-mass stars in Sco-Cen determined by Kouwenhoven et al. (2007), who performed an extensive study of the visual, spectroscopic, and astrometric binaries in the association. They inferred a binary fraction of 100 per cent, a mass ratio distribution fq(q) ∝ q−0.4, and a semi-major axis distribution of the form fa(a) ∝ a−1 (in the range 5 ≥ a/R⊙ ≥ 5 × 106). 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 Kouwenhoven et al. (2007) 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 106 random binaries, providing us with a distribution from which we can sample. We don't consider triple systems because their properties are poorly constrained (Duchêne & Kraus 2013) 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 six free parameters, the central velocity (U0, V0, W0) and velocity dispersion (σU, σV, σW) in each dimension, for which we adopt flat and wide priors of − 100 ≤ v0[km s− 1] ≤ +100 for the central velocities and 0 ≤ σ[km s−1] ≤ 100 for the velocity dispersions. As in Section 3.2, 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 per cent confidence interval used for the 1σ uncertainty. 4.1.2 Results The results are provided in Table 1 and illustrated in Fig. 6, with the central velocities, or bulk motions, also shown in Fig. 4. The bulk motions are different from those presented by de Zeeuw et al. (1999), particularly in U and V, implying that the difference is probably due to the increased availability and quality of RVs for these stars (though the change in Sco-Cen membership by Rizzuto et al. 2011, will have had a small effect). Our bulk motions are in better agreement with those presented by Chen et al. (2011), with most agreeing within 1σ–2σ, and the remaining differences probably due to the change in membership. Figure 6. View largeDownload slide Source distributions (black histograms) for proper motions (left and centre) and RVs (right) for the three subgroups, compared to the best-fitting velocity dispersion model fit results for binary fractions of 100 per cent (green) and 70 per cent (red). The model fit distributions were calculated by sampling the best-fitting forward models 1 000 000 times and plotting the resulting probability distribution functions. Figure 6. View largeDownload slide Source distributions (black histograms) for proper motions (left and centre) and RVs (right) for the three subgroups, compared to the best-fitting velocity dispersion model fit results for binary fractions of 100 per cent (green) and 70 per cent (red). The model fit distributions were calculated by sampling the best-fitting forward models 1 000 000 times and plotting the resulting probability distribution functions. The best-fitting velocity dispersions vary from 0.5 to 2.5 km s−1 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. de Bruijne (1999b) found that the 1D velocity dispersions for all the subgroups are ≤1.0–1.5 km s−1, while Madsen et al. (2002) measured velocity dispersions of 1.1–1.3 km s−1 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−1 (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-Cen5 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). Kouwenhoven et al. (2007) find that the observations of Sco-Cen are best fit with a binary fraction of 100 per cent, though slightly smaller fractions are possible. We explored using binary fractions of 70–100 per cent 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 per cent along that axis for each 10 per cent by which the binary fraction is reduced. Fig. 6 shows the model fit results for a 70 per cent binary fraction, with little difference in the fit quality compared to that for 100 per cent binaries. 4.2 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 5.3). In this section we measure the convergent points of the three subgroups using both the refurbished Jones (1971) method (de Bruijne 1999a) and directly from the 3D bulk motions. 4.2.1 Method Jones's method iterates over a hemispherical grid of points to find the convergent point coordinates that minimizes the proper motions perpendicular to the great circle joining each star to the convergent point. The best-fitting convergent point is that which minimizes the quantity $$\sum t^2_\perp$$, 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 de Bruijne (1999a), who changed this quantity to $$t_\perp = \mu _\perp / \sqrt{\sigma ^2_\perp + \sigma ^2_v}$$, recognizing 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_\perp ^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_\perp ^2$$ will exceed the observed value of $$t_\perp ^2$$ by chance even for a correct model. Therefore, if the value of $$t_\perp ^2$$ is too high, the probability ε will be below the threshold value, which de Bruijne (1999a) recommends setting to a value of 0.954 (giving an ∼5 per cent 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_\perp ^2$$ value. 4.2.2 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 4.1.6 For US a convergent point for the entire subgroup was found with a χ2 probability of 53 per cent, rising to 95.4 per cent 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−1 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−1. Fig. 7 shows this convergent point compared to those determined previously by de Zeeuw et al. (1999) and Madsen et al. (2002) 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σ. Figure 7. View largeDownload slide Convergent points for the three subgroups of Sco-Cen determined using the Jones (1971) method (red circle and contours showing 1σ, 2σ, and 3σ uncertainties) and from the bulk motions (green circles with 1σ error bars typically equal to or smaller than the size of the symbols). Yellow circles show previously published convergent points from de Zeeuw et al. (1999, DZ99) and Madsen, Dravins & Lindegren (2002, M02), with 1σ error bars shown for the latter (uncertainties for the former are expected to be several degrees along the great circle connecting the convergent point with the subgroup). Dashed lines with black markers show the predicted convergent points derived from the mean space motions at a range of expansion ages. Figure 7. View largeDownload slide Convergent points for the three subgroups of Sco-Cen determined using the Jones (1971) method (red circle and contours showing 1σ, 2σ, and 3σ uncertainties) and from the bulk motions (green circles with 1σ error bars typically equal to or smaller than the size of the symbols). Yellow circles show previously published convergent points from de Zeeuw et al. (1999, DZ99) and Madsen, Dravins & Lindegren (2002, M02), with 1σ error bars shown for the latter (uncertainties for the former are expected to be several degrees along the great circle connecting the convergent point with the subgroup). Dashed lines with black markers show the predicted convergent points derived from the mean space motions at a range of expansion ages. For UCL a convergent point for the entire subgroup was found with a χ2 probability of 98 per cent 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 per cent, rising to 95.4 per cent after removing the 10 stars with the highest |t⊥| values. A higher velocity dispersion of 1.5 km s−1 would result in an immediate fit with a χ2 probability of >99 per cent 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−1 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. 5 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 (Blaauw 1964b). 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. 5.1 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 Blaauw (1964b) linear expansion model individual RVs, vrad, are predicted to follow the relation   $$v_\mathrm{rad} = v^\prime \, \mathrm{cos} \, \lambda ^\prime + \kappa d + K,$$ (3)where λ΄ is the angular separation between star and convergent point (calculated using just the proper motions; see Section 4.2), 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^\prime = \frac{ \mu _v A }{ \varpi \, \mathrm{sin} \, \lambda ^\prime },$$ (4)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). Fig. 8 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 (see Hogg, Bovy & Lang 2010) and then marginalizing over this parameter to obtain the uncertainties on the resulting fit parameters. Figure 8. View largeDownload slide Distance to sources in each subgroup of Sco-Cen plotted against the difference between the observed RV and that predicted from the moving group method (v΄ cos λ). The error bars for individual sources show the 16–84 per cent confidence intervals, calculated from Monte Carlo simulations. The solid lines show the best-fitting linear relationship between the plotted quantities, with 1σ uncertainties shown with dashed lines. The best-fitting slopes, κ, and uncertainties are noted in each panel. If the subgroups are not expanding, we would expect a slope of zero, while if they expanding, we would expect a positive slope. 3σ outliers on the ordinate, shown with non-coloured symbols, were excluded from the fit. Figure 8. View largeDownload slide Distance to sources in each subgroup of Sco-Cen plotted against the difference between the observed RV and that predicted from the moving group method (v΄ cos λ). The error bars for individual sources show the 16–84 per cent confidence intervals, calculated from Monte Carlo simulations. The solid lines show the best-fitting linear relationship between the plotted quantities, with 1σ uncertainties shown with dashed lines. The best-fitting slopes, κ, and uncertainties are noted in each panel. If the subgroups are not expanding, we would expect a slope of zero, while if they expanding, we would expect a positive slope. 3σ outliers on the ordinate, shown with non-coloured symbols, were excluded from the fit. 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 Pecaut et al. (2012). 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 Pecaut & Mamajek (2016) 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. 5.2 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. Fig. 9 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. Figure 9. View largeDownload slide Positions versus velocities along the three Galactic Cartesian axes XYZ for each of the three subgroups of Sco-Cen (US on the left, UCL in the centre, and LCC on the right). 1σ error bars are shown for all sources. The solid lines show the best-fit linear relationships between the plotted quantities, with 1σ uncertainties shown with dashed lines. The best-fitting slopes, κ, and uncertainties are noted in each panel. If the subgroups are not expanding we would expect slopes of zero, while if they are expanding we would expect a positive slope (a negative slope implies contraction along that axis). 3σ outliers on the ordinate, shown with non-coloured symbols, were excluded from the fits. Figure 9. View largeDownload slide Positions versus velocities along the three Galactic Cartesian axes XYZ for each of the three subgroups of Sco-Cen (US on the left, UCL in the centre, and LCC on the right). 1σ error bars are shown for all sources. The solid lines show the best-fit linear relationships between the plotted quantities, with 1σ uncertainties shown with dashed lines. The best-fitting slopes, κ, and uncertainties are noted in each panel. If the subgroups are not expanding we would expect slopes of zero, while if they are expanding we would expect a positive slope (a negative slope implies contraction along that axis). 3σ outliers on the ordinate, shown with non-coloured symbols, were excluded from the fits. 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 significance 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. 5.3 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 3) and compare these with the proper motion convergent points. Fig. 7 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 (Pecaut & Mamajek 2016). 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 Pecaut & Mamajek (2016). 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. 5.4 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) = X0 + γUt, 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. Fig. 10 shows the 1σ dispersions (16th–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. Figure 10. View largeDownload slide 1σ dispersions (16th–84th percentiles used to minimize the impact of outliers) of the size of each subgroup in each axis of the XYZ coordinate system (σX shown with a dashed line, σY with a dotted line, and σZ with a dash-dotted line) as well as the quadrature sum of all three dimensions (solid line). Black lines show the results of traceback analysis performed with linear trajectories and the red lines show the analysis performed using the epicycle orbit approximation. Shaded areas around the 3D sum lines show the 90 per cent confidence interval in the quadrature sum dispersion of each subgroup determined from Monte Carlo simulations exploring the impact of uncertainties in the velocities. Figure 10. View largeDownload slide 1σ dispersions (16th–84th percentiles used to minimize the impact of outliers) of the size of each subgroup in each axis of the XYZ coordinate system (σX shown with a dashed line, σY with a dotted line, and σZ with a dash-dotted line) as well as the quadrature sum of all three dimensions (solid line). Black lines show the results of traceback analysis performed with linear trajectories and the red lines show the analysis performed using the epicycle orbit approximation. Shaded areas around the 3D sum lines show the 90 per cent confidence interval in the quadrature sum dispersion of each subgroup determined from Monte Carlo simulations exploring the impact of uncertainties in the velocities. For a more accurate traceback of the stellar motions, we use the epicycle approximation, employing the orbital equations from Fuchs et al. (2006). We use the Oort A and B constants from Feast & Whitelock (1997), the local disc density from Holmberg & Flynn (2004), the local standard of rest velocity from Schönrich, Binney & Dehnen (2010) and a solar Z distance above the Galactic plane of 17 pc (Karim & Mamajek 2017). Fig. 10 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 Fig. 10 we show the 90 per cent 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 Brown et al. (1997) 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 Brown et al. (1997) however. 5.5 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 (referred to as virtual expansion /contraction; Brown et al. 1997). 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 Brown et al. (1997), using the median RV of stars in each subgroup (to minimize the impact of binarity). The corrected proper motion vector maps for the three subgroups are shown in Fig. 11. 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 6, that these patterns are due to the kinematic substructure within the subgroups. Figure 11. View largeDownload slide Proper motion vector maps for the three subgroups of Sco-Cen projected on to inverted H α images from Finkbeiner (2003). The proper motions have been corrected for radial streaming motion according to the equations in Brown et al. (1997) and using bulk (median) RVs of −6.2, 2.8, and 13.0 km s−1 for the three subgroups. Points show the current positions of stars and vectors show the proper motions over 0.5 Myr with the bulk motion of each subgroup subtracted to show the motion in the reference frame of each subgroup. The vectors have been colour-coded based on their direction of motion (see colour wheel in lower-right corner) to highlight kinematic substructure. Note that because the kinematic reference frame shown for each subgroup is different, the colours of the proper motions in each subgroup do not represent the same absolute proper motion as stars in other subgroups. Figure 11. View largeDownload slide Proper motion vector maps for the three subgroups of Sco-Cen projected on to inverted H α images from Finkbeiner (2003). The proper motions have been corrected for radial streaming motion according to the equations in Brown et al. (1997) and using bulk (median) RVs of −6.2, 2.8, and 13.0 km s−1 for the three subgroups. Points show the current positions of stars and vectors show the proper motions over 0.5 Myr with the bulk motion of each subgroup subtracted to show the motion in the reference frame of each subgroup. The vectors have been colour-coded based on their direction of motion (see colour wheel in lower-right corner) to highlight kinematic substructure. Note that because the kinematic reference frame shown for each subgroup is different, the colours of the proper motions in each subgroup do not represent the same absolute proper motion as stars in other subgroups. To quantify the expansion of each subgroup, we follow the method of Wright et al. (2016) and calculate, for each star, the amount of kinetic energy in the radial and transverse directions (using the subgroup centres calculated in Section 3.2), 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 per cent (US), 66 per cent (UCL), and 45 per cent (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 per cent (US), 67 per cent (UCL), and 90 per cent (LCC) of the kinetic energy in the direction of expansion. This suggests that, while there is not evidence for coherent expansion patterns in the correct proper motions of stars in each subgroup, they do appear to be expanding. 6 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 (e.g. the Gamma Velorum cluster shows two distinct kinematic components; Jeffries et al. 2014) or as smaller, more numerous substructures (as observed in Cygnus OB2, Wright et al. 2016). 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 (e.g. Larson 1981; Wright et al. 2016). In this section we explore the evidence for kinematic substructure in Sco-Cen. Fig. 12 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 (as per Wright et al. 2016). 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. Figure 12. View largeDownload slide Proper motion vector map for the entire Sco-Cen association projected on to an inverted H α image from Finkbeiner (2003). The proper motions were corrected for radial streaming motion according to the equations in Brown et al. (1997) and using bulk (median) RVs of −6.2, 2.8, and 13.0 km s−1 for the three subgroups. Points show the current positions of stars and vectors show the proper motions over 0.5 Myr with the bulk (median) proper motion of the entire association subtracted to show the relative motion. The vectors have been colour-coded based on their direction of motion (see colour wheel in lower left) to highlight kinematic substructure. Figure 12. View largeDownload slide Proper motion vector map for the entire Sco-Cen association projected on to an inverted H α image from Finkbeiner (2003). The proper motions were corrected for radial streaming motion according to the equations in Brown et al. (1997) and using bulk (median) RVs of −6.2, 2.8, and 13.0 km s−1 for the three subgroups. Points show the current positions of stars and vectors show the proper motions over 0.5 Myr with the bulk (median) proper motion of the entire association subtracted to show the relative motion. The vectors have been colour-coded based on their direction of motion (see colour wheel in lower left) to highlight kinematic substructure. 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 (e.g. David & Hillenbrand 2015) and therefore suggest that these are field stars that should be removed from the membership of Sco-Cen. Fig. 11 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 5.2. 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. 6.1 Quantifying kinematic substructure 6.1.1 Method In an attempt to quantify this substructure we follow Wright et al. (2016) and use Moran's spatial correlation index I (Moran 1950), 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 I0 = −1/(n − 1), with I > I0, indicating positive spatial correlation and I < I0, indicating negative spatial correlation. We use this instead of Geary's C index (Geary 1954), since Moran's I statistic is a global measure that can identify and quantify large-scale spatial correlation (as appears to be evident in Fig. 11), 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 Wright et al. (2016). 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 on to 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 4.1). We then calculated Moran's I index for the simulated data set 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. 6.1.2 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 Fig. 11 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 (Wright et al. 2016). 7 DISCUSSION 7.1 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   $$\sigma _{3{\rm D}}^2 = \frac{G M_{{\rm vir}}}{2 r_{{\rm vir}}}$$ (5)where σ3D is the 3D velocity dispersion, G is the gravitational constant, Mvir is the virial mass, and rvir is the virial radius. We substitute the parameter η = 6rvir/reff (e.g. Portegies Zwart et al. 2010), where reff is the effective (or half-light) radius, giving   $$M_{{\rm vir}} = \eta \, \frac{\sigma _{3D}^2 \, r_{{\rm eff}}}{3G}$$ (6)where the parameters η and reff can be calculated from the Elson et al. (1987) parameters γ and a (see Section 3.2). For the large γ values calculated for the subgroups we find η values of 9.1 (US), 9.0 (UCL), and 9.1 (LCC),7 while the effective radii, reff 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 1), and η values, we calculate virial masses of $$8.2^{+1.5}_{-1.2} \times 10^4$$ (US), $$8.3^{+1.8}_{-1.5} \times 10^4$$ (UCL), and $$5.1^{+2.9}_{-1.3} \times 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. Preibisch & Mamajek (2008) estimate US to have a total mass of ∼2000 M⊙, while Mamajek et al. (2002) 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⊙ (using the same initial mass function and binary fraction as Preibisch & Mamajek 2008). 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 de Geus (1992) estimates to have a mass of 8 × 104 M⊙. The author argues 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 × 105 and 1 × 105 M⊙, respectively, are more massive than their virial masses, implying that the subgroups might have initially been in a subvirial state. 7.2 Expansion of the Sco-Cen subgroups Ever since Ambartsumian (1947) 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 (e.g. Blaauw 1964a, 1991; Brown et al. 1999). These ideas lead to suggestions that OB associations may be the expanded remnants of compact star clusters (Lada & Lada 2003), disrupted by processes such as residual gas expulsion (e.g. Hills 1980; Lada et al. 1984). In Section 5 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 5.1 and 5.2), while others better account for the Galactic potential (Section 5.4) or provide an overall view of the stellar motions (Section 5.5). 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 5.1 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 5.2) 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, Mamajek & Bell (2014) 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 time-scales involved are probably too short (Dobbs & Pringle 2013, estimate a shear time-scale of ∼70 Myr in the vicinity of the Sun). 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 5.4), 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 Fig. 11 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 (Mamajek & Bell 2014) and Tuc-Hor (Makarov 2007), 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 (e.g. Pfalzner 2009). This is particularly clear in the case of the three subgroups of Sco-Cen, which Pfalzner & Kaczmarek (2013) 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 Fig. 11 suggest a more complex expansion pattern consisting of multiple expanding substructures within each subgroup. 7.3 Kinematic substructure and the formation of Sco-Cen Molecular clouds are known to be highly substructured, both spatially (Falgarone, Phillips & Walker 1991) and kinematically (Hacar et al. 2013), with this hierarchical or fractal structure passed on to the formed stars (Elmegreen & Elmegreen 2001) and often evident in their initial spatial distribution (Gutermuth et al. 2008). Dynamical interactions between groups of stars can erase this substructure (e.g. Scally & Clarke 2002; Parker et al. 2014), a process that acts on the dynamical time-scale of a region. If the dynamical time-scale is longer than the age of the region then the primordial substructure can be preserved (e.g. Goodwin & Whitworth 2004), though dynamical time-scales can change as a region expands or collapses (so an instantaneous measure of the dynamical time-scale is not 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 time-scales 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 (e.g. Parker et al. 2014). The three subgroups of Sco-Cen each display evidence for kinematic substructure, as has been observed in other OB associations (Wright et al. 2014) and across large star-forming complexes (e.g. Fűrész et al. 2006, 2008; Tobin et al. 2015). 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 do not 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 Pecaut & Mamajek (2016), 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 Blaauw (1964a) and extended by Preibisch & Zinnecker (1999), the impact of a high-velocity cloud on the disc (Lepine & Duvert 1994), to its formation in a large ring-like structure containing other young groups of stars known as the Gould Belt (Lindblad et al. 1973). 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 (e.g. Sartori, Lépine & Dias 2003; Preibisch & Mamajek 2008). 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 does not necessarily disprove it. 8 CONCLUSIONS Our kinematic study of the Sco-Cen OB association using Gaia DR1 parallaxes and proper motions has led 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 that 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−1 (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 Blaauw (1964b) 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 (Wright et al. 2016), 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 Pecaut & Mamajek (2016) 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 (e.g. Makarov 2007; Mamajek & Bell 2014; Wright et al. 2016) 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 (Pfalzner 2009) and that star cluster disruption mechanisms, such as residual gas expulsion (Kroupa et al. 2001), 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 discs. Acknowledgements NJW acknowledges an Science and Technology Facilities Council 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 European Space Agency 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 Centre de Données astronomiques de Strasbourg, Strasbourg, France), Astropy (Astropy Collaboration 2013) and TOPCAT (Tool for OPerations on Catalogues And Tables, Taylor 2005). NOTE ADDED IN PROOF At the time this paper was accepted two papers came to our attention that, respectively, find evidence for expansion of the OB associations Per OB1 and Car OB1 (Mel'nik & Dambis 2017) and no evidence for OB association expansion (Ward & Kruijssen 2018), both based on Gaia proper motions. Footnotes 1 Murphy, Lawson & Bento (2015) note that some of the RVs used by Rizzuto et al. (2011) 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. 2 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. 3 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. 4 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 3.2 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. 5 We repeated these fits, using the de Zeeuw et al. (1999) membership list, and find that the change in sample to Rizzuto et al. (2011) has caused the velocity dispersion to increase by ∼10 per cent. 6 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. 7 The dependence of η on γ is very weak at large γ values and therefore the uncertainties on η are almost zero. REFERENCES Adams F. C., Proszkow E. M., Fatuzzo M., Myers P. C., 2006, ApJ , 641, 504 https://doi.org/10.1086/500393 CrossRef Search ADS   Ambartsumian V. A., 1947, Stellar Evolution and Astrophysics . Armenian Acad. Sci. Anderson E., Francis C., 2012, Astron. Lett. , 38, 331 https://doi.org/10.1134/S1063773712050015 CrossRef Search ADS   Arenou F. et al.  , 2017, A&A , 599, A50 CrossRef Search ADS   Astraatmadja T. L., Bailer-Jones C. A. L., 2016, ApJ , 832, 137 https://doi.org/10.3847/0004-637X/832/2/137 CrossRef Search ADS   Astropy Collaboration, 2013, A&A , 558, A33 CrossRef Search ADS   Bailer-Jones C. A. L., 2015, PASP , 127, 994 https://doi.org/10.1086/683116 CrossRef Search ADS   Baumgardt H., Kroupa P., 2007, MNRAS , 380, 1589 https://doi.org/10.1111/j.1365-2966.2007.12209.x CrossRef Search ADS   Blaauw A., 1946, Publi. Kapteyn Astron. Lab. Groningen , 52, 1 Blaauw A., 1964a, ARA&A , 2, 213 CrossRef Search ADS   Blaauw A., 1964b, in Kerr F. J., ed., IAU Symp., Vol. 20, The Galaxy and the Magellanic Clouds . Australian Academy of Science, Canberra, p. 50 Blaauw A., 1991, in Lada C. J. Kylafis N. D., eds., NATO ASIC Proc. 342: The Physics of Star Formation and Early Stellar Evolution . Kluwer, Dordrecht, p. 125 Google Scholar CrossRef Search ADS   Breitschwerdt D., Feige J., Schulreich M. M., Avillez M. A. D., Dettbarn C., Fuchs B., 2016, Nature , 532, 73 https://doi.org/10.1038/nature17424 CrossRef Search ADS PubMed  Bressert E. et al.  , 2010, MNRAS , 409, L54 https://doi.org/10.1111/j.1745-3933.2010.00946.x CrossRef Search ADS   Brown A. G. A., Dekker G., de Zeeuw P. T., 1997, MNRAS , 285, 479 https://doi.org/10.1093/mnras/285.3.479 CrossRef Search ADS   Brown A. G. A. Blaauw A. Hoogerwerf R. de Bruijne J. H. J. de Zeeuw P. T., 1999, in Series C, Lada C. J., Kylafis N. D., eds., NATO Advanced Science Institutes (ASI) Vol. 540, p. 411 Carpenter J. M., Mamajek E. E., Hillenbrand L. A., Meyer M. R., 2006, ApJ , 651, L49 https://doi.org/10.1086/509121 CrossRef Search ADS   Chen C. H., Mamajek E. E., Bitner M. A., Pecaut M., Su K. Y. L., Weinberger A. J., 2011, ApJ , 738, 122 https://doi.org/10.1088/0004-637X/738/2/122 CrossRef Search ADS   Cottaar M., Meyer M. R., Andersen M., Espinoza P., 2012, A&A , 539, A5 CrossRef Search ADS   Dahm S. E., Slesnick C. L., White R. J., 2012, ApJ , 745, 56 https://doi.org/10.1088/0004-637X/745/1/56 CrossRef Search ADS   David T. J., Hillenbrand L. A., 2015, ApJ , 804, 146 https://doi.org/10.1088/0004-637X/804/2/146 CrossRef Search ADS   de Bruijne J. H. J., 1999a, MNRAS , 306, 381 https://doi.org/10.1046/j.1365-8711.1999.02643.x CrossRef Search ADS   de Bruijne J. H. J., 1999b, MNRAS , 310, 585 https://doi.org/10.1046/j.1365-8711.1999.02953.x CrossRef Search ADS   de Geus E. J., 1992, A&A , 262, 258 de Geus E. J., de Zeeuw P. T., Lub J., 1989, A&A , 216, 44 de Zeeuw P. T., Hoogerwerf R., de Bruijne J. H. J., Brown A. G. A., Blaauw A., 1999, AJ , 117, 354 https://doi.org/10.1086/300682 CrossRef Search ADS   Dobbs C. L., Pringle J. E., 2013, MNRAS , 432, 653 https://doi.org/10.1093/mnras/stt508 CrossRef Search ADS   Duchêne G., Kraus A., 2013, ARA&A , 51, 269 CrossRef Search ADS   Elmegreen B. G. Efremov Y. Pudritz R. E. Zinnecker H., 2000, Protostars Planets IV , 179 Elmegreen B. G., Elmegreen D. M., 2001, AJ , 121, 1507 https://doi.org/10.1086/319416 CrossRef Search ADS   Elson R. A. W., Fall S. M., Freeman K. C., 1987, ApJ , 323, 54 https://doi.org/10.1086/165807 CrossRef Search ADS   ESA ed., 1997, ESA Special Publication, Vol. 1200, The HIPPARCOS and TYCHO Catalogues. Astrometric and Photometric Star Catalogues Derived from the ESA HIPPARCOS Space Astrometry Mission . ESA. Available at: https://www.cosmos.esa.int/documents/532822/552851/vol1_all.pdf/99adf6e3-6893-4824-8fc2-8d3c9cbba2b5 Falgarone E., Phillips T. G., Walker C. K., 1991, ApJ , 378, 186 https://doi.org/10.1086/170419 CrossRef Search ADS   Feast M., Whitelock P., 1997, MNRAS , 291, 683 https://doi.org/10.1093/mnras/291.4.683 CrossRef Search ADS   Finkbeiner D. P., 2003, ApJS , 146, 407 https://doi.org/10.1086/374411 CrossRef Search ADS   Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP , 125, 306 https://doi.org/10.1086/670067 CrossRef Search ADS   Fuchs B., Breitschwerdt D., de Avillez M. A., Dettbarn C., Flynn C., 2006, MNRAS , 373, 993 https://doi.org/10.1111/j.1365-2966.2006.11044.x CrossRef Search ADS   Fűrész G. et al.  , 2006, ApJ , 648, 1090, https://doi.org/10.1086/506140 CrossRef Search ADS   Fűrész G., Hartmann L. W., Megeath S. T., Szentgyorgyi A. H., Hamden E. T., 2008, ApJ , 676, 1109 https://doi.org/10.1086/525844 CrossRef Search ADS   Gaia Collaboration, 2016a, A&A , 595, A2 CrossRef Search ADS   Gaia Collaboration, 2016b, A&A , 595, A1 CrossRef Search ADS   Geary R. C., 1954, Incorp. Stat. , 5, 115 https://doi.org/10.2307/2986645 Gontcharov G. A., 2006, Astron. Lett. , 32, 759 https://doi.org/10.1134/S1063773706110065 CrossRef Search ADS   Goodwin S. P., Whitworth A. P., 2004, A&A , 413, 929 CrossRef Search ADS   Guarcello M. G. et al.  , 2016, preprint (arXiv:1605.01773) Gutermuth R. A. et al.  , 2008, ApJ , 674, 336 https://doi.org/10.1086/524722 CrossRef Search ADS   Hacar A., Tafalla M., Kauffmann J., Kovács A., 2013, A&A , 554, A55 CrossRef Search ADS   Hills J. G., 1980, ApJ , 235, 986 https://doi.org/10.1086/157703 CrossRef Search ADS   Hogg D. W. Bovy J. Lang D., 2010, preprint (arXiv:1008.4686) Holmberg J., Flynn C., 2004, MNRAS , 352, 440 https://doi.org/10.1111/j.1365-2966.2004.07931.x CrossRef Search ADS   Hoogerwerf R., de Bruijne J. H. J., de Zeeuw P. T., 2000, ApJ , 544, L133 https://doi.org/10.1086/317315 CrossRef Search ADS   Høg E. et al.  , 2000, A&A , 355, L27 Janson M., Lafrenière D., Jayawardhana R., Bonavita M., Girard J. H., Brandeker A., Gizis J. E., 2013, ApJ , 773, 170 https://doi.org/10.1088/0004-637X/773/2/170 CrossRef Search ADS   Jeffries R. D. et al.  , 2014, A&A , 563, A94 CrossRef Search ADS   Jones D. H. P., 1971, MNRAS , 152, 231 https://doi.org/10.1093/mnras/152.2.231 CrossRef Search ADS   Karim M. T., Mamajek E. E., 2017, MNRAS , 465, 472 https://doi.org/10.1093/mnras/stw2772 CrossRef Search ADS   Kharchenko N. V., Scholz R.-D., Piskunov A. E., Röser S., Schilbach E., 2007, Astron. Nachr. , 328, 889 https://doi.org/10.1002/asna.200710776 CrossRef Search ADS   Kouwenhoven M. B. N., Brown A. G. A., Portegies Zwart S. F., Kaper L., 2007, A&A , 474, 77 CrossRef Search ADS   Kraus A. L., Cody A. M., Covey K. R., Rizzuto A. C., Mann A. W., Ireland M. J., 2015, ApJ , 807, 3 https://doi.org/10.1088/0004-637X/807/1/3 CrossRef Search ADS   Kroupa P., Petr M. G., McCaughrean M. J., 1999, New Astron. , 4, 495 https://doi.org/10.1016/S1384-1076(99)00038-X CrossRef Search ADS   Kroupa P., Aarseth S., Hurley J., 2001, MNRAS , 321, 699 https://doi.org/10.1046/j.1365-8711.2001.04050.x CrossRef Search ADS   Lada C. J., Lada E. A., 2003, ARA&A , 41, 57 CrossRef Search ADS   Lada C. J., Margulis M., Dearborn D., 1984, ApJ , 285, 141 https://doi.org/10.1086/162485 CrossRef Search ADS   Larson R. B., 1981, MNRAS , 194, 809 https://doi.org/10.1093/mnras/194.4.809 CrossRef Search ADS   Lepine J. R. D., Duvert G., 1994, A&A , 286, 60 Lindblad P. O., Grape K., Sandqvist A., Schober J., 1973, A&A , 24, 309 Lindegren L. et al.  , 2016, A&A , 595, A4 CrossRef Search ADS   Lodieu N., 2013, MNRAS , 431, 3222 https://doi.org/10.1093/mnras/stt402 CrossRef Search ADS   Madsen S., Dravins D., Lindegren L., 2002, A&A , 381, 446 CrossRef Search ADS   Makarov V. V., 2007, ApJS , 169, 105 https://doi.org/10.1086/509887 CrossRef Search ADS   Mamajek E. E., Bell C. P. M., 2014, MNRAS , 445, 2169 https://doi.org/10.1093/mnras/stu1894 CrossRef Search ADS   Mamajek E. E., Meyer M. R., Liebert J., 2002, AJ , 124, 1670 https://doi.org/10.1086/341952 CrossRef Search ADS   Marks M., Kroupa P., 2012, A&A , 543, A8 CrossRef Search ADS   Maschberger T., 2013, MNRAS , 429, 1725 https://doi.org/10.1093/mnras/sts479 CrossRef Search ADS   Mel'nik A. M., Dambis A. K., 2017, MNRAS , 472, 3887 CrossRef Search ADS   Michalik D., Lindegren L., Hobbs D., 2015, A&A , 574, A115 CrossRef Search ADS   Moran P. A. P., 1950, Biometrika , 37, 17 https://doi.org/10.1093/biomet/37.1-2.17 CrossRef Search ADS PubMed  Murphy S. J., Lawson W. A., Bento J., 2015, MNRAS , 453, 2220 O'dell C. R., Wen Z., 1994, ApJ , 436, 194 https://doi.org/10.1086/174892 CrossRef Search ADS   Odenkirchen M., Grebel E. K., Dehnen W., Rix H.-W., Cudworth K. M., 2002, AJ , 124, 1497 https://doi.org/10.1086/342287 CrossRef Search ADS   Olczak C., Pfalzner S., Eckart A., 2008, A&A , 488, 191 CrossRef Search ADS   Parker R. J., Goodwin S. P., 2012, MNRAS , 424, 272 https://doi.org/10.1111/j.1365-2966.2012.21190.x CrossRef Search ADS   Parker R. J., Quanz S. P., 2012, MNRAS , 419, 2448 https://doi.org/10.1111/j.1365-2966.2011.19911.x CrossRef Search ADS   Parker R. J., Wright N. J., Goodwin S. P., Meyer M. R., 2014, MNRAS , 438, 620 https://doi.org/10.1093/mnras/stt2231 CrossRef Search ADS   Pecaut M. J., Mamajek E. E., 2016, MNRAS , 461, 794 https://doi.org/10.1093/mnras/stw1300 CrossRef Search ADS   Pecaut M. J., Mamajek E. E., Bubar E. J., 2012, ApJ , 746, 154 https://doi.org/10.1088/0004-637X/746/2/154 CrossRef Search ADS   Pfalzner S., 2009, A&A , 498, L37 CrossRef Search ADS   Pfalzner S., Kaczmarek T., 2013, A&A , 559, A38 CrossRef Search ADS   Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A , 48, 431 CrossRef Search ADS   Preibisch T. Mamajek E., 2008, The Nearest OB Association: Scorpius-Centaurus (Sco OB2), Handbook of Star Forming Regions , p. 235 Preibisch T., Zinnecker H., 1999, AJ , 117, 2381 https://doi.org/10.1086/300842 CrossRef Search ADS   Rizzuto A. C., Ireland M. J., Robertson J. G., 2011, MNRAS , 416, 3108 https://doi.org/10.1111/j.1365-2966.2011.19256.x CrossRef Search ADS   Rosotti G. P., Dale J. E., de Juan Ovelar M., Hubber D. A., Kruijssen J. M. D., Ercolano B., Walch S., 2014, MNRAS , 441, 2094 https://doi.org/10.1093/mnras/stu679 CrossRef Search ADS   Sartori M. J., Lépine J. R. D., Dias W. S., 2003, A&A , 404, 913 CrossRef Search ADS   Scally A., Clarke C., 2001, MNRAS , 325, 449 https://doi.org/10.1046/j.1365-8711.2001.04274.x CrossRef Search ADS   Scally A., Clarke C., 2002, MNRAS , 334, 156 https://doi.org/10.1046/j.1365-8711.2002.05503.x CrossRef Search ADS   Schönrich R., Binney J., Dehnen W., 2010, MNRAS , 403, 1829 https://doi.org/10.1111/j.1365-2966.2010.16253.x CrossRef Search ADS   Taylor M. B., 2005, in Shopbell P. Britton M. Ebert R., eds., ASP Conf. Ser., Vol. 347, Astronomical Data Analysis Software and Systems XIV , p. 29 Tobin J. J., Hartmann L., Fűrész G., Hsu W.-H., Mateo M., 2015, AJ , 149, 119 https://doi.org/10.1088/0004-6256/149/4/119 CrossRef Search ADS   Tutukov A. V., 1978, A&A , 70, 57 van Leeuwen F., 2007, A&A , 474, 653 CrossRef Search ADS   Ward J. L., Kruijssen J. M. D., 2018, MNRAS , preprint (arXiv:1801.03938) Wright N. J., Drake J. J., Drew J. E., Guarcello M. G., Gutermuth R. A., Hora J. L., Kraemer K. E., 2012, ApJ , 746, L21 https://doi.org/10.1088/2041-8205/746/2/L21 CrossRef Search ADS   Wright N. J., Parker R. J., Goodwin S. P., Drake J. J., 2014, MNRAS , 438, 639 https://doi.org/10.1093/mnras/stt2232 CrossRef Search ADS   Wright N. J., Bouy H., Drew J. E., Sarro L. M., Bertin E., Cuillandre J.-C., Barrado D., 2016, MNRAS , 460, 2593 https://doi.org/10.1093/mnras/stw1148 CrossRef Search ADS   © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations