THE SPITZER SPACE TELESCOPE SURVEY OF THE ORION A AND B MOLECULAR CLOUDS. I. A CENSUS OF DUSTY YOUNG STELLAR OBJECTS AND A STUDY OF THEIR MID-INFRARED VARIABILITY

The IRAC observations were taken as part of the Guaranteed Time Observation (GTO) programs PID 43, 50, 30641, and 50070. Table 1 contains a complete summary of the IRAC observations. The maps were made using a triangular grid of pointing centers in celestial coordinates. Since this grid was fixed to the celestial sphere, the map could be adjusted to follow the contours of the molecular cloud. The triangular grid was designed so that the observations could be executed at the full range of possible spacecraft orientations without holes in the coverage. Each field was observed in two epochs scheduled during consecutive visibility periods; this ensured that the detector arrays were rotated by as much as 180° between the epochs. The rotation between the two epochs effectively swapped the positions of the two spatially offset IRAC field of views (FOVs; i.e., the 3.6/5.8 μm and the 4.5/8 μm FOVs), thus providing full four-band coverage over most of the survey area. The rotation also facilitated the removal of bright source artifacts in the data.

The 12 s high dynamic range (HDR) mode was used; this produced two frames with integration times of 0.4 s and 10.4 s (corresponding to frame times of 0.6 and 12 s). At every map position, two dithered observations were obtained during each epoch. Thus, for most positions, the total integration time is at least 41.6 s per band for the long HDR frames and 1.6 s for the short HDR frames. However, due to the rotation of the two FOVs between the two epochs, many positions at the edges of the map have an integration time of 20.8 s per band for the long HDR frames and 0.8 s for the short HDR frames.

In Figure 1, the surveyed FOVs are shown overlaid on an extinction map of the Orion clouds from Gutermuth et al. (2011). The Orion cloud complex consists of two distinct giant molecular clouds, Orion A and Orion B, as shown in ¹²CO and ¹³CO surveys of this region (Bally et al. 1987; Miesch & Bally 1994; Wilson et al. 2005). The Orion A cloud was surveyed in one contiguous 4.04 deg² map extending from NGC 1977 down through the L1641 dark cloud. The survey field was constructed to follow the distribution of ¹³CO emission in the Bell Labs survey of Orion A (Bally et al. 1987). The survey was later extended to fully map the ONC. Subsequent extinction maps of the Orion clouds showed the Orion A cloud extending to a declination of −10° (κ Ori extension in Figure 1; Carpenter 2000); the survey was expanded again to include this southern extension of the cloud. In addition to the L1641 dark cloud, the 5.9 deg² map of Orion A also contains the L1640 dark cloud (part of the ONC region) and the L1647 dark cloud (which we include as part of the L1641 region). The Orion B map was broken into three fields centered on the regions of strong ¹³CO emission (Miesch & Bally 1994): a 1.27 deg² field covering the NGC 2024 and NGC 2023 nebulae in L1630, a 1.1 deg² field covering the NGC 2068 and NGC 2071 nebulae in L1630, and a 0.36 deg² field covering a region of bright ¹³CO emission in between the NGC 2068/2071 and NGC 2024/2023 fields (hereafter Orion B filament). In addition, a 0.23 deg² field covering the dark cloud L1622 was observed; in this paper, we consider L1622 to be part of the Orion B cloud (Bally et al. 2009).

Zoom In Zoom Out Reset image size

In program PID 50, the NGC 2024 nebula and Orion Nebula were observed with 1.2 s integration time frames (corresponding to a frame time of 2 s) to obtain high-signal-to-noise ratios in regions toward which the 10.4 s frames were saturated by bright nebulosity. Toward the Orion Nebula, the map consisted of a 3 × 4 rectangular grid in spacecraft coordinates with 270'' offsets between pointing centers. Toward NGC 2024, a 2 × 3 rectangular grid in spacecraft coordinates was made with the same offsets. In both maps, 24 dithered 2 s frames were obtained per map position for total integration times of 28.8 s per band.

A total of eight reference fields were also observed; these covered a total of 0.41 deg² in all four bands (Figure 1). For each field, a single epoch observation with four dither positions was executed; the typical integration time was 41.6 s per band. Five of the fields were observed in 2005 October (PID 43); these consist of four fields on either side of the ONC and one field targeting a small cloud east of Orion A detected by Bally et al. (1987). The remaining fields, all of which are found to the east of the cloud complex, were observed in 2007 October (PID 30641).

The MIPS observations were obtained from GTO programs PID 43, 47, 58, 30641, and 50070. The observations are tabulated in Table 2 and displayed in Figure 1. The Orion A cloud was mapped in six distinct fast scan maps: one covering the ONC (Table 2), one extending the ONC map to the north, and three covering L1641 and one covering the part of the cloud adjacent in the sky to the star κ Ori. The cross-scan steps were 148'' or 160'', resulting in a typical integration time of 30 s per pixel. The combined map covered 11.43 deg². The Orion B cloud was mapped using four individual scan maps. The NGC 2068/2071 region was mapped using medium-rate scans and 160'' cross-scan offsets, giving a typical integration time of 80 s per pixel. The NGC 2024/2023 region and Orion B filament were mapped with fast scans and either 148'' or 160'' cross-scan steps; the resulting integration time was 30 s. The result was two contiguous maps: a 2.58 deg² map including NGC 2068/2071 and the filament, and a 2.04 deg² map of NGC 2024/2023. Finally, a field covering 0.31 deg² toward L1622 cloud was mapped using medium-rate scans and 302'' cross-scan offsets with a typical integration time of 40 s.

The data were mosaicked using the custom IDL package Cluster Grinder (Gutermuth et al. 2009). The input images were basic calibrated data (BCD) data processed with the Spitzer Science Pipeline. The BCD data pipeline versions included s14, s17, and s18. Cluster Grinder first prepares the BCD images for mosaicking by matching the background of adjacent images, executing cosmic-ray masking using overlapping images, and ameliorating artifacts due to column pull-down, banding, muxbleed and residual bandwidth effects (see the chapter on IRAC in the Spitzer Observation Manual for a discussion of these artifacts).

The cosmic-ray rejection is done by comparing overlapping pixels in the four or more frames covering a given mosaic pixel. Pixels which differ by more than 10 median absolute deviations from the median value are masked out of the mosaic. On the outskirts of the mosaics there are sometimes only two overlapping pixels due to the FOV rotation between epochs; in these cases we increase the redundancy by including multiple pixels from each overlapping image. Specifically, we include pixels that are in the same row or column and are directly adjacent to each overlapping pixel resulting in 10 pixels used in the cosmic-ray rejection for each mosaic pixel. Pixels that deviate by more than 20 median absolute deviations from the median value of the 10 pixels are then masked. For the 1.2 s imaging of the ONC and NGC 2024, the high level of redundancy will translate into improved cosmic-ray rejection compared to the longer integration time mosaics.

The data units were converted from the native BCD units of MJy sr⁻¹ into DN by dividing by the conversion factor given in the BCD header ("FLUXCONV") and then multiplying by the integration time. Cluster Grinder then uses the refined WCS information in the BCD image headers to register the images to the mosaic grids and interpolate the BCD images onto the mosaic grid, which had an identical pixel scale to the BCD data. The mosaic pixels were given the value of the average interpolated BCD images weighted by a bad pixel mask.

