YOUNG STELLAR OBJECT VARIABILITY (YSOVAR): LONG TIMESCALE VARIATIONS IN THE MID-INFRARED

In this section, we first discuss our criteria for picking targets and then review basic properties of each cluster.

2.1.1. Overview of Target Selection

2.1.2. Cluster Properties

To better manage the observation planning and data downloading for the 12 clusters we observed, we separated each cluster into an individual observing program. The individual observing programs and clusters are listed in Table 2; one can download the data from the Spitzer Heritage Archive using these program numbers. We will deliver all of our extracted photometry to the Spitzer Science Center (SSC) and Infrared Science Archive (IRSA²⁹) for public access via IRSA tools, and through those tools, the Virtual Observatory.

Most of the observations were IRAC mapping mode Astronomical Observation Requests (AORs) using full-array 12 s high-dynamic-range (HDR) mode (which is defined such that a 0.4 and 10.4 s exposure is obtained at each pointing—see the IRAC Instrument Handbook, available at the SSC/IRSA Web site³⁰). In some cases, as noted in Table 2, the AORs were staring AORs, meaning that they were continuous or semi-continuous observations without dithering or mapping. These staring data will be discussed in other YSOVAR papers. The present paper is limited to the HDR mapping observations (and is largely further restricted to the "fast cadence" observations; see Section 2.4).

Roughly half of the original YSOVAR time allocation was devoted to observations of Orion. The Orion Spitzer campaign included a ∼0.9 deg² region centered on the Trapezium cluster in Orion; this is far larger than can be obtained in a single AOR at a single epoch. The observed area was thus broken into five segments with a central region of ∼20' × 25' and four flanking fields. The central part was observed in full array mode with 1.2 s of exposure time and 20 dither positions to avoid saturation by the bright nebulosity around the Trapezium stars. The remaining four segments of the map were observed in 12 s HDR mode, and four dither positions. More details on these original YSOVAR Orion observations appear in MC11. We reallocated time within the YSOVAR time budget to follow up on some eclipsing binaries in Orion (Morales-Calderón et al. 2012). We also obtained time in Cycle 7 to follow up on some AA Tau analogues in Orion presented in MC11. These follow-up observations did not map the entire Orion region.

The rest of the original YSOVAR time allocation was divided among the rest of the SFRs. In contrast to Orion, for the 11 other clusters, the monitoring observations are one or at most a few IRAC FOVs; see the last column of Table 2. As noted above, we sometimes refer to these 11 other clusters as "the smaller field clusters."

The observations are typically spread over weeks to months, usually at a cadence of about twice per day, but with an interval between observations that varied both by design (to reduce aliasing problems) and due to constraints imposed by other Spitzer programs or infrastructure operations that were executed during the same campaigns (further details on the cadence follow below in Section 2.4).

NGC 2264 was included as a smaller-field cluster to be monitored as part of the original YSOVAR program, and while this original program was still executing, CSI 2264 was approved. We continued to execute the original YSOVAR small-field observations in this region (program 61027 in Table 2), and these small-field regions are discussed here, since they resemble the rest of the original YSOVAR programs more closely than they resemble CSI 2264.

Two of our clusters, IC 1396A and Serpens Main, were monitored in the cryogenic era with Spitzer with the primary intention of monitoring changes in these objects in the Spitzer bands (as opposed to removing artifacts). In YSOVAR, we re-observed these clusters in the same way so that the same objects continue to be monitored, and we re-reduced all of these data in the same fashion as discussed here. In the context of this paper, we are not particularly focused on these cryogenic-era monitoring data; however, these data will be included in the YSOVAR cluster-specific papers.

Because of the nature of Spitzer and IRAC, the footprints of our observations and how they vary with time are both complicated issues. We now discuss these issues as they pertain to YSOVAR observations.

Figures 1–13 present the outline (footprint) of the YSOVAR observations. The original YSOVAR Orion observations (Figure 3) and the CSI 2264 observations (Figure 7) cover a substantially larger area than those for the other smaller-field clusters. Chandra footprints are included in these figures for reference and discussed below in Section 2.10.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The focal plane of the IRAC camera is such that the 3.6 and 4.5 μm FOVs are not the same; data are obtained in both FOVs at once, with one field placed on the target of interest, and the other obtaining serendipitous data in a non-overlapping ∼5' × 5' field with a center offset ∼65 from the target field; see the IRAC Instrument Handbook for more details. The placement of these FOVs also changes with time. For all Spitzer observations, as discussed in the SpitzerSpace Telescope Handbook (also available at the SSC/IRSA Web site³¹), the ecliptic latitude of the target defines when one can observe the target and for how long. At any one time, Spitzer can observe in an annulus defined by the operational pointing zone, which can be conceptually summarized as "neither too close nor too far away from the Sun." It is about 40° wide, and rotates with the Sun at a rate of about a degree a day. An object near the ecliptic plane can thus only be observed for a period of about 40 days twice a year; objects near the ecliptic pole can be observed at any time during the year for as long as needed. Related to this, the IRAC FOV, as projected onto the sky for any given object on the ecliptic equator, is at an essentially constant angle with respect to north for the duration of the ∼40 day observing window, and, ∼6 months later, is flipped by 180° but then also essentially constant for that ∼40 day observing window. However, the FOV for an object at the ecliptic pole rotates by about a degree a day.

For each of the figures showing the IRAC footprints (Figures 1–13), the projected outline of the observation covered by the entire YSOVAR data set is indicated on top of a 4.5 μm observation obtained (in most cases) during the cryogenic era. The different regions observed by the 3.6 and 4.5 μm cameras are identified; the nominal target of the observation is covered in both FOVs. In addition to the target of the observations, there are serendipitous data obtained offset from the target area, as seen in the figures. The footprint of the serendipitously obtained data in each band can change considerably depending on the time of observation and the ecliptic latitude of the target. For targets at higher ecliptic latitudes, where field variation with time is important, a single epoch of observation is indicated in the figure in addition to the area covered by the entire YSOVAR data set.

The approximate ecliptic latitude for each target is included in Table 2 as a rough indication of how long the observing window was and how much field rotation one can expect to see in these figures. For Orion, we designed our observations to be less sensitive to rotation; that plus the fact that the Orion map contains many FOVs means that the rotation is much less of a factor in terms of how much of the sky may be included in both channels (as opposed to only one). For the smaller-field cluster observations containing only one to four IRAC FOVs, a primary target field is monitored in both IRAC bands, but substantial serendipitous monitoring data have also been obtained, typically in just one band. For most of the smaller-field YSOVAR cluster targets, such as NGC 1333 (Figure 2), the central region has complete two-band coverage and most of the objects in the region of the single-band coverage north and south of the target field have data over most, if not all, of the campaign. For L1688 (Figure 8), with an ecliptic latitude of −3°, the field rotation is essentially zero during one ∼40 day window. However, it was monitored over more than one window, so the observations flip by 180° for the next ∼40 day window. (The targets of observation in Figure 8 are the three regions with brightest stars and nebulosity in the centers of the three "stripes" of coverage.) In contrast, for IC 1396A or Ceph C (ecliptic latitudes of ∼60°), the central 5' × 5' region is covered throughout the monitoring window, but the rapid field rotation creates a "fan" of coverage such that the serendipitous coverage in the outer periphery of those targeted regions have data in just one band, and only for a fraction of the campaign. Some objects are monitored in 3.6 μm at the beginning of the campaign, and in 4.5 μm at the end of the campaign (see Figures 12 and 13).

If one considers only objects with light curves in both channels for the smaller-field clusters, even for the fields that do not rotate, one loses a significant fraction of the available data; while the nominal targets of our observations are the regions with two-band coverage, there are still good light curves for cluster members in the regions with coverage in only one band. Generally, for the smaller-field clusters, we did not require data in both channels, and instead retained all light curves for our analysis, even if obtained in only one channel (Section 3.2). Note also that the highest ecliptic latitude among our targets is ∼65°; no higher ecliptic latitude targets were included, despite the possibility of much longer observing windows, at least in part because we could not find any suitable SFRs at these latitudes at the time we planned our observations. The field rotation for these targets at ∼65° is already substantial.

In general, because of the field rotation, as any given object moves into (or out from) the mapped region, the first (or last) few epochs may be unusable as the object may appear in only a fraction of the dithers at that position and the photometry may thus be compromised by edge effects. Substantial field rotation with time can have a significant effect on faint stars near bright stars. Diffraction spikes of bright stars rotate with time within the IRAC FOV, so if a target star falls near a diffraction spike, this can introduce a source of false longer timescale variability; see also Section 2.5.

We note here that the low ecliptic latitude of L1688 (and specifically the 180° flip between observing windows) enables Günther et al. (2014) to use these observations to characterize possible artifacts in the light curve related to the position in each frame. These tests confirm that the statistical criteria for selecting variable sources (see Section 5 below) is restrictive enough that differences between positions within the frame cannot be a primary contributor for the sources identified as variable.

More than 95% (524 hr) of our original YSOVAR Spitzer Cycle 6 program observing time was devoted to a "fast cadence" mode, designed to be sensitive to timescales from ∼0.15 to ∼40 days, consistent with the known timescales of YSOs due to accretion-related flickering, rotational modulation of star spots, and other effects. The upper limit to this range of timescales was set by the typical duration of Spitzer visibility windows near the ecliptic plane, as discussed above in Section 2.3.

Sampling a light curve at evenly spaced intervals introduces period aliasing, and a bias toward the detection of false periods at integer fraction multiples of the sampling interval (see, e.g., Plavchan et al. 2008b or Dawson & Fabrycky 2010). This period aliasing is common at integer fraction multiples of one day in ground-based photometric surveys due to the natural day-night cycle of the Earth's rotation. Given typical YSO rotation periods of one to several days, we therefore chose to execute our program with a cadence designed to be compatible with the Spitzer scheduling, and minimize period-aliasing while retaining sensitivity to variability on a broad range of timescales.

The fast cadence mode for each region had a total time baseline of ∼40 days, varying depending on the actual duration of the visibility window, with additional days on either end of the visibility window for scheduling flexibility. The actual total duration of the fast cadence observations for each cluster is listed in Table 3. In particular, for AFGL 490, Mon R2, GGD 12-15, IRAS 20050+2720, IC 1396A, and Ceph C, we originally planned a 42 day timespan; for Orion, Serpens Main, and Serpens South, we planned a 38.5 day timespan; and for NGC 1333, NGC 2264 and L1688, we planned a 35 day timespan. We divided the time baseline into a repeating set of 3.5 day sub-cadences to ease scheduling. Within each 3.5 day period, we made eight visits, with time steps between visits of 4, 6, 8, 10, 12, 14, and 16 hr, respectively. By using a linearly increasing time step, we were able to evenly sample in Fourier space a range of higher frequency (shorter timescale) photometric variability than we would have been able to otherwise. Additionally, we were also able to minimize the total amount of necessary observing time to sample these high frequencies and to minimize the period aliasing. For comparison, splitting 8 visits evenly across a 3.5 day cadence would have removed our sensitivity to variability timescales shorter than ∼0.5 day.

To further accommodate scheduling flexibility, we also included a window of ±2 hr in which to obtain a given epoch of observation for our fast cadence observations. This ±2 hr window effectively resulted in a randomization of the cadence about our desired times of observation, which had the additional benefit of further reducing period aliasing, in particular aliases with periods of 3.5 days. We also allowed for interruptions due to data downlink transfers and higher priority observations such as staring mode observations of transiting exoplanets. We are grateful to the SSC scheduling staff for accommodating this fairly complex cadence over the two years of the survey program.

Overall, we obtained ∼40–100 epochs per cluster, which are visualized in Figure 14. Table 2 lists all the clusters, with information about the observations (including the total number of epochs); Table 3 summarizes only the fast cadence observations, including typical values of the minimum and maximum size of the timestep between epochs, and the total time baseline. Our unevenly spaced epochs enable a characterization of different types of observed variability on a variety of timescales; for further discussion, see below (Section 5.4). While periodic analysis tools take advantage of such unevenly spaced data, Findeisen et al. (2014) point out that for auto-correlation functions (ACFs) of non-periodic light curves, timescale sensitivity may be limited by the largest (not the smallest) adjacent timestep in the time series.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For three SFRs (IC 1396A, L1688, and Ceph C) with longer visibility windows, or that we observed in multiple visibility windows, we also observed in a "slow cadence." Besides the 42 day fast cadence, we executed visits with time steps between visits of 1, 2, 3, 4, 5, etc. days to cover the remaining duration of the visibility window. These observations constituted ∼5% of our Cycle 6 observing program (25 hr), and these observations allowed us to be sensitive to variability timescales of ∼40 days to ∼2 yr, e.g., to objects such as KH15D (e.g., Winn et al. 2006) and WL4 (Plavchan et al. 2008a), as well as long term trends on timescales of a year identified in the NIR (e.g., Parks et al. 2014). For L1688 and Ceph C, we also executed this slow cadence observing mode during the three other visibility windows during our survey. Additionally, Orion had some slower cadence follow up, and NGC 2264 had the CSI 2264 program at a different cadence.

