YSOVAR: MID-INFRARED VARIABILITY IN THE STAR-FORMING REGION LYNDS 1688

3.1.1. Spitzer Observations

Three fields in L1688 were observed with Spitzer in four observing windows from 2010 April 12 to 2010 May 16 (visibility window 1), 2010 September 22 to 2010 October 27 (visibility window 2), 2011 April 20 to 2011 May 23 (visibility window 3), and 2011 October 1 to 2011 November 6 (visibility window 4). These windows are consecutive visibility periods dictated by the Spitzer orbit (Werner et al. 2004). They are shown visually in Figure 2. In the first visibility window, the sampling is much denser in time than in the later visibility windows. For most sources, about 70 observations with irregular time intervals to aid period detection were obtained in visibility window 1. In the first visibility window, we use a repeating pattern of eight observations every 3.5 days. Within the 3.5 day period the time step increases from 2 to 16 hr. In the last three visibility windows, the time steps increase linearly with step lengths of roughly 1, 2, 3, ... days to again sample multiple variability frequencies equally. This irregular sampling minimizes aliasing with the observation frequency. Less than 10 data points per source are taken in each of the later three visibility windows. In total, there are 108 observations with a total mapping time of 30.7 hr. Table 1 lists the time of each observation. They can be found under Program Identification number (PID) 61024 in the Spitzer Heritage Archive. Each observation consists of six dithers in IRAC mapping mode using the high-dynamic-range data-acquisition mode, which obtains a 0.4 s and a 10.4 s exposure for each pointing.

Zoom In Zoom Out Reset image size

Table 1. Observation Log

AOR	Start Time	End Time
29267200	2010 Apr 12 11:00:41	2010 Apr 12 11:18:44
29266688	2010 Apr 12 17:02:37	2010 Apr 12 17:28:01
29266176	2010 Apr 12 22:48:49	2010 Apr 12 23:09:28
29265664	2010 Apr 13 08:44:31	2010 Apr 13 09:05:35
29265408	2010 Apr 13 20:09:17	2010 Apr 13 20:26:33

Only a portion of this table is shown here to demonstrate its form and content. Machine-readable and Virtual Observatory (VO) versions of the full table are available.

Download table as:  Machine-readable (MRT)Virtual Observatory (VOT)Typeset image

The three fields chosen were observed with the IRAC 1 and IRAC 2 channels (effective wavelengths 3.6 μm and 4.5 μm). Both channels operate simultaneously, but their fields of view are nonoverlapping. Thus, each target field is observed in two consecutive pointings, one for each channel. A secondary field is observed in the secondary channel, and the primary channel is observing the target field. In visibility windows 1 and 3, sources south of the main fields have only IRAC 1 data, and those to the north only have IRAC 2 data. In visibility windows 2 and 4, the situation is reversed. Not all sources in the central fields have usable data in both bands because they might be too bright or too faint in one channel or fall on the edge of the map. Additional sources with two-band coverage are found where the northern-side field of one target field overlaps with the southern-side field of another target field.

3.1.2. Spitzer Data Reduction

3.1.3. Instrumental Effects Remaining in IRAC Light Curves

Despite the careful data reduction described above, some residual artifacts remain in the Spitzer light curves. In this section, we search for artifacts that are related to the position on the detector. Compared with other clusters in the YSOVAR project, L1688 is particularly suited to discover these kinds of effects because the fields observed in L1688 have almost no rotation within one visibility window. The spacecraft orientation flips between visibility windows, so most instrumental artifacts produce light curves that have one level in visibility windows 1 and 3 and another level in visibility windows 2 and 4. Figure 3 shows sources with a magnitude between 8 and 15 in 4.5 μm that might fall in this category. The relative difference in magnitude is smaller for brighter sources, but because the photometric errors are also smaller, these instrumental effects can still be significant.

Zoom In Zoom Out Reset image size

We visually inspected every frame for a sample of sources with light curves similar to those in Figure 3 to identify the cause for the artifacts. We found that several different effects can cause these steps in the light curve. Some sources are close to a detector edge such that the background is not well determined, some are in the wings of the point-spread function (PSF) of a bright neighbor, and some show residuals from hot pixels. Proximity to a detector edge or hot pixel affects only one of the spacecraft orientations; even the PSF wings change with the orientation because the PSF is not circular.

Sufficiently strong intrinsic variability can mask this offset between visibility windows, and for any individual source, this instrumental effect cannot be distinguished from intrinsic variability with a period of 1 yr. Thus, we use a statistical approach to quantify the number of light curves that suffer from this problem. For each channel, we have about 200 light curves with data points in all four observing visibility windows (the total number of light curves is larger, but only for sources in the three target fields do we have data for all four visibility windows). If the mean magnitude in a visibility window depends on the detector position of the source, this source will on average be brighter in visibility windows 1 and 3 and fainter in visibility windows 2 and 4 (or vice versa). We calculate the mean magnitude in each visibility window and test whether the two brightest mean magnitudes belong to visibility window [1, 3] or [2, 4] or any other combination ([1, 2], [1, 4], [2, 3], [3, 4]).

We use the following abbreviated notation. When we calculate the mean magnitudes in each visibility window for a sample of n = a + b + c sources, then a:b:c means that a sources have the brightest mean magnitudes in visibility windows [1, 3], b sources in visibility windows [2, 4], and c sources have their two brightest mean magnitudes in any other combination of visibility windows. If the difference in the mean magnitude between visibility windows is unrelated to the position on the detector, we expect the ratio 1/6: 1/6: 2/3 = 17%:17%:67%. In contrast, the observed light curves have the ratio 47:66:75 = 25%:35%:40% for IRAC 1 and 52:59:95 = 25%:29%:46% for IRAC 2. This is incompatible with the expected multinomial distribution (the probability to observe a distribution at least as far from the expected 1/6: 1/6: 2/3 by chance is <10⁻⁸). In each case, about 50 light curves, a quarter of the sample, need to be shifted from the first two bins to the last bin to make the observed distribution compatible with the expected distribution. This implies that about a quarter of all light curves suffer from the artifacts described above. However, both IRAC channels use independent detectors, and thus the chance that both channels are affected for the same source and that the effect goes the same way (bright-faint versus faint-bright) is low.

Below in Section 4 we identify variable sources using a Stetson and a χ² test. The limits in those tests are designed to be conservative, and indeed we find a distribution of 5:7:24 = 14%:19%:67% for IRAC 1 and 5:5:27 = 14%:14%:73% for IRAC 2 when we restrict the sample to those light curves that will be classified as variable below. Both ratios are fully compatible with the expected multinomial distribution 1/6: 1/6: 2/3. This shows that the limits we apply are conservative enough that the sample of stars we identify as variable has no or only few sources where the variability is not due to intrinsic source variability. Therefore it is not necessary to remove any source based on magnitude jumps between visibility windows. Only one of the four light curves shown in Figure 3 (SSTYSV J162727.53-242611.2) will be identified as variable below. We visually inspected all light curves that are marked as variable, and, apart from the example shown above, we did not see light curves where the variability seems to be due to the pattern discussed in this section.

