MID-INFRARED VARIABILITY FROM THE SPITZER DEEP WIDE-FIELD SURVEY

A substantial fraction of objects in the universe change in brightness with time. This apparent variability can be due to physical periodic changes in the objects (e.g., Cepheids, RR Lyrae), motion (e.g., eclipsing binaries, microlensing, rotating stars), explosive events (e.g., supernovae, novae), and accretion (e.g., cataclysmic variables, active galactic nuclei (AGNs)). In a high-latitude extragalactic survey, the most common variable sources are AGNs (e.g., Sesar et al. 2007). Optical variability is one method for identifying quasars²² (e.g., Eyer 2002; Geha et al. 2003; Rengstorf et al. 2004), although it is only recently that a fully quantitative approach to variability selection has been developed (Kozłowski et al. 2010).

The emission from quasars generally has three components. The UV to near-IR radiation is dominated by a hot accretion disk extending from an inner edge of a few gravitational radii from the black hole outward with, in simple thin disk theory, a temperature profile T ∝ R^−3/4 (Shakura & Sunyaev 1973). Near the inner edge, there is a corona of hotter gas that produces the non-thermal X-ray emission (Haardt et al. 1994). On scales where the temperature is below the dust sublimation temperature (∼2000 K), dust absorbs radiation from the disk and reradiates the energy in the mid-IR and far-IR (Barvainis 1987). The overall spectrum typically shows a minimum near 1 μm, with the emission from the disk rising toward the UV and the emission reradiated by dust rising toward the far-IR (e.g., Sanders et al. 1989). There is an increasing evidence that many physical relations between AGNs, galaxies, and their large-scale clustering have to be taken into account if their formation and evolution is to be understood (e.g., Tasse et al. 2008).

We know a great deal about the optical variability of quasars both from large studies of the variability seen in ensembles of sparsely monitored quasars and from detailed studies of individual quasars. Ensemble studies (e.g., Vanden Berk et al. 2004; de Vries et al. 2005) have shown that variability increases with decreasing optical wavelength, decreasing luminosity, and potentially decreasing black hole mass. The structure function of the ensemble variability is a power law with smaller variability amplitudes on short timescales, with some evidence for saturation on timescales of order a few decades. Until recently, there were few studies of individual quasars (e.g., Cutri et al. 1985; Clavel et al. 1989; Hook et al. 1994), but this has changed dramatically in the last year. The light curves of individual quasars are well modeled by a damped random walk, a stochastic process described by the amplitude of the random walk and a damping timescale for returning to the mean luminosity (Kelly et al. 2009; Kozłowski et al. 2010; MacLeod et al. 2010). While preliminary indications from Kelly et al. (2009), based on ∼100 quasars, suggest that these two process parameters are related to the quasar luminosity and black hole mass, MacLeod et al. (2010) used ∼9000 Sloan Digital Sky Survey (SDSS) Stripe 82 quasars to find a number of clear trends. For example, the asymptotic variability on long timescales decreases with increasing luminosity and rest-frame wavelength, and is correlated with black hole mass. The timescale for returning to the mean luminosity increases with wavelength and also with increasing black hole mass, but remains constant with redshift and luminosity.

Far less is known about the near-IR and mid-IR variability of quasars. Since AGN luminosities do vary in time at optical wavelengths, we expect to find some variability in the IR, but it could be heavily smoothed by averaging the response over the large scales of either the cooler parts of the disk or the absorbing dust. Glass (2004) found that the majority (39/41) of the surveyed low-luminosity nearby Seyferts do vary in the near-IR and mid-IR, and that the variability is most apparent at longer wavelengths due to diminishing flux from the host galaxy. The most extensive work to date on near-IR and mid-IR variability of quasars is that of Neugebauer & Matthews (1999). These authors studied 25 low-redshift quasars (z ≃ 0.1), which were observed in five bands (1.3–10.6 μm) for ∼30 yr. The amplitude of variations for radio-loud (quiet) quasars is ∼0.3 mag (0.1 mag) on rest-frame time lags (time differences between any two epochs) of 2 yr at 10.6 μm. For one quasar, PG 1226+023, they measured structure functions in five bands, all showing ∼0.1 mag variations for a rest-frame time lag of 2 yr. Another example is the work of Suganuma et al. (2006), who found a delayed response between the K- and V-band variability in four Seyfert 1 galaxies. These authors showed that K-band light curves have a fairly tight time-lag–luminosity relation, τ ≈ L^1/2, while broad-line region (BLR) time lags are more scattered. The interpretation is that the K-band emission is mainly coming from dusty clouds bounding the BLR, and that the size scale of these dust clouds is very dependent on the AGN luminosity, since they have to be colder than the dust sublimation temperature.

In this paper, for the first time, we investigate the variability of a truly large number of mid-IR sources, using the four-epoch Spitzer Deep Wide Field Survey (SDWFS; Ashby et al. 2009) of the Boötes field of the NOAO Deep Wide Field Survey (NDWFS; Jannuzi & Dey 1999). The two goals are to characterize the mid-IR variability of a large sample of AGNs and to use the variability to identify additional AGNs. While mid-IR quasar selection (Stern et al. 2005) works very well when the AGN dominates the luminosity (even for quasars behind dense stellar fields such as LMC/SMC; Kozłowski & Kochanek 2009), the method begins to fail as the luminosity becomes comparable to that of the host (Gorjian et al. 2008; Assef et al. 2010b). With mid-IR variability, we should be able to identify such sources and extend the sample of mid-IR selected quasars to lower luminosities. We also make use of the extensive wavelength coverage of the field, including X-ray (XBoötes; Murray et al. 2005), the earlier Infrared Array Camera (IRAC) Shallow Survey that became the first epoch of SDWFS (Eisenhardt et al. 2004), the 24 μm survey of the field (see Houck et al. 2005), and the FIRST (Becker et al. 1995) and WSRT (de Vries et al. 2002) radio surveys. Based on these data, the AGN and Galaxy Evolution Survey (AGES; C. S. Kochanek et al. 2010, in preparation) selected galaxies and AGNs at all these wavelengths for spectroscopic follow up using the Hectospec instrument on the MMT (Fabricant et al. 2005), providing roughly 23,500 redshifts in the field. In Section 2 we describe the data analysis, and in Section 3 we describe our variability selection criteria. In Section 4, we present the results of this study. The mid-IR structure function analysis is presented in Section 5 followed by a summary in Section 7. Throughout this paper, we use a standard ΛCDM model with (Ω_Λ, Ω_M, Ω_k) = (0.7, 0.3, 0.0) and h = H₀/100 = 0.70. Where needed, we use the galaxy and AGN templates of Assef et al. (2010b) to describe the mid-IR colors of various populations and to make any K-corrections or estimates of absolute luminosities.

