Astrometric Calibration and Performance of the Dark Energy Camera

The Dark Energy Camera (DECam; Flaugher et al. 2015) is one member of a new generation of high-throughput imagers combining a large field of view (FOV) (3 deg² for DECam) with a large telescope aperture (the 4-meter Blanco telescope). In the post-Gaia era (The Gaia Collaboration 2016), when positions, proper motions, and parallaxes are expected to be available with $\lt 1$ milliarcsecond (mas) accuracy for the 10⁹ stars with magnitude $G\lesssim 20$, what need do we have for accurate astrometry from these large ground-based cameras? There remain strong scientific motivations to obtain the most accurate possible positions for sources fainter than Gaia's limit, for transient sources, and for solar system bodies. Ideally these general-purpose, large-format imagers would be capable of obtaining astrometric measurements limited by the unavoidable shot noise and atmospheric fluctuations. Motivation and practice of astrometry from large-format ground-based CCD cameras have been discussed by Anderson et al. (2006), Platais et al. (2006), Bouy et al. (2013), and Magnier et al. (2016), among others. Accurate astrometry underlies many of the science goals of the Large Synoptic Survey Telescope (LSST) now under construction (LSST Science Collaborations 2009).

In addition, one of the motivators for construction of DECam is the measurement of weak gravitational lensing distortions of galaxies. Success in this pursuit requires the ability to register multiple exposures of every galaxy to an accuracy of ≈10 mas or better, otherwise the blur induced by misregistration in combining images could be mistaken for a coherent weak-lensing distortion. Searches for transient sources also benefit from precise image registration to improve subtraction of static sources.

One thing we do not need our wide-field imagers to do is determine absolute positions, as the preliminary Gaia catalogs are sufficiently dense to yield thousands of stars in the spatial and dynamic range overlap between Gaia and most DECam exposures. These suffice to determine the absolute pointing and any low-order astrometric distortion terms across the DECam FOV. In this work, we therefore focus on establishing relative astrometry with DECam on scales $\lesssim 1^\circ$. Indeed, one might ask why to bother at all with the effort making an astrometric model for DECam instead of simply interpolating all positional errors from Gaia stars. First, many of the detector-level effects occur on angular scales too small for Gaia stars to sample. Second, if our model removes discontinuities in the astrometric errors between CCDs, we can interpolate using reference stars from the whole field rather than being confined to those on a single device. Furthermore, the Gaia proper motion catalog is not yet available, so the reference catalog is not yet at mas accuracy. Lastly, many DECam exposures may have dynamic range which does not overlap well with the Gaia catalog.

In Section 2, we describe our method of deriving the DECam astrometric map and its error properties by forcing internal agreement among stellar positions in a series of offset exposures of rich star fields. In Sections 3, 4, and 5, we describe the data used to characterize DECam astrometry, the model applied to it, and the static residuals to this model, i.e., those which repeat from exposure to exposure. Section 6 characterizes the stochastic residuals, i.e., those uncorrelated between consecutive exposures and presumably due to atmospheric fluctuations. Section 7 characterizes the changes in the astrometric model from night to night and over the first 4 years of DECam observations. Section 8 investigates how much of the stochastic distortion can be removed by interpolation from a set of reference stars of a given density.

Our goal will be to model any astrometric distortion that contributes more than ≈1 mas rms error that is correlated between stars separated by $\gt 10^{\prime\prime}$. To put this scale in context, note that the mean scale of a 15 μm DECam pixel is 264 mas, so 1 mas corresponds to 0.004 pixel = 60 nm, or about 100 atoms in the silicon lattice. The DECam science array consists of 62 deep-depletion CCDs, each 2048 × 4096 pixels, and the array spans a roughly hexagonal area of diameter 2°. Thus 1 mas is 1.4 parts in 10⁷ of the DECam FOV. Scale changes due to stellar aberration, air pressure variations, and atmospheric refraction are far larger than this, so we will clearly need to allow each exposure an independent overall linear transformation across the FOV to approach mas accuracy. Indeed the nonlinear portion of atmospheric refraction is expected to have peak-to-peak amplitude of $(14\,\mathrm{mas})\times {\sec }^{2}z\tan z$ across the DECam FOV (where z is the zenith angle), so we must allow at least quadratic freedom to our solution atop any static instrumental model.

RMS positional errors reported in this paper refer to the sum of E-W and N-S components, unless noted otherwise.

The astrometric solution for DECam is a parametric model for the celestial (world) coordinates ${{\bf{x}}}^{w}=({x}^{w},{y}^{w})$ of an object given its observed pixel coordinates ${{\bf{x}}}^{p}=({x}^{p},{y}^{p})$ and some set of observing circumstances, C, which might include the object's color, c, plus discrete variables such as the date, exposure, filter, and individual CCD on which the measurement was obtained. The solution is found by straightforward ${\chi }^{2}$-minimization over the values of the model parameters ${\boldsymbol{\pi }}$. The index, i, ranges over all position measurements used to constrain the solution, and we assume a measurement error, ${\sigma }_{i}^{p}$, that is the same for both positional components of ${{\bf{x}}}_{i}^{p}$. We index the distinct objects on the sky by α, let ${\alpha }_{i}$ be the object targeted by measurement i, and denote by $i\in \alpha$ the subset of measurements with ${\alpha }_{i}=\alpha$. We define

$$\begin{eqnarray}&&{\chi }^{2}\equiv \displaystyle \sum _{i}{w}_{i}{| {{\bf{x}}}^{w}({{\bf{x}}}_{i}^{p},{C}_{i},{\boldsymbol{\pi }})-{\bar{{\bf{x}}}}_{{\alpha }_{i}}| }^{2},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\bar{{\bf{x}}}}_{\alpha }\equiv \displaystyle \frac{{\displaystyle \sum }_{i\in \alpha }{w}_{i}{{\bf{x}}}^{w}({{\bf{x}}}_{i}^{p},{C}_{i},{\boldsymbol{\pi }})}{{\displaystyle \sum }_{i\in \alpha }{w}_{i}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{w}_{i}^{-1}=({\sigma }_{i}^{2}+{\sigma }_{\mathrm{sys}}^{2}){\left|\displaystyle \frac{d{{\bf{x}}}^{w}}{d{{\bf{x}}}^{p}}\right|}_{i},\end{eqnarray} \tag{ 3 }$$

where the last term is the determinant of the Jacobian. In Equation (3) we introduce ${\sigma }_{\mathrm{sys}}$ to prevent very high weights from being assigned to high signal-to-noise ratio (S/N) measurements. We may consider ${\sigma }_{\mathrm{sys}}$ to represent the expected stochastic position errors beyond those arising from image noise that are included in ${\sigma }_{i}.$ In practice, we find that stochastic atmospheric distortions dominate the astrometric residuals, so we set ${\sigma }_{\mathrm{sys}}$ near the typical rms atmospheric distortion in our data (∼10 mas, cf. Section 6). This fairly arbitrary choice appropriately equalizes the weights assigned to individual measurements, but it does mean that our final ${\chi }^{2}$ values should not be expected to follow a ${\chi }^{2}$ distribution: it serves only as the quantity used to optimize ${\boldsymbol{\pi }}$.

Note that measurements from DECam can be freely mixed with other instruments' position measurements in Equation (1). Internal constraints—that multiple DECam observations of a source yield the same world coordinates—are combined with external constraints that DECam match a source of a priori assigned world coordinates of these objects. We denote as a reference catalog any catalog of objects whose map to celestial coordinates has no free parameters to adjust. The canonical case is an external astrometric catalog such as Gaia, for which the catalog directly specifies ${{\bf{x}}}_{i}^{w}$ and errors ${\sigma }_{i}^{w}$ in celestial coordinates, i.e., the function ${{\bf{x}}}^{w}({{\bf{x}}}^{p})$ is simply the identity.

Our strategy for DECam calibration is to produce very strong internal constraints by taking a series of $\approx 20$ consecutive exposures of fields at modest Galactic latitudes, where stellar sources are abundant but not crowded. The pointings of these exposures are shifted by anywhere from 10'' to the FOV diameter, so that a given star is imaged at many places on the array. In this scheme, the reference catalog serves mainly to break degeneracies in overall position and linear scaling of the astrometric map (see Section 2.3). These sequences of exposures are called star flats (as they are also used to calibrate photometric response). Since DECam was installed in 2012, star-flat sequences in all filters have been executed several times per year, usually during bright time. These data, described in Table 1, are the ones used in this paper to derive the DECam astrometric model.

Table 1.  Star-Flat Observing Sequences and DECam Thermal Events Through 2016 September. Italic Entries Are Camera Thermal Events. Roman Entries Are Observing Sequences

Epocha	Field	D50b	Airmass
20121120 c	0640–3400	209	1.04
20121223	0730–5000	104	1.06
2012 Dec 30	Camera warmup
20130221	1327–4845	112	1.06
2013 May 12	Camera warmup
2013 July 22	Camera warmup
20130829	1900–5000	110	1.07
2013 Oct 15	Camera warmup
20131115	0640–3400	141	1.09
2013 Nov 30	CCD S30 fails
20140118	1327–4845	133	1.33
2014 May 12	Camera warmup
20140807 d	1327–4845	143	1.32
20141105	0640–3400	128	1.01
2014 Dec 1	Camera warmup
20150204	1327–4845	088	1.31
2015 May 25	Focal plane temperature drop
2015 June 25	Partial camera warmup
2015 July 25	Camera warmup
2015 Aug 9	Camera warmup
2015 Aug 25	Camera warmup
20150926	2040–3500	119	1.01
2015 Nov 26	Focal plane temperature drop
20160209	0730–5000	125	1.07
2016 Feb 19	Camera warmup & corrector lens cleaning
20160223	1327–4845	110	1.24
20160816	1900–5000	108	1.06

Notes.

^aThe local date at start of the night when the star-flat exposures were taken or event occurred. ^bMedian half-light diameter of the point-spread function for the i-band exposures in the sequence. ^czY star flats were taken on the following night. ^dzY star flats were taken on November 10.

Download table as:  ASCIITypeset image

In the remainder of this section, we detail the algorithmic and coding choices made in defining the maps ${{\bf{x}}}^{w}({{\bf{x}}}^{p})$ and in the minimization of ${\chi }^{2}.$ These choices are realized in C++ code with an executable program called WcsFit. A reader uninterested in the implementation details can skip to Section 4.

2.1. Terminology

2.2. Available Maps

We must specify a functional form (and free parameters) for a map ${{\bf{m}}}_{k}({{\bf{x}}}^{p},c)$ from pixels to world coordinates for each extension k. Here c is the object color, and the other elements of the observational circumstances C are specified by the extension index. WcsFit allows each map ${{\bf{m}}}_{k}$ to be specified as the composition of a series of $j=1,2,\ldots ,{N}_{k}$ "atomic" coordinate transformations ${{\bf{m}}}_{k}^{(i)}$:

$$\begin{eqnarray}&&{{\bf{x}}}^{(0)}={{\bf{x}}}^{p},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{\bf{x}}}^{(i)}={{\bf{m}}}_{k}^{(i)}({{\bf{x}}}^{(i-1)}),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\bf{x}}}^{w}={{\bf{x}}}^{({N}_{k})}.\end{eqnarray} \tag{ 6 }$$

