THE SPATIAL STRUCTURE OF YOUNG STELLAR CLUSTERS. I. SUBCLUSTERS

Young low-mass stars have strong X-ray emission compared to main sequence stars, arising mostly from magnetic reconnection flares from the surface of pre-main sequence stars (Güdel & Nazé 2009; Feigelson 2010). X-ray luminosities are typically L_X ∼ 10⁻⁴–10⁻³ L_bol (Preibisch et al. 2005), several orders of magnitude higher than X-ray emission from main sequence stars like the Sun with L_X ∼ 10⁻⁶ L_bol. Thus, X-ray observations can be used to discriminate young stars both disk-bearing and disk-free in star-forming regions from older Galactic field stars. Although most point sources detected in the Chandra observation discussed here are young stars, some X-ray sources are non-member contaminants that include active galactic nuclei and foreground/background main sequence stars (Getman et al. 2011). Here, members and non-members are distinguished by Bayesian supervised learning using photometric, spectroscopic, and variability properties of the sources in our X-ray/IR data (Broos et al. 2013).

X-ray surveys are also effective at overcoming challenges to studying stars in the complex environments of MSFRs. Optical surveys are limited by absorption, which may exceed 100 mag in the V band, but hard X-rays are effective at penetrating the interstellar medium. Optical and IR studies of stars are also limited by high background levels from nebular emission of the H ii region, which is largely absent in the X-ray. Often young stellar clusters are centered near the regions with the highest IR nebulosity, so the X-ray observations can be more effective for identifying low-mass stars in areas with strong gas and dust emission near OB stars. In addition, the brightness contrast between T-Tauri and OB stars is not so dramatic in the X-ray as in optical and in IR bands (e.g., Cohen et al. 2008). Thus, low-mass stars can be detected in the X-ray, even when they are projected near OB stars in the dense young stellar cores. Finally, the low number of resolved field stars and the high resolution of the Chandra mirrors reduce X-ray source confusion in the cluster centers to a low level for the MYStIX MSFRs.

For X-ray selected complex members, source detection is limited by apparent photon flux in the 0.5–8.0 keV band (Kuhn et al. 2013a; Townsley et al. 2014). In the full MPCM list, sensitivity will vary over a considerable range in different parts of the field of view. This non-uniformity in sensitivity across a MYStIX field is due to two effects: differing integration times at different locations in a Chandra mosaic, and telescope effects such as off-axis mirror vignetting and degradation of the point-spread function. Finally, comparison between MYStIX regions must judiciously consider the different Chandra exposure times and distances to the regions. In the present study, we compensate for the intra-region variations in sensitivity. Inter-region differences in sensitivity are treated in the forthcoming Paper II (M. A. Kuhn et al., in preparation). Readers are cautioned about comparing stellar populations in (sub)clusters of different MYStIX regions from the analysis provided here.

The increased sensitivity to X-ray sources near the center of a Chandra pointing can produce artificially clustered sources; this has been called the "egg-crate effect", when many fields are combined into a mosaic. This problem can be seen clearly in the X-ray source distribution of the Chandra Carina Complex Project (Townsley et al. 2011, their Figure 4) and is quantified by Broos et al. (2011, their Figure 7). Sensitivity is reduced by a factor of ∼four in the outermost portions of a Chandra field compared to on-axis portions. For our spatial analysis, we flatten these inhomogeneities in sensitivity by selecting 11,798 X-ray selected sources that have apparent 0.5–8.0 keV band photon flux greater than the completeness limit for the observation, defined as the minimum apparent photon flux at which a young star has a nearly full chance of being detected anywhere in the field of view (Feigelson et al. 2005).

The completeness limit is calculated empirically with the assumption that the underlying X-ray photon flux distribution follows a power law (e.g., Broos et al. 2011). For the Chandra pointing in the ACIS-I mosaic that has the shortest exposure time, X-ray selected MPCMs are stratified by off-axis angle at 03, 51, 63, 75, 83, and 123. Their apparent 0.5–8.0 keV band photon flux distributions are examined to find a photon flux to truncate the sample such that the photon-flux distributions in all strata look similar. The imposed flux limits, log F_{X, limit}, are given by Column 3 of Table 1, and the number of X-ray selected MPCMs remaining in the flattened sample is given in Column 5. Typically, ∼50% of the X-ray selected MPCMs (Broos et al. 2013, their Table 4) are included in the flattened sample.⁷

The uniform X-ray photon flux limit allows structures to be characterized independently of the Chandra pointing history. However, other spatial biases are present. When variations in gas absorption are present, a fixed X-ray photon flux limit corresponds to a higher intrinsic luminosity limit in high-absorption regions. Deeply embedded regions may be deficient in X-ray MPCMs, but often this is compensated for by the higher fraction of IR-excess sources. Another bias may arise due to mass segregation. X-ray luminosity is correlated to stellar mass (e.g., XEST; Güdel et al. 2007), so detection is not independent of mass and MPCMs may more closely reflect the spatial distributions of the more massive population. A similar bias arises because X-ray observations are less sensitive to stars with circumstellar disks (e.g., Getman et al. 2009). Finally, the use of the spatial prior in the classification of MPCMs (Broos et al. 2013) will produce a bias against discovery of widely distributed MPCMs.

Dusty disks and/or infalling envelopes emit IR radiation above that produced by the young stellar photosphere. The spectral energy distributions (SEDs) may show an IR excess in the K band, but it is strongest in the Spitzer mid-IR bands. Povich et al. (2013) catalog the MYStIX InfraRed Excess Sources (MIRES) with SEDs consistent with disk-envelope-bearing young stars. Povich et al. provide several filters to reduce likely background sources; non-member contaminants include reddened stellar photospheres, star-forming galaxies, active galactic nuclei, (post-)asymptotic giant branch stars, shock emission, unresolved knots of PAH emission, and background young stars. We use here only MIRES stars designated $\mathrm{SED\_FLG}=0$ as likely young stars with protoplanetary disks. For the Orion Nebula and Carina Nebula fields, published lists of probable IR-excess members are used (Megeath et al. 2012; Povich et al. 2011b).

IR-selected MPCMs are not biased by interstellar absorption and are detected down to very low masses. Therefore, the IR-excess selected MPCMs may be more effective than the X-ray selected MPCMs at revealing highly absorbed subclusters of young stars. The importance of the IR-only members to our multi-wavelength survey is illustrated by Getman et al. (2012), where IR-excess stars dominate the stellar population in a triggered star-formation globule on the periphery of an H ii region.

Spitzer IRAC observations of MSFRs suffer from high non-uniform nebulosity, primarily from PAH emission lines in the 3.6, 5.8, and 8.0 μm bands, and from source confusion near the center of rich clusters (Kuhn et al. 2013b). In contrast, Chandra's sensitivity is not affected by nebulosity as it has several times better on-axis spatial resolution than Spitzer. These mid-IR effects may lead to severely decreased MIRES sensitivity in regions with the highest nebular contamination and stellar density. We do not attempt to correct these effects by removing weak sources, as we do for the X-ray selected MPCMs, because any attempt to apply a uniform flux limit across the field of view would leave only the few brightest stars.

The IR-excess selected MPCMs are also biased against widely distributed stars because of the spatial filtering used to remove probable contaminants. This may enhance slight overdensities in source counts, and also interferes with attempts to characterize the distributed populations (Feigelson et al. 2013).

Spectroscopic catalogs of OB stars are available for many regions (e.g., Skiff 2009 and references therein), and all known OB stars in MYStIX regions are included in our spatial analysis because of the large role they play in clustered star formation. However, these lists may be incomplete due to limited spatial coverage by the spectrographic surveys. In addition, OB stars may hide in regions with high absorption and avoid detection in optical surveys (Povich et al. 2011b). Most of the OB stars have X-ray counterparts with luminosities comparable to or higher than the most X-ray luminous T-Tauri stars (e.g., Stelzer et al. 2005). The MYStIX survey uncovers a moderate number of candidate new X-ray emitting OB stars (H. Busk et al., in preparation).