The SDWFS is described in detail by Ashby et al. (2009). The images were taken with the IRAC instrument (Fazio et al. 2004) on the Spitzer Space Telescope, covering wavelengths of 3.6, 4.5, 5.8, and 8.0 μm, hereafter referred to as the [3.6], [4.5], [5.8], and [8.0] bands, respectively. A single IRAC frame is 256 × 256 pixels. We work with full-field SDWFS mosaics constructed from the ∼20,000 IRAC images/epoch. Each approximately 13600 × 16700 pixel mosaic covers 10.2 deg², centered at (R.A., decl.) = (14:31:27, +34:10:43), with a pixel scale of 0.84 arcsec. Counting the original IRAC Shallow Survey (Eisenhardt et al. 2004) as epoch 1, the SDWFS data span a period of four years between 2004 and 2008, where the intervals between the epochs are approximately 3.5 yr, 0.5 yr, and 1 month (observations were done on 2004 January 10–14, 2007 August 8–13, 2008 February 2–6, and 2008 March 6–10). For each individual observation of a patch of the sky (256 × 256 pixels), there is only a ∼40 s offset between the channel 1/3 and 2/4 observations. For a given SDWFS epoch, each patch of sky was observed ∼3 times over a 1–2 day period. Each SDWFS epoch of the entire 10 deg² area took ∼5 days, and we treat it as a single image with a single Julian Date. This does not affect our final conclusions, as we show later that there is little mid-IR AGN variability on these timescales.

To study the variability of objects in the SDWFS field, we used the Difference Image Analysis (DIA) method (Alard & Lupton 1998; Woźniak 2000). While the variability itself is the subject of this paper, the DIA method also proved to be a useful diagnostic tool for constructing artifact-free mosaics (see Ashby et al. 2009). The basic idea behind the DIA method is to match one image, called the template or reference image, both astrometrically and photometrically to another image of the same area of the sky. Once they are astrometrically aligned, the template image (which is chosen to be of better quality) is convolved with a kernel function to match the point-spread function (PSF) of the analyzed image. The two images are then subtracted from each other (Figure 1) to leave only variable objects, noise, and systematic problems (cosmic rays, artifacts, etc.). After some experimentation, we concluded that the [5.8] and [8.0] maps were too shallow to perform DIA and obtain reliable light curves, so our variability analysis is restricted to the [3.6] and [4.5] bands.

Zoom In Zoom Out Reset image size

Spitzer PSFs have a "triangular" shape with extended diffraction-like patterns (more precisely they are "multiplexer bleed" or "muxbleed," "column pull-down," and "banding"²³), and the PSF rotates between the epochs. DIA works best on simple quasi-Gaussian PSFs, so we pre-filtered the images to have a simpler PSF structure (Figure 1). For each region, we scaled and median combined the brightest stars to build an empirical model of the PSF, B. We then linearly filtered the image I to create a filtered image defined by the combination of Fourier transforms $\tilde{I}_F=\tilde{G}\tilde{I}/\tilde{B}$, that transforms the PSF from B to G, where G is a Gaussian broad enough to avoid noise amplification (i.e., $\tilde{G}\rightarrow 0$ faster than $\tilde{B}\rightarrow 0$ for small scales/large wave numbers). The filtering process conserves fluxes and was carried out on a 3 × 4 grid of 4040 × 4040 pixel regions that overlap by 20 pixels along their edges. To save some computation time, each of these 4040 × 4040 regions was then cut into smaller 2040 × 2040 pixel sub-images (again each with 20 pixels of margin) for the DIA analysis. We combined all four epochs to build the template images and these templates were scaled and subtracted from each epoch using DIA (Woźniak 2000). The brightest stars in the field produce artifacts due to subtraction residuals, particularly from the "diffraction" spikes in the PSF. We therefore masked a region around all stars with K < 11 mag in the Two Micron All Sky Survey (2MASS) catalog of point sources (Cutri et al. 2003). The masking radius was 260(1 − K/11) arcsec, which is large enough to include the spikes. We used the 2MASS catalogs because these bright stars were saturated and not included in the SDWFS catalogs. The masked regions constitute ∼3% of the area.

The [3.6], [4.5], [5.8], and [8.0] Vega magnitudes, and the coordinates (J2000.0) are taken from the SDWFS catalogs (Ashby et al. 2009). Throughout this paper, we use 4 arcsec diameter aperture magnitudes with aperture corrections to 24 arcsec. The original 3.6 μm full SDWFS catalog contains 677,522 objects (catalog DR1.1), but we cropped it to 515,055 objects (hereafter called the "cropped sample") as variations in coverage between epochs reduced the usable area for variability analysis by about 20%, to 8.1 deg². Using DIA's PSF fitting routines, we could construct light curves (at least three epochs) for 492,420 objects simultaneously in both [3.6] and [4.5]. The final number of objects used in this analysis after masking regions around bright stars is 474,179. While we include the results for all objects from the cropped sample in Table 1, in the analysis that follows we consider only the sources with at least three epochs. Table 2 provides the DIA light curves for all SDWFS sources from the cropped sample.

Table 1. SDWFS-DIA Variability Catalog

Object Name	R.A.	Decl.	[3.6]	[4.5]	[5.8]	[8.0]	N	v([3.6])	v([4.5])	C12	r	σ12	Mask
	(deg)	(deg)	(mag)	(mag)	(mag)	(mag)		(mag)	(mag)	(mag2)
