NGC 6231: Cluster Structure

The Structure of the Young Star Cluster NGC 6231. II. Structure, Formation, and Fate

[ Millennium Institute of Astrophysics, Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile Instituto de Fisica y Astronomía, Universidad de Valparaíso, Gran Bretaña 1111, Playa Ancha, Valparaíso, Chile Konstantin V. Getman Department of Astronomy & Astrophysics, 525 Davey Laboratory, Pennsylvania State University, University Park, PA 16802, USA Eric D. Feigelson Department of Astronomy & Astrophysics, 525 Davey Laboratory, Pennsylvania State University, University Park, PA 16802, USA Millennium Institute of Astrophysics, Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile Alison Sills Department of Physics & Astronomy, McMaster University, 1280 Main Street West, Hamilton, ON L8S 4M1, Canada Mariusz Gromadzki Warsaw University Astronomical Observatory, Al. Ujazdowskie 4, 00-478 Warszawa, Poland Millennium Institute of Astrophysics, Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile Instituto de Fisica y Astronomía, Universidad de Valparaíso, Gran Bretaña 1111, Playa Ancha, Valparaíso, Chile Nicolás Medina Millennium Institute of Astrophysics, Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile Instituto de Fisica y Astronomía, Universidad de Valparaíso, Gran Bretaña 1111, Playa Ancha, Valparaíso, Chile Jordanka Borissova Millennium Institute of Astrophysics, Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile Instituto de Fisica y Astronomía, Universidad de Valparaíso, Gran Bretaña 1111, Playa Ancha, Valparaíso, Chile Radostin Kurtev Millennium Institute of Astrophysics, Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile Instituto de Fisica y Astronomía, Universidad de Valparaíso, Gran Bretaña 1111, Playa Ancha, Valparaíso, Chile
June 11, 2017October 3, 2017
June 11, 2017October 3, 2017

The young cluster NGC 6231 (stellar ages 2–7 Myr) is observed shortly after star-formation activity has ceased. Using the catalog of 2148 probable cluster members obtained from Chandra, VVV, and optical surveys (Paper I), we examine the cluster’s spatial structure and dynamical state. The spatial distribution of stars is remarkably well fit by an isothermal sphere with moderate elongation, while other commonly used models like Plummer spheres, multivariate normal distributions, or power-law models are poor fits. The cluster has a core radius of  pc and a central density of 200 stars pc. The distribution of stars is mildly mass segregated. However, there is no radial stratification of the stars by age. Although most of the stars belong to a single cluster, a small subcluster of stars is found superimposed on the main cluster, and there are clumpy non-isotropic distributions of stars outside 4 core radii. When the size, mass, and age of NGC 6231 are compared to other young star clusters and subclusters in nearby active star-forming regions, it lies at the high-mass end of the distribution but along the same trend line. This could result from similar formation processes, possibly hierarchical cluster assembly. We argue that NGC 6231 has expanded from its initial size but that it remains gravitationally bound.

stars: kinematics and dynamics; stars: massive; stars: pre-main sequence; stars: formation; open clusters and associations: individual (NGC 6231); X-rays: stars
journal: the Astronomical Journal\correspondingauthor

Michael A. Kuhn

0000-0002-0631-7514]Michael A. Kuhn

1 Introduction

NGC 6231 is a young star cluster (2–7 Myr) at the center of the Sco OB1 association, on the near side of the Sagittarius spiral arm ( kpc). This complex is above the Galactic Plane, projected in front of the Southern Bar of the Milky Way Galaxy, so the line of sight to the cluster is very complex, with numerous field stars. The basic geometry of NGC 6231 and its environs is shown in the mid-infrared image from the WISE all-sky data release (Cutri et al., 2012) presented in Figure 1. NGC 6231 has a substantial population of O-type stars, which power the H ii region Gum 55, covering several square degrees on the sky as shown on the figure. The cluster is larger than the field of view of the Chandra X-ray Observatory (CXO) (outlined in green) that we discuss in this paper. However, like many other very young clusters (cf. the MYStIX study; Feigelson et al., 2013; Feigelson, 2017), a significant fraction of cluster members are concentrated in a dense central region that is the focus of our investigation – the cluster core is shown by the yellow ellipse.

The molecular cloud from which NGC 6231 formed has dissipated, but molecular clouds still exist around the periphery of the H ii region, including the Large Elephant Trunk to the north-west and IC 4628 to the north-east (visible in the WISE image). The lack of molecular cloud material implies that star-formation has ceased in NGC 6231 itself, but these nearby regions may be sites of ongoing star-formation, possibly triggered by NGC 6231 (Reipurth, 2008). In addition to the cluster members of NGC 6231, O, B, and pre-main-sequence stars are distributed throughout the Sco OB1 association, many of which are part of a loose subcluster, Tr 24, centered north-east of NGC 6231 (Perry et al., 1991).

Figure 1: WISE mosaic of NGC 6231 and its surroundings in the 3.4 m (blue), 12 m (green), and 22 m (red) bands, with a logarithmic color scale. The X-ray field of view is shown by the green polygon, and a yellow ellipse marks the cluster core region measured in this paper (labeled ). The Gum 55 H ii region is outlined by the dashed white line, and several other nebulae and star clusters associated with Sco OB1 are marked. The coordinate grid shows Galactic coordinates.

The large young stellar population in NGC 6231 makes it an ideal testbed for early cluster evolution at a time just after star formation has finished. The star cluster has lost its gas, placing it at a critical stage in its evolution, where it may either disperse as an unbound association or remain as a bound open cluster. The molecular clouds that give rise to clusters like NGC 6231 are depleted both through conversion of gas to stars and by dispersal of the cloud, with dispersal accounting for most of the cloud material (Lada & Lada, 2003). There are a variety of reasons clouds are dispersed, including ultraviolet radiation pressure (Dale & Bonnell, 2008; Dale et al., 2013) and stellar winds from O stars (Townsley et al., 2011), outflows from low-mass stars (Li & Nakamura, 2006), and supernovae (Dekel & Krumholz, 2013).

NGC 6231 has likely had at least one supernova explosion occur in the cluster 3 Myr ago, producing the run-away high-mass X-ray binary (HMXB) HD 153919 4 from the cluster (Ankay et al., 2001). It is probable that the progenitor of this supernova was one of the first O stars formed in NGC 6231 because the cluster has 20 main-sequence members with masses above the limit for supernova explosions and only one Wolf-Rayet star in the field of view (Kuhn et al., 2017, and references therein).

Stars born in massive star-forming complexes are often initially gravitationally bound to the complex (e.g., Jeffries et al., 2014; Kuhn et al., 2015a; Mapelli et al., 2015; Rigliaco et al., 2016). However, gravitationally bound groups of stars can be disrupted by cloud dispersal (Tutukov, 1978; Hills, 1980; Lada et al., 1984) or by tidal interactions with external giant molecular clouds (Kruijssen, 2012). If groups of stars are formed with sufficiently high star-formation efficiencies and are sufficiently massive, they may survive as open clusters. Nevertheless, most groups will disperse as unbound associations (Lada & Lada, 2003), and even surviving open clusters may lose a large fraction of their stars (Kroupa et al., 2001). Investigation into cluster survival has focused on both star-formation efficiency and cluster structure (e.g., Kroupa et al., 2001; Bastian & Goodwin, 2006; Goodwin & Bastian, 2006; Pfalzner, 2009, 2011; Pfalzner & Kaczmarek, 2013a, b; Pfalzner et al., 2014, 2015; Kruijssen, 2012; Gregorio-Hetem et al., 2015; Banerjee & Kroupa, 2017).

It has been previously hypothesized that NGC 6231 is in the process of evolving into an unbound association (e.g., Saurin et al., 2015). We will use the structural properties of NGC 6231, revealed by a more-complete census of its cluster members, to help to address the fate of this young star cluster.

This paper is the second in a two-paper investigation of NGC 6231 using observations from the Chandra X-ray Observatory (CXO) and the VISTA Variables in the Vía Lactéa (VVV) survey (Minniti et al., 2010). The first paper (viz., Kuhn et al., 2017, henceforth Paper I) obtains a new census of the stellar population while this paper addresses cluster structure. Section 2 describes the membership catalogs from Paper I, and Section 3 provide estimates of total populations from these incomplete catalogs. Section 4 discusses the spatial distribution of cluster members, and Section 5 models their surface density distribution. Section 6 studies mass segregation in NGC 6231. Section 7 investigates gradients in stellar ages. Section 8 examines the spatial distribution of stars outside the Chandra field of view. Section 9 discusses the implications of the observed cluster properties on the cluster’s formation and fate. And, Section 10 provides the conclusions.

Many of the data reduction and analysis methods used in this study were developed for the Massive Young Star-Forming Complex Study in Infrared and X-ray (MYStIX) project (Feigelson et al., 2013, and references therein) – a comparative study of 20 young star clusters in nearby massive star-forming regions. These include methods for identifying and modeling subclusters of stars (Kuhn et al., 2014), methods for analyzing the intrinsic densities of stars in clusters (Kuhn et al., 2015b), methods for analyzing mass segregation (M. A. Kuhn et al., in preparation), and methods for identifying gradients in stellar ages (Getman et al., 2014a, b). Many alternative methods for these types of analysis are found in the literature. However, by using MYStIX methods, it is easier to make comparisons between different young star clusters observed with Chandra.

For this study we adopt a distance of 1.59 pc for the cluster and a median age of 3.2 Myr for the stellar population, but with star-formation activity going back at least 6.4 Myr (Paper I).

2 Catalogs of Probable Cluster Members

The study of cluster structure is largely based on the catalog of 2148 probable cluster members in the Chandra field of view from Paper I. This CXOVVV catalog includes X-ray sources, classified based on X-ray properties and optical/near-infrared counterparts, spectroscopic OB stars obtained from the literature, and near-infrared variables from the VVV survey. The initial X-ray catalog is estimated to include 13030 unrelated fields stars. Likely field stars were filtered out on optical color-magnitude diagrams if they appeared too far above or below the locus of cluster members; however, some may remain in the final sample.

Most of the sources in the CXOVVV catalog are low-mass pre-main-sequence stars (Paper I). The known high-mass stellar content includes 1 Wolf-Rayet star, 13 O-type stars, and 82 B-type stars. In Paper I, stellar masses are estimated using near-infrared photometry assuming Siess et al. (2000) pre-main-sequence evolutionary models or using early-type stars’ spectral classifications. Kuhn et al. (2010) found that a similar mass-estimation method for pre-main-sequence stars in Taurus yielded typical errors of 0.15–0.30 dex. Stellar age estimates are obtained with two methods: the vs. color-magnitude diagram () and relations between X-ray luminosity and -band luminosity (; Getman et al., 2014b).

2.1 Outside the Chandra Field of View

Methods for selection of cluster members with only optical/near-infrared data, available in the region outside the Chandra field of view, yield samples with much higher rates of non-member contamination. However, even these catalogs may be useful for visualizing spatial clustering on larger spatial scales. We use both variability and color-magnitude diagrams for selection.

A total of 295 near-infrared variables are found in a box surrounding NGC 6231 (VVV tiles “d148” and “d110”). Near-infrared variability in pre-main-sequence stars may be produced by star spots, accretion from a circumstellar disk, and variable extinction from the circumstellar disk (e.g., Joy, 1945; Herbst et al., 1994). In a representative study of high-amplitude variables in the VVV survey by Contreras Peña et al. (2017a, b), it was estimated that approximately 50% of near-infrared variables were pre-main-sequence stars.

Damiani et al. (2016) note that candidate O–B and A–F stars can be selected using the optical photometry alone, and that many of these stars will be cluster members. We use photometry from the VST Photometric H Survey of the Southern Galactic Plane and Bulge (VPHAS+; Drew et al., 2014) to identify candidate stars with spectral types of F or earlier in the portion of the Sco OB1 survey covered by VPHAS+. The magnitude and color cuts,  mag and 1.5 mag, are designed to approximate the selection rules from Damiani et al. (2016), which used the Johnson-Cousins and bands instead. From the spatial distribution of these stars (§8.3), it is clear that there is also a large unclustered population, which are likely to be field stars.

Figure 2: Maps of surface density [stars pc] of X-ray selected stars in NGC 6231 (center) compared to the Orion Nebula Cluster (left) and NGC 1893 (right). All three regions are shown to the same physical scale (a 5-pc length scale is shown) and have been corrected for differences in X-ray sensitivity (a color-bar shows densities scaled to the full IMF). The fields are oriented with north up and east to the left.

3 Intrinsic Stellar Population

Even with the deep Chandra exposure, most stellar-mass cluster members are not detected. The X-ray catalog is only complete for sources with X-ray photon fluxes greater than [photon s cm] in the 0.5–8.0 keV band. For pre-main-sequence stars, there is a positive correlation between stellar mass and X-ray luminosity (e.g., Telleschi et al., 2007), so higher-mass pre-main-sequence stars are more likely to be detected while lower-mass pre-main-sequence stars are less likely to be detected.