The statistical sample of stars used for the analysis of spatial distributions contains 16,608 MPCM out of 31,754 MPCMs in the entire MYStIX project. The resulting samples range from a few hundred to a few thousand young stars for each MYStIX region (Table 1, Column 4). The sample sizes reflect the intrinsic richnesses of the stellar populations, but are affected by ancillary effects such as region distance, obscuration, nebulosity, and observation exposure times. These observational selection effects will be treated in Paper III; the present study is confined to the observed MPCM population.

Nearly all regions have more X-ray selected stars than IR-excess stars. The major clusters are typically dominated by the X-ray selected sources, but peripheral, highly obscured groupings of stars may be dominated by IR-excess sources. Even for W 40 and DR 21, where IR-excess sources outnumber X-ray sources in the whole region, the densest parts of the clusters are dominated by X-ray selected members.

In this study, we adopt a parametric model of young stellar clusters based on equilibrium configurations of self-gravitating stars in the form of isothermal ellipsoids¹⁰ (Chandrasekhar 1942). The isothermal sphere distribution function may be used as the subcluster model. The isothermal sphere surface density profile can be approximated with the Hubble model

$$\begin{equation} \Sigma (r)=\frac{\Sigma _0}{1+(r/r_c)^2}, \end{equation} \tag{ 1 }$$

where Σ(r) is the surface density of stars at a distance r from the cluster center, which is parameterized by the coordinates of the cluster center, the central surface density Σ₀, and the core radius r_c. This approximation is accurate to <7% for r < 4r_c (Binney & Tremaine 2008). The central volume density ρ₀ is related to the central surface density by

$$\begin{equation} \rho _0 = \Sigma _0 / 2 r_c, \end{equation} \tag{ 2 }$$

and the number of stars N projected within a radius r is

$$\begin{equation} N({<} r)=\pi r_c^2 \Sigma _0 \ln \left(1+r^2/r_c^2 \right). \end{equation} \tag{ 3 }$$

However, many young stellar clusters have elliptical shapes on the sky, for example Hillenbrand & Hartmann (1998) find that the ONC has elliptical contours of constant stellar surface density with ellipticity ≃ 0.30. A more general surface density profile is produced by stretching the Hubble model to model ellipsoidal clusters. This fitting function introduces two new parameters, the ellipticity = (a − b)/a and the ellipse orientation ϕ. The surface density for the isothermal ellipsoid becomes

$$\begin{eqnarray} && \Sigma _\mathrm{ell}(\mathbf {r};\Sigma _0,r_0,r_c,\phi ,\epsilon) \nonumber\\ && = \Sigma _0 \left[\! 1\,{+}\, \left| \left(\begin{array}{@{}cc@{}}(1-\epsilon)^{-\frac{1}{2}} & 0 \\ 0 & (1-\epsilon)^{1/2} \\ \end{array} \right) \hat{R}(\phi)(\mathbf {r}-\mathbf {r_0}) \right|^2\! \left /\right. \! r_c^2 \!\right]^{-1}\! , \nonumber\\ \end{eqnarray} \tag{ 4 }$$

where the matrix $\hat{R}(\phi)$ rotates the vector r − r₀ = (x − x₀, y − y₀) by angle ϕ. A single isothermal ellipsoid component has six parameters: central right ascension x₀, central declination y₀, Σ₀, r_c, , and ϕ.

In addition to the subcluster components, a spatially flat unclustered component, Σ_U(r) = Σ_U, will be included to model non-clustered young stars in the region (and any remaining background contaminants). Distributed populations of young stars have been found in NGC 6334 and Carina, where they may represent a previous generation of star formation (Feigelson et al. 2009, 2011).

The finite mixture model for a star-forming region will be the sum of the ellipsoid models for all k subclusters (plus the unclustered component) multiplied by mixture coefficients, a_i, given by

$$\begin{eqnarray} \Sigma _\theta (\mathbf {r}) &=& \sum _{i=1}^{k+1}a_i\Sigma _{i}(\mathbf {r};\theta _i) \nonumber\\ & =& \sum _{i=1}^{k}a_i\Sigma _\mathrm{ell}(\mathbf {r};x_{0,i},y_{0,i},r_{c,i},\phi _i,\epsilon _i)+a_{k+1}\Sigma _{\mathrm{U}}, \end{eqnarray} \tag{ 5 }$$

where θ = {a₁, x_{0, 1}, y_{0, 1}, r_{c, 1}, ϕ₁, ₁, ..., a_k, x_{0, k}, y_{0, k}, r_{c, k}, ϕ_k, _k, a_{k + 1}} denotes the model parameters. Our finite mixture model is only defined over the field of view, so it is integrable and can be used as a probability density function.

We adopt the method of maximum likelihood estimation (MLE) that, following the principles laid out by Fisher (1922), seeks the most likely model given the data. The log-likelihood function of the model parameters θ is given by the equation

$$\begin{equation} \ln L(\theta ;\mathbf {X})=\sum _{i=1}^{N} \ln \Sigma _{\theta }(\mathbf {r}_i)-\int _W \Sigma _{\theta }(\mathbf {r}^\prime){d}\mathbf {r}^\prime , \end{equation} \tag{ 6 }$$

under the assumption of a Poisson point process $\mathbf {X}=\lbrace \mathbf {r}_i\rbrace _{i=1}^N$ containing N points (stars), which are spatially distributed with the surface density of the finite mixture model Σ_θ in a window (field of view) W. The integral in the second term is the expected number of stars within the field of view, and the model Σ_θ is normalized so that this term equals the total number of observed stars N.

MLE of mixture models for cluster analysis is described by Everitt et al. (2011, Chapter 6). We seek the mode (highest value) of the likelihood function. This can be computationally challenging; for example, if five subclusters are present in a MYStIX field, then the parameter space to be searched to maximize the likelihood has 6 × 5 = 30 dimensions. Commonly used computational methods, including the EM Algorithm (McLachlan & Krishnan 2008) and Markov Chain Monte Carlo (MCMC) procedures (Brooks et al. 2011), will find good solutions but not necessarily the global maximum.

Our problem is simpler than many problems where MCMC procedures are needed. As we are seeking physical associations of stars in the two-dimensional (R.A., decl.) subspace of a high-dimensional parameter space, any acceptable subcluster in the final model must be evident as a subcluster in the (R.A., decl.) diagram. We therefore start the MLE computation with a superset of possible clusters formed by visual examination of the MPCM spatial distribution and its adaptively smoothed surface-density map. A Nelder–Mead optimization algorithm, implemented in R's function "optim", is applied to iteratively find the global maximum likelihood solution. Visual examination of a smoothed map of residuals between the data and the model shows that the MLE model solutions are quite satisfactory in explaining the MPCM spatial structure using isothermal ellipses.

The R code for model fitting is presented in the Appendix, using the NGC 6357 MSFR as an example. As is often the case for complex model fits, multiple fitting attempts are needed to find the best-fitting model, with attention paid to individual parameters.

The above steps require that the number of subclusters k in the model be specified. Model selection is the process by which k is varied and the "best" value of k is chosen that balances optimal fit to the data and some concept of model parsimony. This requires that the maximum likelihood for each value of k be "penalized" by some function of k. The most common choices are the BIC and the AIC:

$$\begin{eqnarray} {\rm BIC} &=& {-}2 \ln L + (6k+1) \ln N \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} {\rm AIC} &=& {-}2\ln L + 2(6k+1), \end{eqnarray} \tag{ 8 }$$

where 6k + 1 is the number of parameters for k subclusters, and N is the total number of stars in the MPCM sample. The optimal k value is chosen to minimize the BIC or AIC. When models are found with AIC values very close to the best-fit model, the final model is chosen to minimize maximum excursions in the smoothed residual maps. The BIC places a stronger emphasis of model parsimony (simplicity) than the AIC, whereas the AIC sometimes overestimates the number of components. The choice of AIC, BIC, or other model selection approach is widely debated (Lahiri 2001; Konishi & Kitagawa 2008). Burnham & Anderson (2002) generally favor the AIC, while the BIC has more theoretical support (Kass & Wasserman 1995). Everitt et al. (2011, Section 6.5.2) discuss the issue in the context of mixture models for cluster analysis.