SDWFS J143040.36+322721.8	217.668148	32.456048	17.666	17.488	16.792	16.070	4	0.051	0.062	0.001	0.178	0.894	1
SDWFS J143034.08+322729.0	217.642013	32.458058	17.937	17.466	17.142	16.928	4	0.067	0.185	0.002	0.127	0.758	1
SDWFS J143033.71+322730.3	217.640453	32.458420	18.725	18.352	17.576	16.869	4	0.181	0.223	0.029	0.718	0.392	1
SDWFS J143047.11+322729.5	217.696285	32.458194	17.121	16.995	16.898	14.735	4	0.036	0.096	0.002	0.511	0.295	1
SDWFS J142951.47+322726.0	217.464455	32.457230	18.361	17.875	17.538	16.447	4	0.135	0.102	0.009	0.627	0.734	1
SDWFS J143015.91+322733.8	217.566284	32.459402	19.732	20.200	18.262	99.999	1	99.999	99.999	99.999	99.999	99.999	1

Notes. This is an abridged version of the full table, presenting the variability statistics for the combined channels [3.6] and [4.5] only. The electronic table includes additional columns showing the variability statistics for the separate channels [3.6] and [4.5], and the photometric uncertainties for the IRAC channels. Here, quantities v([X]) and N are, respectively, the variability standard deviation in channel [X] and the common number of epochs in both channels used in the calculation of the joint variability. With C₁₂, r, and σ₁₂ we denote the covariance and the correlation between the channels 1 and 2, and the joint variability significance, respectively. The mask is 0 (1) if there is (is no) bright star with the 2MASS K-band magnitude within 260(1 − K/11) arcsec affecting the measurement. The Vega magnitudes presented in the table are measured in 4 arcsec apertures (see Ashby et al. 2009). The error code for magnitudes, reflecting no measurement, is 99.999. The error code for v, C₁₂, r, and σ₁₂ is 99.999.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

In order to understand the nature of the variable sources, we matched the SDWFS sources to several of the available data sets for the Boötes field, focusing on information related to AGN activity. Thus, we matched to the X-ray and radio catalogs, to the subset of 24 μm sources with high probabilities of being AGNs, and to the AGES spectroscopic data. We used the X-ray catalogs from the XBoötes Chandra survey covering 9.3 deg² of the Boötes field in the 0.5–7 keV band (Murray et al. 2005; Kenter et al. 2005; Brand et al. 2006; Gorjian et al. 2008). We matched a version of the Murray et al. (2005) catalogs extended downward in flux from 3293 sources with 4 or more counts to 4642 sources with 2 or more counts. While using 2 count (about 4 × 10⁻¹⁵ erg cm⁻² s⁻¹) X-ray sources sounds peculiar, the backgrounds in XBoötes are so low that almost all such sources with mid-IR identifications will be real, particularly if the match is to one of the relatively rare mid-IR variable sources. Using the Bayesian matching method of Brand et al. (2006), we identified 3363 X-ray sources with at least two epochs of SDWFS data. We used the same procedures to match to the 20 cm FIRST (Becker et al. 1995) and 1.4 GHz WSRT (de Vries et al. 2002) radio catalogs. Restricting the matching to likely radio point sources in the FIRST and WSRT surveys, we found 322 and 2118 matches, respectively, to SDWFS sources with at least two epochs of data. We know from AGES that most 24 μm, optical point sources satisfying the criterion that I>18 − 2.5log(F₂₄/mJy) are quasars (Brown et al. 2006; Eckart et al. 2010; Assef et al. 2010a), where the optical criterion is used to eliminate stars. Finally, we matched to the AGES redshift survey, finding 19,885 (18,288, 1387, 497) SDWFS sources (with redshifts z>0, 1, 2). The modified SDSS pipeline used for the AGES redshifts (HSRED; Cool et al. 2008) provides an estimate of the best-fit spectral template used for the fits that can be used to classify objects. More generally, we have all the galaxy and AGN target selection codes used by AGES to select and prioritize spectroscopic targets.

For each light curve i (from Table 2) with at least three epochs in both the [3.6] and [4.5] bands, we computed the standard deviation v([X])ⁱ for that light curve as a function of (Vega) apparent magnitude m([X])ⁱ, as shown in Figure 2. We also calculated the variability covariance Cⁱ₁₂ and Pearson's correlation rⁱ coefficients between the two channels, where

$$\begin{equation} v([X])^i = \left({\frac{1}{N_{\rm epochs}-1}\sum _{j=1}^{N_{\rm epochs}}\left(m([X])^i_j-\langle m([X])^i \rangle \right)^2}\right)^{1/2}, \end{equation} \tag{ 1 }$$

$$\begin{eqnarray} C_{12}^i &=& \frac{1}{N_{\rm epochs}-1}\sum _{j=1}^{N_{\rm epochs}}\left(m([3.6])^i_j-\langle m([3.6])^i \rangle \right) \nonumber \\ && \times \left(m([4.5])^i_j-\langle m([4.5])^i \rangle \right), \end{eqnarray} \tag{ 2 }$$

$$\begin{equation} r^i = \frac{C_{12}^i}{v([3.6])^i v([4.5])^i}, \end{equation} \tag{ 3 }$$

and 〈m([X])ⁱ〉 is the average magnitude for the ith light curve in channel X. The correlation coefficient is limited to the range −1 ⩽ r ⩽ 1, where r = 1 (−1) means that the [3.6] and [4.5] light curves are perfectly correlated (anti-correlated). The distribution of v([X])ⁱ is dominated by contributions from noise (and any systematic errors) rather than variability, so we need to define a threshold for selecting candidate variables from these distributions. Let v_m([X]) be the median of the v([X])ⁱ as a function of magnitude and σ([X]) be the standard deviation of the v([X])ⁱ values around this median. Sources are more likely to be true variables as their distance, (v([X])ⁱ − v_m([X]))/σ([X])>0, from the median increases, and if the [3.6] and [4.5] variabilities are strongly correlated. A reasonable metric for the joint significance of the variances in the light curves is

$$\begin{eqnarray} \sigma _{12}^i &=& \Bigg(\bigg({v([3.6])^i-v_m([3.6]) \over \sigma ([3.6])}\bigg)^2 \nonumber \\ && +\bigg({v([4.5])^i-v_m([4.5]) \over \sigma ([4.5])}\bigg)^2 \Bigg)^{1/2}, \end{eqnarray} \tag{ 4 }$$