In summary, about a quarter of all light curves are affected by detector position-dependent artifacts. In extreme cases, the associated jumps reach 0.5 mag for a 14 mag source and up to 0.05 mag in bright sources around 8 mag. However, we show statistically that the definition of variability we use is so conservative that this effect does not contribute a significant number of objects to our sample of variable sources.

Disk-bearing YSOs can be identified from Spitzer data alone, but information in other spectral bands is required to find the other cluster members. X-ray observations are one way to identify diskless (class III) YSOs. In the IR the SEDs of those sources are indistinguishable from a main-sequence field star, but because of their rapid rotation, YSOs are much brighter in X-rays than are field stars (Feigelson & Montmerle 1999).

L1688 was observed by Chandra on 2000 April 13 for 100 ks exposure time in the FAINT mode with the ACIS instrument (ObsId 635). We reprocessed this exposure with the ANCHORS pipeline (Wolk & Spitzbart 2007) using a recent calibration; see discussion in Paper I. These data has been analyzed in detail to study the distribution of X-ray properties in CTTS and to identify brown dwarfs in this star-forming region (Imanishi et al. 2001a, 2001b).

During the observation, five ACIS chips were operational, four from the central ACIS-I imaging array as well as one ACIS-S chip. The PSF degrades significantly for sources located off-axis, and thus the coordinates of the outer sources are less reliable. To cross-match X-ray sources with Spitzer sources, we used a matching radius of 1'' for X-ray sources within 3' of the optical axis of Chandra, 15 for sources between 3' and 6' away from the optical axis, and 2'' for all sources located farther than 6' from the optical axis. The observed Chandra field overlaps about two-thirds of the area covered in the Spitzer monitoring. However, the variable PSF leads to a sensitivity that varies over the observed field; thus the absence of an X-ray detection for Spitzer sources does not necessarily imply the absence of X-ray emission. The ACIS detector has an intrinsic energy resolution, and we use the net flux to characterize the X-ray properties. For sources with more than 20 counts, we also fit an absorbed single-temperature APEC model (Foster et al. 2012) with abundances fixed at 0.3 times the solar value from Anders & Grevesse (1989).

Sources are extracted down to a very low significance. In total, there are 315 detected X-ray sources, but only 31 of them match an object with a Spitzer light curve. We disregard all unmatched sources; for sources with a light curve and an X-ray counterpart, the X-ray properties are given in Table 2. To estimate the number of spurious matches, we multiply the fraction of the total survey area that is included in the positional error circles of the X-ray sources with the number of Spitzer sources with light curves. The result is the average number of spurious matches. We expect at most two to four Spitzer sources to be matched to a spurious X-ray source. X-ray sources that are cross-matched successfully are marked in Figure 1.

The star-forming region L1688 has been the target of intense study over the past decades, and a wealth of additional information exists in the literature. In particular, we refer the reader to two reviews (Wilking 1992; Wilking et al. 2008). The latter review compiles a list of objects with a high probability of membership from a variety of published sources. The membership criteria employed are (1) X-ray emission, which–at the distance of L1688 is detectable only from young and thus active stars; (2) optical spectroscopy, with Hα in emission or Li in absorption; (3) a location above the main sequence in the HR diagram; or (4) IR emission that is indicative of a circumstellar disk.

L1688 was also observed with Spitzer in the cryogenic mission phase with all four IRAC channels and the 24 μm channel of the Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004). Objects classified as YSOs from these data (Padgett et al. 2008) are already contained in the membership list of Wilking et al. (2008). We augment our own Spitzer data reduction with values from the catalog published by the c2d project (c2d = "From Cores to Disks"; Evans et al. 2003). If we did not obtain a photometric value for an IRAC band or the 24 μm MIPS but a value with the quality specifier A, B, or C is present in the c2d catalog, then we use that value. The data are given in Table 2, which specifies if a data point is taken from our own data reduction ("G09": using the pipeline from Gutermuth et al. 2009)¹⁷ or the c2d database.

Near-IR data are taken from 2MASS (Skrutskie et al. 2006) and cross-matched by the cluster grinder pipeline (Gutermuth et al. 2009). Additionally, we take detections from the UKIRT Infrared Deep Sky Survey (UKIDSS) Galactic cluster survey, data release 9. The UKIDSS project is defined in Lawrence et al. (2007). UKIDSS uses the United Kingdom Infrared Telescope (UKIRT) Wide Field Camera (Casali et al. 2007) and a photometric system described in Hewett et al. (2006). The pipeline processing and science archive are described in Hambly et al. (2008). We only retain detections with a mergedClass flag between −3 and 0. Sources brighter than m_K = 10 mag can be saturated, and we discard the UKIDSS data for any source that lies above this threshold either in UKIDSS or 2MASS. This leaves us with 2MASS data in the range m_K = 5–16 mag and UKIDSS m_K = 10–19 mag. The luminosity function for both data sets is almost identical from m_K = 10 to m_K = 15. For fainter sources, 2MASS is incomplete; UKIDSS is incomplete for m_K > 17 mag.

The YSOVAR data is also cross-matched with data from the SIMBAD service to provide an identification with known objects from the literature.

In all cases the matching radius is set to 1''. If a catalog contains multiple entries within 1'' of a YSOVAR source, we match it to the closest catalog entry. In some cases the best cross-match is not obvious. Those sources are discussed in Appendix A.

The UKIDSS has a better spatial resolution than our IRAC data. There are nine sources where more than one UKIDSS source is found within the size of the aperture that we use for IRAC photometry. In six cases (WSB 52, ISO-Oph 152, ISO-Oph 131, SSTYSV J162728.13-243719.6, SSTYSV J162718.11-244814.1, SSTYSV J162718.25-244955.8), the second source is visible in the K band, the UKIDSS band that is closest in wavelength to the IRAC data, so it is likely that both sources contribute to the observed IRAC emission. In the remaining three cases (ISO-Oph 28, ISO-Oph 57, SSTYSV J162741.14-242038.3) the second source is not visible at the K band.

Table 2 contains the position, the designation, the flux densities of each source, and the properties of their light curves. The properties of the light curve will be discussed in detail in the remainder of this article. Most properties of the light curve, e.g., mean, minimum, and maximum, appear twice. They are calculated once over the entire available light curve and once for the first visibility window only; the fast-cadence sampling is available uniformly for all clusters in the YSOVAR project, and thus values calculated over the fast cadence only can be compared between clusters (see, e.g., Paper I). A subset of the properties of light curves with a mean magnitude <15 in IRAC 1 or IRAC 2 is shown in Table 3.

