Instrumental Response Model and Detrending for the Dark Energy Camera

The purpose of astronomical imaging is to infer the surface brightness distribution ${I}_{\star }(\theta ,\phi ,\lambda ,t)$ of light across celestial coordinates $(\theta ,\phi )$, wavelength λ, and epoch of observation t (we will take t as a discrete index labeling exposures). ${I}_{\star }$ has units of power per (area × solid angle × wavelength). For photons that are properly focused by the telescope, the astrometric calibration of the instrument maps celestial coordinates into indices ${\bf{x}}$ of pixels on the camera array. We use boldface to denote quantities that are vectors over the array pixels, i.e., they have an implicit discrete argument ${\bf{x}}$, and arithmetic operations on these vectors are assumed to be element-wise. The desired sky signal is ${{\bf{I}}}_{\star }(\lambda ,t)$. The astrometric solution also yields the solid angle subtended by each pixel, which we take to be a nominal value ${{\rm{\Omega }}}_{0}$ times a relative scaling vector ${\boldsymbol{\Omega }}$ over the array. For wide-field cameras, it is essential to distinguish between images representing the brightness of the sky at a given pixel versus images representing the flux received in a pixel; these differ by a factor of ${\boldsymbol{\Omega }}$. We will ignore the polarization state of the incident light and the polarization sensitivity of the instrument in this analysis.

If the focused sky photons were the only source of signal on the detector, then rate of production of photocarriers in the pixels could be written as

$$\begin{eqnarray}&&{\bf{Rate}}(t)={{\rm{\Omega }}}_{0}\int d\lambda \ \displaystyle \frac{\lambda }{{hc}}{{\bf{A}}}_{\mathrm{eff}}(\lambda ,t){{\bf{I}}}_{\star }(\lambda ,t){\boldsymbol{\Omega }}(\lambda ,t),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{{{\rm{\Omega }}}_{0}}{{f}_{1}\times 1\,{\rm{s}}}\int d\lambda \,{{\bf{I}}}_{\star }(\lambda ,t){r}_{\mathrm{ref}}(\lambda ){\bf{r}}(\lambda ,t){\boldsymbol{\Omega }}(\lambda ,t).\end{eqnarray} \tag{ 2 }$$

Here, ${{\bf{A}}}_{\mathrm{eff}}$ is the collecting area times the total transmission of the atmosphere and optics, times the quantum efficiency of the detector pixel. In the second line, we absorb the typical value of the product ${{\bf{A}}}_{\mathrm{eff}}\lambda /{hc}\times 1\,{\rm{s}}$ into an overall constant $1/{f}_{1}.$ A dimensionless reference bandpass ${r}_{\mathrm{ref}}(\lambda )$ is chosen to peak near unity and represent the spectral response of the system for a typical pixel and observing conditions, and the function ${\bf{r}}$ whose deviation from unity describes spatial and temporal deviations in the spectral response from the reference or "natural" bandpass of the camera. A source with flux f₁ when integrated over the reference bandpass generates one photocarrier per second under typical conditions. The units of f₁ are power/area.

For a point source, the incident surface brightness is spread over the array by a point spread function (PSF) which we normalize such that

$$\begin{eqnarray}&&{{\bf{I}}}_{\star }(\lambda ,t)=\displaystyle \frac{f}{{{\rm{\Omega }}}_{0}}{F}_{\star }(\lambda ){\bf{PSF}}(\lambda ,t),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{{\bf{x}}}{\bf{PSF}}(\lambda ,t){\boldsymbol{\Omega }}(\lambda ,t)=1,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\int d\lambda \,{F}_{\star }(\lambda ){r}_{\mathrm{ref}}(\lambda )=1.\end{eqnarray} \tag{ 5 }$$

Because imagers like DECam count photons without regard to their energy, a single observation can constrain only the amplitude f of the stellar flux, not the spectral shape ${F}_{\star }(\lambda ).$ The shape ${F}_{\star }$ has units of inverse wavelength.

Consider first the limited goal of determining f for sources which are known to share a reference source spectrum ${F}_{\mathrm{ref}}(\lambda ).$ We want to obtain the same result for f regardless of the time or focal plane position of the observation. To obtain f if given the ${\bf{Rate}}$ image, we first define the dimensionless reference flat

$$\begin{eqnarray}&&{\bf{r}}(t)\equiv \int d\lambda \,{F}_{\mathrm{ref}}(\lambda ){r}_{\mathrm{ref}}(\lambda ){\bf{r}}(\lambda ,t)\end{eqnarray} \tag{ 6 }$$

and construct the fluence image, with units of counts:

$$\begin{eqnarray}&&{\bf{Fluence}}(t)\equiv T\,\displaystyle \frac{{\bf{Rate}}(t)}{{\bf{r}}(t)},\end{eqnarray} \tag{ 7 }$$

where T is the exposure time. If we make the assumption that ${\bf{PSF}}$ and ${\boldsymbol{\Omega }}$ are independent of λ across the filter's bandpass and the width of the PSF, we can estimate the source flux by this sum over pixels spanning the object's image:

$$\begin{eqnarray}&&\left(\displaystyle \frac{1\,{\rm{s}}}{T}\right)\displaystyle \sum _{{\bf{x}}}{\bf{Fluence}}(t)\\ &&=\,\displaystyle \frac{f}{{f}_{1}}\displaystyle \sum _{{\bf{x}}}\displaystyle \frac{\displaystyle \int d\lambda {F}_{\mathrm{ref}}(\lambda ){r}_{\mathrm{ref}}(\lambda ){\bf{r}}(\lambda ,t){\bf{PSF}}(t){\boldsymbol{\Omega }}(t)}{\displaystyle \int d\lambda {F}_{\mathrm{ref}}(\lambda ){r}_{\mathrm{ref}}(\lambda ){\bf{r}}(\lambda ,t)},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{f}{{f}_{1}}\displaystyle \sum _{{\bf{x}}}{\bf{PSF}}(t){\boldsymbol{\Omega }}(t),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{f}{{f}_{1}}.\end{eqnarray} \tag{ 10 }$$

Thus, the reference flat is the quantity we need to homogenize the photometry of the survey. Common practice is to estimate the reference flat by the values in an image of a source of near-uniform surface brightness, i.e., a flat field. We will describe below why this is inaccurate—for DECam, dome flats deviate by up to ±5% from the reference flat—and in Section 3.4, we describe our method for estimating the reference flat.

Determination of ${\bf{r}}$ for the entire survey would constitute global relative calibration as it enables accurate determination only of the ratio of fluxes of two objects with the reference spectrum. A further absolute calibration step is needed to determine f₁ if we wish to place the fluxes on a known physical scale.

The output of the DES pipeline for exposure t is an estimate of the ${\bf{Fluence}}$ image, which, per Equation (8), is an estimate of the number of photons incident on the pixel during the exposure for a source with the reference spectrum. The sum over pixels yields a flux (and magnitude) estimate for each object to enter into object catalogs. The DES images and cataloged magnitudes are hence correct only for objects sharing the reference spectrum ${F}_{\mathrm{ref}}(\lambda ).$ For practical reasons, the reference spectrum for DES is taken to be that of the F8IV star C26202 from the HST CalSpec standard set.¹³ Note that DES conforms to the convention that calibrated images represent the flux (or, more precisely, fluence) in each pixel, not the surface brightness.

Consider now the aperture photometry for a source known to have a more general spectrum ${F}_{p}(\lambda )$ parameterized by some value(s) p. The flux estimate will be