which quantifies the degree to which the object deviates from the median variances in both bands. Objects with positive excess variances (large σ₁₂ and v>v_m) and strong correlations r ≃ 1 between the two bands are the best candidates for true variables.

Zoom In Zoom Out Reset image size

Figure 3 shows the distribution of correlation strengths r for sources with σ₁₂ exceeding a range of thresholds. We see that there is a peak for highly correlated sources r ≃ 1 whose strength increases as the variability significance σ₁₂ increases. This demonstrates that the variability in the two bands is tightly correlated for real variables and that we are likely to produce a relatively clean sample of variables (the false positive rate is 6%–7% for AGNs in the Photometry Group 3 defined below) if we restrict ourselves to r>0.8 and σ₁₂>2. Variability selection is always a tradeoff between completeness and contamination, so we will consider sources in three different variability levels, all restricted to have correlations r>0.8 between the bands. To avoid confusion, the variability levels correspond directly to the values of σ₁₂. Variability Levels 2, 3, and 4 correspond to sources with excess variances and significances exceeding σ₁₂>2, 3, and 4, respectively, with Variability Level 2 being the most complete and the most contaminated, and Variability Level 4 being the least complete and least contaminated.

Zoom In Zoom Out Reset image size

We estimate the fraction of false positives by assuming that the density of uncorrelated sources (r < 0.5) in Figure 3 represents the density of false positives amidst the highly correlated (r>0.8) sources. If the number of identified variables is N_r>0.8, that is sources with σ₁₂>2 (or 3, 4) and r>0.8, and the number of sources with σ₁₂>2 (or 3, 4) and r < 0.5 is N_r<0.5, then we estimate that the number of false positives is F = (0.2/1.5)N_r<0.5, where 0.2/1.5 is the ratio of the intervals r>0.8 and −1 ⩽ r ⩽ 0.5. The fraction of the identified variables that are false positives is then F/N_r>0.8.

Our ability to characterize sources also depends on the magnitudes of the object, because this determines how well we can measure the mid-IR colors. We divide sources into three photometry groups based on their mid-IR magnitudes.

• 1.  
  Photometry Group 1 consists of all objects that after masking have three or four epochs of data, with no constraints on either their magnitudes or magnitude uncertainties. There are 474,179 Photometry Group 1 objects after masking, of which 5107, 2208, and 1071 meet the Variability Level 2, 3, and 4 criteria.
• 2.  
  Photometry Group 2 is the subset of the Photometry Group 1 sources for which we can measure the [3.6]−[4.5] color. These sources must satisfy [3.6] < 19.7 mag and [4.5] < 19.3 mag with uncertainties in these bands smaller than 0.1 mag. There are 213,594 Photometry Group 2 objects, of which 1557, 668, and 359 meet the Variability Level 2, 3, and 4 criteria.
• 3.  
  Photometry Group 3 is the subset of Photometry Group 2 objects for which we can also measure the [5.8]−[8.0] color. These sources are required to have uncertainties smaller than <0.2 mag for the 5.8 and 8.0 μm bands. There are 39,522 Photometry Group 3 objects, of which 775, 428, and 260 meet the Variability Level 2, 3, and 4 criteria.

Table 1 provides these significance statistics for the SDWFS sources. Note, however, that our catalogs are limited to the objects detected on the template image—any transient which does not produce a detectable source in the template image, which is just the average of four epochs, will not be identified. The most variable, highly correlated (r>0.8) sources in the field are AGNs, peaking at σ₁₂ ≈ 20.

We will make extensive use of the [3.6]−[4.5] versus [5.8]−[8.0] mid-IR color–color distribution to characterize the sources. Figure 4 shows the distribution of all Photometry Group 3 SDWFS sources in this space. Eisenhardt et al. (2004) suggested that for bright sources the red plume in [3.6]−[4.5] color was likely to be dominated by AGNs. This was confirmed by Stern et al. (2005) based on the first year of spectroscopic observations by the AGES survey. Stern et al. (2005) defined a selection region, the "AGN wedge," for mid-IR AGN selection that is also shown in Figure 4. For these Vega colors, stars lie at color 0 unless they have dusty winds or disks to produce a mid-IR excess. Using the spectral energy distribution (SED) models of Assef et al. (2010b), which are derived by fitting the SEDs of Boötes sources, we can illustrate the locations of various extragalactic populations as a function of redshift. Early-type galaxies at low redshifts have colors similar to stars and then become redder in [3.6]−[4.5] for z>1 at roughly fixed [5.8]−[8.0], tracking along the left edge of the AGN wedge. They become red in [3.6]−[4.5] due to the flux "peak" near rest-frame wavelength of 1.6 μm, which passes through the 4.5 μm band at a redshift of 1.5–2.1. Late-type galaxies at low redshift have red [5.8]−[8.0] colors at low redshift due to the polycyclic aromatic hydrocarbon (PAH) emission in the 8 μm band, then move horizontally to bluer [5.8]−[8.0] along the bottom of the AGN wedge as the feature is redshifted out of the 8 μm band, and finally follow the early-type galaxies up along the left edge of the selection region. Pure AGNs remain inside the AGN wedge except near z ∼ 4.5 where they make an excursion out of the AGN wedge when the Hα line lies in the 3.6 μm band. If we use an equal combination of an AGN and a host galaxy, defined by the total emission from 0.1 to 30 μm, the combined SED is somewhat bluer. With further increases in the host luminosity relative to the AGN, the source drops below the bottom of the AGN wedge.

Zoom In Zoom Out Reset image size