For the primary statistical analysis discussed in this paper, we use the "standard statistical sample," e.g., only the fast cadence data (where there are at least five epochs in the light curve; see Section 3 below). For about half of the clusters, this is the entirety of the data set; see Tables 2 and 3, and the red points in Figure 14. A detailed analysis of the ensemble of all of these observations for each cluster will be presented in papers in preparation, customized to each cluster.

For completeness, we note here the following specifics about those clusters with data beyond the YSOVAR fast cadence. L1688, IC 1396A, and Ceph C all have slow cadence data. Orion has slower cadence data over a second year. Data for Ceph C over its fast cadence window include those originally obtained as part of the original YSOVAR project, as well as data taken using a very similar distribution of time steps as part of a YSOVAR-related Chandra program; see Table 2. For consistency with the maximum timescales sampled in other clusters, we define the Ceph C fast cadence light curve to be a window spanning from 2010 September 19 through 42 days later, which is more representative of the fast cadence window in the other clusters. We note, however, that the sampling within this window is nonetheless enhanced relative to other smaller-field clusters, due to the inclusion of extra Spitzer observations obtained to support the coordinated Chandra observations. Two other clusters (IC 1396A and IRAS 20050) also sample longer timescales than the more typical ∼40 day window (see Figure 27 below for an enlargement of the fast cadence). We did not truncate the last few slower cadence points from the IRAS 20050 time series because there are no other slower cadence observations, and truncating it would effectively "orphan" the last few points. The IC 1396A cadence does not provide a clean breakpoint, so it also has points covering a slightly longer window than "typical" included in its fast cadence (see Table 3 for the total time in the fast cadence for all clusters). IC 1396A was one of the first SFRs to be monitored intensively with Spitzer (Morales-Calderón et al. 2009), so it has additional cryo-era monitoring not shown in Figure 14. Serpens Main also has some cryo-era monitoring. GGD 12-15 and L1688 also have staring data, not included here.

We started with the IDL package Cluster Grinder (G09), which has been used by several other projects (e.g., Megeath et al. 2012), and then made custom modifications to make it suitable for this time series data set. We started with the basic calibrated data (BCD) images released by the SSC. The BCD images used correspond largely to software version S18.18 (with about 20% S19.0 and S19.1); earlier reductions such as those in MC11 used a mix of S18.12, 18.14, and 18.18. Each BCD frame was processed for standard bright source artifacts and combined into mosaics (at 086 per pixel grid resolution) by exposure time (long or short), epoch, and bandpass. Cosmic ray hits and other transient artifacts are flagged during mosaic construction using redundant data at each pixel position in the final mosaic grid as a reference for the nominal value and allowable internal variation.

Since the HDR mode obtains a long and a short frame at each pointing, we need to combine data from the two exposures. HDR mosaics are merged together on a pixel-by-pixel basis, with appropriately scaled short-frame HDR values above a set threshold (individual pixel values corresponding to sources ∼10th mag, dependent on background level) replacing significantly compromised pixels in the long-frame HDR mosaic.

Point source detection and aperture photometry are then performed for each channel and each AOR/epoch. The aperture radius was 24, with a sky annulus having inner and outer radii of 24 and 72, respectively. We adopted Vega standard magnitude zero points of 19.30 and 18.64 for 1 DN s⁻¹ total flux at 3.6 and 4.5 μm, respectively. These values include standard corrections for our chosen aperture and sky annulus sizes (Reach et al. 2005; Carey et al. 2010).

Extensive tests were run on aperture photometry obtained in three different ways. Our initial approach was to obtain photometry from the mosaic created for each AOR, with the assigned time for that point being the average time for the AOR. We also obtained photometry on each BCD frame individually, with the individual, specific time of that observation being assigned to that point. Finally, we took the set of all individual BCD measurements for a given source within a given AOR, and took the mean brightness and the mean time for that set of observations. This latter approach proved to have the best quality photometry. It facilitated the treatment of two well-known systematic effects (residual gain and pixel phase effects) while keeping our signal-to-noise ratios (S/Ns) close to ideal. It provided the highest S/N light curves, even if the time resolution is slightly lower than what might be theoretically possible. Effectively, this choice means that a S/N ∼ 5 is required for a viable source in a given BCD, such that the net S/N of the source over the entire AOR is at least ∼10. If a source fell on the edge of the BCD such that the two native pixel radius was not entirely within a given BCD, no measurement is obtained. Measurements were also excluded as non-viable if the measurement (in the aperture or annulus) included a masked IRAC pixel such that obtaining a finite measurement was not possible, or if the source was faint enough (or the sky variation high enough) that the net measured flux was <0. If there were three or more valid detections of a given source in a given AOR, we could perform outlier rejection, and outliers were flagged; the rest of the measurements were combined by uniform weight arithmetic mean. The uncertainties were added in quadrature and divided by the included measurement count.

Source matching between epochs and wavelengths was performed by position, taking the nearest source within 1''. Final positions for each source were the mean of the individual measurements. Our astrometric residuals as compared to 2MASS suggest an uncertainty of <200 mas root-mean-square (rms), after a single iteration refinement of the astrometry using a first set of mosaics constructed blindly with the BCD-delivered world coordinate system (WCS). Objects with fewer than five viable detections over separate epochs in a single band were not retained as valid light curves.

We note here that there are some residual instrumental column pulldown effects (see the IRAC Instrument Handbook for more information) in some areas near bright sources. This is a more significant problem for clusters with a large number of very bright sources (e.g., Mon R2; see brightnesses in Figures 1–13). The locations of these artifacts can change with time, particularly in fields at high ecliptic latitude and thus with significant field rotation during a given visibility window (see Section 2.3), and can have particularly significant effects on sources fainter than ∼12th mag. Objects falling in such regions that are faint and that had non-repeating structure in their light curves should be particularly scrutinized (including visual inspection of the input frames) for the validity of their light curves. Points obviously compromised by these effects should additionally be identified as non-viable points and rejected. Discussions of this process for individual objects will appear in the clusters papers to follow.

Even after all of the processing to this point, some occasional outliers in the light curves remain. To identify outliers within the fast cadence data, we initially identified points several σ away from the mean for each light curve. We used this approach in MC11, with a different (much earlier) processing of the original Spitzer data. However, our more recent, improved processing described here reduced the number of outliers. We experimented with algorithms for rejecting outliers from each light curve, but were unable to identify a technique that rejected visually spurious points while preserving apparent bona fide variability. An approach as we used in MC11 would, for example, exclude points at the beginning and end of a light curve with a long trend. Our estimates suggest that implementing an outlier rejection strategy affects ≲10% of the light curves tagged as intrinsically variable (as per methods described further below), so we retained all points in the light curves. Outliers may still be identified in a few specific individual cases; those outliers are then effectively excluded from the light curve, and are discussed when relevant.

Many tests of variability (e.g., a χ² test; see Section 5.2) depend on accurate error estimates for each data point in the time series. Initial estimates of uncertainties were derived by combining three terms in quadrature (see also the discussion in Megeath et al. 2004): shot noise in the aperture, shot noise in the mean background flux per pixel integrated over the aperture, and the standard deviation of the sky annulus pixels to account for the influence of non-uniform nebulous background. Such error values proved to be a slight underestimate of the true uncertainty, and do not take into account a "floor" in the errors from calibration errors, primarily the intrapixel gain variation.

We wanted to find a value for this floor that worked for all clusters (since all of the mapping observations discussed here were obtained in the same way). We also wanted to be conservative in that we wanted the set of objects we selected as variable at the end of this process to be reliably variable (with very few, if any, non-variables selected as variable), even if it meant that our sample would not be complete in that some legitimately (lower level) variable objects would not be selected. In other words, we wished to err slightly on the side of overestimating errors. We investigated the rms error for each light curve for each object versus magnitude and versus the median photometric uncertainty for all sources. In MC11, we estimated the true error as a function of magnitude to be single-valued as a function of magnitude bin. In general, this means that the errors will be correct on average, but there will be specific objects (e.g., those in high background regions) for which they are too low.

Here, we aimed to improve on this estimate by finding a value for the noise floor that can be added in quadrature to the error estimates for each point, preserving the errors obtained for each point, localized in time and space. For each channel, for every object in each cluster having at least 20 epochs in its light curve so as to ensure well-determined values, we plotted the measured rms (σ) of the light curve against its median photometric uncertainty; see Figure 15. Orion alone represents about one-quarter of the light curves shown in Figure 15, so it is useful in some contexts to separate Orion from the analysis of the rest of the ensemble. Data for an individual smaller-field cluster are thus shown in the figure for reference.

The largest-amplitude sources in Figure 15 are true variables. To determine the noise floor, we are interested in how the observed rms for non-variables over many epochs compares to the formal errors. The bulk of the distribution of points in Figure 15, e.g., the non-variable and less noisy sources in the ensemble of points, fall in the same region in each plot; the distributions run smoothly down to an asymptotic value. However, there are three components contributing to the overall scatter of the regions with more objects seen in Figure 15, where they are not necessarily seen as distinct features. As we move to brighter objects (generally those with smaller median photometric uncertainty), we have a statistically larger chance that the objects will be cluster members, e.g., legitimately young and therefore expected to have a greater likelihood of being highly variable. There is a transition between the short and long HDR frames (at ∼10th mag) which affects the measured error for sources that are faint in the short frames but not in the long frames. There are variations in depth (∼16th–17th mag) reached across clusters due to variable and nebulous backgrounds (the variations in depth can also be seen in the faint limits of Figures 19 and 20, discussed below) such that a given median photometric uncertainty need not necessarily refer to stars of roughly the same brightness among the clusters. Figure 15 includes an individual cluster (Ceph C) to give an indication of this cluster-dependent variation in the floor; it can be seen that fewer points fall below the red line on the left side of each panel in Ceph C than they do in the Orion plots on the bottom row.

Close inspection of Figure 15 reveals that the asymptotic noise floor value is on average slightly higher for the 3.6 μm channel and slightly lower for the 4.5 μm channel. This is consistent with expectations that the dominant source of error in the asymptotic noise floor is from residual intrapixel gain variations, since the intrapixel gain effect is lower in I2 than in I1. There are also more points with significantly larger σ in the 4.5 μm channel. The 4.5 μm channel has more contributed noise in any given light curve—shocks and scattered light, which are very common in our target regions, contribute to the overall scatter in σ at 4.5 μm. There is also some lower-level contribution from polycyclic aromatic hydrocarbon (PAH) features (also common in these regions) in 3.6 μm (see, e.g., Flagey et al. 2006). These contributions to both channels mean that there is some cluster-to-cluster variation in the uncertainties associated with individual sources at a given magnitude.

We determined that a systematic error of 0.01 mag in I1 and 0.007 mag in I2 (both indicated in Figure 15) provides the best fit to the photometric rms seen in each channel in that it is a conservative representation of the distribution as the distributions approach the asymptote. These values are also quite comparable to the magnitude of the intrapixel gain effect in each of these channels, and it is not expected that the errors in IRAC photometry from mapping observations would be less than 0.1%. Very similar values are obtained via a slightly different approach for the CSI 2264 mapping in Cody et al. (2014), where the populations of members and non-members are better understood. We added our empirically derived systematic values in quadrature to each of the individual errors obtained by our pipeline.

Figure 16 plots the light curve rms (σ) derived from the corrected light curves against the mean magnitude for that light curve. Here, again, one can see more, and more significant, variation among the brighter sources. Relatively few sources have much larger uncertainties because they are close to bright neighbors, chip edges, or affected by other instrumental artifacts; many of the sources with large rms are legitimately variable. For sources brighter than 13th mag, the total uncertainty is dominated by the error floor. The overall noise increases significantly fainter than about 14th magnitude, and the rapid falloff of the number of objects detected in the survey beyond 16th magnitude is immediately apparent (see also Figures 19 and 20 below). There are small variations between clusters, which are caused by different levels of diffuse emission and a different degree of instrumental artifacts from, e.g., residual pull-down. A sample individual cluster is included in Figure 16 as an indication of what is seen in the individual smaller-field clusters.

Zoom In Zoom Out Reset image size