We build two L1688 membership lists based on different criteria. The first is defined in Paper I and is applied uniformly for all clusters in the YSOVAR project. Sources are treated as cluster members if they fulfill at least one of the following criteria: (1) they are classified as YSOs by Gutermuth et al. (2009) based on their IR excess in cryogenic mission Spitzer data or (2) they are detected as X-ray sources in Chandra imaging and have a spectral slope compatible with a stellar photosphere (SED class III; see Section 3.6). At the distance of L1688 (we use 120 pc from Lombardi et al. (2008), but see also the discussion in Paper I), cluster members that are young and thus still magnetically active stars can be detected in X-rays. A total of 57 sources fulfill one or both conditions. In Table 2, these objects are marked as "member (YSOVAR)". The main biases in this sample are that the IR criterion selects only those members with disks and not class III sources, and the Chandra criterion suffers from incomplete spatial coverage and may include late-type foreground stars. The different biases are a common problem in multiwavelength studies of star-forming regions (see, e.g., Preibisch et al. 2011; Wolk et al. 2011; Feigelson et al. 2013; Naylor et al. 2013).

The second membership list is taken from Wilking et al. (2008). On this list, 74 of our 884 sources with Spitzer light curves are cluster members; 51 of those 74 are also included in the YSOVAR membership list because of their X-ray or IR emission. When we compare properties of members and nonmembers below without referring to a specific set, then we mean membership as defined by the standard YSOVAR criteria.

Six sources in the YSOVAR standard member set are not part of the Wilking et al. (2008) list. More details on those sources are given in Appendix B. On the other hand, all of the sources in Wilking et al. (2008) that are not part of the YSOVAR standard set were selected based on an X-ray detection from observations other than Chandra (for the specific references, see Wilking et al. 2008).

Figure 4 shows a histogram of the observed mean magnitudes in 3.6 μm. Almost all bright sources are cluster members. The fraction of members drops below 11th mag, and no cluster member is found below 14th mag. We expect background sources to be fainter than cluster members because of their larger distance and because they are seen through the cloud, but also fainter sources cannot be reliably classified. The thickness of the cloud is inhomogeneous. The three primary fields cover roughly A_V = 20–35 mag in the extinction maps from the c2d project (Evans et al. 2003). This corresponds to A_[3.6] = 1.3–2.2 mag (Cardelli et al. 1989; Indebetouw et al. 2005). At the distance of L1688, a K5 star with little or no extinction has m_[3.6] = 12.2 mag according to the evolutionary tracks of Siess et al. (2000) and the YSO colors from Pecaut & Mamajek (2013). Thus, even late-K or early-M star members in the cloud are bright enough to be detected.

Zoom In Zoom Out Reset image size

Just over 10% of the 884 sources with light curves are classified by Gutermuth et al. (2009) at all because this classification scheme is conservative, representing a set of sources that can be dereddened and classified to a high degree of confidence. However, for many sources, we lack the required spectral coverage. Of those sources that can be classified in this way, 16 are class I sources, 36 are class II sources, and 59 are class III candidates, i.e., sources with a weak or absent IR excess in their SED (in this case the IR SED does not provide sufficient information to decide if a star is a YSO or a field star). Ten of those 59 have an X-ray counterpart in our Chandra data.

In order to classify all sources, including those not classified in Gutermuth et al. (2009), Paper I defines a simpler approach, which uses the observed colors only. In this scheme, we fit the spectral slope α = dlog λF_λ/dlog λ, where λ is the wavelength and F_λ the flux density per unit wavelength interval at that wavelength. We make use of all measured flux densities (no upper limits) in the 2–24 μm range, which corresponds to the range from the K_s or K filter, which we take from 2MASS or UKIDSS, respectively, to the MIPS 24 μm channel. We use the flux densities for all Spitzer bands from the observations presented in Gutermuth et al. (2009; and we include all sources extracted using this pipeline, including those that are not published in Gutermuth et al. 2009 because they are not classified as YSO in that work) or c2d and summarize the new light curves by calculating the mean and standard deviation for the new IRAC1 and IRAC2 data. Thus, sources detected both in the cryogenic mission and in the new data set will have two data points for IRAC1 and IRAC2. If a source is detected in both 2MASS and UKIDSS, we use both values for the fit because they are independent measurements. We perform a least-squares fit and call sources with 0.3 < α class I, −0.3 < α < 0.3 flat-spectrum (F), −1.6 < α < −0.3 class II sources and α < −1.6 class III candidates. The classifications are given in Table 2 in two columns. The s1SED presents the derived SED class using all available literature data as described above and the mean of the light curves from visibility window 1 for comparison with Paper I; SED uses the same literature data but the mean flux density for the light curves calculated for all visibility windows. In this scheme, a significantly reddened main-sequence star, which has an intrinsic slope α < −1.6, may appear as class II object, so this observational classification cannot be translated directly into the evolutionary stage of an object. Of the remaining sources, we find 110 class I, 78 flat-spectrum, and 455 class II sources and 221 class III candidates. We cannot classify 20 sources because they are seen in one band only.

Not only the evolutionary status but also other properties of the individual source such as inclination or stellar mass influence the value of α. Massive stars emit more energetic radiation and can thus change the structure of the accretion disk. However, L1688 is a region without massive YSOs (Wilking et al. 2008), and we define α using wavelengths longward of the K band far from the peak of the stellar SED, so the shape of the photospheric SED has only negligible influence on the total SED. Robitaille et al. (2006) simulated different YSOs of low mass. Their results show that SEDs depend mostly on the evolutionary stage except for very extreme cases such as stars with an edge-on disk. Thus, α provides a good proxy for the evolutionary state of the YSOs in L1688.

Comparing with Gutermuth et al. (2009), the resulting classification is very similar, particularly for the class I and flat-spectrum sources. Most sources with a spectral slope <0.3 according to our slope-fitting are class I sources in both classification schemes. With two exceptions, our flat-spectrum sources are either class I or II in Gutermuth et al. (2009). Of the 52 sources with a spectral slope between −1.6 and −0.3, 21 are also called class II by Gutermuth et al. (2009), but 30 are class III in that paper, indicating that a significant fraction of what we call class II might indeed be reddened background stars. The general agreement between the more complex classification scheme and the observed spectral slope is also found for other star-forming regions (Paper I).

Sources in the target fields are observed in IRAC 1 and IRAC 2 almost simultaneously (within a few minutes). For those sources, we calculate the Stetson index s with points weighted evenly (Stetson 1996):

$$\begin{equation} s = \frac{1}{\sqrt{N(N-1)}}\sum _{k=1}^{N} sig(P_k)\sqrt{|P_k|}, \end{equation} \tag{ 1 }$$

where the sum is taken over all N pairs of observations in IRAC1 and IRAC2 with observed magnitudes a_k and b_k and uncertainties $\sigma _{a_k}$ and $\sigma _{b_k}$; sig denotes the sign of P_k. The variable P_k is the product of the normalized residuals in both bands:

$$\begin{equation} P_k = \left(\frac{a_k-\bar{a}}{\sigma _{a_k}} \right) \times \left(\frac{b_k-\bar{b}}{\sigma _{b_k}} \right)\protect. \protect \end{equation} \tag{ 2 }$$