In fact, the distribution of the IRAC colors for the XBoötes X-ray sources does extend below the Stern et al. (2005) AGN wedge (Gorjian et al. 2008). Based on Gorjian et al. (2008), we define a "modified AGN wedge" that encompasses most X-ray sources. This modified AGN wedge is defined by [3.6]–[4.5] < 2.5 × ([5.8]–[8.0]) − 0.6 mag (left side), [3.6]–[4.5]>2.5 × ([5.8]–[8.0]) − 3.5 mag (right side), and [3.6]–[4.5]>0.2 × ([5.8]–[8.0]) + 0.08 mag (bottom). Because we are identifying AGN candidates based on their mid-IR variability, we should be able to identify them in this [3.6]–[4.5] blueward extension of the Stern et al. (2005) selection region despite any host contamination. The caveat, however, is that this blueward extension is dominated by false positives (galaxies) if used to select AGNs based solely on their mid-IR colors without additional information such as variability or X-ray emission. We also define a "stellar" box which should contain most normal stars. The edges of this box are defined by −0.15mag < [3.6]–[4.5] < 0.15 mag and −0.15mag < [5.8]–[8.0] < 0.15 mag. The coldest stars, reaching into the brown dwarf regime, become red in [3.6]–[4.5] due to methane absorption in the 3.6 μm band, but have [5.8]–[8.0] colors near zero and thus stay to the left of the AGN wedge (Stern et al. 2007; Eisenhardt et al. 2010).

Beyond these mid-IR color diagnostics, we also use the extensive multi-wavelength and spectroscopic data for the field to characterize the sources. We know that the variable sources must be dominated by AGNs and a small fraction of stars, so the basic check on our selection criteria is that our selection limits on σ₁₂ and r preferentially select sources with AGN signatures, without introducing large numbers of contaminating sources.

We will now characterize the variable sources using three approaches. First, we simply examine the distribution of objects in mid-IR color. Second, we examine the variability of AGN candidates identified at other wavelengths. Finally, we consider the spectral properties of the objects with known redshifts.

Figure 5 and Table 3 characterize the distribution of objects in mid-IR color. The Photometry Group 3, Variability Level 2 candidates are concentrated in either the modified AGN wedge, with a distribution very similar to that of X-ray sources (Gorjian et al. 2008; see Figure 4), or in the stellar box. While only 0.5% of the Photometry Group 3 objects outside these two regions are Variability Level 2 sources, 7.4%, 4.6%, and 1.1% of the sources in the old AGN, modified AGN, and stellar regions are Variability Level 2 sources. The estimated false positive rates are 44%, 6%, 7%, and 58% for these four regions (outside, old AGN wedge, modified AGN wedge, stellar box) in the color–color space, respectively. This simply means that, for example, of the 7.4% variable sources in the old AGN wedge 6% are false positives, so in fact 7.4% × (1 − 0.06) = 7.0% of all objects are real variables and 7.4% × 0.06 = 0.4% of all objects are falsely identified as variable. Increasing the variability significance level reduces the variable fractions roughly in the same proportion, suggesting that even the Variability Level 2 candidates are dominated by true variables. The distribution also validates our assumption that the variability is dominated by AGNs—75% of the candidates lie in the modified AGN region for Variability Level 4. The selection remains pure even for Variability Level 2 because the requirement that the [3.6] and [4.5] light curves are tightly correlated (r>0.8) eliminates most of the contamination even when we only require σ₁₂>2 (see Figure 3). Most variable sources in the modified AGN wedge are known to be AGNs based on other evidence. Only 16%, 10%, and 7% of the Photometry Group 3, Variability Level 2, 3, and 4 sources lack photometric or spectroscopic signs of AGN activity. Most of these, however, are likely to be false positives given our estimated false positive rates of 7%, 4%, and 3%. For the Photometry Group 2 objects, where we lack accurate [5.8]−[8.0] color, we can only investigate the distribution in [3.6]−[4.5] color, and the side panel to Figure 5 shows the fraction of variable sources as a function of this color. Not surprisingly, the fraction is high (∼10%) when the color is red, corresponding to AGNs, and then drops dramatically for colors consistent with normal galaxies. The fraction of variables in the other two photometry groups are dramatically lower than for Photometry Group 3 for two reasons. First, a source must be intrinsically more variable in order to meet any variability criterion because the intrinsic measurement errors are larger. Second, with less well-measured colors, objects are scattered from their true location in the color–color space, and the differential scattering of galaxies into the AGN regions and AGN out of these regions dilutes the selection criterion.

Zoom In Zoom Out Reset image size

The 86 (Photometry Group 3, Variability Level 2) variable objects inside the stellar box are red, with optical colors of R − I ≳ 1 (i.e., V − I ≳ 2), so they must be dominated by red giants or main-sequence stars. The red and asymptotic giant branch stars include the OGLE Small Amplitude Red Giants (OSARGs; e.g., Soszyński et al. 2004) with I-band amplitudes of ∼0.1 mag and periods between 10 and 1000 days, the Long Period Variables (LPVs; e.g., Soszyński et al. 2009) with amplitudes of up to ∼3 mag and periods of 10–1000 days, and the Long Secondary Period variables (LSPs; e.g., Soszyński 2007) with periods of 30–1000 days. However, given their absolute magnitudes of M_I ≈ −4 mag, they would have to be at a distance of ∼200 kpc to match the observed R, I, and IRAC magnitudes. Since this scenario seems unlikely, the variable objects inside the stellar box could be the main-sequence flaring M dwarfs (e.g., Kowalski et al. 2009) with absolute magnitudes of M_I ≈ 7 mag, and distances of up to 1 kpc. Based on Kowalski et al. (2009), we expect to find ∼10 M dwarf flares in our data. The R − I and I − [3.6] colors of the variable objects inside the stellar box are inconsistent with the galaxy templates from Assef et al. (2010b). We have not pursued these sources in any detail because they are almost certainly quotidian variable stars and false positives. The estimated false positive rate of 58% is very high.

Next, we examine the variability of sources whose X-ray, MIPS 24 μm, or radio fluxes indicate that the source is an AGN. The color distributions of these sources are shown in Figure 6 and the variability statistics are summarized in Table 4. Here we focus on Photometry Group 1 since we are not dependent on having mid-IR colors to interpret the results. First, we see that a significant fraction of the X-ray sources show mid-IR variability, with 12% (=12.1% × (1 − 0.04) accounting for false positives) being variable even with no limits on the mid-IR photometric uncertainties. The higher fractions in Photometry Groups 2 and 3 are largely due to our improved ability to detect a given level of variability as we require smaller photometric uncertainties. The MIPS quasar candidates, which are found to be either AGNs (∼80%) or optically unresolved star-forming galaxies upon spectroscopic examination, also show high variability rates, with 11% (=11.7% × (1 − 0.05)) of the Photometry Group 1 objects showing variability at Variability Level 2. The variable fraction again increases as we go to the brighter photometric groups. The two sets of radio sources show significantly lower levels of variability. This is consistent with the results of Hickox et al. (2009) and Griffith & Stern (2010) who show that, as compared to X-ray and mid-IR selected AGNs, radio-selected AGNs have lower Eddington ratios and are much less likely to be AGN dominated. Thus, variability in radio-selected AGNs will be more challenging to discern relative to the photometric noise due to the host galaxy.