Early in the Spitzer mission, each of our 12 regions was observed, often as part of the guaranteed time observations or the original Cores-to-Disks (c2d) Legacy program (Evans et al. 2003, 2009b); these observations are summarized in Table 4. In most cases, the data were obtained at two epochs. For clusters near the ecliptic plane, asteroids that happen to be passing through the region at the time of observation appear as long wavelength sources that can resemble embedded low-mass YSOs. Thus, the two epochs were often planned to be separated by a time of the order of hours (e.g., the earliest epochs of Serpens, much of L1688) to allow the moving solar system objects enough time to move, and the corresponding pixels to be rejected as outliers when the individual frames were combined. For some observations (e.g., Orion), to facilitate bright source artifact mitigation, the epochs were separated by ∼6 months such that the orientation of the arrays in the later observation is 180° from that in the early observation (see Section 2.3 above). For the most robust detections, the cryogenic data can be combined into a single early epoch, or one can compare the two epochs to constrain variability on timescales of months (e.g., Megeath et al. 2012 for Orion at IRAC bands) or hours (e.g., Rebull et al. 2007 for Perseus at 24 μm).

To retain information from the individual epochs, the cryogenic data for all of our clusters were reduced identically to the YSOVAR monitoring data, except using cryogenic calibrations (e.g., the zero point and threshold between the short and long HDR frames). Source detection was performed independently in each band for each epoch. However, for reasons we now describe, measurements derived from individual epochs from the cryogenic era cannot be used generally over the entire region, and the identical approach adopted for the rest of the YSOVAR data cannot be used generally for the cyrogenic data.

For the YSOVAR data, each epoch covered the region of interest with sufficient redundancy, so that data reduction can be performed on a per-epoch basis as described in Section 2.5 above. For some of the cryogenic era data, the observations were designed to have sufficient redundancy in the region of interest only once all of the cryogenic era observations were combined. Specifically, for some of those clusters that were observed in two epochs, for moving object identification, the observations from a single epoch were not designed to robustly measure objects at that epoch; it was envisioned that the two epochs would be combined to create a net mosaic of the inertial targets. For example, one epoch might have two frames with another two frames at a later epoch. Combining those after the fact provides four frames of dithered observations to remove artifacts and moving objects. However, a single epoch, having only two frames, does not necessarily have enough data for robust photometry. A complication comes from the combination of individual frames to compose a map at a given epoch. The maps at any one epoch are constructed from multiple pointings, with overlap between each pointing. Therefore, in those thin regions between adjacent pointings, there are more than two frames at a given epoch. Because our pipeline requires at least three frames per position, it only derives (and archives) high-quality photometry at either of those two epochs for that thin region between adjacent pointings in the map, not the entire map at that epoch. Therefore, if one downloads the earliest epoch data from our archive for those clusters, one obtains a "grid-shaped web" of viable data, not complete coverage of the region. We have retained all of these individual epoch data in the data we deliver to IRSA because it may be useful to future researchers to have individual measurements of sources serendipitously falling in those regions. However, for users not aware of this "feature," it will seem strange to only have cryogenic data over a web of coverage in those early epochs.

In the context of the present paper, therefore, all of the earliest cryogenic-era Spitzer observations were combined into one mosaic per epoch per channel, which was then run through Cluster Grinder. The time assigned to this early epoch of observation is the average time of all of the observations that went into the mosaic, summarized in Table 4.

As previously stated, IC 1396A and Serpens Main were explicitly monitored in the cryogenic era with Spitzer, and these data were also reprocessed, retaining individual epochs of observation.

The data were bandmerged across Spitzer bands by position, within a search radius of 1''. (Recall that the IRAC photometry apertures we used were 24.)

Gutermuth et al. (2008b, 2009, 2010) present methodology for identifying YSOs from the cryogenic catalog. The details of the selection process appear in those papers, but in summary, multiple cuts in multiple color–color and color–magnitude diagrams are used to identify YSO candidates, as distinct from, e.g., extragalactic and nebular contamination. This color selection is part of Cluster Grinder, and thus we have this classification, based on the cryogenic data, for all of our clusters; we used it in Section 2.1.2 and in part of Table 1. We have adopted this YSO selection mechanism as part of one of the primary YSOVAR sample definitions; see Section 3.

For completeness, we note here that, in addition to the new IRAC data, we often also obtained contemporaneous optical and NIR observations from the ground using a variety of telescopes. The details of each telescope/camera and set of accompanying observations is beyond the scope of this paper, and is (or will be) provided in each paper discussing the results from each cluster. For example, MC11 discussed observations obtained from four other ground-based telescopes obtained in 2009–2010, contemporaneous with the first epoch of Orion observations; Cody et al. (2014) discussed observations obtained contemporaneously with the CSI 2264 campaign. Because these ancillary monitoring data are different for each cluster, we omit them here for clarity.

We also included 2MASS JHK_s data for all of our clusters. For two clusters (NGC 1333, L1688), deeper than survey 2MASS data (the 6× data) are available from the 2MASS archive, and are included where available. Such detections were also bandmerged to the rest of the catalog by position, within a search radius of 1''.

For three of our clusters (L1688, Mon R2, Serpens Main), deeper NIR data from the UKIRT Infrared Deep Sky Survey (UKIDSS; Lawrence et al. 2007) broadband data were publicly available, but not necessarily at all or even most of the UKIDSS bands (Z, Y, J, H, K); in all three cases, there were at least K-band data, which is important for our calculation of SED class (see Appendix B). Where available, these data are also bandmerged to the rest of the catalog by position, within a search radius of 1''. Details of which bands are available and to what depth will be included in the individual papers associated with the relevant clusters.

For two of our clusters, AFGL 490 and Ceph C, data in two Spitzer bands are available from one of the warm portions of the Spitzer program called the Galactic Legacy Infrared Mid-Plane Survey Extraordinaire (Benjamin et al. 2003). Where available, those single-epoch, shallow data were also bandmerged to the rest of the catalog by position, within a search radius of 1''. Since those measurements are independent measures of these objects at [3.6] and [4.5], these measurements were retained as distinct points from the cryo [3.6] and [4.5], or the means of our light curves.

For three other clusters, NGC 1333, L1688, and Serpens Main, data are available from the c2d (Evans et al. 2003, 2009b) program data deliveries. The data used for these deliveries are typically the same BCDs as were used in the cryogenic data (Section 2.7). As such, then, they are not independent measurements, and these data were only used to supplement our cryogenic-era catalog if a band was missing, e.g., because G09 had identified it as having insufficient S/N in that reduction.

Naturally, each cluster has different amounts of additional data in the literature, and the details of exactly which data are included (and how it was merged to the rest of the catalog) will appear in the corresponding YSOVAR cluster paper.

X-ray observations are an excellent complement to mid- and near-IR observations in regions of star formation. Due to the strong correlation between X-ray luminosity and age, X-ray emission is particularly effective in identifying YSOs that are free of IR excess (Class IIIs), such that X-ray surveys provide samples of young stars unbiased by their IR characteristics.