For the data taken in HDR mode, two mosaics were produced for each wavelength band: one from the 10.4 s integration time data and the other from the 0.4 s integration time data. These are then combined into a single mosaic by replacing all saturated pixels in the 10.4 s mosaic with the pixel values from the 0.4 s mosaic scaled to a 10.4 s integration time. For the 1.2 s data, a single mosaic was produced for each wavelength band.

The point sources on the HDR combined mosaics were identified using PhotVis v1.10 (Gutermuth et al. 2008). PhotVis convolved the mosaics with a sunken Gaussian with an FWHM similar to that of a point source (2 pixels). It then identified all sources in the convolved image where the peak was seven times the rms of the pixel values in a 10 × 10 region centered on the peak. Since regions with bright structured nebulosity have higher rms values, this approach dramatically reduced the number of nebular knots misidentified as stars. The point sources identified on each mosaic were then inspected by eye using PhotVis and all obvious artifacts misidentified as point sources remaining in the data were removed at this stage. In addition, a small number of sources that were missed by PhotVis due to confusion with nebulosity were added. Photometric magnitudes were extracted with PhotVis using a 2 pixel radius aperture and an annulus from 2 to 6 pixels to estimate the sky contribution. The zero points were 19.6642, 18.9276, 16.8468, 17.3909 (for units of DN s⁻¹) and aperture corrections were 1.213, 1.234, 1.379, and 1.584 for the 3.6, 4.5, 5.8, and 8 μm bands, respectively (Reach et al. 2005). The uncertainties for each of the four IRAC bands were calculated by the IDL Astronomy Library routine aper.pro by combining the estimated shot noise from the source, the uncertainty in the background estimation, and the measured standard deviation in the sky annulus (PhotVis incorporates routines from the IDL Astronomy Library of Landsman 1993).

At this point, an initial point source catalog was generated using the bandmerging routine within Cluster Grinder. The routine employs the following procedure. First, each 3.6 μm point source was associated with the nearest point source in the longer wavelength bands as long as they were separated by ⩽1''. If the nearest point source was beyond 1'', no association was made and the source would have a detection only in the 3.6 μm band. This procedure was repeated in turn for the 4.5, 5.8, and 8 μm point sources that were not previously merged with a shorter wavelength band. The point sources were then associated with the nearest 2MASS point source as long as the 2MASS source was within a separation of 12. The final R.A. and declination pair chosen was that given by the first wavelength band with an uncertainty ⩽0.25 mag in the following progression: 4.5 μm, 3.6 μm, 8 μm, 5.8 μm, and 2MASS. The end product was a catalog of point sources with R.A., declination, photometric magnitudes for all seven bands, and the uncertainties in all seven bands.

Photometric magnitudes of stars observed at multiple positions on flat-fielded BCD images show variations of 5% (see IRAC chapter of Spitzer Observing Manual). To correct for these variations and to search for photometric variability, final photometric magnitudes were then derived from the individual BCD images using the coordinates in the point source catalog derived from the mosaics. These final magnitudes were obtained for all point sources identified in the mosaics with uncertainties ⩽0.25 mag in one of the four IRAC bands. The steps of this process were the following. First, each frame containing a given source was identified, divided by the conversion factor of MJy sr⁻¹ over DN s⁻¹ (i.e., "FLUXCONV"), and then multiplied by the integration time to convert the units into DN. These frames were then corrected for column pull-down, banding, muxbleed, and residual bandwidth effects using custom IDL programs developed for Cluster Grinder. Finally, the R.A. and declination pair of each source, as determined from the mosaics, was converted to pixel coordinates in the BCD images using the refined WCS information in the BCD headers. Aperture photometry was extracted at these positions using aper.pro with an aperture of 2 pixels and a sky annulus of 2–6 pixels. The zero points and aperture corrections were identical to those used for the mosaics. The photometry was then further corrected using the array location dependent photometric corrections given by the Spitzer Science Center;¹¹ these corrections account-for the 5% variation in point source photometry observed across the BCD images. We note that point response function (hereafter PRF) fitting was not used due to the difficulty in taking into account the variations of the IRAC PRF across the focal plane. (The PRF is the point-spread function convolved with the pixels response function.)

The photometry was extracted independently from both the 10.4 s and 0.4 s frames. Thus, each point source was typically measured eight times (four dithers and two integration times), although in many cases the sources were not detected at the shorter integration time. For the ONC and NGC 2024 observations with 1.2 s frames, each point source was measured independently for all 24 dithers. The photometry of each point source was collated as a function of Julian date for all three integration times. Magnitudes were calculated for all sources that were detected in two or more images of the same integration time. For each combination of source, wavelength band, and integration time, the median magnitude was tabulated using the average of the middle two magnitudes when an even number of photometry values were present (if there were only two photometry values, the average of those two values was used). The use of the median value eliminates the impact of cosmic-ray strikes as long as there were three or more dithers and only one dither out of the three was affected. Each of the median magnitudes was assigned an uncertainty equal to (∑^N_{i = 0}σ_i²/N²)^1/2 where σ_i was the uncertainty of the ith magnitude and N was the total number of magnitudes for a given source, band, and frame time combination. To construct the final table of photometry, the 10.4, 1.2, and 0.4 s data were compared, and the highest signal-to-noise photometry achieved without saturation was adopted. The integration time at which a given point source would saturate was determined from the peak pixel value of that point source in the HDR combined mosaics.

Comparison of the photometry determined directly from the mosaics with that measured from the individual BCD frames showed a systematic offset. The median values of the differences of the magnitudes, m_mosaic − m_BCD, were 0.02, 0.03, and 0.03 mag for the 3.6, 4.5, and 8 μm bands, but was −0.05 mag for the 5.8 μm band. This resulted in slightly red [5.8] − [8] colors even for pure, Vega-like photospheres. To remove this offset, we corrected the magnitudes in the 5.8 μm band using the average difference between the 5.8 μm band offset and the offsets in the other bands. The correction was applied with the equation m(5.8 μm) = m_BCD(5.8 μm) − 0.077 mag.

A final photometric catalog of IRAC and 2MASS photometry was then created by replacing the photometry from the mosaics with the median value of the individual frame photometry. Note that a photometric magnitude in a given band would be replaced only if the uncertainty for the mosaic photometry was ⩽0.25 mag in that band. Otherwise, the magnitude for that band was set to −100 and the uncertainty to 10. Although the 0.25 mag limit might be considered only a 4σ detection, all of the sources were previously found by PhotVis by showing a peak pixel value 7σ above the noise. Furthermore, we apply more stringent limits to the photometry when we search for YSOs.

