The Spatial Structure of Young Stellar Clusters. I. Subclusters
Michael A. Kuhn, Eric D. Feigelson, Konstantin V. Getman, Adrian J. Baddeley, Patrick S. Broos, Alison Sills, Matthew R. Bate, Matthew S. Povich, Kevin L. Luhman, Heather A. Busk, Tim Naylor, Robert R. King
AAccepted, 13 March, 2014
The Spatial Structure of Young Stellar Clusters. I. Subclusters
Michael A. Kuhn ∗ , Eric D. Feigelson , Konstantin V. Getman , Adrian J. Baddeley ,Patrick S. Broos , Alison Sills , Matthew R. Bate , Matthew S. Povich , Kevin L.Luhman , Heather A. Busk , Tim Naylor , Robert R. King ABSTRACT
The clusters of young stars in massive star-forming regions show a wide rangeof sizes, morphologies, and numbers of stars. Their highly subclustered structuresare revealed by the MYStIX project’s sample of 31,754 young stars in nearby sitesof star formation (regions at distances < * [email protected] Department of Astronomy & Astrophysics, Pennsylvania State University, 525 Davey Laboratory, Uni-versity Park PA 16802 School of Mathematics and Statistics, University of Western Australia, 35 Stirling Highway, CrawleyWA 6009, Australia Department of Physics, McMaster University, 1280 Main Street West, Hamilton ON, L8S 4M1, Canada Department of Physics and Astronomy, University of Exeter, Stocker Road, Exeter, Devon, EX4 4SB,UK Department of Physics and Astronomy, California State Polytechnic University, 3801 West Temple Ave,Pomona, CA 91768 a r X i v : . [ a s t r o - ph . GA ] M a r Subject headings: methods: statistical; open clusters and associations: general;stars: formation; stars: pre-main sequence; H ii regions; ISM: structure
1. Introduction
Massive star-forming regions (MSFRs) are a major, and perhaps the dominant, modeof star formation in the Galaxy (Lada & Lada 2003; Fall et al. 2009; Chandar et al. 2011).Most of the young stars in these massive complexes are clustered (Clarke et al. 2000; Lada& Lada 2003), but a substantial, spatially distributed population of young stars also exists(Evans et al. 2009; Feigelson et al. 2011). Open questions on how young stellar clustersform and evolve in time in these harsh environments are outlined in the overview of theMassive Young Star-Forming Complex Study in Infrared and X-ray (MYStIX Feigelson etal. 2013). This multifaceted project starts by characterizing the young stellar populations of20 OB-dominated MSFRs at distances d < . . Hillenbrand & Hartmann (1998) found that a King (1962)profile describes the radial surface density profile of the ONC, with ρ = 2 − × M (cid:12) pc − and core radius r = 0 . − . χ Per, and young stellar clusters in the LargeMagellanic Clouds (Wang et al. 2008; Kuhn et al. 2010; Wang et al. 2011; Sung & Bessell2004; Harfst et al. 2010; Bragg & Kenyon 2005; Mackey & Gilmore 2003). Pfalzner et al.(2012) note that these clusters tend to have large ratios of half-mass radius to core radius,so the truncated King profile is a better model than the Plummer sphere.In addition to the main subclusters in regions, smaller groups of stars can be seen outsidethe main clusters, or even as subclusters within the main clusters. It is unclear whether arelaxed surface density profile is a good model for subclusters, in part because young stellarclusters are generally not expected to have had time to dynamically relax through two-bodyinteractions and molecular gas may contribute significantly to the gravitational potential.Nevertheless, similar surface density profiles have been used by Smith et al. (2011), who usePlummer sphere subclusters in their simulations, and found by Maschberger et al. (2010) inthe radial density profile of simulated subcluster mergers. The definitions of “young stellar cluster” in the star-formation literature are divergent (for an alternativedefinition cf. Portegies Zwart et al. 2010). Here, we use “cluster” in a statistical sense (subclusters are thecomponents of the finite mixture model in Section 3).
2. Stellar Samples: MYStIX Probable Complex Members
The regions included in this study, listed in Table 1, are nearby sites of high-mass starformation that are included in the MYStIX sample (Feigelson et al. 2013). Each of theseregions has a population of hundreds or thousands of young stars in the MYStIX sample.Some fields appear relatively simple with a single, dominant cluster, like the ONC or W 40;while others have complex structures, like the Rosette nebula which contains the dominantNGC 2244 cluster plus smaller, embedded clusters within the Rosette molecular cloud. Thefield of view must also be taken into consideration: for example, a single
Chandra exposureof the Orion Nebula captures only one major cluster on a ∼ Chandra exposures of the more distant NGC 6357 complex captures three richclusters on a ∼
20 pc scale. The molecular clouds vary from region to region as well —NGC 2362 has no molecular cloud left, others like Orion or NGC 1893 have blown a cavityin their molecular cloud around the star cluster, while others like NGC 2264 and DR 21 stillhave deeply embedded clusters. V -band extinction typically ranges from 2 to 10 mag, butthe obscuration may be highly non-uniform with A V >
100 in the densest molecular coresand filaments. There may also be significant optical and IR nebulosity, and the polycyclic-aromatic-hydrocarbon (PAH) emission (seen with the
Spitzer Space Telescope ) can trace theinterface between the ionized H II region and the molecular cloud.These regions are covered in the near-IR by the UKIRT Infrared Deep Sky Survey(King et al. 2013) or Two-Micron All Sky Survey (Skrutskie et al. 2006), in the mid-IRby
Spitzer
IRAC observations (Kuhn et al. 2013b), and in the X-ray by the ACIS-I arrayon the
Chandra
X-ray observatory in the 0.5-8.0 keV band (Kuhn et al. 2013a; Townsleyet al. 2014). Datasets are also obtained from projects by Getman et al. (COUP; 2005),Benjamin et al. (GLIMPSE; 2003), Kuhn et al. (2010), and Majewski et al. (2007). Themultiwavelength point-source catalogs are combined with a list of spectroscopic OB starsfrom the literature, and sources are classified probabilistically to obtain a list of MYStIXProbable Complex Members (MPCM; Broos et al. 2013) Three types of sources make up theMPCMs: objects that are detected in the X-ray, objects with IR excess consistent with youngstellar object (YSO) disk/envelope emission, and the spectroscopic OB stars. The MPCMsinclude sources with disks and/or envelopes and disk-free complex members, high-mass andlow-mass members, and clustered and unclustered members.In the MYStIX project, data reduction and classification are performed as uniformlyas possible across the study to facilitate comparative analysis between different regions.The fields of view of the study are the intersections of the near-IR, mid-IR, and X-raydata. Nevertheless, differences among the observations (e.g. exposure times, coverage) anddifferences among the regions (e.g. distance, obscuration, nebulosity, and crowding) and 6 –biases in the selection of the different types of MPCMs are present. To improve spatialuniformity of the study, particularly variations in the X-ray sensitivity across a region, ouranalysis here is restricted to a subset of MPCMs.Table 1 lists the MYStIX regions analyzed here. Column 3 gives the X-ray flux limitwith uniform X-ray coverage. Columns 4-7 give the number of MPCMs used for the spatialanalysis, including the overlapping subsets selected by each criterion. 7 –Table 1. MYStIX Targets for Subcluster AnalysisNumber of Sources in SampleRegion Distance log F X, limit Total X-ray IR-excess OB(kpc) (ph s − cm − ) (stars) (stars) (stars) (stars)(1) (2) (3) (4) (5) (6) (7)Orion Nebula 0.414 − − − − − − − − − − − − − − − − − . − . 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¨udel & Naz´e 2009; Feigelson 2010). X-ray luminosities are typically L X ∼ − − − 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 ∼ − L bol . Thus, X-ray observations can be used todiscriminate young stars both disk-bearing and disk-free in star-forming regions from olderGalactic field stars. Although most point sources detected in the Chandra observation dis-cussed here are young stars, some X-ray sources are non-member contaminants that includeactive galactic nuclei and foreground/background main-sequence stars (Getman et al. 2011).Here, members and non-members are distinguished by Bayesian supervised learning usingphotometric, 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 complexenvironments of MSFRs. Optical surveys are limited by absorption, which may exceed100 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 nebularemission of the H II region, which is largely absent in the X-ray. Often young stellar clustersare centered near the regions with the highest IR nebulosity, so the X-ray observations canbe more effective for identifying low mass stars in areas with strong gas and dust emissionnear OB stars. In addition, the brightness contrast between T-Tauri and OB stars is not sodramatic in the X-ray as in optical and in IR bands (e.g. Cohen et al. 2008). Thus, low-massstars can be detected in the X-ray, even when they are projected near OB stars in the denseyoung stellar cores. Finally, the low number of resolved field stars and the high resolutionof the Chandra mirrors reduce X-ray source confusion in the cluster centers to a low levelfor the MYStIX MSFRs.For X-ray selected complex members, source detection is limited by apparent photonflux in the 0 . − . Chandra mosaic, and telescope effects such asoff-axis mirror vignetting and degradation of the point-spread function. Finally, comparisonbetween MYStIX regions must judiciously consider the different
Chandra exposure times anddistances to the regions. In the present study, we compensate for the intra-region variationsin sensitivity. Inter-region differences in sensitivity are treaded in the forthcoming PaperII (M. A. Kuhn et al., in preparation). Readers are cautioned about comparing stellar 9 –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 and canproduce artificially clustered sources; this has been called the “egg-crate effect” when manyfields are combined into a mosaic. This problem can be seen clearly in the X-ray sourcedistribution 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 ∼ Chandra field compared to on-axis portions. For our spatialanalysis, we flatten these inhomogeneities in sensitivity by selecting 11,798 X-ray selectedsources that have apparent 0 . − . Chandra pointing in the ACIS-I mosaic that has the shortest exposure time, X-ray selected MPCMsare stratified by off-axis angle at 0.3, 5.1, 6.3, 7.5, 8.3, and 12.3 arcmin. Their apparent0 . − . F X, limit , are given by Column 3 of Table 1, and the number of X-ray selectedMPCMs remaining in the flattened sample is given in Column 5. Typically, ∼
50% ofthe X-ray selected MPCMs (Broos et al. 2013, their Table 4) are included in the flattenedsample .The uniform X-ray photon flux limit allows structures to be characterized independentlyof the Chandra pointing history. But other spatial biases are present. When variations ingas absorption are present, a fixed X-ray photon flux limit corresponds to a higher intrinsicluminosity limit in high-absorption regions. Deeply embedded regions may this be deficientin X-ray MPCMs, but often this is compensated by the higher fraction of IR-excess sources.Another bias may arise due to mass segregation. X-ray luminosity is correlated to stellarmass (e.g. XEST; G¨udel et al. 2007), so detection is not independent of mass, and MPCMsmay more closely reflect the spatial distributions of the more massive population. A similarbias 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 The young stars used in this study include MPCMs with F X > F X, limit , IR-excess, or spectral classifica-tion as OB stars. If an X-ray selected MPCM is excluded from the “flattened sample”, it may be includedin the spatial analysis if it meets either of the other criteria.
10 –(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 theyoung stellar photosphere. The spectral energy distributions (SEDs) may show an IR-excessin the K -band, but it is strongest in the Spitzer mid-IR bands. Povich et al. (2013) catalog theMYStIX InfraRed Excess Sources (MIRES) with SEDs consistent with disk-envelope-bearingyoung stars. Povich et al. provide several filters to reduce likely background sources; non-member contaminants include reddened stellar photospheres, star-forming galaxies, activegalactic nuclei, (post-)asymptotic giant branch stars, shock emission, unresolved knots ofPAH emission, and background young stars. We use here only MIRES stars that theydesignate SED FLG = 0 as likely young stars with protoplanetary disks. For the OrionNebula 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 tovery 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 importanceof 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 globuleon the periphery of an H II region.
Spitzer
IRAC observations of MSFRs suffer from high, non-uniform nebulosity, primar-ily from polycyclic aromatic hydrocarbon 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 it has several-times better on-axis spatialresolution than
Spitzer . These mid-IR effects may lead to severely decreased MIRES sen-sitivity in regions with the highest nebular contamination and stellar density. We do notattempt to correct these effects by removing weak sources, as we do for the X-ray selectedMPCMs, because any attempt to apply a uniform flux limit across the field of view wouldleave only the few brightest stars.The IR-excess selected MPCMs are also biased against widely distributed stars, becauseof the spatial filtering used to remove probable contaminants. This may enhance slight over-densities in source counts, and also interferes with attempts to characterize the distributedpopulations (Feigelson et al. 2013). 11 –
Spectroscopic catalogs of OB-stars are available for many regions (e.g. Skiff 2009, andreferences therein), and all known OB stars in MYStIX regions are included in our spatialanalysis because of the large roll they play in clustered star formation. However, these listsmay 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 luminositiescomparable 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 emittingOB stars (H. Busk et al., in preparation).
The statistical sample of stars used for the analysis of spatial distributions contains16,608 MPCM out of 31,754 MPCMs in the entire MYStIX project. The resulting samplesrange 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 areaffected by ancillary effects such as region distance, obscuration, nebulosity, and observationexposure times. These observational selection effects will be treated in Paper III; the presentstudy is confined to the observed MPCM population.Nearly all regions have more X-ray selected stars than IR-excess stars. The major clus-ters are typically dominated by the X-ray selected sources, but peripheral, highly-obscuredgroupings 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 ofthe clusters are dominated by X-ray selected members.
3. Modeling Spatial Distributions of Young Stars with Isothermal Ellipsoids
The star-forming regions investigated here include many famous star clusters, but theMYStIX catalogs still reveal new structures in the spatial distributions of the young stellarmembers, which are complicated and multimodal. The structure is sometimes dominatedby a single monolithic cluster, sometimes a dominant cluster with substructure, sometimesseveral clusters, and sometimes a dominant cluster with secondary subclusters. To model allof the spatial structures in an unbiased fashion, we need a statistical method that can findall of the clusters and subclusters in an objective fashion. 12 –A variety of methods exist to perform unsupervised cluster analysis including non-parametric methods (such as kernel density estimation, the k -nearest neighbor algorithm,and graph-based algorithms) and the parametric methods that we use for this study (e.g.Everitt et al. 2011). The nonparametric methods require the choice of a critical value forthe algorithm; for example, the number of nearest neighbors (Rom´an-Z´u˜niga et al. 2008),threshold minimum spanning tree branch length (Gutermuth et al. 2009), and a thresholdsurface density in a kernel density estimator (Feigelson et al. 2011). Nevertheless, the “NoFree Lunch Theorem” states that there is no general way to determine whether a particu-lar cluster solution is better than another found by a different method in the multivariateunsupervised classification problem (Wolpert & Macready 1997). The advantages and dis-advantages of our parametric method are discussed in Section 3.6.Our statistical procedure for identifying subclusters in MPCM stars is based on finitemixture models (McLachlan & Peel 2000). Here, we assume that the MPCM populationis made up of subpopulations that can each be described by a parametric surface-densitydistribution, and the mixture model is the sum of these subcluster surface densities. Mixturemodels are commonly used in parametric modeling of point processes, though usually withGaussian function (rather than an astrophysically motivated function) components (McLach-lan & Peel 2000; Fraley & Raftery 2002). Our procedure is based on several well-establishedsteps: subcluster properties are obtained through maximum likelihood parameter estimation;the number of subclusters is determined by model selection; and subcluster assignments foryoung stars are determined using posterior probabilities from the fitted models (Everitt et al.2011). For example, these steps performed by the cluster-finding program mclust (Fraley& Raftery 2012), which uses normal mixture models and performs model selection usingthe penalized likelihood Bayesian Information Criterion (BIC). Here, we use a subclustermodel that has a roughly flat core and power-law wings, which is generally a better fit toyoung stellar clusters, and we perform model selection using the Akaike Information Criterion(AIC).Figure 1 shows the projected spatial distribution of stars in the Orion Nebula region, inboth individual stars and smoothed stellar surface density, with the clusters found from thefinite mixture model algorithm indicated by black ellipses. In this field, this cluster-findingmethod identifies layered, overlapping clusters and clusters with different sizes and surfacedensities. The statistical analysis is performed using software from the R package spatstat (Baddeley & Turner 2003). This environment facilitates advanced statistical analysis of
13 –spatial point processes within irregularly shaped windows.
In this study, we adopt a parametric model of young stellar clusters based on equilibriumconfigurations of self-gravitating stars in the form of isothermal ellipsoids (Chandrasekhar1942). 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Σ( r ) = Σ r/r c ) , (1)where Σ( r ) is the surface density of stars at a distance r from the cluster center, which isparameterized by the coordinates of the cluster center, the central surface density Σ , andthe core radius r c . This approximation is accurate to <
7% for r < r c (Binney & Tremaine2008). The central volume density ρ is related to the central surface density by ρ = Σ / r c , (2)and the number of stars N projected within a radius r is N ( < r ) = πr c Σ ln(1 + r /r c ) . (3)However, many young stellar clusters have elliptical shapes on the sky, for example Hil-lenbrand & Hartmann (1998) find that the ONC has elliptical contours of constant stellarsurface density with ellipticity (cid:15) (cid:39) .
30. A more general surface density profile is producedby stretching the Hubble model to model ellipsoidal clusters. This fitting function introducestwo new parameters, the ellipticity (cid:15) = ( a − b ) /a and the ellipse orientation φ . The surfacedensity for the isothermal ellipsoid becomesΣ ell ( r ; Σ , r , r c , φ, (cid:15) ) = Σ (cid:34) (cid:12)(cid:12)(cid:12)(cid:12)(cid:18) (1 − (cid:15) ) − /
00 (1 − (cid:15) ) / (cid:19) ˆ R ( φ )( r − r ) (cid:12)(cid:12)(cid:12)(cid:12) (cid:44) r c (cid:35) − , (4)where the matrix ˆ R ( φ ) rotates the vector r − r = ( x − x , y − y ) by angle φ . A single isother-mal ellipsoid component has six parameters: central right ascension x , central declination y , Σ , r c , (cid:15) , and φ . This distribution function is similar to King’s profile, but it lacks the truncation radius where densitygoes to zero. The King profile’s truncation radius is derived from the assumption of a tidal cutoff, which isunlikely to apply to young stellar clusters. For the datasets in our analysis, it is unlikely that model fittingwould be able to constrain any truncation radius due to the limited field of view, the complex distributionsof subclusters that overlap at large radii, and the presence of an unclustered component.
14 –Fig. 1.— Left: MPCMs in the Orion star-forming region (IR-excess selected stars are redpoints; X-ray selected stars without IR-excess detections are green points.) The clustersfound by our statistical method are shown as black ellipses, indicating the cluster coreregion. Overlapping clusters are resolved, indicating the layered structure of the region.Right: The adaptively smoothed stellar surface density is shown, with the ellipses over-plotted. Subclusters from the literature that are identified in the finite mixture model arelabeled (black arrows), including the BN/KL cluster (Becklin & Neugebauer 1967; Kleinmann& Low 1967) and stars embedded in the southern end of the OMC-2/3 molecular filaments(Megeath et al. 2012); the OMC-1S cluster is not resolved in the model (gray arrow). TheONC is modeled by an inner and outer ellipsoid, which may be interpreted a core-halo clusterstructure. 15 –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 remain-ing background contaminants). Distributed populations of young stars have been found inNGC 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 modelsfor all k subclusters (plus the unclustered component) multiplied by mixture coefficients, a i ,given by, Σ θ ( r ) = k +1 (cid:88) i =1 a i Σ i ( r ; θ i ) = k (cid:88) i =1 a i Σ ell ( r ; x ,i , y ,i , r c,i , φ i , (cid:15) i ) + a k +1 Σ U , (5)where θ = { a , x , , y , , r c, , φ , (cid:15) , . . . , a k , x ,k , y ,k , r c,k , φ k , (cid:15) k , a k +1 } denotes the model pa-rameters. Our finite mixture model is only defined over the field of view, so it is integrableand can be used as a probability density function. We adopt the method of maximum likelihood estimation (MLE) that, following theprinciples 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 equationln L ( θ ; X ) = N (cid:88) i =1 ln Σ θ ( r i ) − (cid:90) W Σ θ ( r (cid:48) )d r (cid:48) , (6)under the assumption of a Poisson point process X = { r i } Ni =1 containing N points (stars),which are spatially distributed with the surface density of the finite mixture model Σ θ in awindow (field of view) W . The integral in the second term is the expected number of starswithin the field of view, and the model Σ θ is normalized so that this term equals the totalnumber of observed stars N .MLE of mixture models for cluster analysis is described by Everitt et al. (2011, Chpt.6). We seek the mode (highest value) of the likelihood function. This can be computationallychallenging; for example, if 5 subclusters are present in a MYStIX field, then the parameterspace to be searched to maximize the likelihood has 6 × R ’s function “optim” is applied to iteratively find the global maximumlikelihood solution. Visual examination of a smoothed map of residuals between the dataand the model shows that the MLE model solutions are quite satisfactory in explaining theMPCM spatial structure using isothermal ellipses.The R code for model fitting is presented in Appendix A, using the NGC 6357 MSFRas an example. As is often the case for complex model fits, multiple fitting attempts areneeded 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 thatbalances optimal fit to the data and some concept of model parsimony. This requires thatthe maximum likelihood for each value of k be “penalized” by some function of k . The mostcommon choices are the Bayesian Information Criterion (BIC) and the Akaike InformationCriterion (AIC): BIC = − L + (6 k + 1) ln N (7) AIC = − L + 2(6 k + 1) , (8)where 6 k + 1 is the number of parameters for k subclusters, and N is the total numberof stars in the MPCM sample. The optimal k value is chosen to minimize the BIC orAIC. When models are found with AIC values very close to the best fit model, the finalmodel is chosen to minimize maximum excursions in the smoothed residual maps. The BICplaces a stronger emphasis of model parsimony (simplicity) than the AIC, whereas the AICsometimes overestimates the number of components. The choice of AIC, BIC or other modelselection 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, sec 6.5.2) discuss the issue in the context ofmixture models for cluster analysis.Experimentation with MYStIX samples showed that the BIC model selection criterion 17 –produces final models with fewer subclusters, missing sparser clusters that are readily seenin visual inspection of the MPCM spatial distribution. For a typical MYStIX field with N (cid:39) × ln 1000 (cid:39)
40 in the likelihood. However, a subclusterthat improves the log likelihood by (cid:39)
12 will be acceptable to the AIC. We therefore havechosen the AIC as our model selection criterion in order to capture a greater dynamic rangeof rich and poor clusters.This model selection has two further advantages. First, it obviates the need to arbitrarilyspecify a critical parameter value in the clustering procedure, such as k in k -nearest neighborsor a special branch length in the minimal spanning tree. Second, it reduces features in theresidual 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 identifyingwhere 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 imple-mented by the “diagnose.ppm” function in spatstat . There will be both positive and negativeresiduals, and the residuals over the entire region will integrate to zero. A better fit willhave residuals with lower amplitude that are randomly distributed throughout the regionwithout any coherent structures in the residuals. The residuals for MYStIX clusters will belarger where the surface density is higher; however, there should be roughly equal numbersof positive and negative residuals that are not strongly correlated with the locations of thesubclusters.To smooth the residual maps, we used kernels with bandwidths that are four times themedian nearest-neighbor distance in a region, which produce residual maps for the differentregions that have similar statistical balances between variance and bias. Kernels that aretoo small produce lumpy structures due to stochastic fluctuations, while kernels that aretoo large blur distinct structures together. The residual maps often show coherent positiveresiduals that cannot be eliminated by adding additional subclusters to the model. Theseare typically curved or irregular structures that cannot be well modeled by the isothermalellipsoids used here.In addition to examination of the smoothed residual map, model goodness of fit maybe determined using a χ test based on the number of counts in arbitrarily chosen spatial 18 –bins. If the field of view is divided into n regions A i with approximately equal number ofcounts, then the expected number of points in the i th regions is E ( N i ) = (cid:82) A i Σ( u )d u for asurface-density model Σ. For a sufficiently large number of counts, we can test goodness offit using Pearson’s statistic X = n (cid:88) i =1 [ N i − E ( N i )] E ( N i ) . (9)For this test, the field of view is tessellated using adaptive polygonal cells that typicallycontain ∼
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 particularsubcluster. This is a “soft classification” where probabilities are calculated based on therelative contribution of the different components of the posterior model to the total surfacedensity at the location of a star. This stands in contrast to “hard classification” methodsthat places every object into a single “best” class.However, approximate lists of stars belonging to each subcluster are very useful forastronomical study, so we generate hard classifications based on the mixture model softclassifier. One possibility for assigning stars to a particular subcluster (or to the unclusteredcomponent) is to assign them to the model component that is most probable; this strategy isused by the mclust algorithm (Fraley & Raftery 2012). However, to avoid including the mostambiguous stars in lists for specific subcusters, we additionally require that the probabilityfor the assigned subcluster be better than 30%. We further require that stars assigned to asubcluster be within an ellipse four times the size of the cluster core. (These decision rulesplay 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 thescientific goals of analysis. The MPCM data have variations in surface density of several The degrees of freedom of the problem may not be well-defined when the data are used for both modelfitting and decisions on binning, affecting the interpretation of the X statistics (Feigelson & Jogesh Babu2012).
19 –orders of magnitude and complex, overlapping structures. Our algorithm resolves thesestructures into individual components, and constructs a model for a star-forming regionbased 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 considerationsled us to choose the parametric cluster-analysis methods. However, a this approach mayencounter a variety of challenges or disadvantages.
A model form must be assumed
It is not clear how well finite mixture models workwhen the component distribution models are poor descriptions of the data. Here weassume that real star clusters have an isothermal ellipsoid distribution function, whichis clearly not always true. However, deviations from the isothermal structure can leadto interesting inferences about particular regions. Overall, the clean appearance of theresidual maps and excellent global χ values for the best-fit models (Table 4) showthat the assumption of isothermal ellipsoids is satisfactory. Guidance on the model selection criteria is lacking
Statisticians have not settled ona universal penalty to maximum likelihood estimation for model selection. The AIChas a strong theoretical basis in information theory and thermodynamics, the BIC hasa strong theoretical basis in Bayesian inference, and both are related to minimum- χ fitting under some circumstances. Finite mixture model analysis using the BIC is lesssensitive at identifying small clusters in fields of view that also contain rich clustersbecause the model selection criterion involves improvements to the global likelihood,and does not measure localized surface density enhancements. The AIC used here ismore effective at finding the small clusters. But, the theory does not exist to interpretchanges in the AIC in terms of statistical significance. Finding the best mixture model is computationally difficult
Fitting takes place ina high-dimensional parameter space, even though the point pattern data are two-dimensional. Additional problems may be encountered with infinite likelihoods, mul-timodality of the likelihood (Kalai et al. 2012), and slow convergence of parameterestimators (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 modelfitting and selection rely on penalized likelihood techniques widely used across statis- 20 –tics. In contrast, nonparametric methods, such as the friends-of-friends and the k -nearest neighbor algorithms, provide no mathematical guidance to estimate the num-ber of clusters. They often exhibit weak performance in the presence of noisy orhigh-dynamic range structures (Everitt et al. 2011). A critical scale or surface density is not assumed
Cluster size and central density arefree parameters for each subcluster in the finite mixture models, so structures at vastlydifferent spatial scales can be fit within a single region. This is in sharp contrast tomost cluster algorithms, which rely on a single spatial threshold (e.g. critical distanceor density) applied throughout the region. Given that it is an open area of researchwhether or not subclusters of young stars form with a uniform surface density, applyinga single spatial threshold for cluster detection risks biasing the result with a biasedassumption. Decision rules that play the role of a critical scale are applied only whenthe stellar membership of a given subcluster needs to be defined.
Distributed populations are treated
The finite-mixture model does not insist that theclusters include every star, but rather allow relatively isolated stars to be assignedto a distributed population with uniform surface density. In this way, our approachis similar to density based clustering algorithms, rather than standard nonparametricclustering algorithms that require that every object lie in a cluster.
Probabilistic cluster membership assignments are used
The finite-mixture model isa soft classification method, rather than a hard classification method. Each pointhas 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 prob-ability thresholds can be used for different applications that require different levels ofreliability. The hard classification methods do not include this information, and theambiguous membership will be definitively assigned into a cluster.
The models are flexible and effective
The parametric model chosen can be modified tosuit the situation being modeled. Other methods may require incorrect assumptions(e.g. normal distributions) or lack any explicit assumptions. An elliptical isothermalmodel permits a wide array of different morphologies. If the model is astrophysicallyreasonable, then the parametric clustering method is better optimized for detectionstructures than a nonparametric method. 21 –
4. Results4.1. Surface Density Maps
Surface density maps of young stars are created using an adaptive density estimator onprojected star positions. The density estimator is based on the Voronoi tessellation (Ogata etal. 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 ofview into cells that tend to be smaller in dense regions and larger in less dense regions. Therest 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 ofsurface 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 factorsof 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 densitiesthat depend on the distances and telescope exposure times for each region. The valuestherefore cannot be directly compared between regions.Surface densities vary up to a factor of ∼
100 within a single regions, and many regionsare multimodal with local maxima differing in surface density by more than an order ofmagnitude. An effective cluster finding algorithm must be able to identify these differentsubclusters despite the high contrast in surface density. 22 – - : Right ascension D ec li n a t i on
30 20 10 5:42:00 50 40 30 20 41:10 - : - : Right ascension D ec li n a t i on - : - : Right ascension D ec li n a t i on
50 45 40 35 8:59:30 25 20 15 10 05 - : Right ascension D ec li n a t i on Fig. 2.— The smoothed projected stellar surface density is shown, with a color bar in unitsof observed stars pc − . Brown contours show increase in surface density by factors of 1.5.Black ellipses mark the core regions of the isothermal ellipsoid subclusters. Left to right andtop to bottom: Orion, Flame, W 40, RCW 36. 23 – : : Right ascension D ec li n a t i on - : Right ascension D ec li n a t i on
30 20 10 7:19:00 50 40 30 20 10 18:00 - : - : Right ascension D ec li n a t i on Fig. 2.— Left to right and top to bottom: Rosette, Lagoon, NGC 2362. 24 – : Right ascension D ec li n a t i on
30 9:00:00 30 8:59:00 58:30 - : Right ascension D ec li n a t i on
30 21:00 30 17:20:00 19:30 - : - : Right ascension D ec li n a t i on
30 26:00 30 17:25:00 30 24:00 - : Right ascension D ec li n a t i on Fig. 2.— Left to right and top to bottom: DR 21, RCW 38, NGC 6334, NGC 6357. 25 – - : Right ascension D ec li n a t i on - : Right ascension D ec li n a t i on
10 18:03:00 50 40 30 20 10 02:00 01:50 - : - : Right ascension D ec li n a t i on
30 20 10 5:23:00 50 40 30 20 22:10 : Right ascension D ec li n a t i on Fig. 2.— Left to right and top to bottom: Eagle, M 17, Trifid, NGC 1893. 26 –
40 30 20 10 6:41:00 50 40 30 20 40:10 : : Right ascension D ec li n a t i on - : - : Right ascension D ec li n a t i on Fig. 2.— Left to right: NGC 2264, Carina. 27 –
Table 2 is a catalog of the 142 subclusters found in the best-fit models of 17 MYStIXMSFRs. For each subcluster (letter designations are assigned from west to east) the tableprovides celestial coordinates of the subcluster center, semi-major and semi-minor axes of thecore ellipse in parsecs, ellipticity, orientation, and total number of stars within an ellipse fourtimes the size of the core — all of which are obtained from model fitting ( § J − H and median X-ray M E values. In low-density subclusters, the number of stars inferred fromthe model, N , obs , may not precisely equal the number of stars assigned to that subcluster ifthere is spatial overlap with neighboring subclusters leading to ambiguities in assignments.One can see some strong successes in the decomposition of stars in to isothermal ellip-soids, but also some limitations. In the Orion region (Figure 1) the subclusters disentanglelayered 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 corre-sponding to a dense “core” and the other corresponding to a less dense “halo.” This resultagrees with Henney & Arthur (1998), who find that Orion proplyd and non-proplyd starsare strongly peaked around θ Ori C with statistically-significant excess above the surfacedensities expected from the King profile fit by Hillenbrand & Hartmann (1998). On theother hand certain structures are not modeled: no distinction is made between the unab-sorbed ONC members and the Orion Molecular Cloud stars (OMC; 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 smallwith few stellar members; for example, subclusters B, E, F, J, K in Rosette and subcluster J 28 –in M 17 have fewer than 5 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 whichis associated with a well-studied embedded ultracompact H II region Kleinmann & Wright(1973). Nevertheless, these subclusters often correspond to clumps of stars that appearreal to visual inspection. On the other hand, occasionally model components contain largenumbers of stars that do not appear to be discrete subclusters. This includes subcluster Gin Carina and subcluster E in NGC 6357.Getman et al. (2014a, 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 sparserstructures are coincident with molecular cloud cores, providing corroborating evidence thatthey 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). 29 –Table 2. MYStIX Subcluster Properties
Subcluster α δ r c, major r c, minor (cid:15) φ N , obs J − H M E (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.6Orion B 05 35 15.7 -05 23 23 0.06 0.04 0.30 28 73 0.87 1.6Orion C 05 35 16.7 -05 22 34 0.31 0.16 0.49 5 834 1.05 1.6Orion D 05 35 17.8 -05 16 35 0.23 0.04 0.84 12 48 1.17 1.4Flame A 05 41 42.5 -01 54 14 0.15 0.10 0.37 146 219 1.79 2.8W 40 A 18 31 26.7 -02 05 39 0.17 0.16 0.04 107 187 2.10 2.5RCW36 A 08 59 27.1 -43 45 22 0.20 0.13 0.36 122 196 1.63 2.3RCW36 B 08 59 27.3 -43 45 26 0.03 0.01 0.61 25 24 2.14 2.8NGC 2264 A 06 40 31.5 +09 49 52 0.09 0.08 0.14 136 11 0.83 1.1NGC 2264 B 06 40 37.1 +09 47 31 0.04 0.02 0.44 124 5 0.60 1.1NGC 2264 C 06 40 40.4 +09 51 06 0.03 0.03 0.11 8 6 0.73 1.1NGC 2264 D 06 40 45.8 +09 49 03 0.13 0.11 0.11 194 9 0.65 1.1NGC 2264 E 06 40 59.1 +09 52 22 0.30 0.16 0.47 83 54 0.63 1.1NGC 2264 F 06 40 59.2 +09 53 59 0.07 0.03 0.54 69 12 0.61 1.0NGC 2264 G 06 40 59.4 +09 36 12 0.10 0.07 0.31 96 31 2.11 2.3NGC 2264 H 06 41 02.1 +09 48 44 0.19 0.16 0.19 19 15 0.65 1.0NGC 2264 I 06 41 04.5 +09 35 57 0.11 0.05 0.55 14 31 1.16 1.6NGC 2264 J 06 41 06.3 +09 34 09 0.16 0.12 0.25 63 62 1.22 1.7NGC 2264 K 06 41 08.2 +09 29 53 0.25 0.11 0.55 0 71 0.75 1.3NGC 2264 L 06 41 12.8 +09 29 10 0.03 0.03 0.13 83 16 1.84 3.5NGC 2264 M 06 41 14.9 +09 26 42 0.09 0.06 0.32 20 20 0.73 1.2Rosette A 06 30 57.1 +04 57 57 1.20 0.87 0.28 90 58 0.79 1.4Rosette B 06 31 20.5 +04 50 14 0.16 0.16 0.00 0 4 0.73 1.4Rosette C 06 31 32.0 +04 50 58 0.29 0.09 0.71 129 7 0.77 1.4Rosette D 06 31 55.4 +04 56 39 0.14 0.04 0.72 87 11 1.17 1.4Rosette E 06 31 59.3 +04 54 50 0.91 0.83 0.08 111 244 0.76 1.3Rosette F 06 32 05.5 +04 48 20 0.20 0.09 0.51 0 4 0.83 1.4Rosette G 06 32 46.5 +04 45 35 0.16 0.16 0.00 0 2 0.76 1.0Rosette H 06 33 07.2 +04 46 57 0.79 0.15 0.82 93 34 0.84 1.3Rosette I 06 33 10.2 +04 31 04 0.25 0.11 0.56 42 5 0.94 1.6Rosette J 06 33 15.1 +04 35 08 0.18 0.09 0.53 156 4 2.18 2.4
30 –Table 2—Continued
Subcluster α δ r c, major r c, minor (cid:15) φ N , obs J − H M E (J2000) (J2000) (pc) (pc) (deg) (stars) (mag) (keV)(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)Rosette K 06 33 20.1 +04 37 01 0.09 0.09 0.00 90 2 2.00 2.6Rosette L 06 34 10.7 +04 25 06 1.16 0.57 0.51 167 158 1.10 1.5Rosette M 06 34 12.9 +04 19 07 0.65 0.25 0.61 71 81 2.21 2.3Rosette N 06 34 31.4 +04 19 09 0.18 0.16 0.08 119 10 1.62 1.9Rosette O 06 34 37.1 +04 13 02 0.22 0.05 0.75 151 9 1.76 2.2Lagoon A 18 03 23.8 -24 15 19 0.35 0.16 0.55 126 30 0.84 1.4Lagoon B 18 03 40.1 -24 22 40 0.07 0.05 0.28 134 34 1.22 1.8Lagoon C 18 03 46.3 -24 22 01 0.25 0.12 0.52 109 27 0.85 1.4Lagoon D 18 03 51.3 -24 21 08 0.10 0.06 0.41 135 9 0.78 1.3Lagoon E 18 04 07.6 -24 25 53 0.68 0.27 0.60 119 62 0.83 1.3Lagoon F 18 04 13.3 -24 18 27 1.14 0.73 0.36 3 243 0.80 1.3Lagoon G 18 04 20.1 -24 22 51 0.09 0.05 0.40 156 17 0.83 1.3Lagoon H 18 04 23.3 -24 21 13 0.26 0.21 0.20 6 81 0.79 1.3Lagoon I 18 04 28.3 -24 22 46 0.34 0.31 0.10 17 126 0.83 1.3Lagoon J 18 04 39.6 -24 23 20 0.25 0.24 0.05 4 47 0.87 1.3Lagoon K 18 04 50.5 -24 26 19 0.46 0.25 0.45 39 86 1.03 1.4NGC 2362 A 07 18 37.8 -24 53 58 0.21 0.13 0.38 54 20 0.65 1.2NGC 2362 B 07 18 42.9 -24 57 44 0.43 0.38 0.10 172 126 0.62 1.09DR 21 A 20 38 51.7 +42 18 52 0.20 0.06 0.69 60 12 1.97 2.3DR 21 B 20 38 57.7 +42 17 51 0.07 0.05 0.34 248 5 2.04 2.8DR 21 C 20 39 00.2 +42 18 47 0.07 0.04 0.41 144 11 3.00 3.5DR 21 D 20 39 00.4 +42 19 46 0.17 0.07 0.61 168 22 2.92 4.0DR 21 E 20 39 00.5 +42 22 39 0.22 0.12 0.44 147 75 2.36 3.7DR 21 F 20 39 00.9 +42 24 43 0.06 0.03 0.55 50 5 2.55 4.0DR 21 G 20 39 03.8 +42 16 51 0.14 0.11 0.22 149 11 2.51 3.3DR 21 H 20 39 03.9 +42 25 34 0.11 0.05 0.53 167 22 2.96 3.3DR 21 I 20 39 05.4 +42 21 18 0.13 0.09 0.32 154 15 2.66 3.0RCW 38 A 08 58 47.2 -47 30 52 2.57 2.07 0.19 177 314 1.36 2.2RCW 38 B 08 59 05.0 -47 30 43 0.12 0.08 0.36 54 117 1.06 2.6RCW 38 C 08 59 16.4 -47 27 50 0.42 0.09 0.78 165 15 1.62 2.5NGC 6334 A 17 19 57.9 -35 54 02 0.15 0.10 0.31 126 17 1.44 2.1
31 –Table 2—Continued
Subcluster α δ r c, major r c, minor (cid:15) φ N , obs J − H M E (J2000) (J2000) (pc) (pc) (deg) (stars) (mag) (keV)(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)NGC 6334 B 17 19 58.5 -35 56 24 0.42 0.31 0.25 66 42 1.44 2.0NGC 6334 C 17 20 02.7 -35 58 23 0.06 0.05 0.15 107 11 1.10 1.8NGC 6334 D 17 20 14.9 -35 54 43 0.08 0.07 0.09 35 10 1.88 2.8NGC 6334 E 17 20 19.0 -35 54 59 0.27 0.21 0.24 28 44 1.88 3.1NGC 6334 F 17 20 23.2 -35 57 00 0.22 0.19 0.17 25 18 1.28 1.8NGC 6334 G 17 20 25.2 -35 44 04 0.13 0.09 0.26 49 21 1.61 2.5NGC 6334 H 17 20 31.2 -35 54 14 0.18 0.13 0.26 77 23 1.38 1.8NGC 6334 I 17 20 35.3 -35 59 20 0.14 0.09 0.32 52 12 1.01 1.6NGC 6334 J 17 20 39.2 -35 49 27 0.48 0.17 0.65 40 75 2.27 3.2NGC 6334 K 17 20 48.4 -35 42 59 0.14 0.09 0.31 75 9 3.00 3.1NGC 6334 L 17 20 54.4 -35 45 43 0.27 0.23 0.16 50 50 2.44 3.2NGC 6334 M 17 20 57.3 -35 39 41 0.20 0.16 0.20 167 18 2.03 2.7NGC 6334 N 17 21 32.5 -35 40 27 0.34 0.24 0.30 28 25 2.20 1.6NGC 6357 A 17 24 43.7 -34 12 07 0.28 0.22 0.22 113 151 1.26 1.9NGC 6357 B 17 24 46.7 -34 15 23 0.58 0.37 0.36 127 133 1.30 1.9NGC 6357 C 17 25 33.3 -34 24 43 0.30 0.24 0.19 145 128 1.29 1.9NGC 6357 D 17 25 34.3 -34 23 10 0.06 0.04 0.38 176 27 1.27 2.0NGC 6357 E 17 25 47.9 -34 27 12 1.06 0.18 0.83 14 60 1.33 1.9NGC 6357 F 17 26 02.2 -34 16 42 0.37 0.19 0.50 4 162 1.41 2.1Eagle A 18 18 39.9 -13 47 41 0.12 0.12 0.01 139 46 0.88 1.4Eagle B 18 18 42.2 -13 47 03 0.90 0.45 0.50 148 451 0.98 1.5Eagle C 18 18 52.8 -13 46 43 0.25 0.10 0.60 103 20 0.96 1.4Eagle D 18 18 57.3 -13 45 23 1.65 0.65 0.61 100 275 1.10 1.6Eagle E 18 19 07.9 -13 36 28 0.08 0.07 0.07 13 20 2.40 2.6Eagle F 18 19 12.8 -13 26 03 0.45 0.40 0.13 123 37 1.61 2.3Eagle G 18 19 13.1 -13 33 45 0.20 0.10 0.51 102 20 2.59 3.7Eagle H 18 19 15.0 -13 39 25 0.11 0.08 0.26 51 10 1.57 1.6Eagle I 18 19 19.2 -13 36 34 0.44 0.16 0.64 37 18 1.64 2.4Eagle J 18 19 29.7 -13 23 08 0.30 0.18 0.41 30 12 1.51 2.1Eagle K 18 19 30.7 -13 45 30 0.06 0.06 0.04 50 6 1.42 2.0Eagle L 18 20 02.8 -13 48 26 0.04 0.02 0.31 143 7 2.01 3.0
32 –Table 2—Continued
Subcluster α δ r c, major r c, minor (cid:15) φ N , obs J − H M E (J2000) (J2000) (pc) (pc) (deg) (stars) (mag) (keV)(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)M 17 A 18 20 18.0 -16 14 11 0.07 0.07 0.03 41 10 1.92 3.2M 17 B 18 20 19.5 -16 13 28 0.03 0.02 0.10 80 5 2.49 2.9M 17 C 18 20 21.6 -16 10 35 0.14 0.13 0.09 135 31 1.70 2.5M 17 D 18 20 22.5 -16 08 23 0.29 0.27 0.07 173 119 1.53 2.4M 17 E 18 20 22.5 -16 09 52 0.09 0.07 0.18 68 13 1.56 2.5M 17 F 18 20 22.8 -16 12 23 0.05 0.04 0.03 15 7 1.70 4.4M 17 G 18 20 25.0 -16 11 35 0.04 0.04 0.05 94 6 1.55 3.9M 17 H 18 20 25.6 -16 11 17 0.11 0.11 0.07 151 31 1.25 2.2M 17 I 18 20 26.0 -16 09 34 0.17 0.17 0.02 24 48 1.45 2.1M 17 J 18 20 27.8 -16 09 54 0.01 0.01 0.02 137 4 1.56 2.0M 17 K 18 20 28.5 -16 10 58 0.19 0.16 0.14 177 78 1.41 2.2M 17 L 18 20 29.9 -16 10 46 0.14 0.13 0.05 34 95 1.82 2.8M 17 M 18 20 30.8 -16 03 16 0.20 0.16 0.21 94 31 1.83 2.5M 17 N 18 20 31.4 -16 09 39 0.15 0.11 0.27 4 32 1.33 2.1M 17 O 18 20 31.5 -16 11 22 0.09 0.07 0.18 57 26 1.42 2.4Carina A 10 43 41.9 -59 35 49 0.19 0.09 0.55 115 11 1.03 1.6Carina B 10 43 54.2 -59 33 02 0.70 0.45 0.35 191 193 0.97 1.5Carina C 10 43 56.4 -59 32 54 0.22 0.18 0.17 80 106 0.94 1.4Carina D 10 44 32.9 -59 33 42 0.60 0.30 0.50 56 46 0.90 1.4Carina E 10 44 34.0 -59 44 08 0.28 0.27 0.05 137 36 0.84 1.4Carina F 10 44 37.4 -59 26 03 0.64 0.40 0.36 108 38 0.94 1.5Carina G 10 44 39.7 -59 44 26 10.15 1.88 0.81 177 741 0.91 1.5Carina H 10 44 41.8 -59 22 05 0.30 0.23 0.23 209 49 0.83 1.3Carina I 10 44 45.3 -59 20 07 0.23 0.18 0.24 70 18 0.80 1.3Carina J 10 45 02.4 -59 45 50 0.50 0.46 0.09 163 65 0.96 1.5Carina K 10 45 06.2 -59 40 21 0.35 0.25 0.28 94 36 0.89 1.5Carina L 10 45 11.1 -59 42 46 0.49 0.36 0.27 229 58 0.91 1.5Carina M 10 45 13.7 -59 57 58 0.48 0.33 0.32 122 52 0.95 1.6Carina N 10 45 36.1 -59 48 20 0.39 0.38 0.03 84 37 1.51 2.1Carina O 10 45 53.4 -59 56 53 0.18 0.17 0.09 110 25 1.15 1.8Carina P 10 45 54.4 -60 04 32 1.17 0.59 0.50 60 118 0.87 1.5
33 –Table 2—Continued
Subcluster α δ r c, major r c, minor (cid:15) φ N , obs J − H M E (J2000) (J2000) (pc) (pc) (deg) (stars) (mag) (keV)(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)Carina Q 10 45 55.2 -59 59 51 0.51 0.33 0.35 92 27 0.87 1.4Carina R 10 46 05.4 -59 50 09 1.06 0.51 0.52 83 46 0.99 1.6Carina S 10 46 52.7 -60 04 40 0.44 0.26 0.41 24 24 0.90 1.5Carina T 10 47 12.5 -60 05 58 0.53 0.46 0.14 215 48 0.92 1.5Trifid A 18 02 05.8 -23 05 53 1.88 0.36 0.81 103 87 1.53 1.3Trifid B 18 02 23.1 -23 01 50 0.14 0.11 0.20 158 25 0.92 1.4Trifid C 18 02 31.0 -23 00 40 0.65 0.52 0.20 92 82 0.87 1.3Trifid D 18 02 39.7 -22 49 18 0.87 0.28 0.68 64 27 1.27 1.4NGC 1893 A 05 22 40.3 +33 22 15 0.43 0.32 0.26 98 46 0.81 1.5NGC 1893 B 05 22 46.6 +33 25 17 0.36 0.29 0.21 22 110 0.81 1.5NGC 1893 C 05 22 47.4 +33 24 11 0.23 0.20 0.14 44 20 0.85 1.5NGC 1893 D 05 22 49.7 +33 24 42 0.03 0.03 0.09 13 5 0.75 1.4NGC 1893 E 05 22 49.7 +33 26 51 0.15 0.05 0.67 129 9 0.80 1.4NGC 1893 F 05 22 53.6 +33 25 40 0.17 0.14 0.16 58 18 0.81 1.4NGC 1893 G 05 22 53.8 +33 30 54 0.28 0.19 0.33 25 46 0.89 1.6NGC 1893 H 05 22 57.5 +33 26 34 0.31 0.28 0.10 28 60 0.82 1.4NGC 1893 I 05 23 03.3 +33 28 15 0.48 0.32 0.33 25 114 0.83 1.5NGC 1893 J 05 23 08.5 +33 28 32 0.04 0.03 0.37 82 11 0.97 1.7Note. — The ellipsoid components from the selected finite-mixture model. Column 1: The componentname, labeled from east to west. Columns 2-3: Celestial coordinates (J2000) for the ellipsoid center.Column 4-5: Semi-major 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 starsestimated 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 measureof absorption.
34 –Table 3. Membership of MYStIX Subclusters
Designation R.A. (J2000) Dec. (J2000) X-ray IR-excess OB ME J − H log L X, tc MSFR Subclusterdeg deg (keV) (mag) (erg s − )(1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11)182025.46-160943.4 275.1060920 -16.1620560 1 0 0 2.10 1.38 30.78 M 17 I182025.49-161213.8 275.1062390 -16.2038430 1 0 1 3.70 -9.99 30.66 M 17 U182025.49-161053.8 275.1062440 -16.1816170 0 0 0 2.28 1.31 31.65 M 17 H182025.52-161016.0 275.1063340 -16.1711260 1 0 0 1.98 0.96 30.44 M 17 X182025.52-160906.8 275.1063380 -16.1518940 1 0 0 1.55 1.15 29.79 M 17 I182025.52-160558.5 275.1063600 -16.0996060 1 0 0 3.70 -9.99 31.04 M 17 U182025.53-161004.6 275.1063810 -16.1679530 0 0 0 1.76 1.11 30.79 M 17 I182025.53-161157.7 275.1063840 -16.1993820 1 0 0 3.31 3.18 31.08 M 17 X182025.55-160931.2 275.1064880 -16.1586940 0 0 0 1.93 1.24 30.66 M 17 I182025.56-161124.7 275.1065330 -16.1901970 0 0 0 3.26 -9.99 30.48 M 17 HNote. — The source properties and subcluster assignments for all MPCMs from Broos et al. (2013) within the spatial analysis fields of view.This table thus contains stars that are not in the “flux flattened” sample used to define the subclusters. Column 1: MYStIX MPCM designation(Broos et al. 2013). Columns 2–3: MPCM celestial coordinates. Columns 4–6: Binary variables indicating inclusion (1) or exclusion (0) in theflux flattened X-ray sample, IR-excess selection, or OB selection. Column 7: X-ray median energy in the 0.5–8.0 keV band, when sufficient countsare available. Column 8: J − H color index. Column 9: Inferred absorption corrected X-ray luminosities in the 0.5–8.0 keV band, when available.Column 10: The name of the MYStIX MSFR. Column 11: Subcluster assignment; “A”–“T” indicates the assigned subcluster, “U” indicates theunclustered stellar population, and “X” indicates uncertain assignment. This table is available in its entirety in the electronic edition of the journal.A portion is shown here for guidance regarding its form and content.
35 –Fig. 3.— The MPCM stars in the Carina Nebula field of view are shown as colored pointsindicating subcluster assignments: 20 colors for 20 subclusters, light blue points for the non-clustered stellar population, and gray for the unassigned stars. Black ellipses represent the 20subclusters. All MPCMs are shown (not only the statistical sample defined in Section 2.4),so the “egg-create” effect is evident. Diagrams of the other 16 regions are included in anonline-only figure set. 36 –Orion, Flame, W 40, RCW 36Rosette, Lagoon, NGC 2362, DR 21RCW 38, NGC 6334, NGC 6357, EagleM 17, Trifid, NGC 1893, NGC 2264Fig. 3.— electronic figure set 37 –
The properties of the best-fit model for each region are provided in Table 4, includingthe number of subclusters, the observed surface density of the unclustered young stars, theminimum and maximum of the fitting error residuals and the kernel size used to calculatethese residuals, and a χ goodness-of-fit statistic.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 ofstar formation throughout the Galaxy makes it difficult to choose a sample field of view thatadequately characterizes a region. Often clusters are located near each other and it is unclearwhether they should be treated separately or as subclusters, as is the case with the three richclusters in NGC 6357 and the Tr 14-15-16 clusters in the Greater Carina Nebula. Selectioneffects in our survey of young stars also influence the number of inferred subclusters.The goodness-of-fit statistics for the tessellated observed vs. model MPCM counts inTable 4 suggest that, in most cases, the finite mixture models are neither over-fitted ( χ (cid:28) d . o . f . ) nor poorly fit ( χ (cid:29) d . o . f . ). Poor fits are an indication that structures are not simpleisothermal ellipsoids, which is not surprising in elaborate, complex structures seen in regionslike 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 andnegative residuals in units of stars pc − , respectively. For Carina the maximum residual islocated south-west of Tr 14, between Subclusters A and B – this indicates an asymmetryto Tr 14 due to an excess of stars at the edge of the Northern Cloud in this star-formationcomplex. Otherwise, the positive and negative residuals are low and more-or-less evenlydistributed.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 inselect regions are described here:
Orion Nebula
The minimum residual is located south of the ONC cluster, and may indicatesome deviation from a true isothermal ellipsoid shape. There is also a large negativepatch 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 north-west of the cluster, whichcorresponds to a group of highly absorbed stars. 38 –
W 40
The residuals are nearly zero in the cluster core, but there is a positive residual to thesouth-east, composed mostly of IR-excess selected stars, with a population too smallto be accepted as a subcluster by the AIC statistic.
Rosette Nebula
The residuals are higher amplitude than in other MYStIX regions, andmay 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 longpositive residual stretching from the west to the south of NGC 6530. This residualcoincides with the edge of the molecular cloud and a bright rimmed cloud; it mayrepresent an area of triggered star formation.
DR 21
The large blue residual along the chain of subclusters in DR 21 indicates problemswith this fit. This is also implied by the poor χ statistic for this region. These resultssuggest that the deeply embedded clusters in DR 21 are distributed along the molecularfilament 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 northof 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 complexof stars (G353.1+0.6; Townsley et al. 2014). The best AIC is obtained when threeellipsoids were used to fit this complex, although three distinct components are notobvious 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 sub-clusters (A and B). There are large positive residuals at the edges of the cluster; never-theless, 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 nonparametricadaptive smoothing (x-axis) is plotted against the surface density obtained from the para-metric finite mixture mode (y-axis) where each point corresponds to a pixel in the residualmap. These generally agree to within 0.1 dex. The local maxima from the surface densitymaps are typically somewhat higher in the models than in the smoothed maps; this is anartificial effect of the smoothing algorithm. For Carina, the agreement is good at surfacedensities above ∼ − . Deviations at lower surface densities may reveal differencesbetween the modeling of unclustered young stars with a constant surface density and thetrue surface-density variations of the unclustered stars. 39 –Table 4. Best-Fit Models Residual Map Goodness of FitRegion N subclus Σ U kernel min max χ d.o.f. P χ (stars pc − ) (pc) (stars pc − ) (stars pc − )(1) (2) (3) (4) (5) (6) (7) (8) (9)Orion 4 30 0.076 290 (4%) -270 (4%) 47 60 0.23Flame 1 4 0.12 80 (4%) -60 (3%) 21 23 0.78W 40 1 20 0.17 38 (3%) -26 (2%) 24 17 0.22RCW 36 3 7 0.11 96 (3%) -57 (2%) 17 17 0.82NGC 2264 13 4 0.28 16 (3%) -20 (4%) 50 46 0.62Rosette 15 0.7 0.72 3.8 (7%) -2.5 (5%) 126 62 < < < < χ value, degrees of freedom, and the probability, respectively.
40 –Fig. 4.— Left: Residual surface density for the Carina Nebula. Negative residuals are shownin blue and positive residuals are shown in red as indicated by the color bar in units ofstars pc − . Right: Log surface density obtained from the finite mixture model is plottedagainst log surface density from the adaptively smoothed surface density maps. The y = x lines is shown in red. Analogous plots for the other sixteen MYStIX regions are provided bythe electronic version of this article. 41 –Fig. 4.— Electronic Figure Set (left to right and top to bottom: Orion, Flame, W 40,RCW 36, NGC 2264, Rosette, Lagoon, NGC 2362, DR 21, RCW 38, NGC 6334, NGC 6357,Eagle, M 17, Carina, Trifid, NGC 1893) 42 –
5. Four Classes of Star-Forming Complexes
The 17 MSFRs examined here have a wide diversity of structures, but they generallycan be sorted into four classes. Figure 5 compares the cluster structure (indicated by thepatterns of core ellipses) for 15 of the regions on the same parsec scale. The Carina andRosette Nebulae complexes are omitted because they can be viewed as composites of thefour classes, as described below.
Simple Isolated Cluster Structure
Some regions are dominated by a simple, unimodalcluster that is well fit by a single ellipsoid model. Regions with this structure are themost likely to be dynamically relaxed. This class includes the Flame Nebula, W 40,the principal clusters in the Trifid Nebula and NGC 2362, the two northern subclustersA and F in NGC 6357, and the main cluster NGC 2244 in the Rosette Nebula. TheFlame Nebula is slightly elongated along the direction of the filamentary molecularcloud that passes through the region, while W 40 is round. NGC 2362 is composed oftwo overlapping subclusters, but both are well fit by ellipsoid models. Although thestructure appears dynamically relaxed may suggest that these regions tend to be older,regions in this class are not always the oldest MYStIX regions. Some have molecularmaterial and show signs of active star formation. The two-body relaxation time formany regions is probably longer than the age of the cluster.
Core-Halo Structure
Some clusters appear unimodal with little substructure around amain cluster, but still be poorly fit by a single isothermal ellipsoid model because ofan excess of stars near the center of the cluster. The finite-mixture-model algorithmoften breaks these up into two ellipsoids, one used to model the cluster core and theother to model the cluster halo. Examples include the ONC, RCW 36, RCW 38,Eagle Nebula subclusters A-B, and Tr 14 in the Carina Nebula. These results confirmcore-halo structures reported in the ONC and Tr 14 by Henney & Arthur (1998) andAscenso et al. (2007), respectively. In the cases of the Orion Nebula and RCW 36, thefitted core radius of the central cusp is extraordinarily small, ∼ . − .
02 pc. Thecase of RCW 38 appears different because the core subcluster is larger and off-center.Although, the simple isolated clusters and core-halo clusters may be surrounded byperipheral subclusters or unclustered stars, more than half of individual MPCM starsfrom the flattened sample lie within 5 core radii of the principal clusters.
Clumpy Structure
In some of the MYStIX regions, it is difficult to model the surfacedensity distribution into a few discrete isothermal ellipsoids because of clumpy struc-tures within the main cluster. Regions with such clumpy structure include M 17, theLagoon Nebula, and the Eagle Nebula. These three cases all have central dominant 43 –subclusters, but a substantial fraction of the stars lie in secondary subclusters thatappear contiguous with the central concentration. The embedded subclusters in theSouth Pillars of the Carina Nebula or in the Rosette molecular cloud may also beconsidered clumpy, but these cases lie far from the dominant clusters.
Linear chains of clusters
This structure consists of subclusters that are roughly arrangedin a line. Five MYStIX clusters exhibit such an arrangement: DR 21, NGC 2264,NGC 1893, NGC 6334, and the Carina Nebula. The length of clusters in the MYStIXsample ranges from ∼
10 pc to 30 pc, and in some cases can be as narrow as ∼
6. Properties of MYStIX Subclusters6.1. Subcluster Sizes
The size of subcluster cores varies from 0.01 pc to > p > . ∼ . ∼ M E and J − H for the cluster. A negative relation canbe seen between these quantities, which is emphasized by the gray LOWESS non-parametricregression. LOWESS is a well-established regression technique that combines robust leastsquares 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 medianenergy range corresponds to A V < A V (cid:39)
50 mag (Getman et al. 2010). Althoughthere is much scatter in these relations, nonparametric correlation hypothesis tests using theSpearman rank coefficient (and confirmed using bootstrap resampling) give significant levels p < .
001 for the radius–
M E relation and p < .
01 for the radius– J − H relation.Some of the highly absorbed subclusters with small radius are well-known embedded starforming 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 48 – log Core Radius (pc) N u m be r −2.0 −1.5 −1.0 −0.5 0.0 0.5 Fig. 7.— Histograms of subcluster core radii, for all subclusters (white) and embeddedsubclusters (gray) with a lognormal fit. 49 –median
M E > . A V >
17 mag) are rarely larger than 0.2 pc. Thisphenomenon is more evident in the histogram of core radii, which shows that clusters withmedian X-ray median energies greater than 2.5 (shaded gray region in Figure 8) typicallyhave smaller radii than less embedded clusters. The peak of their lognormal distribution isat 0.08 pc, and the Anderson Darling test shows no statistically significant deviation from alognormal distribution ( p > . J − H index for absorption can serve as a proxyfor age — this is directly demonstrated by Getman et al. (2014a). Thus, the negative size–absorption relation can be viewed as a positive size–age relationship suggesting that clustersexpand as they age. Gieles et al. (2012) find that the combination of multiple clusters atdifferent stages of expansion can produce a distribution of stellar surface densities similar towhat is seen in the MYStIX regions.This effect may explain some of the differences in cluster sizes between different starformation regions. For example, DR 21 and NGC 1893 are both composed of linear chainsof subclusters, but the subclusters in NGC 1893 are larger. This makes sense in light of thesize–absorption relation because the younger DR 21 subclusters are deeply embedded in amolecular filament, while the older NGC 1893 subclusters are lying in a bubble and are lightlyabsorbed. Marks & Kroupa (2012) identify a typical half-mass radius for young embeddedclusters of less than a few tenths of a parsec; for MYStIX subclusters with M E > .
04 pc (Figure 8).We note that the typical size of star clusters is similar to the width of filaments foundby
Herschel which appear to be ubiquitous in MSFRs and involved in the cluster formationprocess (Andr´e et al. 2010). The variation in size and absorption of subclusters within asingle star forming region can be interpreted as an indication of non-coeval clustered starformation. 50 – − . − . − . . J−H (mag) l og c o r e r ad i u s ( p c ) − . − . − . . X−ray median energy (keV) l og c o r e r ad i u s ( p c ) Fig. 8.— Dependence of cluster size on absorption measures. Left: Subcluster core radiusvs. median J − H mag. Right: Subcluster core radius vs. median X-ray median energy inthe 0.5 − Figure 9 (left) shows a histogram of MYStIX cluster ellipticity. The measured elliptici-ties will be affected by projection, so they are lower limits on intrinsic ellipticity. Ellipticitiesrange from 0 to ∼ (cid:15) (cid:39) . Chandra field of view.The ellipticity distribution for MYStIX subclusters has some resemblance to the distri-bution of ellipticities in simulated star-forming regions found by Maschberger et al. (2010).There, ellipticites range from 0.0 to 0 .
82, lower ellipticites are more common, and the peakellipticity is 0 .
33. In the simulated clusters, significant ellipticities are an effect of subclustermergers, and this may also be the case for the MYStIX clusters. The frequent presenceof high ellipticities in MYStIX subclusters suggests that many have not reached dynamicalequilibrium. 52 –
Ellipticity N u m be r −2.0 −1.5 −1.0 −0.5 0.0 0.5 . . . . . . Core radius (pc) E lli p t i c i t y Fig. 9.— Left: Histogram of all 142 subclusters’ ellipticities is shown. The ellipticities rangefrom ∼
0, or nearly circular in projection, to 0.9, a 10 to 1 aspect ratio, with low ellipticitysubclusters favored. Right: The scatter plot shows subcluster core size vs. core radius. 53 –
7. Discussion and Conclusions
The MPCM samples have been designed to improve the statistical characterization ofstellar populations in MSFRs (Feigelson et al. 2013). Using a rich statistical sample of 16,608MPCMs, “flattened” to remove X-ray spatial sensitivity variations, and the finite-mixture-model technique to identify subclusters among these young stars, new details about thespatial structure of young stellar clusters emerge. In the diverse MSFRs surveyed by theMYStIX project, we see a picture of star-forming regions with multiple non-coeval subclustersthat probably formed with similar small sizes, but have expanded to a wide range of sizes.A similar finding of cluster expansion in a different sample is studied by Pfalzner (2009).The catalog of 142 MYStIX subclusters derived here, combined with inferred intrinsic stellarpopulations (Paper II) and age estimates for the subclusters (Getman et al. 2014a), will helpprovide powerful links between observational and theoretical analysis star-cluster formation.The finite mixture model, using isothermal ellipsoids, is a novel cluster analysis methodfor MSFRs that can reveal information that is hard to obtain from other cluster-findingmethods. Pfalzner (2011) finds that multimodal properties in regions with varying clustersurface density and richness are difficult to determine from surface density distributionsalone. But our method is an objective and mathematically founded procedure to separatemultiply-clustered patterns into distinct isothermal ellipsoidal subclusters. The procedure isflexible, finding statistically significant structures on all scales simultaneously, even treatingcases such as small clusters reside inside large ones. This is important because, even inrelatively simple MSFRs like the Orion Nebula, layered subcluster structures are presentand can affect inferred subcluster properties. Furthermore, our algorithm does not requirea single threshold in surface-density or separation to divide subclusters, which can helpavoid bias when subclusters with very different properties exist in the same star-formingregion. However, the finite-mixture-model method requires additional assumptions aboutregion structure and subcluster shape.We validate the subcluster results using residual maps, χ tests of source counts, andcomparison with results from previous studies of these regions. The models reproduce the ob-served surface densities reasonably well, although residual maps reveal some extra structure.The improved stellar census in dense regions improves the accuracy of subcluster radius andpeak surface density maps. On larger scales, we caution that the MPCM sample coverageof widely distributed young stars is limited, and isothermal ellipsoid fitting is not sensitiveto any truncation of cluster sizes that may be present. 54 – The numbers of subclusters and their arrangements vary significantly from region toregion. Heuristically, the clustering morphologies are divided into four classes: isolatedclusters, core-halo clusters, clumpy clusters, and linear chains of clusters (Section 5). Thedominant subcluster in a region (i.e. the subcluster with most stars, N , obs ) often sits at thecenter 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 nosubcluster that is particularly dominant. While, in contrast, regions with simple or core-halostructure may be much more centrally concentrated.In DR 21, NGC 2264, and NGC 6334, the linear chain cluster structure maps ontomolecular filaments in the regions, and it is clear that the star cluster structure was at leastpartially inherited from the filamentary chain of molecular cloud clumps. More detail aboutcorrespondence between subcluster locations and molecular cloud clump locations will begiven in a later MYStIX study. Thus, regions with linear chain structure have not evolvedmuch from their initial state at the time of the birth of stars. The more centrally concentratedclumpy, core-halo, and simple cluster structures are more likely to have dynamically evolvedand, in cases where the subclusters are well fit by single isothermal ellipsoids, may haveachieved some dynamical relaxation.The morphological classification suggests an evolutionary progression from linear chainstructure to clumpy structure to core-halo structure to simple cluster structure. Hierarchicalcluster formation, often in molecular filaments, would evolve into dynamical relaxed uni-modal structures. However, the MYStIX regions indicate that this progression is not alwaysfollowed. 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. Alognormal distribution of subcluster size peaks at 0.17 pc, and many subclusters have el-lipticities of (cid:15) ≤ .
3, but few have ellipticities greater than 0.5. No relationships betweenellipticity and size are seen.An additional property of subclusters, their absorption, can help link these propertiesto cluster formation history and evolution. Line-of-sight absorption is measured both by thenear-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 representslocal cloud cores within which the youngest clusters reside. In complex MYStIX star-formingregions, a wide range of absorption can be seen for the different subclusters, indicating thesesubclusters are not coeval. A statistically significant anti-correlation between absorption and 55 –size of a subcluster is found, most easily explained if both absorption and size are relatedto age. This issue is pursued in the accompanying study by (Getman et al. 2014a) wherespatial-age gradients in MYStIX regions are found based on a new age indicator for individualpre-main sequence stars. A similar relationship was found for subclusters within the RosetteMolecular Cloud by (Ybarra et al. 2013) who infer a gas removal e -folding timescale of0.4 Myr.Cluster expansion is expected to produce a wide range of stellar surface densities seenin many star-forming regions (Gieles et al. 2012). Early cluster expansion arises principallyfrom the loss of the gravitational potential of the molecular material, but stellar dynamicscan also produce stellar halos around dense cores through three-body interactions with hardbinaries and gravothermal core collapse associated with mass segregation (Giersz & Heggie1996). This issue will be discussed further in Paper II when intrinsic populations, ratherthan observed population subject to different sensitivity effects, are derived. The anti-correlation between radius and absorption provides evidence that subclusterexpansion is important part of subcluster dynamical evolution. From numerical simulations,Moeckel & Bate (2010) find that even the very dense simulated cluster produced by Bate(2009a), with a half-mass radius of ∼ > ∼ ∼ ∼ ∼ ∼ ∼ < ∼
30% using theclusters in Lada & Lada (2003). 56 –The ellipticities of MYStIX subclusters could be inherited from the structure of theparental molecular cloud or could arise from subcluster mergers. The subcluster ellip-ticity distribution resembles those found by numerical simulations of merging subclusters(Maschberger et al. 2010, their Figure 10b) in which merger products have a distribution ofellipticities, including a tail at high ellipticities.The cluster morphological classes may correspond to structures seen in star formationsimulations, such as those performed by Bate et al. (2003); Bate (2009a) or Parker & Meyer(2012). Simulations often show multiple pockets of star formation occurring along gas fil-aments as the cloud collapses, which can end up merging and dynamically relaxing as thesimulation progresses. Clumpy clusters may indicate early stages of partial merging, core-halo structures may be an intermediate stages of mergers, and isolated clusters with relaxedstructure a final stage. Getman et al. (2014b) add an important constraint to such a sce-nario: in at least two MYStIX regions, the dense core of the cluster has younger stars thanthe outer regions. This observational result requires some special behavior, such as continualfeeding of gas towards the cluster center or dispersal of older core stars into the halo.The simple clusters may represent those that have completed subcluster merging anddynamical 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 fractionof volume containing subclusters. When this filling factor becomes high, tidal interactionsaccelerate subcluster mergers and formation of a simple unimodal cluster is promoted. Thisphenomenon 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 subclustermergers. Bate (2009a) report a hydrodynamical and N-body simulation of a MSFR resultingin a cluster with a half-mass radius of 0.05 pc and a halo of stars extending out ∼ Chandra
ACIS Team contract SV4-74018(G. Garmire & L. Townsley, Principal Investigators), issued by the
Chandra
X-ray Center,which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASAunder contract NAS8-03060. M.A.K. also received support from NSF SI2-SSE grant AST-1047586 (G. J. Babu, PI). We thank Leisa Townsley for providing advice on this paper andfor the development of
Chandra tools. This research made use of data products from the
Chandra
Data Archive. This work is based on observations made with the
Spitzer SpaceTelescope , obtained from the NASA/IPAC Infrared Science Archive, both of which are oper-ated by the Jet Propulsion Laboratory, California Institute of Technology under a contractwith the National Aeronautics and Space Administration. This research has also made use ofSAOImage DS9 software developed by Smithsonian Astrophysical Observatory and NASA’sAstrophysics Data System Bibliographic Services.
A. Finite Mixture Model with R Code
The procedure for fitting finite mixture models is illustrated here for the MYStIX regionNGC 6357 with the code from the comprehensive R public domain statistical software system(R Core Team 2013). The R code below makes use of online source code and machinereadable tables. Two packages from the Comprehensive R Analysis Network (CRAN) areused: plotrix (Lemon 2006) and spatstat (Baddeley et al. 2005).The structure of the example is as follows: Step 1 loads the required R functions anddata; Step 2 creates a first guess of the ellipsoid parameters; Step 3 conducts the first iterationof model fitting using the Nelder-Mead simplex algorithm; Step 4 displays the results of thefirst iteration; Step 5 refines the model with further fitting iterations; and Step 6 presentsthe results of the fit, producing figures analogous to Figures 2 and 4a in the paper. Thesecommands may be pasted into an R session following the instructions below.The input data consist of the polygonal boundaries of the fields of view and coordi-nates of the stars in the statistical sample in machine readable table format ( fov.mrt and stars.mrt ), which are available in the online version of this article in addition to the libraryof R functions sourcecode.R . The library commands below load the necessary libraries. The These libraries can be installed from the R session by install.packages("spatstat",dependencies=T)
58 –function star.ppp creates a point pattern object (an internal format of the spatstat package)for the region NGC 6357 assuming a distance of 1.7 kpc. Finally, the function unique removesany duplicated positions in the catalog.
To illustrate the analysis, we fit a model with k = 6 subclusters and the uniform den-sity component. The data structure “param2” specifies: the number of model components(6 + 1 = 7), the model that will be used for each component (ellipsoid or constant), andthe number of additional parameters for each component (5 or 0 for ellipsoid or constant,respectively). The function ell.model is an implementation of Equation 4 for the isother-mal ellipsoid in Section 3.1. The array “param.init” contains the initial guess parametersfor all components of the finite mixture model — we can inspect these parameters moreeasily by running the param2ellipse command, which returns a data-frame object with k rows (subclusters) and 6 columns: “x” – ellipsoid x coordinate; “y” – ellipsoid y coordinate;“core” – core radius in parsecs; “theta” – ellipse orientation in radians; “b” – axis ratio;and “mix” – the log of the mixture coefficient corresponding to the relative peak surfacedensities for each subcluster (the log mixture coefficient for the first component is assumedto be 0.0 and not included in the “param” array.) The final entry in “param” is the logmixture coefficient for the unclustered component. and install.packages("plotrix",dependencies=T) .
59 – -2.50, -4.50, 0.3, 0, 1, -0.7415462,-1.25, -2.50, 0.3, 0, 1, 0.8510107,-4.00, 0.75, 0.3, 0, 1, 0.1332489,4.00, 4.00, 0.3, 0, 1, 0.1667037,-2.0);param2ellipse(param.init);
The function multi.model is the finite mixture model, which is specified by the variables“k,” “param,” and “param2.” The function model.lik calculates log-likelihood of this modelusing Equation 6. Model fitting is done using the Nelder-Mead simplex algorithm via optim .This R function takes the star cluster data, the model form and parameters, and the log-likelihood function. It is unlikely that the absolute maximum of the likelihood will be foundon the first attempt, so this step must be performed iteratively — with certain parametersheld fixed in turn while other are free to vary — until the likelihood converges on themaximum. The commands below plot the spatial point pattern object ( plot ) and draws ellipses( draw.model ) indicating the initial guess parameters (red) and new parameters (green). Thelog-likelihood “L” and the AIC are calculated — these values will improve with each iteration.
It is necessary to iteratively converge on the maximum likelihood. The mask.freeze func-tion allows only certain parameters to be fit, while the other parameters are held constant.This can be used to concentrate on improving certain parameters without having to search 60 –the entire (6 k +1)-dimensional parameter space. The commands below improve, sequentially,subcluster positions, subcluster sizes, subcluster orientations, subcluster shapes, subclusterpositions again, and, finally, all parameters again. After each iteration, the choice of thenext fitting function to run depends on inspection of the results of the previous iteration.Nevertheless, each iteration will always result in an AIC value that is greater or equal to theprevious iteration’s AIC values, so the order of these commands is not particularly importantas long as there are enough iterations for the AIC to converge. (AIC is calculated again atthe end of this fitting sequence.) The time for these steps to run depends on the numberof points and the dimensionality of the parameter space being searched, so fits with frozenparameters will run more quickly than fits where all parameters are free to vary. Now that the model has been fit, we can visualize the results. The adaptive.density command produces the adaptively smoothed surface density maps described in Section 4.1.(Here, we only run 10 iterations rather than 100 used to produce Figure 2.) The function make.fig2 produces an image of subcluster ellipses overlaid on the non-parametric surfacedensity map similar to Figure 2. The function make.fig2 produces an image of the spatial fitresiduals similar to Figure 4.
61 – make.fig4(param, model=multi.model, clust=clust,param2=param2, bandwidth=0.38);
A similar procedure is applied for different numbers of ellipsoids (with the initial guessparameters “n”, “param”, and “param2” adjusted accordingly) to find the number of sub-clusters that optimizes the AIC value. It is useful to try the fitting procedure starting withseveral different initial guesses in order to ensure that the results are independent of theinitial parameters used.
REFERENCES
Allison, R. J., Goodwin, S. P., Parker, R. J., et al. 2009, ApJ, 700, L99Allison, R. J., Goodwin, S. P., Parker, R. J., et al. 2009, MNRAS, 395, 1449Alves, J., & Bouy, H. 2012, A&A, 547, A97Andr´e, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102Ascenso, J., Alves, J., Vicente, S., & Lago, M. T. V. T. 2007, A&A, 476, 199Baddeley, A. J., Turner, R. 2003, Journal of Statistical Software, 12, 1Baddeley, A. J., Turner, R., Moller, J. and Hazelton, M. 2005, Journal of the Royal StatisticalSociety, Series B 67, 617Baddeley, A. J., Moller, J. and Pakes, A.G. 2008, Annals of the Institute of StatisticalMathematics 60, 627Baddeley, A. J., Turner, R., Mateu, J., & Bevan, A. 2013, Journal of Statistical Software,55, 1Baddeley, A. 2007, in
Statistical Challenges in Modern Astronomy IV (eds. G. J. Babu & E.D. Feigelson), Astro. Soc. Pacific Conf. Ser. 371, 22Banerjee, S., & Kroupa, P. 2013, ApJ, 764, 29Bate, M. R., Bonnell, I. A., & Bromm, V. 2003, MNRAS, 339, 577Bate, M. R. 2009, MNRAS, 392, 590 62 –Bate, M. R. 2009, MNRAS, 397, 232Battersby, C., Bally, J., Ginsburg, A., et al. 2011, A&A, 535, A128Becklin, E. E., & Neugebauer, G. 1967, ApJ, 147, 799Benjamin, R. A., Churchwell, E., Babler, B. L., et al. 2003, PASP, 115, 953Binney, J., & Tremaine, S. 2008,
Galactic Dynamics , 2nd ed., Princeton University PressBonnell, I. A., Smith, R. J., Clark, P. C., & Bate, M. R. 2011, MNRAS, 410, 2339Bontemps, S., Andr´e, P., K¨onyves, V., et al. 2010, A&A, 518, L85Bragg, A. E., & Kenyon, S. J. 2005, AJ, 130, 134Bressert, E., Bastian, N., Gutermuth, R., et al. 2010, MNRAS, 409, L54Brooks, S., Gelman, A., Jones, G. & Meng, X.-L. (eds.) 2011,
Handbook of Markov ChainMonte Carlo , Chapman & Hall/CRCBroos, P. S., Townsley, L. K., Feigelson, E. D., et al. 2011, ApJS, 194, 2Broos, P. S., Getman, K. V., Povich, M. S. et al., 2013, ApJS, in pressBurkert, A., & Hartmann, L. 2004, ApJ, 616, 288Burnham, K. P. & Anderson, D. (2002)
Model Selection and Multi-Modal Inference , 2nd ed.,SpringerCartwright, A., & Whitworth, A. P. 2004, MNRAS, 348, 589Chandar, R., Whitmore, B. C., Calzetti, D., et al. 2011, ApJ, 727, 88Chandrasekhar, S. 1942,
Principles of Stellar Dynamics , University of Chicago PressChen, J. 1995, Ann. Statist., 23, 221Clarke, C. J., Bonnell, I. A., & Hillenbrand, L. A. 2000, in
Protostars and Planets IV (V.Mannings, eds.), 151Cleveland, W. S. 1979, Journal of the American Statistical Association 74, 368, 829Cohen, D. H., Kuhn, M. A., Gagn´e, M., Jensen, E. L. N., & Miller, N. A. 2008, MNRAS,386, 1855 63 –Elmegreen, B. G., & Elmegreen, D. M. 2001, AJ, 121, 1507Elmegreen, B. G. 2000, ApJ, 530, 277Evans, N. J., II, Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321Everitt, B. S., Landau, S., Leese, M. & Stahl, D. 2001,
Cluster Analysis , 5th ed., ArnoldFall, S. M., Chandar, R., & Whitmore, B. C. 2009, ApJ, 704, 453Feigelson, E. D., Getman, K., Townsley, L., et al. 2005, ApJS, 160, 379Feigelson, E. D. 2010, Proc. NAS, 107, 7153Feigelson, E. D., Getman, K. V., Townsley, L. K., et al. 2011, ApJS, 194, 9Feigelson, E. D., & Jogesh Babu, G. 2012, Modern Statistical Methods for Astronomy, byEric D. Feigelson , G. Jogesh Babu, Cambridge, UK: Cambridge University Press,2012,Feigelson, E. D., Martin, A. L., McNeill, C. J., Broos, P. S., & Garmire, G. P. 2009, AJ,138, 227Feigelson, E. D., Townsley, L. K., Broos, P. S., et al. ApJS, in pressFisher, R. A. 1922,
Phil. Trans. Royal Society, A , 222, 309Fraley, C. & Raftery, A. E. 2002,
J. Amer. Stat. Assn. , 97, 611Fraley, C., Raftery, A. E. Murphy, T. B., Scrucca, L. 2012,F˝ur´esz, G., Hartmann, L. W., Megeath, S. T., Szentgyorgyi, A. H., & Hamden, E. T. 2008,ApJ, 676, 1109Getman, K. V., Flaccomio, E., Broos, P. S., et al. 2005, ApJS, 160, 319Getman, K. V., Feigelson, E. D., Luhman, K. L., et al. 2009, ApJ, 699, 1454Getman, K. V., Feigelson, E. D., Broos, P. S., et al. 2010, ApJ, 708, 1760Getman, K. V., Broos, P. S., Feigelson, E. D., et al. 2011, ApJS, 194, 3Getman, K. V., Feigelson, E. D., Sicilia-Aguilar, A., et al. 2012, MNRAS, 426, 2917Getman, K. V., Kuhn, M. A., Feigelson, E. D., et al. 2013, ApJ, in press, arXiv:1403.2741 64 –Getman, K. V., Feigelson, E. D., Kuhn, M. A. 2013, ApJ, in press, arXiv:1403.2742Gieles, M., Moeckel, N., & Clarke, C. J. 2012, MNRAS, 426, L11Giersz, M., & Heggie, D. C. 1996, MNRAS, 279, 1037Goodwin, S. P., & Bastian, N. 2006, MNRAS, 373, 752Gutermuth, R. A., Megeath, S. T., Myers, P. C., et al. 2009, ApJS, 184, 18G¨udel, M., & Naz´e, Y. 2009, A&A Rev., 17, 309G¨udel, M., Briggs, K. R., Arzner, K., et al. 2007, A&A, 468, 353Harfst, S., Portegies Zwart, S., & Stolte, A. 2010, MNRAS, 409, 628Hartmann, L., Megeath, S. T., Allen, L., et al. 2005, ApJ, 629, 881Henney, W. J., & Arthur, S. J. 1998, AJ, 116, 322Hillenbrand, L. A., & Hartmann, L. W. 1998, ApJ, 492, 540Kalai, A., Moitra, A. & Valiant, G. 2012, Commun. ACM, 55, 113Kass, R. E. & Wasserman, L. 1995,
J. Amer. Stat. Assoc. , 90, 773Kauffmann, J., & Pillai, T. 2010, ApJ, 723, L7King, I. 1962, AJ, 67, 471King, R. R., Naylor, T., Broos, P. S., Getman, K. V. & Feigelson, E. D. 2013, ApJS, in pressKleinmann, D. E., & Low, F. J. 1967, ApJ, 149, L1Kleinmann, D. E., & Wright, E. L. 1973, ApJ, 185, L131Konishi, S. & Kitagawa, G. (2008)
Information Criteria and Statistical Modeling , SpringerKuhn, M. A., Getman, K. V., Feigelson, E. D., et al. 2010, ApJ, 725, 2485Kuhn, M. A., Getman, K. V., Broos, P. S., Townsley, L. K. & Feigelson, E. D. 2013a, ApJS,in pressKuhn, M. A., Povich, M. S., Luhman, K. L., Getman, K. V., Busk, H. A. & Feigelson, E. D.2013b, ApJS, in pressLada, C. J., & Adams, F. C. 1992, ApJ, 393, 278 65 –Lada, C. J., Muench, A. A., Lada, E. A., & Alves, J. F. 2004, AJ, 128, 1254Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57Lahiri, P., ed. 2001,
Model Selection , Inst. Math. Statistics Lecture Notes
The EM Algorithm and Extensions , 2nd ed., WileyMcLachlan, G. & Peel, D. 2000,
Finite Mixture Models , Wiley-InterscienceMcMillan, S. L. W., Vesperini, E., & Portegies Zwart, S. F. 2007, ApJ, 655, L45Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2012, AJ, 144,
Model Selection (P. Lahiri, Ed.), Institute of MathematicalStatistics Lecture Notes