Zoom In Zoom Out Reset image size

Finally, we compare the variable sources with the template types used to estimate the source redshifts in the AGES survey, as shown in Figure 7 and summarized in Table 5. This is somewhat redundant with the previous two tests because the objects targeted as AGNs are chosen to be in the IRAC AGN wedge, to show X-ray or radio emission, or to be optical point sources with non-stellar 24 μm fluxes. The success rate for the spectroscopic identification of candidate AGNs in AGES is quite high. The sample of stars with spectra is very heterogeneous. Roughly half of the stars are F and G stars observed by AGES as spectrophotometric calibrators. Most of the remainder are K and M stars which were flagged as candidate AGNs either because of their X-ray fluxes or the presence of a mid-IR excess. We again see that the fraction of variable sources is highest for the AGNs, with roughly 11% meeting the Level 2 variability criterion.

Zoom In Zoom Out Reset image size

In summary, small relatively shallow mid-IR variability surveys like SDWFS are a workable but not a highly efficient means of identifying AGNs. Roughly speaking, the efficiencies are of order 10% for variability selection criterion providing relatively low false positive rates. It does, however, overcome the weaknesses of mid-IR color selection to identify AGNs outside the Stern et al. (2005) selection region. With difference imaging, the host galaxy flux that moves the source outside this color region ceases to be a problem other than through its effects on the photometric noise.

The structure function is a measure of the expected changes in an object's flux as a function of the time between the measurements. The time lag τ is defined as the time difference between any two data points in a light curve, where the observed time lag is τ = |t_i − t_j| and the rest-frame time lag is τ = |t_i − t_j|/(1 + z). For each time lag, one wants to find all pairs of points in a light curve and calculate the variance of the magnitude differences, Δm = m_i − m_j. Such a variance can also be calculated by averaging over many objects of the same class (i.e., AGNs) for each time lag between any two epochs, as we do in this paper.

Given the depth of the AGES redshift survey, all z>1 objects are spectroscopically confirmed AGNs (see Figure 2 in Stern et al. 2005). We measured the mid-IR structure function of z>1 Photometry Group 3 AGNs to fixed magnitude limits of [3.6] < 16, 17, and 18 mag. Only ∼10% of these objects were actually detected as significantly variable sources (see Table 5), so most of the observed variance in the magnitudes is at the noise level. This is a rather different regime from most optical studies, where the variance is dominated by variability, and makes it critically important to correctly characterize the contribution to the variance from the measurement errors. We do this by matching each AGN with four unresolved galaxies of nearly identical observed flux. We chose these galaxies to be Photometry Group 2 and 3 objects with colors [3.6]–[4.5] < 0.1 mag and 0.5mag < [5.8]–[8.0] < 1.5 mag, thus avoiding AGNs. We then compute the structure function as

$$\begin{eqnarray} S(\tau) &=& \bigg[ { 1 \over N_{\rm qso}(\tau) } \sum _{i<j} (m(t_i) - m(t_j))^2_{\rm qso} \nonumber \\ && - { 1 \over N_{\rm gal}(\tau) } \sum _{i<j} (m(t_i) - m(t_j))^2_{\rm gal}\bigg]^{1/2}. \end{eqnarray} \tag{ 5 }$$

In addition to the simple calculation over the 160, 690, and 1100 AGNs with [3.6] <16, 17, and 18 mag, we also calculated the structure functions for 1000 bootstrap resamplings of both the AGN and galaxy lists. These give consistent results, and we generally report the median result of the bootstrap samples and the region encompassing 68.3% of the samples about the median. The small values at short lags are very sensitive to any systematic problems in the noise estimate, and while this empirical approach to correctly characterizing the noise appears to work well, our estimates of the amplitude of the structure function on the longer timescales are probably more robust than our estimate of the slope of the structure function.

Figure 8 shows the resulting observed and rest-frame structure functions for the [3.6] and [4.5] bands as a function of magnitude limit. We generally see more variability on longer timescales. We fit our structure functions with a power law

$$\begin{equation} S(\tau) = S_0 \left({ \tau \over \tau _0 } \right)^\gamma \end{equation} \tag{ 6 }$$

using τ₀ = 4 yr (2 yr) for the observed (rest frame) estimates. Table 6 presents the resulting fits. The [3.6] and [4.5] structure functions are consistent with each other, but also significantly steeper than the i-band structure function from an ensemble of SDSS quasars (Vanden Berk et al. 2004), with γ = 0.303 ± 0.035, and higher amplitude, with S₀ = 0.20 mag for τ₀ = 2 yr.

Zoom In Zoom Out Reset image size

The steeper mid-IR slope could be explained by the emission being dominated by larger physical scales than the optical emission. At these redshifts, the observed frame mid-IR corresponds to the rest-frame near-IR, where the emission has contributions from both the disk and hot dust. If we focus on the more compact disk emission, the characteristic scale of an accretion disk for emission at wavelength λ, defined by the point where the disk temperature matches the photon energy, is

$$\begin{equation} R_\lambda \simeq 10^{16} \left({ \lambda _{\rm rest} \over \mu \hbox{m} }\right)^{4/3} \left({ M_{\rm BH} \over 10^9 \,M_\odot } \right)^{2/3} \left({ L \over \eta L_E }\right)^{1/3} \rm cm, \end{equation} \tag{ 7 }$$

