CSI 2264: SIMULTANEOUS OPTICAL AND INFRARED LIGHT CURVES OF YOUNG DISK-BEARING STARS IN NGC 2264 WITH CoRoT and SPITZER—EVIDENCE FOR MULTIPLE ORIGINS OF VARIABILITY^*

We now narrow our focus to the NGC 2264 sample of stars with infrared excesses indicative of circumstellar disks. Our main approach is to measure the slope, α, of the spectral energy distribution (SED) at near- and mid-infrared wavelengths, and to compare to the expectation of a bare stellar photosphere. Several different disk classification schemes exist in the literature. The most recent analysis of NGC 2264 members by Sung et al. (2009) involved the five classes I, II, III, pre-transitional, and transitional, based primarily on the loci of photometric points on infrared color–color diagrams. Since the definition of transitional disks has been a subject of debate and the set of infrared photometry is not complete for all objects, we prefer the SED slope definition (see Evans et al. 2009). These are systems that have outer disk properties of normal T Tauri stars but with inner disk holes or gaps.

Following Wilking et al. (2001) with some guidance also from Lada & Wilking (1984), Greene et al. (1994), and Bachiller (1996), we define α = dlog λF_λ/dlog λ for flux F_λ as a function of wavelength, λ. We classify α > 0.3 for Class I, 0.3 to −0.3 for flat-spectrum sources, −0.3 to −1.6 for Class II, and <−1.6 for Class III. The notation II/III is reserved for transitional type disks, as described below. For all objects in our Spitzer sample, we performed a least-squares linear fit to all available photometry (not including upper or lower limits) as observed between 2 and 24 μm. We obtained data at 5.8 μm and longer wavelengths from cryogenic Spitzer observations. To maintain consistency with the approach taken by the YSOVAR project (L. Rebull 2014, in preparation), we obtained archival IRAC and Multiband Imaging Photometry for Spitzer (MIPS) data (program IDs 3441, 3469, 50773, and 37) and performed our own photometric reduction with Cluster Grinder (Gutermuth et al. 2009). We produced photometry for all sources in the NGC 2264 region with a signal-to-noise ratio of at least 5. For objects in the Sung et al. (2009) sample, our absolute photometry agrees with theirs to within 1%. We note that formal errors on the infrared points are so small as to not affect the fitted SED slope. The fit was performed on the observed SED, with no reddening corrections. Accounting for reddening is unlikely to change the SED class, since all of our targets with CoRoT data are in low-extinction regions of the cluster. Variability could also alter the shape of the SEDs slightly, since not all of the photometry was obtained simultaneously. However, we find that an allowance of 15% uncertainty does not change the classification appreciably. We provide this slope assessment, as opposed to that of Sung et al. (2009), to ensure internal consistency and enable comparison with other clusters in the YSOVAR sample (see L. Rebull 2014, in preparation).

For objects with MIPS data (131 in the combined CoRoT/Spitzer sample), we compared the α value and resulting disk class derived with and without the 24 μm point. Disagreement occurred for around one-quarter of these, primarily due to disjoint 24 μm photometry compared to the shorter wavelength SED. In the majority of these cases, the disk is transitional. However, in others, it is possible that nebulosity is causing a false excess at long wavelengths, and we have adopted the α value from 2–8 μm data. For all disagreements, we visually inspected the SEDs to determine the most appropriate wavelength range for slope determination. We identified a total of 140 disk-bearing objects among NGC 2264 members targeted for observation with both Spitzer/IRAC and CoRoT.

There are additional cases where misclassifications may have occurred, such that disk-bearing stars are labeled as diskless. Several objects display high-amplitude variability, but their mid-infrared SED slopes of less than −1.6 result in Class III status. To assess whether they might have weak disks, we examined their [3.6]–[8.0] IRAC colors. Analysis of infrared photometry in other clusters (e.g., Cieza & Baliber 2006; Cody & Hillenbrand 2010) has shown that the requirement [3.6]–[8.0] > 0.7 is a fairly robust disk selection criterion, whereas [3.6]–[8.0] < 0.4 selects diskless stars with high accuracy. Using the 0.7 cutoff, we identified 23 additional stars in the Spitzer and CoRoT fields with weak or transitional disks, which are confirmed by visual inspection of the SEDs, revealing infrared excess above the expected photospheric flux level predicted from the stellar spectral type. We adopt the disk class "II/III" for these objects and present them along with the rest of the disk-bearing sample in common with CoRoT observations in Table 3. We also list the 2MASS and CoRoT cross matches here.

The full sample of ∼1500 cluster members in NGC 2264 includes spectral types ranging from M5 to A7 (Walker 1956; Makidon et al. 2004; Dahm & Simon 2005), corresponding to masses of ∼0.1 M_☉ to several M_☉. Restricting this to the 162 disk-bearing objects among the IRAC and CoRoT targets, the spectral types run from M5 to G0. We present the distributions of brightness in CoRoT and IRAC bands for this disk-bearing sample in Figure 2. Both the spectral type and R magnitude distributions are representative of the cluster sample as a whole. However, the collection of 3.6 μm points is skewed considerably toward brighter values than compared to the available infrared data set. This is because the requirement of CoRoT observations eliminates most faint and embedded objects from the set. We have also produced color–magnitude diagrams to compare the 162-object sample with the distributions of field stars and other cluster members; these are presented in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The pixel sizes for Spitzer/IRAC and CoRoT are 12 and 23, respectively. We are therefore unable to separate some visual binaries. However, thanks to the higher spatial resolution of previous ground-based data sets (i.e., Sung et al. 2008), we have a fairly complete list in Table 3 of blended objects. The effect of binarity on variability will need to be assessed once a more complete multiplicity census is available.

Spitzer has been operating in its Warm Mission mode since the exhaustion of cryogen in mid-2009, and it now observes exclusively in the IRAC 3.6 μm and 4.5 μm channels. Observations of a ∼08 × 08 region of NGC 2264, including the more embedded Cone Nebula and Christmas Tree Cluster regions, were carried out from 2011 December 3 to 2012 January 1 with Warm Spitzer/IRAC, under program 80040 (PI: J. Stauffer). The field center was R.A. = 06^h40^m450, decl. = +09°40'40''. Since the cluster is considerably larger than the camera's 52 × 52 FOV for each channel, we visited targets ∼12 times a day with a series of mapping fields. For each such Astronomical Observation Request (AOR), we executed four pointings. Starting at the first position, two exposures were acquired with an 8'' dither to mitigate any detector artifacts or cosmic-ray hits. Total frame time for each was 12 s. Data were taken simultaneously in the two IRAC bands, but due to the non-overlapping IRAC FOVs, individual objects were monitored in one channel at a time. A spatial shift of just over half the field size was then performed, providing a second pair of dithers for half of the targets and a first dither pair for an equivalently sized set of targets in the newly observed section of the FOV. After three telescope offsets, the FOV has advanced by 72, such that targets in the newly observed section of the "trailing" IRAC dectector's FOV will have been observed by the "leading" IRAC detector less than a minute earlier. Accordingly, most of our targets have two-color IRAC data, but for the approximately 40% of objects near the edges of the mapping region, there is only single-channel coverage.

Since stars brighter than magnitude 9.5 saturate the IRAC detector in 12 s exposures, we also acquired 0.6 s integrations of each field. This high dynamic range (HDR) mode was employed during every sixth round of mapping. The resulting observing cadences were every ∼101 minutes for mapping data and ∼12 hr for the 0.6 s HDR mode data.

To capture light-curve behavior on even shorter timescales, we operated IRAC in staring mode for four blocks of 20, 26, 16, and 19 hr, respectively, toward the beginning of the run (see Table 1 for dates). The channel 1 and 2 FOVs were fixed on a region near the center of NGC 2264, and exposures with 6 s integration times were acquired repeatedly, with no dithering or HDR mode images. After taking readout and data volume restrictions into account, the cadence for staring mode data was ∼15 s. A total of 549 stars were observed in the staring fields. We display these regions, along with the full 3.6 and 4.5 μm FOVs, in Figure 1.

For both the mapping mode and the staring mode data, we used the IDL package Cluster Grinder (Gutermuth et al. 2009) to generate light curves from the BCD images released by the Spitzer Science Center (SSC) pipeline, version S19.1. Each flux-calibrated BCD was processed for standard bright source artifacts. Cluster Grinder delivers mosaics for each AOR, point-source locations, and photometric measurements from the mosaics ("by-mosaic" light curves). We then re-performed photometry on the BCD images using the Cluster Grinder source list, with a 5σ source detection threshold. We applied array-position-dependent systematic corrections for residual gain and pixel-phase effects (although the treatment of this was modified for staring data; see the Appendix). The pipeline also computes a variety of flagging information, including maximum pixel value for saturation detection, spatial coverage, and outlier rejection.

Aperture photometry was obtained from both the BCDs ("by-BCD" photometry) and the mosaics ("by-mosaic" photometry) using an aperture radius of 24 (2 pixels) and a sky annulus from 24 to 72 (2–6 pixels). We also combined the resultant photometric by-BCD measurements into "by-AOR" photometric products, mimicking the cadence of the by-mosaic photometry. We have found that the photometric precision achieved with this approach is higher than that of the by-mosaic data, presumably because the latter does not allow for easy treatment of array-position-dependent systematics, such as residual gain and pixel-phase effects. Hence, we present only the by-BCD light curves (for staring data) and by-AOR light curves (for mapping data) here.

We flagged by-BCD photometry for saturation, low signal-to-noise ratio (S/N < 5), and outlier status (5σ above the median flux of other points in the same AOR). Affected points were omitted from by-AOR photometry, and unaffected points were combined via an unweighted mean. By-AOR photometry values resulting from fewer than three by-BCD measurements were flagged as well. Following these procedures, the majority of mapping mode light curves contain 300–320 valid data points.

For pairs of 0.6 s and 12 s exposures taken during the HDR mode, we selected the latter only if it did not exceed a count level of 20,800 DN (Warm IRAC saturates at ∼30,000 DN); otherwise, the measurements from the 0.6 s exposures were swapped in. For saturated objects brighter than a magnitude of ∼9.5, we retained only the twice daily short HDR exposures, resulting in by-AOR light curves with ∼50 points.

We adopted the standard zero points derived from the official Warm Mission FLUXCONV header values: 19.30 and 18.67 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).

Uncertainties for the by-BCD photometry were derived by combining three terms in quadrature: 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). Finally, the BCD-level photometric uncertainties were combined in quadrature to yield uncertainties for the by-AOR photometry reported here.

To weed out extremely faint sources and artifacts, we required that each light curve contain at least five unflagged data points and have a mean photometric uncertainty of less than 10%. Applying these cuts, we generated a total of 13,856 mapping mode light curves in the 3.6 μm band and 12,186 light curves in the 4.5 μm band, 9541 of which include data in both bands. For our eventual comparison of behavior in the optical and infrared (e.g., Section 7), we performed further cuts and retained only IRAC light curves with at least 30 unsaturated points covering a minimum time span of 20 days. For the subset of observations executed in staring mode, we have additional high-cadence light curves for 290 objects with 3.6 μm band data and 259 objects with 4.5 μm band data. We have also produced binned staring light curves with smaller error bars at 2.5 minute cadence, which we describe below. Since the staring fields were observed simultaneously, there is no overlap between the 3.6 and 4.5 μm sets.