We generically refer to the input of each transformation ${{\bf{m}}}_{k}^{(i)}$ as its "pixel" coordinates and the output as its "world" coordinates, even though the intermediate variables are in fact neither. In our application, the chain of component maps is divided into those of the instrument solution followed by those of the exposure solution.

WcsFit follows the definitions in Section 2.1 by making each coordinate transformation or element thereof an instance of an abstract C++ base class PixelMap. Each has a type, a unique name string, and has a number $\geqslant 0$ of free parameters controlling its actions. PixelMap instances can be (de-)serialized (from) to ASCII files in YAML format, easily read or written by humans. The WcsFit user specifies the transformations to be fit to the data by giving the program such a YAML file as input; the parameters are assigned default starting values if none are specified. Anywhere that the strings BAND, INSTRUMENT, EPOCH, or DEVICE appear in these input files they are replaced with the values appropriate to the extension, allowing a generic model to be specified compactly. The WcsFit user can also specify the names of any PixelMaps whose parameters should be held fixed at their input values. The primary output of WcsFit is another YAML file specifying all of the maps and their best-fit parameters.

The types of PixelMaps available for use are

• The Identity map, which leaves ${\bf{x}}$ unchanged, and has no free parameters;
• Constant maps have ${{\bf{x}}}^{w}={{\bf{x}}}^{p}+{{\bf{x}}}_{0},$ with the two components of ${{\bf{x}}}_{0}$ as parameters;
• Linear maps have ${{\bf{x}}}^{w}=A{{\bf{x}}}^{p}+{{\bf{x}}}_{0},$ with six parameters in ${{\bf{x}}}_{0}$ and the components of the matrix A;
• Polynomial maps have their free parameters as the coefficients of two polynomials of specified degrees d_x and d_y in the ${{\bf{x}}}^{p}$ components that produce x^w and y^w, respectively; and
• Template maps apply transformations based on lookup tables. One has the option of x, y, or radial transformations:
  $$\begin{eqnarray}&&{x}^{w}={x}^{p}+{s}\,{f}({x}^{p}),\end{eqnarray} \tag{ 7 }$$
  $$\begin{eqnarray}&&{y}^{w}={y}^{p}+{s}\,{f}({y}^{p}),\mathrm{or}\end{eqnarray} \tag{ 8 }$$
  $$\begin{eqnarray}&&{{\bf{x}}}^{w}={{\bf{x}}}^{p}+s\displaystyle \frac{{{\bf{x}}}^{p}-{{\bf{x}}}_{c}}{| {{\bf{x}}}^{p}-{{\bf{x}}}_{c}| }f(| {{\bf{x}}}^{p}-{{\bf{x}}}_{c}| ).\end{eqnarray} \tag{ 9 }$$
  The first two cases each operate in only a single cartesian direction. In the third (radial) case, the center ${{\bf{x}}}_{c}$ of the distortion is specified. There is a single free parameter, the scaling parameter s. The template function f is defined as linear interpolation between values v_j at nodes ${a}_{0}+j\,{\rm{\Delta }}a$ for $0\leqslant j\leqslant N$.
• Piecewise maps are functionally identical to the Template map, except that the nodal values v_j are the free parameters, and the scaling is fixed to s = 1.
• A Color term is defined by
  $$\begin{eqnarray}&&{{\bf{x}}}^{w}={{\bf{x}}}^{p}+(c-{c}_{\mathrm{ref}})[{\bf{m}}({{\bf{x}}}^{p})-{{\bf{x}}}^{p}],\end{eqnarray} \tag{ 10 }$$
  where ${c}_{\mathrm{ref}}$ is a reference color and ${\bf{m}}$ is an instance of any of the above forms of transformation. The parameters of the Color map are those of the map it scales.
• Reprojection maps have no free parameters: they merely move coordinates from one projection of the sphere to another.
• Composite maps realize Equations (4)–(6) for a specified sequence of any of maps (including other Composite maps). The parameters of the composite are the concatenation of those of the component maps.

A PixelMap, in combination with a specification of the projection in which the ${{\bf{x}}}^{w}$ maps to the celestial sphere, forms a complete world coordinate system (WCS).

2.3. Degeneracies

When minimizing ${\chi }^{2}$, we must be aware of degeneracies whereby ${\boldsymbol{\pi }}$ can change while ${\chi }^{2}$ is invariant. Such degeneracies will lead to (near-)zero singular values in the normal matrix ${\bf{A}}$ used in the solution for ${\boldsymbol{\pi }}$ (Section 4), and failures or inaccuracies in its inversion. There are several such landmines that we must clearly avoid. We assume in this discussion that the astrometric model for each extension is a device-based instrumental function ${{\bf{x}}}^{T}=D({{\bf{x}}}^{p})$ from pixel to "telescope" coordinates, followed with an exposure-based function ${{\bf{x}}}^{w}=E({{\bf{x}}}^{T}).$

2.3.1. Shift

The simplest degeneracy is a shift in all stellar positions, $E\to E+{\rm{\Delta }}{\bf{x}}$ for every exposure (in the flat-sky limit; more generally, the degeneracy is a rotation of the celestial sphere). Each star α has its derived sky position ${\bar{{\bf{x}}}}_{\alpha }^{w}$ shifted as well, but because ${\chi }^{2}$ is differential, there is no effect on ${\chi }^{2}.$ This degeneracy is broken by having a reference catalog for which ${{\bf{x}}}^{w}$ is fixed. The reference catalogs do not need to be very dense or precise to break this degeneracy.

2.3.2. Color Shift

A color-dependent shift $E\to E+c\,{\rm{\Delta }}{\bf{x}}$ is also undetectable in the differential ${\chi }^{2}$. This degeneracy is broken if colors are known for reference stars over a finite range of color.

2.3.3. Linear

In the flat-sky limit, consider the case where the exposure component for exposure k is an affine transformation ${E}_{k}\,={{\bf{A}}}_{k}{{\bf{x}}}^{T}+{{\bf{x}}}_{k}$ with linear rescaling ${{\bf{A}}}_{k}$ and offset ${{\bf{x}}}_{k}$, the latter corresponding to the pointing of the telescope at exposure k. An object α with world coordinates ${{\bf{x}}}_{\alpha }^{w}$ will be observed at telescope coordinate ${{\bf{x}}}_{k\alpha }^{T}={{\bf{A}}}^{-1}({{\bf{x}}}_{\alpha }^{w}-{{\bf{x}}}_{k}).$ For any non-degenerate matrix ${\bf{B}}$, there is an alternative solution:

$$\begin{eqnarray}&&{{\bf{A}}}_{k}\to {{\bf{BA}}}_{k},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{\bf{x}}}_{k}\to {{\bf{B}}}^{-1}{{\bf{x}}}_{k},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{{\bf{x}}}_{\alpha }^{w}\to {{\bf{Bx}}}_{\alpha }^{w},\end{eqnarray} \tag{ 13 }$$

which leaves ${\chi }^{2}$ unchanged. The shift degeneracy is a special case of affine degeneracy. This degeneracy is also broken by the existence of a sample of reference stars. There is also a color-dependent variant of this degeneracy.

If the exposure solution E_k has freedom to be altered by some global polynomial function B of order n, then there is generalization of this degeneracy in which each E_k is shifted by a polynomial of order $n-1$. Again, the solution is to have a reference catalog of even modest density and accuracy.

2.3.4. Colony Collapse Disorder

WcsFit is accelerated by calculating the weight of each observation in Equation (3) just once at the start of fitting, using the determinant of the starting WCS system to convert the pixel errors into world-coordinate errors. This opens the door to a pseudo-degeneracy in which all output ${{\bf{x}}}^{w}$ values are scaled by some matrix ${\bf{B}},$ sending ${\chi }^{2}\to | {\bf{B}}| {\chi }^{2}$. If $| B| \to 0$, the solution appears to approach perfection while collapsing the output map. This is countered by an increase in ${\chi }^{2}$ contributed by the reference stars, which are not collapsing; however, if the total weight of the reference stars is too low, the solution will tend toward collapse. The collapse becomes complete if the reference stars are then flagged as outliers and removed by our σ-clipping step. WcsFit includes a parameter to scale the weights of the reference catalogs, which can be used to prevent this collapse solution if the reference catalog is sparse. Over-weighting the reference catalog could lead to ill-conditioned solutions, but we are not in this regime.