where L/L_E is the Eddington ratio and $L= \eta \dot{M} c^2$ defines the radiative efficiency η. Taking L/L_E = 1/3, η = 0.1, and z = 1, the optical emission for M_BH = 10⁹ M_☉ comes from a region only a light day (0.003 light yr) in size, while the mid-IR scale is approximately 10 light days (0.03 light yr). Thus, our shortest lags correspond to the fastest possible response time given the size of the emission region since the typical AGES quasar has M_BH ≃ 10^8.5 M_☉ (Kollmeier et al. 2006). The orbital, thermal, and viscous timescales are even longer, and Kelly et al. (2009; also see Kozłowski et al. 2010) find that the optical variability can be characterized by an amplitude and a timescale that might be associated with thermal timescales. In any case, the mid-IR emission is coming from scales so large that it is hard to respond to any driving process on the timescale of our shortest observed lags. This should suppress variability on these short timescales compared to the optical emission, which can respond more easily, leading to a steeper mid-IR structure function.

Next we subdivided the sample into two redshift bins, 1 < z < 2 and 2 < z < 3, and then divided the objects in these bins at the median value of their absolute magnitudes. Figure 9 shows the structure functions of these subdivided samples as compared to the best fits to the full samples. At least for the well-determined values at longer lags, there is a pattern that shorter rest-frame wavelengths vary more than longer rest-frame wavelengths. Also for 1 < z < 2 AGNs, the lower luminosity systems vary more than the higher luminosity systems. This is consistent with earlier findings in the optical part of the spectrum (e.g., Hook et al. 1994; Hawkins 2002; Vanden Berk et al. 2004; de Vries et al. 2005). In general, the trend of more variability with longer lags also holds, although there are some exceptions in the [4.5] results.

Zoom In Zoom Out Reset image size

We also subdivided the data into eight smaller redshift– absolute magnitude bins and performed a global fit that includes a dependence on wavelength in the following form (see also Ivezic et al. 2004):

$$\begin{equation} S(\tau) = A\left[1+B\left(\langle M_J \rangle +25\right)\right]\left(\frac{\tau _{\rm rf}}{\rm 2\; yr}\right)^C\left(\frac{\langle \lambda \rangle }{1\; \mbox{$\mu $m}}\right)^{-D}, \end{equation} \tag{ 8 }$$

where A, B, C, and D are parameters of the fit, 〈M_J〉 is the average absolute J-band magnitude in a bin (the rest-frame M_J magnitudes are taken from Assef et al. 2010a), τ_rf is the rest-frame time lag, and 〈λ〉 is the average rest-frame wavelength in each bin. We either fit (1) all the data or (2) exclude the short time lags (τ < 0.1 yr) as the structure function contains little or no signal in this regime. For case (1) we found A(1) = 0.13 ± 0.01, B(1) = 0.19 ± 0.03, C(1) = 0.41 ± 0.05, D(1) = −0.04 ± 0.12, while for case (2) we found A(2) = 0.13 ± 0.01, B(2) = 0.19 ± 0.03, C(2) = 0.46 ± 0.06, D(2) = −0.02 ± 0.13. We find no detectable dependence of the structure function on wavelength, even though the data cover the redshift range of z = 1–3, the absolute magnitude range of M_J = −24.5 to −26.5 mag, and the rest-frame wavelengths λ_rf = 1–2 μm. As compared to optical studies (e.g., MacLeod et al. 2010), the mid-IR structure function appears to have stronger dependence on the luminosity (parameter B), and no dependence on wavelength (parameter D). Note that the large uncertainty in D (σ_D = 0.13) can potentially hide, within the 2σ–3σ, the wavelength dependence of D ≃ 0.3 observed in optical studies. Another possibility is that the mid-IR wavelength dependence is much weaker than that observed in optical.

Table 7 shows the variability amplitudes at the observed time lag of 4 yr for various classes of objects from Figure 4 assuming a fixed slope of τ^1/2. Objects with z>1 are AGNs, and show uniform and consistent variability amplitudes of ∼0.09 mag as a function of magnitude in both channels for fixed γ = 0.5. The X-ray selected sources inside the modified AGN wedge are expected to be AGNs, and not surprisingly, they also show similar variability amplitudes of ∼0.10 mag to the z>1 objects. The mid-IR counterparts of the 24 μm (MIPS)-selected sources inside the modified AGN wedge, again, have almost identical variability amplitudes. The radio-selected objects, however, have a drastically different distribution in their mid-IR colors (top right panel in Figure 4) from that of AGNs, suggesting that a large fraction of them may not be AGNs. The structure functions based on the radio-selected objects inside the modified AGN wedge may be driven by non-AGN objects, and we find somewhat smaller variability amplitudes of S₀ = 0.05 ± 0.01 mag and S₀ = 0.06 ± 0.02 mag in [3.6] and [4.5], respectively.

Given the structure function, we can now estimate the expected statistical properties of the AGNs in the current survey or any extension by adding new epochs. The structure function S(τ) is related to the correlation function by C(τ) = V²_∞ − S(τ), where V_∞ is the variance as τ → ∞. This makes it straightforward to compute the mean χ² of an AGN 〈χ²〉 = 〈(m_i − 〈m〉)²/σ² since 〈m_im_j〉 = σ²δ_ij + C(|t_i − t_j|), where σ is the measurement uncertainty. Because we fit for the mean, 〈m〉, the exact value of V_∞ is irrelevant. For our observed frame power-law structure function, the result for the SDWFS data and cadence is 〈χ²〉 = 3 + 1.6S²₀/σ², where S = S₀(τ/τ₀)^0.52, S₀ ≃ 0.1 mag, and σ is the typical photometric uncertainty. The constant 3 is simply the number of degrees of freedom after fitting the mean, so a non-variable source would have 〈χ²〉 = 3. This expression makes it clear why our AGN completeness is so low and drops rapidly with magnitude (increasing σ, see Figure 2).

We can, in fact, simulate the completeness. We can decompose the correlation matrix C_ij = C(|t_i − t_j|) into eigenvectors v_i and eigenvalues λ²_i. If G(x) is a Gaussian deviate of dispersion x, then random realizations of light curves consistent with the structure function are

$$\begin{equation} {\bf m } = \sum _i {\bf v}_i G(\lambda _i) + G(\sigma) \end{equation} \tag{ 9 }$$