Here $\bar{a}$ denotes the error-weighted mean of all IRAC 1 mag and $\bar{b}$ the error-weighted mean of all IRAC 2 mag. The Stetson index is very robust to observational errors because those are unlikely to affect both bands in the same way. Following Paper I, we define a source as variable if s > 0.9. In Paper I, this is calculated over the first visibility window only; in this article we use the light curve from all four visibility windows. The Stetson index calculated for the first visibility window only can be found in Table 2 in column stets1 and for the entire light curve in column stet. Figure 5 shows examples of light curves that are classified as variable according to the Stetson index. In the first visibility window, we find 34 sources to be variable, that is, 18% of all light curves where we can calculate the Stetson index. The number is similar (38 sources, 18%) if we consider all observing visibility windows. This is not surprising, given that the Stetson index is designed to be a "robust" statistic that is little influenced by a small number of data points, and compared to the first visibility window, the other visibility windows only contribute a small number of data points. We caution that Equation (1) will introduce a systematic bias in such a way that a source that is variable in only one band will not be recognized as variable.

Zoom In Zoom Out Reset image size

For sources outside the primary target fields or for sources that are too bright or too faint in one IRAC channel, the Stetson index cannot be calculated. For those sources we rely on the reduced chi-squared value to detect variability:

$$\begin{equation} \chi ^2_{{\rm red}} = \frac{1}{N-1}\sum _{k=1}^{N}\frac{(a_k-\bar{a})^2}{\sigma _k^2}\protect. \protect \end{equation} \tag{ 3 }$$

Instrumental uncertainties lead to a non-Gaussian error distribution, so we use a conservative cutoff and mark sources as variable only if $\chi ^2_{{\rm red}} > 5$ (Paper I). Using this metric, we find 22 sources that exhibit variability in IRAC1 and 18 that exhibit variability in IRAC2. For comparison, we note that of the 38 sources classified as variable according to their Stetson index in the last section, 34 are also variable according to the χ² test for their IRAC1 light curve and 37 according to their IRAC2 light curve.

The $\chi ^2_{{\rm red}}$ calculated for the first visibility window only can be found in Table 2 in column chi2s1-36 and chi2s1-45 for IRAC I and IRAC 2, respectively, and for the entire light curve in columns chi2-36 and chi2-45, again for IRAC 1 and IRAC 2, respectively. Examples are shown in Figure 6. If we fit a linear slope instead of comparing to the mean (i.e., a constant), the results are the same. Many sources have significantly larger values of $\chi ^2_{{\rm red}}$ for the whole observation series than within one visibility window, indicating a constant luminosity over 40 days but a change over timescales between 6 months and 2 yr.

Zoom In Zoom Out Reset image size

Last, we search for periodicity in the fast-cadence (visibility window 1) light curves using a Lomb–Scargle periodogram (Lomb 1976; Scargle 1982); again, see Paper I for details. We require a false alarm probability <0.03, and to further reduce the number of false positives, we additionally run a period detection of all light curves on the NASA Exoplanet Archive Periodogram Service.¹⁸ This service employs several algorithms because each algorithm is particularly suited to a different signal shape (see Paper I or the website for details on the other algorithms). We also calculate the autocorrelation function of each signal. We find that the Lomb–Scargle periodogram is advantageous for light curves with <100 points, as is the case for most of our data. We require that at least one of three supplementary algorithms that we run retrieves a periodicity with a similar timescale (see Paper I for details). We only search for periods between 0.1 and 14.5 days so that at least three periods fit in each visibility window. Allowing longer periods leads to the detection of many long-term trends, where the data do not show that these trends are actually periodic. Finally, all algorithmically detected periodicity is vetted by eye. If periods are found in multiple bands, we report the period in IRAC 1, which is generally the most reliable due to the lower measurement uncertainties. If IRAC 1 does not reveal a periodicity, we report the value for the IRAC 2 light curve, and as a last resort the period found in the [3.6]–[4.5] color includes objects where a periodic signal is overlaid by a long-term trend so that it is undetectable in each individual channel but might be visible in the color. The final list of adopted periods is given in Table 4. It contains sources from all SED classes.

We find 14 sources that show a periodic behavior in the first visibility window, but this periodicity is not stable over more than 1 visibility window. The data points taken in the first visibility window follow the folded period much better than do the data from the later visibility windows, indicating that the period is not stable for longer than a few months. In ROXN 41, the data points of the later visibility windows, folded with the same period, seem to follow a different, yet clearly defined, light curve with a similar period.

Ten periodic sources have information in both IRAC bands. Eight of those are already classified as variable by the Stetson index. Three out of four sources with data in one band only are variable according to the $\chi ^2_{{\rm red}}$ test. The remaining two sources with information in both IRAC bands show periodicity but fail the Stetson index test because a larger variability amplitude is required for a significant detection in Stetson index, which does not make any assumption about the form of the variability compared with the Lomb–Scargle periodogram, which only detects periodic signals. Equally, the remaining source with data in only one band fails the $\chi ^2_{{\rm red}}$ test because the $\chi ^2_{{\rm red}}$ test also requires a larger variability amplitude for a significant detection than the Lomb–Scargle periodogram.

The number of data points in visibility windows 2, 3, and 4 is too low to search for periodicity in those visibility windows alone. Figure 7 shows the phase-folded light curve in one band for six periodic sources where the periodicity is significant in a single band. Figure 8 shows the phased light curve for CFHTWIR-Oph 29, where the period is seen only in the color term [3.6] [4.5]. Most of the periods found in Table 4 are in the range of three to seven days. Only two sources have longer periods of 11 and 14 days. The largest amplitude is around 0.3 mag, and the smallest around 0.05 mag.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

All periodic sources except ROXN 41, SSTYSV J162636.08-242404.2, 2MASS J16271881-2448523, and SSTYSV J162622.19-242352.2 are cluster members according to our membership criteria. ROXN 41 and SSTYSV J162636.08-242404.2 have class I and II SEDs, respectively. All but CFHTWIR-Oph 29 (see Appendix B), SSTYSV J162636.08-242404.2 (see Appendix B), and 2MASS J16271881-2448523 (no information beyond the JHK magnitudes is available in the literature) are cluster members according to Wilking et al. (2008).

Of the 56 cluster members, 10 are periodic; in total there are 14 periodic sources out of 60 sources that are variable in the first visibility window. In comparison, Carpenter et al. (2001) find 18% of all variable stars to be periodic in JHK_s monitoring of the Orion A molecular cloud with a time coverage similar to what we have for L1688. Parks et al. (2014) find one-third of all variable stars in L1688 to be periodic; because of their longer time baseline, Parks et al. (2014) are sensitive to different periods, but only four out of their 32 periods have values that are outside of the range we could detect; see their table. We are sensitive to periods of up to 14.5 days here.