The MIPS 24 μm data were reduced, calibrated, and mosaicked using the MIPS instrument team's Data Analysis Tool (Gordon et al. 2005). The resulting MIPS mosaics have a pixel size of 1245 and units of DN s⁻¹ det pix⁻¹ (where det pix is a detector pixel). Point source identification and aperture photometry was first performed with PhotVis. As was done for the IRAC data, PhotVis convolved the mosaics with a sunken Gaussian having a 4.5 mosaic pixel FWHM and then identified all sources in the convolved image where the peak pixel value was seven times the rms of the pixel values in a 10 × 10 region centered on the peak. An aperture of 5 pixels (6225) and a sky annulus extending from 12 (1494) to 15 (18675) pixels were used. We adopted a zero-point magnitude of 16.05; this was calculated using a mosaic pixel that was half the diameter of a detector pixel, a calibration factor of 6.4 × 10⁻⁶ Jy per DN s⁻¹ det pix⁻¹ (again, det pix is a detector pixel), and a zero magnitude flux in the MIPS 24 μm band of 7.17 Jy. In addition, an aperture correction was included in the zero-point magnitude, which consisted of two first parts. First, we multiplied the signal in the aperture by 1.146; this is the correction from a 12 pixel aperture to infinity (Engelbracht et al. 2007). Furthermore, we applied the aperture correction between the 5 pixel aperture used in our photometry and the 12 pixel aperture. This was determined by comparing the photometry in 5 and 12 pixel apertures for a sample of bright, isolated, and unsaturated stars in our sample; the resulting correction was −0.428 mag.

Given that that the MIPS PSF has an FWHM more than twice that of the IRAC PRF, it was necessary to revise the photometry using PSF-fitting photometry to minimize the effect of confusion with other sources and nebulosity. The extraction of the MIPS photometry using PSF-fitting is described in Kryukova et al. (2012). The PSF photometry was calibrated using the aperture photometry obtained for the stars used to construct the reference PSF. Thus, it will have the same calibration as the aperture photometry described above.

There are two regions of spatially extended saturation in the MIPS 24 μm survey where point sources could not be detected. The largest region of saturation is found in the Orion Nebula at the center of the ONC. The saturation encompasses the Trapezium as well as the region surrounding BN/KL. This region has been recently mapped at 6–37 μm with FORCAST on SOFIA by Shuping et al. (2012). The other region of extended saturation is located in the NGC 2024 nebula and contains the exciting star of NGC 2024, IRS 2b (Bik et al. 2003).

Finally, the MIPS point sources were merged with the combined 2MASS and IRAC point source catalog. The R.A. and declination of point sources in the IRAC point source catalog were compared to the positions of the MIPS 24 μm point sources. The MIPS magnitudes and uncertainties were then assigned to the nearest IRAC point sources if they were separated by ⩽25. If there was no IRAC photometry within this limit, the MIPS coordinates, magnitudes, and uncertainties were added to the point source catalog as a new source.

The final point source catalog contains 298,405 sources that are detected in at least one of the IRAC bands or in the MIPS 24 μm band. The magnitude uncertainties for the point source photometry are displayed in Figure 2. The uncertainties in all eight bands follow the expected upward sloping curve with the uncertainties increasing with increasing magnitude; however, the uncertainties in the IRAC bands show a large scatter for a given magnitude. The reason for this spread is that the uncertainties depend not only on the magnitude of the source, but also their location. Sources located toward bright, structured nebulosity or in regions crowded with other point sources show higher uncertainties than sources with equal magnitudes in regions without bright nebulosity or crowding.

Zoom In Zoom Out Reset image size

To assess the dependence of the measured uncertainties on the presence of nebulosity and nearby point sources, we determine the root median square deviation (hereafter RMEDSQ) of the pixel values in an annulus centered on the point source. The RMEDSQ for a given star is given by

$$\begin{equation} {\rm RMEDSQ}(i_0,j_0) = \sqrt{{\rm median}[(S_{ij}-{\rm median}[S_{ij}])^2]}, \end{equation} \tag{ 1 }$$

where i₀ and j₀ are the pixel coordinates of the source, and i, j are the pixels that are found in an annulus centered on the source. For the IRAC photometry, the annulus typically extends from 6 to 11 pixels (72–132). For bright stars, fluctuations in the PRF may dominate the RMEDSQ in this annulus. To ensure that the PRF itself does not contribute to the RMEDSQ, the radius of the annulus is extended to larger radii until variations in the PRF make a negligible contribution to the RMEDSQ. The maximum allowed size of an annulus is from 45 to 50 pixels (54''–60''). Since the units of our mosaics are DN, the resulting RMEDSQ are in units of DN. For the MIPS photometry, the annulus also extends from 6 to 11 pixels (747–13695). Again, the radius of the annulus is extended for bright stars to larger radii until variations in the PSF are negligible with a maximum allowed size of 45–50 pixels (56025–6225). The units for the MIPS data are DN s⁻¹ det pix⁻¹.

The resulting RMEDSQ is a measurement of the spatial varying signal surrounding each point source. These spatial variations are due to highly structured nebulosity, nearby point sources, and background photon noise. To simplify the analysis, we do not distinguish between the source of the variations and include any variation in the signal in the RMEDSQ regardless of the cause. The RMEDSQ is similar to a standard deviation except that it uses median values instead of mean values. Consequently, RMEDSQ is less sensitive to pixels with extreme values than the standard deviation calculated from the same pixels. This moderates the influence of nearby stars, which only affect a small region of the annulus surrounding the star. Note that the RMEDSQ gives the same value as the median absolute deviation.

In Figure 2, we show the uncertainties of the point sources segregated by their RMEDSQ values in the 8 μm band. The 8 μm band is used since the spatially structured nebulosity apparent in all four of the IRAC bands is strongest in that band. By measuring the RMEDSQ in the 8 μm band, we are examining the dependence of the uncertainties on the nebulosity. In all four of the displayed IRAC bands, there is a clear progression to higher uncertainties with increasing RMEDSQ. The same progression is seen in the MIPS 24 μm band, although the spread in uncertainties for a given magnitude is much lower than that observed in the IRAC bands. In contrast, the 2MASS bands do not show the same progression. The reason is that the bright nebulosity that typically dominates the RMEDSQ values at wavelengths >3 μm is much fainter relative to stellar sources in the 1–2.5 μm bands. Thus, the near-IR J, H, and K_s bands are not strongly affected by nebulosity.

The catalog exhibits spatially varying completeness due to confusion with nebulosity and, to a lesser extent, the crowding of point sources in dense clusters and groups. An assessment of the completeness of the point source catalog must take into account this spatially varying background. Instead of measuring completeness for each position in the IRAC maps, we determine completeness as a function of RMEDSQ. The completeness for a given source magnitude at any position in the Orion survey can then be determined by measuring the RMEDSQ at that position.

As described in the Appendix, we added artificial stars to the field containing the ONC and then recovered these stars using the same techniques described in the previous sections. For the IRAC bands and the MIPS 24 μm band, we determined the fraction of artificial stars recovered as a function of the stellar magnitude and the RMEDSQ in the selected wavelength band. The fraction of recovered stars versus RMEDSQ are shown for different source magnitudes in Figure 3 (displaying the four IRAC bands) and Figure 4 (showing the MIPS 24 μm band). These show a strong dependence of completeness on both magnitude and the RMEDSQ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The capability of Spitzer to identify YSOs with infrared excesses from dusty disks and envelopes is now well established in the literature (e.g., Allen et al. 2004, 2007; Megeath et al. 2004; Gutermuth et al. 2004, 2008; Muzerolle et al. 2004; Whitney et al. 2004; Hartmann et al. 2005; Winston et al. 2007; Muench et al. 2007; Robitaille et al. 2007). We have employed the methodology developed in Gutermuth et al. (2009) whereby YSOs exhibiting infrared excesses are identified and classified using color–color diagrams constructed from the combined 2MASS, IRAC, and MIPS photometry. We have augmented this approach by using criteria for identifying and classifying protostars first developed in Megeath et al. (2009) and further refined in Kryukova et al. (2012). We have also adopted new criteria to identify a small number of YSOs that were missed in the criteria of Gutermuth et al. (2009) and Kryukova et al. (2012).