$$\begin{eqnarray}&&\hat{f}={f}_{1}\times \displaystyle \frac{1\,{\rm{s}}}{T}\times \displaystyle \sum _{{\bf{x}}}{\bf{Fluence}}(t)\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&=\,f\displaystyle \sum _{{\bf{x}}}\displaystyle \frac{\int d\lambda {F}_{p}(\lambda ){r}_{\mathrm{ref}}(\lambda ){\bf{r}}(\lambda ,t){\bf{PSF}}(\lambda ,t){\boldsymbol{\Omega }}(\lambda ,t)}{\int d\lambda {F}_{\mathrm{ref}}(\lambda ){r}_{\mathrm{ref}}(\lambda ){\bf{r}}(\lambda ,t)},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&=\,f\times {\bf{C}}(t;p),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\bf{C}}(t;p)\equiv \displaystyle \frac{\int d\lambda {F}_{p}(\lambda ){r}_{\mathrm{ref}}(\lambda ){\bf{r}}(\lambda ,t)}{\int d\lambda {F}_{\mathrm{ref}}(\lambda ){r}_{\mathrm{ref}}(\lambda ){\bf{r}}(\lambda ,t)}.\end{eqnarray} \tag{ 14 }$$

The last line defines the color correction that must be divided into the cataloged flux in order to obtain the correct amplitude f of the source spectrum ${{fF}}_{p}(\lambda )$. In obtaining Equation (13), we assumed that ${\bf{C}}$ is constant across the pixels that comprise the object image, and again the ${\bf{PSF}}$ and ${\boldsymbol{\Omega }}$ are independent of wavelength across the object.

Determination of ${\bf{C}}(t;p)$ most generally requires that we know not only ${F}_{p}(\lambda )$, but also the full response ${\bf{r}}(\lambda ,t)$ of the array over time, space, and wavelength. This general knowledge cannot be obtained solely from internal calibrations of on-sky DECam data.

The ${\bf{C}}$ corrections are, of course, the smallest and smoothest in p if we have made the sensible choice of reference response ${r}_{\mathrm{ref}}(\lambda )$ to be the system response under typical conditions, so that ${\bf{r}}(\lambda ,t)$ is near unity and the cataloged fluxes are in the "natural" system of the camera. In this case, the ${\bf{C}}$ corrections are doubly differential: the cataloged fluxes produced with the reference flat via Equation (8) are correct (${\bf{C}}=1$) if either the source has the reference spectrum or the instrument has the reference response. The departure of ${\bf{C}}$ from unity scales as the product of the (small) variations in camera spectral response times the (not always small) deviation of the source spectrum from the reference.

The color correction simplifies considerably if the sources are stars with color $g-i\lt 2$ (spectral type M0) and we assign the parameter $p\equiv g-i-{(g-i)}_{\mathrm{ref}},$ where ${(g-i)}_{\mathrm{ref}}\,=0.44$ is the color of our reference star C26202 in the AB-normalized natural bandpass of DECam. In this case, synthetic photometry using stellar spectral libraries (Li et al. 2016) indicates that, for the range of ${\bf{r}}(\lambda ,t)$ expected from variation in atmospheric and DES instrumental transmission, the color correction will be fit to $\lt 1$ mmag accuracy by the form

$$\begin{eqnarray}&&-2.5{\mathrm{log}}_{10}{\bf{C}}(t;p)=p\times {\bf{c}}(t).\end{eqnarray} \tag{ 15 }$$

Thus, for the population of not-too-red stars, the photometric response for a given filter is fully specified by the absolute calibration f₁, the reference flat ${\bf{r}}(t)$, and the color term ${\bf{c}}(t)$. We show in Section 3.4 that the last two of these are nearly fully constrained by internal calibration methods.

For sources with highly structured spectra such as cool stars, quasars, and supernovae, the ${\bf{c}}(t)$ maps are useful constraints on ${\bf{r}}(\lambda ,t)$ but not sufficient to determine ${\bf{C}}(t;p)$ to desired accuracy. In Section 3.4.3, we describe how DES models the full response ${r}_{\mathrm{ref}}(\lambda ){\bf{r}}(\lambda ,t)$ by combining internal comparisons of on-sky data, narrowband flat-field observations, and atmospheric monitoring. This allows estimation of ${\bf{C}}$ for an arbitrary specified source spectrum. This process is described fully in Burke et al. (2017). This full-bandpass ${\bf{r}}(\lambda ,t)$ synthesis is currently the primary method for color calibration in DES. The internally derived ${\bf{c}}(t)$ can be used as a crosscheck for the full synthesized ${\bf{r}}(\lambda ,t).$

The focused photons from astronomical sources ${{\bf{I}}}_{\star }$ comprise only some of the light detected by the camera. For ground-based imaging, these are outnumbered by light focused from background surface brightness ${{\bf{I}}}_{\mathrm{bg}}(\lambda ,t).$ The distinction between signal and background is in the eyes of the beholder: diffuse astrophysical emission from dust reflections or intracluster light might be considered a nuisance by those conducting photometry of more compact objects. In this section, we consider as background¹⁴ just the atmospheric emission (airglow) plus scattering of sunlight, moonlight, and terrestrial light in the atmosphere and from zodiacal dust. In cloudless conditions, these sources are very smooth across the sky. We cannot assume they are strictly constant in a DECam exposure, because the distances from the zenith, the Sun, and the Moon all vary enough across the 2-degree diameter of the DECam FOV to produce easily detected background gradients. Clouds can produce more highly variable background emission and absorption (Burke et al. 2015); we will discuss only clear-night processing in this paper.

At least several percent of the photons reaching the DECam focal plane have arrived after unwanted reflections or scattering. We will denote as ${{\bf{S}}}_{\star }(\lambda ,t)$ the pattern of photocarrier production by photons originating in extraterrestrial sources (other than the Sun) and arriving at the focal plane out of focus. More precisely, our stellar photometry will be normalized to the flux measured within a 6'' aperture. Any photons from astrophysical sources whose path brings them to the focal plane outside such an aperture relative to their true direction of origin should not be counted in our response functions and will be assigned to ${{\bf{S}}}_{\star }$. There are various mechanisms producing such misguided photons. Annular "ghosts" appear where light from compact sources arrives after stray specular reflections, particularly light that reflects from the focal plane and returns after a second reflection from a filter or lens surface. Large diffuse scattering halos also surround all sources, and photons from sources outside the FOV can scatter from insufficiently baffled telescope surfaces. Diffraction spikes can extend beyond our aperture. For bright stars, these effects become detectable and interfere with accurate measurement of focused light. We will collectively refer to ${{\bf{S}}}_{\star }$ as the "ghost" signal, and consider it the job of the analysis software, rather than the detrending process, to deal with these signals. In the current DESDM pipeline, regions of detectable ghost signals are usually just flagged as invalid.

More important to the detrending are the photocarriers produced by light scattered from the nearly Lambertian background sources. This signal can be thought of as the integral of all the ghosts and scattered light obtained from sources uniformly distributed across the sky. It will be present in every exposure and scale with the (nearly uniform) brightness of the background.

Our full model for the rate of photocarrier production in the camera is hence,

$$\begin{eqnarray}&&{\bf{Rate}}(t)=\displaystyle \frac{{{\rm{\Omega }}}_{0}}{{f}_{1}\times 1\,{\rm{s}}}\displaystyle \int d\lambda \,{\rm{\Omega }}\,{r}_{\mathrm{ref}}(\lambda ){\bf{r}}(\lambda ,t)\\ &&\times \,\{{{\bf{I}}}_{\star }(\lambda ,t)+{\tilde{{\bf{I}}}}_{\mathrm{ghost}}(\lambda ,t)+{{\bf{I}}}_{\mathrm{bg}}(\lambda ,t)[1+{{\bf{s}}}_{\mathrm{bg}}(\lambda ,t)]\}.\end{eqnarray} \tag{ 16 }$$