In Paper I, we use both the initial mass function (IMF) and X-ray Luminosity Function (XLF) to extrapolate the number of stars missed in the observations (e.g., Kuhn et al., 2015b). The IMF method from Paper I gives a total of 5700250 stars projected within the Chandra ACIS-I field of view, using the masses calculated with a 3.2-Myr isochrone. However, if an older age of 6.4 Myr were assumed, the estimated number of stars would increase to 7500360. The XLF analysis gives 6000530 stars assuming that the change in shape of the XLF is small during the first 5 Myr. The uncertainties reported above reflect statistical errors alone. These are calculated using bootstrap resampling with replacement (1000 draws), where normally distributed random errors were added to or for each draw – a factor of 2 for stellar mass or the reported uncertainty for X-ray luminosity.

Systematic uncertainties due to stellar-mass estimation and assumptions about the IMF and XLF are difficult to quantify and may be larger than the statistical uncertainties. For example, Kuhn et al. (2015b, their Figure 4) compared IMF and XLF estimates for the MYStIX subclusters and found discrepancies of 0.25 dex. This may be regarded as a pessimistic estimate of the systematic uncertainty on number of stars in NGC 6231. However, the estimates for NGC 6231 may be more accurate than for MYStIX subclusters because of NGC 6231’s higher number of observed stars, its lower extinction and lack of strong variation in extinction, and its lack of infrared nebulosity. The IMF and XLF estimates of 5700250 and 6000530 stars in NGC 6231 agree within their estimated statistical uncertainties.

For an analysis of cluster structure, it is essential that inhomogeneities in detector sensitivity not affect star counts. Chandra’s sensitivity is greatest near the optical axis of the telescope, leading to a larger number of faint sources being detected near the cluster center (Broos et al., 2011). Thus, we remove sources with X-ray fluxes below the photon-flux completeness limit from the study, as has been done for previous work (e.g., Feigelson et al., 2011; Kuhn et al., 2014). This leaves 826 X-ray sources, combined with all known O and B-type stars, for a total sample of 885 objects out of a total of 2148 probable cluster members. If we assume that the total number of stars projected in the field of view is 5700, a correction factor of 6.5 would need to be applied to any star counts or surface densities measured from this sample to obtain an intrinsic astrophysical value. (The correction factor would be 8.5 if the older age were assumed.)

The 50% mass completeness limit found by the IMF analysis is 0.5  (Paper I). Source detection probability rolls off gradually as a function of stellar mass, due to the scatter in the relation, so the limit we report for mass is where a pre-main-sequence stars has a 50% chance of being detected; but most pre-main-sequence stars and OB stars above this limit will be detected. Late-B and A stars are expected to be missing from X-ray surveys; however, some of these stars do have X-ray emitting pre-main-sequence companions (Paper I).

4 Cluster Morphology

The projected surface density of stars in NGC 6231, normalized to the full 5700-member population, is displayed in Figure 2 (center panel), along with similarly normalized stellar surface-density maps of two MYStIX star-forming regions, the Orion Nebula Cluster (left panel) and NGC 1893 (right panel), obtained by Kuhn et al. (2015b). These maps were generated by adaptively smoothing the spatial point patterns of stars using the algorithm adaptive.smoothing from the R software package spatstat (Baddeley et al., 2005b, 2015). This algorithm subdivides the field of view using the Voronoi tessellation, using a fraction () of points to generate the tiles and a fraction () of points to estimate surface density in each tile. This procedure is repeated 500 times, and the results are averaged to produce the smoothed maps. The values in the surface-density maps are then multiplied by correction factors (§3) to estimate intrinsic surface density.111Use of the correction factor requires the assumption that X-ray properties of stars are not correlated with position. This assumption does not necessarily hold because the X-ray photon flux is correlated with stellar mass, and the cluster is likely to be mass segregated. Nevertheless, the star-counts are dominated by lower mass stars which are not strongly mass segregated. A complete description of how the maps were generated for the two MYStIX regions is given by Kuhn et al. (2015b).

The three clusters in Figure 2 are shown using the same spatial scale (a 5 kpc arrow is shown) and same surface-density scale (indicated by the color bar) to allow their morphologies to be directly compared. Although all three regions have similar ages and numbers of stars (approximately 2.5 Myr and 2600 stars for the Orion Nebula, 3.2 Myr and 5700 stars for NGC 6231, and 2.6 Myr and 4600 stars for NGC 1893), there are significant differences in cluster structure. The Orion Nebula Cluster is surrounded by a dense molecular cloud, but the center of the cluster is partially evacuated; in contrast, the bubbles around NGC 6231 and NGC 1893 are much larger. NGC 6231 is much less dense than the Orion Nebula Cluster, with a central surface density of 300 stars pc compared to 10,000 stars pc, and has a larger physical size. NGC 6231 is similar in size to NGC 1893; however, NGC 1893 is significantly more clumpy and elongated, being divided into 10 subclusters. The individual subclusters in NGC 1893 also appear physically smaller than the NGC 6231 cluster, although they have similar peak surface densities.

When compared to 17 star-forming regions from MYStIX, shown by Kuhn et al. (2015b, their Figure 5), NGC 6231 appears atypical in having a simple structure with a much larger size than the typical subcluster, rather than a collection of denser subclusters. This suggests that NGC 6231 has expanded considerably from its initial size. NGC 6231 may be most similar to NGC 2244, an expanded cluster of stars that is part of the Rosette Nebula star-forming region. However, NGC 2244 in Rosette does not show mass segregation (Wang et al., 2008) whereas NGC 6231 does (§6).

5 Surface Density Models

Several families of spherically symmetric models have been used to fit the (surface) density profiles of star clusters. Models include the isothermal sphere, the King profile, and the Plummer Sphere, which all represent approximations to density profiles of virilized, gravitationally bound groups of stars in a quasi-equilibrium state (Binney & Tremaine, 2008). None of the cluster profiles would necessarily be expected to provide a good model for young clusters ( Myr), which are not expected to be in dynamical equilibrium due to their young age and would show an imprint from formation in their natal molecular clouds. The lack of kinematic data for cluster members means that we cannot assess the thermodynamic state of the star cluster, so we use these models as empirical descriptions of the surface density.

Other possible models include multivariate normal distributions and a power-law radial-density gradients. If stars have a Maxwell-Boltzmann velocity distribution and are allowed to freely expand from a point, then their resulting spatial distribution would be a multivariate normal distribution. A radial power-law surface-density distribution has also been proposed as a possible model for young star clusters by Cartwright & Whitworth (2004). The power-law model is scale invariant, unlike models that require a critical length scale .

We fit the projected spatial distribution of stars using several functional forms based on the above models.


In these equations is the projected distance from the center of the cluster, is the surface density at the center of the cluster, is a characteristic radius called the “core radius” in the case of the Hubble model and the Plummer sphere, is a scaling constant, and is a power-law index. The isothermal sphere has both singular (density diverges at ) and non-singular solutions. The non-singular solution can be approximated out to several core radii by the Hubble model (Binney & Tremaine, 2008) above. However, for , asymptotically approaches a power-law dependence on , with an index of , not as suggested by the Hubble model. The singular isothermal sphere is a power-law model with .

Often the distributions of stars show elongation; this includes the Orion Nebula Cluster (Hillenbrand & Hartmann, 1998) and many of the subclusters in the MYStIX star-forming regions (Kuhn et al., 2014). To account for elongation of subclusters, Kuhn et al. (2014) introduced additional parameters to the model for ellipticity, , and orientation, , in projection on the sky. Equations 14 can be redefined to describe ellipsoidal distributions if we redefine as


where and are the distances of a point from the cluster center along - and -axes. Following Kuhn et al. (2014), we refer to the distribution described by the transformed version of Equation 1 as the “isothermal ellipsoid,” although this is an empirical model rather than one derived from physics.

We use the maximum likelihood method to fit these models to the spatial distributions of cluster members, using software provided by Kuhn et al. (2014, their Appendix A) written in the R programming language. The log-likelihood is


where is the set of points, are the model parameters, and is the field of view. The Nelder-Mead algorithm is used to optimize , while the Hessian matrix of at the maximum likelihood parameters is used to estimate uncertainty in model parameters.

The resulting cluster parameters include coordinates of the cluster center, the core radius , the ellipticity and orientation of the cluster, and the central density for the isothermal ellipsoid, the Plummer ellipsoid, and the multivariate normal distribution. For the singular isothermal sphere and power-law model, the parameter is not needed, and for the power-law model the index is also fit.

5.1 Single-Cluster Models

The distribution of stars in NGC 6231 was fit with the isothermal ellipsoid, the Plummer ellipsoid, the multivariate normal, the singular isothermal ellipsoid, and the power-law models (Figure 3). We omit the King model, a modification of the isothermal sphere with an outer truncation radius, because the cluster is truncated by the field of view. In the figure, the panels in the left column show green ellipses, representing model parameters, overplotted on spatial maps of fit residuals. These ellipses indicate where the cluster is centered, the ellipticity and orientation of the model, and, for models with a characteristic radius, the size of the cluster core. The panels in the right column show a one-dimensional (constant declination) slice through the best-fit parametric models (gray lines) and the adaptively smoothed data (black lines). The -axes are plotted with logarithmic values to allow a high contrast in surface density to be shown.

Goodnesses of fit can be evaluated using the residual maps (Figure 3, left column). The residual maps are produced by the function diagnose.ppm in spatstat. This tool smooths both the model and the points with a kernel ( pc), and the residual () map is shown by the color scale in units of stars pc. The mathematical foundations for this analysis are given by Baddeley et al. (2005a, 2008). In these residual maps, a good fit is indicated by a residual value of 0 (white), while positive and negative values (red or blue) indicated differences between the data and the model. Thus, patches of dark red may represent subclusters of stars not accounted for by the models. For the singular isothermal ellipsoid and the power-law models, which diverge at , the number of stars remains finite when , so the kernel will smooth over the divergence at the center.

A description of the best fits is given below.

Isothermal ellipsoid

There is close agreement between the model and the adaptively smoothed data in both the cluster center and in the outer region of the cluster, out to  pc (most of the field of view). In particular, in the cluster wings, the slope of the model is a good match to the smoothed data. Residuals are small (25 residual stars pc which is 10% of the peak surface density) except for a peak of 120 residual stars pc to the north-west of the cluster center. The deviation in the outer part of the cluster may be partially due to the inadequacies of the Hubble model as an approximation for the isothermal sphere, but the numbers of stars in these regions are small, so the differences could also be due to stochastic fluctuations in counting statistics.

Plummer ellipsoid

There is a significant difference between the model and the adaptively smoothed data at the cluster center. This arises because the Plummer model has a flatter core region than the isothermal ellipsoid, as well as having a steeper decrease in surface density with radius outside the cluster core. The underestimate in number of stars in the cluster center can be seen as a positive residual at this location in the smoothed map.

Multivariate normal

The central density is significantly underestimated, and the rapidly decreasing wings of the normal distribution do not match the shape of the smoothed data, which is more gradually decreasing. The 1-sigma ellipse of the distribution is significantly larger than the core radius found by the isothermal ellipsoid and the Plummer ellipsoid models because the multivariate normal distribution must be enlarged to fit the broad wings of the distribution of stars.

Singular isothermal ellipsoid

The cusp at the center of this model strongly overestimates the number of stars at the cluster center. Furthermore, the slope (given by ) is too flat in the outer regions of the model producing a strong negative residual surrounded by a ring of positive residuals.

Power-law model

Allowing to be a free parameter decreases the number of stars in the central cusp, decreasing the negative residual. Nevertheless, the power-law index of the model () is even flatter in the outer regions of the cluster than the singular isothermal ellipsoid model, poorly matching the adaptively smoothed data.

Figure 3: Residual maps (left column) and surface density profiles (right column) for the 5 cluster models. Left: Residuals were calculated using a Gaussian kernel with  pc. Negative residuals are blue and positive residuals are red. Possible subclusters appear as deep red spots. The green ellipses indicate model ellipticity and orientation, with solid lines indicating core radius and dashed lines showing a representative radius. Right: Surface density along a constant constant declination (5000) for the adaptively smoothed data (black line) and the model (gray line).

Thus, the isothermal ellipsoid model provides a remarkably good empirical approximation of the observed surface density distribution, while all other models are clearly inadequate. This is good motivation for the use of the isothermal ellipsoid form for analysis of NGC 6231. The good fit found using the isothermal ellipsoid model is consistent with the results from the MYStIX clusters RCW 38, the Flame Nebula Cluster, and M 17, even though these clusters are much denser (Kuhn et al., 2014). The physical cluster parameters based on the isothermal ellipsoid model are provided in Table 1, labeled “Model 1.”