X-ray observations can penetrate up to A_V ∼ 500 mag into a star forming cloud with very deep integrations (Grosso et al. 2005). However, even with shallow integrations, one can reach as deep or deeper in the X-rays than many NIR surveys in the JHK bands. The Chandra X-ray Observatory, with its high angular resolution mirrors and low-noise detectors, is particularly effective in resolving crowded fields down to 05 scales. Furthermore, in X-rays, OB stars are often not much brighter than pre-main-sequence stars, so close companions, even when associated with OB stars, can be identified (Stelzer et al. 2005). Chandra data (or any other X-ray data, e.g., from the European Space Agency's XMM-Newton observatory) and Spitzer data combined have been shown to yield a more complete survey of members; there are many examples of this in the literature. For some of the clusters in our target list, an incomplete listing of such work could include Winston et al. (2007, 2010) for NGC 1333 and Serpens Main, respectively, Günther et al. (2012) for IRAS 20050, Getman et al. (2006) for IC 1396A, or Imanishi et al. (2001) for L1688.

Identification of additional cluster members is the primary purpose for which we include X-ray data, where available, for the YSOVAR clusters. We can then compare variability characteristics for the X-ray-detected sample with that from, e.g., the IR-selected sample. With the exception of AFGL 490,³² all of the YSOVAR clusters have been observed by Chandra using the Advanced CCD Imaging Spectrometer for wide-field imaging (ACIS-I), and nearly all of them have already been published. The date, duration, and aimpoints of these observations are listed in Table 5, along with citations to the literature where possible; footprints of the observations are included in Figures 2–13.

In order to create a more unified set of detection criteria, and in order to reach fainter sources, X-ray data for the nine smaller-field clusters with X-ray data were reprocessed to achieve internal consistency and maximize source detection. We used Chandra pipeline DS 8.4.5. Detailed X-ray data and analysis have already been published for the very deep observation of the ONC (Feigelson et al. 2005; Getman et al. 2005), and for two shallower fields north and south of the ONC (Ramirez et al. 2004b); similarly, for NGC 2264, Flaccomio et al. (2006) and Ramirez et al. (2004a) have published deep Chandra data (our region here is primarily covered by Flaccomio et al. 2006). Because those two regions include X-ray data substantially deeper than those for the rest of the clusters, we have not reprocessed these X-ray data here, but instead taken those source lists from the literature. (A primary reason we reprocessed the nine smaller-field clusters is to reach fainter sources, which the deeper integrations in Orion and NGC 2264 already accomplish.) Additional X-ray data obtained contemporaneously with Spitzer monitoring will be discussed in the appropriate cluster-specific papers.

The region with ACIS-I data is usually smaller than the full region observed with Spitzer during the cryogenic era, but often covers the region monitored for YSOVAR; see Figures 2–13. Aside from AFGL 490 (which has no X-ray data), the region with the least spatial Chandra fractional coverage is Orion; the Orion region with IR light curves is by far the largest of our sample. The other 10 clusters have complete or nearly complete X-ray coverage of the region with two-band IRAC light curves; X-ray coverage of the regions with light curves in only one-band varies. The Chandra sensitivity varies across the FOV, with higher sensitivity in the center of the pointing, decreasing toward the edges.

Source detection from the archival Chandra data is performed with the wavdetect algorithm in CIAO (Chandra Interactive Analysis of Observations; Fruscione et al. 2006), versions 4.5 and 4.6. The details for the detection process used here are identical to those used in Günther et al. (2012). Our chief goal is identification of all X-ray sources, since even a weak detection—provided it is coincident with an IRAC source—has a high likelihood of being a legitimate source, though this does not necessarily prove cluster membership. Table 6 lists the number of X-ray sources detected above the 2σ significance level in the entire ACIS-I field (a single ACIS-I FOV is ∼16' × ∼ 16'). For comparison, Table 6 also indicates the number of detections reported in the literature, often for exactly the same data. We recover essentially all the previously reported sources, and add a significant number of weaker ones due to the lower threshold employed here. Note that without the additional criterion imposed here of requiring a match between an X-ray source and an IRAC source, the 2σ significance level used here would be too low and would result in a high number (about one-third) of spurious sources. The requirement of a positional match with an IRAC source allows us to use such a low threshold. We found that below 2σ, we do not find matches between X-ray and IR sources in excess of matches expected due to random chance.

To determine if a given X-ray source is coincident with an IRAC source, we matched X-ray sources to IRAC source lists generated for Spitzer cryogenic-era observations, as described in Section 2.7 above. Due to the highly spatially dependent Chandra point-spread function, three matching radii were used. For sources within 3' of the aimpoint, the matching radius was 1''. For sources more than 6' off-axis from the aimpoint, the matching radius was 2''. In between, the matching radius was 15. We note that IRAS 20050+2720 and IC 1396A are both exceptions since the fields were each observed multiple times, at different roll angles. Since IRAS 20050+2720 covered a wide range of rotation angles, in that case, the source positions are composites of different point spread functions; matches were made more carefully in these observations. For IC 1396A, the range of roll angles was smaller, so the same approach was used as for the rest of the clusters. For each cluster, in Table 6, we list the total number of X-ray sources with a match to an IRAC source within the appropriate radius. We note here for completeness that there are some very bright X-ray sources without IR counterparts; for our purposes within YSOVAR, we exclude these sources.

If an X-ray source is bright enough for spectral fitting, one can determine four key characteristics: X-ray temperature, gas absorption, mean flux (F_x), and variability. We fit a single Astrophysical Plasma Emission Code (APEC) thermal emission model (Smith et al. 2001) for each source detected above 30 net counts, using C-statistics and leaving the absorbing column density (N_H), the temperature (T), and the volume emission measure as free parameters to determine the first three of the above characteristics. The metallicity is fixed at 30% of the solar abundances (Anders & Grevesse 1989), in keeping with typical values from coronal emission from late-type stars. Details of the extraction and fitting processes are similar to those discussed by Winston et al. (2010), which follows the procedure for automated processing laid out for the ANCHORS (AN archive of CHandra Observations of Regions of Star formation) pipeline (Spitzbart et al. 2005). The key point is that all fits are done assuming that the source is a star with a thermal one-temperature spectrum. We calculate luminosities (L_x) for each source using the fluxes (F_x) from 0.3 keV to 8.0 keV and line-of-sight absorptions, assuming the distance to each cluster given in Table 1. Detailed flux errors due to the fit were not calculated separately for each source, since systematics are likely to dominate; instead, global flux errors were determined using the CIAO tool dmextract. This process calculates a simple error estimate based on photon statistics and the mean value of the exposure map in the source region. These errors are about 4% at 2000 counts, about 35% at 100 counts, and the errors reach 100% below about 50 counts. The error budget is dominated by photon counts and uncertainty in N_H. There is no evidence that such errors are markedly biased by the C-statistic.

Since errors on L_x can be dominated by systematic effects, we have reprocessed all sources (in those nine smaller-field clusters with Chandra data), even those with previously published fluxes, to ensure uniform spectral fitting methodology. For faint sources with fewer than 30 counts, no fit can be performed. In this case, we determine a median photon energy, which, when combined with the count rate, leads to an approximate flux determination. All spectral properties presented here are effectively time averaged over all observations.

To determine the fourth of the characteristics above (variability), we tested the light curves for variability using the Gregory–Loredo method (GL-vary; Gregory & Loredo 1992). This method uses maximum-likelihood statistics and evaluates a large number of possible break points from the prediction of constancy. It assigns an index to each light curve—the higher the value of the index, the greater the variability. Index values greater than 7 indicate >99% variability probability. Values of the GL-vary index >9 usually indicate flares. GL-vary is not reliable below about 30 raw counts, the same limit we used for performing spectral fits. The value of this index will be provided where relevant on a source-by-source basis in the individual cluster papers.

There are five broad classes of X-ray sources in the field.

• 1.  
  First, background active galactic nuclei (AGNs)—these will be numerous; up to 50 of these per 16' × 16' Chandra field are expected, depending on the depth of exposure (Getman et al. 2006). However, they tend to be faint in both X-rays and mid-IR, or are not matched in the IRAC bands at all.
• 2.  
  Second, background starburst galaxies—while much less common than AGNs, they are brighter than AGNs, and have colors similar to Class II YSOs. They tend to be fainter in the IR than typical Class II YSOs, and may be tentatively identified via the faintness of the IR counterpart.
• 3.  
  Third, compact objects such as white dwarfs—these tend to be very faint and usually undetected in the IRAC bands.
• 4.  
  Fourth, YSOs—those with disks have already been identified via their IR excesses. We identify the probable disk-free objects that are detected in X-rays by their star-like IR colors and magnitudes consistent with membership. In addition to matching an IRAC source, the star-like colors are required to ensure that any newly revealed sources are probable Class III objects and not distant starburst galaxies.
• 5.  
  Fifth and finally, active late-type field stars—both foreground and background stars can appear with X-ray fluxes comparable to our targets. Based on a study of IC 1396, Getman et al. (2006) estimate that there could be ≲ 10 of these per 16' × 16' Chandra field. Since the contaminants are a mixture of foreground and background objects of comparable fluxes to our targets (with less foreground and more background contamination likely for closer clusters), without optical spectra, they are very hard to discern from Class III objects. It is necessary, then, that these remain in our sample of candidate members selected via X-rays and represent a source of contamination.

We note that Getman et al. (2012) estimate (again for IC 1396, a region in the Galactic plane) a 22% probability that any X-ray source with any IRAC counterpart is not a member of the cluster. That rate drops to about 10% in IC 1396 if one eliminates sources without 2MASS JHK_s detections. Because the contamination rate is affected by absorption, exposure time, and the depth to which one extracts sources, it may be different for the other clusters.

We can improve our inventory of YSOs in these clusters by identifying objects with X-ray detections, IRAC counterparts, and SEDs that are consistent with those of stars. We have adopted this YSO selection mechanism as the other main component of the primary sample definitions; see Section 3. Note that we define SEDs consistent with stars to be those with a fitted SED Class III (see Appendix B), but that there is room to create an augmented membership list (Section 3.3) to include objects that have an X-ray detection, an IRAC counterpart, were not identified as a YSO from the IR alone, but that have an SED consistent with a YSO. With regard to foreground or background stellar contamination, because any appropriate brightness cutoff is a function of cluster distance (and A_V), we defer any detailed exclusion of likely foreground or background stars from the member sample to the individual cluster papers.

Each cluster has an IR-selected sample of member candidates defined by the Gutermuth et al. (2008b, 2010) and G09 selection algorithm, run on the cryo-era catalogs created anew as per the methodology above (Section 2.7). A detailed description of the YSO candidate color selection algorithm can be found in Appendix A of G09. A sample of YSOs selected this way is thought to be a statistically well-defined sample (see statistical discussions in Gutermuth et al. 2008b and G09), comprised almost entirely of members, though some contamination from background galaxies or asymptotic giant branch stars is always possible. Very few of these IR-selected members have spectroscopic follow-up (or spectra in the literature pre-dating the Spitzer observations), and, as such, many should technically be thought of as YSO candidates, though we include them all in the set of members.

All clusters except AFGL 490 have X-ray data, and thus also an X-ray-selected sample of YSOs. (As with the IR-selected sample, few of these have spectra, so technically they are YSO candidates.) Since we should have identified most of the YSOs with disks (at least disks detectable at 3.6–24 μm) in the IR selection process, we use the X-ray data to identify additional young stars without disks (e.g., Class IIIs). Thus, we add to the set of IR-selected members the X-ray-selected sample defined by the algorithm described above in Section 2.10, which can be summarized as objects having an X-ray detection above a 2σ significance threshold, having a match to an IRAC object in the cryo-era Spitzer catalog, and having an SED shape consistent with it being a disk-free YSO or a star (e.g., Class III; see Appendix B.2). This sample is specifically designed to find and add to our set of members those members without disks. However, it should be noted that (1) the Galactic contamination rate is likely to be higher in this X-ray-selected sample than in the IR-selected sample; and (2) specifically because of the contamination rates, members with disks are identified as members from the IR excess, not the X-ray flux, though, of course, members with disks can also have measured X-ray fluxes in our database.

We have thus defined our "standard set of members" to be the union of all IR-selected members with disks and X-ray-selected members without disks. There are provisions for adding additional objects; see Section 3.3 below. Note that this definition can be applied independently of whether or not there is a light curve, but of course in the context of this discussion of YSOVAR data we require a light curve. Note also that the IR selection requires four bands of IRAC, and thus both very faint and very bright previously identified members may be omitted from the standard set; objects such as these known to be members via some other approach in the literature may be added in the augmented sample (Section 3.3). Finally, note that because the data that go into our selection of cluster members are of various depths, and because the clusters are at a variety of distances, the effective mass limit reached by each set of standard members varies from cluster to cluster.

All of the original YSOVAR light curves were obtained with very similar HDR mapping observations (see Section 2.2) in a fast cadence (see Section 2.4), though the length varies, and some clusters have additional slow cadence observations and/or staring observations. The time sampling and total length of light curves obtained within a single cluster field can also vary as the FOV changes with time (Section 2.3). We have attempted to remove all instrumental effects from the input data, and only retained photometry where there were valid measurements on at least three BCD frames (Section 2.5).

We now define the "standard set of data for statistics" as follows. Since the fast cadence is the most common (and most similar) among the YSOVAR clusters, we used only these mapping fast cadence data, for those light curves that have at least five viable epochs (with each epoch obtained from at least three BCDs per epoch that are not obviously compromised by instrumental effects or cosmic rays), to calculate statistical quantities such as mean, median, etc., as well as Stetson index and χ² (discussed below; see Section 5). We have defined statistical values calculated on the fast cadence data as the "standard set of statistical values" for each cluster and employ them to identify variables and compare values across clusters. Finally, for stars fainter than [3.6]–[4.5] ∼ 16, noise tends to dominate the light curves (Sections 4.2 and 5). We have retained these faint sources, but objects this faint are considered individually where relevant.

Note that the standard set for statistics is thus defined as all light curves with at least five points, just in the fast cadence. This is independent of whether or not the target is identified as a member.

Elsewhere in this paper, we refer to "all objects with a light curve"—this means anything with a light curve in the standard set for statistics. (Essentially no YSOVAR-classic sources have only points outside of the fast cadence.)

When available, additional epochs of data can be included for additional calculations on a per-cluster basis in the corresponding papers, and will clearly be indicated as such where relevant.

We identified members above using IR and X-ray data, which implicitly relies on the shape of the SED between 2 and 24 μm. That sample is still the best set of members that we will use to compare across clusters, because that set of members is defined as similarly as possible across all clusters. However, additional young objects may be identified in the literature, and additional members may be suggested based on our own data. The "augmented sample of members" is where these additional likely members can be included.

Some clusters (e.g., NGC 1333, L1688) have considerable literature discussion of members, and others (e.g., AFGL 490, Mon R2) have far less. We therefore cannot rely exclusively on the literature to select members, but neither should we ignore members identified in the literature and not selected above. Thus, each cluster paper may include in the augmented sample the literature-identified sources that are not already found using our IR or X-ray methods above.

We can use our own data to identify new cluster members. While only 1%–2% of the field population may be variable, YSO variability is the rule, not the exception. It is therefore possible to identify new cluster members from light curve properties alone; one could identify all variables as new members, or one could take only those with certain properties such as amplitude above a threshold. We could identify cluster members from either the standard set for statistics (the fast cadence), or, for those clusters with longer cadences, from those additional data.

The set of statistically selected variables (see Section 5) are those with Stetson index greater than 0.9, and/or with χ² greater than 5, and/or with a significant period, calculated over only the standard set for statistics (the fast cadence). Often, these variables should also be identified as cluster members, but these individual objects will be discussed on a per-cluster basis (because, for example, they can be background eclipsing binaries; Morales-Calderón et al. 2012). Variables identified using data beyond the standard set for statistics (data beyond the fast cadence data) will also be included on a per-cluster basis, and will be discussed in the cluster papers. Those newly identified cluster members may also be included in the augmented sample of members.

While variability-identified objects will make an important contribution to our understanding of the complete membership of each cluster, they should not be used in calculations of, e.g., variability fractions; that sample should be selected on the basis of a parameter distinct from variability, such as disk excess or X-ray emission. This is why the new candidate members we identify from our data are included in the augmented sample of members and not in the standard set of members.

This augmented set of members is only used (where it is clearly identified) in the individual cluster papers, not in the remainder of this paper.

In order to compare results among clusters, it is useful to be able to place clusters in some sort of relative order that could, in the most useful (though hypothetical) case, be tied to age. There are various ways of parameterizing the evolutionary state of these clusters, and we considered several, all aimed at capturing the relative numbers of sources in various SED class bins (see Appendix B for a brief definition of SED classes and our placement of sources therein). Such ratios in some sense capture the relative "degree of embeddedness" of sources in these regions, perhaps with an ultimate (though undefined here) link to cluster ages. Formally, we are binning by SED slope, so the parameterizations are, strictly speaking, relative fractions of sources with a given SED slope, indicative of the amount of circumstellar material, e.g., how self-embedded a given source may be. We interpret clusters with more sources with large SED slopes to, on average, contain more sources that are more embedded.

G09 chose to use the ratio of Class II to Class I sources. These ratios were all obtained internally consistently, e.g., sources selected and categorized according to the same series of IR color cuts and data reduction. Because G09 was working with Spitzer data over a relatively large region for each cluster, this Class II to Class I ratio could be calculated for the whole region and for subsections of the region. The Class II to Class I ratios as calculated for the cores of these clusters (e.g., Table 6 in G09) appear in our Table 1 and in the discussion in Section 2.1.2.

In YSOVAR, we have clusters not included in G09, and moreover, even for the clusters included in G09, we generally have light curves for only a small subset of the sources G09 considered because we observed a smaller region. To calculate a Class II to Class I ratio for the portion of each cluster sampled by the YSOVAR monitoring, we performed the same calculation for objects in the standard set for statistics (e.g., having YSOVAR light curves) by reducing the cryogenic data in the same way and performing the same series of G09 color cuts and classification; see Section 2.7. Then, we recalculated the ratio of Class II to Class I sources specific to the YSOVAR data using the classes assigned via the G09 algorithm only for those sources with light curves. Those Class II to Class I ratios also appear in Table 1 and in the discussion in Section 2.1.2. The values are provided in the present work in part as a link to the G09 analysis; note that they are calculated in the same way as G09, e.g., with the IR-selected sources alone, not on the standard set of members per se.

The G09 parameterization is based only on IR-selected sources. For most of our clusters, we have X-ray data as well (see Section 2.10), so we at least have some information on the Class III population. It is, however, true that we do not always have complete Chandra coverage of our fields (further discussion of coverage appears in Section 2.10 and Figures 1–13), and even for clusters where we have Chandra coverage, the Chandra sensitivity is a strong function of location on the array. Nonetheless, we would like to include the information we have, and simply using the Class II to Class I ratio does not incorporate information about the Class III population.

We explored several alternative parameterizations of the relative fractions of embedded sources, all of which involved various ratios of classes (or groups of classes) to the total or other classes (or groups of classes). Histograms of the relative fractions of the SED classes for the standard set of members (Section 3.1) for objects with light curves in the standard set for statistics (Section 3.2) in each cluster appear in Figure 17. AFGL 490 can be seen to be deficient in a complete sample of Class III objects because it has no X-ray data; the Class III objects identified here have small IR excesses (or sufficient reddening at 2 μm) and thus SED slopes consistent with Class III. By many metrics, Serpens South has the highest fraction of sources with positive SED slopes, which we take to be most embedded sources. Similarly, both Orion and IC 1396A have the highest fraction of sources with negative SED slopes, which we take to be less embedded sources.

Zoom In Zoom Out Reset image size

For further analysis here, we have settled on the fraction of Class I sources to the total number of members, for objects with light curves (objects in the standard set for statistics). These ratios are also included in Table 1. They include (as part of the total number of members) objects selected via X-rays. However, there are still the fundamental, systematic uncertainties inherent in the classification approach, in the selection of members without spectroscopic follow-up, in the completeness of the surveys involved (in both area and depth), and the requirement that there be a light curve (e.g., bright enough in the IRAC channels) as well as in the relative paucity of Class I objects overall. Additionally, for AFGL 490, there are no X-rays to be used, so the ratio is calculated with solely the IR-identified members. However, the Class III objects appear only in the denominator, combined with all the other classes. We also note that our inventory of Class I sources must be incomplete, since we lack the long wavelength coverage that would be needed to find the most embedded sources (e.g., Stutz et al. 2013), but such extremely embedded sources will also not have a YSOVAR light curve.

This parameterization using the Class I/total ratio should be related to the G09 Class II/Class I ratio determined for the objects having light curves. Figure 18 plots these parameterizations against each other. The error bars are derived assuming uncorrelated Poisson statistics, because it is difficult to quantify the additional systematic errors described above. The best-fit slope of the line fit to this relation (taking into account errors in both directions on each point) is −0.024 ± 0.005. The correlation coefficient, Pearson's r, calculated for these two parameters is −0.87, and a probability that the parameterizations are not correlated of only 0.04%. We assume based on the statistics and our underlying physical intuition that these values are indeed correlated.

Zoom In Zoom Out Reset image size

Several individual points in Figure 18 merit additional discussion. Despite having no X-ray data, AFGL 490 is consistent with the trend shown in Figure 18. (If one assumes that there might be about as many Class IIIs as Class IIs in this region, then the point could move down to about 0.1–0.15, which would still be broadly consistent with the trend.) The NGC 2264 point is below the trend. The X-ray data obtained for NGC 2264 is deeper than the X-ray data for the other clusters, except for the central Orion region. This results in more of the fainter sources being included in the total number of YSOs, and pushes the Class I/total ratio toward lower numbers, as seen. For both parameterizations, Serpens South is selected as the cluster with the most embedded sources, as expected from Figure 17. It is well above the fitted line in Figure 18; perhaps an exponential decay rather than a simple line would be a better fit to use, but, in the absence of additional very embedded clusters to constrain the most embedded end of the distribution (or, indeed, spectroscopic vetting of the members and the reduction of other such uncertainties), a line is the simplest fit to use. IC 1396A, based on Figure 17, should be one of the clusters with the fewest embedded sources. It is identified as the least embedded using the Class I/total ratio; there is a large uncertainty on the Class II/Class I ratio, and it is consistent with being the least embedded within 1σ. The Class II/Class I ratio formally identifies Orion as the least embedded. However, aside from AFGL 490 where there are no X-ray data, Orion is the cluster with the poorest match between the YSOVAR-monitored region and the existing X-ray data coverage. Orion has a very deep X-ray pointing, but only in the central ONC region (Feigelson et al. 2005; Getman et al. 2005). There are two shallower pointings that contribute X-ray data (see Section 2.10 and Figure 3), but the region of sky in Orion for which we have IR light curves has the least fractional coverage in X-rays of all our clusters (aside from AFGL 490). To investigate the degree to which the uneven X-ray coverage affects the placement of Orion in this diagram, Figure 18 also includes points for Orion when broken into "North" (declination > −05:05:25), "South" (declination < −05:33:15°), and "ONC" (between those two limits) fields. There is scatter, clearly, in these points, but when those points are used instead of the single Orion point, the fit is functionally indistinguishable from that using the single Orion point.

We use the Class I/total parameterization in subsequent discussions in this paper, most notably Section 6.

Figures 19 and 20 show the distribution of J, [3.6], and [4.5], in units of magnitude and the percent of the sample, for the standard set for statistics (all objects with at least five points in the YSOVAR fast-cadence light curve, including cluster members and field objects). The total number of objects portrayed in the figures ranges from ∼100 (J: L1688 or Serpens South) to ∼7200 ([3.6] and [4.5]: Orion). The J data are generally from 2MASS; notably, NGC 1333 is deeper than the other clusters in J because additional data from the 6 × 2MASS survey have been included (Section 2.9). Generally, more objects are available at [3.6] or [4.5] than at J, largely due to the greater effective depth of Spitzer. In several cases (most notably L1688, Serpens South, and NGC 2264), a relatively small fraction of the objects with Spitzer light curves have J counterparts, since these clusters are, on average, generally more embedded than the others. For most of the clusters, for most of the objects, JHK_s data are not available for objects with [3.6] or [4.5] fainter than about 15th mag.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The Orion maps extend out beyond the edges of the cluster, and include a higher proportion of field stars and other contaminants than do the other smaller-field clusters. This can be seen in the structure in the Orion [3.6] and [4.5] histogram, which is double-peaked; the brighter peak is likely dominated by the cluster members, and the fainter peak is likely dominated by contaminants. For similar reasons, if the UKIDSS data (Section 2.9) are included, the Serpens Main J histogram extends to J ∼ 20 and is also double-peaked, but the L1688 J histogram is not so obviously double-peaked, likely due to the higher obscuration levels of the background population.

Mon R2 seems to be different from the other clusters in that the fainter end of the [3.6] and [4.5] histograms are considerably flatter. This is most likely a symptom of the difficulty of obtaining Spitzer light curves for faint sources in the presence of high and spatially variable background; there are some very IR-bright sources in the IRAC FOV (see Figure 4), and the scattered light is substantial, coupled with intrinsically bright outflows and PAH features. No UKIDSS J mag in this region were in the public archive at the time we checked (in 2013 September).

As can be seen from the turnover at about 16th mag in the [3.6] and [4.5] histograms in Figures 19 and 20, our data do not extend much fainter than [3.6]–[4.5] ∼ 16 mag. By inspection, all of these faint objects appear to be legitimate point sources on the images. For stars fainter than this, noise tends to dominate the light curves (see Section 5). It is hard to completely reject these fainter sources—we could exclude those whose cryogenic-era [3.6] > 16, though the YSOVAR epochs could vary above and below that boundary; or we could discard those where the mean during the YSOVAR campaign is >16, but there are some objects for which the mean [3.6] <16 but the mean [4.5] >16. As noted in Section 3, we have retained these faint sources, but, in general, the large uncertainties associated with their Spitzer photometry preclude strong statements about their variability; objects this faint are considered individually where relevant.

Similar histograms for only the standard set of members (identified through IR excess and X-ray emission; see Section 3) are considerably less populated, and brighter; the peaks are around [3.6]–[4.5] ∼ 12 mag. This is consistent with the location of the brighter peak in Orion, and several of the tails of the distributions seen in Figures 19 and 20.

Figure 21 contains histograms of the log of the X-ray luminosities (log L_x, where L_x is in erg s⁻¹) for those objects with light curves (standard set for statistics) and bright enough in flux (F_x) to have a calculated L_x (Section 2.10). The Orion and NGC 2264 histograms reach fainter values in F_x than those of the other clusters because those integrations were considerably deeper. Moreover, essentially the entire NGC 2264 field considered here has X-ray data, whereas the fractional X-ray coverage of the Orion field is relatively low compared to NGC 2264 or the other clusters here (aside from AFGL 490, where there is no X-ray data). Because this figure has incorporated distance (distances are listed in Table 1) to the clusters in the calculation of L_x, the two closest clusters (NGC 1333 and L1688) have histograms reaching the faintest L_x.

Zoom In Zoom Out Reset image size

The first way we identify variables is the Stetson index (Stetson 1996), which quantifies correlation of variability in two (or more) bands. The Stetson variability index is computed for each object as:

$$\begin{equation} S=\frac{\sum _{i=1}^N g_i \times sgn(P_i) \times \sqrt{|P_i|}}{\sum _{i=1}^N g_i}, \end{equation} \tag{ 1 }$$

where N is the number of pairs of observations for a star taken at the same time,³³P_i = δ_j(i)δ_k(i) is the product of the normalized residuals of two observations, and g_i is the weight assigned to each normalized residual. In our case, the weights are all equal to one. The normalized residual (δ) for a given band is computed as:

$$\begin{equation} \delta _i=\sqrt{\frac{N}{N-1}}\frac{{\rm mag}_i-\overline{{\rm mag}}}{\sigma _i}, \end{equation} \tag{ 2 }$$

where N is the number of measurements used to determine the mean magnitude and σ_i is the photometric uncertainty. Objects with larger values of the Stetson index are typically taken to be variable. Since errors are included in the calculation, light curves that are just noisy are not identified as variable. Objects with variability in different bands that is not correlated will not be identified via this method; physically, we expect most YSOs to have similar variations in the two IRAC channels since it is difficult to imagine processes that would make one IRAC channel vary without the other. However, this method will not find variables in cases where one IRAC channel is compromised, e.g., due to instrumental effects, and the other is not.

If there is correlated noise between the two channels used for the Stetson index, especially at the faint end, one would expect the Stetson index to be correlated with source brightness, as seen in, e.g., Plavchan et al. (2008b) or CHS01. Figure 22 shows the distribution of all calculated Stetson indices against mean [3.6] (for the standard statistical sample, e.g., all objects with light curves, over all clusters, fast cadence only). Unlike the analogous figures found in, e.g., Plavchan et al. (2008b) or CHS01, here we have no substantial change in the bulk of the distribution of the Stetson index toward fainter [3.6] magnitudes, so we do not appear to have correlated noise between the two channels. There is an increase in frequency of large Stetson index values for brighter [3.6], but that is likely because brighter objects are more likely to be legitimately young cluster members, and, as such, are more likely to be variable (with variability correlated between the two IRAC channels). Essentially no large values of the Stetson index are identified for objects fainter than [3.6] ∼ 16; the intrinsic error on each point is sufficiently large that these objects do not have a large Stetson index. (This is consistent with the discussion in Section 4.2 regarding the number of objects in our data set falling off rapidly fainter than [3.6]–[4.5] ∼ 16.)

Zoom In Zoom Out Reset image size

The specific location of the cutoff between variable and non-variable can be unique to each data set, as it is affected by the sampling length and rate of the light curves. This is the primary reason behind our decision to calculate statistics over only the fast cadence window. We now discuss how we chose this cutoff value for the Stetson index. The left panel of Figure 23 shows a histogram of the Stetson indices for all objects in the standard set for statistics having at least five points in both I1 and I2. The bulk of the distribution about 0 are the non-variables, and that part of the distribution can be reasonably well-fit by a Gaussian. There are substantial deviations from Gaussianity toward higher values of the Stetson index, as expected for a population of identified variable stars. There is a change in the distribution of the Stetson index above and below 0.9; from where the distribution deviates from a Gaussian to a Stetson index of 0.9, the slope in the middle panel of Figure 23 is ∼ − 1.0, and from a Stetson index of 0.9 to ∼3, the slope is ∼ − 0.3. Based on this, we take 0.9 as the cutoff for variability in our data set. The value of 0.9 corresponds to about 6σ for the Gaussian fit to the distribution.

Zoom In Zoom Out Reset image size

We note that in the analogous histograms for each individual cluster, each typically has a small gap in the Stetson index distribution at ∼0.9. Orion, however, does not, and Orion contributes about half of the ∼11,000 viable two-band light curves for which the Stetson index appears in Figure 23.

To check whether the Stetson index cutoff of 0.9 is sensible, we conducted a series of Monte Carlo tests. For random light curves (with a Gaussian distribution of points) using the same time sampling as the real data, the distribution of Stetson indices is well-described by a Gaussian with a width typically of 0.1–0.2, comparable to the left-hand side of Figure 23. For Orion, where we have a reasonably well-defined set of members and non-members (from MC11; not just disked and non-disked, but confirmed membership lists), we can compare the distributions of Stetson indices for the members and non-members. In Orion, the distribution for non-members is generally fairly well-described by a Gaussian centered on 0, but there is a small "shoulder" asymmetry toward the larger Stetson index (likely legitimate field variables or as of yet unidentified members). The distribution for members, in contrast, is not well-described by a Gaussian. It is asymmetric with a substantial excess of objects with high Stetson values. This is as expected, since members are more likely to have large amplitude, correlated variability.

Similarly, we can examine the distribution of the Stetson index for our standard set of members that are also in the standard set for statistics (and having sufficient points in both bands) and compare it to the Stetson index distribution for the remaining objects not selected as members (but still in the standard set for statistics, and having sufficient points in both bands). We obtain a similar result in the right panel of Figure 23; the Stetson distribution for members crosses that for non-members at about 0.9, or perhaps a little below that level, suggesting that this division is the dominant cause of the break in the entire distribution at about that level. Even if our separation between members and non-members is imperfect (which it certainly is), these distributions are consistent with our selection of 0.9 as the cutoff. We conclude that a Stetson index cutoff of 0.9 is a sensible boundary for our data set.

While the objects with Stetson indices of ≲0.4 have a very low chance of having legitimate correlated variability, and objects with Stetson indices >0.9 have a high likelihood of correlated variability, there is a continuum between these values. Objects with Stetson indices ≳0.4 and ≲0.9 have low-confidence for correlated variability. In MC11, we took a different Stetson index of 0.55 as the cutoff, based on the distribution of Stetson indices for that particular data set. As such, some of the identified low-confidence variables may have changed between this and the initial analysis. We also note that the cutoff in Stetson index for CSI 2264 is very different (Cody et al. 2014), but that program has a substantially different observing cadence than the YSOVAR-classic data discussed here. In general, the appropriate Stetson index cutoff must be determined for the individual data set, and there is no universal value.

A second method to identify variables is a chi-squared test (χ²), which, for a given band, is given by

$$\begin{equation} \chi ^2=\frac{1}{N-1}\sum \limits _{i=1}^{N} \frac{({\rm mag}_i-\overline{{\rm mag}})^2}{\sigma _i^2}, \end{equation} \tag{ 3 }$$

where σ_i is the estimated photometric uncertainty (corrected as per our discussion of the YSOVAR noise floor in Section 2.5).

This test is used to identify objects with uncorrelated variability, or variability in only one band (perhaps because data exist in only one band). This makes the χ² test more susceptible to instrumental issues affecting only one band. However, to demonstrate that it generally does a good job of recovering variables with large Stetson indices, Figure 24 shows the distributions of $\chi ^2_{I1}$ and $\chi ^2_{I2}$ as a function of Stetson index. For those objects in the standard set for statistics where it is possible to calculate both χ² and the Stetson index, the values are reasonably well correlated for the unambiguously variable objects. Using this plot, we find that a limit of $\chi ^2_{I1}$ or $\chi ^2_{I2}\sim 5$ is an appropriate conservative cutoff for potential variability in those cases where only one χ² can be calculated (e.g., where monitoring in only one band is available).

Zoom In Zoom Out Reset image size

For our largest data set (Orion), there are thousands of light curves that meet the requirement imposed by the Stetson index of having data in both IRAC channels. However, in the 11 smaller-field YSOVAR-classic data sets, imposing such a two-band restriction typically means that more than half the viable light curves would be discarded. Thus, the χ² test is particularly useful in these cases where only one band is available.

Figure 25 shows histograms of $\chi ^2_{I1}$ and $\chi ^2_{I2}$ for the standard sample for statistics, as well as for the subset of objects for which a Stetson index can be calculated (e.g., the sample used in the prior figure), and the smaller subsample of objects identified as variable using the Stetson index (Stetson index >0.9). We fit a Gaussian to the sharply peaked distributions, and found a 3σ value of χ² ≲ 4.5 in all cases. The bulk of the χ² distribution for which the Stetson index is >0.9 also has χ² > 5. To be conservative, we thus set a limit of $\chi ^2_{I1}$ or $\chi ^2_{I2}=5$ for identifying a candidate variable object. In the ideal case where there are multiple bands of monitoring data, one could assess each light curve with a large χ² but small Stetson index for physical plausibility. To attempt to avoid identifying false variability due to instrumental effects, the light curves for each of these potentially variable objects will be examined by hand in the context of the individual cluster analysis to come in separate papers.

Zoom In Zoom Out Reset image size

Finally, as in MC11, there are still legitimately variable sources within our standard statistical sample that fail both the Stetson index test (perhaps because only one band is available) and the χ² test (perhaps because the amplitude of variability is small and our limits for identifying variability were conservative by design). There are many mathematical tools available for identifying periodic behavior in an unevenly sampled time series. The last test for variability we run here is a periodogram analysis using the NASA Exoplanet Archive Periodogram Service³⁴ (Akeson et al. 2013). This service provides period calculations using Lomb-Scargle (LS; Scargle 1982), Box-fitting Least Squares (BLS; Kovács et al. 2002), and Plavchan (Plavchan et al. 2008b) algorithms. These methods are varyingly more or less sensitive to periodic behavior shaped like sinusoids or flat-bottomed transits, and/or may be less sensitive to periodic behavior appearing in addition to other behavior, such as a period superimposed on a long-term trend. The expected periodic variability in our sample includes anything repeated, from a sinusoidal-like signal originating from hot or cool spots on a photosphere, to signals characteristic of close binaries, to repeated dips in the signal (like "dippers" or AA Tau; see, e.g., MC11), or even pulsations (e.g., Morales-Calderón et al. 2009).

Specifically, because of the variety of expected light curve shapes, and the weaknesses inherent in any of these methods for finding periodicity, and noting the approach used by (and results from) McQuillan et al. (2013a, 2013b), we also calculated the ACF for each light curve as a check on repeated patterns. We linearly interpolated the light curve onto evenly spaced times, and then calculated the ACF using the following expression where L is a lag in days, and x is the light curve (with elements x_k):

$$\begin{equation} {\rm ACF}_{x}(L)={\rm ACF}_{x}(-L)=\frac{\sum _{k=0}^{N-L-1} (x_k-\overline{x})(x_{k+L}-\overline{x})}{\sum _{k=0}^{N-1}(x_k-\overline{x})^2}. \end{equation} \tag{ 4 }$$

We experimented with several different timescales as obtained from the ACF, and settled on the location of the first peak, providing that the peak was above an ACF value of 0.2. For those objects with significant periods, this coherence time should be well-matched to the period.

We looked for periods in light curves where we had at least 20 points (more restrictive than our standard statistical sample); we ran all four methods (LS, BLS, Plavchan, and ACF) on not just the 3.6 μm and 4.5 μm light curves, but also, where possible, the [3.6]–[4.5] light curves. In some cases, a long-term trend (astrophysical, not instrumental) is present in the individual I1 or I2 light curve, masking a periodic signal, but the color exhibits the periodic signal. We looked for periods between only 0.05 and 15 days, given the overall sampling of our data, and we require at least two complete periods over the typically ∼40 day window of our observations. We investigated phased light curves for those periods calculated using all of these methods. Based on these many thousands of results, we concluded that LS is the best, for our data set, for finding reliable, plausible periods. BLS and the Plavchan algorithm, while they look for a wider variety of shapes of signals, struggle with light curves that typically have less than 100 points, as ours do; the ACF approach finds only the strongest signals. Typically, if the LS algorithm found a reliable period, those other three approaches found comparable periods.

Thus, we filtered first on the LS results. We excluded candidate periodic objects if the calculated false alarm probability (FAP; see, e.g., Scargle 1982) was >0.03, or if the recovered period was >14.5 days and the FAP for that period was >0.01, or if the period was <0.1 day (slightly larger than the lower limit over which we searched), or if the calculated period was exactly 15.0 days (by inspection of the light curves, input and phased, a returned period exactly equal to the upper limit of our search window was usually indicative not of a true periodicity, but instead of a long-term trend in the data). For each of the remaining objects, we investigated the phased light curve. Perhaps unsurprisingly, given our overall FAP cutoff of 0.03, about 3% of the surviving candidate periodic light curves did not produce physically plausible phased signals. Those objects were omitted from the final set of periods, and will be identified as such in the corresponding cluster papers. The planned individual cluster papers may include a few additional periodic objects not automatically identified due to the presence of outlying photometry points which mask periodicity unless removed by hand.

We proceed to include all of the periodic objects identified via the LS algorithm (excluding the candidate periodic objects as described) in the set of variable objects, even if they fail the other (Stetson, χ²) variability tests. Individual objects will be discussed in the papers dedicated to each individual cluster, but anywhere from 1 to 15 objects, typically ⩽5, were added to the list of likely variables for each of the smaller-field clusters. For any given object, we wish to assign a single period to that object. We take any period derived from the [3.6] data first ([3.6] is less noisy than [4.5]), then, only if there is no [3.6] period of sufficient power, we take the period derived from [4.5], and finally, if no other period of sufficient quality is available, then we take that derived from [3.6]–[4.5].

A preliminary list of periodic objects in Orion appeared in MC11. The approach we are now using to search for periods is more stringent than that in MC11. Our current approach recovers the bulk of the objects that MC11 reported as periodic, but does not, for example, recover the objects reported as having P > 15 days. A complete list of the ∼800 periodic variables in Orion as derived from this YSOVAR data reduction will appear in a later paper.

Above, we described how variable sources are identified using the Stetson test, the χ² test, and a search for periodic variability. The cut-off values for those tests were chosen to yield a conservative list of variable sources and to reject those sources where the variability stems mostly from instrumental artifacts. We now present Monte-Carlo simulations to quantify how much variability is required to meet those criteria.

Several different physical effects may contribute to the observed variability and this can lead to very complex patterns in the light curve. In the absence of a theoretical model to explain the different contributions, we concentrate on simple analytical prescriptions for light curves so as to gain a sense of the sensitivity of our statistical tools. First, we consider a source that has two states, a bright and a faint state. The light curve switches randomly between those two states. We simulate light curves for a different fraction of time spent in the upper state (0%, 10%, 20%, 30%, 40%, and 50%) and we expect that variability is more easily found for sources with a larger amplitude between the two states and an equal chance to find the source in each state. Second, we simulate sinusoidal light curves with different periods—0.1–2 days (in steps of 0.1 day) and 2–20 days (in steps of 2 days). We add Gaussian noise to each light curve and vary the ratio of the signal amplitude and the noise (from 0 to 10 in steps of 0.1, where a relative amplitude of 0 means a constant light curve with noise only). The sensitivity of the Stetson test, the χ² test, and the Lomb–Scargle periodogram is independent of the magnitude of a source (that is not noise-dominated); only the relative amplitude of the signal and the noise plays a role. Thus, it is equally possible to detect strong variability in a weak source (with large photometric uncertainties) as weak variability in a bright source (with small photometric uncertainties).

For each grid point in relative amplitude and fraction of time in the upper level or period, for each band, we simulate 10,000 light curves for each grid point for each band. The light curves are sampled at the time intervals of the actual YSOVAR observations. Figure 26 shows the histogram of the time steps between observations for each cluster (the typical min and max Δt were given above in Table 3). For all clusters, the sampling is non-uniform to avoid aliasing for a specific period. However, the histograms fall in two groups, with AFGL 490, Mon R2, and NGC 2264 having a lower overall sampling rate than the other clusters. Figure 27 is another representation of the sampling rate. In the style of Figure 14, it visually represents the sampling rates for the fast cadence monitoring. Here, AFGL 490, Mon R2, and NGC 2264 can also be seen to have a lower overall sampling rate than the other clusters. We ran all Monte Carlo simulations using the actual time sampling from L1688, one of the clusters with a high sampling rate and Mon R2, a representative cluster with a lower sampling rate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 28 shows results from the Monte Carlo simulations using the sampling of the L1688 cluster. The left panel presents the detection efficiency for light curves from an object with two distinct luminosities. If the star is found in each state half the time, the Stetson test will identify it as variable in almost all cases (99%) if the step size is at least three times larger than the noise level. Since the χ² test uses data from one band only, the step size must be larger (five times the noise level) to reach the same detection efficiency. If the star spends less than 20% of the time in either state, there is a reasonable chance that the variability will not be found, even for larger step sizes, since the sampling might find only few data points in this state. For a periodic light curve (middle panel), variability is again found more easily in the Stetson test than in the χ² test, and the period of the variability does not influence the detection efficiency, since those two tests do not consider the time ordering of the observed data. Using the LS approach, we are sensitive to periods (P) between about 1 and 15 days. Even periods where the amplitude (a) of the light curve asin (2π/Pt) is only twice as large as the noise level are easily detected, almost independent of the period in the range 1–15 days. Such weak signals would not necessarily show up as variable in the Stetson or χ² test (see the middle panel). In general, the region where the tests detect variability in some light curves but not in others is fairly narrow.

Zoom In Zoom Out Reset image size

Figure 29 shows the same plots as Figure 28 but for the time sampling of Mon R2, one of the clusters with a lower cadence in the observations. The general shape of the regions in the parameter space where variability can be detected is the same, but, due to the lower number of observations and the larger time span between observations, a larger amplitude is required, and we are less sensitive to periods below about two days.

Zoom In Zoom Out Reset image size

To explore more complicated light curves, we combined several effects, e.g., a sinusoidal periodicity overlaid on a long-term trend. In these cases, the strongest effect generally determines how the variability will be seen. If the magnitude of the trend is large and the amplitude of the sine wave is small, then the light curve will be marked as variable, but the periodicity might not be detected.

Given the possible complexity of real light curves, it is not possible to cover the entire parameter space with Monte Carlo simulations. However, the scenarios presented here show that, in general, we can detect periodicity with an amplitude just twice the level of the noise, variability at 3–5 times the level of the noise with the Stetson test and about 6–10 times the level of the noise for single-band light curves with the χ² test. Similar results for the relative sensitivity of the Stetson and the χ² test were also found by Flaherty et al. (2013), although they used different cut off levels than this work.

As seen above in Figure 16 and in MC11, there is a reasonably strong correlation between the mean magnitude of a source and its mean error. This overall relation affects our ability to find variability or periodicity. For sources brighter than 13th mag, the total uncertainty is dominated by the error floor introduced in Section 2.5. For those brighter sources, we can detect periodicity if the amplitude is larger than about 0.02 mag and variability if it is larger than about 0.03–0.1 mag (with the exact number depending on the signal shape).

In addition to the variability probed on the YSOVAR monitoring timescale of ∼40 days, we also are interested in the evidence for longer-timescale variations between observations of these same clusters in the Spitzer cryogenic epoch and the post-cryo (YSOVAR) epochs. To identify the variables in this case, for every object with a light curve, we can compare the average measurement from the earliest cryo era (Table 4), and the mean for that object over the YSOVAR standard statistical sample (Section 3.2).

The process we used to identify the long-term variables in each cluster is shown for AFGL 490 in Figure 30. We plot the difference between the cryo and post-cryo measurements for each object as a function of the cryo value and determine, for each cluster, the brightness at which photometric noise clearly dominates. This faintness limit is close to 16th mag, consistent with what we noted in Figures 19–22. We select this limit separately for each cluster and consider only objects brighter than this limit. We fit a Gaussian to the histogram of the difference between the cryo and post-cryo measurements, allowing the zero point as well as the height and width of the Gaussian to be free parameters. We classify as long-term variables all objects with cryo to post-cryo offset further than 3σ from the peak of the (fitted) distribution. All objects from the standard statistical sample, not just the members (or only the variables), go into this process of defining the width of the distribution. The fraction of objects in each field that are classified as variable, and the subset of variables that are also cluster members, are both identified after the long-term variables are identified. A summary of the important parameters in these analysis steps to search for variables over this long-term baseline is in Table 7; the values we used for Δt, the time lapse between the cryo and post-cryo observations, are included in Table 4.

Zoom In Zoom Out Reset image size

Our YSOVAR map in Orion is far larger than the other cluster maps, which focus on the most embedded (possibly youngest) objects in these clusters. For these embedded objects, it is not generally possible to obtain a rotation period from ground-based optical or NIR data due to extinction. Moreover, for stars with more significant disks, it is less likely that the IR light curve will be strictly periodic. Many distinct processes can contribute to a YSO's mid-IR variability, and it often results in a stochastic light curve (Cody et al. 2014).

For the objects for which we can derive a period, we would like to compare our values to those from the literature as a check on our methodology. Since our clusters are, for the most part, very embedded, there are not very many known periods in the literature. As discussed in Section 2.1.2, there are many periods for Orion and NGC 2264, typically obtained in the optical, but with some values from the NIR. Parks et al. (2014) reports on NIR periods from objects in our region of L1688; there are other literature values for rotation periods of objects elsewhere in L1688, beyond our monitored region. We can roughly compare to the MIR timescales reported in Morales-Calderón et al. (2009) for IC 1396A.

Of the objects with periods in the literature, there are ∼200 that also have periods derived here in Orion (excluding those from MC11, since those were derived from the same observations we use here), and an additional ∼15 from NGC 2264, L1688, and IC 1396A. About 75% of those have period measurements that match to better than 10%, so we have confidence that our period-finding approach is at least well-matched to those in the literature. Figure 31 plots the YSOVAR-determined period against the literature period for those objects where it is possible. The clusters that are not Orion (NGC 2264, L1688, and IC 1396A) are plotted separately simply because Orion dominates the statistics, and it is useful to see if there are good matches outside Orion as well as within Orion. Three of the four objects from the three smaller-field clusters that are not well-matched to the YSOVAR-determined period are close to likely harmonics, and the periods are of comparatively low quality. The one that is most discrepant is from NGC 2264, SSTYSV 064101.40+093408.1, and is being compared to a period from Lamm et al. (2004). Our phased light curve looks correct (for our wavelength and epoch of observation). Of the ∼50 Orion periods that do not agree to 10% (out of ∼200 Orion period comparisons total), it is predominantly the case (by a ratio of 3 to 1) that the YSOVAR period is longer than the literature period. About 60% of these have [3.6]–[8] >0.8. All but five of those Orion periods were optically determined. Because we are working in longer wavelengths, it is possible that, particularly in those cases, we may be sampling a different location in the star-disk system, e.g., further away from the photosphere, where Keplerian rotation periods are longer.

Zoom In Zoom Out Reset image size

A relation has already been found between IR excess and rotation rate for young stars, suggesting that IR excess and rotation rate are related; out of our 12 clusters, this relation has been found in Orion (Rebull et al. 2006) and NGC 2264 (Cieza & Baliber 2007). In the 11 smaller-field clusters (i.e., all but Orion), there are ∼350 stars with periods measured from YSOVAR light curves, but only ∼250 of those also have cryogenic Spitzer measurements at 3.6 and 8 μm from which we can get a clear indication of the IR excess in these systems. There are ∼800 stars in Orion with measured YSOVAR periods, but only ∼430 have [3.6] and [8] measurements.

Figure 32 shows the relationship between IR excess (specifically [3.6]–[8]) and YSOVAR-derived IR periodicity for these sources. In both cases, there is a gap near [3.6]–[8] ∼ 0.8, which divides the disk candidates (above that cutoff) from the non-disk candidates (below). There is also different behavior to the left and right of log (P) ∼ 0.25, or P ∼ 1.8 days—excesses do not necessarily imply longer periods, but a star with a longer period is more likely than those with shorter periods to have an IR excess. Figure 33 shows the cumulative distributions of [3.6]–[8] for the same two panels as in Figure 32, for three different bins of log(P), divided at log(P) = 0.25 and 0.75 (1.78 days and 5.62 days, respectively). According to Kolmogorov–Smirnov (K-S) tests, the distributions of [3.6]–[8] are significantly different within each panel. The two distributions that are the most similar are the full (all 12 clusters) distributions for 0.25 < log(P) ⩽ 0.75, and 0.75 < log(P); the probability that those populations were drawn from the same distribution is 4%. The probability that the populations were drawn from the same distribution for log(P) ⩽ 0.25 and 0.75 < log(P) (again, for all 12 clusters) is ∼10⁻¹³; for log(P) ⩽ 0.25 and 0.25 < log(P) ⩽ 0.75, the probability that the populations were drawn from the same distribution is ∼10⁻¹⁰. Similarly, Figure 34 shows the cumulative distributions of log(P) for the same two panels as in Figure 32, for two different bins of [3.6]–[8], divided at [3.6]–[8] = 0.8. Again, according to K-S tests, the distributions of log(P) are significantly different within each panel; the probability that either of the populations were drawn from the same distribution is <10⁻¹⁷.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The plots in Figure 32 are very similar to that obtained for Orion by Rebull et al. (2006), and that by other investigators in other clusters, despite the fact that optically determined periods were used there. The periods derived from our IR YSOVAR data may be photospheric rotation rates, pulsation rates (see, e.g., Morales-Calderón et al. 2009), or inner disk rotation rates (see, e.g., Artemenko et al. 2012); on the other hand, for those clusters where there are periods available, we match the literature reasonably well, and the literature for the most part is using optical data to obtain periods. Therefore, it seems that we are, in most cases, not sampling much different locations in the star-disk system. However, an exception could be that the optical and IR observations are sampling two separate places whose movements are locked together, such as starspots on the photosphere and stellar-magnetosphere-driven disk disturbances at the corotation radius.

It is also surprising that the results for the aggregate set of clusters are so similar to that for Orion, because the clusters should be for the most part substantially younger than Orion. This could imply that disk locking may be in effect at even these young ages, or that accretion-powered stellar winds are the dominant mechanism to slow these objects. However, it is likely that we can obtain viable periods more easily for relatively unobscured stars, e.g., with these periods, we are also sampling the older end of the young star distributions in these clusters. Little is known about many of the objects outside of Orion shown here; additional study of the individual objects will help clarify matters. Individual objects will be discussed in the corresponding YSOVAR cluster paper.

There have been recent reports of debris disks undergoing significant short-term changes; both Meng et al. (2012) and Melis et al. (2012) report on systems that change significantly at wavelengths >10 μm over timescales of years. Our objects are much younger (a few million years rather than a few tens or hundreds of million of years) and our monitoring wavelengths are considerably shorter. However, Rice et al. (2012) also note that nine (36%) of the stars in their sample of young stars in Cygnus OB7 (comparable in age to our sample) have a transient NIR (JHK) excess. We can use this first look at our data to constrain the degree to which IR excesses in our sample vanish (or appear) on the timescales of years, namely between the cryo-era observation and that of our post-cryo observations.

Irrespective of whether or not objects have been identified as variable above, we compared the cryo-era [3.6]–[4.5] color with the maximum and minimum [3.6]–[4.5] color obtained during our YSOVAR (fast-cadence) monitoring (standard set for statistics, the subset of which have measurements in both channels). Out of ∼11,000 objects (cluster members as well as background objects included) for which we have [3.6]–[4.5] color light curves, there are at most 15 objects that seem to have legitimate substantial changes to the [3.6]–[4.5] color (changes of a size that might be consistent with big changes to a disk), and these are all relatively faint ([3.6] > 12 mag) objects. At most, two of those cases have an IR excess that appears to be possibly transient on these timescales, so at most <0.02% frequency of occurrence. For the remaining 13 objects with plausibly real changes in color, the disk is still clearly present, but the brightness and color have changed substantially. Individual objects will be discussed in the cluster papers.

In the literature (e.g., Giannini et al. 2009; Antoniucci et al. 2014), constraints have been placed on the timescales for brightening or fading by comparing how many sources are found to be getting brighter or getting fainter (e.g., for each source, given the two epochs, for how many cases is the second epoch brighter than the first, and for how many cases is the first epoch brighter than the second). If there are random fluctuations in brightness, the same number of sources should get brighter as get fainter. If, instead, there are more fading sources than brightening, then the type of variability may be characterized by a short rise and a long fall.

We can make a similar comparison among our long-term variable sample. To reduce scatter from noise, we consider only the standard set of members (Section 3.1), and consider only those objects identified as long-term variables independently in both [3.6] and [4.5] (Section 5.5). Figure 35 compares the numbers of these remaining sources that become brighter at both [3.6] and [4.5] with those that become fainter at both [3.6] and [4.5] for 11 of the clusters; errors are approximated by simple Poisson counting statistics. A significantly larger number of variables are found in Orion, with essentially equal numbers of sources getting brighter as getting fainter. For the smaller-field clusters, the numbers of sources that brighten is less consistently equal to the number of sources that fade, but there are also far fewer sources to count. Summing up the 11 smaller-field clusters, there are 107 (±10) sources that brighten and 88 (±9) sources that fade; these numbers are within 2σ of each other, but both numbers have to be extended ∼1σ toward each other. Fitting a line to the points in Figure 35, including Orion, results in a line of slope 1.05 ± 0.12, and an intercept of −2.8 ± 3.3. This is consistent with no difference between the brightening and fading sources, but with only a ∼1σ possibility that there are slightly more brightening sources. With or without Orion, there are very similar numbers brightening and fading.

Zoom In Zoom Out Reset image size

In the case of several other papers in the literature, they were looking for much more significant variability than we are finding (e.g., FUors and EXors). They found internally consistent patterns of sharp rises and long falls in brightness. For our sample, we are apparently finding more varieties of variability, even on the long-term, such that the timescales average out over six to seven years to have no significant bias toward brightenings or fadings.

Having identified the long-term variables above (Section 5.5) for each cluster, we can look at the fraction of objects that are variable on the longest timescales we sample. Figure 36 shows the fraction of variables for all objects in each field (the standard set for statistics), including both members and likely field objects. The x axis is our parameterization of the relative fractions of embedded sources (see Section 4.1), the ratio of Class I/total number of members. Figure 36 thus is somewhat incongruous in that the standard set of members is used in the x axis, but the y axis includes everything in the FOV. About 10%–20% of all objects with light curves are identified as variable in the long-term. However, as can be seen in the figure, the variability fraction of everything in the field is not a strong function of the Class I/total ratio—the slopes of the best-fit lines in Figure 36 are both small. (The slopes of these lines are given in the figure caption.) Calculating Pearson's correlation coefficient also suggests that there are no significant correlations shown in Figure 36. We expect that the fraction of stars that are variable should be a function of Galactic latitude (see approximate Galactic latitudes listed in Table 1), because the fraction of sources that are background/foreground stars will be higher in the Galactic plane, so the fraction of cluster members will be higher out of the Galactic plane, and since any young star (cluster member) is more likely to be variable, there should be a higher fraction of variable objects at higher Galactic latitude. Indeed, the Pearson's correlation coefficient suggests that the fraction of long-term variables for all objects in each field is strongly correlated with the absolute value of the Galactic latitude.

Zoom In Zoom Out Reset image size

Figure 37 recasts Figure 36 by using the long-term variability fraction of only the standard set of members, rather than everything in the field. This time, there is a significant correlation; the higher the fraction of embedded members, the higher the fraction of long-term variables. If a cluster has more sources that are embedded and likely to be actively accreting and interacting with their circumstellar material, it also has more sources that are variable on timescales of years. The slopes of the best-fit lines shown in Figure 37 are given in the figure caption. Pearson's correlation coefficient suggests that there is a correlation here, and it is stronger in I1 than I2 (consistent with the calculated slopes).

Zoom In Zoom Out Reset image size

Serpens South has the highest fraction of the most embedded sources; to test if it is providing a significant "lever arm" on this fit, we fit the remaining points omitting Serpens South, which does indeed weaken the correlation, as can be seen in the figure.

We tested the correlations seen in Figures 36 and 37 using a variety of other SED class-based parameterizations, and we found similar results—there is no significant correlation of the long-term variable fraction of all sources in the field with any SED class-based parameterization, and there is a correlation of the long-term variable fraction of the members with any parameterization chosen that uses fractions of various SED classes (or groups of classes). This correlation between the fraction of members that are variable on the longest timescales and the parameterization of "embeddedness" seems robust. We expected that young stars were more likely to be variable than field stars. Moreover, we have found that within the category of young stars, for a higher fraction of embedded sources (a higher fraction of presumably younger sources), we find a higher fraction of long-term variables. This is consistent with what has been found in individual clusters (e.g., NGC 2264 (Cody et al. 2014) and L1688 (Günther et al. 2014)).

The difference between fractions of embedded objects among these clusters is not the only potentially significant factor for this set of observations. Here, we sample timescales of ∼4.5 to ∼7.5 yr. If the amplitude of variability of the members increases as the time baseline increases, we expect more members to be selected as variable in the long-term, and thus a higher long-term variability fraction as the time baseline increases. However, Figure 38 shows this relationship between the variability fraction and the time between epochs of observation (from Table 4). The best-fit lines and correlation coefficients are consistent with no significant effect on the variability fraction as a function of timescale sampled. Serpens South, because it was observed (indeed, discovered; Gutermuth et al. 2008b) comparatively late in the Spitzer mission, samples the shortest timescales. If the Serpens South point is omitted from the fit, the slopes become slightly steeper, but there is still no significant correlation from the correlation coefficients.

Zoom In Zoom Out Reset image size

Scholz (2012), working in the K band, finds that the longer one monitors a cluster, the larger the amplitude of variability, on timescales of years. That work specifically investigated the amplitude of the change in magnitude, for just one cluster at a time, and looked at the change of the range of that distribution of amplitudes as a function of time step. For the times we sampled on these longest timescales, analogous plots do not show a significant change in the median or the top quartile or the 90th percentile. We tried using the most stringent set of objects (only those in the standard set of members, standard set for statistics, that had light curves in both IRAC channels); we still did not find this effect. For this analysis, Scholz (2012) was only working in ρ Oph, over a larger area than we were, using slightly different wavelengths, and a larger range of timescales, the maximum of which (∼2000 days) is comparable to the minimum Δt (∼4.5 yr ∼1600 days) we consider here. These could account for the observed differences.

In the context of understanding the evolution of young stars from a very embedded state to a less embedded state, the community has chosen to parameterize objects based on the slope of the SED (see, e.g., Wilking et al. 2001). This classification is tied to the empirical shape of the SED. Very embedded (presumably very young) objects will have SEDs that peak at long wavelengths; as the object sheds its natal cocoon, the peak of the SED moves to shorter wavelengths. Objects with substantial circumstellar disks will emit more energy in the IR (and longer wavelengths) than in the optical. Objects with less substantial disks will have most of their energy emitted in wavelengths shorter than the NIR, but there will be additional energy contributions in the IR (and longer wavelengths) from the circumstellar dust and/or debris, i.e., the SED has an IR excess. The most embedded phase is referred to as Class 0, then proceeding (based on SED shape) through Class I (rising SED), Flat (flat SED), Class II (falling SED with an IR excess), and finally Class III SEDs, which are objects with photospheric or near-photospheric SEDs, with little or no IR excess, typically identified as young via other means such as X-rays or Hα emission. Sometimes additional classes such as transition disks are added near the end of this sequence. Classical T Tauri stars (CTTSs) are often identified with Class IIs, and weak-lined T Tauri stars (WTTSs) are often identified with Class IIIs. Nomenclature is difficult and inconsistent across the literature; see Evans et al. (2009a) for a discussion of terminology.

We need to establish at least an internally consistent definition of the SED classes so that we can investigate trends as a function of SED class as a proxy for age. The reliability of the translation between SED class and age has been discussed at length in the literature (and will continue to be discussed in the future); other factors such as inclination and multiplicity may play a large role. In the context of our work, we wish to establish an internally consistent approach that can be calculated for all sources in the YSOVAR fields.

In order to assemble our SEDs for each object, we include all the data described in Section 2 between U and 25 μm. For the SED classes, only the IR bands are relevant.

In order to define an internally consistent placement of the YSOVAR objects into SED classes, in the spirit of Wilking et al. (2001), we define the near- to mid-IR (2–24 μm) slope of the SED, α = dlog λF_λ/dlog λ, where α > 0.3 for a Class I, 0.3 to −0.3 for a flat-spectrum source, −0.3 to −1.6 for a Class II, and < −1.6 for a Class III. For each of the objects in our sample, we performed a simple least squares linear fit to all available photometry (only detections, not including upper or lower limits, but including archival and literature data) as observed between 2 and 24 μm, inclusive. We included the mean obtained from the standard set for statistics of the light curve in the SED, in addition to the cryo-era Spitzer points. (This means that there could be more than one point contributing to the fit at 3.6 μm, one from the cryo era and one from the mean of the YSOVAR light curves.) Formal errors (either from individual single-epoch measurements or the mean and standard deviation from the IRAC light curves) on the infrared points are so small as to not affect the fitted SED slope. However, if the mean is calculated over more of the light curve than our standard set for statistics (fast cadence) sample, the mean may be different enough, in the cases of sparse SEDs, that the class may change to an adjacent class; these cases will be noted in the cluster papers where relevant. The linear fit is performed on the observed SED, e.g., no reddening corrections are applied to the observed photometry before fitting.

In the literature, the precise definition of α can vary, or the distribution of slopes and classes can vary (e.g., Sung et al. 2009 lists Class I, II, II/III, pre-transition disk, transition disk categories for NGC 2264), which may result in different classifications for certain objects. Classification via our method is provided specifically to enable comparison within this paper (and to other YSOVAR papers) via internally consistent means. Our classification, since it is based on observed SED, is possible for all objects (not just those identified as YSO candidates). The formal classification definition puts no lower limit on the colors of Class III objects (thereby including those with SEDs resembling bare stellar photospheres, and allowing for other criteria such as X-ray brightness to define youth). The SED slopes and classes for all of the sources of interest will appear in the individual cluster papers. Histograms of the relative fractions of each class for the standard set of members for each cluster appear in Figure 17.

In terms of aggregate statistics over all data from all 12 clusters (members and non-members together), we can constrain the fraction of objects that change class depending on whether or not the 24 μm point is included in the fit. There are about 21,000 objects with light curves over all 12 clusters for which an SED slope can be fit between 2 and 8 μm. Out of those ∼21,000, there are only about 1760 with MIPS-24 detections (not limits), so only about 8% of the sources are affected. Admittedly, those sources that are detectable at 24 μm are the ones with rising SEDs and thus are statistically more likely to be true cluster members than a source selected at random from the map. Out of the ∼1760 with MIPS-24 detections, ∼72% of them do not change class when the 24 μm point is included in the SED fit. The class bins of slopes are relatively large and thus relatively insensitive even to a point at the far red end of the SED. Of the ∼28% that do change class, ∼85% move only one step, to an adjacent class. As expected, there is a bias, when including the 24 μm point, to move the objects to a more positive slope, e.g., toward the more embedded end of the sequence; of the ones that change class at all, ∼61% move one step earlier (toward more embedded, not necessarily in age) in the sequence, ∼24% of those move one step later (toward less embedded, not necessarily in age) in the sequence.

We conclude that it does not make a significant difference for the overwhelming majority of the sources if one uses the slope between 2 and 8 μm or the slope between 2 and 24 μm. (Note that this is fitting all available points between these values, not a simple comparison of the two end points.) For any sources of interest in which the class might change depending on the inclusion of the 24 μm point, they will be noted in the individual cluster papers.

G09 provides placement into SED classes as part of the data presented there, and the same algorithm has been applied to our entire cryogenic-era catalog (Section 2.7). These classes are identified based on dereddened colors, and are only provided for objects thought to be young stars. Objects that are not thought to be young stars are identified as other things such as "PAH emission source"; objects missing bands such that the classification scheme cannot be run are also not classified. The SED classes we are using here are based on fits to the observed SED, not the dereddened SED, and are obtained for any source, regardless of its true underlying nature. To constrain the degree to which we might be introducing a bias by fitting the observed SED, we can compare, for some objects, the class obtained by G09 and by our mechanism above.

It is important to note that the classes provided by G09 are different than the classes we use here. G09 has Class 0, I, II, and II/III, but no Flat class. Here, we do not have Class 0; we have Class I, Flat, Class II, and Class III. Over all 12 clusters, there are ∼3500 objects for which both classifications are available; 75% of those do not change class, even with the different bins that are defined. Out of the ∼3500, about 14% are identified as being from a later (less embedded) class than G09, and about 11% are identified as being from an earlier (more embedded) class than G09.

G09 also dereddens the SEDs before placing them in classes; we are fitting observed SEDs. G09 does this to avoid a reddening bias toward youth—as discussed in Muench et al. (2007), the K_s −[3.6] color can be affected by A_V ∼ 40, though it takes A_V ≳ 200 to affect the [5.8]–[24] color. However, in our case, because we are fitting all available points between K_s and [24], not just subtracting two points in the SED, the influence of reddening on this overall slope in the best case (where there are K_s, four bands of IRAC, and one MIPS band), simulations suggest that we would need A_J ∼ 11 or A_V ∼ 40 before a Class III object would be misclassified as a Class II. Admittedly, this is a best-case scenario, where the SED is well-populated. For a source without MIPS 24 but just K_s through [8], we find that we need A_J ∼ 6 or A_V ∼ 21 before a Class III object would be misclassified as a Class II.

Our objects discussed in YSOVAR have to be bright enough to get good quality light curves at 3.6 and 4.5 μm, which is effectively brighter than 16th mag (see Figures 19 and 20), which means that our sources cannot generally have extremely high extinction. Out of the objects for which we have a Gutermuth-derived value for A_K and a light curve, ∼6% have A_K ≳ 2 (which corresponds to A_J ≳ 6), and ∼1% have A_K ≳ 3.5 (which corresponds to A_J ≳ 11). We conclude that any bias toward young objects is likely not substantial in our data set.

We have opted to use most often our SED class definition as described above for internal comparison and consistency. However, when discussing the Class II/Class I ratio reported in G09, we are using those classes from G09. In Sections 2.1.2 and 4.1, we discuss the Class II/Class I ratio derived in a very similar fashion as G09, using the G09 classes, but only for those objects with light curves in the standard set for statistics. This enables at least some comparison back to G09 and related works.