Experimentation with MYStIX samples showed that the BIC model selection criterion produces final models with fewer subclusters, missing sparser clusters that are readily seen in visual inspection of the MPCM spatial distribution. For a typical MYStIX field with N ≃ 1000, the addition of an additional subcluster is accepted to minimize the BIC only if it improves the log likelihood by 6 × ln 1000 ≃ 40 in the likelihood. However, a subcluster that improves the log likelihood by ≃ 12 will be acceptable to the AIC. We therefore have chosen the AIC as our model selection criterion in order to capture a greater dynamic range of rich and poor clusters.

This model selection has two further advantages. First, it obviates the need to arbitrarily specify a critical parameter value in the clustering procedure, such as k in k-nearest neighbors or a special branch length in the minimal spanning tree. Second, it reduces features in the residual map is similar to model selection based on minimizing the "final prediction error," an approach favorably presented by Rao & Wu (2001).

Maps of model residuals are valuable for both evaluating a model fits and identifying where new subclusters could be added to improve a model. The kernel smoothed "raw" residuals for point-process models are defined by (Baddeley et al. 2005, 2008) and implemented by the "diagnose.ppm" function in spatstat. There will be both positive and negative residuals, and the residuals over the entire region will integrate to zero. A better fit will have residuals with lower amplitude that are randomly distributed throughout the region without any coherent structures in the residuals. The residuals for MYStIX clusters will be larger where the surface density is higher; however, there should be roughly equal numbers of positive and negative residuals that are not strongly correlated with the locations of the subclusters.

To smooth the residual maps, we used kernels with bandwidths that are four times the median nearest-neighbor distance in a region, which produce residual maps for the different regions that have similar statistical balances between variance and bias. Kernels that are too small produce lumpy structures due to stochastic fluctuations, while kernels that are too large blur distinct structures together. The residual maps often show coherent positive residuals that cannot be eliminated by adding additional subclusters to the model. These are typically curved or irregular structures that cannot be well modeled by the isothermal ellipsoids used here.

In addition to examination of the smoothed residual map, model goodness of fit may be determined using a χ² test based on the number of counts in arbitrarily chosen spatial bins. If the field of view is divided into n regions A_i with approximately equal number of counts, then the expected number of points in the ith regions is $\mathbb {E}(N_i)=\int _{A_{i}} \Sigma (u){d}u$ for a surface-density model Σ. For a sufficiently large number of counts, we can test goodness of fit using Pearson's statistic

$$\begin{equation} X^2=\sum _{i=1}^{n}\frac{\left[N_i-\mathbb {E}(N_i)\right]^2}{\mathbb {E}(N_i)}. \end{equation} \tag{ 9 }$$

For this test, the field of view is tessellated using adaptive polygonal cells that typically contain ∼20 points to ensure adequate counting statistics. We calculate p values assuming X² is χ² distributed with n degrees of freedom.¹¹

The finite mixture modeling gives probabilities that any given star is part of a particular subcluster. This is a "soft classification", where probabilities are calculated based on the relative contribution of the different components of the posterior model to the total surface density at the location of a star. This stands in contrast to "hard classification" methods that places every object into a single "best" class.

However, approximate lists of stars belonging to each subcluster are very useful for astronomical study, so we generate hard classifications based on the mixture model soft classifier. One possibility for assigning stars to a particular subcluster (or to the unclustered component) is to assign them to the model component that is most probable; this strategy is used by the mclust algorithm (Fraley & Raftery 2012). However, to avoid including the most ambiguous stars in lists for specific subclusters, we additionally require that the probability for the assigned subcluster be better than 30%. We further require that stars assigned to a subcluster be within an ellipse four times the size of the cluster core. (These decision rules play an analogous role to the choice of k in k-nearest neighbor clustering algorithms.)

The choice of algorithm for cluster finding is related to the data that are used and the scientific goals of analysis. The MPCM data have variations in surface density of several orders of magnitude and complex, overlapping structures. Our algorithm resolves these structures into individual components, and constructs a model for a star-forming region based on a sum of simple astrophysically reasonable cluster shapes (isothermal ellipsoids) that can be used as a basis for studying subcluster ages, relationship to molecular material, physical conditions, and evolutionary state in later MYStIX papers. These considerations led us to choose the parametric cluster-analysis methods. However, a this approach may encounter a variety of challenges or disadvantages.

•   
  A model form must be assumed. It is not clear how well finite mixture models work when the component distribution models are poor descriptions of the data. Here we assume that real star clusters have an isothermal ellipsoid distribution function, which is clearly not always true. However, deviations from the isothermal structure can lead to interesting inferences about particular regions. Overall, the clean appearance of the residual maps and excellent global χ² values for the best-fit models (Table 4) show that the assumption of isothermal ellipsoids is satisfactory.
•   
  Guidance on the model selection criteria is lacking. Statisticians have not settled on a universal penalty to MLE for model selection. The AIC has a strong theoretical basis in information theory and thermodynamics, the BIC has a strong theoretical basis in Bayesian inference, and both are related to minimum-χ² fitting under some circumstances. Finite mixture model analysis using the BIC is less sensitive at identifying small clusters in fields of view that also contain rich clusters because the model selection criterion involves improvements to the global likelihood and does not measure localized surface density enhancements. The AIC used here is more effective at finding the small clusters. However, the theory does not exist to interpret changes in the AIC in terms of statistical significance.
•   
  Finding the best mixture model is computationally difficult. Fitting takes place in a high-dimensional parameter space, even though the point pattern data are two-dimensional. Additional problems may be encountered with infinite likelihoods, multimodality of the likelihood (Kalai et al. 2012), and slow convergence of parameter estimators (Chen 1995).

Finite mixture models also have a number of advantages compared to other cluster-finding methods.

•   
  Estimators are based on well-accepted statistical procedures. Finite-mixture model fitting and selection rely on penalized likelihood techniques widely used across statistics. In contrast, nonparametric methods, such as the friends-of-friends and the k-nearest neighbor algorithms, provide no mathematical guidance to estimate the number of clusters. They often exhibit weak performance in the presence of noisy or high-dynamic range structures (Everitt et al. 2011).
•   
  A critical scale or surface density is not assumed. Cluster size and central density are free parameters for each subcluster in the finite mixture models, so structures at vastly different spatial scales can be fit within a single region. This is in sharp contrast to most cluster algorithms, which rely on a single spatial threshold (e.g., critical distance or density) applied throughout the region. Given that it is an open area of research whether or not subclusters of young stars form with a uniform surface density, applying a single spatial threshold for cluster detection risks biasing the result with a biased assumption. Decision rules that play the role of a critical scale are applied only when the stellar membership of a given subcluster needs to be defined.
•   
  Distributed populations are treated. The finite-mixture model does not insist that the clusters include every star, but rather allow relatively isolated stars to be assigned to a distributed population with uniform surface density. In this way, our approach is similar to density based clustering algorithms, rather than standard nonparametric clustering algorithms that require that every object lie in a cluster.
•   
  Probabilistic cluster membership assignments are used. The finite-mixture model is a soft classification method, rather than a hard classification method. Each point has a membership probability associated with each cluster all of which sum to one, and often one probability is much greater than the others. Therefore, different probability thresholds can be used for different applications that require different levels of reliability. The hard classification methods do not include this information, and the ambiguous membership will be definitively assigned into a cluster.
•   
  The models are flexible and effective. The parametric model chosen can be modified to suit the situation being modeled. Other methods may require incorrect assumptions (e.g., normal distributions) or lack any explicit assumptions. An elliptical isothermal model permits a wide array of different morphologies. If the model is astrophysically reasonable, then the parametric clustering method is better optimized for detection structures than a nonparametric method.

Surface density maps of young stars are created using an adaptive density estimator on projected star positions. The density estimator is based on the Voronoi tessellation (Ogata et al. 2003; Ogata 2004; Baddeley 2007). A randomly selected set of stars containing fraction f of the stars is used to create a Voronoi tesselation of the region, which subdivides the field of view into cells that tend to be smaller in dense regions and larger in less dense regions. The rest of the stars are used to estimate the surface density in each cell. We choose a fraction f = 1/10 and repeat this randomized procedure 100 times producing multiple estimates of surface density that are averaged to create the smoothed map.