Whether or not young star clusters have a critical length scale has been an open question (e.g., Elmegreen & Elmegreen, 2001; Cartwright & Whitworth, 2004). For NGC 6231, the core radius for the isothermal ellipsoid model appears to be such a length scale. The core radius, , is inconsistent with a value of 0. In addition, both scale invariant models—the singular isothermal ellipsoid or the power-law model—poorly fit the data.

When the distribution of stars is fit with only one cluster component, the residual maps show an overdensity of stars north west of the cluster center, at coordinates 165400.6 414807. This residual represents anisotropic structure that could indicate another subcluster. We investigate this possibility quantitatively in Section 5.2.

For the model fitting presented above, all stars of different masses are all treated equally. To investigate whether the combination of stars of different masses affects the functional form of the models, we redo model fitting for 768 low-mass stars, excluding stars with  . The results are nearly the same as when using the full mass range. The isothermal ellipsoid provides the best fit to both the cluster core and cluster wings, while the normal distribution underestimates the number of stars in the core and the scale-invariant models both yield results that are too cuspy in the center and have slopes that are too flat in the outer regions. However, for the smaller sample, it is more difficult to distinguish between the Plummer ellipsoid and the isothermal ellipsoid. A residual corresponding to the candidate subcluster is still apparent for all models.

5.2 Two-Subcluster Model

Beyond the single-cluster model, it is possible to model multiple cluster components using a statistical method known as mixture models (see review by Kuhn & Feigelson, 2017). A mixture model is a probabilistic model in which the probability density function (e.g., surface density) for a set of points is composed of the sum of multiple probability density functions for subpopulations (e.g. subclusters of stars).

For NGC 6231, each component has the form of an isothermal ellipsoid model, which may be used to model different groups of stars in the region. This method of cluster analysis was used by Kuhn et al. (2014) to identify subclusters in MYStIX. For NGC 6231, two subclusters are suspected: the main cluster and the possible subcluster to the north-west – the smaller subcluster is designated Subcl. A following the notation of Kuhn et al. (2014).

Determining whether a population of stars is one or more clusters can be viewed as a model selection problem. Penalized likelihood methods are commonly used for mixture model problems, including the Akaike informaion criterion (AIC; Akaike, 1974) and Bayesian information criterion (BIC; Schwarz et al., 1978), defined by the formulas


where is the log-likelihood, is the number of parameters in the model, and is the number of points. For the isothermal ellipsoid clusters, is equal to 6 times the number of clusters, and we select the that minimizes the AIC or BIC. The BIC generally favors simpler models than the AIC due to its larger penalty for the inclusion of additional parameters. The choice of AIC, BIC, or other model selection approaches has been widely debated (Lahiri, 2001; Burham & Anderson, 2002; Konishi & Kitagawa, 2008).

The best-fit two-component model includes a main cluster, with properties similar to the single-cluster model, and a small subcluster, coincident with the over-density of stars to the north-east of the main cluster. The log-likelihood increases from 1599 to 1618 with the addition of the subcluster (log likelihood will always be greater for the model with more components). Both the BIC and AIC favor the model with two components. The difference in BIC values, , is far above the threshold for confident preference of the 2-cluster model (Kass & Raftery, 1995). Integrated over the entire field of view, the ratio of the number of stars in the main cluster to the number of stars in the subcluster is 24:1.

The cluster and subcluster are shown in Figure 4, where the elliptical contours marking the cores for each component are plotted on an adaptively smoothed surface density map and on a map of smoothed residuals. The addition of the second cluster significantly decreases the amplitude of the residuals. The core radius of the subcluster is clearly much smaller than that of the main cluster (although not too dissimilar from many subclusters found in MYStIX star-forming regions).

Figure 4: Left: Surface density estimate for stars in NGC 6231 with the cluster-core contours overplotted (black ellipses). The two ellipses correspond to the two-component model in Table 1. Right: The residual plot for the two-component model. The color-scale is the same as in Figure 3 (bottom row), but the positive residual to the north-west of the main cluster is gone.

5.3 Properties of the Main Cluster and Subcluster

Table 1 presents the best-fit models for the single-component model and the two-components model. This includes parameters and uncertainties estimated directly from model fitting, including component centers, core radius, ellipticity, and orientation. Numbers of stars and central surface and volume density also include corrections for incompleteness. For the Hubble model, the central volume density, is related to the central surface density, , by the relation . This equation will also apply to the ellipsoidal model if we assume that the three-dimensional ellipsoid has a harmonic-mean core radius of .

The mixture model provides precise and reliable celestial coordinates of the centers of the main cluster and subcluster (Table 1). The modeling finds the centroid of the star counts belonging to each cluster, correcting for truncation of the field of view and overlapping of the clusters. For the rest of the work we will use the coordinates 165415.414959 as the definition of the center of the main cluster. This is offset by 0.5 arcmin to the east of the location of the peak of the adaptively-smoothed surface-density map.222 The difference between the centers is an indication of the slight asymmetry in the distribution of stars. The precise position of the peak density may depend on the smoothing algorithm that is used.

The subcluster is located at 16541.6 414813. Although there is considerable uncertainty in some subcluster properties from the model fit, the number of stars in the subcluster, integrated over the whole field of view,  stars, is less strongly dependent on the core-radius model parameter. The O9.7Ia+O8V system HD 152234 is projected at the location of the subcluster. However, this alignment could be coincidental due to the large number of massive stars in the cluster.

Another anisotropy in the spatial distribution of stars of the main cluster is a shoulder in surface density, 0.5 pc east of the cluster peak. This can be seen in both the smoothed surface-density maps of both Figure 3 (bottom row) and in the surface-density profile in Figure 3 (top row). This asymmetry does not correspond to its own mode in surface density, and it could be merely a statistical fluctuation in the distribution of stars.

5.4 Characterizations of Cluster Size and Total Mass

Estimates of cluster radius and mass must be carefully defined for young star clusters in star-forming complexes like Sco OB1. In theoretical studies of cluster evolution, both the cluster mass and the half-mass radius of a cluster (i.e. the radius of a sphere that encompasses half the mass of the cluster) have been important quantities for characterizing clusters. Neither of these can be directly measured for the isothermal ellipsoid model of NGC 6231 because the cluster is truncated by the field of view. Instead, we report two values that are well defined by our model: , the number of stars (corrected for incompleteness) within a region 4-times larger than the cluster core, and the corresponding radius . In general there is no fixed ratio between the core radius and half-mass radius, but for the clusters in the Portegies Zwart et al. (2010) sample,  dex, so is likely within a factor of several of the half-mass radius. For and , we provide the uncertainty on the number of stars within the reported radii, including only the statistical uncertainties on number of stars. For and , uncertainty from estimation of core radius is also included.

Model 1 Model 2
Main    Main        Subcl. A
R.A. (J2000) 16 54 14.7 [0.5] 16 54 15.9 [0.6] 16 54 1.6 [0.1]
Decl. (J2000) 41 49 47 [11] 41 49 59 [12] 41 48 13 [6]
(pc) 1.20.1 1.20.1 0.0480.031
(pc) 4.60.4 4.60.4 0.190.12
(stars) 1400100 1300100 1710
(stars) 5600250 5400240 7220
(stars) 5700250 5500240 24040
(stars pc) 47080 46060 3400:
(stars pc) 20050 20050 35000:
0.340.05 0.330.05 0.660.46
(degrees) 1594 1625 1659
1599 1618
BIC 3158 3194
AIC 3186 3223

Note. – The best-fit parameters for the isothermal ellipsoids used to model the projected surface density. The results for a one-component model (single cluster) and a two-component model (main cluster subcluster) are shown. Rows 1-2: The coordinates of the ellipsoid centers, with uncertainty given in brackets. Row 3: The harmonic-mean radius of the isodensity ellipse enclosing the cluster core. Row 4: A characteristic radius four times as large as the core radius. Rows 5-7: The number of stars assigned to each component within 1 core radius, within 4 core radii, and within the field of view. Rows 8-9: The surface density and the volume density at the center of the ellipsoids. Row 10: Ellipticity. Row 11: Ellipsoid position angle in degrees east from north. Rows 12-14: Log likelihood, BIC, and AIC of the model. Entries with units of ”stars” have been corrected for sample incompleteness and represent the intrinsic stellar population down to 0.08 . The reported uncertainties represent statistical uncertainty in both estimation of number of stars and model fitting, but exclude possible systematic error discussed in §3.

Table 1: Best-fit Cluster Models

The total number of stars in a cluster described by an isothermal ellipsoid model depends on the outer radius of the cluster. However, the outer radius can be difficult to constrain because cluster members will be most thinly distributed near the outer edge of the cluster. A cluster like NGC 6231 subtends a large area on the sky, outside the fields of view of the various X-ray observations, so large surveys would be needed to determine the outer edge of the cluster. Saurin et al. (2015) use the 2MASS catalog to estimate where the spatial over-density of stars meets the background level, and they report a radius of 36.2 pc. Our investigation in §8.3 shows that these stars are not distributed isotropically around NGC 6231, being mostly distributed to the north and east of the main cluster. With only spatial data, it is impossible to determine whether these stars are part of NGC 6231 or part of other clusters or associations in Sco OB1.

There are several methods to describe the number of stars in NGC 6231 in a way that can be compared to observations of other clusters. The number of stars projected in the Chandra field of view is  stars; the number of stars projected in an ellipse with characteristic radius is  stars; and the number of stars within a three-dimensional ellipsoid with characteristic radius is 4300 stars. These subtle geometric differences in definitions do not strongly affect results—the main point being that, in each case, a large number of cluster members may exist outside the region being considered.

The total mass of stars in the cluster depends on the binary fraction of its low-mass stars, which is not yet well constrained by observation. We follow Maschberger (2013) to estimate the mean mass of single stars and multiple-star systems based on the Chabrier (2003) IMFs. Equation 25 from Maschberger (2013) gives an average mass for the mass range 0.08–150  and an average mass of for systems. When this uncertainty is combined with the statistical uncertainty on total number of stars, the total mass of stars projected within the 4-core radius ellipse, down to the hydrogen-burning limit, is in the range 3300–4200 .

6 Segregation of Stars by Stellar Mass

The spatial distribution of stars, with their masses indicated, is shown in Figure 5 (left). From this figure it can be seen that stars of various masses are mixed together, with both high-mass stars and low-mass stars concentrated toward the center of the cluster. It appears that the O and B stars are relatively more likely to be found in the cluster center compared to low-mass stars. However, several high-mass stars are also located outside the cluster center, including the O9.5 III star HD 152076 (central distance of 13) and the O9 III+O9.7 V system HD 152247 (central distance of 11). Only stars above the mass completeness limit of 0.5  are shown, so as to avoid effects of insensitivity to lower mass stars, either from instrumental effects or from crowding of stars. A mass-complete sample is important for establishing the reality of mass segregation because incompleteness could masquerade as an erroneous signature of mass segregation (Ascenso et al., 2009).

6.1 Statistical Tests for Mass Segregation

Figure 5 (right) shows an estimate of the expected value of as a function of position of a star in the field of view. Interpolations of properties of points to generate a map of expected values is a well-known method of statistical analysis (e.g., Olea, 2000). Due to large differences in surface density, we perform adaptive kernel smoothing, where the kernel bandwidth at a point is set to the distance to the 100th nearest star. The resulting map shows the highest average at the center of the cluster, as would be expected if the cluster were mass segregated. The mean value of this peak is (2.8 ). This peak value is relatively low because of the large population of low-mass stars relative to high-mass stars, even near the cluster center.

The statistical significance of the peak in the map can be judged through Monte Carlo simulations. For these simulations, stellar masses are randomly permuted among the stars to simulate a case in which stellar mass is independent of position, and a surface density map is generated in the same way. One thousand simulations are performed, and the maximum peak in the simulated maps is recorded. Based on the distribution of simulated peak values, the observed value peak of has a -value of 0.01.

Comparison of radial-distance distributions of stars in different mass strata (or comparison of stellar mass distributions in radial bins) has been a common method of testing for mass segregation (Hillenbrand & Hartmann, 1998; de Grijs et al., 2002; Stolte et al., 2006; Kuhn et al., 2010). Figure 6 shows cumulative distributions for three mass strata, the low-mass stars (), intermediate-mass stars (), and high-mass stars (). The radial distributions of stars in these groups are compared using the two-sample Anderson-Darling test (Stephens, 1974). Only the stratum containing the high-mass stars shows any difference in radial distribution (with -value of 0.001 and 0.03 when compared to low- and intermediate-mass strata, respectively), while the low- and intermediate-mass strata have radial distributions that are very similar to each other. This finding agrees with observations of NGC 6231 by Raboud & Mermilliod (1998) that only the most massive stars are segregated, but intermediate stars are well mixed. In contrast, the study of mass segregation in 17 MYStIX star-forming regions by Kuhn et al. (in preparation) reveals many cases where even low-mass stars do appear to be strongly segregated by mass; but this is not the case for NGC 6231.