In this equation, we define ${\tilde{{\bf{I}}}}_{\mathrm{ghost}}$ to be the the spurious brightness signal that scattered source photons cause in detrended images, essentially ${{\bf{S}}}_{\star }/{\bf{r}}.$ We will henceforth ignore this term as the responsibility of the analysis software.

The term ${{\bf{s}}}_{\mathrm{bg}}(\lambda ,t)$ is the ratio of charge production from scattered background light to that from focused background light at a given detector pixel, wavelength and time. It can be only crudely predicted from the optics models for DES. We need to determine this signal empirically to subtract it from the ${\bf{Rate}}$ signal to desired accuracy.

The presence of background and scattered light in the images causes two complications. First, such signals must be removed by a sky-subtraction algorithm to isolate the ${{\bf{I}}}_{\star }$ signals we wish to measure. Less appreciated is that their presence invalidates the standard method of using dome flats to estimate the reference flat ${\bf{r}}(t)$.

2.3.1. Dome Flats

Common practice for visible and NIR image processing is to estimate the reference flat ${\bf{r}}(t)$ by measuring the camera response to a nearly Lambertian source of light observed on the same date: either a twilight sky, or the median night-sky signal, or an illuminated screen on the dome. We generically call this the dome flat. If we optimistically assume that the target source has perfectly uniform illumination with spectrum ${F}_{\mathrm{dome}}(\lambda ),$ Equation (16) gives the resultant signal as

$$\begin{eqnarray}&&{\bf{Dome}}(t)\propto {\boldsymbol{\Omega }}\int d\lambda \,{r}_{\mathrm{ref}}(\lambda ){\bf{r}}(\lambda ,t){F}_{\mathrm{dome}}(\lambda )[1+{{\bf{s}}}_{\mathrm{bg}}(\lambda ,t)].\end{eqnarray} \tag{ 17 }$$

There are several reasons why we should not expect this to yield the same pattern as the reference flat in Equation (6):