The largest group of light curves in the sample of L1688 does not have any significant variability. The lower limit where variability is detected depends on the brightness of the object—fainter sources need a stronger relative variability because of the larger measurement uncertainties. Paper I presents Monte Carlo simulations to show that the Stetson test finds variability when the amplitude is a few times larger than the average uncertainty; the exact number depends on the signal shape. For example, the variability in a source that switches between two discrete levels will be detected when the amplitude is at least twice the uncertainty (2 × 0.015 mag for a star with magnitude 14). If data from only one band are available, the step size must be more than four times the uncertainty to be found by the χ²-test (4 × 0.015 mag for a star with magnitude 14). Given the observing cadence, the Monte Carlo simulations show that we are sensitive to periods between 0.1 and 14.5 days (Paper I).

In summary, we call a source variable if it is either periodic with a low false-alarm probability (14 sources) or fulfills one of the following conditions: sources with simultaneous data in 2 bands need to have a Stetson index >0.9 (34 sources in the YSOVAR standard set fast-cadence data and 38 sources in total), and sources without simultaneous data need to have $\chi ^2_{{\rm red}} > 5$ (22 sources in the YSOVAR standard set fast-cadence data and 29 sources in total).

Figure 9 contains the variability fraction sorted by SED class (see Section 3.6). Within each class, the population is subdivided into previously identified members: bright (<14 mag in 3.6 μm or 4.5 μm) stars not previously identified as members and faint (>14 mag) stars not previously identified as members. Sources are considered variable if they have a low false-alarm probability for periodicity or a Stetson index >0.9 (if simultaneous data in IRAC 1 and IRAC 2 exist) or $\chi ^2_{{\rm red}}>5$ (if no simultaneous data exist). Note that the classes shown here are based on the observed SED, and thus background nonmember stars seen through the cloud might appear red (like a class II source) even if their intrinsic SED slope is compatible with SED class III. No star classified as a member is fainter than 14 mag.

Zoom In Zoom Out Reset image size

The observed fraction of variable stars as a function of SED class is the best guess for the probability of finding variability for a star of a given SED class. We calculate the uncertainty for the probability based on the observed number of variable stars v in a bin with N total sources as follows. If the true probability for a source in the bin to be variable is p, then the probability of observing v variables in N stars is given by p_obs, which follows a binomial distribution:

$$\begin{equation} p_{{\rm obs}}(v,N,p) = {N \choose v} p^v (1-p)^{N-v}\protect. \protect \end{equation} \tag{ 4 }$$

The peak of this distribution is p_obs(p) = v/N. We then calculate the boundaries of the confidence interval p₁..p₂ for p such that p_obs(p₁) = p_obs(p₂) and the area under the curve includes 68% of the probability. Those uncertainty ranges are shown in Figure 9.

Almost all known cluster members are variable; the fraction of variables is decreasing slightly from class I to II. There are five class III members, and none of them shows variability. This is compatible with the variability fraction observed for class III nonmembers. By definition, sources with a class III SED do not have a large IR excess over a stellar photosphere; thus they are not expected to harbor a substantial disk or show disk-related variability. However, even field stars can show significant variability, at least in the near-IR. Wolk et al. (2013b) find 1.6% of all field stars to be variable in JHK. This includes eclipsing systems, stars with unusually large photospheric spots, and other unidentified variables. All effects that influence the JHK light curve are likely also visible in the mid-IR. Our fraction of variable objects with a class III SED is similar (1.8%) to what Wolk et al. found. In the optical and near-IR, the fraction of variable field stars is lower than that value (Pietrukowicz et al. 2009, 2012; Wolk et al. 2013b).

Spitzer observations of other young clusters also find that younger stars are more variable but with a lower variability fraction. Morales-Calderón et al. (2011) find only 70% of all stars with disks in their ONC sample to have detectable variability. Flaherty et al. (2013) find 60% in IC 348, and Cody et al. (2014) find >90% of all members in NGC 2264 to be variable in the IRAC bands. All three clusters are located at larger distances, and thus observations are less sensitive to small variations; the observations of Cody et al. (2014) are more densely sampled and thus provide a better signal. In L1688, a K5 star with little extinction will have a magnitude of 12.2 in 3.6 μm using the evolutionary tracks of Siess et al. (2000) and the YSO colors from Pecaut & Mamajek (2013). In this magnitude range, variability down to 0.05 mag can be detected (Paper I). In contrast, for the same star in IC 348, the variability has to be about twice as large to be detected and three times as large to be detected in the ONC. Of course, even our observations miss the faint end of the distribution. Thus, it seems likely that essentially all YSOs show substantial mid-IR variability on timescales of days to weeks.

Sources with a class III SED are not contaminated by disk-bearing stars because reddening by the cloud can only increase the spectral slope α but never hide an existing IR excess. The larger fraction of nonmember variable objects in the other classes in Figure 9 thus shows that the membership lists are incomplete. One caveat here is that the class III sources in the sample of members are selected in a different way than the other classes. If true class III cluster members with X-ray emission (that are included in our member sample) are systematically less variable than class III members where we do not detect X-ray emission (that are therefore not included in the member sample), that would also lead to a lower observed variability fraction in class III member sources. However, given the size of the observed effect in variability, such a bias seems unlikely to be the sole reason for the observed distribution.

The observed probability that a source is variable for all classes of bright nonmembers (blue bars) is consistent with around 5%. The fact that this barely changes with the SED class, quite unlike the distribution for the cluster members, indicates that unidentified cluster members cannot make up a large fraction of the nonmember sample because they would bias the observed variability fraction to higher values for the earlier classes. The vastly different number of sources in the different bins makes it highly unlikely that each bin is contaminated by the same fraction of class I–II cluster members. The sample of faint stars has almost the same fraction of variable sources in every SED class. The low fraction can be explained by two effects. First, variability of faint sources cannot be reliably detected unless the amplitude is exceptionally large. Second, because most cluster members are brighter than 14 mag, the faint sample contains fewer unrecognized cluster members than the sample of brighter sources.

Figure 10 shows that there is a wide distribution of amplitudes in YSOs of class I, F, and II. The error bars in this figure represent the standard deviation of the mean amplitudes within one class. Amplitude and standard deviation in the figure are measured in magnitudes; a median amplitude of 0.1 mag means that the flux in 4.5 μm varies by 10% for a typical source. For members, the variability amplitude is larger in class I sources than in flat-spectrum and class II sources, where the amplitude is still larger than in class III sources. Also, the spread of the amplitudes in class I and F is much larger than for class II and III. For nonmembers, the observed amplitudes seem to follow a similar but less pronounced trend. This suggests that the nonmember category still includes not only reddened background objects with an intrinsic class III SED, but also some unidentified cluster members. Indeed, we propose that all variable class I, F, and II sources are members, and Table 5 lists those new members that are contained neither in our YSOVAR standard membership set (see Table 2) nor in Wilking et al. (2008). However, for consistency, we do not modify our sample of members at this point and continue to treat the sources in the table as nonmembers for the remainder of the analysis. Two of the objects in the table, CFHTWIR-Oph 21 and CFHTWIR-Oph 74, were suggested as ρ Oph substellar candidate members by Alves de Oliveira et al. (2010), and the remaining objects have not been studied before.