2.3.5. Exposure/Instrument Trades

For any map F, the transformations

$$\begin{eqnarray}&&{E}_{k}\to {E}_{k}F\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&D\to {F}^{-1}D\end{eqnarray} \tag{ 15 }$$

clearly leave ${\chi }^{2}$ and all ${{\bf{x}}}^{w}$ values invariant. If the functional forms being used for E and D admit such a transformation, then the solution is degenerate. The WcsFit code searches for cases where multiple Constant, Linear, or Polynomial atomic map elements are composited into any exposures' pixel maps and are hence able to trade their terms. This degeneracy can be broken by setting one of the exposure maps E_k to the Identity map. WcsFit will do this automatically if the user's configuration leaves such degeneracies in place. Note that this means the "instrumental" solution ends up including the low-order transient components of one exposure. This is immaterial if the instrumental solution is always applied in series with a term with these low-order degrees of freedom adjusted to fit each exposure's data.

2.3.6. Unconstrained Parameters

Map parameters are, of course, degenerate if there are no stellar observations being affected by them (e.g., if a given exposure did not generate any matched detections or they have all been removed as outliers, then the parameters of the exposure solution are unconstrained).

WcsFit checks the normal matrix ${\bf{A}}$ for null rows that arise when a parameter does not act on any observations. In this case, the diagonal element on this row is set to unity, which stabilizes the matrix inversion and freezes this (irrelevant) parameter in further iterations.

More troublesome is the case where there are a small but non-zero number of observations on an exposure, too few to constrain the model, so that ${\bf{A}}$ is degenerate but without null rows. In this case, WcsFit will fail the attempt to do a Cholesky decomposition of the non-positive-definite ${\bf{A}}$. In this case, WcsFit will perform a singular value decomposition of ${\bf{A}}$, report to the user the eigenvectors associated with near-null singular values, and then quit. It is then up to the user to determine the nature of the degeneracy.

2.4. Algorithms

The WcsFit software suite assumes that we are already in possession of an initial WCS for each extension of sufficient accuracy to allow unambiguous matching of common detections of a source. scamp is routinely run on each DES exposure to generate this starting WCS, with accuracy of $\lt 1^{\prime\prime}$ relative to Gaia or other reference catalog.

2.4.1.  ${\chi }^{2}$ Minimization

Another benefit of having a good starting WCS for each exposure is that we can initialize parameters of the maps that are defaulted on input by fitting them to the starting WCS. WcsFit does this by (1) placing test points on a grid of ${{\bf{x}}}^{p}$ spanning the device and treating these as a catalog of detections; (2) using the starting WCS to map these to positions ${{\bf{x}}}_{0}^{w};$ and (3) finding the parameters ${\boldsymbol{\pi }}$ of our model that minimize the mean $\langle {[{{\bf{x}}}_{0}^{w}-{{\bf{x}}}^{w}({{\bf{x}}}^{p},{\boldsymbol{\pi }})]}^{2}\rangle$ over the grid of points.

The algorithm for minimization of ${\chi }^{2}$ assumes that the minimizing solution is close to the starting solution, i.e., we are doing fine tuning after scamp has done the work of bringing us close. The positions are close to linear in the parameters, so the ${\chi }^{2}$ value should be close to the usual quadratic form

$$\begin{eqnarray}&&{\chi }^{2}\approx {\chi }^{2}({{\boldsymbol{\pi }}}_{0})+2{\bf{b}}\cdot {\rm{\Delta }}{\boldsymbol{\pi }}+{\rm{\Delta }}{\boldsymbol{\pi }}\cdot {\bf{A}}\cdot {\rm{\Delta }}{\boldsymbol{\pi }},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}{b}_{\mu } & \equiv & \displaystyle \frac{1}{2}\displaystyle \frac{\partial {\chi }^{2}}{\partial {\pi }_{\mu }}=\displaystyle \sum _{i}{w}_{i}({{\bf{x}}}^{w}({{\bf{x}}}_{i}^{p},{{\boldsymbol{\pi }}}_{0})-{\bar{{\bf{x}}}}_{{\alpha }_{i}})\\ & & \cdot \left(\displaystyle \frac{\partial {{\bf{x}}}^{w}({{\bf{x}}}_{i}^{p},{\boldsymbol{\pi }})}{\partial {\pi }_{\mu }}-\displaystyle \frac{\partial {\bar{{\bf{x}}}}_{{\alpha }_{i}}}{\partial {\pi }_{\mu }}\right)\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{A}_{\mu \nu }\equiv \left(\displaystyle \frac{\partial {{\bf{x}}}^{w}({{\bf{x}}}_{i}^{p},{\boldsymbol{\pi }})}{\partial {\pi }_{\mu }}-\displaystyle \frac{\partial {\bar{{\bf{x}}}}_{{\alpha }_{i}}}{\partial {\pi }_{\mu }}\right)\cdot \left(\displaystyle \frac{\partial {{\bf{x}}}^{w}({{\bf{x}}}_{i}^{p},{\boldsymbol{\pi }})}{\partial {\pi }_{\nu }}-\displaystyle \frac{\partial {\bar{{\bf{x}}}}_{{\alpha }_{i}}}{\partial {\pi }_{\nu }}\right).\end{eqnarray} \tag{ 18 }$$

Note that the weights w_i are being assumed independent of ${\boldsymbol{\pi }}$, i.e., the world-coordinate errors ${\sigma }_{i}^{w}$ of each exposure are held fixed at the values implied by the starting WCS. Also note that WcsFit does not treat the true positions ${{\bf{x}}}_{\alpha }^{w}$ of the sources as free parameters. Instead, the dependence of the mean of the measurements ${\bar{{\bf{x}}}}_{\alpha }$ upon the parameters is propagated directly into the normal equation.

The calculation of ${\bf{b}}$ and ${\bf{A}}$ is the most computationally intensive part of WcsFit. The summation for matches is distributed across cores using OpenMP calls. Each match is dependent upon the limited subset of the parameters ${\boldsymbol{\pi }}$ which appear in the pixel maps for the extensions in which the object is observed, hence the updates to ${\bf{A}}$ are sparse, though the final matrix is dense.

WcsFit first attempts the Newton iteration

$$\begin{eqnarray}&&{\boldsymbol{\pi }}\to {\boldsymbol{\pi }}-{{\bf{A}}}^{-1}{\bf{b}}.\end{eqnarray} \tag{ 19 }$$

The solution time scales as the cube of the number of free parameters, and is executed using a multi-threaded Cholesky decomposition after preconditioning ${\bf{A}}$ to have unit diagonal elements. Despite the cubic scaling, this step is usually faster than the calculation of the normal matrix. If the decomposition fails due to a non-positive-definite ${\bf{A}}$, WcsFit performs a singular value decomposition on ${\bf{A}}$ and informs the user which parameters dominate the degenerate vectors.

The Newton step is iterated until ${\chi }^{2}$ no longer decreases by more than a chosen fraction. Should ${\chi }^{2}$ increase during an iteration or fail to converge within a selected number of steps, then the minimization process is re-started using a Levenberg-Marquart algorithm based on the implementation by Press et al. (2003).

2.4.2. Outlier Rejection

The WcsFit solutions must be robust to astrometric measurements perturbed by unrecognized cosmic rays or defects on the stellar images and by stars with proper motion or binary partners that alter the photocenter by amounts exceeding measurement errors. We do not at this time fit for proper motion or parallax within WcsFit.

Outlier rejection is done using standard σ-clipping algorithms. A clipping threshold, t, is specified at input. After each ${\chi }^{2}$ minimization, a rejection threshold is set at $t\sqrt{{\chi }^{2}/\mathrm{DOF}}$. For N detections grouped into M distinct matches and P free parameters in ${\boldsymbol{\pi }}$, we have ${DOF}=2\ast (N-M)-P.$ Detections whose residual to the fit (in units of σ) exceeds the threshold are discarded. At most one outlier per match is discarded at each clipping iteration.

Outlier clipping is alternated with ${\chi }^{2}$ minimization until the clipping step no longer reduces the ${\chi }^{2}$ per degree of freedom by a significant amount.

2.4.3. Procedure

The steps in the astrometric solution process are as follows:

• 1.  
  A preparatory Python program reads an input YAML configuration file specifying the desired input catalog files, plus the definitions of the fields, epochs, and instruments. It then collects from all the catalogs and their headers any information necessary to construct tables of extensions, devices, exposures, and instruments. This includes extracting the serialized starting WCS, usually as produced by scamp and stored in the headers of the FITS catalog extensions.
• 2.  
  A second preparatory program reads all the detections from the input catalogs, applying any desired cuts for S/N and stellarity, and then runs a standard friends-of-friends algorithm to identify all matching detections. Any match that includes multiple detections from the same exposure is discarded. The IDs of all groups of matching detections are then stored in another FITS table.
• 3.  
  WcsFit starts by ingesting the input FITS tables and creating the structures defining instruments, devices, exposures, and extensions.
• 4.  
  The YAML file specifying the pixel maps to be applied to each extension is parsed, and a PixelMap is created with specified or defaulted parameters. Any of the map elements may have its parameters frozen by the user, and the remainder are the free parameters of our model.
• 5.  
  WcsFit checks the map configuration for degeneracies: is there reference catalog in each field? Are there are any exposure/instrument degeneracies? If so, WcsFit will attempt to break the degeneracies by setting one or more exposures' maps to Identity.
• 6.  
  All exposures in a field are reprojected to a common gnomonic system about the field center.
• 7.  
  Any parameters of PixelMaps that were set to defaults have their values set by a least-squares fit to the starting WCS. Any degeneracies halt the program.