• 1.  
  The spectrum ${F}_{\mathrm{dome}}$ does not match ${F}_{\mathrm{ref}}.$ It is possible to mitigate this problem by using twilight flats and assigning a G-star spectrum times Rayleigh scattering spectrum as the reference (though this would not correspond to any real source's spectrum). For DECam, however, twilight flats are prohibited since the detectors are susceptible to permanent damage from over-exposure. Night-sky flats are very poor matches to any realistic ${F}_{\mathrm{ref}}$ in redder bands where airglow lines generate fringes because of Fabry–Perot oscillations in ${\bf{r}}(\lambda )$.
• 2.  
  Dome flats do not fill the telescope pupil in exactly the same way as distant sky sources, causing subtle differences in both the focused and scattered light patterns.
• 3.  
  The dome flat contains a factor of pixel area ${\boldsymbol{\Omega }}$, which we do not want in our reference flat if we are performing aperture photometry.
• 4.  
  The dome flat is contaminated by the scattered light signal ${{\bf{s}}}_{\mathrm{bg}}$.
• 5.  
  The polarization state of the dome illumination will not be the same as the true atmospheric and zodiacal backgrounds, and the instrument may have polarization sensitivity.

Even if the first two problems—spectral and pupil mismatches—could be mitigated, the last three problems are still present. It is impossible to distinguish pixel-area fluctuations and scattered light from true QE variations solely from an image of a Lambertian source. In Section 3.4, we describe our process of combining dome flats with internal calibration of stellar catalogs to generate better estimates of the reference flat ${\bf{r}}$.

2.3.2. Sky Subtraction

We need a sky subtraction operator H that we can apply to the ${\bf{Rate}}$ image of Equation (16) that removes the background terms, leaving us with an image obeying Equation (1) that can be used in Equation (8) to measure object fluxes. A good linear operator would split the parts of ${\bf{Rate}}$ that arise from sources versus background:

$$\begin{eqnarray}&&{\bf{H}}[{{\bf{Rate}}}_{\mathrm{bg}}]=0,\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{\bf{H}}[{{\bf{Rate}}}_{\star }]={{\bf{Rate}}}_{\star },\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{{\bf{Rate}}}_{\mathrm{bg}}\equiv {\boldsymbol{\Omega }}\int d\lambda \,{r}_{\mathrm{ref}}(\lambda ,t){\bf{r}}(\lambda ,t){{\bf{I}}}_{\mathrm{bg}}(\lambda ,t)[1+{{\bf{s}}}_{\mathrm{bg}}(\lambda ,t)],\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{{\bf{Rate}}}_{\star }\equiv {\boldsymbol{\Omega }}\int d\lambda \,{r}_{\mathrm{ref}}(\lambda ,t){\bf{r}}(\lambda ,t){{\bf{I}}}_{\star }(\lambda ,t).\end{eqnarray} \tag{ 21 }$$

The usual approach to sky subtraction is morphological, based on the assumption that ${{\bf{I}}}_{\mathrm{bg}}$ varies slowly with sky position, while no (interesting) part of the source signal does so. A high-pass filter or robust fitting of low-order polynomial to the image can then be used for ${\bf{H}}$. While ${{\bf{I}}}_{\mathrm{bg}}$ can be safely assumed to have little spatial structure, the function ${\bf{r}}(\lambda ,t)$ can have small-scale structure due to Fabry–Perot behavior, producing fringes in the ${{\bf{Rate}}}_{\mathrm{bg}}$ image. Likewise, the ${\boldsymbol{\Omega }}$ function of DECam (and most deep-depletion CCD cameras) has structure on scales of interest in the source image due to stray electric fields, e.g., the "tree rings" analyzed by Plazas et al. (2014). Furthermore, the scattered-light function ${{\bf{s}}}_{\mathrm{bg}}(\lambda )$ can have sharp features, e.g., when a stray reflection path places an image of the pupil nearly at the focal plane, as happens with the KPNO Mosaic camera (Jacoby et al. 1998).

A partial remedy is to divide the ${\bf{Rate}}$ image by a dome flat before applying the high-pass filter. This cancels out the ${\boldsymbol{\Omega }}$ and $(1+{{\bf{s}}}_{\mathrm{bg}})$ terms in ${{\bf{Rate}}}_{\mathrm{bg}},$ rendering the sky smooth again:

$$\begin{eqnarray}&&{\bf{H}}[{\bf{Rate}}]={\bf{HighPass}}\left[\displaystyle \frac{{\bf{Rate}}}{{\bf{Dome}}}\right]\times {\bf{Dome}}.\end{eqnarray} \tag{ 22 }$$

This is imperfect, as the dome flat spectrum does not exactly match the night-sky spectrum. For example, the dome illumination does not excite the fringe pattern like real airglow. A better background subtraction scheme is to replace ${\bf{Dome}}$ with a median night-sky flat in Equation (22), which by definition contains the same fringes and other structures as the ${{\bf{Rate}}}_{\mathrm{bg}}$ we are targeting.

This common approach is essentially equivalent to fitting and subtracting a scaled version of the median sky signal from the image. This too fails, however, if the background signal pattern varies in time, which it does, for example, as the relative contributions of airglow (and hence fringes), zodiacal light, moonlight, etc. change.

We choose to dispense with the morphological approach entirely, and instead exploit the knowledge that the background signal is determined by a handful of physical variables, such as the phase of the moon, the excitation state of the airglow, and the relative positions of the telescope, zenith, moon, and Sun. The charge collection rate should be expressible as

$$\begin{eqnarray}&&{\bf{Rate}}(t)={{\bf{Rate}}}_{\star }+\displaystyle \sum _{k}{b}_{{tk}}{{\bf{Sky}}}_{k}+{\bf{Noise}}(t),\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&=\,({\bf{Sparse}})+({\bf{Low}}\,{\bf{Rank}})+({\bf{Noise}}),\end{eqnarray} \tag{ 24 }$$

where k runs over a small number K of sky template patterns. The second line restates the sky-subtraction problem as the standard problem of robust principal component analysis: the camera output, expressed as a matrix of dimension (pixel number) × (exposure number), is expected to consist of a sparse term (the astronomical objects, which cover a minority of the pixels), a low-rank term (the low-dimensionality background signal), and shot noise. Many algorithms exist to identify the low-rank components in this situation; we use a standard principal components analysis (PCA) coupled with outlier rejection iteration. Within DES, we have the luxury of thousands of exposures taken in each filter, from which we can derive the series of sky templates ${{\bf{Sky}}}_{k}.$ The sky subtraction process then consists of identifying the coefficients b_tk of the templates that best describe exposure t. The processes for doing so in the presence of the noise and sparse sky signal are detailed in Appendix A. Once we have these, the sky subtraction is simply

$$\begin{eqnarray}&&{\bf{H}}[{\bf{Rate}}]={\bf{Rate}}-\displaystyle \sum _{k}{b}_{{tk}}{{\bf{Sky}}}_{k}.\end{eqnarray} \tag{ 25 }$$

The common technique of subtracting a scaled median night-sky flat is equivalent to the special case K = 1. In practice, we find that four components are sufficient to describe the background in nearly all clear-weather DES exposures. The PCA method will fail if a majority of the pixels in the DES field of view are covered by astronomical objects, though the Galactic plane and cores of the Magellanic Clouds are likely the only sources visible to DECam that are large and crowded enough to cause this problem. The three to four brightest stars in the DES footprint cause large enough ghosts to confuse the PCA algorithm, but these images are largely irretrievable anyway. The morphology-based sky subtraction algorithm of the SExtractor code (Bertin & Arnouts 1996) that runs during the cataloging phase of the DESDM pipeline removes the sky variations from patchy clouds.

Dark current is negligible in DES images, so it is not included in the model.

2.4.1. Brighter-fatter Effect

The detector and signal chain alter the ${\bf{Rate}}$ values in several more ways before producing the output digital signal. When charge is collected by the CCDs, it is repelled by the charge that has been collected earlier in the integration. This has the effect of pulling the pixel boundaries inward to the charge, and may also increase diffusion of incoming charge near highly- occupied pixels. The net effect is to shift charge so as to broaden the PSF for brighter stars, hence the name "brighter-fatter effect." A model for this effect was proposed by Antilogus et al. (2014) and the parameters of this model were fit to DECam data by Gruen et al. (2015). Hence, we have an operator ${\bf{BF}}$ such that the charge actually collected in each pixel of exposure t is

$$\begin{eqnarray}&&{\bf{Charge}}(t)={\bf{BF}}[{\bf{Rate}}(t)\times T].\end{eqnarray} \tag{ 26 }$$

Gruen et al. (2015) also demonstrate that an algorithm that calculates and reverses the charge shift from the Antilogus model leads to at least a tenfold reduction of the effect both in theory and in practice, so that we have a practical inverse operator ${{\bf{BF}}}^{-1}.$

2.4.2. Gain, Nonlinearity, Crosstalk, and Bias

The photocarrier count is converted to ADU by two output amplifiers per CCD. The process is slightly nonlinear. Fortunately, the process appears to be local, meaning that the output value for a pixel depends only upon the charge collected in that pixel, not the values in any other neighboring pixels. The only exception is that there is some crosstalk between amplifiers at the level of ∼10⁻³ for intra-CCD amplifier pairs and ≲10⁻⁴ for inter-CCD pairs. The values output by the camera (the "raw image") are modeled as

$$\begin{eqnarray}&&{\bf{Raw}}(t)={\bf{XTalk}}\left[{\bf{Nonlin}}\left[\displaystyle \frac{{\bf{Charge}}(t)}{{\bf{Gain}}}\right]\right]+{\bf{Bias}}(t).\end{eqnarray} \tag{ 27 }$$

Gain is a constant per amplifier giving the number of photocarriers per output ADU in the linear response regime. The XTalk function takes the form of a matrix operation among groups of pixels that are read simultaneously through the 124 CCD amplifiers. We take the gain values, the nonlinearity function ${\bf{Nonlin}}$ for each amplifier, and the XTalk matrix elements to be fixed in time, as there is currently no evidence for change over time.¹⁵ Temporal changes in gain can, in any case, be absorbed into the reference flat ${\bf{r}}(t)$.

Charge transfer inefficiency (CTI) introduces non-local errors into CCD images. Fortunately, DECam is ground-based, hence shielded by Earth's atmosphere from cosmic rays that create charge traps in the device. Furthermore, DES is an exclusively broadband survey and all science exposures are 45–330 s duration, so the background level is high enough to render CTI unimportant. Observers using u-band or narrowband filters with DECam, or those using short exposures, may need to revisit CTI.

${\bf{Bias}}(t)$ is determined by the standard CCD techniques: a line-by-line subtraction of the mean overscan values is done first. We also create a Zero image from the robust mean of a series of zero-length exposures and subtract this from each image to remove any residual spatial structure. The Zero signal in the DECam CCDs is $\lt 15e$ and smoothly structured.

The nonlinearity of each amplifier is determined from a series of dome exposures of varying duration. Exposure times T between 0.05 and 35 seconds are selected, shuffled into random order, then interleaved with exposures of fixed duration T = 2 s. The median bias-subtracted ADU level of pixels read through a given amplifier is calculated at each exposure. The 2 s exposures are used to track a small linear drift in the dome illumination level, which is then divided out of all the signal levels. Nonlinearities in amplifier response are then manifested as nonlinearities in the plot of illumination-corrected ADU level vs exposure time. A ${{\bf{Nonlin}}}^{-1}$ function is chosen to map the ADU level back to a value proportional to exposure time.

The DECam results in Figure 1 show two classes of nonlinearity before the onset of saturation. All devices show a high-light-level nonlinearity consistent with a quadratic response term. A subset of devices show poorly understood nonlinear behavior at very low illumination levels in which some of the first few tens of electrons "disappear." One amplifier (#31B) shows unstable nonlinearity and is masked from all images as unsuitable for quantitative analyses. The worst remaining amplifier "loses" $\approx 20e$ of charge. The ${{\bf{Nonlin}}}^{-1}$ function does linearize this response, though one should be cautious about precision photometry at background levels below a few hundred e because the effect is poorly understood.

Zoom In Zoom Out Reset image size

Crosstalk is nominally a linear leakage from "source" amplifier to "victim" channel, hence there are 140² leakage coefficients to determine, as leakage may occur between any of the channels being read out synchronously with the valid science channels (including five non-science-quality amplifiers and 16 focus/alignment channels). Each coefficient was evaluated by plotting the median victim amplifier output versus source amplifier signal for a large number of sky images. Fortunately, only a handful of inter-CCD amplifiers show detectable ($\gtrsim {10}^{-5}$) crosstalk, and always at $\lt 2\times {10}^{-4}$ level. Crosstalk between science channels and the focus/alignment channels is small enough to ignore, so we have only the 124 science amplifiers in the crosstalk matrix.

Figure 2 shows the victim versus source trends for all of the intra-CCD amplifier pairs, which have crosstalk coefficients up to 10⁻³. It is seen that the crosstalk becomes nonlinear when the source amplifier exceeds its saturation level. These nonlinearities are tabulated and incorporated into the ${{\bf{XTalk}}}^{-1}$ operator. Nonlinearity is not significant for the already small inter-CCD crosstalk.

Zoom In Zoom Out Reset image size

The model for the brighter-fatter effect, the determination of its parameters, and the correction algorithm are fully described in Gruen et al. (2015). To summarize, the left, right, top, and bottom ($d\in \{L,R,T,B\}$, respectively), edges of pixel ij on the detector are shifted outwards by a distance determined by the charges q in neighboring pixels:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{ij}}^{d}=\displaystyle \sum _{{mn}}{a}_{{mn}}^{d}{q}_{i+m,j+n}.\end{eqnarray} \tag{ 28 }$$

The charge in pixel ij is altered as

$$\begin{eqnarray}{q}_{{ij}} & \to & {q}_{{ij}}+{{\rm{\Delta }}}_{{ij}}^{T}\displaystyle \frac{{q}_{{ij}}+{q}_{i,j+1}}{2}+{{\rm{\Delta }}}_{{ij}}^{B}\displaystyle \frac{{q}_{{ij}}+{q}_{i,j-1}}{2}\\ & & +{{\rm{\Delta }}}_{{ij}}^{L}\displaystyle \frac{{q}_{{ij}}+{q}_{i-1,j}}{2}+{{\rm{\Delta }}}_{{ij}}^{R}\displaystyle \frac{{q}_{{ij}}+{q}_{i+1,j}}{2}.\end{eqnarray} \tag{ 29 }$$

As the a coefficients are small, the effect can be reversed to good approximation by simply changing the + signs in (29) to minus signs. The correction is fast to calculate because each Δ array is a convolution of the Charge image with the kernel(s) a_mn. The convolution is executed by fast Fourier transform; then the charge shifting is a straightforward vector operation on the pixels.

The a coefficients are the parameters of the model. These are determined by measuring the covariance induced by the brighter-fatter effect between the noise on nearby pixels of dome flat images. These covariances are measured to high precision from 1200 r-band dome flats taken over a full season. The covariances (and the a coefficients derived from them) are observed to be constant in time and identical for the A and B amplifiers on a given CCD, as expected for an effect deriving from properties of the wafer and the gate structures.

Another manifestation of the brighter-fatter effect is the suppression of variance in dome flat images. This invalidates the usual method of determining the amplifier ${\text{Gain}}$, which is to note that the noise variance V expected for a pixel with expectation value μ (measured in ADU² and ADU, respectively) should be the sum of Poisson noise and read noise:

$$\begin{eqnarray}&&V=\displaystyle \frac{\mu }{{G}{ain}}+{(\mathrm{RN})}^{2}.\end{eqnarray} \tag{ 30 }$$

${\text{Gain}}$ is typically determined by regressing V against μ for a series of dome flats of varying exposure time. The brighter-fatter effect, however, introduces a term $\sim a{\mu }^{2}$ into this equation. The gain is overestimated by 5% to 10% if the data are fit without this term in the case of DECam CCDs. We determine a single gain per amplifier by a fit that includes the brighter-fatter effect.

We choose to factor the reference flat image as

$$\begin{eqnarray}&&{\bf{r}}(t)={\bf{Dome}}\times {\bf{SFlat}}\times {\bf{GRC}}(t).\end{eqnarray} \tag{ 31 }$$

The star flat ${\bf{SFlat}}$ is a parametric function of array coordinates that serves to correct the errors in the dome flat. ${\bf{SFlat}}$ is held constant for a given filter and observing epoch. The time dependence of ${\bf{r}}(t)$ within an epoch is captured by the global relative calibration function ${\bf{GRC}}(t)$, that is nearly constant across the field of view, but which can vary for each exposure, correcting for atmospheric extinction, mirror-coating aging, etc. There are many fewer parameters in ${\bf{SFlat}}$ and ${\bf{GRC}}$ than there are pixels in the array, and we can solve for these parameters using purely internal constraints.

3.4.1. Dome Flat Usage

We noted that the ${\bf{Dome}}$ image is a poor approximation to the response ${\bf{r}}(t)$ of the array to focused light from a source with the reference spectrum, as it confuses pixel-area variation and scattered light with the desired measure of response to focused light. Not surprisingly, the best way to calibrate the camera's response to focused starlight is to measure focused starlight, not diffuse light. Early recognition of this principle led to scanning a single source across the (small) arrays of the time (McLeod et al. 1995; Manfroid 1995), but we can solve for ${\bf{r}}(t)$ much more efficiently using the many sources detected by wide-field arrays. If a population of objects with the reference spectrum were spread across the sky, we could determine the reference flat empirically by finding a model of ${\bf{r}}(t)$ that produces uniform fluxes for all exposures of a given reference source. The DES photometric calibration relies primarily upon this internal calibration method (dubbed "ubercalibration" by Padmanabhan et al. (2008)). Note that the internal calibration does not actually require us to know either the reference source spectrum ${F}_{\mathrm{ref}}(\lambda )$ or the reference bandpass ${r}_{\mathrm{ref}}(\lambda )$ to determine ${\bf{r}}(t)$, nor does it offer information on either.

The success of this method depends upon choosing an observation sequence that can constrain, without degeneracy, a model for ${\bf{r}}(t)$ that adequately describes its variations. Unfortunately, the stellar data in DES are insufficient to estimate ${\bf{r}}$ at all 500 million DECam pixels, never mind doing so along with time dependence. Even with a huge number of stellar measurements, the fact that stellar images span several pixels precludes the use of stellar data to map response down to single-pixel scales. Therefore, we choose to rely on the dome flats to provide small-scale response information for which a parametric form is not readily constructed.

When the small-scale dome flat variations are truly QE effects, this is a good idea, such as for the "spots" visible in Figure 3. However, it is detrimental in that pixel-to-pixel variations in subtended sky area, due to imperfections in the CCD lithography, are mistakenly interpreted as QE variations. The latter probably lead to $\approx 1$ mas rms astrometric errors for DECam. In fact, it may be more accurate to replace the ${\bf{Dome}}$ factor in Equation (31) with a constant, except for a few regions with blemishes that are clearly QE features. For future cameras, these small-scale variations should be determined by laboratory characterization of the devices before deployment (e.g., Baumer & Roodman 2015, for LSST).

Zoom In Zoom Out Reset image size

3.4.2. Star Flats

The choices of ${\bf{SFlat}}$ parameterizations, and the assignment of parameters, are described in detail in G. M. Bernstein et al. (2017, in preparation). We give a short summary here.

The ${\bf{SFlat}}$ parameterization is expressed as a sum of terms in magnitude space, i.e., multiplicative in ${\bf{r}}.$ We begin with independent polynomial variation across each CCD in the mosaic, which essentially replace the large-scale behavior of ${\bf{Dome}}$ with patterns constrained by stellar data. The ${\bf{SFlat}}$ model also includes terms to compensate for known pixel-area variations at the edges of the devices ("glowing edges") and in annular patterns related to inhomogeneities in the silicon wafers ("tree rings") (Plazas et al. 2014).

The parameters of ${\bf{SFlat}}$ are chosen to optimize internal agreement of stellar photometry in a series of $\approx 20$ exposures targeted on a rich star field. The exposures are dithered with spacings from 10'' up to the 1° radius of the DECam field, to generate constraints on ${\bf{SFlat}}$ across this range of spatial scales, using several $\times {10}^{5}$ stellar measurements per filter. These "star flat" sequences are taken in each filter every few months during bright time. The ${\bf{SFlat}}$ for the star flat sequence is taken as the instrumental response for the observing epoch. For processing of the star flat sequences, we model ${\bf{GRC}}(t)$ with a single free extinction constant per exposure.

The internal calibration procedure can, if applied to all stars within the well-behaved range of color p, also derive the spatial and temporal variation of the color term ${\bf{c}}(t)$. As with the reference flat ${\bf{r}}$, the color term is modeled with a static function of array coordinates, plus a per-exposure term c_G(t) that is constant across the array. All stars in the well-behaved range are used to constrain ${\bf{SFlat}}$ and ${\bf{c}}$ simultaneously.

The images are first flattened using ${\bf{Dome}}$, and fluxes for each star are extracted using aperture photometry. The flux ${f}_{\alpha t}$ of star α in image t is converted to ${m}_{\alpha t}={m}_{0}-2.5{\mathrm{log}}_{10}\ {f}_{\alpha t}$, with some nominal zeropoint m₀. If ${\sigma }_{\alpha t}$ is the measurement uncertainty on ${m}_{\alpha t}$, then we seek to minimize

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{\alpha }\displaystyle \sum _{t}\\ &&\quad \times \,\displaystyle \frac{{[{m}_{\alpha t}-S({{\bf{x}}}_{\alpha ,t})-{p}_{\alpha }c({{\bf{x}}}_{\alpha t})-G(t)-{p}_{\alpha }{c}_{G}(t)-{m}_{\alpha }]}^{2}}{{\sigma }_{\alpha t}^{2}}.\end{eqnarray} \tag{ 32 }$$

We define $S=-2.5{\mathrm{log}}_{10}\,{\bf{SFlat}},$ and likewise the scalars G(t) and c_G(t) are the logarithms of ${\bf{GRC}}$ and a color term for each exposure. The observed pixel position, ${{\bf{x}}}_{\alpha t}$, and color, p_α, of star α are known, and its true magnitude, m_α, is a free parameter. There are $O({10}^{3})$ free parameters in the S and c functions and the G and c_G values.

Figure 3 shows the dome flat, star flat, and their product, the reference flat, for one of the DECam filters in one of the observing epochs. The second row of the figure indicates that the effects measured by the star flat—scattered light, pixel-area variations, and color mismatches between dome lamps and stars—are generally smooth across the array and dominated by a scattered light "doughnut" with a peak-to-peak amplitude of 4%. The structure remaining in the reference flat is dominated by variations in QE from CCD to CCD, with $\lt 0.01$ mag variation within each device. The zoomed images show that most of the "tree-ring" pattern in the domes is suppressed by the star flat, as expected from pixel-size fluctuations, although we do not understand why the tree rings still remain to some extent in the reference flat. The corner "tape bump" features that remain in the reference flat arise at the locations of small spacers that define the thickness of a glue layer between the CCD and its carrier. It is thought that stresses in the CCD induce electric fields that produce pixel-size variations in these regions, but our star flat parametric models do not have enough freedom to constrain these variations or remove them from the reference flat. We currently have no means to determine whether the remaining fine-scale features in the reference flat are from QE variations or pixel-size variations, although we suspect primarily the latter.

The optimized ${\chi }^{2}$ value is invariant under the transformations $G\to G+{G}_{0},$ ${c}_{G}\to {c}_{G}+{c}_{0}$ and ${m}_{\alpha }\to {m}_{\alpha }+{G}_{0}\,+{p}_{\alpha }{c}_{0},$ reflecting our inability to determine the absolute magnitude or color scale with purely internal calibrations. External absolute calibrators are needed to set the flux and color scales. There are other degeneracies in the internal calibration process; further details on the derivation of star flats are in G. M. Bernstein et al. (2017, in preparation).

3.4.3. Global Relative Calibration and Spectral Response Function

The above processing of specialized "star flat" observation sequences determines the ${\bf{SFlat}}$ vector for a given epoch and filter. When multiplied by the dome, it yields the ${\bf{r}}$, which homogenizes the measurements of stellar fluxes across the array at the time the star flat data were acquired. In the DES pipeline, this ${\bf{r}}$ is used to flatten all images obtained during the epoch, and stellar fluxes are derived from these images.

The next step is to derive the global relative calibration function, ${\bf{GRC}}(t)$, which can be applied to these cataloged fluxes to homogenize fluxes over all t, ensuring that a given star would yield the same measured flux when placed at any surveyed sky location at any time during the survey. Starting with Y3A1 reductions, this is done for DES survey exposures by the forward global calibration module (FGCM) described in Burke et al. (2017), which in fact yields a model for the chromatic response function ${\bf{r}}(\lambda ,t)$ by combining a static model of the telescope/camera spectral response with a time-varying model of the atmospheric transmission, ${S}_{\mathrm{atm}}(\lambda ,t)$. The atmospheric transmission is in turn modeled with nightly values for the parameters for the major atmospheric optical processes: absorption by water vapor, ozone, and well-mixed molecules, plus Rayleigh scattering and aerosol scattering.

The atmospheric modeling is attempted only for cloudless nights. For data taken through clouds, ${\bf{GRC}}(t)$ is derived by finding a distinct zeropoint shift and color term that bring each CCD of each exposure into agreement with overlapping exposures taken on clear nights. This procedure can benefit from the on-site thermal-infrared all-sky imager RASICAM (Reil et al. 2014), which definitively indicates which nights are free of clouds.

As proposed by Stubbs & Tonry (2006), FGCM splits the system spectroscopic response into an atmospheric part (anything outside the dome) and an instrumental function. The atmospheric portion is taken as constant across the field on clear nights and the instrumental function is assumed constant for a given filter during an observing epoch. Both functions are characterized using a mix of auxiliary instrumentation and internal calibrations.

The instrumental response is measured by the DECal narrowband flat-field illumination system (Marshall et al. 2016). DECal directs the output of a monochromator to the dome flat screen via optical fibers and a projector system, which are designed to ensure that the illumination pattern is independent of wavelength. Furthermore, the pupil is close to uniform. A calibrated photodiode at the top ring of the telescope measures the brightness of dome illumination as the monochromator is swept across each filter.

Let ${\bf{DECal}}(\lambda )$ be the DECam output divided by the photodiode signal, so it gives the response to diffuse light at the top ring. If we multiply this by the atmospheric transmission ${S}_{\mathrm{atm}}(\lambda ,t)$, we should obtain the response of the instrument to celestial diffuse light. Because the DECal illuminators do not fill the telescope pupil in the exact same way as the night sky does, the DECal signals differs from the celestial one by a spatial function D. Specializing Equation (17) to a monochromatic source spectrum and inserting the atmospheric transmission to get the total system response, we get

$$\begin{eqnarray}&&{S}_{\mathrm{atm}}(\lambda ,t){\bf{DECal}}(\lambda )\propto {r}_{\mathrm{ref}}(\lambda ){\bf{r}}(\lambda ,t)[1+{{\bf{s}}}_{\mathrm{bg}}(\lambda )]\times {\boldsymbol{\Omega }}\times {\bf{D}},\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&\Rightarrow {r}_{\mathrm{ref}}(\lambda ){\bf{r}}(\lambda ,t)\propto {S}_{\mathrm{atm}}(\lambda ,t)\times \displaystyle \frac{{\bf{DECal}}(\lambda )}{[1+{{\bf{s}}}_{\mathrm{bg}}(\lambda )]\times {\boldsymbol{\Omega }}\times {\bf{D}}}.\end{eqnarray} \tag{ 34 }$$

DECal is designed to have the wavelength dependence of D be small, as is that of ${\boldsymbol{\Omega }}$. If we ignore the wavelength dependence of the scattered-light fraction ${{\bf{s}}}_{\mathrm{bg}}(\lambda )$ across the band of each filter, then the product of ${\bf{DECal}}$ and ${S}_{\mathrm{atm}}$ give the spectral response ${\bf{r}}(\lambda ,t)$ up to some time-independent function of pixels. The unknown function is determined by the star flat measurement of ${\bf{r}}(t).$

The scattered-light function, ${{\bf{s}}}_{\mathrm{bg}}(\lambda )$, is problematic because we have no means to distinguish scattered from focused light in the DECal flats, and unfortunately there is reason to suspect that it does vary within a band: stray light from filter reflections requires (at least) one transmission and one reflection by the filter, so scales as $T(1-T)$ if T is the filter transmission. This signal should peak dramatically as DECal scans across the filter shoulders. In the future, we may be able to model this effect and refine the instrument spectral models. Because the transition from 10% to 90% of peak transmission occurs within $\leqslant 12 \%$ of the bandwidth for each DES filter, we do not expect wavelength of variation of ${{\bf{s}}}_{\mathrm{bg}}(\lambda )$ to be significantly altering the color corrections, which are already small.

The nightly parameters of the atmospheric transmission, ${S}_{\mathrm{atm}}(\lambda ,t)$, are constrained by forcing agreement on magnitudes and colors between observations taken on different nights. The atmospheric constituents are further constrained by external data: the barometric pressure and airmass for every observation are recorded and determine the optical depths of the well-mixed molecular components. Synthetic photometry indicates that ozone variations will contribute less than 1 mmag to ${\bf{C}}(t;p)$ in the DECam bands, hence a typical value is assumed for all observations (Li et al. 2016). The ATMCam instrument (Li et al. 2014) continuously monitors a bright star through four narrow filters selected to constrain the remaining aerosol and water vapor components. The latter is also monitored by continuously by a high-precision GPS receiver.¹⁶ These instruments (which are not always in service) are supplemented by internal color constraints to estimate the remaining atmospheric values for each night.

With the ${\bf{DECal}}$ scans and the nightly atmospheric model in hand, we can integrate any source spectrum family to derive the spectral flux correction factor ${\bf{C}}(\lambda ,t;p)$. Recall that the star flat process yields an empirical color term, ${\bf{c}}(t)$, for stellar spectra. Stellar spectra are well known, so ${\bf{c}}$ is a known special case of ${\bf{C}}(t;p)$, hence ${\bf{c}}$ is a consistency check. For example, DECal scans detect a 6 nm center-to-edge change in the blue cutoff of the DES i-band filter (Marshall et al. 2016). Synthetic stellar photometry predicts a radial color term that agrees to better than 1 mmag with the values empirically measured from the star flats (Li et al. 2016).

The absolute flux scale f₁ of the DES natural system in each filter is determined through observations of spectrophotometric standard stars. Current estimates of f₁ are made using the HST CalSpec standard star C26202, which lies in one of the DES supernova fields (which are observed ≈ weekly) and is within the dynamic range of our imaging. A sample of DA white dwarf stars in the DES footprint is being measured and modeled for use as absolute calibrators. A more thorough treatment will be presented in a later publication.

The background subtraction technique requires derivation of the component templates ${{\bf{Sky}}}_{k}.$ This is done through a principal-components analysis (PCA) of $\approx 1000$ images taken in a given filter in a given observing epoch. Appendix A describes this procedure.

Figure 4 plots the four sky templates derived for the z band. Here, it is apparent that the dominant template (essentially the median sky signal) contains signatures of fringing as well as the scattered-light "doughnut," and a color difference from the dome flat. Weaker components capture gradients in the sky illumination, and changes in the ratio of airglow to continuum background.

Zoom In Zoom Out Reset image size

Determination of the map $(x,y,t)\leftrightarrow (\theta ,\phi )$ between pixel and sky coordinates is accomplished by an internal-calibration procedure nearly completely analogous to the photometric methods, using the same star flat exposure sequences. Degeneracies in the internal calibration are resolved by including stellar positions from the Gaia DR1 catalog in the analysis (The Gaia Collaboration 2016). The full DECam astrometric model and its derivation are described in Bernstein et al. (2017).

For each epoch, a static bad pixel mask (BPM) is created by noting pixels in the array which have unusually high, low, or noisy levels in the bias and/or dome flat exposures. We also mask the 15 pixels closest to each edge of each CCD, because distortions of the electric field near the device edges cause changes in pixel effective area that are too large and potentially background-dependent to be correctable to the desired science accuracy. Pixels with any of these flags set are ignored in all processing.

The BPM additionally flags as suspect some other array regions that are useful for most analyses, but that have subtle quirks that preclude calibration to the same accuracy as the bulk of the array. The suspect regions include the areas 16–30 pixels from the detector edge; the "tape bump" features on the CCD that have small but highly structured pixel-size variations, and some regions that acquire excess noise when transferring through traps or hot pixels on their way to the readout amplifier. Suspect pixels are detrended and carried through to final images in the normal way, but the flag persists so that the most precise analyses (e.g., weak gravitational lensing, parallax, and proper motion) can ignore these data. E. Morganson et al. (2017, in preparation) document the BPM flags and the mask generating process in more detail.

Our operator D should yield a decimated $\tilde{{\bf{R}}}$, which is altered by all of the physical parameters expected to alter the background. We expect changes in the large-scale pattern of focused and scattered background light to arise from changes in lunar phase and the positions of the Sun, Moon, and telescope. We divide the image into ${N}_{D}\times \ {N}_{D}$ pixel subarrays, and produce as $\tilde{{\bf{R}}}$ the medians of each subarray. Typically, we take N_D = 128. The resultant "mini-sky" image is not only smaller, but has greatly reduced noise and is more sparse than the original image. Objects smaller than N_D pixels are filtered away by the median, and large disturbances to the background (such as ghosts of bright stars, reflection nebulae, or intracluster light) occupy the same fraction of the compressed image's pixels as in the original.

A concern is that airglow fringing, one of the most important backgrounds to remove, will be largely suppressed by the median at scale ${N}_{D},$ such that PCA of $\tilde{R}$ cannot recover the elements of U that depend on the strength and excitation state of the airglow. We therefore take an 8 × 8 block median of one of the CCDs (which still resolves the fringes), and append to the data vector $\tilde{{\bf{R}}}$ a 2 × 2 subsampling of this image. The final feature vector $\tilde{R}$ has ${N}_{f}\approx {10}^{4.5}$ elements and the linear algebra operations are much faster than the disk access times for the exposures.

With a compressed $\tilde{R}$ that should manifest all changes to the background, we proceed to execute a robust PCA decomposition via the following algorithm.

• 1.  
  Begin a rank-one approximation of $\tilde{R}\approx {{\bf{uv}}}^{T}$ with vectors ${\bf{u}}$ and ${\bf{v}}$ by initializing v_k = 1.
• 2.  
  Set ${u}_{t}={{\rm{median}}}_{i}({\tilde{R}}_{{ti}}/{v}_{i}).$
• 3.  
  Set ${v}_{i}={{\rm{median}}}_{t}({\tilde{R}}_{{ti}}/{u}_{t}).$ Iterate with step 2 to convergence of ${\bf{u}}$ and ${\bf{v}}.$
• 4.  
  Rescale ${\tilde{R}}_{{ti}}\leftarrow {\tilde{R}}_{{ti}}/({u}_{t}{v}_{i})$ so all elements are near unity.
• 5.  
  Conduct a standard PCA of $\tilde{R}$ and retain the ${N}_{\mathrm{keep}}\sim 20$ most significant components in a trial decomposition $\tilde{R}\approx U{\tilde{V}}^{T}.$
• 6.  
  For each exposure t, define ${\sigma }_{t}$ as the rms value of the residual mini-image ${\tilde{{\bf{R}}}}_{t}-{\sum }_{k}{U}_{{tk}}{\tilde{{\bf{V}}}}_{k}.$
• 7.  
  For some chosen σ-clipping threshold g, replace all pixels having residual $\gt g{\sigma }_{t}$ with the value of the rank-${N}_{\mathrm{keep}}$ approximation. This step excludes super pixels that remain contaminated by sources and ghosts. Exposures with very large ${\sigma }_{t}$ are discarded from the process as ill-described by the low-rank decomposition.
• 8.  
  Iterate to step 5 until U converges, typically three to four iterations.
• 9.  
  Rescale U and $\tilde{V}$ with ${\bf{u}}$ and ${\bf{v}}$ to make them represent the low-rank decomposition of the original $\tilde{R}$ matrix. We further rescale U and $\tilde{V}$ such that the elements of each row of $\tilde{V}$ (which are decimated sky templates) have rms value of unity. The U_tk elements then give the typical amplitude (in count rate) of the contribution of template k to the background of exposure t.
• 10.  
  Select the number K of templates to retain in the sky model by noting where the amplitudes U_tk of the corresponding corrections become insignificant.

With U in hand, Equation (43) separates into an independent robust linear regression for each pixel i with unknowns V_ki:

$$\begin{eqnarray}&&{R}_{{ti}}={S}_{{ti}}+\displaystyle \sum _{k=1}^{K}{U}_{{tk}}{V}_{{ki}}+{N}_{{ti}}.\end{eqnarray} \tag{ 46 }$$

We solve this system for each i using a standard linear regression with aggressive σ-clipping of outliers, as we expect a substantial fraction of the ${N}_{\exp }\approx {10}^{3}$ samples at this pixel to be "contaminated" by object flux. This is the most time-consuming portion of the analysis because there are 5 × 10⁸ such solutions to execute per DECam filter. However, this operation is trivially parallelized and the computation time is insignificant compared to the DESDM pipeline execution time.

Once the sky template matrix V is known for a given filter during a given epoch, we can execute the sky subtraction on any exposure's data ${{\bf{R}}}_{t},$ whether or not it was among the ${N}_{\exp }$ used to derive the templates. This is again a robust linear regression operation, which we conduct on a decimated version of the exposure using the decimated templates $\tilde{V}:$

$$\begin{eqnarray}&&{[D({{\bf{R}}}_{t})]}_{i}={\tilde{S}}_{t}+\displaystyle \sum _{k=1}^{K}{u}_{{tk}}{\tilde{V}}_{{ik}}+{\tilde{N}}_{{ti}}.\end{eqnarray} \tag{ 47 }$$

We solve for the K coefficients u_tk which best satisfy this equation, once again using a σ-clipping iteration to exclude super pixels that are perturbed by large-scale objects or image artifacts.

This fit to the mini-sky image yields valuable diagnostic information as well as the background coefficients u_tk: the fraction ${f}_{\mathrm{clip}}$ of super pixels that were clipped away; the rms variation σ of those super pixels that remained; and the residual of the mini-sky image to the best sky model. A high ${f}_{\mathrm{clip}}$ usually indicates the presence of a large spurious-light source in the image, such as an airplane trail, or ghosts of a bright star. A high σ value indicates an aberrant background pattern, such as can occur for images taken through patchy clouds, or with a very bright, nearby moon. These artifacts are readily identified by visual inspection of the mini-sky residuals.

For DECam, we find that $\gt 99 \%$ of the background signal is described as a linear combination of $K\leqslant 4$ sky components (except for exposures taken through patchy clouds). Figure 4 shows the first four principle components of background in the DECam z filter. The top row shows the "mini-sky" decimated images for each component row of $\tilde{V}$. The second row shows a closeup of a single CCD to reveal the presence of fringing signals in the full-resolution V. All of the signals here are shown relative to the dome flat pattern.

The dominant PC1 is basically the median night-sky signal. The upper left panel reveals that the background signal differs from the dome flat in several respects: first, there is a significant slope, suggesting that there is a gradient in the dome illumination system; second, the doughnut structure suggests that the scattered light amplitude differs between dome and sky illumination; and we also see CCD-to-CCD step functions, which are likely due to a spectral differences between dome and night sky combined with varying spectral response of the CCDs. The lower left panel reveals the fringing signal, which we expect to be present in the night sky but not in the domes. In the bluer filters there is less fringe signal, but one can see in the sky templates other small-scale deviations between dome and sky signals, which are likely attributable to color differences. A morphology-based sky subtraction would fail to identify such issues.

In all of the filters we find that PC2 and PC3 are essentially two orthogonal linear slopes of sky flux across the array. These are readily associated with gradients in brightness expected across the FOV as the distance to the Moon or Sun varies, or gradients in the zodiacal brightness. Note there are no fringes visible in PC2 or PC3: the airglow is not involved in these gradients. PC2 and PC3 coefficients are found to be usefully diagnostic of problems such as the occultation of the telescope beam by the observatory dome.

Fringes reappear in PC4 for z band, as does a doughnut pattern. Our hypothesis is that PC4 is generated by a change in the ratio of airglow intensity to the intensity of backgrounds that have nearly solar continuum spectra (twilight, moonlight, and zodiacal light). PC4 basically allows the fringe amplitude to vary independently of the large-scale patterns, but we see from the upper right panel that this is accompanied by a change in the overall illumination pattern. Apparently, the airglow/solar ratio change is manifested as an overall color shift of the night sky as well as in the fringing patterns. In fact, we found that the large-scale sky pattern is fully predictive of the fringe amplitude for this filter; we recover the fringe pattern in PC4 even if we omit from our decimation operator D the secondary small-scale feature vector described above. In practice, therefore, we have found the DECam PCA can be conducted just with the 128 × 128 median-filtered "mini-sky" images.

Principal components at $k\gt 4$ are found to contain smooth quadratic and increasingly higher-order polynomial variation across the focal plane, with no visible fringe or small-scale components. As such patterns are, at most, a few e amplitude in normal images, and are getting into the realm where they may be excited by astrophysical sources such as reflection nebulae or bright star ghosts, we truncate our templates at K = 4 and leave higher-order sky variation to the sky subtraction algorithm executed at the cataloging step.

In the $g,r,$ and i bands the fringes are nearly absent, but other small-scale structures (such as tree rings) are visible in PC1. The DECam detectors have low fringe amplitudes because their thick, deep-depletion bulk and good anti-reflection coatings yield high quantum efficiency to the silicon bandgap cutoff. This makes them weak Fabry–Perot resonators for the airglow lines. Other detectors with stronger fringe patterns may benefit even more than DECam from this PCA-based sky subtraction technique, because their fringe patterns will be stronger and may also exhibit detectable pattern shifts as excitation conditions of the ionospheric radicals change.