To identify and classify YSOs, we utilized all photometric magnitudes with uncertainties less than 0.1 mag for the J, H, K_s, 3.6, and 4.5 μm bands and 0.15 mag for the 5.8 and 8 bands. We did not apply an uncertainty cutoff for the 24 μm band, which in cases of saturated sources that are fit to the wings of PSFs may exceed 0.2 mag (Kryukova et al. 2012). The sources satisfying these criteria are displayed in the color–color and color–magnitude diagrams shown in Figures 5–7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.7.1. The Reddening Law

2.7.2. YSO Identification

The Gutermuth et al. (2009) methodology is divided into multiple phases designed to classify sources with the most sensitive diagram for the available photometric bands. We first applied the Phase I criteria to sources with detections in all four IRAC bands, 3.6, 4.5, 5.8, and 8 μm, as described in Appendix A.1 of Gutermuth et al. (2009). These criteria isolate sources with excess emission in the Spitzer bands relative to a pure stellar photosphere, and they reject probable active galactic nuclei (AGNs), galaxies with bright polycyclic aromatic hydrocarbon (PAH) emission, sources contaminated by structured PAH emission, or knots of shock-heated emission. The selection process is displayed in Figure 5. Sources lacking detections in all four IRAC bands, but with magnitudes in the H, K_s, 3.6 μm, and 4.5 μm bands were then applied to the Phase II criteria, as shown in Appendix A.2 of Gutermuth et al. (2009; see Figure 5 of this work). Again, we excluded all sources with colors and magnitudes consistent with galaxies, bright PAH emission, AGN, or knots of shock-heated emission from the Phase II sample using the criteria in Gutermuth et al. (2009). Note that in contrast to Gutermuth et al. (2009), we group all of the Class 0/I, flat spectrum sources and Class II sources (i.e., all protostars and pre-main-sequence stars with disks) together as YSOs with IR-excesses. The classification of the YSOs is performed later during the protostar selection (see Section 2.7.3).

We identified additional YSOs by applying criteria designed to identify YSOs that are not detected in the wavelength bands combinations required by the criteria of Gutermuth et al. (2009). The first of these criteria examines stars detected only in the J, H, and 4.5 μm bands. We followed the approach of Gutermuth et al. (2004) and Winston et al. (2007) and identified as YSOs sources satisfying the inequality:

$$\begin{equation} H-[4.5]+\sigma _{H-[4.5]} \ge 0.921 (J-H-0.55+\sigma _{J-H}). \end{equation} \tag{ 2 }$$

This selection criterion is illustrated in Figure 6. In addition, for sources with detections in the 3.6, 4.5, and 5.8 μm bands that are missing 8 μm band detections or detections in the H and K bands, we identified YSOs that satisfy the inequalities:

$$\begin{equation} [3.6]-[4.5] \ge 0.5 + \sigma _{[3.6]-[4.5]}, [4.5]-[5.8] \ge 0.25 + \sigma _{[4.5]-[5.8]}. \end{equation} \tag{ 3 }$$

These criteria are displayed in Figure 5. To reduce the contamination from background galaxies, we further required that YSOs identified by the above two criteria satisfy the inequality:

$$\begin{equation} [4.5]-A_{4.5} \le 14, \end{equation} \tag{ 4 }$$

where A_4.5 is the total cloud extinction determined from the A_V map.

The final phase in Gutermuth et al. (2009) is the identification of IR-excess sources through their MIPS photometry. At this step, we diverged from the Gutermuth et al. (2009) criteria, and hence we describe our criteria explicitly. We began by identifying sources with weak excesses in the IRAC bands, but a strong excess in the MIPS 24 μm band. Such objects include transition disks with large inner disk holes (Muzerolle et al. 2010). We took all sources that satisfied our uncertainty limits in the four IRAC bands, but which did not show an IR-excess in the IRAC bands. We then identified as YSOs sources that satisfied the condition:

$$\begin{equation} [8]-[24] > 2.5 + \sigma _{[8]-[24]}, \end{equation} \tag{ 5 }$$

where the 2.5 mag threshold comes from Gutermuth et al. (2009). The identified sources can be seen in Figures 5–7. Note that we only mark sources as transition disk in these figures if they have detections in all four IRAC bands, have no apparent IR-excess in the IRAC bands, and their dereddened [3.6] − [4.5] colors are less than 0.2 mag. Next, we identified sources that showed strong excesses at 24 μm but were not identified in all four IRAC bands. These sources satisfied one of the following four inequalities:

$$\begin{eqnarray} &&[3.6]-[24] \ge 2.5\quad {\rm or} \quad [4.5]-[24] \ge 2.5\nonumber\\ && \quad {\rm or} \quad [5.8]-[24] \ge 2.5\ {\rm or}\ [8]-[24] \ge 2.5. \end{eqnarray} \tag{ 6 }$$

Again, the 2.5 mag thresholds come from Gutermuth et al. (2009). For a magnitude to be used in one of the above inequalities, it must satisfy the uncertainty limits for that band. In addition, we require that the sources satisfy the criterion:

$$\begin{equation} [24] \le 7. \end{equation} \tag{ 7 }$$

This criterion rejects the large population of faint, red 24 μm sources apparent in Figure 7. Unlike the YSOs, which are concentrated within the boundaries of the cloud, these faint sources are evenly dispersed throughout the mapped regions and are most likely extragalactic sources (Kryukova et al. 2012; Megeath et al. 2009). Note that the [24] ⩽ 7 inequality is not required for the transition disk criterion in Equation (5). The reason is that the transition disk criterion is applied only to the sources with detections in all four IRAC bands that have already been subject to criteria designed to eliminate extragalactic sources.

After applying all the additional criteria described above, we further reduced contamination from galaxies by eliminating all sources with colors and magnitudes consistent with galaxies, bright PAH emission, AGN, or knots of shock-heated emission using the criteria of Gutermuth et al. (2009).

2.7.3. Protostar Identification

The protostar identification serves two purposes: finding protostellar YSOs that were not found using the criteria overviewed in the previous section, and classifying the likely protostars among the entire sample of YSOs. We implemented the criteria of Kryukova et al. (2012), who used colors constructed from the IRAC and MIPS bands to find sources with flat or rising spectra between the 4.5, 5.8, or 8 μm IRAC bands and the 24 μm MIPS band. These criteria were derived from the spectral index criteria of Greene et al. (1994), with the adopted color thresholds equivalent to a spectral index criterion α ⩾ −0.3, where α = dlog (λF_λ)/dlog (λ). To minimize extragalactic contamination, we adopted the magnitude cutoff that Kryukova et al. (2012) determined from the Orion cloud such that all protostars must satisfy [24] ⩽ 7.