• 8.  
  The ${{\bf{x}}}_{i}^{p}$ and ${\sigma }_{i}^{p}$ of all detections that are part of useful matches are extracted from their source catalogs. For any detections whose maps include color terms, we require a measurement from a color catalog to be matched to the same object. The color catalog is read at this point.
• 9.  
  A requested fraction of the matches are excluded from the fit at random. These reserved matches can be used later to validate the fit.
• 10.  
  Any exposures containing insufficient detections are removed from the fit.
• 11.  
  The iteration between ${\chi }^{2}$ minimization and σ-clipping begins. At each iteration, ${\bf{A}}$ is checked for null rows as noted in Section 2.3, which are altered so as to freeze the associated parameter. If ${\bf{A}}$ is not positive-definite, WcsFit reports the nature of the associated degenerate parameters, then exits.
• 12.  
  The best-fit astrometric model is written to an output YAML file.
• 13.  
  After completion of the fit, the best-fit map is applied to both the fit and reserved matches. The σ-clipping algorithm is applied iteratively to the reserved matches.
• 14.  
  The rms residual and ${\chi }^{2}$ statistics are reported for the unclipped detections on each exposure.
• 15.  
  The input, output, and best-fit residual for every detection are written to an output FITS table for further offline analyses.

2.5. Performance

The astrometric solution is derived from multiple epochs of the star-flat observations described above. Table 1 lists the dates and conditions of the star-flat sequences during the first four years of DECam operations for which there were neither clouds nor instrument anomalies.³⁵ Exposures are usually 30 s long, with 25–30 s dead time for readout and repointing, so the star-flat sequence for five filters consumes about 2 hours of clock time. Figure 1 shows a typical star-flat pointing sequence of 22 exposures. There are ≈10⁵ Gaia reference stars available in the area covered by the star-flat exposures of a given field.

Zoom In Zoom Out Reset image size

Each exposure is run through the standard DES data reduction pipeline, including linearization of images; crosstalk removal; correction for the "brighter-fatter effect" (Gruen et al. 2015); debiasing and division by dome flats; and subtraction of sky and fringe signals. Sources are detected and measured using SExtractor (Bertin & Arnouts 1996). For the following analyses, we filter the catalogs for sources with no SExtractor flags set, no defective, saturated, or cosmic-ray-flagged pixels within the isophote, with MAGERR_AUTO < 0.01, indicating signal-to-noise ratio $S/N\gtrsim 100$, and with $| {\mathtt{SPREAD}}\_{\mathtt{MODEL}}| \,\lt {\mathtt{0}}.{\mathtt{003}}$ to select only stellar sources. The flag cut removes objects that overlap detected neighbors.

The windowed centroids $({\mathtt{XWIN}}\_{\mathtt{IMAGE}},{\mathtt{YWIN}}\_{\mathtt{IMAGE}})$ are used for centroid positions, ${{\bf{x}}}^{p}$, as they have been demonstrated to be robust to the details of the point-spread function (PSF) while approaching the accuracy of ideal PSF-fitting astrometry. Our focus on astrometric errors that correlate over space and/or time means we will not investigate the vagaries of centroid measurement, e.g., pixel-phase errors.

The density of useful stellar positions varies with field, filter, seeing, and sky conditions, but is usually 200–400 per CCD, or more than 10,000 per exposure and $\gt {10}^{6}$ per star-flat epoch.

Our goal is to produce an astrometric model that maps the ${{\bf{x}}}^{p}$ of a source to ICRS sky coordinates, such that any coherent errors are at ≲1 mas rms. Coherence applies here to both time and space, meaning that the error should persist across more than one star and more than one exposure. Note that we are not attempting to model the following effects:

• Shifts in the centroids of individual detector pixels due to variation in lithography of the gate structures. There is not enough on-sky stellar data to calibrate this for the 500 megapixels in DECam, but astrometric errors due to pixel-to-pixel variations will behave as noise in individual stars' positions, and will not correlate between stars. In a well-designed use of DECam, a given star will be exposed on different parts of the array in each exposure, and hence this error will not correlate across time either. In any case, the rms variation in DECam pixel sizes is estimated (from flat-field behavior) to be at a few parts per thousand, or <1 mas. Stellar position errors will be even lower because they average over a PSF containing O(10) pixels.
• Stochastic atmospheric distortions on $\lt 1^\circ$ scale. Such distortions are not coherent between exposures, but they are ≫1 mas and dominate the astrometric error budget for high-S/N detections. The characteristics of these fluctuations are investigated in Section 6.
• Other sub-mas effects.

The DECam astrometric model was constructed through careful examination of the residual astrometric errors in the star-flat data. The final choice of model is given in Table 2. Here we describe each element of the model in more detail, tracing backwards from the collected charge in the pixel well back to the top of the atmosphere.

Table 2.  Components of the DECam Astrometric Model

Description	Name	Type	Max. Size
Tree-ring distortion	$\langle {band}\rangle /\langle {device}\rangle /{\mathtt{rings}}$	Template (radial)	$\approx 0\buildrel{\prime\prime}\over{.} 05$
Serial edge distortion	$\langle {band}\rangle /\langle {device}\rangle /{\mathtt{lowedge}}$	Template (X)	$0\buildrel{\prime\prime}\over{.} 03$
Serial edge distortion	$\langle {band}\rangle /\langle {device}\rangle /{\mathtt{highedge}}$	Template (X)	$0\buildrel{\prime\prime}\over{.} 03$
Optics	$\langle {band}\rangle /\langle {device}\rangle /{\mathtt{poly}}$	Polynomial (order = 4)	$\gg 1^{\prime\prime}$
Lateral colora	$\langle {band}\rangle /\langle {device}\rangle /{\mathtt{color}}$	Color × Linear	$\approx 0\buildrel{\prime\prime}\over{.} 04$
CCD shift	$\langle {epoch}\rangle /\langle {device}\rangle /{\mathtt{ccdshift}}$	Linear	$\approx 0\buildrel{\prime\prime}\over{.} 1$
Exposure	$\langle {exposure}\rangle$	Linear	$\gg 1^{\prime\prime}$
Differential chromatic refraction	$\langle {exposure}\rangle /{\mathtt{dcr}}$	Color × Constant	$\approx 0\buildrel{\prime\prime}\over{.} 05$

Note.

^aThe lateral color correction is set to Identity transformation for izY bands.

Download table as:  ASCIITypeset image

4.1. Tree Rings

In $g,r,$ and i bands, photons generate holes near the DECam CCD surface and then have to drift the 250 μm thickness of the device before being collected in the pixels. As described in Plazas et al. (2014), any electric field components transverse to the surface will cause the charge carriers to drift sideways before collection and induce an apparent astrometric shift. The DECam CCDs are known to have two significant sources of such stray fields. The first are "tree rings," which arise from fluctuations in the impurity density of the silicon boules from which the CCD wafers were cut. The zone refining of the boules results in approximate circular symmetry about the boule axis and the wafers are cut perpendicular to this axis, so the astrometric distortions are realized as an irregularly oscillating pattern of rings. For some DECam devices, the ring centers are on the device, for others the centers are off their edges. Because the distortions also produce oscillations in the solid angle of sky received by each pixel, they are readily apparent in the flat-field images. The nearly circular symmetric pattern in the flat fields implies that the astrometric distortions share this symmetry and are directed radially toward (or away from) the ring center. As described in Plazas et al. (2014), we locate the ring center for each CCD by visual inspection of the flat-field images, and then create templates of the expected astrometric distortion about this center from a high-pass-filtered, azimuthally averaged profile of the flat-field signal. Figure 2 plots the template derived for a representative device.

Zoom In Zoom Out Reset image size

In WcsFit, the tree-ring signal is realized as a Template map, with both the variation and the displacement expected to be purely radial to the rings. We have a single free parameter for each device/filter combination, which is a multiplicative scaling of the distortion predicted by the template. We do not allow for any time variation of the tree-ring signal because the effect is literally built into the device. We do allow for a band dependence, however, as photons in the z and Y bands penetrate well into the device and are therefore expected to suffer less deflection before collection, on average. Figure 3 plots the best-fit template coefficients for all devices and filters. We do not know why the tree-ring distortions are seen to be only 80%–90% of the values predicted from the r-band flat-field images, but the scaling of these coefficients with filter band hews closely to the values calculated from the absorption versus wavelength characteristics of silicon. Figure 2 plots the azimuthally averaged residual position for all detections from the gri exposures of a representative device. The rms of this residual is at our goal level of ≈1 mas.

Zoom In Zoom Out Reset image size

4.2. Edges

The electric field in the CCD also develops a substantial transverse component near the device edges. The subsequent astrometric distortion and pixel size variation is readily apparent in the flat-field images as a "glowing edge." It is found that the flat-field (photometric) edge behavior is not a good predictor of the astrometric distortions, so we derive a template for edge behavior entirely from the stellar astrometry. We assume throughout that the edge distortion is directed in the x direction (parallel to the serial register on the short edge of the device) and is constant along y at each edge. We first fit the star-flat data to a model with a Piecewise displacement term with a free node position every eight pixels within 180 pixels of each edge.

Note that the 25 (15) pixels of the device nearest to the long (short) edges are completely masked from analysis because the distortion is too large. Thus, we do not have useful stellar centroids closer than ≈30 pixels to the $x=1,x=2048$ boundaries. Any nodal values in these regions are unconstrained and ignored. There are also unusable nodal values near the locations of any defective columns on a device.

Upon examination of the best-fit piecewise solutions at the x edges, we find that all edges of all CCDs in all filters are consistent with a common "master" edge template once we allow for a multiplicative scaling and a shift as large as 12 pixels (0.18 mm). These shifts might arise from the finite precision of the cutter tooling relative to the array during CCD dicing. The master edge template is shown in Figure 4. In the final astrometric fit, we allow each device/filter combination to have a Template pixel map at the high- and low-x edges. The templates are shifted versions of the master template, and the scaling is left as a parameter for WcsFit to optimize.