Figure 5: Left: Stars in flux-complete sample with mass estimates 0.5  are plotted on the Chandra/ACIS-I field of view for NGC 6231. The area of the circles is proportional to the estimated stellar masses, and stars with are plotted with black circles and stars with are plotted with red circles. Right: Adaptively smoothed mean values of . The smoothing uses a Gaussian kernel, with a width equal to the distance to the 10th nearest point. The color scale shows the mean values over a range from 1–5 .
Figure 6: Cumulative distributions of projected distance from the center of the cluster for three different mass strata, which are divided at 1.8  and 7.9  into low- (blue), intermediate- (green), and high-mass (red) strata. The Anderson-Darling two-sample test gives -values of 0.23 (low vs. intermediate), 0.001 (low vs. high) and 0.03 (intermediate vs. high).

6.2 An Empirical Model of Mass Segregation

To further investigate the effect of stellar masses on the distributions of stars, we subdivided the sample by stellar mass and fit the subsamples with the single isothermal-ellipsoid model from Section 5. (The small subcluster to the northwest is ignored here.) Five samples were used, divided at , containing 251, 223, 126, 75, and 41 stars, respectively. Figure 7 shows the plot of subcluster core radius vs. stellar mass. The points have an abscissa value equal to the mean mass of stars in a subsample (the horizontal bars show the range of masses in the subsample) and an ordinate value equal to the core radius (the vertical error bars show the uncertainty on core radius). The gray, dashed lines show relations of for comparison.

Figure 7 shows that, aside from the second mass bin ( ), the core radius decreases monotonically with increasing stellar mass. However, for the range 0.5–3 , the statistical uncertainties on core radius show that core radius is not statistically different for the first three mass bins. Although there is a persistent decrease between 1  and 50 , the difference in core radius only becomes statistically significant for the highest mass bin (stars with  ), which agrees with the results of the Anderson-Darling tests which only show mass segregation for high mass stars. The relation between core radius and mean stellar mass is described by an empirical relation (black line) found using weighted ordinary least-squares regression using the bins shown, where weights are equal to the reciprocal of the estimated uncertainty. Gray dashed lines with a relation are shown for comparison.

Figure 7: The best-fit model-cluster radius is plotted versus stellar mass, based on model fits to 5 sets of stars stratified by stellar mass. The horizontal lines show the range of stellar mass included in each sample, which are divided at , 1, , , and 10 , and the position of the point is the mean mass in each sample. The error bars on the core radius are based on the Hessian matrix at the maximum of the likelihood function, which appear asymmetric on the logarithmic axis scale. Gray dashed lines indicate possible radius–mass relations for a cluster/association with energy equipartition. The solid black line shows an empirical ordinary-least-squares regression to the data.

7 Age Distribution

In Paper I stellar ages are estimated using two independent techniques. The first estimates (denoted ) are based on X-ray and near-infrared photometry, using the method from Getman et al. (2014a, b). The second estimates (denoted ) are based on the optical vs.  color-magnitude diagram. Age estimates from both methods are calibrated using the Siess et al. (2000) models. In addition, the method was based on relations derived for pre-main-sequence stars with ages 5 Myr old, so stars with ages greater than 5 Myr will be assigned 5 Myr as a lower limit. As reported by Paper I, the median for stars in NGC 6231 is 3.2 Myr, which is very similar to the median value of 3.3 Myr. While statistical uncertainties on ages of individual stars may be large, statistical uncertainties on the median ages of sufficiently large sample of stars will be smaller. Thus, median-age estimates can be used to identify spatial age gradients.

In many young star-forming regions, differences in ages of groups of stars are of the order 1 Myr (Getman et al., 2014b). This supports a model of star formation in which all stars do not form at in a monolithic cluster in a single cluster-crossing timescale. Instead, stars form either over multiple free-fall timescales or in multiple, independent subclusters that form at different times. In the two cases from the MYStIX study with sufficient data quality, the Orion Nebula Cluster and NGC 2024, stars within an individual cluster also showed a radial gradient in stellar age, with the youngest stars nearest the cluster center and the older stars in the cluster periphery (Getman et al., 2014a). This was regarded as an unexpected result because stars are typically expected to form first where gas in molecular clouds is densest, which would be in the centers of clusters. Additional cases are reported by K. V. Getman et al. (2017, in preparation) indicating that age gradients are common in young clusters.

The NGC 6231 cluster appears to have an age spread of 2–7 Myr, noted by previous studies (Sana et al., 2007; Sung et al., 2013; Damiani et al., 2016) and Paper I. With a median age of 3.2 Myr, NGC 6231 is older than most MYStIX star-forming regions, including the Orion Nebula Cluster and NGC 2024, and star formation has ended in NGC 6231, so stars cannot be extremely young (e.g., 0.1 Myr). Figure 8 shows the median and values (25%, 75%, and uncertainties on the median are also shown) for stars at several projected distances from the cluster center. The range of ages (median, 1st-quartile, and 3rd-quartile) estimated using both independent techniques are very similar, although the ages have a greater range because they are not limited to 1–5 Myr.

Median ages of the stellar population range from 3 to 4 Myr on a distance baseline ranging from 0 to 4 pc, and no systematic trend is seen. Thus, there is no large-scale radial age gradient. Figure 9 shows an adaptively smoothed map of mean . Here, variations in age throughout the field of view are small, and no global trend is evident. As before, Monte Carlo simulations can be used to determine if structures in a map may be the result of random fluctuations. For the map of mean ages, the low-amplitude structure is consistent with variations for a distribution where stellar age is independent of position. The lack of a radial age gradient implies that, either NGC 6231 never had a radial age gradient, or that a previously existing gradient disappeared with age due to dynamical mixing. This result is clearly different from the younger MYStIX clusters where a radial–age gradient has been observed.

Figure 8: Box-And-Whisker plots of the (left) and (right) distributions for stars in several radial bins. The boxes indicate the 25%, 50% (median), and 75% quartiles for the Age values for stars 0–1 pc, 1–2 pc, 2–3 pc, and 3–4 pc from the cluster center (from the single-ellipsoid model). Notches indicate uncertainty on the sample median. The whiskers are indicate the range of the data, up to 1.5 times the interquartile ratio, and outliers are drawn as open circles. The data on these plots span different ranges because is only estimated for stars with ages 5 Myr and stars that appear older are assigned 5  as a lower limit. Nevertheless, both plots show very similar median ages.
Figure 9: Adaptively smoothed mean ages from the Age method. The smoothing uses a Gaussian kernel, with a width equal to the distance to the 5th nearest point. The color scale shows variations in mean calculated age from 2 to 5 Myr.

8 Larger-Scale Galactic Environment

NGC 6231 extends beyond the Chandra field of view. However, in these regions, spatial distributions of stars can be mapped using the catalog of 295 near-infrared variables from VVV and the candidate early/mid-type stars from VPHAS+. Both these catalogs suffer from more incompleteness and more field-star contaminants than the Chandra-based catalog. Nevertheless, contaminants are expected to be distributed smoothly (except for possible patchy absorption) with some dependence on Galactic latitude, and candidate selection is not strongly affected by position. Thus, clustering of these sources is likely associated with real clusters or associations of stars. These two samples trace different types of populations. The amplitude of variability in the near-infrared decreases with a pre-main-sequence star’s age (Rice et al., 2015), so VVV variables will trace the youngest stellar population, while candidate early/mid-type stars will be less sensitive to age.

8.1 Modeling Clusters of Near-infrared Variables

We use a mixture model approach to identify possible clusterings of variable stars measured in the VVV band. Given that there may be a high number of field stars, one component of the model will account for these objects. Field stars are expected to be smoothly distributed rather than clustered. However, the large field of view means that the projected density of field stars will vary with Galactic line of sight, mostly as a function of Galactic latitude . The contribution of the unclustered field-star component in the mixture model is similar to the use of an unclustered component by Kuhn et al. (2014), but here it is a model with several parameters. We use the flexible model form,


to describe variation in field-star densities, where the variable is the normalization of the model and the polynomial coefficients and are the parameters to be fit.

Figure 10 (left) shows the model residuals when the data have been fit with the unclustered model (the hypothesis that all variables are field stars). The best-fit parameters of the density gradient are  deg and  deg. Several density peaks are not well modeled, leading to high residuals (red patches). The residuals include a peak associated with NGC 6231, several peaks to the south-east of NGC 6231 near the Galactic plane, and several peaks to the north and south of NGC 6231.

Next, we test various mixture models, composed of “isothermal ellipsoid” components and an “unclustered” component. The identification of multiple statistical clusters follows the same method as Kuhn et al. (2014). Models with , 1, 2, 3, 4, 5, and 6 are fit, and the AIC and BIC calculated. For this model, the number of parameters is , so the penalty for each additional component is 12 for the AIC and 34 for the BIC. The AIC values are 1767, 1771, 1717, 1668, 1671, 1668, and, 1666 and the BIC values are 1778, 1804, 1773, 1745, 1770, 1789, and 1810, for , 1, 2, 3, 4, 5, and 6 clusters respectively. Thus, the BIC clearly favors a 3-cluster model, while the AIC is consistent with models with 3–6 clusters. Note that groups of stars are only clusters in a statistical sense, while the physical nature of these groups is mostly uncertain.

Table 2 provides the list of cluster candidates identified from the 6-cluster model. Cluster candidates that are included in both the best AIC model and best BIC model are listed first, followed by cluster candidates that only appear in the best AIC model. These clusters are listed in approximate order of significance. Uncertainties on model parameters related to cluster shape are large. For example, in all cases the radius of the cluster core is poorly constrained. Thus, we do not report cluster core radius, ellipticity, or orientation. The number of stars reported in the table is the total number of stars in the cluster. Figure 10 (right) shows the residual map for a 6-cluster model that is favored by the AIC. Red ellipses show the locations of the clusters.

No RA Dec Unc.
(J2000) (J2000) (arcmin) (stars)
Favored by the AIC and BIC
1 16 59 18 42 34 50 [1.7 1.9] 53
2 16 59 58 42 12 00 [0.9 0.6] 21
3 16 54 28 41 02 50 [1.4 2.0] 30
Consistent with the AIC
4 16 54 33 41 53 20 [2.1 1.6] 33
5 16 53 53 41 18 50 [1.2 0.8] 22
6 16 53 26 42 27 00 [1.2 2.0] 11

Note. – Clusters of VVV -band variables identified using the mixture model analysis. The cluster supported by both AIC and BIC is listed first, followed by clusters supported by only the AIC. Column 1: Cluster number. Columns 2–4: Celestial coordinates of cluster center, and uncertainty on these positions in arcminutes. Column 5: Number of variable stars in the cluster (integrated over the entire field of view). Model properties such as core-radius, ellipticity, and orientation are poorly constrained, so we do not report these value. Note that cluster #4 is NGC 6231.

Table 2: Clusters of VVV Variables
Figure 10: Spatial distribution of near-infrared variables in the VVV fields plotted on a residual map for two models: a smoothly varying model of Galactic field stars (left) and the multi-cluster mixture model (right). When the field stars are modeled, as shown in the left panel, possible clusters of stars stand out as positive (red) residuals. In the right panel, the locations of the clusters favored by the AIC are shown by red ellipses. The residuals at these locations are very small, implying a good fit to the data. The grid lines show right ascension and declination.

8.2 Relation of Clusters of Near-infrared Variables to NGC 6231

VVV tiles “d148” and “d110” cover large angular areas in the Galactic plane, so it is likely that multiple, unrelated clusters of young stars will be identified within the field of view. The densest cluster #1 of VVV variables at coordinates 165918423450 is spatially coincident with the cluster DBSB 176 associated with IRAS 16558-4228 (Dutra et al., 2003; Wang & Looney, 2007). Mid-infrared images of DBSB 176 from the GLIMPSE survey (Benjamin et al., 2003) show significant nebulosity in the region, suggesting active star formation. This cluster contains more VVV variables than NGC 6231, which may be explained if it is younger than NGC 6231. To the north east of DBSB 176 is another strong over-density of VVV variables at coordinates 165958421200.

A cluster of VVV variables is associated with the center of NGC 6231, but this cluster is found in the best AIC model, but not the best BIC model. The relatively low number of VVV variables in NGC 6231 may be related to the relatively low disk fraction in NGC 6231 because Class III pre-main-sequence stars typically have lower near-infrared variability amplitudes than Class 0/I/II young stellar objects. Although the analysis by Saurin et al. (2015) suggested a radius for NGC 6231 of 68 arcmin, in the spatial distribution of VVV variables, there is no radially symmetric over-density of this size. However, the VVV variables include only a small fraction of the cluster members, and may not be sufficient for probing low-density distributions of stars that define cluster’s outer boundary.

Rather than a radially symmetric distribution of VVV variables around NGC 6231, there are clumps with higher densities of variables to the north of NGC 6231, modeled by candidates at 165428410250 and 165353411850. These two groups of stars lie between NGC 6231 and Tr 24. If they do lie at the same distance of the rest of this complex, they could indicate further subclustered structure in the outer northern portions of NGC 6231. However, no information is currently available about their distance.