Comparison of the rms values for all mapping light curves with their mean predicted uncertainties revealed additional systematic errors that limited the photometric accuracy to 1%, as seen in Figure 4. In order to identify variable stars, we need to characterize the systematic effects as a function of magnitude in detail. Because of the transition to shorter exposure times toward the bright end of the sample, there is an upward bump in the rms distribution near 10th magnitude, and it is difficult to assess the systematic contribution here. We therefore computed a different measure of the systematic, "S," which should be independent of exposure time:

$$\begin{equation} S=\sqrt{{\rm rms }^2-\sigma ^2}, \end{equation} \tag{ 1 }$$

where σ is the uncertainty estimate described above, without systematics, and rms is the rms deviation of each light curve omitting flagged points. Plotting S against magnitude, we found that the nonvariable stars traced out an exponential trend, with a larger systematic contribution for fainter objects. To omit as many variable objects as possible while characterizing this effect, we considered only stars lacking firm NGC 2264 membership status (see the Appendix for membership evaluation criteria) and computed the median of the distribution of S values in 0.5 mag bins. We then fit an exponential curve of form

$$\begin{equation} S(m)=e^{(m-m_0)}+S_0, \end{equation} \tag{ 2 }$$

where m is magnitude and m₀ and S₀ are free parameters in units of magnitude.

Zoom In Zoom Out Reset image size

The distribution of S values, as well as the fits for each IRAC band, are shown in Figure 5. We obtained best-fitting values of S₀ = 0.014, m₀ = 19.75 in the 3.6 μm band and S₀ = 0.008, m₀ = 19.28 in the 4.5 μm band, indicating that the photometry includes systematics at the 1% level. Indeed, this is expected from intrapixel sensitivity variations and detector gain variations, both of which are known to result in position-dependent flux measurements (see the Appendix). As shown in Figure 5, S increases exponentially from ∼0.01 mag for 13th magnitude sources to 0.04 mag for 16th magnitude sources. Figure 4 demonstrates that the expected Poisson noise and sky noise increase from ∼0.01 mag to ∼0.1 mag over the same range of magnitudes. Thus, while both random and systematic uncertainties increase toward faint magnitudes, the systematic effects increase more slowly toward the faint end, such that they dominate the error budget for sources brighter than ∼14th magnitude, with random uncertainties becoming increasingly dominant for fainter sources.

Zoom In Zoom Out Reset image size

To fully account for the errors, we added the fitted magnitude-dependent systematic offsets in quadrature to the derived uncertainties and adopted the result for the selection of variables and production of variability statistics, as will be discussed in Section 6. In addition, we have performed similar assessments of the staring data. Since those light curves are not presented extensively here, we refer the reader to the appendix for a discussion of their quality and the mapping/staring data merging process.

NGC 2264 is the only rich 1–5 Myr cluster situated within one of the CoRoT eyes, making it the primary target for young star variability with this telescope. For the CSI 2264 campaign, we monitored a ∼13 × 13 field in NGC 2264 centered at R.A. = 06^h40^m180, decl. = +09°41'4624 (see Figure 1) from 2011 December 1 to 2012 January 9. These observations composed the fifth CoRoT short run ("SRa05"). All stars were placed on the second exoplanet channel CCD (E2), as the first channel ceased to function in early 2009. NGC 2264 was previously observed during the SRa01 short run in 2008, using the first exofield CCD, E1. In each case, targets consisted of confirmed and candidate NGC 2264 members, as well as field objects selected for CoRoT's transiting planet search program. Only data for pre-selected targets are downloaded from the satellite, with photometry consisting of flux values within a pre-defined aperture mask for each star. For SRa05, we monitored 489 confirmed and 1617 candidate NGC 2264 members, along with an additional 2129 field stars.

Light curves are produced by the CoRoT pipeline (Samadi et al. 2006), which performs corrections for gain and zero offset, jitter, and electromagnetic interference, as well as background subtraction. We obtained level N2 reduced data from the CoRoT archive, which consists of fully processed light curves. This includes fluxes and background levels, flagged for outliers and hot pixels. Typically ∼10%–15% of these were flagged and omitted from the light curves.

Exposure times were 512 s, resulting in up to ∼6300 data points per light curve. For a subset of stars with signs of eclipses or transits, a high-cadence mode was triggered with 32 s exposure times. This mode generated light curves with up to 100,850 points and was mainly used for stars brighter than 14th magnitude in R band.

The full range of magnitudes for observed stars was R ∼ 11–17. Light-curve rms values ranged from 0.00055 for the brightest stars to 0.01–0.1 for the faintest objects. There is a substantial spread in rms as a function of magnitude, because of strong systematics in a subset of the light curves, which we address below.

We produced light curves from the CoRoT data by converting raw fluxes to the magnitude scale and subtracting the median after removal of outliers automatically flagged by the CoRoT pipeline. The zero point for CoRoT photometry is not well determined and may vary slightly between runs. To estimate a calibration for the SRa05 data set, we compared the logarithmic mean flux with R-band photometry reported in the literature by Rebull et al. (2002), Lamm et al. (2004), and Sung et al. (2008). A fit of the intercept results in an R-band zero point of 26.74.

The full 2011 CoRoT data set on NGC 2264 contains light curves for a total of 4235 objects, just over half of which are field stars, based on the membership criteria outlined in the Appendix. These were included for the transit monitoring portion of the campaign. A total of 2500 of the CoRoT targets were previously observed during the SRa01 run in 2008. In total, 1303 stars were observed by both CoRoT and Spitzer, about 65% of which are possible or likely field stars. All CoRoT light curves encompassed the full 39 day campaign.

In general, we performed source cross matching with a 1'' radius. For targets with source confusion within the CoRoT mask, we first identified the object in the optical by requiring its CoRoT calibrated R magnitude to match photometry of a known source to within 0.5 mag. We then identified the 2MASS counterpart and used those coordinates to select the appropriate Spitzer source.

CoRoT light curves include a number of well-known systematic effects, including outliers and flux discontinuities due to radiation events or changes in detector temperature (Auvergne et al. 2009). Isolated outliers are automatically detected by the CoRoT pipeline and flagged in the resulting photometry. Other types of systematics can be removed if the light curves are well behaved or if color data are available (e.g., Mislis et al. 2010). However, correction becomes more difficult in a highly variable, mostly monochromatic data set such as ours.

We found that ∼10% of light curves contained abrupt jumps in flux not attributable to stellar variation. Two prominent discontinuities appeared at the same time stamps in many light curves. These are due to detector temperature jumps, the times of which were provided by the CoRoT team. We searched for and corrected discontinuities at these locations by computing a weighted difference between the flux difference on either side of each jump.

Background correction has also introduced time-dependent systematics into the light curves. Background levels are measured as a median across 400 10 × 10 pixel windows placed in star-free regions of the FOV. As the E2 CCD has aged, the level of dark current, as well as its gradient across the detector, has increased. Further complicating background measurements is the intrinsic peak in background due to nebulosity near the NGC 2264 cluster center. As a result, the median background value subtracted from stellar fluxes can be an underestimate for many sources. The combination of these effects, along with the different CCD used to observe NGC 2264 during 2008, has resulted in artificial variability amplitude changes for stars observed in both the SRa01 and SRa05 runs. For well-characterized variables such as eclipsing binaries, we note amplitude differences of up to 10%.

To omit problematic data from variability statistics, we visually inspected all 2011 CoRoT light curves and flagged those with discontinuities that were at least a factor of seven times the rms (and often much more). About 10% of light curves received flags.

To assess the noise levels and fit the rms distribution as a function of magnitude, we have plotted the rms values of the entire ∼5000-object CoRoT data set versus the magnitude values, using the zero point of 26.74; the result is displayed in Figure 6. For bright objects (R < 14), the trend of rms values follows a slope consistent with Poisson noise plus a constant systematic. The distribution expands to lower values for objects beyond R = 14, with some rms values reaching the Poisson limit. This structure in the diagram is related to the change from high to low cadence near R = 14 and is not explained by the flagged objects, as these have a large range of rms values at all magnitudes. As seen in the figure, the rms values of the final light curves ranged from 0.001 mag (at R ∼ 11.5) to 0.04 mag (at R ∼ 17), with larger amplitudes for variable objects.

Zoom In Zoom Out Reset image size

Prior optical and near-infrared YSO variability studies have proposed several different classes of T Tauri star behavior, depending on light-curve morphology, accretion status, and stellar mass. The first general survey of YSO variability was performed by Parenago (1954), who classified T Tauri stars based on histograms of their brightness levels. Objects that were more often bright than faint were referred to as Class I, objects with typical brightness similar to the mean were Class II, and stars that were more often faint than bright were called Class III. Following up on this, Weaver & Frank (1980) explored the relationship between light-curve class and emission line strength, finding that Class I objects were associated with preferentially weak strengths. Later, Herbst et al. (1994) presented a variability classification scheme based on 15 yr of photometric monitoring of 80 young stars. They divided light-curve behavior into four types, including periodic variations caused by cool spots in WTTSs ("type I"), irregular and periodic variations in CTTSs (types "II" and "IIp") attributed to variable accretion luminosity and rotational modulation of hot spots, respectively. In addition, they identified flux changes of 0.8–3 mag in the V band among T Tauri stars earlier than spectral type K1 ("type III"). The mechanism for this behavior was proposed to be obscuration by circumstellar material, since it was not associated with veiling changes. It frequently involves pronounced flux decreases, interpreted as extinction events as in the prototype UX Ori.

A similar scheme was suggested by Alencar et al. (2010) based on higher cadence light curves of T Tauri stars from the first CoRoT short run on NGC 2264 in 2008 ("SRa01"). They divided the variability sample into periodic sinusoidal (spot-like), periodic flat topped (AA Tau-like), and aperiodic (irregular). These correspond roughly to the Herbst et al. (1994) classes, except that the AA Tau-type objects have much shorter timescales (a few to 10 days) than the higher mass UX Ori stars.

Further variability studies (e.g., Carpenter et al. 2001; Lamm et al. 2004; Grankin et al. 2007, 2008; Scholz et al. 2009) continued the strategy of dividing variability into periodic and aperiodic and assessing physical mechanisms by analyzing color changes as a function of magnitude. These photometric campaigns typically involved monitoring at lower cadences (days to years) or with many more gaps in time sampling than our own data set. In the case of NGC 2264, Lamm et al. (2004) classified stars only as periodic or irregular, based on periodograms and χ² tests. They concluded that variability was likely due to rotational modulation of cool spots (if periodic), or hot spots generated by accretion flows onto the star (if irregular).

