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

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 nonmember 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 $2\buildrel{\circ}\over{.} 3\times 1\buildrel{\circ}\over{.} 5$ box surrounding NGC 6231 (VVV tiles "d148" and "d110"). Near-infrared variability in pre-main-sequence stars may be produced by starspots, 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, 2017b), 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, $g\lt 15.5$ mag and $g-i\lt 1.5$ mag, are designed to approximate the selection rules from Damiani et al. (2016), which used the Johnson-Cousins V and I bands instead. From the spatial distribution of these stars (Section 8.3), it is clear that there is also a large unclustered population, which are likely to be field stars.

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 decl.) slice through the best-fit parametric models (gray lines) and the adaptively smoothed data (black lines). The y-axes are plotted with logarithmic values to allow a high contrast in surface density to be shown.

Zoom In Zoom Out Reset image size

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 ($\sigma =0.3$ pc), and the residual ($\mathrm{model}\mbox{--}\mathrm{points}$) 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 R = 0, the number of stars remains finite when $\alpha \gt -2$, 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 the outer region of the cluster, out to $R\sim 4$ 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 northwest 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σ 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 $\alpha =-1$) 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 ($\alpha =0.7$) is even flatter in the outer regions of the cluster than the singular isothermal ellipsoid model, poorly matching the adaptively smoothed data.

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 M17, 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."

Table 1.  Best-fit Cluster Models

Parameter	Units	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]
rcore	(pc)	1.2 ± 0.1	1.2 ± 0.1	0.048 ± 0.031
r4	(pc)	4.6 ± 0.4	4.6 ± 0.4	0.19 ± 0.12
${N}_{\mathrm{core}}$	(stars)	1400 ± 100	1300 ± 100	17 ± 10
N4	(stars)	5600 ± 250	5400 ± 240	72 ± 20
${N}_{\mathrm{FOV}}$	(stars)	5700 ± 250	5500 ± 240	240 ± 40
${{\rm{\Sigma }}}_{0}$	(stars pc−2)	470 ± 80	460 ± 60	3400:
${\rho }_{0}$	(stars pc−3)	200 ± 50	200 ± 50	35,000:
		0.34 ± 0.05	0.33 ± 0.05	0.66 ± 0.46
ϕ	(deg)	159 ± 4	162 ± 5	165 ± 9
${ \mathcal L }$		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 one core radius, within four 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 ${M}_{\odot }$. The reported uncertainties represent statistical uncertainty in both estimation of number of stars and model fitting but exclude possible systematic error discussed in Section 3.

Download table as:  ASCIITypeset image

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 r_c for the isothermal ellipsoid model appears to be such a length scale. The core radius, ${r}_{c}=1.2\pm 0.1$, is inconsistent with a value of 0. In addition, both scale-invariant models—the singular isothermal ellipsoid and 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 northwest of the cluster center, at coordinates 16^h54^m006, –41°48'07''. 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, 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 $M\geqslant 3$ ${M}_{\odot }$. 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.

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 northwest—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 information criterion (AIC; Akaike 1974) and Bayesian information criterion (BIC; Schwarz 1978), defined by the formulae

$$\begin{eqnarray}&&\mathrm{AIC}=-2{ \mathcal L }+2k,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\mathrm{BIC}=-2{ \mathcal L }+k\mathrm{ln}n,\end{eqnarray} \tag{ 8 }$$

where ${ \mathcal L }$ is the log-likelihood, k is the number of parameters in the model, and n is the number of points. For the isothermal ellipsoid clusters, k is equal to 6 times the number of clusters, and we select the k that minimizes the AIC or BIC. The BIC generally favors simpler models than the AIC because of 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 overdensity of stars to the northeast 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, ${\rm{\Delta }}\mathrm{BIC}=36$, is far above the threshold ${\rm{\Delta }}\mathrm{BIC}\simeq 8-10$ for confident preference of the two-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).

Zoom In Zoom Out Reset image size