Zoom In Zoom Out Reset image size

In summary, the probability for any given source to be observed as a mid-IR variable decreases little between class I and II sources, but the mean amplitude and the differences within a class are much larger for sources that have more circumstellar material.

Next, we will look at the timescales of the light curves. The list of periodic sources contains objects of all evolutionary stages, but it is too small to recognize any trends in the period with evolutionary stage. We calculate the discrete autocorrelation function ACF(τ) for each fast-cadence (visibility window 1) light curve to characterize the timescale in all sources, not just the periodic ones. One complication here is that our data are unevenly sampled. In order to calculate the ACF, we linearly interpolate the light curve on a grid with time steps of 0.1 days. This process can change the properties of the light curves on short timescales. However, for most light curves the relevant timescales are longer than the distance between two observations (0.1–0.8 days). We use the following definition of the ACF:

$$\begin{equation} {\rm ACF}(\tau) = \frac{1}{N-\tau }\sum _{k=1}^{N-\tau }\frac{(a_k-\bar{a})(a_{k+\tau }-\bar{a})}{\sigma _a^2}\protect, \protect \end{equation} \tag{ 5 }$$

where N is the total number of points in the light curve a₁, a₂, ..., a_N with a mean of $\bar{a}$ and a standard deviation of σ_a. In the discrete ACF, τ is the number of time steps. By definition, ACF(τ = 0) = 1 and the ACF decay for longer time lags, τ. We take the first value of τ with ACF(τ) < 0.5 as the characteristic timescale for a light curve. A more common definition would be to use the position of the first local maximum in the ACF, but because of the low number of data points in our light curves, the noise in the ACF is large, and this value is often not well defined. Although our definition might not give us the timescale of the physical processes at work, the ACF timescale still provides relative comparisons. We find that the average ACF timescale decreases from four days for class I sources to one day for class III sources (Figure 11). A more detailed analysis of stochastic and quasi-periodic properties of the light curves is beyond the scope of the current paper (for a discussion of these properties in JHK_s light curves, see Parks et al. 2014).

Zoom In Zoom Out Reset image size

The majority of the stars included in our data show no significant variability, and we exclude them from further discussion. The remaining stars are definitely variable, with light curve shapes that display a variety of morphologies. Visual inspection and some quantitative analysis allow us to group these stars into sets with similar light curve morphologies, which can be the first step in attaching physical mechanisms as the cause of the observed variability. Based on a near-IR (JHK) monitoring campaign of cadence and duration similar to ours, Wolk et al. (2013a) identified four light curve morphological classes: (1) periodic; (2) quasi-periodic, which they defined as stars with cyclic brightening and fading, but where the frequency and amplitude of the variations varied from cycle to cycle; (3) long duration, or stars with relatively long-term monotonic changes in brightness over weeks or months, with eventual changes in sign of the variability; and (4) stochastic, which indicated all other variability types, where no obvious pattern to the variability was present.

In another recent paper, Cody et al. (2014) analyzed a 30 + day monitoring campaign for the star-forming region NGC 2264, using optical data from CoRoT and IR data from IRAC, to assign light curve morphology classes to their young star set. Their proposed light curve morphologies were: (1) periodic: stars whose light curves show periodic waveforms whose amplitude and shape are unchanging or only change in minor ways over the 30 day observing window. These light curves were ascribed generally to cold spots on the stellar photosphere; (2) dippers: stars showing a well-defined maximum brightness, upon which are superposed flux dips of variable shape and amplitude. These are subdivided into stars where the flux dips occur at an approximately constant period (designated as quasi-periodic systems) and stars with no obvious periodicity to the dips (designated aperiodic). In previous literature, these two subclasses, or portions of them, have been referred to as AA Tau systems and UX Ori systems, respectively; (3) short-duration bursters: stars with relatively well-defined minimum light curve levels, superposed on which are brief (hours to days) flux increases, attributed to accretion bursts. See Stauffer et al. (2014) for further discussion of this set; (4) quasi-periodic variables: stars lacking a well-defined maximum brightness but showing periodic variability whose waveform changes shape and amplitude from cycle to cycle; (5) stochastic: stars with prominent luminosity changes on a variety of timescales with no preference for "up" or "down." (6) long timescale variability: stars with slow (weeks to months) changes in brightness.

Although there are some clear similarities in the two schemes, their usage of the terms quasi-periodic and stochastic are not the same, and light curves described as belonging to those classes in one scheme would not necessarily be so classified in the other scheme. It will be important to resolve these nomenclature issues in the future, perhaps in the same manner as was done for sorting out similar issues on classifying pre-main-sequence disk SED morphologies (Evans et al. 2009). In the meantime, we must choose which scheme to adopt and be clear that that is what we have done. For this paper, we adopt the Wolk et al. (2013a) scheme.

Periodic light curves are described in Section 4, and examples can be seen in Figure 7. Examples of quasi-periodic light curves are shown in Figure 12. WL 3 shows variability with a timescale of ∼3 days in the first half of the first visibility window, but around MJD 55,316, the flux increases significantly and this rise masks out any periodicity. The dip around MJD 55,325 again has a duration similar to those observed in the beginning of the visibility window. The other two sources shown might be similar, but the signal-to-noise ratio is not as good.

Zoom In Zoom Out Reset image size

Other sources are variable, but no preferred timescale for the variability is discernible. ISO-Oph 50 and 51 (Figure 13) are examples of this aperiodic behavior. In yet other cases, the timescale of the variability is so long that we cannot decide if a feature is a singular event or part of a recurring pattern. Good examples of this are WL 4 and [EDJ2009] 892, which are shown in Figure 15.

Zoom In Zoom Out Reset image size

In some cases, we see short-duration features in the light curve in addition to a longer trend. This can either be a brightening (e.g., MJD 55,303 and 55,311 in ISO-Oph 51, Figure 13) or short dips in the light curve that last between 1 and 10 days ([EDJ2009] 809 in the same figure).

Like Wolk et al. (2013a), we find that even strongly periodic sources are not perfect clocks in the mid-IR in that the amplitude can vary from cycle to cycle. There might also be phase shifts from visibility window to visibility window. In a study very similar to the one we present here, Flaherty et al. (2013) observed IC 348 over a 40 day window with IRAC1 and 2 and find 25% of the variable stars are likely periodic. Here we find a similar division: 10–52 members are strongly periodic.