(see Rybicki & Press 1992). If we assume the variability in both bands is identical but the noise is independent, then we can compute the fraction of such AGNs we would select as a function of the noise σ using Monte Carlo simulations. Figure 10 shows the approximate completeness as a function of magnitude for sources with r>0.8 and either σ₁₂>2 or 4. We use the observed structure function and the normalized measurement uncertainties to track the median values for [3.6] in Figure 2. For simplicity, the [4.5] band is treated as being identical to the [3.6] band. As might be expected from the structure of Figure 2, the completeness is high when the measurement errors are small compared to the amplitude of the structure function and then drops rapidly for fainter sources.

Zoom In Zoom Out Reset image size

In general, the statistics of these simulated light curves agree with observations, indicating that our estimates of the structure function are correct. However, they also show that we are "pushing the envelope" of variability, because for most AGNs the variability amplitudes are merging into the measurement uncertainties, as can be seen in Figure 2.

In this paper, we performed DIA photometry on the four epochs of SDWFS data. The common area of the Spitzer mosaics covers 8.1 deg² of the NDWFS Boötes field, and contains 474,179 mostly extragalactic objects. We concentrate on the deeper 3.6 and 4.5 μm channels, and provide variability catalogs and light curves. These variability data are cross-correlated with the AGES redshift survey, X-ray and radio data, and 24 μm AGN candidates from the same area of the sky to study the variability of various classes of objects.

We study the mid-IR variability of objects based on the significance level of their joint variability (σ₁₂) in both the [3.6] and [4.5] channels focusing on the sources that also show a high correlation (r>0.8) between the two channels. Our standard selection criteria with σ₁₂>2 identify ∼5100 objects as variables, which constitute 1.1% of all detected sources in the field. We find that the majority (76%) of the mid-IR variable sources in the NDWFS Boötes field are AGNs, where we define AGNs to be objects that satisfy any of the following criteria: (1) lie inside the modified AGN wedge, (2) are X-ray, 24 μm, or radio AGN candidates, (3) have z>1, or (4) have an AGN spectroscopic template in any region.

Amongst all the extragalactic classes of objects, AGNs are the primary sources able to produce significant changes in their mid-IR luminosity over a several year timescale. Thus, mid-IR variability may be used as a new tool for selecting AGNs, especially if combined with other methods. However, in our relatively shallow few-epoch survey, only ∼15% of AGNs are sufficiently variable to be detected; therefore, variability-selected samples of AGNs are highly incomplete. This incompleteness is caused by the typical variability amplitude being comparable to the survey photometric errors, and is not due to the small number of available epochs. For our standard criteria, 14% of X-ray sources, 17% of 24 μm-selected AGN candidates, and 15% of spectroscopically confirmed AGNs are detected by their variability.

The variable AGNs primarily occupy the mid-IR AGN selection wedge of Stern et al. (2005), with an extension to bluer [3.6]–[4.5] colors where the host galaxy dominates the mid-IR colors of low-luminosity AGNs (see Gorjian et al. 2008; Assef et al. 2010b). If we redefine the original Stern et al. (2005) selection wedge to encompass X-ray or variability-selected AGNs, we find that variable objects bluer than [3.6]–[4.5]>0.3 mag are likely to be AGNs. This opens a new window for the Spitzer Warm Mission, where the [3.6]–[4.5]>0.3 mag color cut combined with the variability criterion will be a solid indication of AGN activity.

The structure function of z>1 AGNs is well described by a power law with a logarithmic slope of γ = 0.56 ± 0.18 and amplitude S₀ = 0.11 ± 0.02 mag in [3.6] and γ = 0.44 ± 0.14, S₀ = 0.12 ± 0.01 mag in [4.5] for a rest-frame time lag of τ₀ = 2 yr. These structure functions are steeper than those observed in optical bands, and we argue that this is due to the inability of either disk or dust emission to respond on the short timescales seen in optical studies. The mid-IR structure functions go approximately as τ^1/2, which corresponds to the short timescale (τ ≪ τ_model) behavior in the damped random walk model used by Kelly et al. (2009), Kozłowski et al. (2010), and MacLeod et al. (2010). If we examine the variability of various subsamples, we detect more variability in low-luminosity objects observed at shorter rest-frame wavelengths. The mid-IR sources corresponding to X-ray and 24 μm AGN have structure functions identical to that of z>1 AGNs. However, this is not the case for the counterparts of radio sources which have S₀ ≃ 0.05 ±  0.01 mag in both [3.6] and [4.5] bands, as compared to S₀ ≃ 0.10 ±  0.01 for either z>1 objects, X-ray, or MIPS sources. One explanation may be that the radio sample has more contamination by non-AGN (e.g., due to confusion between radio cores and lobes), another is that the radio-selected AGNs are accreting in a different mode (lower accretion rates or efficiencies; e.g., Tasse et al. 2008) and that variability in such a mode is inherently lower.

Expanding the SDWFS survey with new epochs during the Spitzer Warm Mission would represent a significant improvement. Given our structure functions, we can estimate the results from expanding the SDWFS survey with new epochs. We considered two cases of (1) doing one new epoch in 2011 March and (2) doing new epochs in both 2011 and 2012 March. The improvements in the mean χ² are significant, with 〈χ²〉 = 4 + 3.0S²₀/σ² and 5 + 4.3S²₀/σ², respectively, for AGNs, compared to 〈χ²〉 = 4 or 5 for non-variable objects. These represent a significant improvement over the SDWFS baseline, where 〈χ²〉 = 3 + 1.6S²₀/σ² for AGN and 3 for non-variable objects. Figure 10 shows the increases in completeness for these five and six epoch scenarios. Most of the gain for detecting variable sources comes from adding one additional epoch. The greatest astrophysical gain from two additional epochs comes from being able to better map out and subdivide the structure function by luminosity, redshift, and source type to better characterize the physics of AGN variability at these wavelengths. In particular, it would confirm the different slope of the mid-IR structure function from that observed in the optical.

This work is based on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of the Technology under contract with the National Aeronautics and Space Administration (NASA). Support for this work was provided by NASA through award numbers 1310744 (C.S.K. and S.K.), 1314516 (M.L.N.A.), and 1317692 (H.A.S.) issued by JPL/CalTech. C.S.K. and S.K. are also supported by NSF grant AST-0708082. The NDWFS and the research of A.D. and B.T.J. are supported by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under cooperative agreement with the National Science Foundation. Support for M.B. was provided by the W. M. Keck Foundation.