The adaptively smoothed maps for each star-forming regions are presented in Figure 2, with surface density indicated by color and brown contours showing increases by factors of 1.5. The maps show the surface density of the MPCMs from the "flattened sample" (Section 2). These are observed stellar surface densities, not intrinsic stellar surface densities that depend on the distances and telescope exposure times for each region. The values therefore cannot be directly compared between regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Surface densities vary up to a factor of ∼100 within a single regions, and many regions are multimodal with local maxima differing in surface density by more than an order of magnitude. An effective cluster finding algorithm must be able to identify these different subclusters despite the high contrast in surface density.

Table 2 is a catalog of the 142 subclusters found in the best-fit models of 17 MYStIX MSFRs. For each subcluster (letter designations are assigned from west to east), the table provides celestial coordinates of the subcluster center, semimajor and semiminor axes of the core ellipse in parsecs, ellipticity, orientation, and total number of stars within an ellipse four times the size of the core—all of which are obtained from model fitting (Section 3). In Figure 2, subclusters are shown overlaid on the smoothed surface density maps using a black ellipses to mark the core ellipse of each subcluster.

Table 2. MYStIX Subcluster Properties

Subcluster	α	δ	rc, major	rc, minor		ϕ	N4, obs	J − H	ME
	(J2000)	(J2000)	(pc)	(pc)		(deg)	(stars)	(mag)	(keV)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