We present the first mid-IR monitoring of L1688, but the field has been monitored before in the near-IR. Figure 14 compares the K_s-band amplitude found in 2MASS data taken between 1997 and 2001 (Parks et al. 2014) and the K-band amplitude observed at UKIRT between 2005 and 2006 (Alves de Oliveira & Casali 2008) with the amplitude of our 3.6 μm light curves. Given the different bands, it is not surprising that the absolute value of the amplitude differs, but there is a good correlation such that the sources with the largest K- or K_s-band amplitudes in earlier observations also have the largest amplitudes in the IRAC bands 5–10 yr later. In most objects, color changes are small or happen in parallel with luminosity changes (Section 5.5), so the K or K_s band and the [3.6] variability should be strongly correlated. Figure 14 shows that the amplitude of the variability is relatively stable over at least one decade. Parks et al. (2014) used a larger number of observations than Alves de Oliveira & Casali (2008) and observed a larger spread between a typical YSO's brightest and faintest K-band magnitude, indicating that the longest timescale of variability is longer than the time span of the Alves de Oliveira & Casali (2008) observations.

Zoom In Zoom Out Reset image size

Flaherty et al. (2013) noted that roughly 6% of the stars in IC 348 had data that indicated significant changes in the source's mid-IR flux over the three-year interval since the last IRAC observation. Wolk et al. (2013a) also displayed numerous examples of stars that showed continuous change over the course of at least one observing visibility window.

All of this shows that the amplitude of variability of a given source may evolve over very long timescales but is consistent over at least one decade. YSOs apparently do not switch between highly variable and much less variable states.

For sources in the primary target fields, the observations in 3.6 μm and 4.5 μm are separated by only a few minutes. For those sources, we compare the color and the magnitude in the 3.6 μm band in a color–magnitude diagram (CMD, Figures 12 and 15). In some cases, the color is fairly stable within one visibility window but changes between visibility windows; in others the timescale of variability is much shorter, and color changes are seen within one visibility window as well. For nonvariable sources, the CMD forms a point cloud with the size set by the photometric uncertainties, but for variable sources, the shape of the CMD can reveal the physical cause of the variability; for example, if a disk warp or accretion funnel passes in front of a YSO, we expect its color to become redder as it becomes fainter. If the absorber has the same gas and dust properties as the interstellar medium (ISM) and absorbs star and inner disk at the same time, then the data points in the CMD follow a line with the slope of the interstellar reddening law.

Zoom In Zoom Out Reset image size

In each CMD, the observational data (I1, I1-I2) are fit to a line segment using an orthogonal distance regression method that takes the errors in both the x and y directions into account. We calculate the length in magnitudes of the line segment, excluding the 10% of the data that are outliers on either side. The CMD slopes of the 37 sources with a well-defined slope are shown in Figure 16. The x and the y axis in the CMD are correlated because they both depend on the 3.6 μm magnitude, so sources that are noise-dominated would show a slope of 45°. No source is seen in this region.

Zoom In Zoom Out Reset image size

From a total of 38 variable sources with CMDs with at least 10 data points, 34 can be fitted with a formal uncertainty on the best-fit slope below 6°. In three of the four sources with a statistical error on the slope in the CMD >6°, the color variability is dominated by observational uncertainties; the remaining source, ISO-Oph 140, is discussed below.

As can be seen in Figure 15, even in those sources with a well-defined slope in the CMD, the scatter around the best linear fit is larger than the measurement uncertainties; individual visibility windows are systematically above or below the fitted line.

The slopes shown in Figure 16 can be separated (somewhat arbitrarily) into two groups. One group becomes redder when the sources are fainter, with slopes comparable to the ISM reddening. This group contains almost all variables that are not classified as cluster members in Section 3.5. Most sources in this group have class I or flat-spectrum SEDs (right panel).

Sources that bluen (we use the term "bluen" as a verb to mean that a source becomes bluer similar to the common expressions "redden" or "reddening") when they are fainter are mostly class II and flat-spectrum sources. Figure 12 shows the light curves and CMDs for two of the bluening sources, WL 3 and [GY92] 264.

The left panel in Figure 16 quantifies the magnitude of the reddening or bluening observed in each source. For a source where the slope in the CMD is compatible with ISM reddening, this can be interpreted as the A_3.6 of the intervening material. A range of values is observed, but in most sources it is <0.3 mag. No significant difference in vector length is seen between sources that bluen or redden when they dim. We discuss the physical mechanisms that could cause these slopes in Section 5.8.

Comparing the amplitudes of all sources in Figure 16 that either redden or bluen over the full time span of the observations, we do not find significant differences in value or distribution of the variability amplitudes.

Similar to our results, Alves de Oliveira & Casali (2008) find reddening and bluening slopes in the CMD in their HK monitoring, but compared to Figure 16 they see a higher fraction of stars that bluen as they dim (33/49). Only six sources have well-defined slopes in both their study and ours. In the data of Alves de Oliveira & Casali (2008), all of them belong to the group that bluen with lower fluxes. In five cases, our data agree. In contrast, ROXN 44 reddens in our observations, so its behavior changed over the time interval of 5 yr, or it is behaving differently in the near and mid-IR.

ISO-Oph 140 is the only source with a statistical error on the slope in the CMD >6° where the color variability is not dominated by observational uncertainties but where the slope varies substantially over time. Its light curve and CMD for the first visibility window are shown in Figure 17. It has a class II SED slope and is a low-mass YSO (spectral type M1; Luhman & Rieke 1999). There is no single slope in the CMD; instead, the source seems to switch between different modes. Initially, the luminosity drops sharply, while the color becomes slightly bluer. The behavior changes around MJD 55,305. For the next 15 days, the source becomes redder and dimmer, but the slope is less steep than the ISM reddening law. The evolution is not monotonic, but a brightening (e.g., around MJD 55,313) corresponds to a bluer color, and the source reddens as it dims. At the end of this period, the brightness increases again sharply, and the slope in the CMD is comparable to the period MJD 55,295–55,305: the source becomes noticeably redder with higher luminosity. In the last 10 days of the monitoring, the source again has a slope close to the ISM reddening law. Its luminosity continues to increase, and the color bluens until it has a color and luminosity similar to what it had at the beginning of the monitoring. This closes the circle in the CMD. This is the only source in our sample where we observe multiple changes of the reddening slope over time. In this case, variability similar to ISM reddening corresponds to slow changes in luminosity, while the other direction of the slope in the CMD is associated with faster luminosity changes.

Zoom In Zoom Out Reset image size

Now we discuss the light curves with the largest changes in magnitude in our sample. WL 4, shown in Figure 15, has an outburst event consisting of a brightness increase by about 0.6 mag and bluening at the same time. The rise and fall take only a few days, and the entire outburst lasts about a month.

[EDJ2009] 892 (shown in the same figure) presents a slow rise by about 0.2 mag in 4.5 μm and a smooth decay over 2 yr where the source gradually bluens (although more activity between visibility windows cannot be excluded). There is an indication of a slow and smooth increase in brightness again in the last observing visibility window, pointing to a recurring phenomenon. These timescales are much longer than the dynamical time in the inner disk, so they are presumably driven by disk phenomena that originate at larger radii.

Morales-Calderón et al. (2011) identified a few such bursting or fading events in their analysis of YSOVAR data for the ONC, and they have also been observed in the optical (Cody & Hillenbrand 2010; Findeisen et al. 2013). The physical cause is not known yet. Findeisen et al. (2013) find a few light curves that match theoretical predictions for short mass accretion events, but without simultaneous spectroscopy this is hard to prove.