In regions showing bright nebulosity or saturation in the 24 μm bands, we must use alternative criteria for identifying protostars. Following Gutermuth et al. (2009), we use two criteria based on the IRAC and 2MASS photometry. First, we identify sources with detections in the 3.6, 4.5, and 5.8 μm bands guided by the criteria of (Gutermuth et al. 2009)

$$\begin{equation} [3.6]-[4.5] \ge 0.7 + \sigma _{[3.6]-[4.5]}, [4.5]-[5.8] \ge 0.7 + \sigma _{[4.5]-[5.8]}. \end{equation} \tag{ 8 }$$

To minimize contamination to the protostar sample, all sources identified as AGN, galaxies with PAH emission, outflow shock knots, and stars contaminated by PAH emission (as described in Appendix A.1 from Gutermuth et al. 2009) are eliminated. Furthermore, we also require that the protostars satisfy the criterion in Equation (4). For sources without 5.8 μm and 24 μm detections, we use the Phase II criteria from Gutermuth et al. (2009) to find protostars. These criteria require photometry in the H, K_s, 3.6 μm, and 4.5 μm bands.

Once the potential protostars without 24 μm detections are identified, we reject sources where the upper limit to the [4.5] − [24] is <4.761 mag and the lower limit of the 24 μm magnitude is >7 mag. We determine these limits using the 24 μm magnitude at which 90% of the sources are recovered for the local value of RMEDSQ (Figure 4). In areas where the 24 μm data are saturated, we do not reject any sources.

As shown in Figure 7, there are sources that clearly exhibit the red colors of protostars and satisfy all the Phase I criteria of Gutermuth et al. (2009); these sources should be considered dusty YSOs and are unlikely to be contamination. Although they satisfy the color criteria for protostars, they ultimately fail the Kryukova et al. (2012) criteria because they have [24] > 7. We classify these YSOs as faint candidate protostars.

A number of objects are also detected that show MIPS detections with [24] < 7, but no detections in the 4.5, 5.8, or 8 μm bands. To distinguish between protostars and sources such as asteroids or galaxies, these extremely red sources were visually inspected (asteroids would not be detected at the same position by IRAC due to their orbital motion). If they appeared as a distinct point source in a region with high extinction, if they were found in a region of active star formation, or if they were associated with jets apparent in the IRAC images, they were classified as red candidate protostars. A total of 10 protostars were added to the catalog. Far-IR observations with the PACS instrument on board Herschel will be used to confirm their protostellar nature (A. M. Stutz et al., in preparation).

2.7.4. Distinguishing YSO Candidates from Outflow Knots

We classified the YSOs as either protostars, faint candidate protostars, red candidate protostars, or pre-main-sequence star with disks. The protostars and candidate protostars were identified by the criteria in Section 2.7.3. Any YSO with an IR-excess that was not classified as a protostar or candidate protostar was classified as a pre-main-sequence star with a disk. We have not distinguished between transition disks and disks with optically thick inner regions in our catalog; these can be identified using the criteria described in Section 2.7 or with criteria proposed in the literature (e.g., Muzerolle et al. 2010; Cieza et al. 2010; Merín et al. 2010).

Highly reddened stars with disks may be misidentified as protostars; this can occur when a star with a disk is observed through a dense region of the molecular cloud or when a star is observed through a nearly edge-on disk. Kryukova et al. (2012) applied the approach of Gutermuth et al. (2009) to put upper limits on the number of young stars misidentified as protostars due to their nearly edge-on disks, and they used Monte Carlo simulations to estimate the number of young stars with disks that are misidentified as protostars due to extinction from the molecular cloud. They estimated that at most 24% of the protostars are actually misclassified stars with disks. The ongoing Herschel Orion Protostar survey (Fischer et al. 2010; W. J. Fischer et al., in preparation) will provide more stringent limits on this percentage by using Herschel 70 and 160 μm observations to confirm the detection of infalling envelopes around Spitzer-identified protostars.

The total statistics of the number of sources, their classification, and the method used to identify them is given in Table 3. Note that red candidate protostars are not tabulated as IR-excess sources; since they are only detected in one to two bands, we cannot establish whether they have excesses using our criteria. In the table, the number of additional sources is given in each column, starting with the leftmost column and progressing to the right.

The final YSO catalog is given in Table 4. This catalog contains all of the identified Spitzer dusty YSOs, their 2MASS photometry, their photometry in all four IRAC bands, their photometry in the MIPS 24 μm band, and their classification. The uncertainties given are the formal uncertainties from the extraction of the photometry. In addition to these formal uncertainties, the IRAC and MIPS instruments have absolute calibration uncertainties of 0.05 mag.¹² Since the absolute uncertainties affect all the sources equally, we tabulate only the formal uncertainties which varies from source to source. For consistency with the tables in young cluster survey of Gutermuth et al. (2009), we also include the extinction toward the source and the spectral index of the source over the four IRAC bands. The determination of those values is described in Gutermuth et al. (2009).

Background galaxies, and in particular AGN, often show colors similar to those of YSOs (Stern et al. 2005; Gutermuth et al. 2008). Although the YSO criteria described above are designed to minimize contamination due to galaxies, a residual contamination is expected. To assess the residual contamination of galaxies, we take two approaches. First, we use data from the SWIRE Spitzer Legacy program, a wide-field survey of galaxies executed with the IRAC and MIPS instruments (Lonsdale et al. 2003). Second, we use the reference fields observed as part of the Orion molecular cloud survey.

To estimate the contamination from the SWIRE survey, we used the IRAC and MIPS photometry from a circular region in the Elias N2 field with a total area of 2 deg². We merged these data with the J-, H-, and K_s-band photometry from the 2MASS point source catalog. To facilitate a comparison between the Orion survey and the deeper SWIRE data, we first "degrade" the uncertainties of the SWIRE data so that they are comparable to those achieved with the Orion survey. Using the Orion photometry, an empirical curve of uncertainty versus magnitude is computed by calculating the median uncertainty as a function of magnitude for each band. For each point source in the SWIRE survey, the tabulated uncertainty was replaced by the median Orion uncertainty for the equivalent magnitude and wavelength band. We then applied our methodology for identifying YSOs to the SWIRE data. From this analysis, we find a total of 14.5 ± 2.7 galaxies are misidentified as YSOs per square degree. Of these, 2.5 ± 1.1 are misidentified as protostars and 4.0 ± 1.4 are misidentified as faint candidate protostars.

We also imaged eight off-cloud reference fields with the IRAC instrument; the overlapping regions of 3.6 and 4.5 μm cover a total of 029 (Figure 1, Table 1). In these fields, we find a total of three IR-excess objects, for a density of 10.4 ± 6.0 objects per square degree. In addition to these fields, there is one cloud field, the Orion B "filament," in which we detected only a single YSO (Figure 12). We include this 0.36 deg² map as a reference field since it was covered by MIPS and IRAC and since it provides a measurement of the contamination toward the cloud with its higher extinction. If we assume that all the IR-excess sources in the combined reference fields are objects misidentified as dusty YSOs, we find the total contamination is 6.1 ± 3.1 objects per square degree, of which 1.5 ± 1.5 are misidentified as protostars and 2.8 ± 1.5 are misidentified as faint candidate protostars. Note that the density of IR-excess sources in the off-cloud fields is consistent with the density of galaxies misidentified as YSOs in the SWIRE field. However, when we include the reference field toward the Orion B cloud, the contamination drops significantly; this may result from the higher extinction through the on-cloud region. We adopt the value for the combined reference fields, 6.1 ± 3.1 objects per square degree misidentified as YSOs, as the typical source contamination since it is more representative of on-cloud regions targeted in this survey. We note this is similar to the contamination of eight objects per square degree quoted by Gutermuth et al. (2009).