Orion A	05 35 14.6	−05 22 31	0.01	0.01	0.23	85	18	0.96	3.6
Orion B	05 35 15.7	−05 23 23	0.06	0.04	0.30	28	73	0.87	1.6
Orion C	05 35 16.7	−05 22 34	0.31	0.16	0.49	5	834	1.05	1.6
Orion D	05 35 17.8	−05 16 35	0.23	0.04	0.84	12	48	1.17	1.4
Flame A	05 41 42.5	−01 54 14	0.15	0.10	0.37	146	219	1.79	2.8
W 40 A	18 31 26.7	−02 05 39	0.17	0.16	0.04	107	187	2.10	2.5
RCW36 A	08 59 27.1	−43 45 22	0.20	0.13	0.36	122	196	1.63	2.3
RCW36 B	08 59 27.3	−43 45 26	0.03	0.01	0.61	25	24	2.14	2.8
NGC 2264 A	06 40 31.5	+09 49 52	0.09	0.08	0.14	136	11	0.83	1.1
NGC 2264 B	06 40 37.1	+09 47 31	0.04	0.02	0.44	124	5	0.60	1.1
NGC 2264 C	06 40 40.4	+09 51 06	0.03	0.03	0.11	8	6	0.73	1.1
NGC 2264 D	06 40 45.8	+09 49 03	0.13	0.11	0.11	194	9	0.65	1.1
NGC 2264 E	06 40 59.1	+09 52 22	0.30	0.16	0.47	83	54	0.63	1.1
NGC 2264 F	06 40 59.2	+09 53 59	0.07	0.03	0.54	69	12	0.61	1.0
NGC 2264 G	06 40 59.4	+09 36 12	0.10	0.07	0.31	96	31	2.11	2.3
NGC 2264 H	06 41 02.1	+09 48 44	0.19	0.16	0.19	19	15	0.65	1.0
NGC 2264 I	06 41 04.5	+09 35 57	0.11	0.05	0.55	14	31	1.16	1.6
NGC 2264 J	06 41 06.3	+09 34 09	0.16	0.12	0.25	63	62	1.22	1.7
NGC 2264 K	06 41 08.2	+09 29 53	0.25	0.11	0.55	0	71	0.75	1.3
NGC 2264 L	06 41 12.8	+09 29 10	0.03	0.03	0.13	83	16	1.84	3.5
NGC 2264 M	06 41 14.9	+09 26 42	0.09	0.06	0.32	20	20	0.73	1.2
Rosette A	06 30 57.1	+04 57 57	1.20	0.87	0.28	90	58	0.79	1.4
Rosette B	06 31 20.5	+04 50 14	0.16	0.16	0.00	0	4	0.73	1.4
Rosette C	06 31 32.0	+04 50 58	0.29	0.09	0.71	129	7	0.77	1.4
Rosette D	06 31 55.4	+04 56 39	0.14	0.04	0.72	87	11	1.17	1.4
Rosette E	06 31 59.3	+04 54 50	0.91	0.83	0.08	111	244	0.76	1.3
Rosette F	06 32 05.5	+04 48 20	0.20	0.09	0.51	0	4	0.83	1.4
Rosette G	06 32 46.5	+04 45 35	0.16	0.16	0.00	0	2	0.76	1.0
Rosette H	06 33 07.2	+04 46 57	0.79	0.15	0.82	93	34	0.84	1.3
Rosette I	06 33 10.2	+04 31 04	0.25	0.11	0.56	42	5	0.94	1.6
Rosette J	06 33 15.1	+04 35 08	0.18	0.09	0.53	156	4	2.18	2.4
Rosette K	06 33 20.1	+04 37 01	0.09	0.09	0.00	90	2	2.00	2.6
Rosette L	06 34 10.7	+04 25 06	1.16	0.57	0.51	167	158	1.10	1.5
Rosette M	06 34 12.9	+04 19 07	0.65	0.25	0.61	71	81	2.21	2.3
Rosette N	06 34 31.4	+04 19 09	0.18	0.16	0.08	119	10	1.62	1.9
Rosette O	06 34 37.1	+04 13 02	0.22	0.05	0.75	151	9	1.76	2.2
Lagoon A	18 03 23.8	−24 15 19	0.35	0.16	0.55	126	30	0.84	1.4
Lagoon B	18 03 40.1	−24 22 40	0.07	0.05	0.28	134	34	1.22	1.8
Lagoon C	18 03 46.3	−24 22 01	0.25	0.12	0.52	109	27	0.85	1.4
Lagoon D	18 03 51.3	−24 21 08	0.10	0.06	0.41	135	9	0.78	1.3
Lagoon E	18 04 07.6	−24 25 53	0.68	0.27	0.60	119	62	0.83	1.3
Lagoon F	18 04 13.3	−24 18 27	1.14	0.73	0.36	3	243	0.80	1.3
Lagoon G	18 04 20.1	−24 22 51	0.09	0.05	0.40	156	17	0.83	1.3
Lagoon H	18 04 23.3	−24 21 13	0.26	0.21	0.20	6	81	0.79	1.3
Lagoon I	18 04 28.3	−24 22 46	0.34	0.31	0.10	17	126	0.83	1.3
Lagoon J	18 04 39.6	−24 23 20	0.25	0.24	0.05	4	47	0.87	1.3
Lagoon K	18 04 50.5	−24 26 19	0.46	0.25	0.45	39	86	1.03	1.4
NGC 2362 A	07 18 37.8	−24 53 58	0.21	0.13	0.38	54	20	0.65	1.2
NGC 2362 B	07 18 42.9	−24 57 44	0.43	0.38	0.10	172	126	0.62	1.1
DR 21 A	20 38 51.7	+42 18 52	0.20	0.06	0.69	60	12	1.97	2.3
DR 21 B	20 38 57.7	+42 17 51	0.07	0.05	0.34	248	5	2.04	2.8
DR 21 C	20 39 00.2	+42 18 47	0.07	0.04	0.41	144	11	3.00	3.5
DR 21 D	20 39 00.4	+42 19 46	0.17	0.07	0.61	168	22	2.92	4.0
DR 21 E	20 39 00.5	+42 22 39	0.22	0.12	0.44	147	75	2.36	3.7
DR 21 F	20 39 00.9	+42 24 43	0.06	0.03	0.55	50	5	2.55	4.0
DR 21 G	20 39 03.8	+42 16 51	0.14	0.11	0.22	149	11	2.51	3.3
DR 21 H	20 39 03.9	+42 25 34	0.11	0.05	0.53	167	22	2.96	3.3
DR 21 I	20 39 05.4	+42 21 18	0.13	0.09	0.32	154	15	2.66	3.0
RCW 38 A	08 58 47.2	−47 30 52	2.57	2.07	0.19	177	314	1.36	2.2
RCW 38 B	08 59 05.0	−47 30 43	0.12	0.08	0.36	54	117	1.06	2.6
RCW 38 C	08 59 16.4	−47 27 50	0.42	0.09	0.78	165	15	1.62	2.5
NGC 6334 A	17 19 57.9	−35 54 02	0.15	0.10	0.31	126	17	1.44	2.1
NGC 6334 B	17 19 58.5	−35 56 24	0.42	0.31	0.25	66	42	1.44	2.0
NGC 6334 C	17 20 02.7	−35 58 23	0.06	0.05	0.15	107	11	1.10	1.8
NGC 6334 D	17 20 14.9	−35 54 43	0.08	0.07	0.09	35	10	1.88	2.8
NGC 6334 E	17 20 19.0	−35 54 59	0.27	0.21	0.24	28	44	1.88	3.1
NGC 6334 F	17 20 23.2	−35 57 00	0.22	0.19	0.17	25	18	1.28	1.8
NGC 6334 G	17 20 25.2	−35 44 04	0.13	0.09	0.26	49	21	1.61	2.5
NGC 6334 H	17 20 31.2	−35 54 14	0.18	0.13	0.26	77	23	1.38	1.8
NGC 6334 I	17 20 35.3	−35 59 20	0.14	0.09	0.32	52	12	1.01	1.6
NGC 6334 J	17 20 39.2	−35 49 27	0.48	0.17	0.65	40	75	2.27	3.2
NGC 6334 K	17 20 48.4	−35 42 59	0.14	0.09	0.31	75	9	3.00	3.1
NGC 6334 L	17 20 54.4	−35 45 43	0.27	0.23	0.16	50	50	2.44	3.2
NGC 6334 M	17 20 57.3	−35 39 41	0.20	0.16	0.20	167	18	2.03	2.7
NGC 6334 N	17 21 32.5	−35 40 27	0.34	0.24	0.30	28	25	2.20	1.6
NGC 6357 A	17 24 43.7	−34 12 07	0.28	0.22	0.22	113	151	1.26	1.9
NGC 6357 B	17 24 46.7	−34 15 23	0.58	0.37	0.36	127	133	1.30	1.9
NGC 6357 C	17 25 33.3	−34 24 43	0.30	0.24	0.19	145	128	1.29	1.9
NGC 6357 D	17 25 34.3	−34 23 10	0.06	0.04	0.38	176	27	1.27	2.0
NGC 6357 E	17 25 47.9	−34 27 12	1.06	0.18	0.83	14	60	1.33	1.9
NGC 6357 F	17 26 02.2	−34 16 42	0.37	0.19	0.50	4	162	1.41	2.1
Eagle A	18 18 39.9	−13 47 41	0.12	0.12	0.01	139	46	0.88	1.4
Eagle B	18 18 42.2	−13 47 03	0.90	0.45	0.50	148	451	0.98	1.5
Eagle C	18 18 52.8	−13 46 43	0.25	0.10	0.60	103	20	0.96	1.4
Eagle D	18 18 57.3	−13 45 23	1.65	0.65	0.61	100	275	1.10	1.6
Eagle E	18 19 07.9	−13 36 28	0.08	0.07	0.07	13	20	2.40	2.6
Eagle F	18 19 12.8	−13 26 03	0.45	0.40	0.13	123	37	1.61	2.3
Eagle G	18 19 13.1	−13 33 45	0.20	0.10	0.51	102	20	2.59	3.7
Eagle H	18 19 15.0	−13 39 25	0.11	0.08	0.26	51	10	1.57	1.6
Eagle I	18 19 19.2	−13 36 34	0.44	0.16	0.64	37	18	1.64	2.4
Eagle J	18 19 29.7	−13 23 08	0.30	0.18	0.41	30	12	1.51	2.1
Eagle K	18 19 30.7	−13 45 30	0.06	0.06	0.04	50	6	1.42	2.0
Eagle L	18 20 02.8	−13 48 26	0.04	0.02	0.31	143	7	2.01	3.0
M 17 A	18 20 18.0	−16 14 11	0.07	0.07	0.03	41	10	1.92	3.2
M 17 B	18 20 19.5	−16 13 28	0.03	0.02	0.10	80	5	2.49	2.9
M 17 C	18 20 21.6	−16 10 35	0.14	0.13	0.09	135	31	1.70	2.5
M 17 D	18 20 22.5	−16 08 23	0.29	0.27	0.07	173	119	1.53	2.4
M 17 E	18 20 22.5	−16 09 52	0.09	0.07	0.18	68	13	1.56	2.5
M 17 F	18 20 22.8	−16 12 23	0.05	0.04	0.03	15	7	1.70	4.4
M 17 G	18 20 25.0	−16 11 35	0.04	0.04	0.05	94	6	1.55	3.9
M 17 H	18 20 25.6	−16 11 17	0.11	0.11	0.07	151	31	1.25	2.2
M 17 I	18 20 26.0	−16 09 34	0.17	0.17	0.02	24	48	1.45	2.1
M 17 J	18 20 27.8	−16 09 54	0.01	0.01	0.02	137	4	1.56	2.0
M 17 K	18 20 28.5	−16 10 58	0.19	0.16	0.14	177	78	1.41	2.2
M 17 L	18 20 29.9	−16 10 46	0.14	0.13	0.05	34	95	1.82	2.8
M 17 M	18 20 30.8	−16 03 16	0.20	0.16	0.21	94	31	1.83	2.5
M 17 N	18 20 31.4	−16 09 39	0.15	0.11	0.27	4	32	1.33	2.1
M 17 O	18 20 31.5	−16 11 22	0.09	0.07	0.18	57	26	1.42	2.4
Carina A	10 43 41.9	−59 35 49	0.19	0.09	0.55	115	11	1.03	1.6
Carina B	10 43 54.2	−59 33 02	0.70	0.45	0.35	191	193	0.97	1.5
Carina C	10 43 56.4	−59 32 54	0.22	0.18	0.17	80	106	0.94	1.4
Carina D	10 44 32.9	−59 33 42	0.60	0.30	0.50	56	46	0.90	1.4
Carina E	10 44 34.0	−59 44 08	0.28	0.27	0.05	137	36	0.84	1.4
Carina F	10 44 37.4	−59 26 03	0.64	0.40	0.36	108	38	0.94	1.5
Carina G	10 44 39.7	−59 44 26	10.15	1.88	0.81	177	741	0.91	1.5
Carina H	10 44 41.8	−59 22 05	0.30	0.23	0.23	209	49	0.83	1.3
Carina I	10 44 45.3	−59 20 07	0.23	0.18	0.24	70	18	0.80	1.3
Carina J	10 45 02.4	−59 45 50	0.50	0.46	0.09	163	65	0.96	1.5
Carina K	10 45 06.2	−59 40 21	0.35	0.25	0.28	94	36	0.89	1.5
Carina L	10 45 11.1	−59 42 46	0.49	0.36	0.27	229	58	0.91	1.5
Carina M	10 45 13.7	−59 57 58	0.48	0.33	0.32	122	52	0.95	1.6
Carina N	10 45 36.1	−59 48 20	0.39	0.38	0.03	84	37	1.51	2.1
Carina O	10 45 53.4	−59 56 53	0.18	0.17	0.09	110	25	1.15	1.8
Carina P	10 45 54.4	−60 04 32	1.17	0.59	0.50	60	118	0.87	1.5
Carina Q	10 45 55.2	−59 59 51	0.51	0.33	0.35	92	27	0.87	1.4
Carina R	10 46 05.4	−59 50 09	1.06	0.51	0.52	83	46	0.99	1.6
Carina S	10 46 52.7	−60 04 40	0.44	0.26	0.41	24	24	0.90	1.5
Carina T	10 47 12.5	−60 05 58	0.53	0.46	0.14	215	48	0.92	1.5
Trifid A	18 02 05.8	−23 05 53	1.88	0.36	0.81	103	87	1.53	1.3
Trifid B	18 02 23.1	−23 01 50	0.14	0.11	0.20	158	25	0.92	1.4
Trifid C	18 02 31.0	−23 00 40	0.65	0.52	0.20	92	82	0.87	1.3
Trifid D	18 02 39.7	−22 49 18	0.87	0.28	0.68	64	27	1.27	1.4
NGC 1893 A	05 22 40.3	+33 22 15	0.43	0.32	0.26	98	46	0.81	1.5
NGC 1893 B	05 22 46.6	+33 25 17	0.36	0.29	0.21	22	110	0.81	1.5
NGC 1893 C	05 22 47.4	+33 24 11	0.23	0.20	0.14	44	20	0.85	1.5
NGC 1893 D	05 22 49.7	+33 24 42	0.03	0.03	0.09	13	5	0.75	1.4
NGC 1893 E	05 22 49.7	+33 26 51	0.15	0.05	0.67	129	9	0.80	1.4
NGC 1893 F	05 22 53.6	+33 25 40	0.17	0.14	0.16	58	18	0.81	1.4
NGC 1893 G	05 22 53.8	+33 30 54	0.28	0.19	0.33	25	46	0.89	1.6
NGC 1893 H	05 22 57.5	+33 26 34	0.31	0.28	0.10	28	60	0.82	1.4
NGC 1893 I	05 23 03.3	+33 28 15	0.48	0.32	0.33	25	114	0.83	1.5
NGC 1893 J	05 23 08.5	+33 28 32	0.04	0.03	0.37	82	11	0.97	1.7