The best example presented here, WL 4, becomes bluer during the burst, but the opposite happens for WL 3, which otherwise has a similar light curve. Thus, it is likely that different mechanisms cause these burst events.

The reddening sources in Figure 16 have slopes that are roughly compatible with the slope of an interstellar reddening law. Small deviations can be explained by modifications of the dust-to-gas ratio or the grain size and composition. This is commonly observed in individual young stars (Günther & Schmitt 2008) and star-forming regions as a whole (Winston et al. 2010; Günther et al. 2012).

The most prominent example of a class II source where the extinction changes on timescales of days is AA Tau. In this case, the absorption dip is periodic and caused by an inner disk warp partially occulting the star (Bouvier et al. 1999; Schmitt & Robrade 2007). Many examples of comparable light curves have been found in the ONC by Morales-Calderón et al. (2011) and in NGC 2264 by Cody et al. (2014). These authors call objects with absorption events in the light curve "dippers."

We observe reddening on timescales of months to years in about half of all sources with a well-defined slope in the CMD. An example of that is shown in the top right panel of Figure 15. We do not know the inclination of the objects in our sample, but if it is close to edge-on, this might be caused by an asymmetric disk similar to AA Tau. The required warp or local change in scale height that lets some gas and dust protrude above the average disk height could be caused by a low-mass star, brown dwarf, or planetary mass object orbiting in the disk (e.g., Uribe et al. 2011; Baruteau & Masset 2013). Alternatively, vortices in the disk (e.g., Lesur & Papaloizou 2010; Lin 2012) or dust traps (Pinilla et al. 2013) can also change the local scale height.

Changes in the optical brightness and color of a YSO are often attributed to spots, either hot accretion spots or cold magnetic spots, rotating in and out of view or to time-variable extinction (e.g., Carpenter et al. 2001, 2002). If the reddening increases and, at the same time, the visible fraction of the hot spot decreases, then the resulting slope in the CMD would be intermediate between ISM reddening and colorless.

Because spots on the star cannot explain sources that bluen as they dim, this observed phenomenon must be related to the disk. Several parameters of the inner disk could be time variable. The most obvious one is the accretion rate. However, the direct effect of increased accretion is a larger accretion spot, which would make the star bluer and brighter, not fainter. Also, Faesi et al. (2012) find no relation between mid-IR light curves and spectroscopic accretion tracers in their limited sample. This indicates that the relevant parameter can change without affecting the accretion rate. In the mid-IR, we see the optically thick dust at the inner disk edge. Because the disk is heated over its entire vertical height at the inner edge while only the disk surface absorbs radiation at larger radii, it can form a puffed-up inner rim that reaches above the usual disk scale height and casts a shadow on the remaining disk, reducing the luminosity at longer wavelengths (e.g., Wisniewski et al. 2008; Espaillat et al. 2010, 2011). If this rim grows, it intercepts more stellar light. Thus, the emission from the inner wall increases and the emission from the outer disk decreases, which results in bluer colors. Because the total emitting area also decreases, the source can become bluer and fainter. Kesseli et al. (2013) also present radiative transfer models with hot spots and warped disks that reproduce features of the observed light curves.

Alternatively, Flaherty et al. (2013) suggested a nonaxisymmetric model. If the emission from the stellar photosphere is stronger at some longitudes due to a strong accretion spot or a spot from magnetic activity on the star, it irradiates the inner disk like a searchlight beam and causes the dust at this longitude to retreat. The result is a disk where the inner rim is broken at one (or several) positions. This would reduce the emission in 3.6 μm but cause more energy to be emitted at a larger radius, where the dust is cooler on average and thus emits at longer wavelengths such as 4.5 μm. In this way, the source can become redder as it brightens.

There is one group of YSOs, the so-called UX Ori (UXOr) variables, where several members have been observed to bluen when they fade in optical observations (Bibo & The 1990; Waters & Waelkens 1998). This can be explained by a larger fraction of scattered light that contributes to the observed SED in their faint state. A similar scenario for the YSOs discussed here requires, first, scattering by relatively large grains (>1 mm) because only for large grains will radiation at 3.6 μm be scattered more than radiation at 4.5 μm (Draine & Lee 1984; Andersen et al. 2013). Grains this large will be present in the disk, but not at all radii and all disk heights. Second, for this mechanism to work, a large fraction of the IR radiation must be scattered light and not intrinsic emission. In contrast, many bluening sources are class I to flat-spectrum sources, where the IR luminosity is comparable to the total stellar irradiation.

None of the simple models presented so far offers a convincing explanation of the observed light curves and CMDs. However, detailed simulations that take into account turbulent transport processes in the disk and different dust species with a complete chemical network have not yet been performed. The real situation is likely to be much more complex than sketched above. Because different dust species form and sublimate with different speeds and at different temperatures, time-variable irradiation can actually cause a very complex mixture of species with different opacities, which will not react linearly to increased irradiation. Such a complex network might cause a delayed and nonlinear response of the disk, which masks the underlying relation between accretion and disk emission.

The characteristics of the time variability in the light curves could be related to the spectral type of the central star because hotter stars irradiate the disk with harder spectra. We compared the 90% quantiles of the observed 4.5 μm magnitudes for those sources in our sample that have spectroscopically determined spectral types from Erickson et al. (2011). We use the relation between spectral type and T_eff from Pecaut & Mamajek (2013). This gives values for ten objects in our sample, and we do not find a correlation in this set. Similarly, Cody et al. (2014) do not identify any trends in IR variability versus effective temperature among NGC 2264 stars that are a few megayears old. In apparent contrast, Flaherty et al. (2013) find that the variability in [3.6] and [4.5] increases with increasing T_eff, but they show that this trend probably does not reflect a change in disk properties, but is instead due to the lower relative contribution of the photosphere compared to the disk for hotter stars.

In this section, we analyze the subsample of sources with X-ray counterparts. Because the subsample is not selected for its IR properties, it provides a clean sample for calculating the variability fraction. There are 31 sources with detected X-ray emission, of which 20 are variable in the IR. Again, we find high variability fractions in class I (6/7), F (3/4), and II (11/15), and none of the four class III sources is variable. The sample size is small, but this is consistent with our analysis of the IR-selected sample in Section 5.1.

We searched for correlations between the parameters of the X-ray emission (median energy, fitted temperature, and absorbing column density) and the amplitude of the mid-IR variability. No such correlation is apparent within each SED class. For a given absorbing column density, soft X-ray emission is more strongly absorbed than hard X-ray emission. Thus, the observed X-ray spectrum from more embedded sources always appears harder. Because class II sources have more circumstellar matter than class III sources, it is not surprising that they appear harder. Imanishi et al. (2001a) already discuss this for the L1688 X-ray data; Figures 9 and 10 show that class II sources are more variable than class III sources. We see no additional effect of the X-ray emission on the variability characteristics in the mid-IR.