Another form of contaminants are stars misclassified as dusty YSOs due to their coincidence with bright nebulosity. Aperture photometry of point sources superimposed on structured nebulosity may result in some of the nebulosity contributing to the observed magnitudes and creating a fake IR-excess. To assess the number of such misidentified stars, we first identified every star in the survey field with sufficient photometry to detect an IR-excess but without the detection of an IR-excess. From this sample of stars, we excluded likely extragalactic contamination leaving a total of 27,339 sources.

We then added an artificial star 10'' east of every one of the identified stars. The artificial stars closely reproduce the spatial distribution of stars without IR-excesses, which are primarily field stars, and thus eliminate potential biases induced by the adopted spatial distribution. For example, this experiment will correctly account for the fact that most of the field stars are found in dark regions without bright nebulosity. The colors and magnitudes of the artificial stars were taken from the set of stars in the reference fields which also had sufficient photometry to detect an excess, but were not identified as IR-excess sources or galaxies. We ignore the MIPS data due to the small number of additional YSOs identified solely with MIPS and since the MIPS photometry was extracted with PSF-fitting that may reduce, but not eliminate, the effect of nebulosity. Furthermore, all but one of the adopted reference fields used to supply the colors and magnitudes of the artificial stars lacked MIPS photometry.

We next recovered the photometry of the artificial stars. The input and recovered photometry is shown in Figure 8. The IRAC color–color diagrams show a significantly larger spread in the [5.8] − [8.0] color of the recovered data due to the nebulosity in the image. In comparison, the increase in the spread of the [3.6] − [4.5] color is negligible. This can be shown quantitatively by examining the standard deviation; while the standard deviation of the [3.6] − [4.5] color changed from 0.05 mag to 0.07 mag, the standard deviation of the [5.8] − [8] color changed from 0.14 mag to 0.19 mag. Furthermore, a number of objects are scattered to [5.8] − [8] colors exceeding 0.5 mag but maintain small [3.6] − [4.5] colors. This diagram demonstrates the necessity of rejecting candidate YSOs with IR-excesses only in the [5.8] − [8.0] color. Out of the 27,399 artificial stars, only 13 showed infrared excess using our criteria for the IRAC color–color diagram. We note that 11 of those 13 stars have recovered 4.5 μm magnitudes that are more than 0.5 mag brighter than the input magnitudes; these may be artificial sources that have been added on top of actual sources. We conclude that the contribution to our YSO sample by stars with spurious IR-excesses due to confusion with nebulosity and stars is negligible.

Zoom In Zoom Out Reset image size

Since the Orion survey observations were not designed to measure variability, a search for variable sources requires the extraction of photometry from one or two individual data frames with no redundancy to eliminate cosmic rays. To reliably identify variable sources, we use the Stetson index which measures correlated variability in the 3.6 and 4.5 μm bands of IRAC (Stetson 1996). These bands were chosen as they are the least affected by nebulosity and typically show the highest signal-to-noise ratios. The requirement of correlated variability minimizes the chance of mistaking a cosmic-ray strike during one of the epochs as variability.

We have adapted the Stetson analysis outlined in the 2MASS variability study of Carpenter et al. (2001) to the IRAC 3.6 and 4.5 μm bands. The Stetson index S was calculated using the equation:

$$\begin{equation} S = \frac{1}{N_{{\rm epoch}}}\sum ^{N_{{\rm epoch}}}_{k=1} {\rm sgn}(P_k)\sqrt{|P_k|}, \end{equation} \tag{ 9 }$$

where N_epoch is the number of independent epochs and

$$\begin{equation} P_k = \delta _k(3.6 \,\mu {\rm m})\delta _k(4.5 \,\mu {\rm m}). \end{equation} \tag{ 10 }$$

The variable δ_K is calculated for each epoch from for the 3.6 and 4.5 μm magnitudes using the equation:

$$\begin{equation} \delta _k(\lambda) = \sqrt{\frac{n_{{\rm epoch}}}{n_{{\rm epoch}}-1}}\left(\frac{m_{\lambda }-\bar{m}_{\lambda }}{\sigma _\lambda } \right), \end{equation} \tag{ 11 }$$

where $\bar{m}_{\lambda }$ is the average of m_λ over the epochs. In the determination of S, we only used magnitudes with uncertainties less than 0.1 mag and we required that this uncertainty limit was met in at least two epochs. To increase the signal-to-noise ratio and reliability of the photometry, we average observations which were taken within two days of each other. This will smooth out any short-period variability over the two-day interval. The number of overlapping observations can range from one frame (in the case of L1641 where there are separate second and third epochs), two frames (in the case of most of the other fields), to as high as six frames (in the interstices of our mapping grid where frames from three separate positions overlap). In cases where the peak signal of a source in the combined mosaics exceeded the saturation limit of the 10.4 s frame or the magnitude of the short 0.4 s frame was brighter than the 10.4 s frame saturation magnitude, we used the photometry for the 0.4 s frame. A noise floor was set by adopting a minimum value of 0.05 mag for σ_λ.

The distribution of Stetson indices for the entire cloud sample are displayed in Figure 17. We compare the distribution of YSOs with IR-excesses with that for stars that show no evidence of IR-excesses (hereafter pure photospheres). The pure photospheres are a mixture of field stars in the line of sight as well as diskless pre-main-sequence stars. The difference between these two distributions is clearly apparent. While the IR-excess sources show a strong bias toward positive values of S, the pure photospheres are mostly confined to small values of S centered on 0. We show the histograms for sources with m_4.5 ⩽ 14; this magnitude threshold eliminates the sharp increase of S with fainter magnitudes apparent for the pure photospheres. We assume the distribution of S for the photospheres is dominated by noise and systematic uncertainties in the photometry and not by intrinsic variability in the sources. Given the that the distribution is symmetric with S = 0, this appears to be a reasonable assumption. We can then compare this distribution to that for the YSOs with IR-excesses to assess the frequency of variability in both pre-main-sequence stars with disks and protostars. This comparison shows that the YSOs with IR-excesses are clearly biased to positive values of S, indicating the presence of significant variability.

Zoom In Zoom Out Reset image size

In Figure 18, we compare the variability as a function of evolution. (Note: protostars do not include faint candidate protostars in this analysis.) To identify a sample of diskless pre-main-sequence stars, we use the 2MASS variability survey to find sources that show variability in the near-IR bands but show no evidence for a disk in the Spitzer wavelength bands. These objects are thought to be pre-main-sequence stars with large star spots that rotate in and out of view (Carpenter et al. 2001). We use all sources that are not found to show an IR-excess in the Spitzer bands (despite having photometry in a sufficient number of wavelength bands to detect an excess) and which have near-IR Stetson indices ⩾0.55 in the observations of Carpenter et al. (2001).