Table 1 presents the best-fit models for the single-component model and the two-component 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, ${\rho }_{0}$, is related to the the central surface density, ${{\rm{\Sigma }}}_{0}$, by the relation ${\rho }_{0}={{\rm{\Sigma }}}_{0}/2{r}_{c}$. This equation will also apply to the ellipsoidal model if we assume that the three-dimensional ellipsoid has a harmonic-mean core radius of r_c.

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 16^h54^m159, −41°49'59'' 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.⁷

The subcluster is located at 16^h54^m16, −41°48'13''. 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, 240 ± 40 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 Figure 3 (bottom row) and 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.

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 r_eff 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: N₄, the number of stars (corrected for incompleteness) within a region 4 times larger than the cluster core, and the corresponding radius ${r}_{4}=4{r}_{c}$. 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, $\mathrm{log}({r}_{\mathrm{eff}}/{r}_{c})\sim 0.7\pm 0.4$ dex, so r₄ is likely within a factor of several of the half-mass radius. For ${N}_{\mathrm{core}}$ and N₄, we provide the uncertainty on the number of stars within the reported radii, including only the statistical uncertainties on number of stars. For ${{\rm{\Sigma }}}_{0}$ and ${\rho }_{0}$, uncertainty from estimation of core radius is also included.

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 Two Micron All Sky Survey catalog to estimate where the spatial overdensity of stars meets the background level, and they report a radius of 36.2 pc. Our investigation in Section 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 ${N}_{\mathrm{core}}\approx 5700$, the number of stars projected in an ellipse with characteristic radius r₄ is ${N}_{4}\approx 5400$, and the number of stars within a three-dimensional ellipsoid with characteristic radius r₄ is ∼4300. 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 of $\bar{m}=0.61\,{M}_{\odot }$ for the mass range 0.08–150 ${M}_{\odot }$ and an average mass of $\bar{m}=0.78\,{M}_{\odot }$ for systems. When this uncertainty is combined with the statistical uncertainty on total number of stars, the total mass of stars projected within the four-core radius ellipse, down to the hydrogen-burning limit, is in the range of 3300–4200 ${M}_{\odot }$.

The right panel of Figure 5 shows an estimate of the expected value of $\mathrm{log}M$ as a function of position of a star in the field of view. Interpolation 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 (x, y) is set to the distance to the 100th-nearest star. The resulting map shows the highest average $\mathrm{log}M$ at the center of the cluster, as would be expected if the cluster were mass segregated. The mean value of this peak is $\langle \mathrm{log}M/{M}_{\odot }\rangle =0.44$ (∼2.8 ${M}_{\odot }$). 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 $\langle \mathrm{log}M/{M}_{\odot }\rangle =0.44$ has a p-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 ($\mathrm{log}M/{M}_{\odot }\lt 0.25$), intermediate-mass stars ($0.25\lt \mathrm{log}M/{M}_{\odot }\lt 0.9$), and high-mass stars ($\mathrm{log}M/{M}_{\odot }\,\gt 0.9$). 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 p-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 M. A. Kuhn et al. (2017, 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.

Zoom In Zoom Out Reset image size

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 $\mathrm{log}M/{M}_{\odot }=0.0,0.25,0.5,1$, containing 251, 223, 126, 75, and 41 stars, respectively. Figure 7 shows the plot of subcluster core radius versus 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 $1\sigma$ uncertainty on core radius). The gray dashed lines show relations of ${r}_{c}\propto {M}^{-1/2}$ for comparison.

Zoom In Zoom Out Reset image size

Figure 7 shows that, aside from the second mass bin ($0.5\lt M\lt 1.0$ ${M}_{\odot }$), the core radius decreases monotonically with increasing stellar mass. However, for the range 0.5–3 ${M}_{\odot }$, 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 ${M}_{\odot }$, the difference in core radius only becomes statistically significant for the highest-mass bin (stars with $M\gt 10$ ${M}_{\odot }$), 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 ${r}_{c}\propto {M}^{-0.29\pm 0.06}$ (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 ${\sigma }_{v}(m)\propto {m}^{-1/2}$ relation are shown for comparison.

We use a mixture model approach to identify possible clusterings of variable stars measured in the VVV K_s 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 b. 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,

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{unc}.}({\ell },b)=C\exp [{a}_{1}b+{a}_{2}{b}^{2}],\end{eqnarray} \tag{ 9 }$$

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