One of the most common phenomena in our time series is transient optical fading events. These have been noted previously in YSOs, from the early-type HAeBe stars (e.g., Herbst & Shevchenko 1999) to lower mass AA Tau stars (e.g., Bouvier et al. 1999). Based on color trends and an inconsistent relationship of photometric behavior with veiling in accompanying spectra, these episodes have been attributed primarily to extinction by circumstellar material (e.g., Bouvier et al. 1999; Alencar et al. 2010). In some cases, the events recur periodically, although their depth changes. Morales-Calderón et al. (2011) identified similar behavior in infrared and near-infrared time series of YSOs in Orion, likewise attributing "dipper" events to stellar occultations by regions of enhanced optical depth in a surrounding dusty disk. In several cases, periodic fading on longer timescales (months to years) has been traced to obscuration by a circumbinary disk precessing about a double- or triple-star system (e.g., Cohen et al. 2003; Plavchan et al. 2008a; Herbst et al. 2010; Plavchan et al. 2013). If the dips are caused by varying dust extinction along the line of sight, then the light curves should become redder as they grow fainter. Bouvier et al. (2003) indeed found this behavior in the prototype AA Tau itself, noting color slopes similar to the values expected from an interstellar medium extinction curve.

Unsurprisingly, prominent fading events appear in our CoRoT and Spitzer light curves. These are distinguished from other types of variability by both the sharpness and rapidity of the brightness troughs. Outside of these events, the light curves display a "continuum" flux level with less variability. To quantitatively identify and compare dipper behavior in the different bands, we have selected by eye objects that are clearly asymmetric, in that their light curves appear different when flipped upside-down.

We have identified 35 optical and seven infrared dipper stars in the joint CoRoT/Spitzer disk-bearing sample. Two (Mon-000183 and Mon-000566) exhibit clear dips only in the infrared, although with knowledge of the Spitzer data, one can see small dips (2%–3% depth) in the CoRoT light curves at the same times. Many of the optical dippers, however, are accompanied by entirely different infrared light curves, with high-amplitude, long-timescale variability and little to no sign of dips.

Prominent example light curves are shown in Figures 9–12. In Section 6.2, we introduce criteria to distinguish between periodic, quasi-periodic, and aperiodic light curves. The dipper sample here is split nearly evenly among quasi-periodic and aperiodic behavior, although it is not clear how the physics variability mechanisms differ between the two groups. We find that seventeen of the optical dippers are quasi-periodic, compared to four in the infrared; we label them "QPD" in Table 4, and these are synonymous with the periodic AA Tau variables described in Alencar et al. (2010). While these may be explained by a warped disk or orbiting companion, it is more difficult to envision a mechanism for the aperiodic dippers. A possible origin is that magnetic turbulence induces scale height changes in the inner disk, as envisioned by Turner (2013).

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

Also problematic is the total fraction of our data set displaying dipping behavior. With more than 20% of the sample in this category, the assumption that dust obscuration is happening in randomly aligned disks implies that the occulting material lies at scale heights of ∼0.2, in units of distance from the central star. The presence of material this high up in a hydrostatic disk might require large surface densities at this radius.

In several cases, dipper behavior does not persist over the entirety of the time series. The fact that we detect fewer dips in the infrared is not simply a selection bias due to the lower cadence of the IRAC time series; a number of objects display prominent dips only in the optical, whereas any corresponding dipping behavior at longer wavelengths is hidden by high-amplitude, long-timescale behavior originating in the disk.

In general, dipper light curves display infrared amplitudes that are less than the optical amplitudes; this supports the idea that the fading events are caused by enhanced dust extinction. We explore the degree of correlation between optical and infrared flux in dipper light curves, as well as implications for the underlying physics, in Section 7.

Dipper light curves are not the only ones displaying asymmetric behavior with respect to flux. We have identified a separate set of variables in the CoRoT data set that exhibit abrupt (i.e., 0.1–1 day) increases in flux, followed by decreases in flux on similar timescales. Examples are shown in Figures 13 and 14. We refer to these events as short-duration bursts (hereafter "bursts" or "B" in Table 4). We wish to distinguish them from more powerful outbursts such as FU Ori events, as well as the longer bursts identified in young stars by Findeisen et al. (2013). The bursts in our data set can be differentiated from coronal flares by their symmetric shapes (with respect to time), generally much longer durations, and sequential nature (i.e., many bursts occur in a row; they are by no means isolated events). Our examination of the light curves of WTTSs in the CoRoT sample revealed at most a few coronal flares per star over the 40 day observing window. Thus, while some of our bursting events may be due to stellar activity, most require a different explanation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Some YSOs exhibit isolated bursts above a "continuum" flux level, while others display similar behavior superimposed on a longer term trend. The baseline flux level in bursting light curves of some members of this class even displays quasi-periodic variability (i.e., "QPB" in Table 4). Since the burst durations are typically less than 1 day, we have in some cases subtracted out underlying trends to see the bursts more clearly. It should be noted that we are consequently insensitive to bursting behavior in timescales of ∼10 days or longer, unless there is no long-term trend.

We detect a total of 20 cases of short-duration optical bursting behavior in our 162 disk-bearing CoRoT/Spitzer data set, four of which also occur in the infrared (see Table 4). In addition, we find three examples (Mon-000119, Mon-000346, Mon-000185) of bursting in the IRAC data set that appear stochastic in the corresponding CoRoT light curves (as discussed in Section 5.5), as well as one case where the object does not display prominent variability at all in the optical (Mon-000273).

Bursting behavior in the Spitzer data set is much more difficult to identify with the sparser time sampling of ∼12 points per day in the mapping observations. We only identify a few examples exclusively in the IRAC time series but note that in about 50% of cases, there are increases in infrared flux coincident with strong bursts in the optical. It therefore appears that a combination of short cadences and high precision (∼1% or better) is crucial to detecting the full range of timescales involved in bursting, at least as it appears in our data set. Longer timescale analogs may be present in multi-year data sets, such as the one presented by Findeisen et al. (2013).

We believe that the bursts are caused by accretion instabilities, based on their resemblance to the predictions of Romanova et al. (2008) and strong ultraviolet (UV) excesses evident for many of these stars. To our knowledge, these are the first observations of an entire class of objects with such behavior. As such, a more detailed exploration of the bursters, including average durations and strengths, is presented in a companion paper (Stauffer et al. 2014).

We wish to note that while this paper and Stauffer et al. (2014) are closely linked, they were written in parallel and evolved somewhat independently. The two papers used slightly different sets of data—most importantly, Stauffer et al. (2014) included stars with only 2008 CoRoT light curves, whereas such stars were excluded here—and slightly different criteria for defining light-curve classes. This resulted in slightly different sets of stars belonging to each class. In particular, the burst-dominated class in Table 1 of Stauffer et al. (2014) included 23 stars; 19 of those are also listed as burst dominated in Table 4 of this paper. Of the remaining four, one (Mon-000185) was included in Stauffer et al. (2014) based on its 2008 CoRoT light curve and hence was not in the parent sample for this paper. The other three were all classified as stochastic in our Table 4. Inclusion or exclusion of these stars from the burst-dominated class would not appreciably change the conclusions in either paper.

Our formal period search (Section 6.2) uncovered many quasi-periodic light curves in both the CoRoT and Spitzer data sets. Only five optical and four infrared light curves in our infrared-excess restricted data set displayed strictly periodic behavior typical of the spotted WTTSs in the current and previous CoRoT run on NGC 2264 (Affer et al. 2013). This set includes only two stars that are periodic in both bands (Mon-001085 and Mon-001205; see Table 4). This is not entirely surprising considering that our stars were selected specifically to have disks based on their SEDS, and many are actively accreting based on spectroscopy.

The remaining quasi-periodic disk-bearing stars exhibit repeating patterns that either change in shape from cycle to cycle or include lower amplitude stochastic behavior superimposed on a periodicity. We refer to these as quasi-periodic symmetric, or "QPS," in Table 4. Excluding the quasi-periodic dipper class, as well as the handful of quasi-periodic bursters already accounted for above, we detect 27 quasi-periodic symmetric cases in the optical and 22 in the infrared, with 6 appearing in both bands. Examples are shown in Figures 15 and 16. Unlike the bursters and dippers, all of these light curves are relatively symmetric with respect to an upside-down flip. We consider here only variables with periods less than half the total observation baseline; longer timescale behavior is indistinguishable from the "stochastic" class, introduced below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The infrared variables that show repeating patterns do tend to be partially correlated with the optical behavior, even in cases for which we have assigned different morphology classes for the two wavelength bands. This may reflect multiple origins for the infrared flux variations: if a light curve is a superposition of periodic behavior and higher amplitude aperiodic changes, then we may detect the former but assign a different overall classification.

We propose two possible origins for this quasi-periodic symmetric variability. First, it may be a combination of purely periodic variation (such as cool starspot modulation) with longer timescale aperiodic changes (e.g., from accretion). Second, it could reflect a single variability process that is not entirely stable from cycle to cycle. Examples include stellar hot spots for which brightness evolves as a function of stochastic accretion flow (e.g., the type IIp variability highlighted by Herbst et al. 1994), or structural inhomogeneities in the surrounding disk periodically occulting the central star, but with geometries different from the standard dipper orientation. The latter is indicated in cases where a vsini value suggests that the measured period is too long to originate at the stellar photosphere (e.g., Artemenko et al. 2012; Cody et al. 2013). In addition, we cannot rule out that we are seeing cool spots that evolve much faster on disk-bearing stars than on their weak-lined counterparts.

The majority of remaining variables show no preference for fading or brightening but nevertheless exhibit prominent brightness changes on a variety of timescales. Fourier transform periodograms (see Cody & Hillenbrand 2010) for these objects display amplitudes devoid of dominant periodicities and consistent with a "1/f" trend in amplitude; two examples are provided in Figure 17. This is in contrast to typical red or "flicker" noise, which follows a 1/$\sqrt{f}$ trend in amplitude and 1/f in power (Press 1978). Thus, our stochastic objects display more power at low frequency than expected for standard flickering. Similar behavior was detected in a Herbig Ae star by Rucinski et al. (2010). The combination of lack of periodicity and coherence of the light curves suggests that these objects may represent a superposition of variable extinction (as in the dipper stars) and stochastic accretion (as in the bursting stars).

Zoom In Zoom Out Reset image size

Stochastic behavior constitutes one of the largest classes of variability in the CoRoT data set. Figure 18 illustrates prototypical examples of optically stochastic behavior. We focus specifically on short-timescale stochastic YSOs, which have peak-to-peak timescales much less than the 40 day duration of the light curve. There are in fact only two longer duration variables among the optical sample (see Section 6.5); thus, there is a true dearth of optical behavior on timescales longer than a week.

Zoom In Zoom Out Reset image size