Zoom In Zoom Out Reset image size

Figure 5 plots the binned displacement residuals in all filters near all four edges of the CCD after the x-edge template is included in the WcsFit model. The master template reduces rms x residuals to well below 1 mas. Note that we have elected to make no correction at all for the glowing edge effect on the short (y) edges because the displacement is already <3 mas before correction. Because it affects only a small fraction of the focal plane, the rms error is ≪1 mas.

Zoom In Zoom Out Reset image size

The best-fit coefficients to the master template are found to be in the range 0.8–1.2 in the gri bands, with lower values in z and Y as expected again from the deeper photon conversion. We take the edge coefficients to be independent of time, as one would expect for such detector-physics effects.

4.3. Optics Polynomials

4.4. CCD Shifts

4.5. Lateral Color

Any axisymmetric refractive optical system is expected to have color-dependent radial distortion, leading to color terms described by odd-order polynomials in radius. We check this assumption by including in our initial fits a more general color term: a time-independent linear function of coordinates on each CCD. The displacement is assumed to be proportional to

$$\begin{eqnarray}&&c\equiv {(g-i)-(g-i)}_{\mathrm{ref}},\end{eqnarray} \tag{ 20 }$$

where the reference color is chosen to be $0.44,$ the color in the natural DECam $(g-i)$ system of the F8IV star C26202 from the Hubble Space Telescope CalSpec system.³⁶ We restrict the fit to stars with $-0.2\leqslant (g-i)\leqslant 1.8,$ to avoid M stars for which the expected shifts may no longer be linear in c, and assume that the color term is time-independent for a given filter.

Figure 6 plots the best-fit static color solution in the g and r bands, which show the radial patterns (at mas accuracy) and approximate amplitude expected from the optical solution (S. Kent, private communication). The solutions for $i,z,$ and Y bands are, as expected, undetectably small as the chromatic terms of the corrector lenses are weaker, and we disable their color terms for the final WcsFit run. For the g and r bands, we continue by re-fitting to a more restrictive function, namely a Color PixelMap wrapping a radial Piecewise displacement. A fifth-order function of field radius is found to fit the resultant piecewise function with 2 mas/mag rms scatter (Figure 6). We used the fitted polynomials (implemented as radial Template pixel maps) for our final model, with no free parameters.

Zoom In Zoom Out Reset image size

4.6. Exposure Solution

4.7. Differential Chromatic Refraction

The color dependence of atmospheric refraction in the context of wide-field cosmological surveys is studied by Plazas & Bernstein (2012) and Meyers & Burchat (2015), who conclude that in some bands it will be present at O(10) mas and significant for weak gravitational lensing analyses. The atmospheric refraction is very large ($1^{\prime} \times \tan z$ at zenith angle z) compared to our desired accuracy so the chromatic effect is significant. For a given site, it is expected to behave as

$$\begin{eqnarray}&&{\rm{\Delta }}{{\bf{x}}}^{w}={K}_{b}c\tan z\,\hat{{\bf{p}}},\end{eqnarray} \tag{ 21 }$$

where c is the object color (again g − i), $\hat{{\bf{p}}}$ is the unit sky vector toward the zenith (the parallactic direction), and K_b is a constant derivable for each band b from the instrumental bandpass and atmospheric index of refraction.

We test this model by allowing each exposure to have its own constant differential chromatic refraction (DCR) term ${\rm{\Delta }}{{\bf{x}}}^{w}$ when fitting to the star-flat detections. Figure 7 shows these results for g and r bands, along with the model (21) with the best-fit value of K_b. The standard atmospheric model is seen to describe the measured $c\,{\rm{\Delta }}{{\bf{x}}}^{w}$ well, with rms residuals of 2–3 mas/mag. Table 3 gives the DCR amplitudes K_b derived from the star-flat data (supplemented with the supernova-field data described in Section 7 for the i and z bands). The K_b are in good agreement with the predictions of Plazas & Bernstein (2012) and for future use it should suffice to simply fix the DCR term to Equation (21) instead of allowing freedom to each exposure.

Zoom In Zoom Out Reset image size

Table 3.  Differential Chromatic Refraction for DECam

Band	Kb (mas/mag)
g	45.0
r	8.4
i	3.2
z	1.4
Y	1.1

Download table as:  ASCIITypeset image

The first panel of Figure 8 plots the errors ${\rm{\Delta }}{\bf{x}}$ in the stellar positions of a randomly selected single exposure on a randomly selected CCD, relative to the mean positions determined for the same stars from the entire stack of star-flat exposures. The residuals after application of our final astrometric model are dominated by a coherent pattern with rms $| {\rm{\Delta }}{\bf{x}}| \approx 15$ mas. This pattern is found to differ from exposure to exposure and is plausibly attributed to refraction by atmospheric turbulence. Before investigating these stochastic distortions in Section 6, we ask here whether there are any distortion patterns that recur from exposure to exposure. To find them, we will need to average down the stochastic atmospheric signal by stacking and binning the residuals from many exposures.

Zoom In Zoom Out Reset image size

The middle panels of Figure 8 show the result of averaging the residuals from all of the $g,r,$ and i band exposures from a single star-flat epoch and then from all epochs. The amplitude of the residual pattern drops steadily with number of exposures included, although not as quickly as the square root of the number of exposures, as would be expected if all remaining errors were uncorrelated between exposures. The all-gri residuals for this CCD exhibit clear coherence of several kinds:

• A radial pattern just above center contains distortions with amplitude of up to 30 mas. This roughly coincides with the edge of the electrical connector that is soldered and glued to the CCD and protrudes through a hole in the CCD mounting board. The right-hand panel of Figure 8 averages over all 30 functional DECam science CCDs mounted in the same orientation, showing that this pattern is recurrent and exists at both ends of the connector mount and is likely a product of stresses induced in the CCD lattice the connector or the hole in the mounting board.
• There are excess residuals along the long edges of the device, suggesting that the edge distortion is not uniform along the edge of the device. The right-most panel confirms that this is true in a systemic fashion for the devices.
• The right-most panel shows the largest residuals in the corners of devices, plus two patches at the midpoints of the long edges. These six locations are known as "tape bumps" because they are underlaid by thin spacers that define the thickness of the glue layer between the CCD and its carrier. These regions, each 100–200 pixels on a side, exhibit structure in the flat fields that is indicative of stray transverse electric fields induced by lattice stresses. There are clear astrometric disturbances associated with these fields as well. Because these are difficult to model and cover only a small fraction of the focal plan, we do not attempt to remove them: detections occurring on the tape bumps are flagged as having less reliable astrometry, and in fact have been omitted from the characterization and modeling performed in this paper.

The DECam flat fields show evidence for modulations of the pixel size with period ${\rm{\Delta }}x=27.33$ pixels, a behavior seen in many CCDs due to the step-and-repeat accuracy of the mask generator (Anderson & King 1999). A corresponding periodicity is detectable in the astrometric deviations, but with peak-to-peak amplitude <1 mas, as is predicted from the amplitude of the flat-field fluctuations. We ignore this effect for DECam.

A final question we address about the residuals to the model is whether the linear-per-ccd ccdshift terms are sufficient to describe the change of the astrometric solution between epochs. Figure 9 shows a test of this, whereby we plot the mean residual astrometric errors for an entire epoch's gri detections, averaged by position on the array. The residuals are consistent with the expectations of averaging 66 realizations of the stochastic atmospheric pattern. In particular there are no statistically significant discontinuities across CCD boundaries.

Zoom In Zoom Out Reset image size

We conclude that our astrometric model captures the recurrent instrumental distortion pattern to an accuracy of 2–4 mas rms. The only residual distortions that have detectable coherence between the 64-pixel (16'') bins of Figure 8 are those associated with stresses from the CCD mounting, and with inhomogeneities in the transverse electric fields generated at the CCD edges.

After application of the astrometric model, the dominant form of astrometric error is a field that varies from exposure to exposure and has a coherence length of 5–10'. Figure 10 plots the residual vector field for a representative exposure.

Zoom In Zoom Out Reset image size

6.1. Atmospheric Turbulence

Multiple lines of evidence support the hypothesis that these distortions arise from refraction by atmospheric turbulence:

• The patterns are uncorrelated between exposures, thus change on timescales of 1 minute or less. The only physical conditions that should change this quickly are the atmosphere and the settings of the hexapod that fixes the alignment of the camera to the primary.
• The distortion pattern appears to be curl free. The lower panel of Figure 10 suggests that the curl arises from white noise, i.e., errors in stellar positions due to shot noise. This is shown more rigorously in Figure 11, in which the two-point correlation function of astrometric errors is split into E- and B-mode components (curl- and divergence-free, respectively), as explained in Appendix. The latter is seen to be consistent with zero. Curl-free distortion patterns are expected in the ray-optic limit, where the astrometric displacement of each photon is the gradient of the integral of the scalar index of refraction (time delay) along the line of sight to the star.
• The distortion pattern is clearly anisotropic, with a long correlation length in one direction. The preferred direction is roughly, but not exactly, consistent between exposures, as is expected from having the atmospheric turbulence pattern blowing across the FOV during the exposure. The streaky patterns are very similar to the maps of PSF ellipticity in short exposures with the CFHT Megacam presented by Heymans et al. (2012), which they also attribute to wind.
• The power spectrum of the distortion in the cross-wind direction is roughly consistent with that expected of Kolmogorov turbulence.
• The amplitude and correlation length of the distortion are roughly consistent with numerical simulations of Kolmogorov turbulence (J. Peterson, private communication). The simulations suggest that the astrometric perturbations are strongly dependent on the outer scale of the turbulence.