The left panel of Figure 10 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 ${a}_{1}=-1.47$ deg⁻¹ and a₂ = 0.26 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 southeast of NGC 6231 near the Galactic plane, and several peaks to the north and south of NGC 6231.

Zoom In Zoom Out Reset image size

Next, we test various mixture models, composed of G "isothermal ellipsoid" components and an "unclustered" component. The identification of multiple statistical clusters follows the same method as Kuhn et al. (2014). Models with G = 0, 1, 2, 3, 4, 5, and 6 are fit, and the AIC and BIC are calculated. For this model, the number of parameters is $k=6G+3$, 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 G = 0, 1, 2, 3, 4, 5, and 6 clusters, respectively. Thus, the BIC clearly favors a three-cluster model, while the AIC is consistent with models with three to six 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 six-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. The right panel of Figure 10 shows the residual map for a six-cluster model that is favored by the AIC. Red ellipses show the locations of the clusters.

Table 2.  Clusters of VVV Variables

No	R.A.	Decl.	Unc.	${N}_{\star }$
	(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 K_s-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 values. Note that cluster 4 is NGC 6231.

Download table as:  ASCIITypeset image

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 16^h59^m18^s, −42°34'50'' 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 northeast of DBSB 176 is another strong overdensity of VVV variables at coordinates 16^h59^m58^s, −42°12'00''.

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 in 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 overdensity 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 the 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 16^h54^m28^s, −41°02'50'' and 16^h53^m53^s, −41°18'50''. 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 six-cluster model is located at 16^h53^m26^s, −42°27'00''. 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.

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 inhibit the detection of A and F stars there.

Zoom In Zoom Out Reset image size

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 16^h55^m10^s, −39°57'00'', just to the northeast of VVV variable cluster 3 from Table 2.

The core region of NGC 6231 derived in Section 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 overdensity associated with NGC 6231 blends into the overdensities 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 on a larger scale (Hillenbrand & Hartmann 1998; Megeath et al. 2012, 2016).

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 ${\rho }_{0}=200\pm 50$ stars pc⁻³ (or a column density of ${{\rm{\Sigma }}}_{0}=460\pm 60$ stars pc⁻²).
• 3.  
  The total number of stars in the cluster has a lower limit of 5700 ± 250, 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 ${M}_{\odot }$.
• 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 ($\epsilon =0.33\pm 0.05$) elongation of the cluster. The measured isothermal ellipsoid core radius is 1.2 ± 0.1 pc.
• 6.  
  There is a second mode in the surface density map, which corresponds to a minor subcluster of stars with 4% of 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 ${r}_{c}\propto {M}^{-0.29}$.
• 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.

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

9.2.2. Multiplicity of Star-forming Cloud Clumps

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 the 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) have 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 four 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.

Zoom In Zoom Out Reset image size

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 r₄. This cluster is very close in age and radius to NGC 6231; it is slightly less massive and less dense, but it 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 a range of properties. The only massive cluster to significantly deviate from these trends is RCW 38, the blue point with the smallest value of r₄. RCW 38 has ∼50 times more stars than expected for a cluster its size⁸ given the ${N}_{4}\sim {r}_{4}$ 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 line⁹ shown for MYStIX subclusters is ${r}_{4}\propto {\mathrm{age}}^{1.8}$. 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 ${N}_{4}\propto {r}_{4}^{1.4}$. Several effects could produce this relation, including buildup 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 ${\rho }_{0}\propto {r}_{4}^{-2.4}$ and ${{\rm{\Sigma }}}_{0}\propto {r}_{4}^{-1.6}$. 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., ${\rho }_{0}\propto {r}_{4}^{-3}$). 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, to clumpy distributions of stars, to centrally concentrated clusters) may be an evolutionary progression. Thus, we suggest that the common factor between NGC 6231 and MYStIX subclusters is hierarchical assembly.

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