Most infrared light curves are aperiodic, with fewer examples of quasi-periodic or asymmetric flux behavior (as in dippers or bursters) than in the optical. The timescales evident in the mapping data are also somewhat longer. Since the flux is generally disk dominated (although not by a large factor) at 3.6 and 4.5 μm, stochastic infrared variability could have a different origin from that in the optical, although it may simply reflect reprocessed variable starlight. The latter possibility is supported by the fact that we find stochastic behavior to be more than twice as common in the optical band than in the infrared. For long-timescale trends (see Section 5.6), the reverse is true: the behavior is much more common in the infrared. Orbital timescales in the disk emission region range from a couple of days to several weeks, depending on the mass and temperature of the star. It is therefore difficult to produce stochastic infrared behavior without invoking reprocessed optical variability or long timescales, which we consider separately below.

We present a selection of stochastic infrared light curves in Figure 19.

Zoom In Zoom Out Reset image size

A subset of our disk-bearing stars displays variability that grows or declines all the way out to the longest timescale of observation. We select as "long timescale" those objects for which the largest peak-to-peak amplitudes occur on timescales of 15 days or longer. This behavior is much more common among the infrared light curves than those in the optical. Because of the 30–40 day duration of our time series, we are unable to classify the trends morphologically and simply refer to them as "long-timescale" variables. The long duration of flux changes likely reflects disk dynamics beyond the inner edge, where Keplerian timescales are longer than one week. It also involves fairly large amplitude changes, from 0.08 to 0.6 mag among the sample identified here. We present prototypical examples of long-timescale infrared behavior in Figure 20.

Zoom In Zoom Out Reset image size

Notably, the optical behavior of these objects involves much shorter timescales. Only two CoRoT light curves contained trends gradual enough to be classified as "long timescale." The rest fall into a variety of classes, from quasi-periodic dippers to stochastic.

The distribution of stochastic and long-timescale stars in each band suggests that while some infrared variability reflects reprocessed starlight, this is a minor contribution to the flux. Some other process must dominate the long-wavelength variability at long timescales, and this could involve changes in the disk luminosity or shape.

There are several notable cases of rare light-curve types that merit discussion, since they may signify unusual geometry, disk evolutionary states, or variability processes.

We detect purely periodic optical behavior in only five stars (Mon-000335, Mon-001064, Mon-001085, Mon-001205, and Mon-000954), which we show in Figure 21. The latter two are nearly sinusoidal. These light curves are not easily distinguishable from the magnetically active WTTSs, and in fact all but two of them have weak or transitional disk SED classifications ("II/III"), along with weak reported Hα emission equivalent widths and low UV excesses. By some schemes, these sources could have been included with the WTTS set. We therefore believe that the variability here reflects cold spots on the stellar surface in systems where circumstellar dust obscuration and variable accretion are minimal. Nevertheless, it could be generated by very stable accretion flow onto a hot spot in the two objects that have significant infrared excesses (Mon-000335 and Mon-001064).

Zoom In Zoom Out Reset image size

Similar to the optical light curves, we detect only four instances of periodic infrared behavior, shown in Figure 22: Mon-000279, Mon-000427, Mon-001205, and Mon-001085. The latter two of these are also periodic in the optical, likely reflecting the variability induced by rotational modulation of starspots.

Zoom In Zoom Out Reset image size

We also highlight an unusual optical light curve consisting of a periodic morphology interrupted by fading episodes, as shown in Figure 23. The origin of the variability in this K5.5 star, Mon-000378, may be a combination of spots interspersed by circumstellar extinction dip events; this is difficult to verify since the infrared behavior does not correlate with the optical in this case.

Zoom In Zoom Out Reset image size

The fractions of different variability types offer insights into the diversity and relative importance of different physical mechanisms at work in the inner disks, magnetospheres, and stellar photospheres of young stars. We reproduce the schematic of the different variability types in Figure 24, with their fractions overlaid.

Zoom In Zoom Out Reset image size

The fractions for different types of optical variability are also presented in Table 5. A subset of light curves did not fall into any of the morphological classes described above, either because the data were noisy or because the light-curve shape was unusual. We denote these as unclassifiable in Tables 3 and 5.

Table 5. CSI 2264 Morphology Fractions in Disk-bearing Stars

Morphology Class	Optical	Infrared
Bursters	13$^{+3}_{-2}$% (21)	5$^{+2}_{-1}$% (8)
Periodic dippers	10.5$^{+3}_{-2}$% (17)	3$^{+2}_{-1}$% (5)
Aperiodic dippers	11$^{+3}_{-2}$% (18)	1$^{+2}_{-1}$% (2)
Quasi-periodic symmetric	17±3% (27)	13.5$^{+3}_{-2}$% (22)
Stochastic	13$^{+3}_{-2}$% (21)	6$^{+2}_{-1}$% (9)
Long timescale	1$^{+2}_{-1}$%	30$^{+4}_{-3}$% (48)
Periodic	3$^{+2}_{-1}$% (5)	2.5$^{+2}_{-1}$% (4)
Multiperiodic	1$^{+2}_{-1}$% (2)	1±1% (1)
Eclipsing binary	1±1% (1)	1±1% (1)
Unclassifiable	11$^{+3}_{-2}$% (18)	29$^{+4}_{-3}$% (47)
Non-variable	19±3% (30)	9$^{+3}_{-2}$% (15)

Notes. We list the fractions of objects in each variability class, along with the number (out of the 162 member disk-bearing object set) in parentheses. These classifications have been made by eye, whereas we provide statistical support for them in Section 6.4.

Download table as:  ASCIITypeset image

The non-variable stars are those sources that did not pass the Stetson test (infrared) or rms test (optical and single-band infrared). Among the CoRoT light curves, we are only able to distinguish variability down to rms values of 0.02 at R = 14 and 0.06 at R = 16. It is likely that the 19% of optical sources that are not variable according to our statistical measures simply have lower level brightness fluctuations.

Of note, the total fraction of dippers (both periodic and aperiodic) is just over 21%, a value that significantly exceeds the prediction of Bertout (2000) (15%), but is not as large as the ∼30% found by Alencar et al. (2010). These results suggest that the frequency of dipper behavior is surprisingly large, given typical inner disk scale height estimates (e.g., H/R ∼ 0.05).

Assessment of our infrared data set revealed that long-timescale variability was the most common type, and the other variability groups were less populated. The associated fractions for each infrared morphology class are also included in Table 5. Evidently, infrared variability is more common than optical variability, at least at the wavelengths and timescales sampled by these observations of disk-bearing sources.

In order to compare our identified variability types with models, we must first account for any biases introduced by the CoRoT sample selection. Targets for this telescope were selected in advance, whereas we observed all stars in a ∼1 square degree region with Spitzer/IRAC. CoRoT targets were selected primarily with NGC 2264 membership in mind, and priority was given to known CTTSs. Since we have already restricted our light-curve morphology sample using membership and infrared excess, there should be no associated bias. There was, however, knowledge of variability from CoRoT's SRa01 observations in 2008, and one might ask whether this contributed to more variables being observed than in a randomly selected sample of CTTSs. We argue that the answer is no, for the following reasons.

First, we have examined the input target priority list for the 2011 CoRoT observations. It is this set from which the final target list and pixel masks were selected. Of 772 possible high-priority targets, 350 were not previously observed by CoRoT and hence should have no variability bias. Of the 474 targets that were already observed in 2008, none had variability as the sole signature of youth. All input stars had prior membership information consistent with our criteria described in the Appendix and therefore should not be biased by previous knowledge of variability.

To probe further for systematics in our variability class fractions, we have selected a set of disk-bearing NGC 2264 members that were observed with CoRoT in 2008 and Spitzer in 2011. In this way, we assemble a sample similar to the 2011 data set highlighted in the rest of this paper, but with no prior knowledge of variability (apart from some sparser ground-based data, such as that of Lamm et al. (2004)). We classified the optical morphologies into the same categories described in this section and computed their fractions. We find that the fraction of non-variables is actually higher in 2011 than 2008 (19% versus 11%), again refuting the idea that variability knowledge biased the input catalog, thereby affecting the results reported here. On the contrary, we attribute the higher fraction of non-variables in 2011 to a slightly fainter sample. The total number of disk-bearing stars approximately doubled from the 2008 to the 2011 CoRoT run, and many of the newly observed sources were fainter and thus had lower signal-to-noise ratio than the typical SRa01 target.

Comparing the distribution of variable types in 2008 versus 2011, we find that in every one of the categories presented in Sections 5.1–5.6, except for the aperiodic dippers, the fraction of stars in each class for 2008 is the same to within the uncertainties presented above in Section 5.8. The fraction of aperiodic dippers decreased from ∼20% in the 2008 data set to 11% in the 2011 data set. This change is offset by the increase in the fraction of non-variable sources in 2011, leaving all of the other class fractions the same to within the 1σ errors. We speculate that the slightly fainter (and hence lower mass) 2011 sample may bias against the detection of aperiodic dipper behavior. However, it is not clear as to why the periodic dipper fraction is then similar in both samples. Because of this discrepancy in the aperiodic dipper fraction, we suggest that its uncertainty be increased from ∼3% to ∼10% in any theoretical attempt to account for the distribution of variability types.

6.1.1. IRAC Data Set

The high precision and cadence, as well as large sample size, of our Spitzer and CoRoT data sets permit detection of stellar variability down to fairly low amplitudes in comparison to previous photometric monitoring studies. Since our focus here is on disk-bearing stars, we might expect the majority of targets to display non-sinusoidal or aperiodic variability dominated by accretion and circumstellar effects rather than the cool spot rotational modulation that explains well the periodic variability of the WTTSs, as pointed out by Herbst et al. (1994). We have therefore taken several approaches to identifying photometrically variable objects in the photometry.

For IRAC data, we take advantage of the fact that nearly simultaneous light curves are available at both 3.6 and 4.5 μm for the majority of objects (151/162 stars in the disk-bearing CoRoT/IRAC sample). Furthermore, the behavior in these two bands should be similar, since their emission regions in the disk are at comparable radii. This enables identification of correlated variability via a Stetson cross-correlation index (Stetson 1996), the use of which was previously promoted for variability detection in young stars by Carpenter et al. (2001) and Plavchan et al. (2008b). For data in two bands, this index is defined in the following way:

$$\begin{equation} {\rm Stetson}=\Sigma _{i=1}^{N}{\rm sgn}(P_i)\sqrt{|P_i|}, \end{equation} \tag{ 3 }$$

where N is the number of simultaneous pairs of observations in the two bands, sgn is the sign function (i.e., −1 for negative input values, 0 at zero, and 1 for positive input values), and P_i is the product of the normalized residuals of two observations:

$$\begin{equation} P_i=\frac{N}{N-1}\left(\frac{d_{1i}-\langle d_1\rangle }{\sigma _{1i}}\right)\left(\frac{d_{2i}-\langle d_2\rangle }{\sigma _{2i}}\right), \end{equation} \tag{ 4 }$$

where d_1i and d_2i are simultaneous photometric points in bands 1 and 2, respectively, 〈d〉 refers to their means, and σ_1i, σ_2i their uncertainties, including systematics.