If the stochastic distortions do indeed arise from atmospheric turbulence, we expect their amplitude to decrease with the square root of exposure time as we average over phase screens. We cannot verify this with our data because nearly all star-flat exposures were taken with 30 s exposures, save the first three epochs which used 50 s. While these early epochs do show the lowest stochastic distortion (see Figure 12), there is substantial variation from epoch to epoch so we cannot draw any quantitative conclusions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.2. Behavior of the Stochastic Component

For a more quantitative picture of the stochastic/atmospheric distortion field, we produce its two-point correlation function,

$$\begin{eqnarray}&&{\xi }_{+}(r)\equiv {\langle {\rm{\Delta }}{{\bf{x}}}_{i}\cdot {\rm{\Delta }}{{\bf{x}}}_{j}\rangle }_{| {{\bf{x}}}_{i}-{{\bf{x}}}_{j}| =r},\end{eqnarray} \tag{ 22 }$$

where ${{\bf{x}}}_{i}$ is the sky position of detection i, ${\rm{\Delta }}{{\bf{x}}}_{i}$ is the measurement error on this position, and the average is taken over all pairs of detections $i\ne j$ in the same exposure separated by distance r (in practice the "true" position is taken as the average of our many measured positions). The Appendix explains how ξ can be split into two components ${\xi }_{E}$ and ${\xi }_{B}$ which arise from the curl- and divergence-free parts of the vector field, respectively. These are plotted in Figure 11, where it is clear that ${\xi }_{B}$ is insignificant in comparison to ${\xi }_{E},$ as expected for atmospheric refraction. In this case, and if we consider the turbulence to be a Gaussian random field, then ${\xi }_{E}={\xi }_{+}$ fully characterizes the field.³⁷

Before calculating ξ, we subtract from the ${\rm{\Delta }}{{\bf{x}}}_{i}$ the best-fit cubic polynomial function of field coordinates. As discussed earlier, we expect significant linear and quadratic-dependence distortions from normal (homogeneous) atmospheric refraction; turbulent refraction should also have a substantial large-scale component, and indeed we observe $\approx 25 \%$ of the distortion variance to come from this polynomial signal. Because the low-order component is easily determined in practice by fitting to the Gaia catalog, we remove it from our analysis, leaving small-scale distortions. Note that the virtue of using $\xi (r)$ is that it is unaffected by the shot-noise measurement errors of the stellar positions (for $r\gt 0$) and therefore is a pure measure of the astrometric map.

We characterize the astrometric correlation function by two quantities, ${\xi }_{0}$ and ${R}_{\mathrm{corr}}.$ The former is essentially the small-scale (largest) value of ${\xi }_{+},$ defined as

$$\begin{eqnarray}&&{\xi }_{0}\equiv \displaystyle \frac{{\int }_{0}^{\infty }r\,{dr}\,W(r){\xi }_{+}(r)}{{\int }_{0}^{\infty }r\,{dr}\,W(r)}\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&W(r)\equiv {e}^{-{r}^{2}/2{\sigma }^{2}},\qquad \sigma =1^{\prime} .\end{eqnarray} \tag{ 24 }$$

The second quantity ${R}_{\mathrm{corr}}$ is defined as the radius at which ${\xi }_{+}({R}_{\mathrm{corr}})=0.5{\xi }_{0},$ i.e., the smoothing scale that would cut the astrometric variance in half.

Figure 12 plots the values of ${\sqrt{\xi }}_{0}$ and ${R}_{\mathrm{corr}}$ for all the star-flat exposures under analysis, along with the half-light diameter D₅₀ of the PSF in each exposure. It is clear that there are nights when a degradation of seeing is accompanied by an increase in astrometric distortion (e.g., 20140118, 20150926), as one might expect if both are proportional to the amplitude of a strictly Kolmogorov turbulence spectrum. However, there are also cases of anti-correlation, and the mean seeing of a night is a very weak predictor of astrometric accuracy. In particular, the epoch 20140807 is astrometrically awful, exhibiting 30–50 mas rms atmospheric contribution whereas most other epochs are 10–20 mas rms. Yet the seeing on that night was not as poor. Clearly, there are other variables besides Kolmogorov amplitude, such as wind speed or outer scale, that determine the astrometric quality of the night.

The correlation angle is in the range $4^{\prime} \lt {R}_{\mathrm{corr}}\lt 10^{\prime}$ at all epochs, with no apparent relation to the seeing. This suggests that interpolation between stars in Gaia catalog, with source density of $\approx 1$ star per arcmin² at high latitude, could significantly reduce the stochastic atmospheric errors. We investigate this in Section 8.

The DECam measures of stochastic atmospheric astrometric fluctuations are in rough agreement with previous characterizations by Han & Gatewood (1995); Zacharias (1996), and Bouy et al. (2013) (and references therein), in terms of typical amplitude at good sites, and substantial night-to-night variation that is at best tenuously tied to the seeing full width at half maximum (FWHM).

Is the DECam astrometric solution stable over weeks, days, or years? We already know that there is short-term (seconds) variability due to atmospheric turbulence, at a typical level of 10–20 mas rms in a 30 s exposure and 5–10' correlation length, but this should average to zero with longer exposures. We have verified in Section 5 that the astrometric errors within a given star-flat epoch (i.e., a few hours' clock time) are consistent with a single solution, up to the accuracy allowed by the stochastic atmospheric distortions, as long as we allow for expected exposure-to-exposure variations at low order across the focal plane. We are interested in the duration over which a single astrometric solution can otherwise be considered to maintain few-mas accuracy.

Our fit to all of the star-flat data allowed for variations between star-flat epochs in the form of a free linear transformation for each device. Figure 13 depicts the "CCD shift" patterns found to best fit five of the epochs. In these plots, and in all analysis, we removed from the CCD shift coefficients any components consistent with an overall cubic polynomial distortion of the focal plane, since we know that the solution will have time variability of this nature which must be resolved on an exposure-by-exposure basis, e.g., by using Gaia reference stars. This figure makes it clear that there are epoch-to-epoch changes that would dominate the stochastic errors even in a single 30 s exposure. The relative motions of CCDs can be surprisingly large, e.g., over 100 mas, or >6 μm in the focal plane, over a time period of just 14 days, in the case of the last two epochs plotted.

Zoom In Zoom Out Reset image size

The star-flat sequences were taken too infrequently to resolve the temporal behavior of the CCD shifts. Fortunately, the Dark Energy Survey observing program includes repeated visits to 10 fields in a search for high-redshift supernovae. Each field is imaged roughly once per week during the 6-month DES observing season. We examine here the stellar detections in $\approx 1500$ exposures taken of the SNC3 field in the 4 years following camera commissioning. Visits to this field usually comprise 3 × 200 s exposures in g band, 3 × 400 s in r, 5 × 360 s in i, and 11 × 330 s in $z.$ After matching all of the stellar detections in these images, we run WcsFit on the 1123 i and z band exposures from 112 distinct nights during which the photometric solutions indicate absence of clouds and the seeing D₅₀ is predominantly $\lt 1\buildrel{\prime\prime}\over{.} 6$. This WcsFit adopts the astrometric solution derived from the star flats, holding all parameters fixed except

• the linear solution for each exposure,
• the differential chromatic refraction constant for each exposure, and
• the linear CCD shifts, one per CCD per epoch (night) of observations.

From this solution we extract the CCD shifts for each SNC3 epoch, combining them with the CCD shifts for the star-flat data, and projecting out a FOV-wide cubic polynomial fit to each epoch. We analyze only the 59 CCDs that are fully functional over the 4 years.

The first row of Figure 14 plots the difference between each epoch's CCD shifts and the preceding epoch's solution. We quantify this difference by giving the rms displacement between the solutions, averaged over the active regions of the array. We immediately see that the largest changes occur for the first epoch to follow a warming or cooling event for the camera. DECam is cycled to room temperature for various maintenance purposes or when electrical power is lost for long periods. On three occasions, the focal plane temperature dropped from the normal −100 C to −120 C when power to its heaters was lost for several hours. We will refer to these as "camera events" and the periods between them as "camera intervals."

Zoom In Zoom Out Reset image size

We perform a principal components analysis of the 354-element feature vector specifying each exposure's CCD shifts, in hopes of revealing the temporal structure of the largest contributors to astrometric variation. We should be aware, however, that two spurious signals will be present in these data:

• 1.  
  The SNC3 exposures are taken with minimal dithering, and have only $\approx 20$ high-S/N stars per CCD, many fewer than the star-flat data. This means that the CCD shift fits will be pulled by the proper motions of the stars within each device. We should therefore expect to see one or more principal components (PCs) containing a signature that is linear in time for the SNC3 exposures and absent from the star-flat exposures. Parallax motion of the SNC3 stars should be a small perturbation to this.
• 2.  
  The SNC3 CCD shifts have a different reference epoch than the star-flat solution's. Therefore, we should see a static difference between these two in at least PC.

Figure 14 presents the coefficients of the six most significant PCs, plotted against date of the epoch solution. For this purpose, the PCs are scaled to have unit rms displacement across the FOV, so that the coefficients represent typical focal-plane shifts. We immediately identify PC0, the largest contributor, as the expected signature of proper motion in the SNC3 stars (although it is also possible that a smaller, linear-with-time focal plane drift is also in this PC). We also suspect that PC2 contains the second expected spurious signal, the distinct reference epochs for the SNC3 and star-flat solutions.