The least statistically significant cluster candidate in the 6-cluster model is located at 165326422700. This is near enough to NGC 6231 that it is plausible that it is related to the Sco OB1 association, but it may also be an unrelated cluster or association in the Galactic plane.

8.3 Distributions of Candidate Early- and Intermediate-type Stars in VPHAS+

The spatial distribution of the candidate O, B, A, and F-type stars from the VPHAS+ survey is shown in Figure 11. Surface density varies by 1.4 dex. We exclude a region around the center of NGC 6231 from the diagram, because the wings of bright O-stars in NGC 6231 inhibits the detection of A and F stars there.

These stars are also not distributed around NGC 6231 evenly, but concentrated to the north of the cluster instead. The two main peaks in surface density correspond to NGC 6231 and another cluster at coordinates 165510395700, just to the north east of VVV variable cluster #3 from Table 2.

The core region of NGC 6231 derived in §5 is shown as a red ellipse on the map in Figure 11, and an additional ellipse 4 times the size of the core is also shown. It is difficult to define an outer boundary to the cluster because the over-density associated with NGC 6231 blends into the over-densities associated with the larger-scale Sco OB1 region. Without measurements of stellar kinematics or three-dimensional coordinates of stars in the Sco OB1 region, it is impossible to determine whether there is a meaningful astrophysical distinction between the NGC 6231 stellar population and the Sco OB1 population.

The relation between NGC 6231 members and Sco OB1 members may resemble the relation between members of the Orion Nebula Cluster and members of the Orion Molecular cloud. In both cases the distribution of stars resembles a smooth, centrally concentrated cluster on a smaller scale, and a more elongated and clumpy distribution at on a larger scale (Hillenbrand & Hartmann, 1998; Megeath et al., 2012, 2016).

Figure 11: Spatial distribution of candidate O, B, A, and F-type stars in the VPHAS+ fields around NGC 6231. Yellow points mark the candidates selected from the vs.  diagram. The adaptively smoothed surface-density map is shown by the color map. The Chandra field of view is shown as a black polygon, and the 1-core radius and the 4-core radii ellipses for the main NGC 6231 cluster are shown in red. Selection of stars is impeded in the center of NGC 6231 due to bright stars, so this area is not included in the analysis.

9 Discussion:
Cluster Formation and Fate

9.1 Summary of Observational Results