Notes. The ellipsoid components from the selected finite-mixture model. Column 1: the component name, labeled from east to west. Columns 2 and 3: celestial coordinates (J2000) for the ellipsoid center. Column 4 and 5: semimajor and minor axes of the ellipsoid core region in parsec units. Column 6: ellipticity. Column 7: orientation angle of the ellipse in degrees east from north. Column 8: number of stars estimated by integrating the model component out to four times the size of the core. Column 9: median J − H color index as a measure of absorption. Column 10: median X-ray median energy as a measure of absorption.

Download table as:  ASCIITypeset images: 1 2 3

Some subclusters are dominated by X-ray selected sources, while others are dominated by IR-excess selected sources. For example, in the Eagle Nebula region, the main star cluster is dominated by the former while the embedded subclusters to the north and east are dominated by the later. Thus, the combined selection criteria of the MYStIX survey help reveal a fuller picture of structures. But observed star counts for individual subclusters can not be directly compared without correction for detection sensitivity effects.

The assignments of individual MPCM stars to subclusters (or to the unclustered young stellar component) is shown in Figure 3 and listed in Table 3. Table 3 also provides X-ray and near-IR properties—these are used in Table 2 to calculate the median J − H and median X-ray ME values. In low-density subclusters, the number of stars inferred from the model, N_{4, obs}, may not precisely equal the number of stars assigned to that subcluster if there is spatial overlap with neighboring subclusters leading to ambiguities in assignments.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (17 panels)

Tarball of figure set images

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (17 panels)

Tarball of figure set images

One can see some strong successes in the decomposition of stars in to isothermal ellipsoids, but also some limitations. In the Orion region (Figure 1) the subclusters disentangle layered structures in the region, including the BN/KL cluster (Becklin & Neugebauer 1967; Kleinmann & Low 1967), the rich ONC, and stars in the north–south OMC-2/3 filament (Megeath et al. 2012). Furthermore, two ellipsoids are needed to fit the ONC—one corresponding to a dense "core" and the other corresponding to a less dense "halo." This result agrees with Henney & Arthur (1998), who find that Orion proplyd and non-proplyd stars are strongly peaked around θ¹ Ori C with statistically significant excess above the surface densities expected from the King profile fit by Hillenbrand & Hartmann (1998). On the other hand, certain structures are not modeled: no distinction is made between the unabsorbed ONC members and the Orion Molecular Cloud stars (OMCs; Feigelson et al. 2005), the OMC-1S subcluster (Lada et al. 2004) is not identified, and stars from the ι Ori cluster (Alves & Bouy 2012) may contribute to the subclusters and unclustered component.

On occasion, the minimum AIC criterion requires a subcluster that appears very small with few stellar members; for example, subclusters B, E, F, J, K in Rosette and subcluster J in M 17 have fewer than five stars. This can be a consequence of poorly constrained core radii (see discussion in Section 6.1). One of these sparse structures is Subcluster B in M 17 which is associated with a well-studied embedded ultracompact H ii region Kleinmann & Wright (1973). Nevertheless, these subclusters often correspond to clumps of stars that appear real to visual inspection. On the other hand, occasionally model components contain large numbers of stars that do not appear to be discrete subclusters. This includes subcluster G in Carina and subcluster E in NGC 6357.

Getman et al. (2014b, their Figures 7–22) show subclusters superimposed on 500 μm Herschel-SPIRE maps of the molecular clouds in the star-forming region. Some of the sparser structures are coincident with molecular cloud cores, providing corroborating evidence that they are real embedded subclusters. Some examples include, NGC 2264 (Subclusters G, I, J, K, L, and M), Rosette (Subclusters J, M, N, O), Lagoon (Subclusters B, K), DR 21 (Subclusters H, F, E, D, C), NGC 6334 (Subclusters M, L, J, E, D), Eagle (Subclusters E, G, H, L), and NGC 1893 (Subclusters J and G).

The properties of the best-fit model for each region are provided in Table 4, including the number of subclusters, the observed surface density of the unclustered young stars, the minimum and maximum of the fitting error residuals and the kernel size used to calculate these residuals, and a χ² goodness-of-fit statistic.

Table 4. Best-fit Models

Region	Nsubclust	ΣU	Residual Map			Goodness of Fit
			Kernel	Min	Max	χ2	dof	$P_{\chi ^{2}}$
		(stars pc−2)	(pc)	(stars pc−2)	(stars pc−2)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
Orion	4	30	0.076	290 (4%)	−270 (4%)	47	60	0.23
Flame	1	4	0.12	80 (4%)	−60 (3%)	21	23	0.78
W 40	1	20	0.17	38 (3%)	−26 (2%)	24	17	0.22
RCW 36	3	7	0.11	96 (3%)	−57 (2%)	17	17	0.82
NGC 2264	13	4	0.28	16 (3%)	−20 (4%)	50	46	0.62
Rosette	15	0.7	0.72	3.8 (7%)	−2.5 (5%)	126	62	<0.01
Lagoon	11	0.2	0.33	14 (3%)	−9.1 (2%)	61	57	0.61
NGC 2362	2	0.3	0.58	3.5 (7%)	−2.0 (4%)	9	9	0.86
DR 21	9	3	0.39	13 (6%)	−11 (5%)	72	31	<0.01
RCW 38	4	0.08	0.30	13 (1%)	−8.6 (1%)	25	25	0.88
NGC 6334	14	1	0.50	4.9 (4%)	5.9 (5%)	43	37	0.41
NGC 6357	6	0.01	0.38	14 (3%)	−12 (3%)	120	73	<0.01
Eagle	12	1	0.41	14 (6%)	−9.1 (4%)	75	66	0.42
M 17	15	0.2	0.24	25 (3%)	−24 (3%)	75	73	0.80
Carina	20	0.01	0.59	12 (4%)	−3.5 (1%)	200	145	<0.01
Trifid	4	0.3	0.8	1.5 (3%)	−1.2 (3%)	11	12	0.98
NGC 1893	10	0.06	0.64	3.0 (3%)	−7.0 (7%)	51	34	0.06

Notes. Properties of the best-AIC models. Column 1: region name. Column 2: number of subclusters in the model. Column 3: observed surface density of the uniform model component. Column 4: Gaussian kernel size used to smooth residuals. Columns 5–6: maximum and minimum residuals (as percent of peak observed surface density). Columns 7–9: goodness-of-fit test results: Pearson's χ² value, degrees of freedom, and the probability, respectively.

Download table as:  ASCIITypeset image

The number of subclusters per region ranges from 1 (Flame and W 40) to 20 (Carina). The variation in this number is an indication of the vast diversity in spatial structure. However, this number also relates to the size of the field of view. The fractal distribution of star formation throughout the Galaxy makes it difficult to choose a sample field of view that adequately characterizes a region. Often clusters are located near each other and it is unclear whether they should be treated separately or as subclusters, as is the case with the three rich clusters in NGC 6357 and the Tr 14-15-16 clusters in the Greater Carina Nebula. Selection effects in our survey of young stars also influence the number of inferred subclusters.

The goodness-of-fit statistics for the tessellated observed versus model MPCM counts in Table 4 suggest that in most cases, the finite mixture models are neither overfitted (χ² ≪ dof) nor poorly fit (χ² ≫ dof). Poor fits are an indication that structures are not simple isothermal ellipsoids, which is not surprising in elaborate, complex structures seen in regions like Rosette, DR 21, NGC 6357, and Carina.