Since data taken at each wavelength are not exactly simultaneous, we interpolate one of the light curves onto the time stamps of the other. We find that the Stetson index has a roughly constant baseline regardless of stellar brightness. Objects with Stetson values well above this level have a degree of correlation in the two bands that cannot be accounted for by random noise. We determine a threshold for variability by examining the Stetson index distribution of likely field stars that do not meet any of the NGC 2264 membership criteria outlined in the Appendix, as displayed in Figure 25. These have a low probability of being variable and should therefore reflect the intrinsic noise spread in Stetson index. We fit a Gaussian profile to the distribution, finding a 1σ width of 0.07. We therefore adopt a variability threshold of 0.21, for 3σ confidence; this is denoted by the dotted line in Figure 25. The distribution of Stetson indices for likely cluster members (also shown in Figure 25) displays a break at this value, confirming that this cutoff is a reasonable dividing line between variable and non-variable objects.

Zoom In Zoom Out Reset image size

To select infrared variables among the set of targets with only single-band IRAC mapping or staring data, we must rely on another variability criterion. Light-curve rms is a suitable metric, as long as we take the systematic noise contribution into account. We incorporate the noise model determined in Section 2.2 to compare the measured rms values against the expected noise level in each mapping light curve. For staring light curves, we use the estimated uncertainties, since these are a good approximation to the errors. Each of the four staring time series is treated separately.

We determine a cutoff value in the difference between measured and expected log (rms). As with the assessment of systematic errors, we have divided our sample into sets of field stars and confident NGC 2264 members. Stars that are candidates but not confident members are not included in the determination of the cutoff, but they are evaluated for variability afterward. We take the distribution of rms values for the group of field stars to be representative of non-variable behavior, and find it to be approximately Gaussian. The excess of rms beyond expectations is much larger for the second group, as shown in Figure 26. We adopt as our single-band variability criterion a cutoff of 3σ, or 0.15, in the distribution of log (rms/rms_expected). This is indicated by dotted vertical lines in Figure 26. For stars with data in both IRAC bands, we find that the rms test is not as sensitive as the Stetson index in detecting variability.

Zoom In Zoom Out Reset image size

Out of the 162 member disk-bearing data set considered here, 11 do not have multi-band Spitzer data, and nine (82%) of these are variable by the rms criterion. In the subset that do have data in both bands, 137 of 151 (91%) are variable by the Stetson index. In contrast, only 36% of stars without disks detected by Spitzer are variable in the infrared, and many of these display sinusoidal patterns attributable to rotating starspots.

Objects with amplitudes larger than 0.2 mag were inspected on the raw IRAC images to determine whether stray light from a bright neighboring star might be inducing artificial variability via background contamination. This is often the case when an object with both 3.6 and 4.5 μm data is detected as variable based on its rms, but not based on the Stetson index. A handful of candidate variables were eliminated in this way.

6.1.2. CoRoT Data Set

Selection of variables in the CoRoT data set is more difficult than with Spitzer/IRAC given the systematic effects (e.g., jumps, hot pixels) present in some of the light curves. Defects remained even after correction of the two most prominent jumps in the light curves (see Section 3.2).

For variable selection, we first considered the full set of YSOs observed by CoRoT and then narrowed our focus to the 162 disk-bearing stars with simultaneous Spitzer data. We determined the median trend in rms as a function of magnitude, using 0.5 mag bins and omitting objects that were flagged or known NGC 2264 members (implying high likelihood of variability). Despite the unexplained structure at R > 14 (see Figure 6), the median fit in log space is well modeled by the Poisson noise expectation, shifted upward by a constant value of 0.35. We therefore adopted this as the underlying noise distribution. The fit breaks down for R < 12, and here we simply adopt the median trend in rms, which is roughly linear.

We require that variables have rms values at least three times the median noise level for their magnitude. In Figure 27, we display the distributions of log (rms/noise) for likely cluster members and field stars. The cluster members are clearly more variable, and the chosen cutoff at rms/noise = 3.0 selects variables at the ∼3σ confidence level, according to a fit of the non-member distribution. Additional low-amplitude periodic objects that did not meet this rms selection criterion were identified via periodogram and autocorrelation function (ACF) analysis, which we discuss below.

Zoom In Zoom Out Reset image size

We find the chosen threshold to be adequate in that all of the selected objects are clearly variable by eye. However, there is a collection of faint (R > 16) objects in the CoRoT sample that are also clearly variable by eye, but do not meet the selection criterion since their rms values are <0.05 mag. The brightness fluctuations in these cases consist of short-duration upward or downward events that depart significantly from the median value but are too transient to contribute significantly to the overall variance. The coherence of the light curves during the events causes them to stand out in comparison to the artificial fluctuations caused by CoRoT systematics. In a number of cases, this variability is supported by strong correlation with infrared behavior in the corresponding IRAC light curve. We identified a handful of additional such variables by eye, as indicated in Table 4.

Disregarding flagged light curves, we find 218 definitively variable stars in the entire CoRoT data set, circled in Figure 6. Although the WTTSs are not the subject of this paper, we find that ∼45% of them are variable by the rms criterion alone (this fraction would increase were periodicity detection to be employed as well). Considering hereafter only the ones in the 162-object disk-bearing sample with both IRAC and CoRoT data, we are left with 83 variable by the rms criteria. A further 30 faint stars mainly with R > 16 were added to this set based on variability identified by eye. We identified 19 additional variables based on periodicity or quasi-periodicity alone, as explained below. The variability status of all objects in the disk-bearing data set is indicated in Table 4.

Periodic photometric behavior is common among YSOs, although it is mostly associated with spotted WTTSs. Our disk-dominated sample does include a number of stars with quasi-periodic light curves, few of which are sinusoidal. We have developed a criterion below to differentiate between "periodic" and "quasi-periodic" stars. By "quasi–periodic," we refer to light curves that have a stable period but whose shape and/or amplitude changes from one cycle to the next. We use "periodic," on the other hand, to denote stable repeating patterns with shapes that evolve minimally over the 40 days of observation.

To select these objects from the IRAC and CoRoT data sets and identify legitimate variables that did not meet the rms or Stetson criteria, we have developed a period search technique. McQuillan et al. (2013) showed that the ACF is a particularly good tool for selecting the correct period from a non-sinusoidal light curve, by considering the larger of the first two local ACF peaks. The commonly used periodogram, on the other hand, tends to display many peaks corresponding to harmonics and aliases that may be confused with the true signal. We have therefore carried out a preliminary period search by interpolating all light-curve magnitude values onto evenly spaced time grids with interval Δτ and N points, and computing the ACF based on the following equation:

$$\begin{equation} {\rm ACF}(\tau)=\frac{\sum _{i=0}^{N-\tau /(\Delta \tau)-1}(d_i-\langle d \rangle)(d_{i+\tau /(\Delta \tau)}-\langle d\rangle)}{\sum _{i=0}^{N-1}(d_i-\langle d\rangle)^2}. \end{equation} \tag{ 5 }$$

Here d_i are the light-curve data points, 〈d〉 is their mean, τ = nΔτ is the time lag, and N is the total number of points in the interpolated light curve. We let the time lag run from zero to the maximum baseline of the time series; peaks in the ACF indicate lag values for which the light curve is self-correlated. While interpolation may alter the light curves slightly, the sampling is dense enough compared to variability timescales that we expect any resulting inconsistencies to be small. We typically oversample the light curve by a factor of 1.5, corresponding to Δτ ∼ 6 minutes for CoRoT data and Δτ ∼ 1.5 hr for IRAC data. For each ACF, we note all local maxima occurring at time lags greater than zero and less than half the total time baseline (i.e., ∼15 days for IRAC light curves and ∼20 days for CoRoT light curves). We required the amplitude of any such peaks to be greater than 0.05 over the surrounding local minima. We then select the first or second local maximum, depending on which is higher. An example of this process and the succeeding steps are illustrated for Mon-000660 in Figure 28.

Zoom In Zoom Out Reset image size

To check whether the selected ACF peak corresponds to a significant periodicity in the light curve, we then compute a Fourier transform periodogram and search for peaks within 15% of the frequency expected from the period of the ACF peak. Upon identifying the periodogram peak, we phase-fold the light curve around the selected period. Based on the coherence of phased light curves, we find that periods extracted from the periodogram are more accurate than those adopted from the ACF; the latter is more sensitive to long-term trends in the light curve. Initial use of the ACF is nevertheless vital to determining which of multiple peaks is the correct period to phase around, as shown by McQuillan et al. (2013).

Once the light curve is folded, we generate a smoothed phase-folded light curve smoothing over a boxcar with width 25% of the period. This approach is similar to that of Plavchan et al. (2008b), except that we only phase the light curve to a single period, as opposed to a continuum of periods. We overlay the smoothed phase curve on the original light curve, repeating it once per period. Comparison of the phase trend curve with the raw light curve provides an impression of how well a periodic model explains the behavior. We subtract the two curves to produce a residual as a function of time. For strictly periodic light curves, the remaining points should consist of noise, and indeed the residuals are consistent with the uncertainties shown in Figures 2 and 6. However, most of our light curves are better described as quasi-periodic: the amplitude of the residuals is significantly reduced compared to the raw data, but there remain strong trends not attributable to systematic errors. This is particularly the case for objects that display repeating flux dips of varying amplitude.

We have adopted a metric to assess the degree of periodicity in the light curves, by comparing the rms value before and after subtraction of the smooth phased curve. Since we are only considering timescales up to half the data length, we first remove any long-term trends on timescales over 15 days (IRAC) or 20 days (CoRoT) by subtracting a boxcar-smoothed version with a window of 10 days. We compute a periodicity metric, "Q," by assessing how close the light-curve points are to the systematic noise floor before and after the phased trend is subtracted from the light curve:

$$\begin{equation} Q= \frac{\left({\rm rms}_{\rm resid}^2-\sigma ^2\right)}{\left({\rm rms}_{\rm raw}^2-\sigma ^2\right)}, \end{equation} \tag{ 6 }$$

where rms_raw and rms_resid are the rms values of the raw light curve and the phase-subtracted light curve, respectively, whereas σ is the estimated uncertainty including the systematics (e.g., Section 3.3). Testing on sinusoidal light curves typical of WTTSs, we find Q to be a few percent or even negative (i.e., when the uncertainty is an overestimate). However, for light curves that appear to contain multiple sources of variability, the value is larger, ∼0.15–0.60. Light curves with no detectable periodicity have Q values of ∼0.6–1. Examples of different Q values are shown in Figure 29. We will test the dependence of Q on parameters such as period and time sampling in future work.

Zoom In Zoom Out Reset image size

As an independent check on the periodicities, we note whether the corresponding peak in the Fourier transform periodogram exceeds the local noise level by at least 4.0. This criterion was put forth by Breger et al. (1993) as indicating a 0.1% false alarm probability; see Cody & Hillenbrand (2010) for further discussion. All of the periodic sources and many of the quasi-periodic sources have significant periodogram detections. A more extensive analysis of the periodic sources in the 2011 CoRoT sample will be presented in a forthcoming paper in the context of an updated rotation rate analysis for NGC 2264.

