Cosmic Infrared Background Fluctuations and Zodiacal Light

We have performed a specific observational test to measure the effect that the zodiacal light can have on measurements of the spatial fluctuations of the near-IR background. Previous estimates of possible fluctuations caused by zodiacal light have often been extrapolated from observations of the thermal emission at longer wavelengths and low angular resolution, or from IRAC observations of high latitude fields where zodiacal light is faint and not strongly varying with time. The new observations analyzed here target the COSMOS field, at low ecliptic latitude where the zodiacal light intensity varies by factors of $\sim2$ over the range of solar elongations at which the field can be observed. We find that the white noise component of the spatial power spectrum of the background is correlated with the modeled zodiacal light intensity. Roughly half of the measured white noise is correlated with the zodiacal light, but a more detailed interpretation of the white noise is hampered by systematic uncertainties that are evident in the zodiacal light model. At large angular scales ($\gtrsim100"$) where excess power above the white noise is observed, we find no correlation of the power with the modeled intensity of the zodiacal light. This test clearly indicates that the large scale power in the infrared background is not being caused by the zodiacal light.


INTRODUCTION
The study of astronomical backgrounds at various wavelengths allows the examination of sources that are intrinsically diffuse, or individually too faint or too confused to be detected. Over time, improvements in instrumentation may resolve increasingly fainter sources, but very faint and intrinsically diffuse sources always remain in the realm of background studies.
Studies of the cosmic infrared background (CIB) have aimed at measuring the cumulative stellar emission of galaxies across the entire history of the universe. Measurements made by the DIRBE instrument on COBE provided the first space-based measurements of the absolute sky surface brightness at wavelengths from 1.25 to 240 µm with an angular resolution of 0. • 7 (Hauser et al. 1998). However, difficulties in accurately removing foreground contributions from the zodiacal light , and from Galactic stars and interstellar dust  prevented precise detections of the CIB except at the longest wavelengths.
Subsequent measurements and studies of the CIB at near-to mid-IR wavelengths have made use of telescopes and instruments that provide much higher angular reso-1 CRESST/University of Maryland -Baltimore County 2 Science Systems & Applications Inc. 3 NASA lution than DIRBE, but which usually cover only a small fraction of the sky, e.g. 2MASS, Spitzer, HST, IRTS, AKARI, CIBER. The higher angular resolution allows the exclusion or subtraction of stars and bright galaxies from sky brightness measurements, such that the contribution of Galactic stars is minimized, and the CIB that remains does not include the contributions of galaxies brighter than certain limits dependent on depth of the observations. These data sets are often better suited for measuring the spatial fluctuations or structure of the "source-subtracted" CIB, rather than for measuring the mean value of the CIB (Kashlinsky et al. 1996a(Kashlinsky et al. ,b, 2002(Kashlinsky et al. , 2005(Kashlinsky et al. , 2007Kashlinsky & Odenwald 2000;Odenwald et al. 2003;Matsumoto et al. 2005Matsumoto et al. , 2015Matsumoto et al. , 2011Thompson et al. 2007a,b;Cooray et al. 2007;Sullivan et al. 2007;Arendt et al. 2010;Matsuura et al. 2011;Pyo et al. 2012;Zemcov et al. 2014;Donnerstein 2015;Mitchell-Wynne et al. 2015;Seo et al. 2015). These studies seem to show reasonable consistency of fluctuation measurements reported by different experiments and different groups at λ 2 µm, though the picture is less clear at shorter wavelengths.
The origin of the fluctuations is not yet clear. The fluctuations may arise from any or all of: solar system and Galactic foregrounds, nearby extragalactic contributions, and distant extragalactic contributions. In most cases (e.g. Kashlinsky et al. 2005;Arendt et al. 2010;Matsumoto et al. 2011;Zemcov et al. 2014) foregrounds are estimated by extrapolation of measurements at other wavelengths and locations. This is particularly true of zodiacal light contributions. Mid-IR observations using ISO have limited zodiacal light fluctuations to < 0.2% on scales > 3 ′ (Abraham et al. 1997). More recent AKARI measurement set the limit even lower at 0.03% (Pyo et al. 2012).
Direct detection of zodiacal light influences in Spitzer IRAC measurements have been checked in existing deep data sets by constructing A − B difference maps, where A and B represent observations made at two different epochs, typically 6 months and/or 1 year apart (e.g. Kashlinsky et al. 2007). The expectation is that the first-order gradient (and any instrumental imprint) should reverse at 6-month intervals, while any smaller scale, physically distinct structures in the interplanetary dust cloud should not remain fixed in a given field (due to differential rotation of the cloud) and should not appear the same at different epochs separated by 6 months. Excess differences have not been seen in existing deep CIB studies. However, these are generally high latitude fields where the zodiacal light is faint and not strongly modulated. Additionally, at intervals of 6 months and (especially) 1 year, the interplanetary dust cloud may be sufficiently symmetric such that certain structures (e.g. those associated with the earth-resonant ring) may still cancel out in A − B difference maps.
As a more certain test for zodiacal light influences, our new Spitzer IRAC observations have monitored a low latitude field over an entire visibility window to obtain a data set where the zodiacal light is both brightest and most strongly modulated. The observations were planned to be sufficiently deep to detect the reported large scale background structure. Thus these data are uniquely suitable for checking if the zodiacal light intensity has any effect on the reported background fluctuations at large angular scales. This paper reports on this test and the results. The observations and data reduction are described in sections 2 and 3. Section 4 provides the characterization of the power spectra of the background. Section 5 discusses the correlations between the various components of the power spectrum and the zodiacal light intensity. Section 6 summarizes the results. An appendix provides additional detail on temporal variations in the data that are tracked and corrected by the self-calibration procedure that we apply. The appendix also features additional details on the effects of the source model, and on the comparison with previous CIB measurements.
2. OBSERVATIONS Most commonly observed extragalactic fields are chosen to be at high Galactic and ecliptic latitudes to minimize the influence of foregrounds. The observations presented here were designed, proposed, and approved specifically for the purpose of examining the effect of zodiacal light on CIB measurements. The COSMOS field (Scoville et al. 2007;Ashby et al. 2013Ashby et al. , 2015 is suitable for our experiment, as it lies at relatively low ecliptic latitude. We have selected a subregion at ecliptic coordinates (λ, β) = (151.73, −8.63) which is relatively free of bright sources. The patch is observed 5 times across an entire visibility window of the Spitzer spacecraft (Werner et al. 2004;Gehrz et al. 2007), covering the widest possible range of solar elongation, and thus brightness. The size, ∼ 10 ′ × 10 ′ , and depth, ∼ 4 hr per epoch, of the patch were chosen to be sufficiently wide and deep to distinguish the large scale fluctuations above the random white noise in the observations. We use Spitzer's IRAC instrument (Fazio et al. 2004) to collect 3.6 and 4.5 µm data. As IRAC's field of view is 5 ′ × 5 ′ , the observations require mosaicking 2 × 2 fields of view and a total observing time of ∼ 16 hr at each epoch. Table 1 lists the dates, solar elongations, and astronomical observation request (AOR) numbers of the observations. 3. DATA REDUCTION 3.1.

Self-Calibration
We reduce the data using the same self-calibration techniques (Fixsen et al. 2000) that have been previously employed (Arendt et al. 2010). The data reduction began with the IRAC cBCD individual frames and applied a data model of where D i is the measured intensity at pixel i in a single frame (q), S α is the true sky intensity at location α, F p is an offset for pixel i which is constant over all frames, and F q is an offset for frame q which is constant over all pixels (but variable with time). Each of the 5 epochs is self-calibrated separately to provide detector offsets, F p , that are appropriate for each epoch. Sky maps are generated on a pixel scale of 0.6 ′′ (half the size of the detector pixels) using an interlacing algorithm. In addition to regular sky maps, we also create A − B sky maps where all the odd numbered frames are multiplied by −1 before mosaicking the image. This has the effect of removing the contributions of fixed celestial sources and leaving only instrumental noise (and photon shot noise, see Section 4.1). Figure 1 shows the 3.6 and 4.5 µm images (S α determined by self-calibration) of the combined 5 epochs Figure 1. Combined five-epoch 3.6 µm (left) and 4.5 µm (right) mosaics of the study fields, illustrating the source density and the absence of bright stars and large galaxies. Only the two brightest stars in the 3.6 µm image show a small residual of incompletely corrected "column pull-down" artifacts (dark vertical stripes). The images are logarithmically scaled on the range [2.5 × 10 −4 , 1.0] MJy sr −1 after addition of an offset of 0.004 MJy sr −1 to avoid logarithmic scaling of negative data values. The images are 900×900 pixels or 9 ′ × 9 ′ in size (1 pixel = 0.6 ′′ ). Celestial north is at a position angle of 66. • 4 (counterclockwise from vertical). There is only ∼ 25% overlap between the fields. of observations. The images are cropped to show only the roughly uniformly covered region that was used for power spectrum analysis. Linear background gradients have been fitted and subtracted. Figure 2 displays the derived values of the detector offsets, F p , for 3.6 and 4.5 µm at all 5 epochs. At all epochs there are different patterns of light and dark latent images, which are very slowly decaying imprints of very bright stars as they were dithered across the detector in previous observations. In several cases the tracks of bright stars that slewed across the detector between pointings can also be seen. Residual stray light in the cBCD frames is also revealed and removed by the selfcalibration. Diffuse patches of stray light are created by the zodiacal light (and the sum of all other backgrounds) in the upper left and upper right corners of both the 3.6 and 4.5 µm detectors. The BCD pipeline uses estimates of the expected brightness of the zodiacal light to model and remove this stray light component. These self-calibration results show that at 3.6 µm that process works well initially, when the zodiacal light is bright, but has an increasing error at later epochs as the elongation increases and the zodiacal light becomes fainter. The opposite occurs at 4.5 µm, where the standard correction works best at the later high-elongation epochs and less well at the early epochs when the zodiacal light is brighter. All these features visible in the F p maps are systematic effects that are removed by the self-calibration of the data. Figure 3 shows expected general trend of the zodiacal light intensity as a function of time, as estimated by the Spitzer foreground model 4 and given as the ZODY EST keyword in the header of each frame. Note the 4.5 µm intensity is a stronger function of time or elongation than the 3.6 µm intensity, indicating that the color of the zodiacal light is not expected to be constant. The self-calibration model applied in Equation (1) assumes a fixed sky intensity. Thus any variations in brightness due to changing zodiacal light intensity across the span of the data being self-calibrated, or a single epoch, are absorbed by the variable offset term, F q . Figure 4 shows that F q varies by a much larger amount than the expected zodiacal light trends at 3.6 µm, but is similar to the expected zodiacal light trends at 4.5 µm. In the Appendix we show that the differences are real and represent corrections for transient instrumental effects (not fully corrected in the BCD pipeline), and residual linear gradients across the field.

Source Subtraction and Masking
For measurement of the power spectra of the background, resolved sources need to be masked or modeled and subtracted from the images. As in prior studies, we subtract the flux from sources above the noise level using an iterative algorithm described by Arendt et al. (2010). The iterations are halted when the skewness of the intensity distribution of the remaining pixels is zero (independently for each epoch). Because the removal is imperfect, the image is also masked using a mask derived from a surface brightness threshold in the original images for all epochs combined. The masking threshold can equivalently be expressed as a specified surface brightness, a specified maximum outlier in the The fractional decrease at 4.5 µm is larger than that at 3.6 µm. Plotting simply as a function of frame number (bottom) provides a slightly clearer look at the very small oscillations in intensity that are caused by moving up and down the zodiacal light gradient at different pointings.
distribution of surface brightness of unmasked pixels, or a specified fraction of area masked. In this case we chose the last constraint, limiting the masked out regions at both wavelengths to 25% of the image, leaving 75% of the image remaining. At 3.6 and 4.5 µm, this limit corresponds surface brightness thresholds of max(I ν ) − mean(I ν ) = 0.0057 and 0.0044 MJy sr −1 , and to maximum outliers of max(I ν )/σ Iν = 3.2 and 2.6 respectively.

POWER SPECTRA
The power spectra of the source-subtracted images are calculated as described by Arendt et al. (2010). The power on the horizontal (u = 0) and vertical (v = 0) axes in the Fourier domain is omitted when averaging the power in bins at different angular scales, 2π/q. This makes the results less susceptible to certain systematic errors (such as the residual column pull-down seen at 3.6 Comparison between the self-calibration F q offsets (black dots) and the ZODY EST keyword values (red dots) after subtraction of the mean values for each epoch. The lack of perfect correlation between these reflects errors in the zodiacal light model and the presence of additional temporally variable effects in the data (especially the firstframe effect at 3.6 µm, see the Appendix), which are removed by F q . µm in Figure 1), but limits the maximum angular scale to 382 ′′ = 540 ′′ / √ 2 instead of the full 540 ′′ size of the field. The uncertainties assigned to each binned point in the power spectra are the standard deviations of all measurements contributing to each bin. The resulting power spectra, for the five epochs combined, and for each epoch individually are plotted in Figures 5 and 6.
The power spectra are characterized as in Arendt et al. (2010) by fitting a combination of 3 simple components: where the parameters a 0 and a 1 are the amplitude and index of a power-law component that is modulated by the instrument beam or point response func-tion, P PRF (q), the parameter a 2 is the amplitude of the sky shot noise: a white shot-noise component that is also modulated by the beam (e.g. Poisson variation in the number of faint unresolvable sources in the beam at each location), and a 3 is the amplitude of a white (shot noise) component that is not modulated by the beam. As an alternate characterization, we also fit: (3) where the white noise component is replaced by the measured A − B power spectrum with no rescaling allowed. These characterizations of the power spectra are overplotted in Figures 5 and 6 and are tabulated in Table 2.
The values of reduced chi squared (χ 2 ν ) in Table 2 are calculated for each fit using the full data set, i.e. ν = 445 or 446 degrees of freedom for Equations (2) or (3). However, for the 3.6 µm power spectra, we strongly deweight the two data points at 2π/q > 200 ′′ because the estimated uncertainties are clearly inconsistent with the behavior of these data points. The absolute values of χ 2 ν are poor, indicating that the uncertainties may be underestimated or these models do not accurately reflect the power spectra. However, the overall aim of these fits is not to validate a particular physical model of the power spectra, but to reduce the power spectra to a small number of parameters that can be easily compared between power spectra (as we do below), and which provide a convenient approximation to the power spectra.
We note that one could apply a more physical model for characterization of the large-scale component of the power spectra, e.g. a ΛCDM template. However, a simple power law provides a sufficient approximation for a wide possibility of origins, given the angular scales ( 400 ′′ ) and quality of the data analyzed here.   Note-Units for a0, a2, a3, b0, and b2 are nW 2 m −4 sr −1 .

White Noise Component
The white noise component of the power spectrum, a 3 in Equation (2), includes instrumental noise, but it also includes the photon shot noise from celestial sources. In particular, the zodiacal light is the dominant brightness component. Because it is an approximately uniform source of emission, the power spectrum of its photon shot noise is not modulated by the beam. As a noise term it also does not cancel out in the construction of the A − B difference images, and therefore the power spectra of those images also include the photon shot noise of the zodiacal light.
The slight rise in the white noise component at the smallest angular scales is an artifact of mapping the data, sampled on 1. ′′ 2 detector pixels, onto a parallel sky map with 0. ′′ 6 pixels. Given our interlacing mapping algorithm, a slight mismatch in the mean level of any single frame will insert power into the map at the Nyquist frequency of the 0. ′′ 6 pixels of the sky map, i.e. at 1. ′′ 2. The multiplication of the image by a mask then convolves this power in the Fourier domain, spreading it to larger angular scales. Thus the shape of the turn up at the smallest spatial scales is related to the masking, and the amplitude is related to the size of the frameto-frame background errors. Mapping on a grid that is not so well aligned with detector orientation can mitigate the effect somewhat. For example, if the sky map is generated on grid that is rotated by 45 • to the de-tector orientation, then the white noise component does appear flat. Because this rise cannot be well fit with flat white noise, more than half the contribution to χ 2 ν comes from 2π/q < 2 ′′ . Figure 7 shows the very strong correlation between the zodiacal light intensity and the measured level of the white noise power spectrum. Extrapolation to zero intensity of the zodiacal light indicates that white noise power levels in the absence of zodiacal light would be 5.7 × 10 −12 and 5.2 × 10 −12 nW 2 m −4 sr −1 at 3.6 and 4.5 µm respectively. Comparison to the measured white noise power (a 3 in Table 2) indicates that photon shot noise of the zodiacal light contributes ∼ 40 − 60% of the amplitude of this component, depending on the epoch of the observations.

Sky Shot Noise and Power Law Components
The sky shot noise component, characterized by a 2 or b 2 , is not an essential component for fitting the observed power spectra for the five individual epochs. These observations are so shallow that the flat white noise component can dominate to sufficiently large scales where the power law component takes over. At 4.5 µm, the best fits are obtained without a sky shot noise component, though this requires shallower power law indices, a 1 and b 1 , than previously found (Arendt et al. 2010). This is likely caused by an inability to distinguish separate sky shot noise and power law components in these shallow data. Constraining the power law index to be Figure 7. Correlation between the white noise level (a3) and the zodiacal light intensity across the 5 epochs shows that roughly half of the white noise is correlated with the zodiacal light. The formal uncertainties in the white noise values are smaller than the plotted symbols (see Table 2) The Kendall's τ rank correlations (bounded on the interval [-1,1]) for these quantities at 3.6 and 4.5 µm are 0.95 and 1.0. The probabilities of finding these values of |τ | (or larger) given uncorrelated data are 0.02 and 0.01. a 1 = b 1 = 1.0 does result in a weak sky shot noise component that constitutes ∼ 15% of the power at 100 ′′ , but the constraint produces a poorer fit at scales of ∼ 20 ′′ . Results of this constrained fit are listed in the last section of Table 2. At 3.6 µm there is no similar motivation for a constrained fit, as the best fits find non-zero sky shot noise components (a 2 and b 2 ) and power law slopes (a 1 and b 1 ) that are similar to those previously derived (Arendt et al. 2010). Conversely, omitting the sky shot noise component (a 2 and b 2 ) at 3.6 µm causes the power law component to flatten, to better fit the spectrum at 2π/q 10 ′′ . For the all-epochs-combined power spectrum, the power law flattens completely, (a 1 ∼ 0 and b 1 ∼ 0), to become the sky shot noise component, leaving the rising power at 2π/q 10 ′′ poorly fit.
The amplitudes, a 0 and b 0 , of the power law fitting the large scale power at 4.5 µm are much lower than at 3.6 µm, unlike prior results where powers were compara- ble (Arendt et al. 2010). Figure 8 shows no correlation between the amplitude of the power law component and the intensity of the zodiacal light.

DISCUSSION
The correlation between the zodiacal light intensity and the white noise in the power spectra (Figure 7) indicates both a constant component (i.e. the intercept as the zodiacal light intensity goes to 0) and a component that is proportional to the zodiacal light intensity. Instrumental noise terms (e.g. read noise and dark current), and the photon shot noise from the mean intensity of unresolved Galactic and extragalactic backgrounds would contribute to this constant term. In Figure 9 we present a comparison between the measured white noise power as a function of zodiacal light intensity, and the expected noise levels.
Using the equations and parameters presented in Sec- Figure 9. Comparison of white noise predictions and measurements as a function of background intensity. The expected white noise levels calculated from the IRAC Instrument Handbook using "High", "Medium", and "Low" background levels, and with the background set to 0.0, are shown by the black circles. The mean white noise levels derived from the measured standard deviation of individual frames at each of the 5 epochs are shown by red diamonds. The mean white noise levels derived from the measured standard deviation of the self-calibrated sky maps are shown by orange + symbols. The measured white noise levels obtained by fits to the power spectra are indicated by the blue × symbols. See text for explanation of consistencies and discrepancies.
tion 2.5 of the IRAC Instrument Handbook 5 , we can evaluate the expected 1-σ noise level for extended emission for 100s frame time exposures during the IRAC warm mission. These numbers are listed in Table 3 (Note 1), and, after conversion to white noise power, P = σ 2 Ω pix (Table 3, Note 5), are shown in Figure 9. For direct comparison to these expected noise estimates, we calculated the mean of the standard deviations of all (576) individual exposures at each epoch. A robust procedure was used to exclude the effect of stars and other statistical outliers. These numbers (Table 3, note 2 and note 6) are in generally good agreement with the expected uncertainties, with the exception that at 3.6 µm the measured standard deviations show more variation with zodiacal light intensity than expected. This discrepancy may indicate that the zodiacal light model predicts the correct mean intensity of the zodiacal light, but underestimates the modulation of the intensity as a function of time. This effect is very similar to that observed at the north ecliptic pole (NEP) by Krick et al. (2012). The standard deviation in mosaicked sky maps should be reduced by a factor of 1/ √ N , where N is the number of frames contributing at each pixel (Table 3, note 3). In this case there are 144 frames per pointing, but interlaced mapping on 0. ′′ 6 pixels reduces the per pixel coverage by a factor of 4, so N = 36. The actual measured standard deviations of the self-calibrated mosaics ( Table 3, note 4) are ∼ 10% smaller at 3.6 µm and ∼ 24% smaller at 4.5 µm. Converted to a power (Table 3, note 7), these values are plotted as the orange + symbols in Figure 9. They show a similar trend to the white noise levels derived from the power spectra (a 3 , blue × symbols), but are biased upwards because they average the power at all angular scales, rather than just fitting the minimal white noise at the smallest angular scales. The roughly constant size of the bias is a further indication that the large scale structure is independent to the zodiacal light intensity. Ideally, the white noise level extrapolated to a zodiacal light intensity of 0, could be interpreted in terms of instrumental noise plus the photon shot noise from galactic and extragalactic backgrounds. However, systematic errors in the zodiacal light model, as clearly evident at 3.6 µm, will directly affect the intercept. This is the likely reason for the mismatch between the extrapolated intercepts for the measured 3.6 µm power (either from the fitted power spectra or from single frames) and the expected noise power as calculated from the IRAC Instrument Handbook. However, the ultimate origin of the white noise component is not very important for current CIB studies, which normally subtract any "instrumental" or A − B noise term that is constructed to be independent of fixed sources on the sky.
The large-scale component in the source-subtracted CIB power spectrum is the term that is of greatest cosmological interest. Prior studies agree that there is power here in excess of that expected from the faint unresolved galaxies extrapolated from known galaxy populations (Kashlinsky et al. 2005;Sullivan et al. 2007;Helgason et al. 2012;Cooray et al. 2012). The lack of significant correlation between the large-scale power and the zodiacal light intensity (Figure 8) suggests that the zodiacal light is not influencing the large scale power. Additionally, we note that while the zodiacal light intensity is greater at 4.5 µm than at 3.6 µm, the data exhibit weaker large scale power at the 4.5 µm than at 3.6 µm. This also indicates that the zodiacal light is not the main source of the large scale power. Note-(1) Expected extended source surface brightness sensitivity for single 100s frames, as calculated from equations in Section 2.5 of the IRAC Instrument Handbook.
(3) Measured average standard deviation (outliers excluded) of single frames scaled by of 1/ √ N , (4) Average standard deviation (outliers excluded) of self-calibrated sky map (0. ′′ 6 pixels). (5) White noise power derived from expected sensitivity (as in Note 1). (6) White noise power derived from standard deviation of single frames. (7) White noise power derived from standard deviation of sky map. (8) White noise power derived from fit to power spectrum.

SUMMARY
We have performed an experiment specifically designed to measure the impact of the zodiacal light on the estimate of the spatial fluctuations of the CIB. To provide the greatest possible sensitivity to zodiacal light effect, our test monitored a fixed patch in the COSMOS field at low ecliptic latitude as the mean zodiacal light intensity varied over the full accessible range of brightness (or solar elongation). The CIB spatial power spectrum was calculated at 5 epochs over this 5-week interval. The power spectra are characterized as the sum of (a) a white noise component, (b) a sky shot noise component, and (c) a power law component dominating on large angular scales.
We find that approximately half of the white noise component is correlated with the varying mean intensity of the zodiacal light. Photon shot noise of the zodiacal light is expected to be the main contribution to this correlated component. Detailed analysis of the nonzodiacal light portion of the white noise is limited by inaccuracies of the zodiacal light model in predicting the intensity in this direction as seen from Spitzer's location within the interplanetary dust cloud.
The sky shot noise in the angular power spectra of the background is not reliably distinguished in the relatively shallow observations of this experiment. The power law component does not show significant corre-lation with the mean zodiacal light intensity at 3.6 or 4.5 µm. This confirms that observed spatial fluctuations at large scales ( 100 ′′ ) are not being influenced by zodiacal light. Prior observations had been less conclusive because they were usually limited to high ecliptic latitudes, where the zodiacal light is faintest, and to epochs ∼ 6 or 12 months apart, where there should be minimal modulation of the zodiacal light intensity.
We thank the referee, M. Zemcov, for comments that improved this paper. Additional helpful comments were provided by M. Ashby. This work is based on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. This work was funded by JPL under Spitzer Cycle 8 funding contract 1464716, and in part through NASA/12-EUCLID11-0003 "LIBRAE: Looking at Infrared Background Radiation Anisotropies with Euclid". This research has made use of NASA's Astrophysics Data System Bibliographic Services, and the IDL Astronomy Library (Landsman 1993).
Facility: Spitzer (IRAC) Software: IDL elongations are sampled. The corresponding variations in brightness are thus mapped as temporal variations in F q . Additionally, the self-calibration is degenerate with respect to linear gradients across the field, which would have the same effect as the intrinsic zodiacal light gradient. We model spatial gradients in F q as a constant plus linear function of the x and y coordinates. The choice of the coordinate system is irrelevant, as the coefficients will simply adjust appropriately for any chosen system.

A.4. Exponential Decay
The final systematic effect evident in F q is an apparent decaying response (with a negative amplitude) across each AOR. This trend can be fit by a simple exponential decay as a function of frame number in each AOR. However, there are both fast and slow decay terms with e-folding constants of 9 and 70 frames. We model this decay as a linear combination of these two exponential decays. This effect was also noted by Krick et al. (2011), but it appears more cleanly here.
The net model for the temporal variation in F q is thus given by: E 9 e −f rame/9. + E 70 e −f rame/70. Figure A1 shows the residual F q after successive subtraction of each of the components of F q model (with arbitrary offsets for clarity).
At 3.6 µm the first-frame effect is responsible for most of the variance in F q , as seen by comparing the derived F q (black dots) with the derived F q minus the fitted first frame effect (red dots). At 4.5 µm, the zodiacal light trend is clearly dominant. Therefore at 4.5 µm the red dots represent the derived F q minus the fitted zodiacal light trend.
The comparison of the red and orange dots in both panel shows the subsequent subtraction of the zodiacal light trend (3.6 µm) and the first frame effect (4.5 µm). The zodiacal light trend has little effect at 3.6 µm. The influence of the first frame effect is now clear at 4.5 µm, though far smaller in amplitude than at 3.6 µm.
Comparison between the orange and green dots shows the subsequent subtraction of the spatial gradient terms at both wavelengths. The AORs were designed such that the spatial gradients map into oscillations with a period of 24 frames. This makes them distinguishable from the slower monotonic change in the zodiacal light intensity which occurs as a function of time. The amplitudes of the gradient terms are visibly larger in the earliest AOR and decrease across AORs as does the zodiacal light intensity.
The green dots clearly show the exponential decay behavior. At 3.6 µm, the slow 70-frame decay is evident in the first AORs, but there is a gradual transition to the faster 9-frame decay in the later AORs. At 4.5 µm, the amplitude of the decay is much smaller, and the slow 70-frame decay is dominant for all AORs.
After removal of the exponential decays, the blue dots show a fairly random distribution, with a standard deviation Figure A1. The successive removal of identifiable components of F q . For 3.6 µm (left) the black dots show the full F q (as in Figure 4). Each segment of 576 frames corresponds to one AOR. The colored dots show the residual F q after successive subtraction decreasingly smaller components: a first-frame effect (red), the zodiacal light variation (orange), an arbitrary linear gradient (green), and, an exponential recovery term (blue). At 4.9 µm (right), the relative strength of the components differs and sequence shown is after successive subtractions of: the zodiacal light variation (red), a first-frame effect (orange), an arbitrary linear gradient (green), and an exponential recovery term (blue).
that is > 10 times smaller than present in the original F q . These residual dispersions of σ = 4.7 × 10 −4 and 2.5 × 10 −4 MJy sr −1 represent an upper limit on the accuracy of the F q term of the self-calibration. The residual variation may still be accounting for real effects in the data, but the specific nature of such effects has not been identified here.
B. MODEL DEPTH A critical aspect of this analysis is the use of a source model to remove the effects of (a) the emission of extended sources and PSF wings that project beyond the masked areas, and (b) faint sources which cannot be masked without adversely decreasing the fraction of area available for analysis. Figure B2 shows power spectra of the masked images as the depth of the source model (i.e. the number of components subtracted) is linearly increased. Our choice is to use the model depth at which the skewness of the intensities of the unmasked pixel is zero. This assumes that a positive skewness signifies the presence of residual sources in the image, and that a negative skewness indicates that the model has begun subtracting the positive side of the noise distribution rather than actual sources. The zero-skewness models are evaluated independently for each epoch, and are indicated by red lines (with error bars) in Figure B2. The power spectra of the A − B images are indicated by the blue lines. This figure is analogous to Figures 8 -13 of Arendt et al. (2010).
C. COMPARISON TO PRIOR RESULTS Figure C3 compares the power spectra measured here with those measured previously in Arendt et al. (2010). The previous power spectra are indicated by the black error bars. The power spectra measured here for epochs 1 -5 are rainbow colored from red to blue. Magenta symbols indicate all epochs combined. The smallest scale power (dominated by white noise) decreases appropriately when all epochs are combined, whereas the large scale power is not strongly reduced by combining epochs. The large scale power is much higher than that seen in the deeper CDFS and HDFN observations, but is comparable to that of the shallower QSO1700 and EGS fields. Figure C4 shows the same comparisons at 4.5 µm. In this case, the shot noise level of the combined epochs is similar to that observed in the deep CDFS and HDFN fields. The large scale power is also close to the levels measured in these deep fields. Figure B2. Power spectra shown at various model depths for the combined and separate epochs, at 3.6 µm (left) and 4.5 µm (right). The black lines show that the large scale power does drop as the source model is pushed deeper. The red lines show the power spectra at "optimal" model depth when the skewness of the background (excluding outliers) is zero. The blue lines show the A − B power spectra. Figure C3. Comparison of the 3.6 µm power spectra measured here (red, orange, green, blue, violet, and magenta = Epochs 1, 2, 3, 4, 5, and all epochs combined), with the power spectra measured in 6 fields from Arendt et al. (2010) as indicated by the black error bars. Figure C4. As Figure C3, but for 4.5 µm power spectra.