Figure 4 (left) shows a map of the smoothed fitting residuals for the Carina region (residual maps for the other regions are provided by the electronic version of this article) that were smoothed by a 0.59 pc Gaussian kernel. Red and blue colors indicate positive and negative residuals in units of stars pc⁻², respectively. For Carina the maximum residual is located southwest of Tr 14, between Subclusters A and B—this indicates an asymmetry to Tr 14 due to an excess of stars at the edge of the Northern Cloud in this star formation complex. Otherwise, the positive and negative residuals are low and more or less evenly distributed.

In all the MYStIX regions, minimum/maximum error residuals for the fits are typically <10% of the peak stellar surface density in a cluster. Salient features of residual maps in select regions are described here.

•   
  Orion Nebula. The minimum residual is located south of the ONC cluster and may indicate some deviation from a true isothermal ellipsoid shape. There is also a large negative patch east of the ONC, which was noted as a void by Feigelson et al. (2005).
•   
  Flame Nebula. There is a large positive residual to the northwest of the cluster, which corresponds to a group of highly absorbed stars.
•   
  W 40. The residuals are nearly zero in the cluster core, but there is a positive residual to the southeast, composed mostly of IR-excess selected stars, with a population too small to be accepted as a subcluster by the AIC statistic.
•   
  Rosette Nebula. The residuals are higher amplitude than in other MYStIX regions and may indicate that our finite mixture model fit is inadequate for this complex region. The χ² value for this fit is also poor.
•   
  Lagoon Nebula. Residuals are low indicating an overall good fit. However, there is a long positive residual stretching from the west to the south of NGC 6530. This residual coincides with the edge of the molecular cloud and a bright rimmed cloud; it may represent an area of triggered star formation.
•   
  DR 21. The large blue residual along the chain of subclusters in DR 21 indicates problems with this fit. This is also implied by the poor χ² statistic for this region. These results suggest that the deeply embedded clusters in DR 21 are distributed along the molecular filament and not yet fully collected into isothermal ellipsoids.
•   
  RCW 38. The small residuals suggest a good match between the data and the "core–halo" model, especially in the center of the cluster. The maximum residual is located north of the cluster and does not appear to be associated with any identified subcluster.
•   
  NGC 6357. A strong negative residual is located in the center of the southern complex of stars (G353.1+0.6; Townsley et al. 2014). The best AIC is obtained when three ellipsoids were used to fit this complex, although three distinct components are not obvious from a visual examination of the asymmetrical G353.1+0.6 cluster.
•   
  Eagle Nebula. The main cluster in the Eagle Nebula is composed of two overlapping subclusters (A and B). There are large positive residuals at the edges of the cluster; nevertheless, different configurations of the two subclusters do not produce lower residuals, so the structure can not be fully represented with isothermal ellipsoids.

In the right panels of Figure 4, the surface density obtained from the nonparametric adaptive smoothing (x axis) is plotted against the surface density obtained from the parametric finite mixture mode (y axis) where each point corresponds to a pixel in the residual map. These generally agree to within 0.1 dex. The local maxima from the surface density maps are typically somewhat higher in the models than in the smoothed maps; this is an artificial effect of the smoothing algorithm. For Carina, the agreement is good at surface densities above ∼3 stars pc⁻². Deviations at lower surface densities may reveal differences between the modeling of unclustered young stars with a constant surface density and the true surface-density variations of the unclustered stars.

The size of subcluster cores varies from 0.01 pc to >2 pc. To investigate subcluster sizes, we use the geometric mean of semimajor and semiminor axes to define a characteristic radius. For the subclusters from Table 2, the distribution of characteristic radii is shown in Figure 7. The distribution is fit well by a lognormal distribution with a mean 0.17 pc and standard deviation 0.43 dex; that is, 68% of core radii lie between 0.06 and 0.45 pc. The Anderson–Darling test for normality shows no statistically significant deviation from this distribution (p > 0.05). The peak of the distribution is similar to the ONC core radius of ∼0.2 pc found by Hillenbrand & Hartmann (1998).

Zoom In Zoom Out Reset image size

The tails of this distribution are not well-defined because the smallest radii are seen in sparse subclusters without clear core regions, while the largest radii arise from large-scale overdensities that do not appear to be distinct subclusters. There may be a bias against the detection of subclusters with very large or small sizes. Detection of large subclusters is limited by the field of view and its stars may instead be assigned to the flat distributed star component in our statistical model. Very small subclusters will also be difficult to detect if they lack sufficient numbers of stars to significantly increase the likelihood of the model. A trend of subcluster size increasing with MYStIX region distance suggests this bias is present. On larger scales, in regions with where the young star surface density is low, distinct subclusters may be merged together to create larger artificial subclusters; deeper X-ray and IR surveys may result in more subclusters. Nevertheless, the most common subcluster size (∼0.2 pc) is consistent for MYStIX regions at all distances.

Figure 8 is a scatter plot of cluster core radius versus two measures of absorption along the line of sight, X-ray median energy ME and J − H for the cluster. A negative relation can be seen between these quantities, which is emphasized by the gray LOWESS non-parametric regression. LOWESS is a well-established regression technique that combines robust least squares local regression technique with a k-nearest neighbor meta-model (Cleveland 1979). The J − H range corresponds approximately to 2 < A_V < 10 mag, and the X-ray median energy range corresponds to A_V < 1 to A_V ≃ 50 mag (Getman et al. 2010). Although there is much scatter in these relations, nonparametric correlation hypothesis tests using the Spearman rank coefficient (and confirmed using bootstrap resampling) give significant levels p < 0.001 for the radius–ME relation and p < 0.01 for the radius–J − H relation.

Zoom In Zoom Out Reset image size

Some of the highly absorbed subclusters with small radius are well-known embedded star forming regions including the BN/KL cluster (subcluster A in Orion) and the Kleinmann–Wright Object cluster in M 17 (subcluster B in M 17). In particular, subclusters with median ME > 2.5 keV (corresponding to A_V > 17 mag) are rarely larger than 0.2 pc. This phenomenon is more evident in the histogram of core radii, which shows that clusters with median X-ray median energies greater than 2.5 (shaded gray region in Figure 8) typically have smaller radii than less embedded clusters. The peak of their lognormal distribution is at 0.08 pc, and the Anderson–Darling test shows no statistically significant deviation from a lognormal distribution (p > 0.05).

Clusters lose their gas as they age, so the J − H index for absorption can serve as a proxy for age—this is directly demonstrated by Getman et al. (2014b). Thus, the negative size–absorption relation can be viewed as a positive size–age relationship suggesting that clusters expand as they age. Gieles et al. (2012) find that the combination of multiple clusters at different stages of expansion can produce a distribution of stellar surface densities similar to what is seen in the MYStIX regions.

This effect may explain some of the differences in cluster sizes between different star formation regions. For example, DR 21 and NGC 1893 are both composed of linear chains of subclusters, but the subclusters in NGC 1893 are larger. This makes sense in light of the size–absorption relation because the younger DR 21 subclusters are deeply embedded in a molecular filament, while the older NGC 1893 subclusters are lying in a bubble and are lightly absorbed. Marks & Kroupa (2012) identify a typical half-mass radius for young embedded clusters of less than a few tenths of a parsec; for MYStIX subclusters with ME > 3 the mean core radius is considerably smaller than this around 0.04 pc (Figure 8).

We note that the typical size of star clusters is similar to the width of filaments found by Herschel, which appear to be ubiquitous in MSFRs and involved in the cluster formation process (André et al. 2010). The variation in size and absorption of subclusters within a single star-forming region can be interpreted as an indication of non-coeval clustered star formation.

Figure 9 (left) shows a histogram of MYStIX cluster ellipticity. The measured ellipticities will be affected by projection, so they are lower limits on intrinsic ellipticity. Ellipticities range from 0 to ∼0.8, with a higher fraction having the lower ellipticities. Nevertheless, there are a significant number of subclusters with substantial ellipticities. Figure 9 (right) shows no relation between ellipticity and cluster core radius; this is confirmed using Spearman's rank correlation test.