Most striking is that the other PCs appear to be dominated by changes that occur at camera events. The largest, PC1, began on Thanksgiving day 2015, and its associated distortion persisted at normal operating temperatures until the camera was warmed on 2016 February 19–23. Note that this warmup occurs between the last two star-flat epochs plotted in Figure 13. Figure 15 plots the change in focal plane mapping that occurred during this cooldown. Some of the CCDs appear to have moved by up to 100–150 mas, or 6–10 μm. CCDs also show significant rotations, contributing ≈20 mas rms displacement. Scale changes or shears of the CCDs are much smaller (≲3 mas rms), as expected if the cooldown distorted the mounting structure.

Zoom In Zoom Out Reset image size

Concluding that the bulk of the change in astrometric maps occurs during camera events, we adopt a scheme for astrometric calibration of DES data whereby each exposure is assumed to have the same instrumental astrometric map as the nearest-in-time star-flat epoch that lies within the same camera interval.

We test this DES astrometric procedure by re-running WcsFit on the SNC3 iz-band exposures using a model with the CCD shifts assigned per the local star-flat data, and allowing each exposure a free cubic polynomial distortion across the focal plane.³⁸

Figure 16 illustrates the quality of the resulting fits. In the top row, we plot the rms deviation of all stellar residuals from a given night relative to the mean within the night. This is a measure of astrometric map components that vary on timescales of hours or less, which we expect to be dominated by atmospheric effects. At 5–7 mas rms per component, these are roughly consistent with the values seen in Figure 12 after considering the $\sqrt{330/30}\approx 3.3\times$ reduction in atmospheric noise expected from the longer exposure times in the SN field. It is also possible that some of this rms is due to unmodelled instrumental effects that vary across the 10–20'' dithers of the supernova exposures, e.g., the pixel-to-pixel fluctuations in CCD gate lithography.

Zoom In Zoom Out Reset image size

The lower panel of Figure 16 plots the rms deviation of each night's stellar exposures after we average the night's measurements of each star and subtract the rms noise expected from the intra-night variations. This yields an estimate of the rms astrometric error that is coherent through a night, such as might be attributable to unmodelled shifts in the CCD positions or changes in optical alignment. Such errors are seen to be in the 2–4 mas range on most nights. The errors appear larger in the E-W direction than the N-S direction, particularly during early Y1 observations when the E-W errors reach 6 mas rms, perhaps indicating lower quality in the star-flat solution derived for the corresponding camera interval.

Finally, we note that the cubic polynomials we fit to each exposure are often much larger in amplitude than could be ascribed to atmospheric effects. This suggests that changes in optical alignment over time are significant at ∼100 mas level.

Given an astrometric reference catalog with errors at mas scale and $\gt 1$ star in each 5–10' coherence patch, one could measure some fraction of the atmospheric (or other) astrometric errors and add them to the solution, i.e., interpolate the map between reference stars. The Gaia catalog will provide such a reference catalog. The Gaia DR1 secondary catalog (The Gaia Collaboration 2016) does not contain proper motions so falls slightly short of our ideal, but these will appear in the DR2 release scheduled for 2018 April.³⁹

Reference stars could also be obtained by repeated ground-based observations to average atmospheric and instrumental effects. We use this approach for a cursory investigation of the potential of reference catalog interpolation.⁴⁰ We extract as a reference catalog the mean positions over all star-flat observations of a randomly selected set of stars with mean density of 0.75 arcmin⁻². These "truth" positions are used to interpolate the astrometric distortions for individual exposures. Because our star-flat observations span multiple years, the truth positions may also be degraded by proper motions.

We use the scikit-learn implementation of Gaussian process (GP) regression to interpolate the errors in the astrometric model on a given exposure. The GP technique requires a kernel specifying the covariance between the error vectors of two stars separated by ${\bf{x}}$. We take this covariance function to have a white-noise (δ-function) component of amplitude ${(4\mathrm{mas})}^{2}$ plus a Gaussian with amplitude ${(3\mathrm{mas})}^{2}$ in the cross-wind direction and ${(15\mathrm{mas})}^{2}$ in the wind direction. The cross- and along-wind components of ${\bf{x}}$ have independent GP models. The procedure is to

• 1.  
  Randomly select a training set of stars at the chosen density and fit the GP model to these.
• 2.  
  Use the GP to interpolate to the location of each training star, and reject training stars with outlier residuals (e.g., high-proper-motion stars).
• 3.  
  Refit the GP using the retained training stars.
• 4.  
  Interpolate to the positions the validation set of remaining high-S/N stellar detections.
• 5.  
  Remove outlying residuals from the validation set.
• 6.  
  Calculate the two-point correlation functions of the residuals.

We execute this process for 21 z-band exposures in epoch 20130829, the same set plotted in Figure 11. The mean correlation functions before and after interpolation are plotted in Figure 17. The E mode (curl free) remains dominant even though the interpolation process is not designed to conserve E/B behavior. As expected, the interpolation reduces ${\xi }_{E}$ to negligible levels (≲1 mas²) at scales $\gt 3^{\prime}$ where multiple reference stars can contribute to interpolation. The removal of large-scale power reduces the ${\xi }_{E}(r)$ at $r\to 0$ by a factor ≈8 from the pre-interpolation value. The average post-interpolation residuals are <7 mas rms for this epoch, which has typical stochastic signal level ${\xi }_{0}$. The correlation length of the astrometric errors is reduced to 1'. One would expect the amplitude of the post-interpolation residuals to decrease with the square root of integration time until the systematic error floor of either the reference catalog or DECam is reached.

Zoom In Zoom Out Reset image size

This is just an initial investigation; we have attempted to optimize the interpolation procedure neither for accuracy nor speed. Certainly, there is improvement to be had through GP kernel optimization or other approaches, including interpolation schemes that exploit the known absence of B modes in the vector distortion field. Doubling the density of reference stars appears to have little effect on the residual amplitude.

An astrometric model for DECam with errors near mas level requires terms not only for the polynomial optical "plate solution," but also contributions from: stray electric fields near the edges in the detector and from "tree ring" impurity fluctuations; lateral color and differential chromatic refraction in the bluer bands; shifts in the CCD positions primarily accrued during focal-plane temperature excursions; and time-variable low-order (cubic) distortions across the FOV from a litany of instrumental and atmospheric effects.

All of these distortions components are determined reliably by fitting a model to stellar positions measured from dithered DECam exposures. External reference catalogs play little role in this process, being needed only to stabilize some large-scale degeneracies such as the overall pixel scale. The WcsFit software that we created for this purpose is similar to the widely-used scamp code in optimizing the parameters of a model to maximize agreement among multiple exposures of the same star. WcsFit uses simple linearized iterations to minimize a ${\chi }^{2}$, relying on scamp or some other code to have produced an initial solution that maps each exposure to $\lesssim 1^{\prime\prime}$ accuracy. WcsFit complements scamp by: the ability to specify and fit a complex model with many components interlacing their effects among many exposures; enhanced outlier rejection, necessary to achieve precise modeling; and some optimizations for fitting large exposure sets with large numbers of free parameters.

Once this model is fit to an ensemble of DECam exposures, the remaining astrometric errors are dominated by a curl-free stochastic field of atmospheric refraction fluctuations. In a typical 30 s exposure, the stochastic atmospheric distortions are 10–30 mas rms with coherence length of 4–10' and a strongly anisotropic pattern from wind-blown turbulence. Some nights are much worse than this; unfortunately there is no strong connection between seeing FWHM and astrometric quality.

The atmospheric turbulence averages down with longer exposures or through stacking of residuals on many exposures. Doing so reveals weaker but persistent errors in the astrometric model. Fixed patterns in the devices at 2–4 mas rms (0.008–0.015 pixel) are dominated by larger residuals on small regions of the device subject to edge effects and mounting structures. These could be tabulated from the device stacks and added to the model if we acquired even more stellar measurements.

Star-flat exposures sequences taken every few months are used to monitor shifts in CCD positions. Using DES supernova-field observations, we determine that the bulk of the observed shifts occurs when DECam's focal plan warms to room temperature or cools below normal operating temperature. If we apply the CCD shifts measured in the star-flat epochs to the SN data, we find that remaining errors inter-night variation in the solution is 2–4 mas.

The 4'–10' coherence length of the dominant atmospheric distortions suggests that the Gaia reference catalog, with positions and (in the future) proper motions for $\approx 1$ star per arcmin² at high Galactic latitude, can be used to constrain and remove the atmospheric pattern (and, trivially, the low-order polynomial distortions). Indeed we find that a trial of Gaussian-process interpolation using reference stars at this density reduces the correlation function ${\xi }_{+}(r)$ of errors to <1 mas² on scales $r\gt 3^{\prime}$ and reduce the rms value at smaller scales to <7 mas in a 30 s exposure.

The astrometric maps and procedures derived herein will be applied to the DES observations and made available to other users of DECam. We find that errors in the resultant positions are likely to be dominated by (in order of decreasing importance):

• Unavoidable shot noise in the measurement of object centers (for fainter stars).
• Atmospheric turbulence of typical rms amplitude $(15\mbox{--}20)\times \sqrt{30/T}$ mas, where T is the exposure time in seconds, and coherence length $\approx 10^{\prime} .$ This amplitude depends on the weather. Given a Gaia catalog with proper motions, the atmospheric field can be measured and interpolated to random locations, leaving residuals $\approx 3\times$ smaller and with coherence length $\approx 1^{\prime} .$
• Unmodelled night-to-night variations in the DECam astrometric solution at 2–4 mas rms.
• Up to 2–4 mas rms additional errors from static detector effects (mounting holes, edge fields) that we do not yet model, but could potentially include given a much larger set of stellar measurement residuals.