The main cluster properties, as listed below, can be used to place the cluster in a Galactic context and to test theoretical models for cluster formation and evolution.

  1. The median age of the cluster is estimated to be 3.2 Myr, but there may be systematic uncertainty in age. There is evidence of a large age spread with stellar age estimates ranging from 2 to 7 Myr.

  2. The density of stars at the center of the cluster is  stars pc (or a column density of  stars pc.

  3. The total number of stars in the cluster has a lower limit of 5700250 stars, down to the hydrogen burning mass limit. However, the ability to estimate total cluster population is limited by the Chandra field of view. This population corresponds to a total stellar mass of 3300–4200 .

  4. Gas mass does not contribute significantly to the cluster’s gravitational potential.

  5. The radial density distribution of stars resembles an isothermal-ellipsoid distribution with significant () elongation of the cluster. The measured isothermal-ellipsoid core radius is 1.20.1 pc.

  6. There is a second mode in the surface density map, which corresponds to a minor subcluster of stars with 4% the population of the main cluster.

  7. No radial stellar-age stratification is evident.

  8. The spatial distribution of O and B stars shows statistically significant mass segregation, but lower-mass stars shown no sign of mass segregation. The dependence of spatial dispersion on stellar mass is described by a power-law relation .

  9. The cluster follows the “isothermal ellipsoid” surface-density model out to at least 4 times the core radius. However, the distribution of pre-main-sequence stars in a several-square-degree field of view around the cluster is clumpy, rather than radially symmetric.

9.2 Origin of NGC 6231

The initial properties of the molecular cloud and the early dynamical interactions of groups of stars may determine what type of star cluster or association is produced and whether it will survive as an open cluster (Kruijssen, 2012). The lack of remaining cloud material and the dynamically evolved state of NGC 6231 mean that formation characteristics such as the star-formation efficiency cannot be directly measured. However, some properties of the progenitor cloud can be inferred from the observed star cluster.

9.2.1 Filamentary Natal Cloud

Young star clusters are often associated with filamentary molecular clouds on many size scales (e.g., Jackson et al., 2010; Hacar et al., 2013). These elongated clouds may imprint their structure on the star clusters that form within them. For example, in the MYStIX study, several examples of “linear chains of subclusters” are noted (Kuhn et al., 2014, their Figure 5). In some of these, like DR 21, the young stars are deeply embedded in a massive infrared-dark filament. In NGC 1893, multiple subclusters are arranged linearly within a bubble, evacuated of molecular gas.

The distribution of VVV variables around NGC 6231 is not isotropic, but instead concentrated to the north of the cluster in two subclusters. This suggests that there is a population of young stars bridging the gap between NGC 6231 and Tr 24. The concentration of high and intermediate mass stars north of NGC 6231 also shows excess stars to the north of NGC 6231. This spatial distribution suggests that the progenitor cloud for NGC 6231 likely was filamentary with a north-south orientation.

An elongated cloud may also impart ellipticity on the clusters that form. For example, the Orion Nebula Cluster is elongated in a north-south direction (Hillenbrand & Hartmann, 1998) approximately matching the orientation of the molecular filament. Simulations of the collapse of an elongated molecular cloud, starting with the likely initial state of Orion A cloud (Hartmann & Burkert, 2007), do produce an elongated young star cluster (Kuznetsova et al., 2015).

The core of NGC 6231 is also clearly elongated with statistically significant ellipticity (). This is similar to the ellipticity of the Orion Nebula Cluster (–0.5) measured by Kuhn et al. (2014) using the same method. The orientation of the core is close to being north-south, with a moderate offset of 20. Thus, the elongation of the cluster could be explained by cluster formation in a collapsing, elongated cloud, similar to the scenario described by Hartmann & Burkert (2007).

9.2.2 Multiplicity of Star-Forming Cloud Clumps

Molecular clouds are typically clumpy, with fractal like density structures (Stutzki et al., 1998). These can give rise to subclusters of stars (Elmegreen, 2000), which often lie at the locations of dense molecular clumps (Ybarra et al., 2013; Feigelson et al., 2009). The clumpiness of the distribution of stars can vary from region to region, with subclusters of stars formed in individual cloud clumps possibly merging to form more centrally concentrated young star clusters (Fellhauer et al., 2009). Kuhn et al. (2014) find between 1 and 20 subclusters in each of the star-forming regions surveyed by MYStIX.

Our study of NGC 6231 reveals a main cluster accounting for the vast majority of the stars seen in the Chandra field of view. However, a group of stars that is more densely clustered is identified as a statistically-significant subcluster offset from the cluster center. This group of stars most likely formed as a separate density enhancement in the molecular cloud, revealing that the initial cloud was clumpy. The distribution of VVV variables outside the Chandra field of view is also not uniform, suggesting that the progenitor cloud was clumpy on a larger spatial scale as well.

9.2.3 Cluster Assembly

Overall, the structure of NGC 6231 looks quite different from regions with ongoing star formation in the MYStIX study. The regions investigated by MYStIX show significant diversity in the spatial distributions of their stars – some MYStIX rich clusters are smooth with “simple” structures, while others are “clumpy” and/or linear “chains” of subclusters (Kuhn et al., 2014). On one hand, NGC 6231 is very different from clusters like Orion Nebula Cluster and RCW 38 that are much smaller and more concentrated. On the other hand, it is also quite different from clusters like NGC 1893, NGC 6334, and Carina that have complicated clumpy or filamentary morphologies. NGC 6231 is most similar to NGC 2244 in Rosette – both have similar ages, sizes, and numbers of stars, and both are located within evacuated bubbles.

Here we compare the physical properties of NGC 6231 (age, radius, number of stars, density) to the properties of the MYStIX clusters and subclusters. The diversity of these complexes may make it seem as if the (sub)clusters of stars they contain would not be directly comparable. Nevertheless, Kuhn et al. (2015a) has shown that properties of (sub)clusters of stars in MYStIX, even in regions with different global morphology, do follow common trends. A comparison of NGC 6231 to subclusters in MYStIX may make sense if the individual subclusters in a clumpy distribution of stars are the building blocks for more massive clusters, and thus could yield insight into how NGC 6231 formed. In the following discussion, we highlight several examples from MYStIX of more fully-formed clusters that may be better analogs to NGC 6231. These include: RCW 38 (Subcl. B), NGC 2024, W40, Pismis 24 (Subcl. A in NGC 6357), G353.2+0.7 (Subcl. F in NGC 6357), NGC 2362, NGC 6611 (Subcl. B in the Eagle Nebula), and NGC 2244 (Subcl. E in the Rosette Nebula).

Figure 12 shows the cluster properties of the main NGC 6231 cluster (marked by a red square) compared to properties of MYStIX subclusters (black circles) and the highlighted MYStIX clusters (blue squares). The properties include cluster age, cluster core radius, number of stars (within 4 core radii), the projected central stellar density, and the central volumetric stellar density. The relations between these properties were investigated by Kuhn et al. (2015a), who show that this set of subcluster properties is statistically correlated, exhibiting a positive age–radius relation, a positive radius–number of stars relation, and a negative density–radius relation. The uncertainties on the properties of NGC 6231, which result from model fitting, from estimation of completeness, and from age estimation, are shown by the error bars. The regression lines for the relations found in MYStIX are also drawn.

NGC 6231’s properties lie near the regression lines, at the older, larger, more-massive, and less-dense end of the distribution (Figure 12). NGC 2244 also lies at the same end of the parameter distributions. In the scatter plots, the point corresponding to NGC 2244 (Rosette Subcluster E) is the blue point with the largest value of . This cluster is very close in age and radius to NGC 6231; it is slightly less massive and less dense, but is still one of the most massive and least dense clusters in the MYStIX study. In contrast, most of the other rich clusters are also located near the regression lines, but with range of properties. The only massive cluster to significantly deviate from these trends is RCW 38, the blue point with the smallest value of . RCW 38 has 50-times more stars than expected for a cluster its size333A recent analysis of RCW 38 by Mužić et al. (2017) has reported a maximum surface density 6 times lower than reported in MYStIX using near-infrared observations. This difference may arise due to correction for completeness or differences in data-smoothing method to calculate density. The Mužić et al. (2017) central density moves this point closer to the regression line found for the other MYStIX clusters; however, it would still be nearly an order of magnitude more dense than predicted by the regression line. given the regression line, and its unusually high central density has already been noted by Kuhn et al. (2015b).

These observations would suggest that NGC 6231 arose from the same cluster assembly processes that formed the majority of MYStIX subclusters because a different cluster formation process would not necessarily produce a cluster following the same age–radius–mass–density relations. Below we briefly describe the subcluster relations obtained by Kuhn et al. (2015a) and how NGC 6231 relates to each case.

Radius–age relation

The regression line444Several methods exist to obtain linear regression fits to bivariate data. Here we present the slope of the reduced major-axis regression line to the data on a log-log plot. In contrast, Kuhn et al. (2015b) provide the orthogonal regressions, but these are mislabeled as reduced major-axis regressions in their Tables 3 and 4. shown for MYStIX subclusters is . The statistically significant correlation between age and radius was interpreted as cluster expansion, which may be an effect of mass loss (e.g., loss of cloud material), binary stars, subcluster mergers, or a cluster that is initially supervirial or unbound. For an age of 3.2 Myr, NGC 6231 is slightly larger than given by the regression, but well within the scatter observed for MYStIX subclusters. A 6.4 Myr age would place NGC 6231 below the regression line.

Number of stars–radius relation

The regression line shown is . Several effects could produce this relation, including build up of cluster mass while expansion is occurring as suggested by Kuhn et al. (2015a) or a birth relation inherited from the mass–size relation of molecular clumps as suggested by Pfalzner et al. (2016). The coincidence of NGC 6231 with the regression line is quite close. However, given that NGC 6231 has almost certainly expanded from its original size, the second explanation is unlikely in this case.

Density–radius relation

The regression lines shown are and . The slopes of these lines are slightly less steep than the relation that would be expected for a cluster expanding while neither gaining nor losing stars (i.e. ). Kuhn et al. (2015a) proposed that a cluster gaining stars through hierarchical mergers of subclusters could produce such a relation. Alternatively, the distribution could be produced by combined effects of initial conditions and evolution of the cluster. Although NGC 6231 is less dense than most MYStIX subclusters, it has a higher density than expected given its radius.

It is intriguing that a simple, isolated, massive cluster like NGC 6231 would follow the same relations as less massive subclusters in complex regions with ongoing star formation. If this is not a coincidence, then it may suggest that the same processes that govern the properties of subclusters in star-forming complexes also govern the properties of fully-formed young clusters. Many of the MYStIX star-forming regions are likely sites of hierarchical cluster assembly. Fellhauer et al. (2009) show that subclusters still embedded in clouds (like MYStIX regions DR 21 or NGC 6334) rapidly merge. Kuhn et al. (2014) note that the variation in morphology of young star clusters (ranging from linear chains of subclusters, clumpy distributions of stars, centrally concentrated clusters) may be an evolutionary progression. Thus, we suggest that the common factor between NGC 6231 and MYStIX subclusters is hierarchical assembly.

Figure 12: Cluster properties (radius, number of stars, central surface density, central volume density, and age) for NGC 6231 shown in comparison to the properties of subclusters of stars in other star-forming region from the MYStIX study. The red points mark NGC 6231 on these diagrams, and the uncertainties are shown by the red error bars. The MYStIX (sub)clusters are shown as black or blue points, typical uncertainties are indicated by black crosses, and the orthogonal regression lines are shown as a dashed, black line. Gray lines indicate simplistic evolutionary tracks for clusters, whereby a cluster expands at a constant radial velocity and neither gains nor looses stars. For the MYStIX subclusters, only points with no missing values for age, , , and are included. Eight rich clusters in MYStIX, which make the best comparisons to NGC 6231, are highlighted in blue: RCW 38, NGC 2024, W40, Pismis 24, G353.2+0.7, NGC 2362, NGC 6611, and NGC 2244 (from smallest to largest in radius).

9.3 Current State

The CXOVVV catalog can reveal the effects of the cluster’s dynamical state on its spatial structure, but kinematic data are necessary to obtain fundamental properties of the current state, such as the cluster’s total energy and how the energy is partitioned. In the near future, measurements of proper motion are expected to become available from the Gaia survey. For the stars in the CXOVVV catalog, it is estimated that the precision will be 0.1–2.0 km s, although the performance may be degraded near the Galactic plane (de Bruijne et al., 2005). These measurements could be used to test some of the inferences about the cluster’s current state from this paper.

9.3.1 Radial Structure

The projected density of stars in NGC 6231, based on star counts in the CXOVVV catalog, is one of the most accurately known of nearby young star clusters. This has allowed us to test a variety of empirical cluster models. We find that the isothermal ellipsoid model is a remarkably good fit, while other commonly used models are not. We note that the isothermal ellipsoid is a special case of the EFF profile (Elson et al., 1987), which has been found to describe the distribution of light in very massive young clusters in the Milky Way and nearby galaxies. Recently, Grudić et al. (2017) have argued that EFF surface-density profile arises from hierarchical cluster assembly.

Another theoretical implication of this result is that the formation of star clusters is not an entirely scale-invariant process, given that the resulting cluster does have a characteristic length scale, . This is interesting because many of the astrophysical processes of star formation are scale invariant, leading to fractal-like distributions of star-formation activity in the Galaxy. Nevertheless, within individual star clusters, some process must produce the observed length scales.

9.3.2 Cluster Expansion

NGC 6231’s current core radius of 1.20.1 pc is a factor of 15 larger than the typical core radius of highly-absorbed subclusters of stars in MYStIX, and a factor of 7 larger than the typical MYStIX young stellar cluster (Kuhn et al., 2014). Other studies also suggest that typical embedded clusters have half-mass radii of a few tenths of parsecs (e.g., Marks & Kroupa, 2012). Thus, NGC 6231’s size is likely dominated by expansion, not its initial formation size.

The cause of cluster expansion can have an effect on a cluster’s radial structure, including the surface-density profile and the presence or absence of a radial–age gradient. Expansion of young clusters can be driven by mass loss from the dispersal of a molecular cloud, while older clusters can expand due to mass loss from stellar evolution. However, a gravitationally bound cluster will expand even in the absence of these effects, and Gieles et al. (2012) has suggested that this expansion can explain the surface-density distribution of pre-main-sequence stars in our Galactic neighborhood. This expansion is driven by transfer of energy from binaries, mass segregation, and mass loss due to dynamical relaxation, and the resulting expansion is homologous (preserving radial structure) and scale-invariant in time (Giersz & Heggie, 1996).

In contrast, an unbound cluster may have a different structure. If the total energy of a cluster is highly positive, the cluster will evolve toward a structure determined by the velocity dispersion of its stars. If stars have a Maxwell-Boltzmann distribution, this would be a multivariate normal surface density distribution. Alternatively, in a cluster where a fraction of stars are escaping, models by Bastian & Goodwin (2006) show that the escaping stars may lead to an excess number of stars at large radii compared to an EFF profile. Given that NGC 6231 is well described by the isothermal ellipsoid model (Figure 3) and does not appear to have excess stars within the Chandra field of view, it is unlikely that the cluster disintegrating in either of these ways.

The test for a radial gradient in stellar ages can also be used to evaluate whether stars have been escaping from the region from the cluster during the star-formation process. If stars in NGC 6231 were free to drift from their initial location as soon as they were formed, one would expect the first stars that formed to have drifted the farthest, creating a radial–age gradient. The age spread in NGC 6231 has been estimated to be –7 Myr (Sana et al., 2007; Sung et al., 2013; Damiani et al., 2016), and most stellar ages calculated in Paper I range from 2.5 Myr to 9 Myr. If stars have drifted outward over this period, the oldest stars would have on average moved farthest from the center of the cluster, producing a clear radial-age gradient. However, no radial age gradient is seen, and the stars of various ages are well mixed. Thus, the cluster must have been gravitationally bound at least until star formation had ceased. In contrast, we note that the MYStIX cluster NGC 2244 does have hints of a radial–age gradient in the analysis of Getman et al. (2014b, their Figure 7c).

9.3.3 Dynamical Evolution

Dynamical timescales are important for cluster evolution. However, the lack of kinematic data for NGC 6231 means that the timescales can only be approximated. Velocity dispersions are likely to be similar to those in other young clusters, which range from 1–3 km s (e.g., Fűrész et al., 2008; Cottaar et al., 2012; Rigliaco et al., 2016). We also use  pc as a characteristic radius for the cluster, yielding a cluster crossing time of –4.6 Myr. The number of stars within this three-dimensional volume555For these calculations we will ignore the effects of cluster elongation, treating the cluster as a radially symmetric distribution of stars. is  stars (§5.4). The timescale for cluster virialization through two-body interactions is approximated by Binney & Tremaine (2008) as


Thus, the two-body virilization timescale for NGC 6231 would be 60 cluster crossing times, or 100–300 Myr, much longer than the cluster has existed. Nevertheless, in the cluster’s evolution so far, it is likely that more rapid dynamical processes (e.g., violent relaxation) have been dominant as we discuss below.

Some structural properties of NGC 6231 suggest that the cluster has undergone some dynamical evolution. This includes a surface density distribution that has a radial profile similar to what is expected from a kinematically isothermal cluster, the segregation of high-mass stars, and the thorough mixing of stars of different ages. On the other hand, some of the observed features of NGC 6231 would likely have been erased in a cluster that had already reached a quasi-equilibrium state: these include the existence of a small subcluster and the possible asymmetry in the main cluster mentioned in Section 5.1. Subclustered structure would be erased during dynamical relaxation, and Hillenbrand & Hartmann (1998) and Feigelson et al. (2005) attribute an asymmetry in the spatial configuration of the Orion Nebula Cluster to ongoing violent relaxation. However, it is also possible that the subcluster is physically separated from NGC 6231 and is just projected onto the main cluster by chance. The age of the NGC 6231 of 2–7 Myr is much less than the 100–300-Myr dynamical-relaxation timescale.

The spatial distribution of stars in NGC 6231 appears to be smoother on shorter length-scales (within the Chandra field of view) and more clumpy on larger length-scales (outside the Chandra field of view). This may be related to different cluster crossing times and dynamical timescales for different length scales. For example, the cluster crossing time of 1.5–4.6 Myr within a radius of is less than or approximately equal to the median age of stars in the cluster. While, the cluster-crossing timescale outside this region, where the distribution is clumpy can be significantly larger than the median age of the stellar population. However, there is not sufficient information to definitively distinguish subclusters of stars in the Sco OB1 complex from unrelated young star clusters or associations along the same field of view

9.3.4 Mass Segregation

Mass segregation has been observed in a number of young star clusters. For example, mass segregation is seen for high-mass stars in some star-forming regions (e.g., the Orion Nebula Cluster; Hillenbrand & Hartmann, 1998) and low-mass stars in others (e.g., W40; Kuhn et al., 2010). A variety of astrophysical phenomena can produce mass segregation, making it difficult to test any particular model. For example, mass segregation will appear in dynamically relaxed clusters due to two-body interactions in which high-mass stars lose energy and sink toward the centers of clusters. The theoretical models for this process are complicated. Several families of models have been discussed by Hénon (1961), Gunn & Griffin (1979), and Gieles & Zocchi (2015). However, mass segregation can also be induced rapidly (in a single cluster crossing time) in young star clusters through violent relaxation or mergers of subclusters (McMillan et al., 2007; Allison et al., 2009; Moeckel & Bate, 2010). It may also be the case that OB stars can only form in certain regions, leading to primordial mass segregation. Parker et al. (2014) found that dynamical mass segregation will increase with dynamical age as clusters become more radially symmetric.

Although OB stars in NGC 6231 are mass segregated, the mass segregation in this cluster is not as strong as in other star-forming regions. The differences in radial distribution between massive stars and low-mass stars (Figure 6) are not as prominent as in the Orion Nebula Cluster (Figure 6 from Hillenbrand & Hartmann, 1998), and the segregation does not extend down to low-mass stars as it does in W40. On the other hand, some other young star clusters do not show mass segregation (e.g., NGC 2244; Wang et al., 2008). Both NGC 6231 and NGC 2244 are older than Orion and W40, suggesting that mass segregation in newly formed clusters may not be just a simple effect of dynamical age.

Mass segregation can arise from the dissintegration of an unbound cluster if stellar velocities depend on mass because, after a stellar association had expanded sufficiently, the distance of a star from the center of the cluster would be proportional to its velocity. If such a cluster had energy equipartition (), this would result in a power-law relation with index between stellar mass and characteristic radii. Figure 7 shows that the regression line for core radius versus average mass is slightly less steep than a relation.

9.3.5 Comparison to Very Massive Young Clusters

Insight into cluster formation mechanisms can also be gained by comparing NGC 6231 to very massive young star clusters. The review by Portegies Zwart et al. (2010) gives a sample that includes several of the most extreme young star clusters ( ) in the Milky Way and nearby galaxies, hereafter the PZ sample. In contrast, NGC 6231 and the MYStIX star-forming regions are all in our relatively quiet neighborhood of the Galactic Disk; neither the starburst at the Galactic Center nor where the Galactic Bar ends. Figure 13 shows NGC 6231, the MYStIX subclusters, and the PZ clusters on the same plot of radius vs. number of stars.666There is one cluster in common in the two samples: Tr 14. The PZ properties of this cluster come from Ascenso et al. (2007). The measured central densities of stars in the clusters are similar,  stars/pc (PZ) vs. 10 stars/pc (MYStIX), and the measured core radii are similar, 0.14 pc (PZ) vs. 0.2 pc (MYStIX). Nevertheless, the total number of stars reported by Ascenso et al. (2007) is 5 times greater than reported by MYStIX. The explanation for this is that Ascenso et al. (2007) include a much larger fraction of the young stellar population of the Carina Nebula Cluster as part of Tr 14 than is included in MYStIX. Thus, we only include the Tr 14 measurement from MYStIX; however, the reader is cautioned that similar differences in definition of cluster boundary may affect other starburst clusters from the PZ sample. NGC 6231 lies between these samples in mass, more massive than the MYStIX subclusters and less massive than the PZ clusters. The properties compiled for the PZ clusters include a core radius where surface density decreases by a factor of 2 from the central density, an effective radius that contains half the light from the cluster, and a photometric cluster mass. The core radius is directly comparable to our values of for NGC 6231 and the MYStIX subclusters. We do not have a direct measurement of total mass for NGC 6231 or the MYStIX subclusters. However, for the PZ sample, we can approximate by assuming an average stellar mass of   and assuming that approximately half the cluster mass is contained within a radius . (Given that Figure 13 is a log-log plot, errors resulting from these approximations likely have little effect on the overall distributions of points.)

NGC 6231 has fewer stars than every starburst clusters, as tabulated by Portegies Zwart et al. (2010). However, NGC 6231 also has a larger core radius than most of them, including all of the PZ clusters with estimated ages less younger than NGC 6231, and 16 out of 19 PZ clusters with ages less than 10 Myr. This may suggest that NGC 6231 has expanded more than the PZ clusters. NGC 6231 has a lower central density than all but a few of the oldest clusters in the PZ sample.

From Figure 13, it can be seen that the MYStIX mass–radius relation, which appears to hold for NGC 6231, does not hold for massive young star clusters in general. For the PZ clusters, the range of core-radius values is similar to that of the larger MYStIX clusters and subclusters, but their values are greater by 2 to 4 orders of magnitude, which places them all above the MYStIX mass–radius relation. A subset of 12 PZ clusters do lie on the plot to the upper right of NGC 6231, and these are only slightly off the MYStIX mass–radius relation, but these are a minority of the sample.

The different locations of very massive clusters and less massive clusters/associations represented by MYStIX and NGC 6231 on the plot suggests different mechanisms of cluster assembly and/or different types of cluster evolution. For example, Banerjee & Kroupa (2014, 2015) argue for monolithic (or prompt) formation of the massive young cluster R136 (left-most orange triangle). Pfalzner (2009, 2011) suggest distinct cluster evolution for “starburst” versus “leaky” young star clusters. NGC 6231 and the MYStIX subclusters have masses of the more-common leaky clusters, while the PZ clusters have masses of the less-common starburst clusters.

Figure 13: Comparison of NGC 6231 to massive young clusters from Portegies Zwart et al. (2010) (both Galactic and extragalactic) and MYStIX subclusters. NGC 6231 is marked by the red square, while the massive young clusters are marked by symbols indicating their age (orange triangles, green diamonds, and blue circles), and the MYStIX subclusters are shown as black circles. The dashed black line shows the regression line found for MYStIX subclusters from Figure 12. Only clusters from the Portegies Zwart et al. (2010) sample with estimates of core radius are included.

9.4 Fate of NGC 6231

The fate of NGC 6231 is tied to the issue of gravitational boundedness. If the system is unbound, the members may form a coherent moving group for a number of Galactic orbits (e.g., Fellhauer & Heggie, 2005), but the surface density will rapidly diminish to the point where the it is nearly indistinguishable from the field (Pfalzner et al., 2015). We have provided evidence that stars were born gravitationally bound to the system, but that the cluster has likely expanded by a factor of 7–15 since its birth. Here, we investigate whether this expansion indicates that the cluster is already unbound or liable to tidal disruption.

9.4.1 Gravitational Boundedness or Unboundness

We can test whether the velocity of stars required to expand NGC 6231 from an initially compact configuration to its current size is greater than the escape velocity. To calculate this velocity, we assume that velocities are mostly radial, which is a reasonable assumption if stars are escaping. For this discussion, we use the mean velocity of a star from its point of origin to its current location, but a star with an outward trajectory would typically be slowing. We also assume that if stars at central distance are unbound, then stars at larger radii will also be unbound. Integrating the mass within a radius for the three-dimensional Hubble model gives the equation for the escape velocity,


where is the gravitational constant. Cluster expansion most likely started 2–7 Myr ago, so the mean velocity of a star at distance from the cluster center would be to . Figure 14 shows a comparison of the escape velocity curve and to the mean velocity a star would need to travel a distance since the cluster expansion began. The gray regions show the effect of uncertainty on cluster mass and uncertainty on when cluster expansion began, which we assume to be approximately equal to the age of the cluster. If  stars pc and expansion started 3.2 Myr ago, the average outward velocity would be lower than the escape velocity out to 8 pc, much larger than the radius. The peak escape velocity would be 2.75 km s.

The expansion of a cluster core to 1 pc during the first several million years is also not inconsistent with simulations of a bound cluster. For example, Kroupa et al. (2001) and Goodwin & Bastian (2006) find expansion of this magnitude in simulations of bound clusters, with the former describing the simulation as a “Pleiades” analog.

Figure 14: The curved black line shows escape velocity as a function of distance from the cluster center , ignoring stars outside radius . The straight black line shows the velocity that a star would need to travel from the center of the cluster to distance in the time since the cluster started to expand. When , stars are not gravitationally bound, assuming they originated near the cluster center. Since this is only the case at large radii, greater than the cluster core radius or the cluster characteristic radius , it is likely that most stars in the cluster are gravitationally bound. We assume a cluster central density of 200 stars pc and a cluster age of 3.2 Myr, but the gray regions show uncertainties in these lines due to the uncertainty in cluster mass and an age range of 2–7 Myr.

9.4.2 Tidal Effects

A young star cluster must be smaller than its Jacobi radius if it is to survive as a open cluster in the Milky Way without disruption due to tidal stripping. The Jacobi radius is


where is cluster mass, is the gravitational constant, and and are Oort’s constants at the location of the cluster (King, 1962, his Equation 24). The values of the Oort constants at the location of NGC 6231 are  km s and  km s using the equations from Piskunov et al. (2006, 2007). For a total cluster mass of 2200–6800 , the Jacobi radius would be –30 pc, or on the sky. Given that 50% of stars reside with a radius  pc from the cluster center, NGC 6231 is currently much smaller than its Jacobi radius.

Young clusters may also be tidally disrupted by molecular clouds and other young star clusters in their natal environment, in a phenomenon known as the “cruel cradle effect” (Kruijssen, 2012). The remaining massive elements of Sco OB1, besides NGC 6231, include the clouds IC 4628 and the Large Elephant Trunk and the clusters Tr 24, VDBH 211, 207, G342.1+00.9, and various other groupings of stars discussed in §8. IC 4628 is likely to be the most massive of these, with a total mass of more than   (Phillips et al., 1986), but more precise mass estimates are not available from the literature. Nevertheless, at an angular distance of from NGC 6231, IC 4628 would need to be at least 20 times more massive than NGC 6231 to tidally disrupt the cluster.

10 Conclusions

The new CXOVVV catalog of 2148 probable cluster members in NGC 6231 from Paper I can be used as a testbed for star-cluster formation theory. This sample is only complete down to 0.5 ; the total intrinsic stellar population of cluster members within the Chandra field of view is estimated to be 5700 stars (§2-3), and the full population is likely to be significantly larger.

The isothermal ellipsoid (Equation 1) provides a remarkably good empirical fit to the surface density distribution of stars in the cluster. This model is notably better than Plummer sphere, multivariate normal, or power-law models. The cluster has core radius  pc and ellipticity . 4% of the stars are in a small subcluster embedded in, or projected upon, the main cluster (§5). Several additional small young clusters are present within 30 pc within the large Sco OB1 complex (§8). Mass segregation is present with a statistical significance of for   and lower significance for lower-mass stars (§6). The empirical dependence of spatial dispersion stellar mass is given by . The median age of stars in the cluster is around 3-4 Myr with a significant age spread, but no radial age gradient is present (§7). The basic structural properties of the cluster measured in this paper are listed in more detail in §9.1 of the Discussion.

NGC 6231 has physical properties similar to other clusters and association in the nearby Galaxy studied in the MYStIX project. NGC 6231 most closely resembles NGC 2244 in the Rosette Nebula. Both clusters are similar in size, number of stars, and age, and both are located in bubbles from which gas has been removed. In addition, NGC 6231 lies on the empirical relationships between core radius, age, mass, and density found for clusters and subclusters in the MYStIX star-forming regions (§9.2.3). This could be considered to be an unexpected result, given that the spatial distributions of stars in MYStIX star-forming regions with clumpy, subclustered structures (like NGC 6334 or the Carina Nebula) are very different from the smooth distribution of stars in NGC 6231. However, this situation may arise if both subclusters in star-forming regions and more evolved clusters like NGC 6231 are assembled through similar processes. In contrast, extremely massive young star clusters in the Milky Way and nearby galaxies from the PZ sample do not follow this relation, possibly suggesting a different mode of cluster formation for very massive clusters (§9.3.5).

Several lines of evidence indicate that NGC 6231 is gravitationally bound. The lack of a radial gradient in stellar age shows that most stars have not been freely drifting away from the cluster during the several million years (approximately 2 to 7 Myr ago) when stars were forming (§9.3.1). The cluster’s expansion velocity is significantly less than the escape velocity for stars, which would not be true for an unbound cluster. The cluster is also significantly smaller than its tidal disruption radius (§9.4.1).

Future measurements of cluster kinematics will be able to better constrain theoretical modeling. Proper motions from Gaia are expected to have precisions of 0.1–2 km s. When combined with object selection using the CXOVVV catalog, the Gaia dataset will provide four-dimensional kinematic information about the cluster, which will allow the conclusions from this paper, based on the two-dimensional spatial distributions of stars, to be tested.

MAK, EDF, MG, NM, JB, and RK acknowledge support from the Ministry of Economy, Development, and Tourism’s Millennium Science Initiative through grant IC120009, awarded to The Millennium Institute of Astrophysics. MAK was also supported by a fellowship (FONDECYT Proyecto No. 3150319) from the Chilean Comisión Nacional de Investigación Científica y Tecnológica, and RK received support from FONDECYT Proyecto No. 1130140. The scientific results reported in this article are based on data obtained from the Chandra Data Archive, the Vista Variables in Vía Lactéa project (ESO program ID 179.B-2002), and the Two Micron All Sky Survey catalog. This work make use of analysis methods developed at Penn State for the MYStIX project and Chandra data reduction procedures (including ACIS Extract) developed by Patrick Broos and Leisa Townsley. We thank the referee for many useful comments and suggestions. We also thank Anushree Sengupta for useful feedback on the article and Amelia Bayo and Estelle Moraux for helpful discussions about star-forming regions. CXO(ACIS-I), ESO:VISTA(VVV) \software ACIS Extract (Broos et al., 2010), astro (Kelvin, 2014), astrolib (Landsman, 1993), astrolibR (Chakraborty & Feigelson, 2015), celestial (Robotham, 2016), SAOImage DS9 (Joye & Mandel, 2003), SIMBAD (Wenger et al., 2000), spatstat (Baddeley et al., 2005b), TOPCAT (Taylor, 2005), XPHOT (Getman et al., 2010)


  • Akaike (1974) Akaike, H. 1974, IEEE transactions on automatic control, 19, 716
  • Allison et al. (2009) Allison, R. J., Goodwin, S. P., Parker, R. J., et al. 2009, ApJ, 700, L99
  • Ankay et al. (2001) Ankay, A., Kaper, L., de Bruijne, J. H. J., et al. 2001, A&A, 370, 170
  • Ascenso et al. (2009) Ascenso, J., Alves, J., & Lago, M. T. V. T. 2009, A&A, 495, 147
  • Ascenso et al. (2007) Ascenso, J., Alves, J., Vicente, S., & Lago, M. T. V. T. 2007, A&A, 476, 199
  • Baddeley et al. (2008) Baddeley, A., Møller, J., & Pakes, A. G. 2008, Annals of the Institute of Statistical Mathematics, 60, 627
  • Baddeley et al. (2015) Baddeley, A., Rubak, E., & Turner, R. 2015, Spatial point patterns: methodology and applications with R (CRC Press)
  • Baddeley et al. (2005a) Baddeley, A., Turner, R., Møller, J., & Hazelton, M. 2005a, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67, 617
  • Baddeley et al. (2005b) Baddeley, A., Turner, R., et al. 2005b, Journal of statistical software, 12, 1
  • Banerjee & Kroupa (2014) Banerjee, S., & Kroupa, P. 2014, ApJ, 787, 158
  • Banerjee & Kroupa (2015) —. 2015, MNRAS, 447, 728
  • Banerjee & Kroupa (2017) —. 2017, A&A, 597, A28
  • Bastian & Goodwin (2006) Bastian, N., & Goodwin, S. P. 2006, MNRAS, 369, L9
  • Benjamin et al. (2003) Benjamin, R. A., Churchwell, E., Babler, B. L., et al. 2003, PASP, 115, 953
  • Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
  • Broos et al. (2010) Broos, P. S., Townsley, L. K., Feigelson, E. D., et al. 2010, ApJ, 714, 1582
  • Broos et al. (2011) —. 2011, ApJS, 194, 2
  • Burham & Anderson (2002) Burham, K., & Anderson, D. 2002, Model Selection and Multivariate Inference: A Practical Information–Theoretical Approach, Springer, New York
  • Cartwright & Whitworth (2004) Cartwright, A., & Whitworth, A. P. 2004, MNRAS, 348, 589
  • Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  • Chakraborty & Feigelson (2015) Chakraborty, A., & Feigelson, E. D. 2015, astrolibR: Astronomy Users Library for R, v0.1, CRAN.
  • Contreras Peña et al. (2017a) Contreras Peña, C., Lucas, P. W., Minniti, D., et al. 2017a, MNRAS, 465, 3011
  • Contreras Peña et al. (2017b) Contreras Peña, C., Lucas, P. W., Kurtev, R., et al. 2017b, MNRAS, 465, 3039
  • Cottaar et al. (2012) Cottaar, M., Meyer, M. R., Andersen, M., & Espinoza, P. 2012, A&A, 539, A5
  • Cutri et al. (2012) Cutri, R. M., et al. 2012, VizieR Online Data Catalog, 2311
  • Dale & Bonnell (2008) Dale, J. E., & Bonnell, I. A. 2008, MNRAS, 391, 2
  • Dale et al. (2013) Dale, J. E., Ercolano, B., & Bonnell, I. A. 2013, MNRAS, 430, 234
  • Damiani et al. (2016) Damiani, F., Micela, G., & Sciortino, S. 2016, A&A, 596, A82
  • de Bruijne et al. (2005) de Bruijne, J., Perryman, M., Lindegren, L., et al. 2005, Gaia–JdB–022
  • de Grijs et al. (2002) de Grijs, R., Gilmore, G. F., Johnson, R. A., & Mackey, A. D. 2002, MNRAS, 331, 245
  • Dekel & Krumholz (2013) Dekel, A., & Krumholz, M. R. 2013, MNRAS, 432, 455
  • Drew et al. (2014) Drew, J. E., Gonzalez-Solares, E., Greimel, R., et al. 2014, MNRAS, 440, 2036
  • Dutra et al. (2003) Dutra, C. M., Bica, E., Soares, J., & Barbuy, B. 2003, A&A, 400, 533
  • Elmegreen (2000) Elmegreen, B. G. 2000, ApJ, 530, 277
  • Elmegreen & Elmegreen (2001) Elmegreen, B. G., & Elmegreen, D. M. 2001, AJ, 121, 1507
  • Elson et al. (1987) Elson, R. A. W., Fall, S. M., & Freeman, K. C. 1987, ApJ, 323, 54
  • Feigelson (2017) Feigelson, E. D. 2017, ArXiv e-prints, arXiv:1704.08115
  • Feigelson et al. (2009) Feigelson, E. D., Martin, A. L., McNeill, C. J., Broos, P. S., & Garmire, G. P. 2009, AJ, 138, 227
  • Feigelson et al. (2005) Feigelson, E. D., Getman, K., Townsley, L., et al. 2005, ApJS, 160, 379
  • Feigelson et al. (2011) Feigelson, E. D., Getman, K. V., Townsley, L. K., et al. 2011, ApJS, 194, 9
  • Feigelson et al. (2013) Feigelson, E. D., Townsley, L. K., Broos, P. S., et al. 2013, ApJS, 209, 26
  • Fellhauer & Heggie (2005) Fellhauer, M., & Heggie, D. C. 2005, A&A, 435, 875
  • Fellhauer et al. (2009) Fellhauer, M., Wilkinson, M. I., & Kroupa, P. 2009, MNRAS, 397, 954
  • Fűrész et al. (2008) Fűrész, G., Hartmann, L. W., Megeath, S. T., Szentgyorgyi, A. H., & Hamden, E. T. 2008, ApJ, 676, 1109
  • Getman et al. (2010) Getman, K. V., Feigelson, E. D., Broos, P. S., Townsley, L. K., & Garmire, G. P. 2010, ApJ, 708, 1760
  • Getman et al. (2014a) Getman, K. V., Feigelson, E. D., & Kuhn, M. A. 2014a, ApJ, 787, 109
  • Getman et al. (2014b) Getman, K. V., Feigelson, E. D., Kuhn, M. A., et al. 2014b, ApJ, 787, 108
  • Gieles et al. (2012) Gieles, M., Moeckel, N., & Clarke, C. J. 2012, MNRAS, 426, L11
  • Gieles & Zocchi (2015) Gieles, M., & Zocchi, A. 2015, MNRAS, 454, 576
  • Giersz & Heggie (1996) Giersz, M., & Heggie, D. C. 1996, MNRAS, 279, 1037
  • Goodwin & Bastian (2006) Goodwin, S. P., & Bastian, N. 2006, MNRAS, 373, 752
  • Gregorio-Hetem et al. (2015) Gregorio-Hetem, J., Hetem, A., Santos-Silva, T., & Fernandes, B. 2015, MNRAS, 448, 2504
  • Grudić et al. (2017) Grudić, M. Y., Guszejnov, D., Hopkins, P. F., et al. 2017, ArXiv e-prints, arXiv:1708.09065
  • Gunn & Griffin (1979) Gunn, J. E., & Griffin, R. F. 1979, AJ, 84, 752
  • Hacar et al. (2013) Hacar, A., Tafalla, M., Kauffmann, J., & Kovács, A. 2013, A&A, 554, A55
  • Hartmann & Burkert (2007) Hartmann, L., & Burkert, A. 2007, ApJ, 654, 988
  • Hénon (1961) Hénon, M. 1961, Annales d’Astrophysique, 24, 369
  • Herbst et al. (1994) Herbst, W., Herbst, D. K., Grossman, E. J., & Weinstein, D. 1994, AJ, 108, 1906
  • Hillenbrand & Hartmann (1998) Hillenbrand, L. A., & Hartmann, L. W. 1998, ApJ, 492, 540
  • Hills (1980) Hills, J. G. 1980, ApJ, 235, 986
  • Jackson et al. (2010) Jackson, J. M., Finn, S. C., Chambers, E. T., Rathborne, J. M., & Simon, R. 2010, ApJ, 719, L185
  • Jeffries et al. (2014) Jeffries, R. D., Jackson, R. J., Cottaar, M., et al. 2014, A&A, 563, A94
  • Joy (1945) Joy, A. H. 1945, ApJ, 102, 168
  • Joye & Mandel (2003) Joye, W. A., & Mandel, E. 2003, in Astronomical Society of the Pacific Conference Series, Vol. 295, Astronomical Data Analysis Software and Systems XII, ed. H. E. Payne, R. I. Jedrzejewski, & R. N. Hook, 489
  • Kass & Raftery (1995) Kass, R. E., & Raftery, A. E. 1995, Journal of the american statistical association, 90, 773
  • Kelvin (2014) Kelvin, L. 2014, astro: Astronomy Functions, Tools and Routines, v1.2, CRAN.
  • King (1962) King, I. 1962, AJ, 67, 471
  • Konishi & Kitagawa (2008) Konishi, S., & Kitagawa, G. 2008, Information criteria and statistical modeling (Springer Science & Business Media)
  • Kroupa et al. (2001) Kroupa, P., Aarseth, S., & Hurley, J. 2001, MNRAS, 321, 699
  • Kruijssen (2012) Kruijssen, J. M. D. 2012, MNRAS, 426, 3008
  • Kuhn et al. (2015a) Kuhn, M. A., Feigelson, E. D., Getman, K. V., et al. 2015a, ApJ, 812, 131
  • Kuhn et al. (2015b) Kuhn, M. A., Getman, K. V., & Feigelson, E. D. 2015b, ApJ, 802, 60
  • Kuhn et al. (2010) Kuhn, M. A., Getman, K. V., Feigelson, E. D., et al. 2010, ApJ, 725, 2485
  • Kuhn et al. (2017) Kuhn, M. A., Medina, N., Getman, K. V., et al. 2017, AJ, 154, 87
  • Kuhn et al. (2014) Kuhn, M. A., Feigelson, E. D., Getman, K. V., et al. 2014, ApJ, 787, 107
  • Kuznetsova et al. (2015) Kuznetsova, A., Hartmann, L., & Ballesteros-Paredes, J. 2015, ApJ, 815, 27
  • Lada & Lada (2003) Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57
  • Lada et al. (1984) Lada, C. J., Margulis, M., & Dearborn, D. 1984, ApJ, 285, 141
  • Lahiri (2001) Lahiri, P. 2001, in Vol. 38 of IMS Lecture 2V)tes Monograph Series, Institute of Mathematical Statistics
  • Landsman (1993) Landsman, W. B. 1993, in Astronomical Society of the Pacific Conference Series, Vol. 52, Astronomical Data Analysis Software and Systems II, ed. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 246
  • Li & Nakamura (2006) Li, Z.-Y., & Nakamura, F. 2006, ApJ, 640, L187
  • Mapelli et al. (2015) Mapelli, M., Vallenari, A., Jeffries, R. D., et al. 2015, A&A, 578, A35
  • Marks & Kroupa (2012) Marks, M., & Kroupa, P. 2012, A&A, 543, A8
  • Maschberger (2013) Maschberger, T. 2013, MNRAS, 429, 1725
  • McMillan et al. (2007) McMillan, S. L. W., Vesperini, E., & Portegies Zwart, S. F. 2007, ApJ, 655, L45
  • Megeath et al. (2012) Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2012, AJ, 144, 192
  • Megeath et al. (2016) —. 2016, AJ, 151, 5
  • Minniti et al. (2010) Minniti, D., Lucas, P. W., Emerson, J. P., et al. 2010, New A, 15, 433
  • Moeckel & Bate (2010) Moeckel, N., & Bate, M. R. 2010, MNRAS, 404, 721
  • Mužić et al. (2017) Mužić, K., Schödel, R., Scholz, A., et al. 2017, MNRAS, 471, 3699
  • Olea (2000) Olea, R. A. 2000, Technometrics, 42, 444
  • Parker et al. (2014) Parker, R. J., Wright, N. J., Goodwin, S. P., & Meyer, M. R. 2014, MNRAS, 438, 620
  • Perry et al. (1991) Perry, C. L., Hill, G., & Christodoulou, D. M. 1991, A&AS, 90, 195
  • Pfalzner (2009) Pfalzner, S. 2009, A&A, 498, L37
  • Pfalzner (2011) —. 2011, A&A, 536, A90
  • Pfalzner & Kaczmarek (2013a) Pfalzner, S., & Kaczmarek, T. 2013a, A&A, 555, A135
  • Pfalzner & Kaczmarek (2013b) —. 2013b, A&A, 559, A38
  • Pfalzner et al. (2016) Pfalzner, S., Kirk, H., Sills, A., et al. 2016, A&A, 586, A68
  • Pfalzner et al. (2014) Pfalzner, S., Parmentier, G., Steinhausen, M., Vincke, K., & Menten, K. 2014, ApJ, 794, 147
  • Pfalzner et al. (2015) Pfalzner, S., Vincke, K., & Xiang, M. 2015, A&A, 576, A28
  • Phillips et al. (1986) Phillips, J. P., de Vries, C. P., & de Graauw, T. 1986, A&AS, 65, 465
  • Piskunov et al. (2006) Piskunov, A. E., Kharchenko, N. V., Röser, S., Schilbach, E., & Scholz, R.-D. 2006, A&A, 445, 545
  • Piskunov et al. (2007) Piskunov, A. E., Schilbach, E., Kharchenko, N. V., Röser, S., & Scholz, R.-D. 2007, A&A, 468, 151
  • Portegies Zwart et al. (2010) Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431
  • Raboud & Mermilliod (1998) Raboud, D., & Mermilliod, J.-C. 1998, A&A, 333, 897
  • Reipurth (2008) Reipurth, B. 2008, in Handbook of Star Forming Regions, Volume II, ed. B. Reipurth (Astronomical Society of the Pacific Monograph Publications), 401
  • Rice et al. (2015) Rice, T. S., Reipurth, B., Wolk, S. J., Vaz, L. P., & Cross, N. J. G. 2015, AJ, 150, 132
  • Rigliaco et al. (2016) Rigliaco, E., Wilking, B., Meyer, M. R., et al. 2016, A&A, 588, A123
  • Robotham (2016) Robotham, A. S. G. 2016, Celestial: Common astronomical conversion routines and functions, Astrophysics Source Code Library, , , ascl:1602.011
  • Sana et al. (2007) Sana, H., Rauw, G., Sung, H., Gosset, E., & Vreux, J.-M. 2007, MNRAS, 377, 945
  • Saurin et al. (2015) Saurin, T. A., Bica, E., & Bonatto, C. 2015, MNRAS, 448, 1687
  • Schwarz et al. (1978) Schwarz, G., et al. 1978, The annals of statistics, 6, 461
  • Siess et al. (2000) Siess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593
  • Stephens (1974) Stephens, M. A. 1974, Journal of the American statistical Association, 69, 730
  • Stolte et al. (2006) Stolte, A., Brandner, W., Brandl, B., & Zinnecker, H. 2006, AJ, 132, 253
  • Stutzki et al. (1998) Stutzki, J., Bensch, F., Heithausen, A., Ossenkopf, V., & Zielinsky, M. 1998, A&A, 336, 697
  • Sung et al. (2013) Sung, H., Sana, H., & Bessell, M. S. 2013, AJ, 145, 37
  • Taylor (2005) Taylor, M. B. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert, 29
  • Telleschi et al. (2007) Telleschi, A., Güdel, M., Briggs, K. R., Audard, M., & Palla, F. 2007, A&A, 468, 425
  • Townsley et al. (2011) Townsley, L. K., Broos, P. S., Chu, Y.-H., et al. 2011, ApJS, 194, 15
  • Tutukov (1978) Tutukov, A. V. 1978, A&A, 70, 57
  • Wang et al. (2008) Wang, J., Townsley, L. K., Feigelson, E. D., et al. 2008, ApJ, 675, 464
  • Wang & Looney (2007) Wang, S., & Looney, L. W. 2007, ApJ, 659, 1360
  • Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9
  • Ybarra et al. (2013) Ybarra, J. E., Lada, E. A., Román-Zúñiga, C. G., et al. 2013, ApJ, 769, 140
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description