Zoom In Zoom Out Reset image size

Dynamical studies have shown that highly elliptical cloud configurations are unstable to gravitational fragmentation (Burkert & Hartmann 2004). Note that the most extreme ellipticities with ≃ 0.8 in our model fits may not distinct subclusters: for example, subcluster E in NGC 6357 appears to be part of a larger asymmetric complex, subcluster G in the Carina Nebula is a huge structure along the spine of the complex likely including unclustered stars and fitting residuals from other subclusters, and subcluster A in Trifid is one end of a larger filamentary structure that was truncated by the Chandra field of view.

The ellipticity distribution for MYStIX subclusters has some resemblance to the distribution of ellipticities in simulated star-forming regions found by Maschberger et al. (2010). There, ellipticities range from 0.0 to 0.82, lower ellipticities are more common, and the peak ellipticity is 0.33. In the simulated clusters, significant ellipticities are an effect of subcluster mergers, and this may also be the case for the MYStIX clusters. The frequent presence of high ellipticities in MYStIX subclusters suggests that many have not reached dynamical equilibrium.

The numbers of subclusters and their arrangements vary significantly from region to region. Heuristically, the clustering morphologies are divided into four classes: isolated clusters, core–halo clusters, clumpy clusters, and linear chains of clusters (Section 5). The dominant subcluster in a region (i.e., the subcluster with most stars, N_{4, obs}) often sits at the center of a group of subclusters (e.g., subcluster L in M 17), but this is not always true (e.g., subcluster B in the Eagle Nebula). In regions with linear chain structure there may be no subcluster that is particularly dominant. In contrast, regions with simple or core–halo structure may be much more centrally concentrated.

In DR 21, NGC 2264, and NGC 6334, the linear chain cluster structure maps onto molecular filaments in the regions, and it is clear that the star cluster structure was at least partially inherited from the filamentary chain of molecular cloud clumps. More detail about correspondence between subcluster locations and molecular cloud clump locations will be given in a later MYStIX study. Thus, regions with linear chain structure have not evolved much from their initial state at the time of the birth of stars. The more centrally concentrated clumpy, core–halo, and simple cluster structures are more likely to have dynamically evolved and, in cases where the subclusters are well fit by single isothermal ellipsoids, may have achieved some dynamical relaxation.

The morphological classification suggests an evolutionary progression from linear chain structure to clumpy structure to core–halo structure to simple cluster structure. Hierarchical cluster formation, often in molecular filaments, would evolve into dynamical relaxed unimodal structures. However, the MYStIX regions indicate that this progression is not always followed. For example, the linear chain region NGC 1893 is almost entirely unembedded, and the chain of subclusters in NGC 6334 exhibit a wide range of absorptions.

The distribution of subcluster properties show trends common to all the regions. A lognormal distribution of subcluster size peaks at 0.17 pc, and many subclusters have ellipticities of ⩽ 0.3, but few have ellipticities greater than 0.5. No relationships between ellipticity and size are seen.

An additional property of subclusters, their absorption, can help link these properties to cluster formation history and evolution. Line-of-sight absorption is measured both by the near-IR color index J − H and the X-ray median energy indicator (Getman et al. 2010). While absorption may be caused by overlaying molecular material, in most cases it represents local cloud cores within which the youngest clusters reside. In complex MYStIX star-forming regions, a wide range of absorption can be seen for the different subclusters, indicating these subclusters are not coeval. A statistically significant anti-correlation between absorption and size of a subcluster is found, most easily explained if both absorption and size are related to age. This issue is pursued in the accompanying study by Getman et al. (2014b) where spatial-age gradients in MYStIX regions are found based on a new age indicator for individual pre-main-sequence stars. A similar relationship was found for subclusters within the Rosette Molecular Cloud by Ybarra et al. (2013) who infer a gas removal e-folding timescale of 0.4 Myr.

Cluster expansion is expected to produce a wide range of stellar surface densities seen in many star-forming regions (Gieles et al. 2012). Early cluster expansion arises principally from the loss of the gravitational potential of the molecular material, but stellar dynamics can also produce stellar halos around dense cores through three-body interactions with hard binaries and gravothermal core collapse associated with mass segregation (Giersz & Heggie 1996). This issue will be discussed further in Paper II when intrinsic populations, rather than observed population subject to different sensitivity effects, are derived.

The anti-correlation between radius and absorption provides evidence that subcluster expansion is important part of subcluster dynamical evolution. From numerical simulations, Moeckel & Bate (2010) find that even the very dense simulated cluster produced by Bate (2009), with a half-mass radius of ∼0.05 pc, grew to >2 pc over a period of ∼2 Myr. The typical small size that we find for the embedded clusters indicates that the high initial densities from the simulation in Bate (2009) are realistic, which can make competitive accretion and dynamical interaction between protostars much more important than if clusters form ∼1 pc in size.

Simulation of cluster expansion after gas expulsion by Goodwin & Bastian (2006) show a relationship between star formation efficiency and the age–radius relationship. They find that, after 2.5 Myr, core radii will grow by a factor of ∼1.5 for a star-forming efficiency of 60% and a factor of ∼3 for a star-forming efficiencies of 10%–30%. Larger expansion factors are found in models of very rich clusters by Banerjee & Kroupa (2013). For comparison, MYStIX results indicate that subcluster core radii start out at ∼0.08 pc and grow to the typical ∼0.2 pc for the typical unabsorbed MYStIX subclusters (a factor of more than two expansion). The stars in MYStIX regions are generally <5 Myr old, but more precise age determination is necessary to determine which age–radius relation provides the better fit. Pfalzner (2011) find a trend consistent with a star formation efficiency of ∼30% using the clusters in Lada & Lada (2003).

The ellipticities of MYStIX subclusters could be inherited from the structure of the parental molecular cloud or could arise from subcluster mergers. The subcluster ellipticity distribution resembles those found by numerical simulations of merging subclusters (Maschberger et al. 2010, their Figure 10(b)) in which merger products have a distribution of ellipticities, including a tail at high ellipticities.

The cluster morphological classes may correspond to structures seen in star formation simulations, such as those performed by Bate et al. (2003), Bate (2009) or Parker & Meyer (2012). Simulations often show multiple pockets of star formation occurring along gas filaments as the cloud collapses, which can end up merging and dynamically relaxing as the simulation progresses. Clumpy clusters may indicate early stages of partial merging, core–halo structures may be an intermediate stages of mergers, and isolated clusters with relaxed structure a final stage. Getman et al. (2014a) add an important constraint to such a scenario: in at least two MYStIX regions, the dense core of the cluster has younger stars than the outer regions. This observational result requires some special behavior, such as continual feeding of gas toward the cluster center or dispersal of older core stars into the halo.

The simple clusters may represent those that have completed subcluster merging and dynamical relaxation. Subcluster morphology can also influence these rates. Smith et al. (2011) find that the cluster formation rate is related to a filling factor measuring the fraction of volume containing subclusters. When this filling factor becomes high, tidal interactions accelerate subcluster mergers and formation of a simple unimodal cluster is promoted. This phenomenon may be occurring today in crowded, clumpy MYStIX regions like M 17.

The core–halo structures in 5 out of 17 MYStIX regions may also be a result of subcluster mergers. Bate (2009) report a hydrodynamical and N-body simulation of a MSFR resulting in a cluster with a half-mass radius of 0.05 pc and a halo of stars extending out ∼0.4 pc. Smith et al. (2011) found that gravitational interactions between subclusters of young stars in simulated star-forming environments can scatter stars out of the clumps, which can lead to a dense core surrounded by a less dense halo of scattered stars. Maschberger et al. (2010) also identify a halo of scattered low-mass stars around the dominant cluster, albeit with a surface-density–radius relation opposite of the effect we find with concentric ellipsoids. Pfalzner & Kaczmarek (2013) report that two-body encounters can produce a population of cluster stars with high ellipticities which will spend most of the time in a halo. The accompanying paper by Getman et al. (2014a) describes star formation scenarios that may account for core–halo structures. All of these issues relating to cluster formation and MSFR star formation histories will be pursued in MYStIX papers by Getman et al. (2014b) and Paper II, where quantitative measures of subcluster intrinsic populations and cluster ages will be combined with the structural results obtained here.