Star formation across the w3 complex

STAR FORMATION ACROSS THE W3 COMPLEX Carlos G. Román-Zúñiga 1 , Jason E. Ybarra 1 , Guillermo D. Megías 2 , Mauricio Tapia 1 , Elizabeth A. Lada 3 , and Joáo F. Alves 4 1 Instituto de Astronomía, Universidad Nacional Autónoma de México, Unidad Académica en Ensenada, Km 103 Carr. Tijuana–Ensenada, Ensenada 22860, Mexico 2 Facultad de Física. Universidad de Sevilla. Dpto. Física Atómica, Molecular y Nuclear, Sevilla, E-41080, Spain 3 Astronomy Department, University of Florida, 211 Bryant Space Sciences Center, FL 32611, USA 4 Institute of Astronomy, University of Vienna, Türkenschanzstr. 17, A-1180 Vienna, Austria Received 2014 November 24; accepted 2015 June 29; published 2015 August 20 ABSTRACT We present a multi-wavelength analysis of the history of star formation in the W3 complex. Using deep, nearinfrared ground-based images combined with images obtained with Spitzer and Chandra observatories, we identified and classified young embedded sources. We identified the principal clusters in the complex and determined their structure and extension. We constructed extinction-limited samples for five principal clusters and constructed K-band luminosity functions that we compare with those of artificial clusters with varying ages. This analysis provided mean ages and possible age spreads for the clusters. We found that IC 1795, the centermost cluster of the complex, still hosts a large fraction of young sources with circumstellar disks. This indicates that star formation was active in IC 1795 as recently as 2 Myr ago, simultaneous to the star-forming activity in the flanking embedded clusters, W3-Main and W3(OH). A comparison with carbon monoxide emission maps indicates strong velocity gradients in the gas clumps hosting W3-Main and W3(OH)and shows small receding clumps of gas at IC 1795, suggestive of rapid gas removal (faster than the T Tauri timescale)in the cluster-forming regions. We discuss one possible scenario for the progression of cluster formation in the W3 complex. We propose that early processes of gas collapse in the main structure of the complex could have defined the progression of cluster formation across the complex with relatively small age differences from one group to another. However, triggering effects could act as catalysts for enhanced efficiency of formation at a local level, in agreement with previous studies. Key words: infrared: stars –radio lines: ISM –stars: formation –stars: luminosity function, mass function –stars: pre-main sequence –X-rays: stars 1. INTRODUCTION In the current picture of star formation, giant molecular clouds (GMCs)are highly inefficient factories, in which only a small fraction of the available gas is converted into stars. Moreover, star-forming regions in GMC complexes are highly heterogeneous. There is considerable diversity among the stellar aggregations they produce. Differences are evident even in their basic properties, e.g., sizes, numbers, and density structures. Aggregations of young stars in molecular clouds are classified as embedded star clusters. These clusters are born as bound systems of gas and stars. Once gas is removed, embedded clusters may or may not stay bound. They either survive (e.g., open clusters), or dissolve to become part of the field population. For a majority of clusters in the galaxy, the latter is the more common outcome (referred to as the “infant mortality effect”; Lada & Lada 2003). However, recent discussions consider that cluster morphology is more heterogeneous (Bressert et al. 2010; de Grijs 2011). Numerical studies suggest that small groups of stars may be abundant (Adams & Myers 2001)and they may either merge to form clusters or disrupt rapidly into the field. According to various studies (Kruijissen et al. 2012; Banerjee & Kroupa 2015), star clusters may assemble, under certain conditions, from subclusters that merge into larger entities after the gas from which they formed is both disrupted by stellar feedback and torn apart by tidal shocks from the surrounding cloud (the so called “cruel cradle”effect). Unfortunately, the dynamical picture of cluster evolution provided by numerical studies is difficult to directly compare with observations. On one hand, the initial conditions of cluster formation, a key requirement for realistic simulations, are still under debate. On the other hand, young cluster evolution from observations would require additional information on the kinematics of stars (e.g., radial velocities)and spectroscopic age estimations of individual sources to infer the progression of formation. Such studies are only beginning to arise from surveys like APOGEE 5 (e.g., Cottaar et al. 2014)and are limited to nearby ( d 1<kpc)regions. Meanwhile, it is possible to partially reconstruct the history of star formation in a region from photometric information, which can provide evolutive classification and the spatial distribution of young sources, and from molecular gas emission maps, which can provide gas distribution and kinematics. Surveys of GMCs show that embedded star clusters are rarely (if ever)born in isolation. Most embedded clusters form as part of “families”related to a particular complex, defining together a history of star formation from various levels of interaction. For instance, massive stars of one cluster may have influence on the efficiency of star formation of a neighboring cluster in the same cloud. Considering all this, it should be clear that star formation is highly dynamic; the embedded populations we observe are snapshots of a rather convoluted process of evolution and interaction that changes significantly along the star-forming history of a complex. In order to investigate the progression of cluster formation under the influence of the local environment, we considered W3, a prominent cluster-forming region in the “Heart and The Astronomical Journal, 150:80 (18pp), 2015 September doi:10.1088/0004-6256/150/3/80 © 2015. The American Astronomical Society. All rights reserved. 5 Sloan Digital Sky Survey, Apache Point Galactic Evolution Experiment. 1 Soul”molecular complex. W3 is considered a clear example of sequential formation, where cluster formation was induced by the expansion of the giant H II region W4 (Lada et al. 1978; Thronson et al. 1985). One of our goals is to trace the formation of distinct embedded clusters in this region and to attempt to reconstruct the star-forming history of the cloud. The W3 Complex is located at a distance of 2.04 ±0.107 kpc 6 and has hosted at least three major episodes of recent cluster formation. According to spectroscopic studies, IC 1795 formed first, about 3–5 Myr (Oey et al. 2005)ago, followed by the W3 “main” cluster located to its west edge, and the W3(OH)cluster group to the east, both with ages of 2–3 Myr (Navarete et al. 2011; Bik et al. 2012). It has been suggested that IC 1795 triggered the other two episodes in a hierarchical progression (Oey et al. 2005). The Chandra study by Feigelson & Townsley (2008)suggested that the clusters in W3 extend widely and are highly structured, with sources located at relatively large distances from the dynamical centers, including a relatively isolated O star that might have escaped from the main cluster. In two recent studies, Rivera-Ingraham et al. (2011,2013) made use of mid-infrared photometry and far-IR emission mapping from the Spitzer and Herschel 7 space observatories that comprise the entire W3 region. They were able to catalog hundreds of young stellar sources across the complex and determine their spatial distribution. They compiled important evidence that projected distances among YSOs are consistent with cluster-forming clump scales, favoring the cluster-forming mode. They also concluded that small aggregations and distributed populations account for a significant fraction of the recent stellar production in the region. They suggest relatively large age spreads from central to external regions in W3 and proposed a “convergent constructive feedback” scenario, where the gas flows from massive star formation clumps and collects into new dense regions, favoring a progression of formation and the age spread. In this study we emphasize the relative importance of different young stellar groupings in and around IC 1795, W3- Main, and W3(OH), augmenting the level of detail achieved in previous studies. We present a new set of deep, high-resolution near-infrared photometry that we combine with other available data sets, allowing us to increase the number of young star candidates in the complex. We use a K-band luminosity function (KLF)analysis to investigate the presence of an age spread and to help reconstruct the history of star formation in the cloud. Finally, we investigate gas kinematics near the clusters, providing additional information on the gas–star interaction during the early evolution of the complex. 2. OBSERVATIONS AND DATA REDUCTION 2.1. Near-infrared Imaging We obtained near-infrared images of four fields in W3 covering a large area ( 0 .58 0.46´  ). Images were obtained with the Omega 2000 camera at the 3.5 m telescope in the Calar Alto Observatory of the Centro Astronómico Hispano Alemán (CAHA)in Almería, Spain. Omega 2000 provides a 15 15 ¢´ ¢ field of view. Observations were made in J,H, and K(1.1, 1.6, and 2.2 μm, respectively). A list of all fields observed is given in Table 3, which gives the field identification, the center of field positions, observation date, filter, seeing (estimated from the sigma-clipped average FWHM of the stars in each field), and the peak values for the brightness distribution, which is a good estimate of the sensitivity limits achieved. The brightness distribution peaks are in all cases at or above J20.5=mag, H 19.75=mag, and K19.25=mag. 2.1.1. Image Reduction The Omega 2000 images were reduced with modified versions of the FLAMINGOS near-infrared reduction and photometry/astrometry pipelines, which are built in the standard IRAF Command Language environment. One pipeline (see Román-Zúñiga 2006)processes all raw frames by subtracting darks and dividing by flat fields, improving signal- to-noise ratios by means of a two-pass sky subtraction method, and combining reduced frames with an optimized centroid offset calculation. We used dark frames and dome flats obtained within 48 hr of each observation. The final combined product images were then analyzed with a second pipeline (see Levine 2006), which identifies all possible sources from a given field using the SExtractor algorithm (Bertin & Arnouts 1996). We improved the SExtractor detection efficiency in the near-infrared images using a Gaussian convolution filter and maximum deblending (see Bertin & Arnouts 1996). The pipeline then performs Daophot pointspread function (PSF)photometry (Stetson 1987), calibrates observed magnitudes to a zero point, and finds accurate astrometric solutions. After astrometric solutions were found, individual field images were combined into mosaics using Montage. 8 In Figure 1we show an RGB panorama constructed from a combination of the J,H, and Kmosaics. The image represents a complete spatial coverage of the IC 1795, W3-Main, and W3- OH clusters. In the near-infrared the region is transparent to most of the prominent regions of nebulosity and obscuration observable in optical images (e.g., Ogura & Ishida 1976). The near-infrared images at W3(OH)are particularly interesting because they reveal, with unprecedented detail, several small, embedded stellar groups that lie in a small “chain”structure north and east of the W3(OH)cluster (these groups have been previously identified; e.g., Feigelson & Townsley 2008; Navarete et al. 2011). In Figure 2we show a close-up image of this region, where we now combine Hand Kimages with the Spitzer/IRAC 3.6 μm image to enhance illuminated nebulosity features and highlight the most reddened sources. We have labeled the W3(OH)cluster with the letter “A.”The two other most conspicuous groups in the “chain,”both associated with B-type stars (Navarete et al. 2011), have been labeled with letters “B”and “C.”The former group does not have as prominent an overdensity, but as can be seen in the mosaic, it is associated with thick nebulosity, so it is possibly more deeply embedded. There is one more bright source sitting at the center of a cavity between groups B and C, which we suspect could be an additional sub-group. 6 Based on Very Long Baseline Interferometry parallaxes to H 2 O masers in W3(OH)by Hachisuka et al. (2006). 7 Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. 8 http://montage.ipac.caltech.edu 2 The Astronomical Journal, 150:80 (18pp), 2015 September Román-Zúñiga et al. Figure 1. False color mosaic made from J(blue),H(green)and K(red)images from our OMEGA 2000 observations. Figure 2. RGB color mosaic of the W3(OH)region, from Spitzer/IRAC 3.6 μm(red)and Kand HOMEGA 2000 images (green and blue, respectively)showing possible sub-cluster structure in W3(OH). The W3(OH)cluster is labeled as “A,”while two possible sub-clusters are labeled “B”and “C.”These last two labels are placed closer to B-type sources. 3 The Astronomical Journal, 150:80 (18pp), 2015 September Román-Zúñiga et al. 2.1.2. Catalog Preparation The OMEGA 2000 photometry was calibrated relative to the 2MASS using catalogs retrieved from the All Sky Release Point Source databases. Final photometry catalogs were prepared with aid of TOPCAT-STIL (Taylor 2005).We combined individual photometry catalogs from each frame and filter and merged catalogs from all four fields into a master photometry list. In the overlapping areas, we selected the duplicates with the smaller total photometric error across the three bands in either field. Finally, we added 2MASS entries to replace failed measurements from saturated sources. Our final near-infrared catalog contains a total of 72,168 sources. 2.2. Mid-infrared Spitzer Imaging We performed a new reduction and photometric extraction of the W3 mosaics from project 1–127 available in the Spitzer Heritage Archive. Mosaics were constructed by processing both the short and long exposure (0.4 and 10.4 s)Basic Calibrated (BCD)sets from IRAC (all four bands)and the BCD set from MIPS (24 μm only)observations using the MOPEX package. 9 We used the IRAC short exposure data sets to improve extraction of sources down the centers of W3-Main and W3(OH)regions, which present large areas of saturation (increasing with wavelength)due to bright nebulosities in the long exposure sets. However, even at the short exposure images, the bright nebulosity made it very difficult to obtain uniform quality from aperture photometry in the IRAC images, so we performed PSF photometry with Daophot on detection lists obtained with SExtractor. In this case, we improved the detection efficiency using a “Mexican hat”convolution filter and reduced deblending (see Bertin & Arnouts 1996; compare to the near-infrared in Section 2.1.1). Then, our pipeline selects a group of moderately bright and relatively isolated sources in order to make a PSF model for each Daophot run. In all cases we were able to fit the PSF model to a large majority of the detected sources. Final source detection lists were slightly cleaned with the aid of PhotVis (Gutermuth et al. 2004)to remove spurious detections. The resultant photometry lists showed a net reduction of up to 12% in the photometric scatter to the aperture photometry obtained with the APEX software. For the MIPS 24 μm mosaic we used the APEX aperture photometry pipeline, which was good enough for non-saturated regions of the field. In heavily saturated regions, particularly near the centers of W3-Main and W3(OH), missing pixels impeded the ability to make good measurements, and we could not improve these measurements with PSF photometry. Zero points for our Spitzer photometry catalogs were checked by direct comparison of our photometry lists with values from the Galactic Legacy Infrared Midplane Survey Extraordinaire (GLIMPSE)360 Legacy Project catalog 10 and from the photometry tables of YSO candidates in Rivera-Ingraham et al. (2011). The agreement between both GLIMPSE and Rivera-Ingraham et al. (2011)is excellent and zero point corrections were not required. The final photometry catalogs from the IRAC and MIPS 24 μm mosaics were then combined into a mid-infrared catalog with 37267 sources, which we later merged with the near-infrared catalog from our CAHA data set to form a master source catalog. 2.3. Other Data Sets In this study we also make use of other data sets available in the literature and public web-based archives. 1. From the Chandra Source Catalog (CSC; Evans et al. 2010)the positions and X-ray photometry of a total of 611 sources in a box of 3 0 ¢ around the center of IC 1795. 11 2. The 12 CO(2–1)and 13 CO(2–1)maps of the W3 region from the study of Bieging & Peters (2011). 3. The 2.5 level (science grade)SPIRE and PACS mosaics of the W3 region from the Herschel Space Telescope data archive. 4. The Bolocam 1.1 mm mosaic of W3 from the the Canadian Galactic Plane Survey (Aguirre et al. 2011). 3. DATA ANALYSIS AND RESULTS 3.1. Identification of YSOs Using our master infrared photometry catalog we identified Young Stellar Objects (YSOs)in the W3 complex. First, we identified Classes I/0 and II sources from IRAC colors using the criteria applied by Ybarra et al. (2013), which are in turn based on the IRAC color criteria of Gutermuth et al. (2008). Then, using the Chandra CSC catalogs we identified Class III sources as X-ray sources with infrared counterparts and no excess associated with circumstellar material. Finally, we identified additional Classes I and II sources as X-ray CSC sources with infrared excess by applying criteria that combine near-infrared and IRAC colors (see the appendix of Gutermuth et al. 2008). Figure 3shows a H K-versus K4.5[ ] -color–color diagram for all infrared sources that coincide in position with X-ray CSC sources. The colors of the symbols are indicative of the hardness ratio, calculated from the hard and soft X-ray fluxes, Hand S, estimated from aperture photometry measurements, as HSHS ( )( )-+ . It is clear that most of the sources have large hardness ratios, which indicates that in most cases X-ray sources are obscured by large amounts of dust, which reduces the soft X-ray emission. We identified a total of 1008 YSOs in W3, distributed as follows: 60 Class I/0 candidate sources, 780 Class II candidate sources, and 167 Class III candidate sources. Our list contains a significantly larger number of YSO candidates compared to the list of Rivera-Ingraham et al. (2011). We suspect that the differences are mostly due to the addition of X-ray and nearinfrared criteria, which increased the number of candidates in each class. Rivera-Ingraham et al. only list candidates with detections in all IRAC bands, which complicates a direct assessment of completeness. In Figure 4we show the Kand [3.6]band brightness distributions of the three identified YSO classes in our sample. In the case of the [3.6]band, the three distributions show a sharp drop at a similar limit of about 15.5 mag. In contrast, for the Kband, the distribution for Class III sources drops at about 15.0 mag, while the distributions for Classes I and II sources 9 http://irsa.ipac.caltech.edu/data/SPITZER/docs/dataanalysistools/tools/ mopex/ 10 Only the IRAC [3.6]and [4.5]bands are available as this is a warm mission extension of the original GLIMPSE project. 11 Observation ID numbers for W3 ACIS-I data are 446, 611, 7356, 5889, 5890, 5891, 6335, and 6348. 4 The Astronomical Journal, 150:80 (18pp), 2015 September Román-Zúñiga et al. drop near 18.0 mag. These histograms show that due to the high and non-uniform extinction, the brightness distributions for young sources are likely incomplete above the sensitivity limits. However, our near-infrared observations are deep enough to detect young sources across the entire region. For this reason, analysis of the luminosity functions discussed in Section 3.3 were performed using extinction-limited samples. 3.2. Spatial Distribution of Young Stellar Sources In Figure 5we illustrate the spatial distribution of YSO candidates in the W3 region, overlaid on a dust extinction map constructed with the optimized near-infrared excess method NICEST of Lombardi (2009). The YSO sources appear to be very well constrained to the three well known areas of the complex (namely IC 1795, W3-Main, and W3(OH)). Class I sources are mostly confined to areas of high extinction: following the prescription of Lada et al. (2013), we found that over 80% of the Class I sources are located in high extinction regions ( A 7.0 V>mag), and the surface density of Class I sources is linearly correlated with A V within A 3 1 2 V << mag. This is in very good agreement with other cluster-forming regions like Orion A (Lada et al. 2013)and the Rosette Molecular Cloud (Ybarra et al. 2013). Class II sources clearly permeate the entire complex, while Class III sources appear to have a slightly more centrally condensed distribution in the IC 1795 region but are also present in the western edge of W3-Main and surrounding the areas of low extinction to the east of W3(OH). Two aspects are important to note. First, cluster regions are not separated clearly in the layout of the complex, which suggests at least some partial overlap along the line of sight perpendicular to the field. Second, the spatial distribution of young sources suggests that clusters are not single, centrally condensed systems. Instead, they appear to have substructure. In the following we (a)make a consistent separation of the individual cluster populations and (b)provide evidence of the presence of substructure in the cluster population using a surface density analysis. We make the following two definitions. 1. We define a Principal Cluster as a surface density structure that can be associated with a single Gaussian peak with an FWHM of at least 1.0 pc (see Section 3.2.1). 2. We define a substructure as a significant surface density peak, smaller than a principal cluster. We identified as significant all those peaks with a surface density above 1.0 stars per square arcmin (see Section 3.2.2). 3.2.1. Identification of Principal Clusters In order to delimit individual clusters, we employed a Gaussian Mixture Model (GMM)analysis via the mclust package in the Rcode (Banfield & Raftery 1993; Fraley & Raftery 2002; Raftery et al. 2012). The package provides a model-based clustering analysis of any population, based on covariance parametrization and an Expectation-minimiza- tion (EM)algorithm. It can also divide a population in an optimized number of clusters by means of the Bayesian Information Criterion (BIC). The BIC is the value of the maximized log-likelihood with a penalty on the number of model parameters, and allows comparison of models with different parameter sets as well as solutions with different numbers of clusters. We chose the Class II sample as it is the most complete and it is distributed more or less uniformly across the complex. R. A., decl. positions were fed into the code as a sample of bivariate observations. We chose to use circular Gaussians to minimize the number of parameters required for the estimation. Given the complicated layout of the clusters in W3, we opted to run the GMM analysis including noise modeling and rejection Figure 3. H K-vs. K 4.5[]-color–color diagram for CSC X-ray sources coinciding with an source in our infrared master catalog. Open diamond symbols are sources without an infrared excess indicative of a disk. Open square symbols indicate sources classified as Class I. Open triangle symbols indicate sources classified as Class II. The colors of the symbols are coded to the hardness-to-soft X-ray ratio, as indicated in the color bar. Figure 4. Above: K-band brightness distribution for young stellar source candidates identified in the W3 region. Below: [3.6]brightness distribution. 5 The Astronomical Journal, 150:80 (18pp), 2015 September Román-Zúñiga et al. of outliers. 12 For our data set, a model of five clusters provided the best fit for the sample. This is shown in Figure 6. While this finite mixture model shows less detail than—for instance—the nearest neighbor analysis, we see that it succeeds at identifying the three well known clusters in the region (W3-Main, W3 (OH), IC 1795). It also adds two more principal clusters, one coinciding with NGC 896, the compact H II region located at the western edge of the complex, and one located closer to the eastern edge. We labeled this group IC 1795-N. In Table 2we list the result of the GMM analysis, with the centers and radii of the five principal clusters identified (columns 1–4). We emphasize that the purpose of this exercise is to define the principal clusters as consistent samples to study the history of star formation in the complex. The GMM analysis provides such consistent identification, which may differ with respect to other identification methods (e.g., Rivera-Ingraham et al. 2011). 3.2.2. Cluster Substructure in W3. Nearest Neighbors Estimating the surface density of YSOs is often done by applying a variant of the nearest neighbor method (formalized for star cluster detection by Casertano & Hut 1985). This method has been applied successfully in different studies (e.g., Román-Zúñiga et al. 2008; Gutermuth et al. 2010). However, two very important issues regarding this method must be pointed out. First, when calculating distances to the nth neighbor in a distribution of points on a plane, there is a Figure 5. Plot showing the positions of infrared young stellar sources in the entire W3 region. Red symbols indicate sources classified as Class I. Green symbols indicate Class II. Blue symbols indicate Class III sources. The background image and contours correspond to a visual extinction (A V , in units of magnitudes)map of the region constructed with the NICEST method. The map is Nyquist sampled over an 18″pixel grid, smoothed with a Gaussian filter with a 36″FWHM. The contours indicate extinction at A V=7.0, 10.0, 15.0, 20.0, and 25.0 in units of magnitudes. Figure 6. Identification of principal clusters in the W3 region using a Gaussian Model Mixture (GMM). This five cluster model, applied to the projected distribution of the population of Class II YSO candidates in the region, provides the optimal mixture based on a Bayesian Information Criterion (BIC). The circles are placed at the cluster centers determined by the GMM analysis and have radii equal to the Gaussian σ. 12 In our case, we chose to use a first order estimation of noise in the spatial distribution of the sources on the plane based on nearest-neighbors. This way the mclust package applies the hierarchical clustering algorithm to a denoised data set and finally rejects from the final solution those sources that do not belong to any cluster in the final solution (Raftery et al. 2012). 6 The Astronomical Journal, 150:80 (18pp), 2015 September Román-Zúñiga et al. minimum value n=6 to avoid biasing the detection toward small structures formed purely by chance in a two-dimensional distribution of points. Second, when mapping the surface density defined by such distances, the choice of an optimal smoothing kernel that would not dilute genuine substructures or generate spurious ones is not straightforward. These two issues were taken into consideration by Gladwin et al. (1999), who developed two nearest neighbor mapping methods, gather and scatter, both optimized to correctly detect small groups (in our case, possible cluster substructures)in source position maps. The gather algorithm is more adequate for determining hierarchical structure, while the scatter method is more adequate for defining the size of the smallest group in a distribution of points on a map. The gather algorithm was better suited for the task of detecting and separating substructures. We applied gather to the distribution of all Class II YSO candidates in the W3 region, which is the most complete of our samples. Densities were calculated with a square kernel of 36″over a Nyquist sampled R.A.–decl. grid. In Figure 7we show the resultant map, in the form of surface density contours overlaid on an infrared emission image that combines data from the Spitzer MIPS 24μm, HerschelPACS 70 μm, and Herschel SPIRE 250 μm. If we compare the density map to the spatial distribution of sources in Figure 5we see how gather highlights the three main regions (IC 1795, W3-Main, and W3(OH)) as the principal clusters, but also shows a number of peaks that clearly define significant substructure in the spatial distribution of young sources. The positions and equivalent sizes of substructure peaks in the gather map are listed in Table 1. Some of the substructure peaks in the map coincide in position with one or several Class I sources. Those peaks are likely related to the most recent episodes of formation in the complex. Note how none of the peaks near the center, corresponding to IC 1795, are associated with Class I sources, possibly indicating an absence of current formation. We also show the positions of O- and B-type stars listed by Ogura & Ishida (1976), Navarete et al. (2011), and Bik et al. (2012). The largest concentration of massive stars (6 O-type stars, 7 B-type stars)is located at W3-Main, along with a noticeable concentration of protostars, associated with at least four YSO surface density substructures. In IC 1795 there are two O-type stars and one B-type star identified. The central O star in IC 1795 (No. 89 in the Ogura & Ishida 1976 list)is associated with the highest surface density peak. In the W3(OH)region, two substructures define the positions of groups B and C in W3 (OH), and both are associated with B-type stars. It is interesting that the main cluster (group A)in W3(OH)is not associated with any identified massive source. In Figure 7we also placed labels (“RI”)next to regions coinciding with “sub-clusters”identified by Rivera-Ingraham et al. It is worth noting how the gather map highlights the groups A, B, C that were identified by eye in Figure 2. 3.2.3. YSO Ratio Maps It is possible to trace the evolution of a star-forming region through the construction of maps of YSO number ratios. Figure 8shows maps for the Class II to III number ratio, R II:III, and the Class I/0 to II number ratio, R I:I I , both constructed using the method described in Ybarra et al. (2013). For the R II:III map we limited the source candidates to a dereddened brightness H=13.5 mag. This brightness cut avoids biasing the number counts from incompleteness due to patchy extinction. We used the Bayesian estimator described in Ybarra et al. (2013)to calculate the ratios within circular projected regions with radius 3¢ . The map only shows ratios for regions with sufficient number counts such that N NNN10.25 1 1 12 1 2 ()()++< + - , where N 1is the number count of the numerator, and N 2is the number count of the denominator. We found that the whole complex has R 1.0 II:III , which suggests star formation was ubiquitous throughout the complex ∼2–3 Myr ago. We can also trace the most recent star formation using the R I:I I ratio (Figure 8, right). Regions with significant R I:I I are expected to have been actively forming stars within the last Myr. We found that the regions of high R I:I I appear to form an oblong ring around IC 1795. This mirrors the distribution of Class I/0 sources, which also appear to surround IC 1795 (see Figure 7). Interestingly, the high R II:III in IC 1795 suggests that stars were forming there up to 1 Myr ago. Table 1 Substructures in W3 No. Center Coords. Class I a Notes J2000 1 02:27:33.29 61:56:24.2 LL 2 02:27:21.43 61:54:57.4 LW3(OH)-C 3 02:27:16.37 61:54:16.0 LW3(OH)-B 4 02:27:12.43 61:53:08.8 1 L 5 02:27:03.46 61:52:33.3 3 W3(OH)-A 6 02:27:18.20 62:00:04.1 1 R106 7 02:27:08.58 61:56:30.6 1 L 8 02:26:58.74 61:54:53.8 4 L 9 02:26:49.52 61:57:46.0 1 L 10 02:26:48.17 62:03:28.3 LL 11 02:26:40.22 62:00:18.4 LIC 1795 center 12 02:26:41.14 62:04:33.6 LR107 13 02:26:38.87 62:02:01.3 LL 14 02:26:38.87 62:00:42.2 LL 15 02:26:35.50 61:59:15.1 LL 16 02:26:33.25 61:56:25.0 LL 17 02:26:20.01 62:05:01.3 1 L 18 02:26:16.96 61:59:48.7 LL 19 02:26:19.22 61:56:42.8 LR101 20 02:26:03.36 62:07:13.7 LL 21 02:25:59.28 61:58:37.3 LL 22 02:25:57.97 62:08:50.6 1 L 23 02:25:52.42 62:04:21.4 2 R100 24 02:25:51.10 62:00:20.0 LL 25 02:25:47.76 62:10:35.2 LR104 26 02:25:42.74 62:08:10.7 LL 27 02:25:43.36 62:06:18.0 LW3-Main center 28 02:25:32.91 62:06:49.4 LL 29 02:25:32.98 62:04:32.9 2 L 30 02:25:30.88 62:00:05.7 2 NGC 896 center 31 02:25:30.41 61:57:11.6 3 R103 32 02:25:21.90 62:07:06.8 LL 33 02:25:10.23 62:03:11.0 LR105 34 02:25:13.21 61:59:21.6 LL 35 02:25:00.93 62:09:08.7 3 R102 36 02:25:03.58 62:06:52.4 1 L 37 02:25:00.26 62:05:37.0 3 L Note. a Number of Class I candidates within the gather surface density peak contour. 7 The Astronomical Journal, 150:80 (18pp), 2015 September Román-Zúñiga et al. Therefore, the absence of Class I/0 sources in IC 1795 suggests that star formation was recently and possibly abruptly terminated in that region. It is important to note that values for the R II:III ratio in the map can be high (typically larger than 5), particularly because the Class III sources are more spatially scattered. However, R II:III is not the same as the circumstellar disk fraction in the cluster because we are only measuring against weak disk sources and not the whole population (see also Section 3.3). 3.3. KLF of Clusters in the W3 Complex We obtained the KLF for each of the five principal clusters identified with our GMM analysis using the following procedure. Each cluster sample is composed of all K-band sources detected within a “3σ”radius defined as R r3.0 8 ln 2 3cl ()= s, where r cl is defined from the GMM analysis as the half-width at half-maximum of the best fit 2D Gaussian. Defining samples within R 3 s minimizes the effects of cluster overlapping, which could be misleading in the estimation of age and age spreads in cluster-forming regions (e.g., NGC 1980, Alves & Bouy 2012). We obtained an extinction-limited sample, defined as follows. (1)We projected an A 20 V=mag extinction vector from a 1 Myr isochrone (D’Antona & Mazzitelli 1997)until it reached the sensitivity limit of the observations. The A V -limited sample is defined by all sources enclosed within the area defined by the extinction vector, the zero-age main sequence, and the reddened isochrone (see Figure 9, left).(2)We separated the possible foreground population from sources with colors coinciding with the locus of the dwarf and giant sequences in the J H -versus H K-color–color diagram. This defines a foreground contribution. In Figure 9(right panel) we show a color–magnitude diagram for W3-Main showing the A V -limited sample selection. We estimated the fraction of sources with circumstellar disks within R 3 s , for each principal cluster, as the total number of Class I plus Class II sources in the extinction-limited sample, down to a de-reddened brightness H=13.5 mag. (see also Section 3.2.3). The highest fractions of circumstellar disks, close to 30%, are in W3-Main, IC 1795, and W3(OH). Clusters NGC 896 and IC 1795-N have smaller fractions, closer to 15%. Table 2 Main Parameters and Age Estimates for Principal Clusters in W3 Cluster Center Coords. r cl Disk Fraction b Mean Age c Age Spread c J2000 (pc) a (%)(Myr)(Myr) W3-Main 36.368671 62.097485 1.313 27 ±4 2.0 [1.5, 2.5]1.5 [1.0, 2.5] IC 1795 36.652903 61.990908 1.244 29 ±5 3.0 [2.0, 4.0]4.0 [1.0, 5.0] W3(OH)36.798850 61.895349 1.644 27 ±3 3.0 [2.5, 3.5]3.5 [3.0, 5.0] NGC 896 36.388061 61.996607 1.428 14 ±3 3.5 [2.5, 4.5]4.0 [2.0, 5.0] IC 1795-N 36.569302 62.103874 1.618 12 ±2 4.0 [3.5, 5.0]2.0 [1.0, 5.0] Notes. a Assuming a distance d=2.045 kpc to the W3 Complex. b Estimated as the fraction of Class I+II sources within R 3s, in the extinction-limited sample, down to H dered =13.5 mag. c The number pairs inside the brackets indicate the range of model ages within the 2σ(68%)confidence level. The number listed before the bracket is the central value, which we list as the best estimate for the age/age spread of the cluster. Table 3 Near-Infrared Observations of W3 Fields Field ID Date Obs. Center Coords. Filter Seeing BD Peak a J2000 ( ) (mag) OMEGA 2000-CAHA 3.5 m OBSERVATIONS CAHA-W301 2008 Nov 08 36.580583 62.135000 J1.00 21.00 CAHA-W301 2008 Nov 08 36.580625 62.135000 H1.05 20.00 CAHA-W301 2008 Nov 07 36.580625 62.135000 K1.07 19.50 CAHA-W302 2008 Nov 09 37.100042 62.108056 J1.09 20.50 CAHA-W302 2008 Nov 08 37.100042 62.108056 H0.98 20.25 CAHA-W302 2008 Nov 09 37.100042 62.108056 K1.00 19.25 CAHA-W303 2008 Nov 09 36.659000 61.985833 J1.08 20.75 CAHA-W303 2008 Nov 09 36.659000 61.985833 H1.07 19.75 CAHA-W303 2008 Nov 09 36.658958 61.985833 K1.16 19.25 CAHA-W304 2008 Nov 09 37.056875 61.904167 J1.04 20.75 CAHA-W304 2008 Nov 09 36.580583 62.135000 H0.96 20.50 CAHA-W304 2008 Nov 09 37.056875 61.904167 K0.91 19.50 CAHA-W3CF b 2008 Nov 09 33.924208 60.545000 J1.07 21.00 CAHA-W3CF 2008 Nov 09 33.924208 60.545000 H1.14 19.75 CAHA-W3CF 2008 Nov 09 33.924250 60.545000 K1.06 19.50 Notes. a Brightness distribution peak. Expresses the turnover point of the observed magnitude distribution. b Control Field. 8 The Astronomical Journal, 150:80 (18pp), 2015 September Román-Zúñiga et al. The total contribution from the field was estimated from a nearby control field with minimum extinction, observed in the same conditions and to the same depth as the cloud fields (see the appendix). We made a K-band histogram from the control field catalog and scaled it to the same area of the cluster sample. Then we subtracted the foreground contribution from the control field KLF and applied artificial reddening to the remaining distribution down to the same level as the cloud sample. Artificial reddening is applied by constructing a normalized A V extinction distribution from the background stars in each cluster region using the NICER algorithm, and used it as a frequency distribution to convolve the control field K-band histogram, producing a final control field K- band brightness distribution. The NICER algorithm provides estimates of extinction with average errors 0.1 AV s < mag below A 30 V=mag. These errors are way below the size of the K-band brightness bins used to construct the KLF. The final KLF of the cluster sample is obtained by subtracting the final A V -convolved (reddened)control field distribution from the one constructed from the K-band limited sample. The advantages of using a K-band limited sample to construct the KLF are (a)that all samples are cut at a common maximum extinction, (b)that the contribution of background extragalactic sources is reduced to a minimum, and (c)that patchy extinction is less prone to biasing the brightness distribution, especially at the most embedded regions—e.g., the centers of W3-Main and W3(OH). Also, by convolving the raw control field KLF by the “on”field extinction distribution, we are able to obtain a fair estimation of the contribution from field stars on different lines of sight (Muench et al. 2003). 3.3.1. Age Inference from Cluster Model Fits The KLF for each principal cluster sample was then compared to artificial K-band brightness distributions of embedded cluster populations of varying ages, constructed with the model interpolation code of Muench et al. (2000). The method is relatively easy to apply: using a set of pre-main sequence models, the code uses a Monte Carlo generator to draw an artificial population from an assumed initial mass function (IMF)with a total number of members input by the user, between a minimum and a maximum age, defined by model isochrones, properly shifted by the distance modulus of the complex. 13 For this study, we chose the isochrone set of D’Antona & Mazzitelli (1997), with a deuterium fraction DH 210 5 [ ]=´ -. The user has the option to choose the form of the IMF and also has the option to fix or to sample values for a set of parameters that includes an extinction distribution, a binary fraction, and the infrared excess fraction as a function of color ( H K-). The peak and the slopes of the KLF are particularly sensitive to age and age spread, which are the two parameters that were left to vary. We decided to leave the binary fraction as a fixed value of 20% based on estimations by Padgett et al. (1997)suggesting that the binary fraction in the Orion B clusters may be about 15% in separations up to 1000 AU—the separations we are able to resolve in our nearinfrared images. Also, Muench (2002)in a similar experiment showed that, in particular, the mean age and spread estimates Figure 7. RGB image constructed by combining the maps from Spitzer MIPS 24 μm(blue),Herschel PACS 70 μm(green), and Herschel SPIRE 250 μm(red). Contours indicate the surface density of Class II young stellar sources in the W3 region, estimated with the gather method. Contour levels are set as 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, and 10.0 stars per square arcmin. The positions of Class I sources are indicated with red circle symbols. The positions of O- and B-type stars are indicated with magenta and cyan symbols, respectively. Areas in yellow/green are saturated pixels in the infrared images. The labels correspond to previously identified clusters and sub-clusters in the W3 region; labels starting with letters “RI”indicate those identified by Rivera-Ingraham et al. (2011). 13 The 0.107 kpc distance error quoted by Hashisuka et al. (2006)for their maser parallax distances translates into an error of about 0.2 mag in the brightness estimates. This is smaller than the size of our KLF bins (0.5 mag).In terms of mass, this translates into an uncertainty of about 0.05 in a log-mass scale. 9 The Astronomical Journal, 150:80 (18pp), 2015 September Román-Zúñiga et al. rather than a process of induction in a sequence from the older to the younger episodes. Our YSO surface density analysis provides evidence of substructure in the principal cluster populations. Small aggregations of stars around the main cluster events (very evident in W3 (OH)and groups like R102 and R103)were revealed to be significant peaks in the surface density of YSOs, in many cases coincident with the locations of Class 0/I candidates. This is suggestive of a complex scenario in which stellar collections of different sizes may form with small age differences by the compression of gas in various locations within large clumps. Strong evidence of sub-clustering in W3-Main was also provided by the very deep near-infrared KLF by Bik et al. (2014): they pointed out that an age gradient was evident from the edge to center in W3-Main, judging from a higher fraction of circumstellar disks near the center. The large number of disk bearing sources in IC 1795 could be indicative of a similar process, where cluster formation occurs in substructures. If these substructures do not form simultaneously, they certainly may lead to a large age spread, as observed. The projected separations between surface density peaks in the gather map range is 1.4 ±0.4 pc. These separations are comparable to the estimated radii of the principal clusters identified with our GMM model. This scale is shorter than the Jeans length (2.0 pc)for a mean density of 10 cm 23 and T10 K=K, typical of the diffuse gas traced by 12 CO, and may be tracing the primordial fragmentation of the cloud. One possibility is that density peaks in the gather map are subclusters that eventually merge into larger groups, but that will depend on the dynamical stability of the aggregations and the gravitational potential of the corresponding principal cluster. Numerical simulations (Allison et al. 2010; Bonnell et al. 2011) suggest that clusters are initially highly structured, but rapidly evolve by merging into large entities. On the other hand, small groups of stars may disperse rapidly or even remain stable depending crucially on the number of members and the starforming efficiency (Adams & Myers 2001; Lada & Lada 2003). It is also expected that most small, embedded aggregations end up being disrupted early, as proposed by the “infant mortality” and “cruel cradle”models. Figure 14. Cartoon illustrating our proposed scheme for the history of star formation in the W3 complex. The blue regions indicate the extension of the molecular gas cloud. Dotted, empty regions indicate regions of collapse toward cluster formation. Red shaded areas indicate active cluster formation regions. Green shaded areas indicate cluster regions where gas dissipation has started and yellow shaded areas indicate slightly older clusters where gas dispersal has removed a significant amount of the parental gas and the cluster. The rectangular easternmost area represents a region of earlier formation with unknown extension that might contribute today to a larger age spread observed in the W3(OH)cluster complex. Black arrows indicate the suggested progression of cluster formation in the complex. 16 The Astronomical Journal, 150:80 (18pp), 2015 September Román-Zúñiga et al. The formation of small groups that disperse rapidly and merge into a larger central structure could provide a simple explanation for the presence of larger age spreads in a lower surface density population in regions like IC 1795-N or NGC 896, and would be in agreement with age spreads and substructure in the younger principal clusters. Numerical studies like those of Banerjee & Kroupa (2015) suggest that scenarios that lead to the formation of centrally condensed, massive clusters from the merging of smaller subclusters may not be the most common but may exist under certain conditions, particularly that stellar overdensities have ∼1 pc scale separations and may form via the process of local triggering. This is in good agreement with our observations. The strong outflows in the molecular gas shell at W3-Main and W3(OH)confirm a very rapid removal of gas in the complex. This removal may occur in a period comparable to the T Tauri timescale, as evidenced by Class II source candidates in IC 1795, where gas is mostly dispersed. However, it is important to stress that the mass loss rate cannot be constant, and will probably decline within a few 10 5 Myr, if it has not already done so. The estimate we provided, however, may not be far from the average mass loss rate during the main gas dissipation period. 5. CONCLUSIONS We obtained deep, high-resolution ground-based nearinfrared images of the eastern edge of the W3 region in the Perseus Arm. These data were combined with archive imaging data from the Spitzer and Chandra space telescopes to obtain a new catalog of young sources for the W3 Complex. We used our resultant catalog to: 1. Identify and classify young (Classes I–III)embedded sources in the region. 2. Determine the locations and extension of the principal clusters in the complex. 3. Determine the spatial distribution of young sources across the complex. This distribution shows significant substructure. 4. Construct extinction-limited samples that allowed us to obtain clean KLFs of the principal clusters. We compared the resultant functions with those from artificial clusters of ages created from pre-main sequence models. This analysis permitted to estimate mean ages and possible age spreads for the principal clusters of the W3 Complex. We found that IC 1795, the central cluster of the complex, still has a large fraction of Class II sources, despite not being embedded anymore, like clusters W3-Main and W3(OH). The R II:III ratio is large across IC 1795, just as it is in the flanking embedded clusters. The 13 CO map shows some receding molecular gas clumps coincident with the center of IC 1795. This, along with the absence of Class I sources, suggests that parental gas was removed quickly from IC 1795 in a period shorter than the T Tauri timescale. The high mass loss rate estimate at the center of W3-Main strongly supports this. Our GMM analysis clearly identified the well known clusters IC 1795, W3-Main, and W3(OH). In addition, two principal clusters could be identified, one associated with the compact HII region NGC 896 and one located more north of IC 1795. We estimated the ages and age spreads of the principal clusters by comparing their KLFs with those of artificial young clusters. The differences between the ages of the clusters are not significantly large. We also were able to constrain the age spread of the W3 cluster, but we could not well constrain the age spread of the other four clusters, so there might be overlaps in the formation periods of all clusters in the complex. This is confirmed by a large number of Class II and III sources across the whole complex. This suggests a progression of the cluster formation scenario in which formation proceeded as a result of the primordial density structure of the cloud. Using our list of Class II candidates, we constructed a surface density map that shows a significant amount of substructure in the principal clusters of the complex. The large age spread in the clusters may indicate that formation proceeds in episodes. If the observed substructures do not form simultaneously, they certainly may lead to a large age spread, as we observe. We propose a scenario in which star cluster formation progressed over 4 Myr ago from the north and south (clusters NGC 896 and IC 1795-N)toward IC 1795. Later (about 2–3 Myr ago), the progression moved from the center (IC 1795)toward the east and west (W3(OH)and W3-Main).We also suggest that an undefined region east of W3(OH)could form groups of stars before the formation of the multiple embedded groups we see today. The formation of W3(OH)and W3-Main is still in progress at the edge of the giant molecular shell of the complex. Using molecular emission maps, we showed that the two main gas clumps associated with these clusters (which together contain more than M 4 .5 103 ´)show strong velocity gradients and gas outflow activity, spatially coincident with clusters and substructures. This is in good agreement with the convergent flow scenario (Rivera-Ingraham et al. 2013), which contributes to a large massive star formation efficiency in W3. Our analysis suggests that both the progression of formation following the primordial structure of the original GMC and the triggered formation caused by the compression and dispersal of gas (that is, within the cloud)are both important in W3. The W3 complex presents a rather complicated layout. The principal clusters overlap partially along the line of sight and they may still interact with each other. Moreover, the clusters present a significant substructure. If the substructures are bona fide sub-clusters, they pose an interesting problem for future higher quality data. Particularly, space-based mid-infrared imaging with high resolution and deeper photometric limits (e.g., the James Webb Space Telescope)will provide more complete samples of the sub-cluster population. Also, spectroscopic data for individual young candidate members may refine age and age spread estimates. We thank an anonymous referee for comprehensive and useful comments that greatly improved the content of our manuscript. C.R.Z. acknowledges support from CONACYT project CB2010-152160, Mexico and program UNAMDGAPA-PAPIIT IN103014. J.Y. and M.T. acknowledge support from program UNAM-DGAPA-PAPIIT IN101813. G.M.V. acknowledges support from an introduction to research scholarship from the Programa de Introducción a la Investigación de la Junta de Ampliación de Estudios (JAE INTRO) CSIC/Instituto de Astrofísica de Andalucía. E.A.L. acknowledges support from the National Science Foundation through NSF Grant AST-1109679 to the University of Florida. This study is based on observations collected at the Centro Astronómico Hispano Alemán (CAHA)at Calar Alto, operated jointly by the Max-Planck Institut für Astronomie and the 17 The Astronomical Journal, 150:80 (18pp), 2015 September Román-Zúñiga et al. Instituto de Astrofísica de Andalucía (CSIC). We acknowledge the staff at Calar Alto for top of the line queued observations at the 3.5 m with OMEGA 2000. We acknowledge use of data products from 2MASS, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Centre/California Institute of Technology (funded by the USA National Aeronautics and Space Administration and National Science Foundation). This work is partly based on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. The scientific results reported in this article are based to a significant degree on data obtained from the Chandra Data Archive; particularly, we made use of data obtained from the CSC, provided by the Chandra X-ray Center (CXC)as part of the Chandra Data Archive. In this paper we discuss observations from the Bolocam Galactic Plane Survey (BGPS; PI John Bally); the BGPS project is supported by the National Science Foundation through NSF grant AST-0708403. In the present paper we discussed observations performed with the ESA Herschel Space Observatory (Pilbratt et al. 2010), in particular employing Herschelʼs large telescope and powerful science payload to do photometry using the PACS (Poglitsch et al. 2010)and SPIRE (Griffin et al. 2010)instruments. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. This research made use of Montage, funded by the National Aeronautics and Space Administrationʼs Earth Science Technology Office, Computation Technologies Project, under Cooperative Agreement Number NCC5-626 between NASA and the California Institute of Technology. Montage is maintained by the NASA/IPAC Infrared Science Archive. We made use of the pvextractor tool by Adam Ginsburg, which is part of the Radio Astro Tools repository (http://github.com/radio-astro-tools). Facilities: CAO:3.5 m, Spitzer, CXO, CSO, Herschel, ARO:12 m. APPENDIX NEAR-IR OBSERVATIONS In Table 3we list information corresponding to our nearinfrared Calar Alto 3.5 m OMEGA 2000 observations. REFERENCES Adams, F. C., & Myers, P. C. 2001, ApJ,553, 744 Aguirre, J. E., Ginsburg, A. G., Dunham, M. K., et al. 2011, ApJS,192, 4 Allison, R. J., Goodwin, S. P., Parker, R. J., Portegies Zwart, S. F., & De Grijs, R. 2010, MNRAS,407, 1098 Alves, J., & Bouy, H. 2012, A&A,547, A97 Banerjee, S., & Kroupa, P. 2015, MNRAS,447, 728 Banfield, J. D., & Raftery, A. E. 1993, Biometrics,49, 803 Bertin, E., & Arnouts, S. 1996, A&AS,117, 393 Bieging, J. H., & Peters, W. L. 2011, ApJS,196, 18 Bik, A., Henning, T., Stolte, A., et al. 2012, ApJ,744, 87 Bik, A., Stolte, A., Gennaro, M., et al. 2014, A&A,561, A12 Bonnell, I. A., Smith, R. J., Clark, P. C., & Bate, M. R. 2011, MNRAS, 410, 2339 Bressert, E., Bastian, N., Gutermuth, R., et al. 2010, MNRAS,409, L54 Casertano, S., & Hut, P. 1985, ApJ,298, 80 Cottaar, M., Covey, K. R., Meyer, M. R., et al. 2014, ApJ,794, 125 D’Antona, F., & Mazzitelli, I. 1997, MmSAI, 68, 807 de Grijs, R. 2011, in Stellar Clusters, Associations: A RIA Workshop on Gaia, ed. E. Alfaro, T. Gallego & M. R. Zapatero-Osorio (Granada, Spain: IAACSIC), 331 Evans, I. N., Primini, F. A., Glotfelty, K. J., et al. 2010, ApJS,189, 37 Feigelson, E. D., & Townsley, L. K. 2008, ApJ,673, 354 Fraley, C., & Raftery, A. E. 2002, J. Am. Stat. Assoc., 97, 611 Fraley, C., Raftery, A. E., Murphy, T. B., & Scrucca, L. 2012, Tech. Rep. No. 597 (Department of Statistics, Univ. Washington) Ginsburg, A., Bally, J., & Williams, J. P. 2011, MNRAS,418, 2121 Gladwin, P. P., Kitsionas, S., Boffin, H. M. J., & Whitworth, A. P. 1999, MNRAS,302, 305 Griffin, M. J., Abergel, A., Abreu, et al. 2010, A&A,518, L3 Gutermuth, R. A., Megeath, S. T., Muzerolle, J., et al. 2004, ApJS,154, 374 Gutermuth, R. A., Megeath, S. T., Myers, P. C., et al. 2010, ApJS,189, 352 Gutermuth, R. A., Myers, P. C., Megeath, S. T., et al. 2008, ApJ,674, 336 Hachisuka, K., Brunthaler, A., Menten, K. M., et al. 2006, ApJ,645, 337 Kruijssen, J. M. D., Maschberger, T., Moeckel, N., et al. 2012, MNRAS, 419, 841 Lada, C. J., Elmegreen, B. G., Cong, H.-I., & Thaddeus, P. 1978, ApJL, 226, L39 Lada, C. J., & Lada, E. A. 2003, ARA&A,41, 57 Lada, C. J., Lombardi, M., Román-Zúñiga, C., Forbrich, J., & Alves, J. F. 2013, ApJ,778, 133 Levine, J. 2006, PhD thesis, Univ. Florida Lombardi, M. 2009, A&A,493, 735 Muench, A. 2002, PhD thesis, Univ. Florida Muench, A. A., Lada, E. A., & Lada, C. J. 2000, ApJ,533, 358 Muench, A. A., Lada, E. A., Lada, C. J., et al. 2003, AJ,125, 2029 Navarete, F., Figueredo, E., Damineli, A., et al. 2011, AJ,142, 67 Oey, M. S., Watson, A. M., Kern, K., & Walth, G. L. 2005, AJ,129, 393 Ogura, K., & Ishida, K. 1976, PASJ, 28, 651 Padgett, D. L., Strom, S. E., & Ghez, A. 1997, ApJ,477, 705 Peng, T.-C., Wyrowski, F., van der Tak, F. F. S., Menten, K. M., & Walmsley, C. M. 2010, A&A,520, A84 Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A,518, L1 Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A,518, L2 Qiu, K., Zhang, Q., Wu, J., & Chen, H.-R. 2009, ApJ,696, 66 Rivera-Ingraham, A., Martin, P. G., Polychroni, D., & Moore, T. J. T. 2011, ApJ,743, 39 Rivera-Ingraham, A., Martin, P. G., Polychroni, D., et al. 2013, ApJ,766, 85 Román-Zúñiga, C. G. 2006, PhD thesis, Univ. Florida Román-Zúñiga, C. G., Elston, R., Ferreira, B., & Lada, E. A. 2008, ApJ, 672, 861 Ruch, G. T., Jones, T. J., Woodward, C. E., et al. 2007, ApJ,654, 338 Schneider, N., Csengeri, T., Bontemps, S., et al. 2010, A&A,520, A49 Stetson, P. B. 1987, PASP,99, 191 Taylor, M. B. 2005, in ASP Conf. Ser. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton & R. Ebert (San Francisco, CA: ASP),29 Thronson, H. A., Jr., Lada, C. J., & Hewagama, T. 1985, ApJ,297, 662 Wang, J., Feigelson, E. D., Townsley, L. K., et al. 2009, ApJ,696, 47 Ybarra, J. E., Lada, E. A., Roman-Zúñiga, C. G., et al. 2013, ApJ,769, 140 18 The Astronomical Journal, 150:80 (18pp), 2015 September Román-Zúñiga et al.