Once a significant periodic or quasi-periodic behavior has been identified and the phased trend subtracted out of the light curve, we recompute the periodogram to assess whether further periodicities are present. We once again use the criterion of periodogram local signal-to-noise ratio greater than 4.0. We mask out all frequencies associated with harmonics of the first detected signal, as well as the aliases. The frequencies of aliases are determined by overlaying on the primary peak the window function associated with the periodogram. If one or more significant signals remain outside of the masked values, we then repeat the period search as outlined above.

We report the results of our period search in Table 4, noting in the morphology column which objects are significantly periodic, quasi-periodic, and multi-periodic. We identify two optical light curves with multiple periodicities (Mon-000434 and Mon-000164), as well as one infrared light curve with two periods, only one of which is observed in the optical (Mon-001181). These could indicate binary systems in which both stars show spot modulation.

Among the stars that are periodic or quasi-periodic, the distributions of periods in the optical and in the infrared both peak near 5 days, but we find a significant dearth of infrared periods beyond 9 days compared to a steady decline in the number of optical periods out to 15 days. Of note, some of the infrared periods may originate at the stellar surface, if stellar emission dominates the disk flux at these wavelengths (as is the case for very weak disks).

In addition to the cases of quasi-periodicity, we observe that many of the light curves are asymmetric with respect to a reflection along the magnitude axis. Some stars have prominent downward flux dips, while others have abrupt increases (see Section 5). We believe that this behavior is connected with the physical mechanisms causing variability. To quantify the degree of flux asymmetry, we have developed a metric, M. To determine its value, we first select the 10% highest and 10% lowest magnitude values in each light curve, after removal of 5σ outliers. This process is carried out by first smoothing the light curve on 2 hr timescales (CoRoT) or 6 hr timescales (for the more sparsely sampled IRAC data) and subtracting the smoothed trend from the raw data. Outliers are then measured on the residual light curve. After their removal, we compute the mean of the remaining points and compare with the median of the entire outlier-filtered light curve. We define the asymmetry metric via

$$\begin{equation} M=(\langle d_{\rm 10\%}\rangle -d_{\rm med})/\sigma _{d}, \end{equation} \tag{ 7 }$$

where $\langle d_{\rm 10\%}\rangle$ is the mean of all data at the top and bottom decile of light curve, d_med is the median of the entire light curve, and σ_d is its overall rms.

In some cases, there is a clear asymmetry in the light curve, but it is superimposed on a longer timescale trend. Our asymmetry metric is only sensitive to asymmetries on timescales less than about half the light-curve duration, or 15–20 days. We therefore remove trends on longer timescales by subtracting out a smoothed version of the light curve before computing M. We find that a smoothing window of a few days is sufficient; in cases where short-term light-curve features were oversubtracted, we retain the M value computed on the raw data. These cases were identified visually comparing the raw light curve with the long-term trend overplotted. All light curves were also visually inspected to determine which value was most appropriate.

We encounter a range of M values in our disk-bearing data set, from approximately −1 (prominent upward flux peaks) to just over 1 (flux dips). Examples are shown in Figure 30. Examining the light curves by eye, we observe the most obvious asymmetric behavior to occur for values of M greater than 0.25, or less than −0.25.

Zoom In Zoom Out Reset image size

The M metric, in combination with timescale (Section 6.5) and the periodicity measure, Q, enables us to quantitatively retrieve the morphology classes that were first established by eye and that presumably represent different variability mechanisms. We plot Q against M in Figure 31, along with suggested boundaries to divide the different variability types. Examining the classifications made by eye (color-coded points), we find that the selected boundaries in Q (0.11, 0.61) and M (±0.25) are quite successful in separating the variability sample into classes. They appear slightly less useful in the infrared since many more sources in this band vary on long timescales and relatively fewer exhibit dipping or bursting behavior. The lower infrared cadence may play a role here as well. Nevertheless, the fact that two parameters can divide our sample so well into the predetermined groups is promising for future variability classification efforts based on sparser data.

Zoom In Zoom Out Reset image size

The Q–M diagram reveals several facets of the variability in our sources. First, while we have selected boundaries between different classes, the statistics show that there is a continuum of light-curve behavior along both the Q and the M axis. This suggests that sources on the boundaries of variability classes may be characterized by multiple physical mechanisms. Second, the Q–M diagram displays several areas devoid of points. There are relatively few variables that are both bursting and periodic. This may point to the stochastic nature of accretion. Likewise, there are relatively few dipper objects with highly periodic variability. The highest amplitudes of these light curves correlate with unstable variability patterns, such that they are quasi-periodic but not periodic.

Timescale is an important quantifiable aspect of our light curves. Visual inspection of the data reveals that some objects oscillate quite rapidly (i.e., multiple zero crossings in a week), while others show only long-term trends. We wish to define a timescale measure for the purposes of examining correlations with physical parameters. While it is easy to attach a value to the periodic and quasi-periodic light curves, a timescale for aperiodic stars is not so obvious. Tools such as the ACF and periodogram have traditionally been used to assess whether there is a characteristic timescale for aperiodic variability. However, it is not immediately clear what amplitude level one should choose to extract a timescale from the ACF, or how to determine whether there is such a special time at all. Furthermore, the time sampling and baseline may strongly influence the appearance of the ACF.

We have adopted a different strategy, which is to identify the median timescale separating the largest consecutive peaks in a light curve. This process consists of two steps: (1) for each light curve, we calculated how the inferred timescale varies as a function of the amplitude threshold applied in selecting peaks to calculate timescales between; and (2) we collapse this set of timescales into a single value representing the highest amplitude variations within the light curve. For the first step, we calculate timescales for amplitudes as small as the light curve's noise level and as large as the light curve's full range (i.e., the absolute maximum minus the absolute minimum, after filtering for outliers). For each of these values, we run through the light curve, marking maxima and minima that differ from surrounding peaks by more than that amplitude (see the top and middle panels of Figure 32). We start by selecting the absolute maximum and count all peaks preceding it. We then repeat the process to identify peaks succeeding it. From this "PeakFind" algorithm, we estimate a characteristic timescale for each amplitude by computing the median time difference of all consecutive pairs of peaks and then multiplying by two. The normalization by 2.0 ensures that the derived timescale will match the period for periodic objects. This process results in a trend of timescale versus amplitude, which we invert to amplitude as a function of timescale, as shown in the bottom panel of Figure 32. We note that this function is not necessarily monotonic, because of the smaller number statistics in the case of large amplitudes and small numbers of peaks.

Zoom In Zoom Out Reset image size

For many of our light curves, the PeakFind timescale reaches a maximum value that is much shorter than the duration of the light curve. This implies that the full variability amplitude is accounted for by a relatively short timescale phenomenon. To extract a single timescale from the trend—the second step in the process—we note the maximum peak-to-peak amplitude of the light curve and adopt the timescale corresponding to 70% of this (i.e., the vertical lines in Figure 32). If this timescale is comparable to the duration of the time series, or ∼30–40 days, then it is a lower limit on the true maximum variability timescale. We emphasize that this "70% amplitude" timescale is not necessarily a characteristic timescale, since variability may be produced on a continuous spectrum of timescales, and we are not sensitive to those longer than 30–40 days. However, for fairly rapid variability we consider it an approximate upper limit to the timescales on which variability is generated. Simulations using damped random walks show that, for aperiodic signals with characteristic timescales of 0.1–5 days and amplitudes of 0.2 mag or more, the peak-finding timescale is well correlated with the true timescale on average but shows scatter comparable to the true timescale for any individual source (K. Findeisen 2014, in preparation). Therefore, the timescales are best interpreted in an ensemble sense. Using them as guidance, we are able to distinguish between aperiodic light curves that oscillate on day to week timescales and those that wander up and down over the course of a month or more. We find that the distribution of aperiodic stochastic timescales from the PeakFind method centers around 5–10 days, indicating that short timescales are dominant in the optical, at least among the spectral types encompassed by our data set.

With two types of timescales (i.e., quasi-periodic and aperiodic) in hand for the variables in our sample, we can compare them with the rms values in the optical and infrared, as shown in Figure 33. Separating objects into their variability type (aperiodic or quasi-periodic), we do not find any obvious correlations, apart from a subtle rise in variability amplitude from 1 to 5 day timescales. Aperiodic infrared variables show higher amplitudes than the periodic varieties, with many of the former achieving rms values between 0.1 and 0.2 mag. Both optical and infrared variables achieve their highest quasi-periodic amplitudes on timescales between 5 and 10 days. The amplitudes of aperiodic light curves have a much wider distribution of amplitude versus timescale.

Zoom In Zoom Out Reset image size

To measure in detail the degree of correlation between optical and infrared time series, we first interpolate each CoRoT light curve onto the same time stamps as the IRAC 3.6 and 4.5 μm mapping data. Since the CoRoT observations were obtained at such high cadence, the effect of interpolation on the variability is negligible. The result is pairs of optical and infrared light curves with typically ∼300 points.

We have performed a cross-correlation of the CoRoT/IRAC sets using the Stetson index (Stetson 1996), as well as subtraction to generate color trends. We find that many sets of light curves are well correlated for part of the time series but less correlated or varying with a different color slope in other parts. This makes it difficult to quantify correlation via a single parameter.

As an alternative, we have computed a time-dependent Stetson index by comparing 20 point sections (1.7 days) of the IRAC and interpolated CoRoT light curves. The 20 point step size was motivated by the typical duration of short-timescale light-curve fluctuations, which is a few days. Well-correlated light curves should display a positive Stetson index with little dependence on time. Light curves that are correlated at some times but not others will instead display a fluctuating Stetson index, reaching large values at some times and dropping to zero or below at others. To differentiate this behavior from chance correlation between sections of two light curves, we have performed simulations using time series from different stars. We match optical and infrared light curves from different stars, thereby guaranteeing that there should be no true correlation. Only light curves with at least 250 points are allowed, since the more sparsely sampled HDR mode light curves will have a broader Stetson index distribution. One thousand simulations of randomly matched light-curve pairs revealed that 95% of uncorrelated light curves have a median running Stetson index between −0.5 and 0.5 when sampled at every 20 points. Repeating the simulations on HDR light curves for every five points, the value increased to 1.0. We also ran the simulations for a sampling step size of 50 points, but found that the running Stetson index showed more systematics, and its median was less consistent between the two IRAC channels, despite flux behavior being nearly the same in these bands. The distribution of median running Stetson index for CoRoT light curves paired with IRAC light curves is shown in Figure 35.

Zoom In Zoom Out Reset image size

We report the median running Stetson index for all objects in Table 4. Stars with large negative running Stetson indices exhibit anticorrelated behavior in the optical and infrared.