We conclude that the DECam astrometric model, with registration to the Gaia catalog, has rms errors below 10 mas in a typical 30 s exposure, small enough to be negligible for cosmic-shear measurements and likely to be even smaller in the standard 90-second DES exposure. For a general-use, wide-field instrument like DECam to reach the limit of astrometric accuracy imposed by atmospheric turbulence (with Gaia interpolation gains), the best observing scheme is to dither successive exposures so that the few mas of remaining systematic camera-centered distortions are sampled differently for each exposure of the desired targets. Such a strategy is intrinsic to the DES Wide 5000 deg² survey, so we should expect astrometric catalogs from this survey that are limited by the combination of image shot noise and atmospheric turbulence.

The LSST aims to achieve this goal as well, and the DECam results here show that this should be entirely feasible. LSST has a larger field, shorter exposures, and many more stellar detections to work with. The DECam experience perhaps shows the value of regular star-flat observation sequences. A substantial complication for LSST (as well as other modern wide-field imagers/telescopes such as Hyper Suprime-Cam) is its alt-az mounting and consequent need of an instrument rotator. This introduces a degree of freedom to the optical system absent from the equatorial-mounted Blanco telescope, perhaps greatly increasing the number of constraints that must be analyzed to yield a solution valid at all rotator angles. LSST will obtain many more stellar images, so the necessary data will likely exist but pose a bigger computational challenge.

The WcsFit code is part of the public repository https://github.com/gbernstein/gbdes, where further documentation of its function and use may be found.

G.M.B. gratefully acknowledges support from grants AST-1311924 and AST-1615555 from the National Science Foundation, and DE-SC0007901 from the Department of Energy. A.A.P. is supported by the Jet Propulsion Laboratory. Part of the research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. We thank Steve Holland and Greg Derylo for help with interpretations of detector and camera behaviors, and the anonymous referee for improving the clarity of this paper.

The WcsFit code and subsequent analyses make extensive use of the following excellent public software packages: tmv for linear algebra and TreeCorr for fast correlation functions, both by R. M. Jarvis; yaml-cpp by J. Beder for YAML decoding; and cfitsio by W. Pence, fitsio by E. Sheldon, and AstroPy for FITS access in C and Python environments.

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey.

The collaborating institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

The DES data management system is supported by the National Science Foundation under grant number AST-1138766. The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2015-71825, ESP2015-88861, FPA2015-68048, SEV-2012-0234, SEV-2012-0249, and MDM-2015-0509, some of which include ERDF funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya.

We calculate the correlation functions of the curl-free and divergence-free components of the residual astrometric distortion field ${\rm{\Delta }}{\bf{x}}$ on each exposure, given an irregular, noisy sampling of this field by stellar detections. This is closely analagous to the E/B decompositions performed on the spin-2 polarization field of the cosmic microwave background and the spin-2 weak gravitational shear field in many cosmological investigations. We can derive the vector E/B decomposition by a very slight alteration to the shear-field derivation given by Schneider et al. (2002).

We start with a 2d vector field ${\bf{v}}=({v}_{x},{v}_{y})$. It is useful to work with a complex notation $v={v}_{x}+{{iv}}_{y}$ and complex derivatives $\partial ={\partial }_{x}+i{\partial }_{y}.$ We can write an arbitrary vector field as

$$\begin{eqnarray}&&v=\partial \phi ({\bf{x}}),\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\phi \equiv {\phi }_{E}+i{\phi }_{B}.\end{eqnarray} \tag{ 26 }$$

The curl-free E mode of ${\bf{v}}$ is sourced by ${\phi }_{E}$ and the divergence-free B mode by ${\phi }_{B}.$

The two-point correlation functions of ${\bf{v}}$ at separation vector ${\bf{r}}$ are best posed in terms of the quantity ${v}_{\parallel }+{{iv}}_{\perp }={{ve}}^{-i\beta }$, where β is the position angle of ${\bf{r}}$. We define

$$\begin{eqnarray}&&{\xi }_{+}({\bf{r}})={\xi }_{+}(r,\beta )=\langle v({\bf{x}}){e}^{-i\beta }{[v({\bf{x}}+{\bf{r}}){e}^{-i\beta }]}^{\star }\rangle \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{\xi }_{-}(r,\beta )+i{\xi }_{\times }(r,\beta )=\langle v({\bf{x}}){e}^{-i\beta }v({\bf{x}}+{\bf{r}}){e}^{-i\beta }\rangle \end{eqnarray} \tag{ 28 }$$

Taking $\tilde{\phi }({\bf{k}})$ to be the Fourier transform of ϕ, and the generation of ϕ to be a stationary stochastic process, we define the power spectra via

$$\begin{eqnarray}\langle {\tilde{\phi }}_{E}({{\bf{k}}}_{1}){\tilde{\phi }}_{E}^{\star }({{\bf{k}}}_{2})\rangle & = & {(2\pi )}^{2}\delta ({{\bf{k}}}_{1}-{{\bf{k}}}_{2}){P}_{{EE}}({{\bf{k}}}_{1})\\ \langle {\tilde{\phi }}_{E}({{\bf{k}}}_{1}){\tilde{\phi }}_{B}^{\star }({{\bf{k}}}_{2})\rangle & = & {(2\pi )}^{2}\delta ({{\bf{k}}}_{1}-{{\bf{k}}}_{2}){P}_{{EB}}({{\bf{k}}}_{1})\\ \langle {\tilde{\phi }}_{B}({{\bf{k}}}_{1}){\tilde{\phi }}_{B}^{\star }({{\bf{k}}}_{2})\rangle & = & {(2\pi )}^{2}\delta ({{\bf{k}}}_{1}-{{\bf{k}}}_{2}){P}_{{BB}}({{\bf{k}}}_{1}).\end{eqnarray} \tag{ 29 }$$

P_EB must be real if the vector field statistics are invariant under 180° rotation, so we will assume this is true. The real part will vanish as well if the process generating ${\bf{v}}$ is invariant under parity flips. We will leave the real part as a free parameter.

By propagating the derivatives in Equation (25) through a Fourier transform we can express the correlation functions (27) and (28) as

$$\begin{eqnarray}{\xi }_{+}(r,\beta ) & = & \displaystyle \int d\alpha \displaystyle \int k\,{dk}| {k}^{2}| {e}^{{ikr}\cos (\alpha -\beta )}\\ & & \times [{P}_{{EE}}(k,\alpha )+{P}_{{BB}}(k,\alpha )],\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{\xi }_{-}(r,\beta )+i{\xi }_{\times }(r,\beta )\\ &&\quad =\displaystyle \int d\alpha \displaystyle \int k\,{dk}| {k}^{2}| {e}^{-{ikr}\cos (\alpha -\beta )}{e}^{2i(\alpha -\beta )}\\ &&\qquad \times [{P}_{{EE}}(k,\alpha )-{P}_{{BB}}(k,\alpha )+2{{iP}}_{{EB}}(k,\alpha )],\end{eqnarray} \tag{ 31 }$$

where α is the position angle of the wavevector ${\bf{k}}$. Even though the atmospheric distortions are anisotropic, we will concern ourselves only with the angle-average quantities ${\xi }_{+}(r)\,=\langle {\xi }_{+}(r,\beta ){\rangle }_{\beta },$ and the corresponding angle-averaged power spectra ${P}_{{EE}}(k),$ etc. If we average the preceding equations over β, Bessel's first integral yields

$$\begin{eqnarray}&&{\xi }_{+}(r)=2\pi \int k\,{dk}\,{k}^{2}{J}_{0}({kr})[{P}_{{EE}}(k)+{P}_{{BB}}(k)]\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{\xi }_{-}(r)=-2\pi \int k\,{dk}\,{k}^{2}{J}_{2}({kr})[{P}_{{EE}}(k)-{P}_{{BB}}(k)]\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{\xi }_{\times }(r)=-4\pi \int k\,{dk}\,{k}^{2}{J}_{2}({kr}){P}_{{EB}}(k).\end{eqnarray} \tag{ 34 }$$

The last equation tells us that ${\xi }_{\times }$ is produced purely by EB power. Let us define pure-E and pure-B quantities

$$\begin{eqnarray}&&{\xi }_{E,B}(r)\equiv 2\pi \int k\,{dk}\,{k}^{2}{J}_{0}({kr}){P}_{{EE},{BB}}(k)\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{1}{2}\left\{{\xi }_{+}(r)\mp \int r\,{dr}\,{J}_{0}(r)\int {r}^{\prime }\,{{dr}}^{\prime }{J}_{2}({r}^{\prime }){\xi }_{-}({r}^{\prime }).\right\}\end{eqnarray} \tag{ 36 }$$

where the last line combines Equation (32) with the order-2 Hankel Transform of Equation (33). After making use of identities 9.1.27 and 11.4.42 from Abramowitz & Stegun (1965), this can be converted to

$$\begin{eqnarray}&&{\xi }_{E,B}(r)=\displaystyle \frac{1}{2}\left\{{\xi }_{+}(r)\pm \left[{\xi }_{-}(r)-2{\int }_{r}^{\infty }{{dr}}^{\prime }\displaystyle \frac{1}{{r}^{\prime }}{\xi }_{-}({r}^{\prime })\right]\right\}.\end{eqnarray} \tag{ 37 }$$

Equations (34) and (37) allow us to produce measures of pure E, B, and cross-EB power from two-point correlations constructed from all pairs of detections in a given exposure. In Figure 18, we plot ${\xi }_{E}$ and ${\xi }_{B}$ inferred for the astrometric residuals in each of 20 consecutive star-flat exposures (after projecting out a cubic polynomial function of field coordinates from each exposure). It is clear that the astrometric residuals are indeed dominated by E modes, e.g., curl free. Figure 11 plots the mean of ${\xi }_{E},{\xi }_{B},$ and ${\xi }_{\times }$ in another set of exposures, confirming that any divergence modes are very small. In further analysis we assumed ${\xi }_{B}=0$ such that we can more simply take ${\xi }_{E}={\xi }_{+}$.

Zoom In Zoom Out Reset image size