9.3.2. Cluster Expansion

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 to 3 km s⁻¹ (e.g., Fűrész et al. 2008; Cottaar et al. 2012; Rigliaco et al. 2016). We also use r₄ = 4.6 pc as a characteristic radius for the cluster, yielding a cluster-crossing time of ${t}_{\mathrm{cross}}=1.5\mbox{--}4.6$ Myr. The number of stars within this three-dimensional volume¹⁰ is $N\approx 4300$ (Section 5.4). The timescale for cluster virialization through two-body interactions is approximated by Binney & Tremaine (2008) as

$$\begin{eqnarray}&&{t}_{\mathrm{vir}}={t}_{\mathrm{cross}}\times \displaystyle \frac{N}{8\mathrm{ln}N}.\end{eqnarray} \tag{ 10 }$$

Thus, the two-body virialization 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 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 r₄ 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

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 (M ∼ 10⁴–10⁷ ${M}_{\odot }$) 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 in the starburst at the Galactic center nor at the end of the Galactic bar. Figure 13 shows NGC 6231, the MYStIX subclusters, and the PZ clusters on the same plot of radius versus number of stars.¹¹ 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 r_c where surface density decreases by a factor of ∼2 from the central density, an effective radius r_eff that contains half the light from the cluster, and a photometric cluster mass. The core radius is directly comparable to our values of r_c 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 N₄ by assuming an average stellar mass of $\bar{m}=0.61\,{M}_{\odot }$ and assuming that approximately half the cluster mass is contained within a radius ${r}_{4}=4{r}_{c}$. (Given that Figure 13 is a log–log plot, errors resulting from these approximations likely have little effect on the overall distributions of points.)

Zoom In Zoom Out Reset image size

NGC 6231 has fewer stars than every starburst cluster, 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 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 N₄ values are greater by 2–4 orders of magnitude, which places them all above the MYStIX mass–radius relation. A subset of ∼12 PZ clusters do lie on the r₄–N₄ 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 r₄–N₄ plot suggest 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 (leftmost 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.

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 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 Unboundedness

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 r are unbound, then stars at larger radii will also be unbound. Integrating the mass within a radius r for the three-dimensional Hubble model gives the equation for the escape velocity,

$$\begin{eqnarray}&&{v}_{\mathrm{esc}}=\sqrt{\displaystyle \frac{8\pi {{Gr}}_{c}^{3}{\rho }_{0}\bar{m}}{r}\left({\sinh }^{-1}(r/{r}_{c})-\displaystyle \frac{r/{r}_{c}}{{({r}^{2}/{r}_{c}^{2}+1)}^{1/2}}\right)},\end{eqnarray} \tag{ 11 }$$

where G is the gravitational constant. Cluster expansion most likely started 2–7 Myr ago, so the mean velocity of a star at distance r from the cluster center would be ${v}_{\mathrm{mean}}=r/(2\,\mathrm{Myr})$ to $r/(7\,\mathrm{Myr})$. Figure 14 shows a comparison of the escape velocity curve to the mean velocity a star would need to travel a distance r 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 ${\rho }_{0}\,=200$ 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 r₄ radius. The peak escape velocity would be 2.75 km s⁻¹.

Zoom In Zoom Out Reset image size

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.

9.4.2. Tidal Effects

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

$$\begin{eqnarray}&&{r}_{J}=\sqrt[3]{\displaystyle \frac{{GM}}{4A(A-B)}},\end{eqnarray} \tag{ 12 }$$

where M is cluster mass, G is the gravitational constant, and A and B 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 $A\approx 11.8$ km s⁻¹ and $A-B\approx 22.2$ km s⁻¹ using the equations from Piskunov et al. (2006, 2007). For a total cluster mass of 2200–6800 ${M}_{\odot }$, the Jacobi radius would be r_J = 20–30 pc, or $0\buildrel{\circ}\over{.} 7\mbox{--}1\buildrel{\circ}\over{.} 1$ on the sky. Given that ≥50% of stars reside within a radius r₄ = 4.6 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 Section 8. IC 4628 is likely to be the most massive of these, with a total mass of more than 10³ ${M}_{\odot }$ (Phillips et al. 1986), but more precise mass estimates are not available from the literature. Nevertheless, at an angular distance of $1\buildrel{\circ}\over{.} 9$ from NGC 6231, IC 4628 would need to be at least ∼20 times more massive than NGC 6231 to tidally disrupt the cluster.