The fractions of correlated versus uncorrelated and anticorrelated light curves may indicate the percentage of sources in which variability is disk dominated. To determine this number, we assemble the distribution of median running Stetson indices in Figure 35, along with the distribution expected for uncorrelated light curves, as derived in our simulations. Although we have drawn a cutoff of 95% confidence, it is clear that the set of Stetson indices is skewed with respect to the distribution defined by the simulations. We conclude that there is low-level optical/infrared correlation in many of the objects; this is not surprising since stellar variability will generally have a small infrared contribution due to emission from the long-wavelength tail of the star's SED. In general, however, the larger infrared contribution from the dust means that variability originating in the disk will dominate the light curve.

Overall, 38% of our sample shows evidence of correlation at the 2σ level, and 58% shows correlation at the 1σ level. There is generally a correlation between the shapes in the optical and those in the infrared, but in many cases it is quite weak (i.e., it consists of transient dips or bursts matching up). To further illuminate the relationship between optical and infrared behavior, we have plotted a correlation matrix comparing the morphologies assigned in each band. Figure 36 confirms a low-level optical/infrared correlation at best in most sources.

Zoom In Zoom Out Reset image size

We have identified light curves that are correlated in the optical and infrared by selecting those with a large median running Stetson index. Some of these light curves display partial correlation in that the infrared flux trends mimic those in the optical over part of the time series, but at other points there is no resemblance between the two bands.

There are two scenarios for correlation of behavior in the optical and infrared. First, if there is only one variability process affecting the star, then the light curves should exhibit a relatively simple wavelength dependence. Second, if the disk flux dominates in the infrared bands but variability reflects mainly reprocessed stellar light, then we would also expect to see correlated behavior, either in or out of phase. We test these scenarios in Section 8 by comparing the Stetson index to measures of disk to star flux.

We have identified three broad categories of morphological behavior among the well-correlated light curves: optical dippers, bursters and stochastic stars with similar infrared amplitude, and optical stochastic stars with lower amplitude infrared light curves. A selection of these is presented in Figure 37, while optical and infrared light curves for the entire 162 member disk-bearing data set are provided in the Appendix, Figure 49. Some of these have behavior that makes sense in the context of the extinction-dominated variability model. Here, fluctuations in the light curve are caused by changing amounts of dust blocking the star, and the optical/infrared correlation simply reflects the wavelength dependence of extinction, diluted by any flux from the inner disk. Confronting this model are the two objects Mon-000183 and Mon-000566, which display deeper dips in the infrared than in the optical. This behavior is not expected from any standard extinction law and may require unique geometry, such as occultations of the disk by itself (i.e., the front blocks flux from the back wall).

Zoom In Zoom Out Reset image size

The remaining set of well-correlated optical and infrared light curves displays very similar morphology and amplitude in each band. Half of these are members of the optical bursting class, and the rest are stochastic or unknown morphology. Strong optical/infrared correlation for these would be expected if flux changes are reflective of increases in accretion, the radiation from which is absorbed and reemitted by the inner disk.

Some 62% of the variable light curves in our sample display little correlation between the optical and infrared bands. In these cases, there are likely two or more variability mechanisms at work, connected separately with the star and the disk. Alternatively, a single variability mechanism could be at work if it only affects one wavelength region. Since the infrared flux is in most cases not dominated by stellar or accretion emission, disk-driven mechanisms are necessary to explain the often large-amplitude (∼0.1 mag) infrared modulation on week or longer timescales. One possible origin is magnetic turbulence, proposed by Turner (2013) to modulate the inner disk scale height and thereby alter the observed mid-infrared emission.

We have assembled a set of the most prominent uncorrelated optical and infrared behavior in Figure 38. The wide range of morphologies in both bands is evident. Particularly interesting examples of uncorrelated behavior occur in Mon-000185, Mon-000273, Mon-000357, Mon-000876, and Mon-000928, for which there is almost no optical variability, but high-amplitude (0.1–0.3 mag) infrared light-curve excursions on 5–10 day timescales. We have checked the individual images of these objects for erroneous cross matching of CoRoT and Spitzer sources, but in all cases there is no nearby star that could be a better match. This type of behavior occurs in a few percent of our disk-bearing stars.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The selection of objects with uncorrelated behavior is particularly helpful for investigating the distribution of infrared variability timescales, since we do not have to worry about much contamination from optical processes. We have measured the aperiodic infrared timescales of stars with median running Stetson indices less than 0.6. As illustrated in Figure 39, the distribution displays a clump around 5–15 days, although there are many additional objects with infrared timescales at 20 days and longer. The individual examples of rapid infrared changes, as well as the overall distribution, lend support to the idea that disk structural changes may be occurring close to or even faster than the local dynamical timescale in some cases. The latter is between a few days and a few weeks, depending on stellar mass and dust properties.

Zoom In Zoom Out Reset image size

A prediction of YSO variability models for edge-on systems (e.g., Kesseli et al. 2013) is that inner disk regions may receive non-uniform illumination from the central star, thereby reradiating infrared flux in a time-variable manner. An example is an accretion hot spot on the stellar surface, radiation from which interacts only with the region of the inner disk on that hemisphere of the star. In this case, we expect to observe the infrared emission 180 degrees out of phase with the hot spot signature, since we do not view the heated region of the disk wall until it has rotated behind the star.

Surprisingly, we encounter very few such examples of optical/infrared phase shifts in our time series. The only object with a clear inverse relationship between the two bands is Mon-001031. Mon-001132 displays a more subtle phase lag, but this shift does not appear to be constant across the entire observation time. We present the two light curves in Figure 40.

Zoom In Zoom Out Reset image size

The main stellar parameter of interest is the effective temperature, which serves as a proxy for mass. We have estimated this for the 112 stars in our sample that have available spectral types in Table 3. We then plotted the variability properties involving timescale and amplitude against temperature. One might expect infrared variability timescale to scale with temperature since it could reflect dynamics near the sublimation radius, which is in turn dependent on stellar mass. For completeness, we have separated the timescale measurements into quasi-periodic and aperiodic sets, as well as divided them into the two wavelength bands before comparing with temperatures. The correlation diagrams resulting from this exercise are shown in Figures 41. Timescale, whether periodic or aperiodic, does not show any prominent trends as a function of temperature. However, there is a lower envelope to the distribution of aperiodic infrared timescales versus temperature; this could be consistent with the orbital period at the sublimation setting the minimum timescale for disk variability. In addition, the implied lack of a global dependence of the variability timescale on mass does not necessarily mean that there is no relationship between these two parameters. We suspect that there are subclasses of variability that exhibit different behavior as a function of stellar mass. For example, the dipper objects show larger infrared amplitudes in comparison to the optical for preferentially later spectral types. Finally, there may be currently unobservable factors, such as disk inclination, that influence the variability properties.

Zoom In Zoom Out Reset image size

We have also carried out a similar comparison of variability amplitude versus temperature, as shown in Figure 42. Here, we find the surprising result that aperiodic variability amplitudes are significantly lower around the cool stars in our sample. This trend is the opposite of what one would expect from a detection bias, in which it is harder to detect low-amplitude variability around fainter stars. The effect is stronger in the optical, but also seems to appear in the infrared. This is in contrast to the periodic variability, which shows no appreciable mass dependence. The rms/temperature trend in aperiodic variables may reflect a correlation of magnetic field strength and/or configuration with mass. However, this may be contradicted by recent results from Zeeman–Doppler imaging that suggest that the dipole field component weakens with increasing mass (e.g., Gregory et al. 2012).

Zoom In Zoom Out Reset image size

Perhaps the most telling correlation we have noted is the connection between light-curve flux asymmetry (M) and Hα equivalent width. We plot these two parameters in Figure 43. While there is not a one-to-one correspondence between them, it is clear that the more negative M values correspond to stronger Hα emission. This finding supports the idea that the burster class of light curves represents the most strongly and unstably accreting stars in our data set (see also Stauffer et al. 2014 for confirmation of this idea).

Zoom In Zoom Out Reset image size

In addition to searching for global correlations between variability and stellar properties, we can also ask whether variability is well connected to disk properties. We first compared the rms light-curve value for infrared variables with the slope of the SED, α (as calculated in Section 2.1). The result, shown in Figure 44, does not display a strong dependence on morphology type, but there does appear to be a subtle lower envelope, with few high-amplitude variables at low α values. This suggests that disks with earlier classes can achieve higher levels of variability, a result that has also been borne out in other clusters (e.g., Morales-Calderón et al. 2009). It also makes sense given that the Class II/III objects in our sample have less dust to generate variability.

Zoom In Zoom Out Reset image size

We also suspect that the ratio of disk excess to stellar flux plays a significant role in light-curve morphology. As a proxy for this ratio, we have computed the K − [4.5] color of all sources in our sample. We compare this with the median running Stetson index in Figure 45. Surprisingly, there is no strong dependence of optical/infrared correlation on disk strength either. However, the light curves of a number of individual sources do make sense in the context of their SEDs. For example, several stars (e.g., Mon-001094) display well-correlated behavior at similar amplitudes in the two bands. We suspect in these cases that there is a single variability process at work, and it is associated with the stellar surface. The disk, on the other hand, is weak and has little emission compared to the stellar flux in all bands. We see that this is the case from the small K − [4.5] values of these objects. Estimates of their disk to stellar flux ratios in the infrared also support this idea, since the values are less than 10⁻³.

Zoom In Zoom Out Reset image size

There are other objects in our diagram for which the high degree of optical/infrared correlation is not expected to relate strongly to the disk flux strength. These are the dippers, which are suspected to be caused by dust extinction events. This phenomenon is mostly dependent on disk inclination and can occur even for transitional disks, as is the case for prototype Mon-000660 (also known as V354 Mon). Therefore, as with the lack of variability and temperature trends, the fact that there is no global relationship with disk flux is not particularly surprising.

NGC 2264 has been the subject of many young cluster studies, from photometric and Hα censuses (e.g., Rebull et al. 2002; Lamm et al. 2004) to X-ray (e.g., Ramírez et al. 2004a; Flaccomio et al. 2006; Feigelson et al. 2013) and radial velocity surveys (e.g., Fűrész et al. 2006). To cull a reliable membership list from candidates reported in the literature, we required that objects meet at least two of six criteria: (1) photometric data consistent with the V−I or R−I cluster locus defined by Flaccomio et al. (2006) (see their Section 3.2); (2) strong photometric Hα emission, according to the criteria of Sung et al. (2008), or spectroscopic Hα EW larger than 10 Å; (3) X-ray detection at a flux greater than 10⁻⁴ L_bol (Ramírez et al. 2004a, 2004b; Flaccomio et al. 2006); (4) radial velocity consistent with NGC 2264 membership, as classified by Fűrész et al. (2006); (5) mid-infrared excess indicating a disk (i.e., Class I, II, or flat SED, according to selection methods described in Section 4.2); and (6) spatial location coinciding with an A_V > 7 region of NGC 2264, if an object displays an infrared excess or X-ray emission but is not detected in the optical.

We identified highly embedded objects for criterion (6) by plotting X-ray and mid-infrared source locations on the extinction map produced by Teixeira et al. (2012). These obscured stars were added to our overall membership list for further study, but not included in this paper since we are focusing specifically on targets with CoRoT detections.