Zoom In Zoom Out Reset image size

The histograms show a clear progression in the distribution of S with evolution, with the protostars showing the highest values of S, the young stars with disks exhibiting intermediate values, and the diskless young stars showing typically very small values. We note that the IR-excess sources exhibit much higher values of S even though the diskless pre-main-sequence stars were selected by their inherent variability. Nevertheless, the diskless pre-main-sequence stars show a higher incidence of variability than the total sample of pure photospheric sources that includes both field stars and diskless pre-main-sequence objects.

To better quantify the fraction of YSOs that exhibit variability, we tabulate the fraction of sources with S > nσ_S(photo), where σ_S(photo) is the standard deviation in the S for the photospheres (Table 5). We show that the IR-excess sources have systematically larger values than the photospheres. If we use the photospheres as a reference distribution for non-variable sources (i.e., the spread in the values of S result from random and systematic errors), we estimate that ∼50% of the YSOs with IR-excesses exhibit variability.

We now characterize the changes in magnitude during the variations. For the IR-excess sources where S > 2σ_S(photo), Figure 19 displays the change in the magnitude in each of the IRAC bands. In each case, we identified the epoch that showed the minimum and maximum 4.5 μm magnitudes. We then plotted the difference in magnitude for each IRAC band between those two epochs. In each case, we plot the difference in magnitude relative to that observed in the 4.5 μm band. We find that the distribution of maximum deviations is strongly correlated and approximately follows a diagonal in each plot. This indicates that the change in magnitude is approximately the same in all four wave bands. Note that there are several points where variability is seen in one band but not the other. In these cases, the photometry may have been affected by cosmic-ray strikes and are removed from the analysis. Nevertheless, there is a strong correlation between the independent bands: the Pearson's r correlation coefficients are given in Table 6. These show strong correlations in all three cases which can be further strengthened by limiting the analysis to sources showing changes ⩾0.25 mag.

Zoom In Zoom Out Reset image size

We fit a line to each correlation using the functional form $[\lambda ]-[\lambda ]_{4.5_{{\rm min}}} = A + B \times ([4.5]-[4.5]_{4.5_{{\rm min}}})$ after eliminating all the points that are likely cosmic-ray strikes from the fit (Figure 19). We try fits that include all the data points (after rejecting the possible cosmic-ray strikes) and fits limited to those data points showing changes in both bands greater than 0.25 mag. By ignoring sources that show small variations, we reduce the sensitivity of the fits to spatially varying nebulosity. The results of the fits are given in Table 6. In three out of four of the observed correlations, the slopes are within 3σ of unity, indicating that the variability does not change the color of the sources. The only exception is the slope for the Δ3.6 μm versus Δ4.5 μm correlation, which shows a somewhat lower slope.

The distribution of the magnitude variations is shown for all four IRAC bands in Figure 20. We plot the maximum variation in magnitude for each source where S > 2σ_S(photo). Although the variations can be as large as 1.76 mag, most of the variations are smaller than 0.5 mag. The maximum variations in each band as a function of percentiles are given in Table. 7. In addition to having a higher incidence of variability, the protostars exhibit larger variations in magnitudes than the young stars with disks.

Zoom In Zoom Out Reset image size