The above requirements eliminate most field dwarf and extragalactic contaminants from our membership sample. We have avoided selecting members based on variability detection so as not to bias our statistical analysis of the light-curve types.

While not used extensively in this paper, the high-cadence staring data from Spitzer/IRAC offer an important window into short-timescale infrared variability in YSOs. In preparing it for analysis, we performed a series of procedures to remove systematics. The steps described below may be of general interest to other Spitzer/IRAC users.

A.2.1. Staring Data Quality and Correction of Systematics

The pixel-phase effect, or variation of measured flux as a function of sub-pixel position, is a known issue affecting the IRAC detector (see the IRAC instrument handbook, v2.0.3). Intrapixel sensitivity variations introduce two systematic effects to the staring (i.e., continuous, non-dithered) light curves: first, the flux may vary up to 6% as the pointing shifts gradually from one side of a pixel to the other; second, the flux exhibits oscillations at the percent level on <1 hr timescales. The first effect reflects the pointing stability of Spitzer, while the second was determined by the engineering team to be associated with the cycling of a battery heater on board the satellite, causing periodic flexure between the star tracker and the cold focal plane. The result is that the pointing center on the detector oscillates within individual pixels, convolving their sensitivity variations with the intrinsic flux. In 2010 October, the Spitzer engineering team was able to reduce both the amplitude and period of the heater-associated temperature fluctuations, the latter from ∼60 minutes to ∼40 minutes. However, the effect remains prominent at the 0.5% level in both channels, comparable to the white-noise contribution for bright stars in the sample.

Our goal of assessing infrared variability of young stars on short timescales and at low amplitudes necessitates the removal of pixel-phase variations to the extent possible from the data. A pixel-phase correction exists for Warm Spitzer observations, but tests on our light curves revealed that it works for stars on some pixels, but not others. Several techniques have been developed to provide a more consistent correction, many of which involve fitting a polynomial to the flux as a function of detector X and Y position (e.g., Bonfils et al. 2011). Unfortunately, this approach does not work well for many of our targets, as the erratic variability that we seek to study cannot be modeled analytically and thus prevents a robust fit. In Cody & Hillenbrand (2011), we adopted a different technique, involving a Gaussian model of the intra-pixel sensitivity variation, for which free parameters were fit by minimizing the height of the pixel-phase oscillation peak in the periodogram. Maps of the pixel sensitivity were determined independently for each star, and the inferred flux variation as a function of X and Y position was then subtracted from the light curves.

While this method provided satisfactory corrections without compromising intrinsic stellar variability, it was computationally intensive. To apply a similar technique to our ∼500 staring targets, we require an algorithm with fewer free parameters. Fortunately, during our 2011 observations of NGC 2264, object centroids were stable to ∼0.08 pixel (i.e., 01) in the X detector position over the course of each 16–26 hr staring mode observation. Flux variations are much better correlated with the Y centroids, which vary by 0.3–0.4 pixels (03–05). Plots of the normalized flux versus Y for stars lacking obvious variability by eye revealed nearly linear trends, as shown in Figure 46. We therefore chose a basic pixel-phase correction that fits a single slope characterizing the pixel sensitivity as a function of Y position. To determine this slope, A, we express the sensitivity, s, as a linear function of Y position on the detector:

$$\begin{equation} s(y)=A(y-y_{\rm med})+1. \end{equation} \tag{ A1 }$$

Here we have normalized the sensitivity to a median value of 1.0 in the area of the pixel where the data fall. Applying this sensitivity function, the true flux data, d_true(y), will be observed as

$$\begin{equation} d_{\rm obs}(y)=d_{\rm true}(A(y-y_{\rm med})+1). \end{equation} \tag{ A2 }$$

Zoom In Zoom Out Reset image size

Note that we have specified the observed and true data to have the same median. To correct the observed data, we then invert this equation, exploiting the fact that y − y_med ≪ 1:

$$\begin{equation} d_{\rm true}(y)=d_{\rm obs}/(A(y-y_{\rm med})+1)\sim d_{\rm obs}(1-A(y-y_{\rm med})). \end{equation} \tag{ A3 }$$

We wish to determine the value of A so that the corrected light curve is devoid of pixel-related systematics. To do this, we take advantage of the pixel-phase effect, which causes y(t), the detector position, to oscillate rapidly as a function of time. The period of oscillation is approximately 40 minutes. To determine this value more accurately, we compute a Fourier transform periodogram for each chunk of staring data, y(t), first subtracting out its median value y_med to remove low-frequency systematics:

$$\begin{eqnarray} FT(f,y(t_k))&=&\frac{2}{N}\left[\left(\sum _{k=0}^{N}\sin (-2\pi ft_k)(y(t_k)-y_{\rm med})\right)^2\right.\nonumber\\ &&\left.+\,\left(\sum _{k=0}^{N}\cos (-2\pi ft_k)(y(t_k)-y_{\rm med})\right)^2\right]^{1/2},\nonumber\\ \end{eqnarray} \tag{ A4 }$$

where f is frequency in the periodogram, k is the index of points in the time series, and N is the total number of points. We numerically locate the highest peak in this periodogram and note its frequency, f₀. We find a value of ∼37.5 cycles per day; this can be seen with the solid curve in the bottom panels of Figure 47. The corresponding pixel-phase period is 1/f₀ ∼ 38.4 minutes.

Zoom In Zoom Out Reset image size

We will now use the oscillatory behavior of y(t) to determine the best value of the slope A in the pixel sensitivity function s(y). We do this by minimizing the value of the Fourier transform periodogram of the data themselves, d(y(t_k)), at f₀. The value of A resulting from this process can then be used to remove the effects of pixel sensitivity variation from the data without compromising variability. We begin by writing out the Fourier transform periodogram for the true data, again subtracting out the median value, d_med, to eliminate periodogram systematics:

$$\begin{eqnarray} && FT(f,d(y(t_k))=\displaystyle\frac{2}{N}\left[\left(\sum _{k=0}^{N}\sin (-2\pi ft_k)(d_{\rm true}(t_k) -d_{\rm med})^2\right.\right.\nonumber\\ &&\left. +\left(\sum _{k=0}^{N}\cos (-2\pi ft_k)(d_{\rm true}(t_k) -d_{\rm med})\right)^2\right]^{1/2}.\nonumber\\ &&=\! \frac{2}{N}\left[\!\!\left(\!\sum _{k=0}^{N}\sin (-2\pi ft_k)(d_{\rm obs}(t_k)\,{\cdot}\, (1\,{-}\,A(y_k\,{-}\,y_{\rm med})) \,{-}\,d_{\rm med}))\!\!\right)^2\right. \nonumber\\ &&\left.+\!\left(\!\sum _{k=0}^{N}\cos (-2\pi ft_k)(d_{\rm obs}(t_k) \,{\cdot}\, (1\,{-}\,A(y_k\,{-}\,y_{\rm med})) \,{-}\,d_{\rm med})\!\!\right)^2\right]^{1/2}\!\! . \nonumber\\ \end{eqnarray} \tag{ A5 }$$

We next minimize the periodogram with respect to the slope, A, at the peak frequency f₀:

$$\begin{equation} \frac{d{\rm FT}(f_0,d(y(t_k))}{dA}=0. \end{equation} \tag{ A6 }$$

It can be shown that this equation takes the form

$$\begin{equation} 0=AC_1+C_2+AC_3+C_4, \end{equation} \tag{ A7 }$$

where C₁, C₂, C₃, and C₄ are constants, as follows:

$$\begin{equation} C_1=-\left[\sum _{k=0}^{N}\sin (-2\pi f_0t_k)\frac{d_{\rm obs}(t_k)}{d_{\rm med}}\cdot (y_k-y_{\rm med}) \right]^2, \end{equation} \tag{ A8 }$$

$$\begin{eqnarray} C_2&=&\left[\sum _{k=0}^N \sin (-2\pi f_0t_k)\left(\frac{d_{\rm obs}(t_k)}{d_{\rm med}}-1\right) \right] \nonumber\\ &&\times \left[\sum _{k=0}^{N}\sin (-2\pi f_0t_k)\frac{d_{\rm obs}(t_k)}{d_{\rm med}}\cdot (y_k-y_{\rm med})\right], \qquad \end{eqnarray} \tag{ A9 }$$

$$\begin{equation} C_3=-\left[\sum _{k=0}^{N}\cos (-2\pi f_0t_k)\frac{d_{\rm obs}(t_k)}{d_{\rm med}}\cdot (y_k-y_{\rm med}) \right]^2, \end{equation} \tag{ A10 }$$

$$\begin{eqnarray} C_4&=&\left[\sum _{k=0}^N \cos (-2\pi f_0t_k)\left(\frac{d_{\rm obs}(t_k)}{d_{\rm med}}-1\right) \right] \nonumber\\ &&\times \left[\sum _{k=0}^{N}\cos (-2\pi f_0t_k)\frac{d_{\rm obs}(t_k)}{d_{\rm med}}\cdot (y_k-y_{\rm med})\right]. \qquad \end{eqnarray} \tag{ A11 }$$

Operating on the four sections of by-BCD staring light curves, we determine analytically the slope A that minimizes the pixel-phase oscillation signal by computing

$$\begin{equation} A=\frac{-(C_2+C_4)}{C_1+C3}. \end{equation} \tag{ A12 }$$

Finally, we divide out the pixel response from the light-curve sections by calculating d_true as in Equation (3). Because the X and Y positions are well correlated in time, elimination of the variation associated with shifts in Y only should effectively remove flux trends in X as well (see Figure 46).

Application of our algorithm resulted in a substantial reduction in the level of systematic variability in our light curves, as indicated by both visual inspection and measurement of the height of the periodogram peak at the 40 minute pixel-phase oscillation period. An example comparing our modified correction with the standard correction and a raw light curve is shown in Figure 47. For approximately 1% of light curves, the pixel sensitivity distribution was not well modeled by a linear fit and our correction introduced oscillatory features that were not present in the original light curve. In these cases, we retained the raw light curve or the version corrected with the standard Warm Spitzer pixel-phase prescription, depending on which one displayed lower levels of systematics.

The rms values of the corrected 3.6 μm staring light curves ranged from 0.003 mag (at a magnitude of 7.5) to 0.1 mag (at a magnitude of 15.5). The 4.5 μm light curves had similar but slightly lower precision (by a factor of ∼1.5). The rms values as a function of magnitude for one of the staring photometry sections are presented in Figure 48; these are consistent with the predicted uncertainties based on Poisson noise and sky background.

Zoom In Zoom Out Reset image size

A.2.2. Merging the Mapping and Staring Data

We assemble in Figure 49 the entire 162 member disk-bearing data set observed by CoRoT and Spitzer. Small black dots are optical data, light gray points are 3.6 μm data, and dark gray points are 4.5 μm data.

Zoom In Zoom Out Reset image size

View extended figure (8 panels)