We find evidence of variability for 50% of the YSOs with IR-excesses. The magnitude changes are approximately the same in all four IRAC bands and are typically on the order of 0.2 mag. Both the fraction of sources with variability and the changes in magnitude are largest for protostars. Since diskless objects do not exhibit similar variability, it is likely that the observed variability is due to the changes in the emission from the inner disk region where dust emits strongly in the IRAC bands. This emission is thought to originate primarily in the inner rim of the disk; at this inner rim the central stars heats the dust to the sublimation temperature of silicate grains. Dust within this radius will be sublimated, thus eliminating the primary source of opacity in the disk. Since the disk is vertically extended, the emission in the inner rim originates come from an "opacity wall" facing the central star. The infrared emission from the disk traced by the IRAC bands arises mainly in the rim, but it includes contributions from disk radii up to 1 AU, particularly in the 8 μm band (D'Alessio et al. 2006).

The variations in overall disk flux could be due to changes in the heating of the disk, fluctuations in its physical structure, or some combination of the two. The wavelength dependence of the variability, as measured by the slope of Δm_λ versus Δ[4.5], and the timescale of the variations provide constraints on the mechanism for the infrared variability. A physical model for the variability must explain how material over the inner 1 AU could uniformly change its flux over the course of six months. Variations on timescales of less than one year could be due to fluctuations and rotating structures in the inner 1 AU of the disk where the dynamical timescale ranges from months to days with decreasing radius. These variations may result from a variety of reasons such as a rotating disk warp or magneto-rotational instability instabilities that launch dust out of the midplane (Turner et al. 2010). Previous observational studies have invoked scale height perturbations to the dust at the very inner edge of the disk to explain the wavelength dependence and timescale of the variability for evolved transition disks where the inner disk has been cleared. Scale height perturbations at the very inner edge of a disk result in a "seesaw" behavior for the infrared spectral energy distribution (SED), in which the short-wavelength flux increases while the long-wavelength flux decreases with the entire SED appearing to pivot at 8 μm (Flaherty et al. 2011; Espaillat et al. 2011). The observed variability could also result from changes in the heating of the inner disk by the central star. The disk is heated by both the intrinsic luminosity of the star and the accretion luminosity generated by gas funneled down onto the star. Changes in the heating of a disk could arise from fluctuations of the accretion rate onto the central star and the resultant change in accretion luminosity. The accretion luminosity is thought to be concentrated in accretion hot spots; the rotation of these hot spots may also result in variations in the heating of the portions of the inner disk most visible from the Earth (Bouvier & Bertout 1989; Morales-Calderón et al. 2009).

To examine how changes in the structure and heating of the inner disk would affect the observed magnitudes in the IRAC bands, we employ the CGPlus models (Chiang & Goldreich 1997; Dullemond et al. 2001). These models include a two-layer, flared disk with a puffed-up inner wall. While more complicated combinations of grain growth and settling are often needed to fit detailed SEDs, we can use these models to estimate the wavelength dependence due to changes in the disk and star. Our base model consists of a 0.5 M_☉, 2 L_☉, T_eff = 3800 K star surrounded by a 0.01 M_☉ disk that extends from the dust sublimation radius, defined as the point where T_dust = 1500 K, out to 100 AU. From these models, we can directly estimate magnitudes, and by varying different model parameters we can calculate the slope of Δm_λ versus Δ[4.5] in the same manner as was done with the actual data. First, we consider variations due to a puffed-up inner rim by increasing the scale height of the dust and gas of the inner rim. In this case, we find slopes of 1.014, 0.96, and 0.60 for [3.6], [5.8], and [8.0], respectively (Table 8). While the radiative transfer in the shadowed region behind the inner rim is not fully understood, changing the scale height in this way produces a "seesaw" behavior similar to what has been observed in transition disks and thus serves as a useful approximation. The calculated slopes are inconsistent with the observations, suggesting that in the majority of the YSOs, the variability is not due to a change in the scale height of the inner disk rim. We next vary the radial location of the inner rim of the disk by lowering the temperature at which the disk becomes truncated. In this situation, we find slopes of 0.53, 1.468, and 1.499 which are also highly inconsistent with the data (Table 8). This implies that, unlike in the evolved disks, we are not observing a change in the structure of the dust at the very inner edge of the disk. We cannot rule out synchronized changes in the disk scale height over the entire inner 1 AU, which would affect all the IRAC bands equally; however, it is not clear what mechanism could alter the disk in this manner.

Next, we consider changing the stellar flux in the CGPlus models without varying the flux of the disk. This produces slopes of 1.09, 0.861, and 0.687 for [3.6], [5.8], and [8.0], respectively; these values show the decreasing dominance of the stellar photosphere with increasing wavelength. The lack of agreement with the observed slopes confirms the need for a change in the disk flux. If we vary the stellar luminosity while allowing the disk structure to respond to this change in the impinging flux, we find slopes of 0.985, 1.012, and 1.013. If we vary the disk heating without varying the infrared emission from the photosphere, which would be the case if the disk was heated by UV radiation from the accretion of gas onto the star, the resulting slopes of 0.922, 1.053, 1.124 provide the best match to the data (Table 8). In this case, the reduced slope at [3.6] is due to the significant contribution of the photosphere in that band. The thermal timescale of the upper layers of the disk that are responsible for the observed infrared emission are very short (Chiang & Goldreich 1997) and the disk will respond almost instantaneously to changes in the illumination. This suggests that we are observing a change in the overall flux level of the disk by an amount that is relatively constant with wavelength. Variable heating may come from rapid changes in the rate of accretion of gas onto the central stars; the accretion luminosity contributes significantly to the total luminosity of protostars and young stars with disks. Alternatively, Morales-Calderón et al. (2009) constructed models of the variable heating in a protostar due to a rotating accretion hot spot. In this case, the variability is not due to a change in total luminosity, but due to the region of the inner disk rim illuminated by the hot spot rotating in and out of view. These models also produce a wavelength dependence that is similar to what has been observed. We note that the Morales-Calderón et al. (2009) models which add a disk warp rotating in sync with the hot spot predicts a weaker variability at [8.0]; this is inconsistent with our observations.

We conclude that the time variable heating of an extended portion of the disk driven by rotating accretion hot spots or rapid accretion variations is the most likely explanation for the observed variations in the mid-IR magnitudes. This is in distinct contrast to transition disks for which structural perturbations of the puffed-up wall at the dust sublimation radius can explain the observed mid-IR variability (Flaherty et al. 2011; Espaillat et al. 2011). Hence, the mechanism for mid-IR variability may change as young stars objects and their disks evolve. A more definitive assessment of the origin of the variability will require coordinated observations designed to study the time dependence of mid-IR emission, stellar luminosity, and accretion (e.g., Flaherty et al. 2012).

To assess the completeness of the data, artificial stars were added to the ONC field in the four IRAC bands and the MIPS 24 μm bands. As described in Section 2, the Orion data exhibit spatially varying completeness primarily due to confusion with the bright and highly structured nebulosity detected in all wavelength bands. Instead of measuring the completeness for each position in the IRAC and MIPS maps, we determine completeness as a function of RMEDSQ, as defined in Section 2. To do this, we used the ONC field since this field shows the largest range in RMEDSQ values.

For the IRAC artificial stars, we used a PRF that was reconstructed from the observations of an IRAC standard star.¹³ To account for the rotation of the spacecraft during the two epochs of the observation, the artificial PRF was the sum of two rotated PRFs. The artificial star analysis was done separately for "faint," "bright," and "ultra-bright" stars. The input magnitudes of the faint stars were distributed in 0.5 mag increments between 12.5 and 16.5 mag for the 3.6 and 4.5 μm bands, and between 10.5 and 14.5 mag for the 5.8 and 8 μm bands. Bright stars were added in 0.5 mag increments between 8 and 12 mag for the 3.6 and 4.5 μm bands, and between 6 and 10 mag for the 5.8 and 8 μm bands. Finally, the ultra-bright stars covered 3.5–7.5 mag in 0.5 mag increments for the 3.6 and 4.5 μm bands and 1.5–5.5 mag in 0.5 mag increments for the 5.8 and 8 μm bands. In all three cases, the stars were distributed in a rectangular grid in which each added star was separated by 125 from the adjacent artificial stars.

The added stars were then recovered with PhotVis using automated source detection; any source within 15 of the input position, within 0.25 mag of the input source magnitude, and with uncertainties within the limits set for YSO identification (0.1 mag for the 3.6 and 4.5 μm bands and 0.15 mag for the 5.8 and 8 μm bands) was considered recovered. The PhotVis object detection parameters were the same as those used for the point source catalog. The recovered faint stars were binned by their log(RMEDSQ) values using a bin of 0.5. The reason for the logarithmic binning is that increasingly bright nebulous regions cover progressively smaller portions of the image; a logarithmic binning thus provides better statistics for regions where RMEDSQ > 1000. We then calculated the fraction of recovered stars (n_rec) to added stars (n_add) using the equation f = n_rec/n_add. The fractions as a function of input magnitude and log(RMEDSQ) are shown for the IRAC bands in Figure 3.

For the purpose of the correction, it is of value to obtain for each magnitude a continuous curve parameterized by the RMEDSQ. For each input magnitude of the artificial stars, we fit a function to the fraction of recovered stars versus log(RMEDSQ). We fit the equation

$$\begin{equation} x = \log ({\rm RMEDSQ}), f = \frac{1}{2}\left(1-{\rm erf}\left(\frac{x-a}{\sqrt{2}b}\right)\right), \end{equation} \tag{ A1 }$$

where a and b are the value of RMEDSQ at which the completeness declines to 0.5 and the width of the decline, respectively. This equation has desirable features of having a fraction approaching one for small values of log(RMEDSQ) and a fraction approaching zero for high values of log(RMEDSQ). When fitting the data, we adopted an uncertainty

$$\begin{equation} \sigma _f^2 = (f(1-f)/n_{{\rm add}})^2+0.00001+0.000001(1-f), \end{equation} \tag{ A2 }$$

which is the uncertainty from Bolte (1989) with a noise floor added to both high and low values of f. The noise floors were set by using trial and error to find values of σ_f for which the fit would converge. The resulting fits are shown in Figure 3. The values of a and b for the selected magnitudes are shown in Table 9.

The MIPS 24 μm completeness was assessed in a similar way. A grid of point sources were added to the images; the point source was generated by STinyTim for a 50 K source. This was done with two iterations: one for faint stars with magnitudes ranging from 6.5 to 10.5 mag in increments of 0.5 mag and the other for bright stars with magnitudes ranging from 2 to 6 mag in increments of 0.5 mag. The stars were separated by 12.5 pixels. Point sources detected by PhotVis with positions and magnitude within 15 of the input position, 0.25 mag of the input magnitude and with an uncertainty ⩽0.25 mag were considered recovered. The PhotVis object detection parameters were the same as used for the point source catalog and the magnitudes were determined using aperture photometry using the aperture radius, sky annulus, and calibration described in Section 2.5. The number of recovered and added stars were binned in intervals of log(RMEDSQ), and the resulting fraction was fit repeating the steps for the IRAC data (Figure 4). The values of a and b for the MIPS data are shown in Table 9.