Brighter-fatter effect in near-infrared detectors -- I. Theory of flat auto-correlations

TL;DR

This paper develops a formalism to measure the brighter-fatter effect in near-infrared detectors using flat field auto-correlations, accounting for inter-pixel capacitance and enabling more accurate detector characterization for weak lensing studies.

Contribution

It introduces a new method to distinguish and measure the brighter-fatter effect and inter-pixel capacitance in infrared detectors using non-destructive readouts and correlation functions.

Findings

BFE can be measured separately from IPC in IR detectors.

Simulations show accurate extraction of BFE and IPC parameters.

Biases in BFE coefficients are around 12%, indicating higher order effects.

Abstract

Weak gravitational lensing studies aim to measure small distortions in the shapes of distant galaxies, and thus place very tight demands on the understanding of detector-induced systematic effects in astronomical images. The Wide-Field Infrared Survey Telescope (WFIRST) will carry out weak lensing measurements in the near infrared using the new Teledyne H4RG-10 detector arrays, which makes the range of possible detector systematics very different from traditional weak lensing measurements using optical CCDs. One of the non-linear detector effects observed in CCDs is the brighter-fatter effect (BFE), in which charge already accumulated in a pixel alters the electric field geometry and causes new charge to be deflected away from brighter pixels. Here we describe the formalism for measuring the BFE using flat field correlation functions in infrared detector arrays. The auto-correlation of CCD flat fields is often used to measure the BFE, but because the infrared detector arrays are read out with the charge "in place," the flat field correlations are dominated by capacitive cross-talk between neighboring pixels (the inter-pixel capacitance, or IPC). Conversely, if the BFE is present and one does not account for it, it can bias correlation measurements of the IPC and photon transfer measurements of the gain. We show that because the infrared detector arrays can be read out non-destructively, one can compute numerous cross-correlation functions between different time slices of the same flat exposures, and that correlations due to IPC and BFE leave distinct imprints. We generate a suite of simulated flat fields and show that the underlying IPC and BFE parameters can be extracted, even when both are present in the simulation. There are some biases in the BFE coefficients up to 12%, which are likely caused by higher order terms that are dropped from this analysis.

Tables

Table 1. Table 1: Averaged results for the simulations, based on stacks of flat ramps. These values were obtained with advanced characterization with ncycle=3.

Quantity	Units	Flat type, number			Uncert			Notes
		sim,n3	sim,n10	truth	stat.(3)	stat.(10)	sys.(3)
Charge, $I{t}_{n,n+1}$	ke	1.4607	1.4615	1.4604	0.0006	0.0003
Gain $g$	e/DN	2.0606	2.0620	2.0600	0.0008	0.0004
IPC $\alpha$	%	1.6764	1.6793	1.6900	0.0055	0.0025
IPC ${\alpha }_H$	%	1.6809	1.6806	1.6900	0.0039	0.0018
IPC ${\alpha }_V$	%	1.6720	1.6779	1.6900	0.0038	0.0018
IPC ${\alpha }_D$	%	-0.0002	-0.0021	0.0000	0.0027	0.0012
Non-linearity ${\beta }_{\mathrm{ramp}}$	ppm/e	0.5835	0.5782	0.5800	0.0003	0.0001	0.0091
Alternative intervals
Non-linearity ${\beta }_{\mathrm{ramp}}$	ppm/e	0.5862	0.5794	0.5800	0.0006	0.0003	0.0191	Frames 3,7,9
Non-linearity ${\beta }_{\mathrm{ramp}}$	ppm/e	0.5806	0.5801	0.5800	0.0002	0.0001	0.0052	Frames 3,19,21
Non-overlapping correlation function (Method 1)
BFE Coefficients - frames 3,11,13,21, baseline-corrected
${\left[K^2{a}^{\prime }\right]}_{0,0}$	ppm/e	-1.0373	-1.0192	-1.1590	0.0145	0.0064	0.0103	Central pixel
	ppm/e	0.1838	0.1980	0.2034	0.0073	0.0033		Nearest neighbor
	ppm/e	0.0362	0.0428	0.0505	0.0072	0.0032		Diagonal
	ppm/e	0.0155	0.0133	0.0120	0.0074	0.0032
	ppm/e	0.0049	0.0010	0.0027	0.0052	0.0023
	ppm/e	0.0271	0.0179	0.0185	0.0075	0.0033
BFE Coefficients - frames 3,7,9,13 baseline-corrected
${\left[K^2{a}^{\prime }\right]}_{0,0}$	ppm/e	-1.0400	-1.0293	-1.1590	0.0288	0.0130	0.0216	Central pixel
	ppm/e	0.2381	0.2195	0.2034	0.0152	0.0066		Nearest neighbor
	ppm/e	0.0479	0.0392	0.0505	0.0151	0.0066		Diagonal
BFE Coefficients - frames 3,19,21,37 baseline-corrected
${\left[K^2{a}^{\prime }\right]}_{0,0}$	ppm/e	-0.9156	-0.9214	-1.1590	0.0068	0.0031	0.0059	Central pixel
	ppm/e	0.1850	0.1818	0.2034	0.0034	0.0015		Nearest neighbor
	ppm/e	0.0422	0.0450	0.0505	0.0035	0.0016		Diagonal
Mean-variance relation (Method 2)
${\hat{a}}_{0,0,M2}$	ppm/e	-1.3120	-1.3513	-1.3720	0.0383	0.0144	0.0273
$\beta -4\left(1+3\alpha \right){\alpha }^{\prime }$	ppm/e	0.5613	0.5677	0.5800	0.0218	0.0079
${\sum }_a-8\left(1+3\alpha \right){\alpha }^{\prime }$	ppm/e	-0.0445	-0.0211	0.0000	0.0436	0.0158	0.0182
Adjacent pixel correlations (Method 3)
	ppm/e	0.1816	0.1855	0.2034	0.0072	0.0032

Equations

$${\mathcal{A}}_{i,j}={\mathcal{A}}_{i,j}^0\left[1+\sum _{\Delta i,\Delta j}{a}_{\Delta i,\Delta j}Q(i+\Delta i,j+\Delta j)\right],$$

$$a_{\Delta i,\Delta j}=a_{\Delta i,\Delta j}^R+a_{\Delta i,\Delta j}^T+a_{\Delta i,\Delta j}^L+a_{\Delta i,\Delta j}^B,$$

$$\sum _{\Delta i,\Delta j}{a}_{\Delta i,\Delta j}=0.$$

$${\Sigma }_a=\sum _{\Delta i,\Delta j}{a}_{\Delta i,\Delta j}.$$

$$a_{\Delta i,\Delta j}^{\prime }\equiv a_{\Delta i,\Delta j}-{\delta }_{\Delta i,0}{\delta }_{\Delta j,0}{\Sigma }_a.$$

$$a_{\langle\Delta i,\Delta j\rangle}\equiv\frac{1}{8}\Bigl{(}a_{\Delta i,\Delta j}+a_{\Delta j,\Delta i}+a_{-\Delta j,\Delta i}+a_{-\Delta i,\Delta j}+a_{-\Delta i,-\Delta j}+a_{-\Delta j,-\Delta i}+a_{\Delta j,-\Delta i}+a_{\Delta i,-\Delta j}\Bigr{)}.$$

$$S_{\mathrm{initial}}-S_{\mathrm{final}}=\frac{1}{g}\left(Q-\beta Q^2\right),$$

$$S_{\mathrm{initial}}(i,j)-S_{\mathrm{final}}(i,j)=\frac{1}{g}\sum _{\Delta i,\Delta j}{K}_{\Delta i,\Delta j}{Q}_{i-\Delta i,j-\Delta j}+[\mathrm{nonlinear}\mathrm{terms}],$$

$$\sum _{\Delta i,\Delta j}{K}_{\Delta i,\Delta j}=1.$$

$$S_{\mathrm{initial}}(i,j)-S_{\mathrm{final}}(i,j)=\frac{1}{g}\sum _{\Delta i,\Delta j}[K_{\Delta i,\Delta j}+K_{\Delta i,\Delta j}^{\prime }\bar{Q}]Q_{i-\Delta i,j-\Delta j},$$

$$S_{\mathrm{initial}}(i,j)-S_{\mathrm{final}}(i,j)$$

$$C_{abcd}(\Delta i,\Delta j)=\mathrm{Cov}\left[S_a(i,j)-S_b(i,j),S_c(i+\Delta i,j+\Delta j)-S_d(i+\Delta i,j+\Delta j)\right].$$

$${\mathrm{Cov}}_{\mathrm{meas}}\left[\mathcal{O},{\mathcal{O}}^{\prime }\right]=\frac{1}{2}⟨({\mathcal{O}}_A-{\mathcal{O}}_B)({\mathcal{O}}_A^{\prime }-{\mathcal{O}}_B^{\prime })⟩,$$

$${\bar{C}}_{abcd[n]}(\Delta i,\Delta j)=\sum _{\nu =0}^{n-1}{C}_{a+\nu ,b+\nu ,c+\nu ,d+\nu }(\Delta i,\Delta j).$$

$$⟨Q_a(i,j)⟩=I{t}_a$$

$$\mathrm{Cov}\left[Q_a(i,j),Q_b(i+\Delta i,j+\Delta j)\right]=I{t}_{\min (a,b)}{\delta }_{\Delta i,0}{\delta }_{\Delta j,0}.$$

$$C_{abcd}(\Delta i,\Delta j)=\frac{I}{g^2}\left[t_{\min (a,c)}-t_{\min (a,d)}-t_{\min (b,c)}+t_{\min (b,d)}\right]{\delta }_{\Delta i,0}{\delta }_{\Delta j,0}.$$

$$t_{abcd}\equiv t_{\min (a,c)}-t_{\min (a,d)}-t_{\min (b,c)}+t_{\min (b,d)}\mathrm{and}{t}_{ab}\equiv t_b-t_a;$$

$$C_{abab}(0,0)=\frac{I}{g^2}{t}_{ab}\mathrm{and}⟨S_{ab}(i,j)⟩=\frac{I}{g}{t}_{ab}.$$

$$C_{abcd}(\Delta i,\Delta j)$$

$$C_{abcd}(0,0)$$

$$C_{abcd}(\pm 1,0)$$

$$C_{abcd}(0,\pm 1)$$

$$C_{abcd}(\pm 1,\pm 1)=C_{abcd}(\pm 1,\mp 1)$$

$$⟨Q_a(i,j)Q_b(i,j)Q_c(i,j){⟩}_{\mathrm{conn}}=I{t}_{\min (a,b,c)},$$

$$\mathrm{Cov}[Q_a(i,j),Q_b(i,j{)}^2]=2I^2{t}_b{t}_{\min (a,b)}+I{t}_{\min (\mathrm{a},\mathrm{b})}.$$

$$C_{abcd}(0,0)$$

$$T_{abcd}\equiv -(t_a+t_d)t_{\min (d,a)}+(t_b+t_d)t_{\min (d,b)}+(t_a+t_c)t_{\min (c,a)}-(t_b+t_c)t_{\min (c,b)}$$

$$C_{abcd}(0,0)=\frac{1}{g^2}(1-2\beta )I{t}_{abcd}-2\frac{\beta }{g^2}{I}^2{T}_{abcd}.$$

$$[S_{\mathrm{initial}}-S_{\mathrm{final}}](i,j)$$

$$\Delta C_{abcd}(0,0){|}_{\alpha \beta }=\frac{\alpha \beta }{g^2}(16I^2{T}_{abcd}+12I{t}_{abcd})$$

$$\Delta C_{abcd}(\pm 1,0){|}_{\alpha \beta }=-\frac{4{\alpha }_{\mathrm{H}}\beta }{g^2}(I^2{T}_{abcd}+I{t}_{abcd})$$

Full text

instrumentation: detectors

1 Introduction

Weak gravitational lensing (WL) – the distortion of the shapes and sizes of distant galaxies by the curvature of the intervening space-time – is a powerful method for probing the matter distribution in the Universe (for recent results, see e.g. Heymans et al. 2013; Abbott et al. 2016; and Hildebrandt et al. 2017). However, the signal is small and must be measured to a fraction of a percent to meet the science goals of current and future WL surveys. Therefore, WL programs place a very strong emphasis on understanding every systematic error that can occur in measuring the shape of the galaxy – this includes the contribution of the atmosphere, optics, and image motion to the smearing of an image, as well as imprints of the detector system and data processing.

The brighter-fatter effect (BFE; e.g. Antilogus et al. 2014) is one of these subtle effects that has been observed in silicon CCD detectors. This is a non-linear effect in which a brighter point source produces a larger image (as measured by e.g. full width at half maximum) in the CCD than a fainter point source. It is caused by changes in the electric field geometry in the CCD as a well fills up with electrons: if at any instant during the exposure a pixel $(i,j)$ contains more charge than its neighbors, then due to self-repulsion of the electrons, additional photo-electrons generated will be less likely to land in pixel $(i,j)$ and more likely to land in its neighbors. This is described phenomenologically for CCDs by supposing that the pixel boundaries111A “pixel boundary” in an either astronomical CCD or a NIR detector is not a physical barrier between neighboring pixels, but rather is defined by the electric fields and diffusion coefficients that determine which well ultimately collects a photo-electron. move as a function of accumulated charge. In thick CCDs, the BFE has been observed to have a significant range, e.g. in the Dark Energy Camera (DECam) CCDs, the pixel boundary shifts have been measured from charge up to $\sim 10$ pixels away (Gruen et al., 2015). The BFE also manifests itself in correlation properties (variance and correlation function) of flat-field images, where the shifting pixel boundaries break the usual assumption that each photo-electron behaves independently from previous electrons and hence causes non-Poisson correlations in the flat images (Guyonnet et al., 2015). Indeed, this provided one of the early hints to the existence of the BFE (Downing et al., 2006).

The BFE and techniques for modeling it have been well-established in current WL surveys such as the Dark Energy Survey (DES; Gruen et al. 2015) and the Hyper-Suprime Cam (HSC; Coulton et al. 2018). Higher precision will be demanded of the next generation of WL surveys: the Large Synoptic Survey Telescope (LSST), the Euclid space mission, and the Wide-Field Infrared Survey Telescope (WFIRST) space mission. The BFE has been observed and characterized in Euclid development CCDs (Niemi et al., 2015) and in candidate sensors for LSST (Baumer & Roodman, 2015; Lage et al., 2017). In setting requirements on any detector effect, it is important to study all of the ways that effect can enter into the analysis; in the case of the brighter-fatter effect in weak lensing, the stars used for determination of the point spread function (PSF) are much brighter (and have steeper intensity gradients) than either the galaxies or the sky background, and so the stars are most affected by the BFE. This means that the principal effect of the BFE is to make the measured (star-based) PSF larger than the correct PSF for faint galaxies. Subsequent stages of analysis will then over-correct for the smearing out of galaxy ellipticities by the finite size of the PSF, and hence will over-estimate the shear signal. Note that we consider the BFE to be a calibration problem, in the sense that WFIRST will need to develop a model for it; we do not need to eliminate it.

WFIRST plans to measure weak lensing in the near infrared (NIR), and thus silicon CCDs are not an option. Instead, it will use infrared detector arrays: each of the 18 detector arrays will use Teledyne’s H4RG-10 readout integrated circuit ($4088\times 4088$ active pixels, 10 $\mu$m pitch) with mercury cadmium telluride (HgCdTe; $2.5$ $\mu$m cutoff) as the light-absorbing component.222For more background on the H$x$RG series devices, we refer the reader to the overview (Beletic et al., 2008; Blank et al., 2011), and reports on the JWST/NIRSpec H2RGs (Rauscher et al., 2007, 2014) and on the WFIRST H4RG development (Piquette et al., 2014). Just as for CCDs, the boundaries between adjacent pixels are defined by the solution to the drift-diffusion equation rather than a physical barrier, and so a BFE in WFIRST detectors would be physically plausible; however one would expect the details to be very different. It is therefore important to understand whether the BFE, or other non-linear detector effects, are present in WFIRST prototype devices, and if so how the BFE fits into the overall WFIRST calibration plan. A discussion of the physics of the BFE as applied to NIR detectors, as well as some characterization efforts on earlier generations of Teledyne detectors (H1RG and H2RG, 18 $\mu$m pitch) can be found in Plazas et al. (2017). The BFE has been measured in point source illumination tests on an H2RG tested for the Euclid program (Plazas et al., 2018).

In order to develop calibration and data reduction procedures for WFIRST data, it is essential to characterize non-linear effects such as the BFE. Such characterization efforts also provide us with the opportunity to learn about what calibration procedures work and what pitfalls exist. Such studies must be undertaken early in the life cycle of the project, especially in the context of a space mission where late changes to the calibration requirements could be expensive or impossible. Detailed characterization of subtle effects such as the BFE is also outside the scope of the technology development milestones (Spergel et al., 2015). These milestones were recently completed for the WFIRST NIR detectors, but focus on basic performance (e.g. quantum efficiency, read noise, persistence), production yield, and environmental testing (e.g. thermal cycling, vibrational, radiation). However, the studies for these milestones do not explore the detector-related systematic effects that are relevant to control shear errors at the few$\,\times 10^{-4}$ level.

The two major methods of measuring the BFE are (i) measurement of spots projected onto the detector (either with a laboratory spot projector or using real stars observed through a telescope) and (ii) flat field statistics. This paper considers flat fields, since spots were not yet available for the H4RG-10 HgCdTe detectors at the time we began this project.333Jay Anderson (private communication) has presented some results on the BFE in the HST WFC3-IR channel using observations of stars to the WFIRST Detector Working Group. While WFC3-IR is a useful guide to some of the issues WFIRST will encounter, it is an 18 $\mu$m pitch device and so it is important to measure BFE parameters on H4RG-10 HgCdTe devices. The data were acquired as part of general detector characterization and made available to the science teams; the test procedure was not specifically optimized for BFE studies.

The use of flat field statistics for BFE measurements presents some special challenges for infrared arrays, most notably that the flat field auto-correlation function is dominated by the effect of inter-pixel capacitance (IPC), which gives a positive correlation between adjacent pixels (Moore et al., 2004). IPC can have a linear component and a non-linear component, the latter of which is known as non-linear inter-pixel capacitance (NL-IPC, a signal dependent coupling that occurs when one converts from charge to voltage; e.g., Cheng 2009; Donlon et al. 2016, 2017, 2018).444There are also reports of NL-IPC in H2RG detectors (Arielle Bertrou-Cantou, private communication). On the other hand, the H4RG-10 provides a non-destructive read capability, which is very useful for BFE studies as it enables intermediate stages of the image to be observed: the flat field is then a 3D data cube, and correlations across different time slices (“frames”) of the image can be measured. We will find that this capability allows us to simultaneously measure the IPC and BFE using flat fields. We will further find that the gain measurement from the photon transfer curve (e.g., Mortara & Fowler, 1981; Janesick et al., 1985) must be corrected for the BFE in addition to the now-standard IPC correction (Moore et al., 2006; Brown et al., 2006).

The cleanest method for BFE measurements from H4RG-10 flat fields is to cross-correlate two correlated double sample (CDS) images, obtained from non-overlapping parts of the ramp. This eliminates any possible correlations from Poisson noise – including those that couple between neighboring pixels through linear IPC – as well as any read noise correlations that occur within a single frame. It does leave correlations due to classical non-linearity, which must be removed based on the standard non-linearity curve analysis. The non-overlapping correlation function method cannot tell the difference between the brighter-fatter effect (which occurs during the process of collecting charge into a well) from NL-IPC. To distinguish the BFE from NL-IPC, we must resort to correlations of CDS images over the same (or at least an overlapping) time interval, and observe how pixel variances or covariances of adjacent pixels change as one varies the time interval of interest. While these tests mix together many different detector effects, the BFE and NL-IPC hypotheses make distinct predictions. NL-IPC appears in these methods with a factor of 2 different from BFE, because in the BFE each electron555The Teledyne HxRG detectors actually collect holes, rather than electrons, although the statistical techniques in this paper are agnostic to the sign of the charge. We will use the standard nomenclature of charge in “electrons” – this is common usage in the astronomical community and HST WFC3-IR documentation, even though the WFC3-IR detector also collects holes. collected only affects the behavior of subsequent electrons, but in NL-IPC (which acts on collected charge) every electron affects every other electron.666In §3, the mathematics of this is worked out in great detail, but the underlying reason for the factor of 2 is the simple combinatorial effect.

This paper is organized as follows. In §2, we build up our description of the brighter-fatter effect as well as other detector effects relevant to the flat field (IPC and non-linearity) and the formalism for correlation functions among the frames of a flat field. In §3, we work out the theoretical predictions for the 2-point correlation function of the flat fields in the presence of the various effects. In §4, we describe a simulation incorporating IPC, non-linearity, and the BFE that allows us to test the characterization methods in §5. We conclude in §6. The appendix contains some technical material on the covariances of clipped data (Appendix A). The application to laboratory H4RG-10 data – and associated evidence for the BFE – is presented in a companion paper (Paper II).

2 Formalism

2.1 Brighter-fatter effect

Autocorrelation measurements are sensitive to the BFE via changes in the effective pixel area.777In the BFE literature, the operational definition of pixel “area” is that ${\cal A}_{i,j}={\rm QE}_{\rm ref}^{-1}\int p_{i,j}(x,y)\,dx\,dy$, where $p_{i,j}(x,y)$ is the probability that a photon incident at position $(x,y)$ on the detector leads to an electron collected in pixel $(i,j)$. This integral is divided by a reference value of the quantum efficiency, QEref, so that at low flux levels the sum of pixel areas in some region corresponds to the geometrical area. If the probability of collecting a charge in any well $\sum_{i,j}p_{i,j}(x,y)$ remains fixed, then the BFE conserves total pixel area. We suppose that a pixel $(i,j)$ has effective area that changes depending on the charge in neighboring pixels:

[TABLE]

where $i$ and $j$ denote column and row indices (0…4095 for the H4RG), $Q({i,j})$ is the charge (in number of elementary charges) in pixel $(i,j)$, and $a_{\Delta i,\Delta j}$ denotes a coupling matrix. In a general case, one might allow $a$ to also depend on $i$ and $j$, which would correspond to a BFE that varies from one pixel to another. However with flat autocorrelation data and a plausible number of flats, it is possible only to measure averages of the $a$ coefficients in groups of pixels. Therefore we will assume discrete translation invariance; note, however, that the ability to perform autocorrelations with different sets of pixels allows some (limited) ability to test for translation invariance.

While $a_{\Delta i,\Delta j}$ is formally dimensionless, we will normally quote $a_{\Delta i,\Delta j}$ in units of 10*-6* e*-1*, ppm/e, or %/$10^{4}$ e (all of which are equivalent). These units are convenient because $10^{4}$ e is a typical integrated signal level in the central pixel of a PSF star for WFIRST, so a measured value of $a$ in ppm/e maps into the expected order of magnitude of the effect on a star in percent.

In a phenomenological BFE model, one specifies how much of the area change comes from each of the boundaries by writing

[TABLE]

where the superscripts $R$, $T$, $L$, and $B$ refer to the right, top, left, and bottom boundaries respectively. If the quantum efficiency depends on the charge in the well, then we would include an additional term $a_{\Delta i,\Delta j}^{QE}$. An image simulation requires as input all of these components individually, and they can be probed with spot illumination or individual pixel resets; however flat autocorrelations are only sensitive to the total.

A true BFE that works by moving pixel boundaries conserves total area, so we should have

[TABLE]

(We expect that $a_{0,0}$ would be negative, and the $a$ coefficients for the neighbors would be positive.) However, the sum in Eq. (3) is ill-behaved, since the noise diverges as we continue to add pixels. It can therefore be tested only in the context of fitting a model to $a_{\Delta i,\Delta j}$. Moreover, if adding charge to a pixel changes the QE or charge collection probability, then Eq. (3) may be violated. Therefore, it is important to measure all of the $a_{\Delta i,\Delta j}$, without assuming Eq. (3). In general we will define:

[TABLE]

We will find it useful later to define

[TABLE]

By definition, the $a^{\prime}$ coefficients sum to zero. Some of the BFE tests that we will conduct are not sensitive to ${\Sigma_{a}}$, and hence can only measure the $a^{\prime}_{\Delta i,\Delta j}$.

Finally, we define symmetry-averaged functions over coefficients $(\Delta i,\Delta j)$ related by the rotation or reflection symmetries of the pixel grid:

[TABLE]

There is no law of physics requiring the BFE to respect the rotation and reflection symmetries (indeed, in CCDs it does not), so some test results are provided for, e.g. $a_{\Delta i,\Delta j}$ and $a_{\Delta j,\Delta i}$ separately.

2.2 Gains, nonlinearities, and IPC

“Raw” data from the detector arrays are not in electrons but in data numbers (DN), which are voltages quantized as 16-bit integers. As each pixel is exposed to light, the voltage difference across the photodiode decreases; for detectors that collect holes, the voltage on the readout (p-type) side of the diode increases. The observed signal $S$ (units: DN) may increase or decrease depending on the polarity of the analog-to-digital converter.888We have worked with raw data of both polarities, and it is easy enough to switch between them by inserting an optional mapping $S\rightarrow 2^{16}-1-S$ in the input script. In this paper we work with the convention that $S$ decreases during integration. Ideally the relation between the accumulated charge and signal drop would be linear, but in practice it is not. This effect can contain contributions both from the non-linearity of the $p-n$ junction itself as well as any step in the readout chain, and is generically observed in NIR detectors (e.g. Bohlin et al., 2005; Deustua et al., 2010; Hill et al., 2010); see Plazas et al. (2016) for a study of its impact on the WFIRST weak lensing program.

A polynomial model is typically used to describe the non-linearity curve; the most important correction is typically the quadratic term. In flat illumination, where each pixel accumulates charge $Q$, the drop in signal level is given by

[TABLE]

where $g$ is the gain (units: e/DN) and $\beta$ is the leading-order non-linearity coefficient. Note that $\beta$ has the same units as ${a}$, and so it will be convenient to quote it in ppm/e. We define “initial” for the purposes of Eq. (7) to mean immediately following a reset, which we take to be $t=0$. Note that non-linearity parameters will likely depend on the reset voltage.

In thick CCD detectors, it is common to use auto-correlations of the flat field to measure the BFE. However, in infrared detectors the auto-correlation of a flat is instead dominated by inter-pixel capacitance (IPC). IPC is an electrical coupling between neighboring pixels, in which the voltage on one pixel is sensitive to the charge in its neighbors (e.g. Moore et al., 2004, 2006). This coupling increases the apparent size of the image of a star on the detector; see Kannawadi et al. (2016) for a study of its impact on WFIRST. This means a more complex procedure is needed to probe the BFE in infrared arrays. Furthermore, determination of the gain from variance vs. mean plots must be corrected for IPC to obtain meaningful results (e.g. Moore et al., 2006; Fox et al., 2008; Crouzet et al., 2012). Normally the IPC is described by replacing the linear term in Eq. (7) with a kernel describing the capacitive cross-talk among the pixels:

[TABLE]

where the kernel matrix ${\bf K}$ satisfies

[TABLE]

In the case where the IPC only talks to the nearest neighbors and does so equally, we have $K_{0,0}=1-4\alpha$, $K_{0,\pm 1}=K_{\pm 1,0}=\alpha$, and all others are zero. However, asymmetries between the horizontal and vertical directions ($K_{0,\pm 1}\neq K_{\pm 1,0}$) are commonly observed in NIR detectors. Therefore we measure separately $\alpha_{\rm H}=K_{\pm 1,0}$ and $\alpha_{\rm V}=K_{0,\pm 1}$; if these are different then we define $\alpha$ to be their average $(\alpha_{\rm H}+\alpha_{\rm V})/2$. We will also allow for diagonal IPC, $\alpha_{\rm D}=K_{\pm 1,\pm 1}$ (when this notation is used, we will not distinguish between the “northeast-southwest” and “northwest-southeast” directions, although in principle their IPC may be different).

Inter-pixel capacitance in a semiconductor device may depend on signal since such devices do not obey the principle of superposition. This “non-linear inter-pixel capacitance” (NL-IPC) can be phenomenologically similar to the BFE: if the IPC grows with signal level, then this will also lead to brighter stars showing a larger observed FWHM on account of the greater amount of coupling. However, NL-IPC is a different mechanism – it arises in the conversion of charge to voltage, whereas the BFE arises in the collection of charge – and as such there are subtle differences in how it impacts both flat field statistics and science data. Disentangling the two effects proves to be one of the most difficult part of the BFE analysis.

Discussions of NL-IPC are complicated by the fact that NIR detectors both have non-linear charge-to-signal conversion (Eq. 7) and IPC. In the presence of both of these effects, we generically expect some kind of non-linear cross-talk between neighboring pixels of order $\alpha\beta$, and great care is needed to even define a quantitative measure of NL-IPC. In most astronomical data processing, the non-linearity correction is applied to individual pixels as one of the first steps – and certainly before any attempt at IPC correction (if the latter is done at all). This is equivalent to the assumption that all of the non-linearity acts on the signal after IPC. In this paper, we use “NL-IPC” to denote any non-linearity in the charge-to-signal conversion that deviates from this assumption.

For the purposes of flat fields, we parameterize NL-IPC by a mean signal-level-dependent kernel,

[TABLE]

where $\bar{Q}$ is the mean accumulated charge ($It$ in a flat exposure). One can equivalently write this in terms of $\alpha^{\prime}$, $\alpha^{\prime}_{\rm H}$, $\alpha^{\prime}_{\rm V}$, etc.:

[TABLE]

In the flat illumination case, there is no ambiguity of what mean accumulated charge $\bar{Q}$ should be used. In other cases such as spot illumination or pixel reset tests, neighboring pixels can have wells filled to very different levels and $K^{\prime}$ is no longer the appropriate concept. Some studies have indicated that NL-IPC is a function of both contrast and signal level – see, e.g., Donlon et al. (2016, 2017, 2018) – and in this case one should write an IPC coupling constant $\alpha(Q_{i,j},Q_{i+1,j})$ that is a function of charge in both pixels.999If viewed as a capacitor network, voltage in the pixel might be a more fundamental variable than the charge. The flat field test probes the case of $Q_{i+1,j}\approx Q_{i,j}\approx\bar{Q}$, whereas single pixel reset and hot pixel tests measure the case where $Q_{i+1,j}\approx 0$.

2.3 Correlation functions

In a CCD flat, there is only a single read of the detector in each flat exposure. However, in an infrared array flat, one typically obtains $N$ samples up the ramp, and correlation functions can be defined not just between different pixels but between different frames. If one denotes the frames by indices $abc...$ then let us define:

[TABLE]

If one obtained $N=2$ samples and took the autocorrelation of the CDS image $S_{1}-S_{2}$, this would correspond to $C_{1212}(\Delta i,\Delta j)$. This is the procedure that is most similar to a CCD flat autocorrelation. However, as noted above, it contains IPC as well as BFE and therefore cannot distinguish the two. Fortunately, with multiple up-the-ramp samples an infrared array flat is much richer in information than a CCD flat, and the temporal structure ($abcd$ indices) is the key to disentangling the various effects.

In what follows, we will simplify some expressions by writing $S_{ab}(i,j)\equiv S_{a}(i,j)-S_{b}(i,j)$. (Note the sign convention!)

The correlation functions satisfy the trivial properties:

• $\bullet$

$C_{abcd}(\Delta i,\Delta j)=0$ if $a=b$ or $c=d$.

• $\bullet$

$C_{abcd}(\Delta i,\Delta j)=C_{cdab}(-\Delta i,-\Delta j)$.

• $\bullet$

$C_{abcd}(\Delta i,\Delta j)=-C_{bacd}(\Delta i,\Delta j)=-C_{abdc}(\Delta i,\Delta j)=C_{badc}(\Delta i,\Delta j)$.

• $\bullet$

$C_{abcd}(\Delta i,\Delta j)=C_{aecd}(\Delta i,\Delta j)+C_{ebcd}(\Delta i,\Delta j)$ and $C_{abcd}(\Delta i,\Delta j)=C_{abcf}(\Delta i,\Delta j)+C_{abfd}(\Delta i,\Delta j)$.

The last property means that all of the correlation functions can be composed of “elementary” correlation functions $C_{a,a+1,c,c+1}(\Delta i,\Delta j)$.

Equation (12), like all covariance matrices, requires more than one realization to make a measurement. The standard approach, followed here, is to compare a pair of two flats. The “measured covariance” of any two observables ${\cal O}$ and ${\cal O}^{\prime}$ is then

[TABLE]

where the average is taken over pixels $(i,j)$ in the region of interest. The differencing removes small deviations such as imperfect illumination patterns, intrinsic variations in pixel area or QE, etc.

In a flat or dark field with many samples, we may also construct a correlation function averaged in the time direction. If we take the average of $n$ time-translations of the time windows $abcd$, then we find

[TABLE]

Since a flat field is not time-stationary (gain, non-linearity, and possibly other quantities will change as the voltage across the $p-n$ junction decreases), we must keep track of all the time indices $abcdn$ when fitting a model to the time-translation-averaged correlation function.

While in this study we consider individual time samples, space-based infrared surveys are often data rate limited and thus not every sample can be downlinked. Therefore future work should also examine how the BFE and other effects appear in the cross-correlation functions of flat field data with compression along the time axis, e.g., using the first few Legendre coefficients (Rauscher et al., 2019). If compression by linear combinations is used (Legendre coefficients, group averaging, etc.), any such cross-correlation function can be written trivially as an appropriate weighted sum over the $C_{abcd}(\Delta i,\Delta j)$.

2.4 Multiple exposures

When we discuss statistical algorithms, it will be essential to describe operations acting on multiple exposures. Here a specific exposure will be denoted with a $|$ separator, followed by the exposure type and number. For example, we write $S_{a}(i,j|{\rm F}_{k})$ to denote the signal in time step $a$ and pixel $(i,j)$ in the $k^{\rm th}$ flat field, and $S_{a}(i,j|{\rm D}_{k})$ for the $k^{\rm th}$ dark exposure. This formalism could be extended in the future to include other types of tests (besides flats and darks).

3 Theory

We now embark on the main calculation in this paper: the determination of the correlation function $C_{abcd}(\Delta i,\Delta j)$ including IPC, classical non-linearity, and the brighter-fatter effect to leading order. We will summarize the general result in Eq. (51). The measured correlation functions can then be used to simultaneously constrain the gain, IPC, and BFE parameters of a detector array or sub-region thereof – a step we will take on real data in Paper II. We will assume $a<b$ and $c<d$, since these functions contain all the information because of symmetries, but we do not assume anything else about the ordering. In particular, the exposure intervals $a...b$ and $c...d$ may be the same, may overlap, or may be non-overlapping. The “same interval” case $(a,b)=(c,d)$ will be the most familiar to readers who have worked with 2D image products before (e.g., CCD images, or CDS images from infrared arrays). However, the non-overlapping case turns out to be of particular use for measuring the sum of BFE and NL-IPC. It is the combination of these many different cases (itself possible due to the non-destructive read capability) that allows us to constrain so many parameters.

Rather than trying to solve everything at once, we begin this section by considering a “perfect” detector (no IPC, BFE, or any form of non-linearity; §3.1), and then adding layers of physical and mathematical complexity. In particular, we add linear inter-pixel capacitance, which simply introduces a convolution kernel (§3.2). Next is the classical non-linearity (§3.3), where covariances of higher-order moments of the charge appear; then we put the classical non-linearity together with the IPC (§3.4). We introduce the brighter-fatter effect in §3.5; there the buildup of charge in each pixel $Q(i,j;t)$ is a stochastic process with interactions between pixels, and we write and solve a differential equation for the moments of the charge. Non-linear IPC is briefly discussed and included in §3.6. The main result, Eq. (51), is presented in §3.7, and we discuss some special cases in §3.8. A reader not interested in the details of the derivation may skip directly to Eq. (51) and the simplifications in §3.8 (though one will still have to refer to the definitions for the quantities that summarize the time intervals: $t_{abcd}$, $T_{abcd}$, $\sigma_{abcd}$, and $\tau_{abcd}$, defined in Eqs. 18, 25, 46, and 47, respectively).

In this calculation, we suppose that the flat illumination provides current $I$ per pixel (units: e/s) and that the frame $a$ is saved at time $t_{a}$. We also assume that the flat field illumination uses a wavelength long enough for quantum yield effects to be insignificant (i.e., where one photon produces at most one electron-hole pair). At wavelengths blueward of the quantum yield threshold, it is possible for multiple carriers to be produced, and then (by diffusion) end up in separate wells, leading to an additional contribution to the flat field autocorrelation function (e.g., McCullough et al., 2008) as well as errors in gain determination.

3.1 Perfect detector

In a perfect detector, with $\alpha$, $\beta$, and ${a}$ all zero, each pixel operates independently. The mean charge accumulated in pixel $(i,j)$ in frame $a$ is

[TABLE]

and the covariance structure is

[TABLE]

Since the signal difference $S_{ab}(i,j)=g^{-1}[Q_{b}(i,j)-Q_{a}(i,j)]$, we then find a covariance structure:

[TABLE]

It is convenient then to define

[TABLE]

by inspection if $a<b$ and $c<d$, then $t_{abcd}$ is the amount of time in the intersection of the intervals $(t_{a},t_{b})\cap(t_{c},t_{d})$. In this case, we also have $t_{abab}=t_{ab}$.

Equation (17) is behind the usual concept of obtaining a system gain from a variance vs. mean plot: we have

[TABLE]

For a detector with no read noise, the ratio of variance to mean is then $1/g$. In practice $C_{abab}(0,0)$ contains a contribution from read noise, which can be removed by taking the slope of the variance vs. mean plot.

We now consider the various non-ideal detector effects. We begin by considering the effects one at a time, but we also need to consider interactions between the IPC and non-linearity, i.e. effects of order $\alpha\beta$ and $\alpha a$.

3.2 Inter-pixel capacitance

In the presence of IPC, the covariance structure of Eq. (17) is modified via smoothing by the IPC kernel. The IPC kernel acts locally in time, so we may write

[TABLE]

If the IPC kernel is represented by nearest-neighbor parameters $\alpha_{\rm H,V}$, and diagonal-neighbor parameters $\alpha_{\rm D}$, then

[TABLE]

(There are other non-zero terms.) The nearest-neighbor correlations are thus useful for measuring $\alpha_{\rm H}$ and $\alpha_{\rm V}$, and the diagonal-neighbor correlations for $\alpha_{\rm D}$.

Note that regardless of $K$, the IPC-induced correction to $C_{abcd}(\Delta i,\Delta j)$ remains proportional to $t_{abcd}$. Therefore, if IPC is the only non-ideal effect in the detector, the correlation function will be zero if $t_{abcd}=0$. The “disjoint correlation functions” with $t_{abcd}=0$ are therefore diagnostics of other effects – including, as we shall see, the brighter-fatter effect.

3.3 Classical non-linearity

The classical non-linearity – that arising from the nonlinearity of the electrons to data numbers conversion, Eq. (7) – contributes a correction to the correlation function that involves the third moment of the Poisson distribution. In the presence of only classical non-linearity, but no IPC or BFE, the pixels still operate independently, so for simplicity we will consider only one pixel. The connected skewness101010Connected skewnesses are defined by $\langle ABC\rangle_{\rm conn}=\langle\Delta A\Delta B\Delta C\rangle$, where $\Delta A\equiv A-\langle A\rangle$, etc. of charges at different times is

[TABLE]

since the connected skewness of a Poisson distribution is its mean, and all counts received after $t_{\min(a,b,c)}$ are independent of $Q_{\min(a,b,c)}$. This leads, after some algebra, to the ancillary result

[TABLE]

Then we find (suppressing $i$ and $j$ indices to avoid clutter):

[TABLE]

We will define

[TABLE]

(units: s2) so that

[TABLE]

Here the “$1-2\beta$” correction term is of little interest, since the correction is tiny even compared to WFIRST requirements – indeed, it represents the nonlinearity generated by a single electron, and if $\beta\sim{\cal O}(1)$ ppm/e, then this is a correction of order $10^{-6}$. The $T_{abcd}$ term can be much larger.

Note the following special cases of $T_{abcd}$:

• $\bullet$

If $a\leq b\leq c\leq d$, then we have $t_{abcd}=0$ and $T_{abcd}=t_{ab}t_{cd}\geq 0$.

• $\bullet$

If $a\leq c\leq b\leq d$, then we have $t_{abcd}=t_{bc}$ and $T_{abcd}=t_{ab}t_{cd}+(t_{b}+t_{c})t_{bc}\geq 0$.

• $\bullet$

If $a=c$ and $b=d$, then we have $t_{abcd}=t_{ab}$ and $T_{abcd}=2t_{b}t_{ab}$.

3.4 Interdependence of IPC and non-linearity

Because the IPC corrections to flat results are often large (e.g. $\alpha=1.25\%$ leads to an $\sim 10$% correction to the gain!) we need to consider the way in which IPC interacts with the non-linearity curve. This is particularly true given that IPC-non-linearity interactions affect both of the flat auto-correlation measurements of the BFE presented in this document. In particular, we want to capture the order $\alpha\beta$ terms in the flat auto-correlation function.

IPC and non-linearity may interact in a complicated way because the two steps do not in general commute. The approach taken here is the mathematical point of view: one chooses an ordering, and any additional effects – including issues associated with order of operations – are packaged into “non-linear IPC” (§3.6). The ordering we choose here is IPC first and then non-linearity (consistent with “standard” pipelines that treat non-linearity as the last step in the signal chain and thus the first correction implemented in data processing). From the physical point of view, the non-linear capacitance coming from the depletion region in the photodiode and the capacitive links between pixels should be thought of as a non-linear capacitor network that is solved simultaneously (the non-linearities coming from the rest of the signal chain and the analog-to-digital converter really do come later). Further consideration of this physical point of view is deferred to future work.

With these assumptions, the non-linearity interacts with the IPC according to

[TABLE]

where the approximation includes terms of order $\alpha_{\rm H,V}\beta$ but not $\alpha_{\rm D}\beta$ or $\alpha^{2}\beta$.

Our principal interest is in the contributions of order $\alpha\beta$ to the correlation function, which occur at either zero lag $(\Delta i,\Delta j)=(0,0)$ or for nearest-neighbor pixels, $(\Delta i,\Delta j)\in{\cal N}$. In general, the contribution of order $\alpha\beta$ to $C_{abcd}(\Delta i,\Delta j)$ (denoted below as $\Delta C_{abcd}(\Delta i,\Delta j)|_{\alpha\beta}$) has four parts: the covariance of the order $\alpha\beta$ term in $(i,j)$ with the order $1$ term in $(i+\Delta i,j+\Delta j)$ (which we will call the “$\alpha\beta\times 1$” term); the $\alpha\times\beta$ term; the $\beta\times\alpha$ term; and the $1\times\alpha\beta$ term. These can be read off from Eq. (27), and covariances can be computed using the fact that (i) the charges in each pixel are independent, and (ii) the Poisson statistics needed are in Eqs. (15), (16), and (23); this is an algebraically lengthy but straightforward exercise. The result is

[TABLE]

for zero lag,

[TABLE]

for the horizontal neighbors, and

[TABLE]

for the vertical nearest neighbors. (The order $\alpha_{\rm H}\beta$ and $\alpha_{\rm V}\beta$ contributions beyond the 4 nearest neighbor pixels are zero.) Note that we normally have $I^{2}T_{abcd}\gg It_{abcd}$, so that term is dominant.

3.5 Brighter-fatter effect; moving pixel boundaries

The effect of the BFE on pixel correlation functions in infrared arrays can be treated by considering the charge $Q(i,j;t)$ at time $t$ in pixel $(i,j)$ as a stochastic function. The new ingredient is that the charge deposited in each pixel between $t$ and $t+\delta t$ depends on the charge already present at time $t$. The BFE process is Markovian, in the sense that the state at $t+\delta t$ depends on the state at time $t$, but has no further dependence on the state of the system at any earlier times.111111If we also tried to include charge trapping effects such as persistence, then they would not be Markovian. These are not treated in the present formalism. This allows us to write the moments of the charge (here we consider the first two moments, the mean and covariance) at time $t+\delta t$ in terms of those at time $t$; by working to order $\delta t$, we can construct a system of differential equations for the moments, which we solve starting from the initial condition at $t=0$. This is analogous to what has been done for CCDs (see Coulton et al. 2018 for an implementation of a very similar method; see also Astier et al. 2019, who work directly with derivatives rather than very short but finite time steps). We extend this method to apply to unequal-time correlation functions by writing a similar set of differential equations for the covariance between a pixel at time $t_{1}$ and another pixel at time $t\geq t_{1}$, again with $t$ as the independent variable. We can take the covariance at $t_{1}$ as an initial condition.

Let us define the area defect of a pixel at time $t$ to be

[TABLE]

this is close to 1, with deviations controlled by the BFE. Then – given the state of the system $Q(i,j;t)$ at time $t$ – we can find the mean charge in pixel $(i,j)$ at time $t+\delta t$ as

[TABLE]

where $\delta t$ is taken to be small, and the subscript $|_{t}$ denotes that the state of the detector at time $t$ is fixed. Here $IW(i,j;t)\delta t$ is the probability that an electron is collected in pixel $(i,j)$ between times $t$ and $t+\delta t$; we assume $I\delta t\ll 1$ (one electron at a time), and will take the limit as $\delta t\rightarrow 0$ so that this approximation becomes arbitrarily good. Since we are turning the result into a first-order differential equation for each moment, we can drop terms of order $\delta t^{2}$ and higher in what follows. The change in 2nd moment is

[TABLE]

where we have expanded $Q(i,j;t+\Delta t)=Q(i,j;t)+\Delta Q(i,j;t)$, and the four terms on the right hand side correspond to the expectation values of $Q(i,j;t)Q(i^{\prime},j^{\prime};t)$, $\Delta Q(i,j;t)Q(i^{\prime},j^{\prime};t)$, $Q(i,j;t)\Delta Q(i^{\prime},j^{\prime};t)$, and $\Delta Q(i,j;t)\Delta Q(i^{\prime},j^{\prime};t)$ respectively. The last term is only non-zero if the two pixels are identical ($\delta_{ii^{\prime}}\delta_{jj^{\prime}}=1$), since then a single electron can increment both $Q(i,j)$ and $Q(i^{\prime},j^{\prime})$.

It is now possible to solve the above system of equations to first order in ${a}$. Let us first consider Eq. (32). Taking the average of the right-hand side over possible realizations at time $t$, we see that

[TABLE]

Recalling that ${\Sigma_{a}}=\sum_{\Delta i,\Delta j}{a}_{\Delta i,\Delta j}$, and using translation invariance to show that the $Q(i,j;t)$ all have the same expectation value, we see that

[TABLE]

This becomes a differential equation for $\langle Q(i,j;t)\rangle$:

[TABLE]

with solution starting from $\langle Q(i,j;t)\rangle=0$ at $t=0$:

[TABLE]

(The approximation holds to first order in the ${a}$ coefficients.)

The next step is to solve for the covariance matrix. This proceeds in two steps. First, one tracks the full second moment from time 0 to some later time $t_{1}$. Then one tracks a conditional second moment to a later time $t_{2}\geq t_{1}$. The mean of Eq. (LABEL:eq:Q2) is

[TABLE]

This can be turned into a differential equation. The first moment solution from Eq. (37) can be substituted in, and all second order terms in ${a}$ dropped:

[TABLE]

The initial condition is that $\left\langle Q(i,j;t)Q(i^{\prime},j^{\prime};t)\right\rangle=0$ at $t=0$. To solve Eq. (39) to first order in $a$, we use the standard method of first solving the equation at $a=0$ (the zeroth order solution), then substituting this into any term multiplying $a$ (or ${\Sigma_{a}}$) and solving again. This gives

[TABLE]

In our case, however, we need not just equal-time but also unequal-time correlation functions of the charge. This means we need to propagate the second moment at time $t_{1}$ to the covariance between times $t_{1}$ and $t_{2}$ (without loss of generality, $t_{2}>t_{1}$). Since the system is Markovian, if $t_{1}\leq t$ we can take the expectation value in Eq. (32) to be conditioned not only on the state of the system at $t$ but also at $t_{1}$. Then we can multiply Eq. (32) by $Q(i^{\prime},j^{\prime};t_{1})$. Since $Q(i^{\prime},j^{\prime};t_{1})$ is fully determined by the state of the system at $t_{1}$, we can pull it inside the expectation value:

[TABLE]

We next average this over states of the system at $t_{1}$ and $t$ (i.e., we remove the condition, but then get an expectation value on the right-hand side). Turning the result into a differential equation, we have

[TABLE]

with initial condition from Eq. (40), $\left\langle Q(i,j;t_{1})Q(i^{\prime},j^{\prime};t_{1})\right\rangle|_{{a}=0}=I^{2}t_{1}^{2}+It_{1}\delta_{ii^{\prime}}\delta_{jj^{\prime}}$. The solution to first order in $a$ is

[TABLE]

Subtracting out $\left\langle Q(i,j;t)\right\rangle\left\langle Q(i^{\prime},j^{\prime};t_{1})\right\rangle$ from Eq. (37) gives the covariance matrix:

[TABLE]

Recall that this is for $t\geq t_{1}$; for $t<t_{1}$, one can use the symmetry of the covariance matrix to obtain the result. If one considers only the linear response of the detector, this maps directly into the flat autocorrelation function:

[TABLE]

where we define the auxiliary quantities

[TABLE]

and

[TABLE]

Here $\sigma_{abcd}$ and $\tau_{abcd}$ have units of s2 and satisfy the following rules:

• $\bullet$

$\sigma_{abcd}=\sigma_{cdab}$ and $\tau_{abcd}=-\tau_{cdab}$.

• $\bullet$

If $a\leq b\leq c\leq d$, then $\sigma_{abcd}=0$ and $\tau_{abcd}=t_{ab}t_{cd}\geq 0$.

• $\bullet$

If $a=c\leq b=d$, then $\sigma_{abcd}=t_{ab}(t_{a}+t_{b})$ and $\tau_{abcd}=0$.

Note that in Eq. (45), $\tau_{abcd}$ describes the response of a correlation function to the odd part of ${a}$ while the response to the even part is described by $t_{ab}t_{cd}$. The response to the summed effect ${\Sigma_{a}}$ is encoded in $\sigma_{abcd}$.

In the presence of IPC, the BFE contribution to the correlation function should be convolved twice with the IPC kernel:

[TABLE]

3.6 Non-linear inter-pixel capacitance (NL-IPC)

The contribution of NL-IPC to the covariance of signals is, to order $K^{\prime}$,

[TABLE]

where in the first expression the first term comes from the order $K^{\prime}$ contribution to $S_{a}(i,j)$ and the second term from the contribution to $S_{c}(i+\Delta i,j+\Delta j)$. The final expression used the symmetry of $K^{\prime}_{\Delta i,\Delta j}$ under $(\Delta i,\Delta j)\rightarrow(-\Delta i,-\Delta j)$, and has defined the convolution $[KK^{\prime}]_{\Delta i,\Delta j}=\sum_{i_{1},j_{1}}K_{i_{1},j_{1}}K^{\prime}_{\Delta i-i_{1},\Delta j-j_{1}}$. The contribution to the correlation function is

[TABLE]

3.7 Combined correlation function

Putting together all of the combinations – IPC, non-linearity, BFE, NL-IPC, and the leading order interactions – we have the following expression, including corrections of order $\alpha$, $\alpha^{2}$, $\beta$, $\alpha\beta$, $a$, $\alpha a$, $\alpha^{\prime}$, and $\alpha\alpha^{\prime}$:

[TABLE]

where we have defined $[K^{2}]_{\Delta i,\Delta j}$ to be the auto-convolution of $K$, and $[K^{2}a]_{\Delta i,\Delta j}$ to be the convolution of $K^{2}$ and $a$.

3.8 Special cases used in detector characterization

We now turn our focus to the special cases that are used in detector characterization. We consider two special cases of the correlation function: the equal-interval correlation functions ($a=c<b=d$ – most similar to the auto-correlation that one would obtain from a CCD) and the non-overlapping correlation functions ($a<b<c<d$ – which exhibits new features only accessible with a non-destructive read capability). We also consider the mean-variance plot, which is a common diagnostic of the gain of a detector system, with a particular emphasis on how IPC, nonlinearity, and BFE affect the gain measurement.

In what follows, terms of order $\alpha$, $\alpha^{2}$, $\beta$, ${a}$, $\alpha\beta$, and $\alpha a$ are kept. Higher terms in the non-ideal detector effects are dropped.

3.8.1 Equal-interval correlation function

The case of $a=c<b=d$ corresponds to the auto-correlation of a single difference image $S_{a}-S_{b}$. It is therefore most comparable to what one would obtain with a CCD. The contributions at zero lag sum to:

[TABLE]

while the horizontal nearest neighbors are

[TABLE]

where we define ${a}_{\rm H}=({a}_{1,0}+{a}_{-1,0})/2$. A similar equation holds for the vertical nearest neighbors. For the diagonal neighbors, we have

[TABLE]

The equal-interval correlation function, especially but not exclusively at zero lag, contains a large contribution from read noise (from various sources), and this must be removed before interpreting it.

The time-translation-averaged versions of Eqs. (52–54) can be evaluated with straightforward algebra; they are

[TABLE]

for zero lag;

[TABLE]

for the nearest neighbor; and

[TABLE]

for the diagonal neighbor.

3.8.2 Non-overlapping correlation function

The case of $a<b<c<d$ is special, because then $\sigma_{abcd}=t_{abcd}=0$ and the correlation function – including contributions of IPC, classical non-linearity, and BFE – simplifies to

[TABLE]

That is, the non-overlapping correlation function is directly sensitive to the coefficients ${a}_{\Delta i,\Delta j}$, has no sensitivity to linear IPC at order $\alpha$, and only has sensitivity to the classical non-linearity $\beta$ at zero lag. At order $\alpha\beta$, there is a contribution in the nearest neighbors. There is a trivial mapping from the pixel-space lag in the correlation function $(\Delta i,\Delta j)$ to the lag in the BFE kernel ${a}_{\Delta i,\Delta j}$. Thus this should be a “clean” measurement of the inter-pixel non-linear effects (BFE and NL-IPC), insensitive to small errors in the determination of $I$ and $g$. Any source of noise that is uncorrelated across frames is also removed. The reset ($kTC$) noise is also removed, since the correlation function is constructed from correlated double samples. The main drawback is that the method is only sensitive to a combination of BFE and NL-IPC, and cannot distinguish between the two mechanisms.

The one large correction that is necessary is that ${a}_{0,0}$ must be corrected for the classical non-linearity (which is a larger effect than the BFE). Therefore we need to measure $\beta$ from the non-linearity ($t^{2}$ term) of the signal vs. time plot of the flat. Interestingly, this is sensitive to the combination $\beta-\frac{1}{2}{\Sigma_{a}}$. It follows that the flat non-linearity and the non-overlapping correlation functions contain an intrinsic degeneracy where $\beta$ and ${a}_{0,0}$ are both changed but holding the combination $\beta-\frac{1}{2}{a}_{0,0}$ constant. Other correlation functions are needed to break this degeneracy.

A secondary correction is that $C_{abcd}(\Delta i,\Delta j)$ in Eq. (58) is the correlation function of the signal, but the measurement contains signal+noise. Therefore any noise that is correlated across frames must be characterized and removed from $C_{abcd}(\Delta i,\Delta j)$.

3.8.3 Mean-variance slope

A common method to estimate the gain of a system is to determine the ratio of the mean signal in a pair of matched flats to the variance. In practice, since the measured variance contains read noise, one measures the slope of the variance as a function of the mean, e.g.:

[TABLE]

where $M_{ab}=\langle S_{a}(i,j)-S_{b}(i,j)\rangle$ and $V_{ab}=C_{abab}(0,0)$ is the variance of a difference frame. This construction only makes sense for $(a,b)\neq(c,d)$ (a common case is $a=c<b<d$). The mean is

[TABLE]

We obtain the variance from Eq. (52). The mean-variance slope is related to the gain by

[TABLE]

In the special case of $a=c<b<d$ (which will be used herein), we find

[TABLE]

Here, we have used that to order $\alpha\alpha^{\prime}$,

[TABLE]

In Eq. (62), the pre-factor is the traditional IPC correction to the gain. Following this is a non-linear correction term that depends on the “start time” $t_{a}$ of the measurement. Then comes a second non-linear correction term that depends on the “duration” $t_{ad}+t_{ab}$ of the measurement. Both are proportional to $\beta$ (or to the ${a}_{\Delta i,\Delta j}$); they have the same dependence in the special case of $t_{a}=t_{c}=0$. The last term is formally of order $\beta$, but is smaller than the previous two correction terms as it does not contain a factor of accumulated charge $It$.

4 Flat field simulations

We construct simulations for validation and interpretation. This simulated data set contains flats and darks that are designed to resemble the real data cubes from the Detector Characterization Laboratory (DCL) at NASA Goddard Space Flight Center, with an implementation of the key effects described in the earlier sections of this paper. The procedure consists of three main steps, visualized in a flowchart in Fig. 1: first, user inputs and specifications are read from a configuration file (top); second, charge is accumulated via draws from a Poisson distribution and modified by a BFE kernel, if the BFE is turned on (loop on lower left); third, all other effects including linear IPC, classical non-linearity, and noise are applied to each time step of the charge array, which is ultimately converted to a signal and stored in an output fits data cube (lower right). The remainder of this section delves into the specifics of how the simulated flat fields are constructed.

4.1 Details of simulation procedure

The first part of the script sets up the simulation that will be created by ingesting a configuration file and using defaults when selections are not specified by the user. The default settings create a datacube with dimensions of $4096^{2}$ pixels2 with 66 time samples with the bounding 4 rows and columns designated as reference pixels. Substeps set the total number of time slices at which the charge is computed between the stored time slices with the default set to substep=2 (for this default the computation is done for $2\times 66=112$ time steps). This setting exists to ensure convergence when the BFE mode is turned on. The user specifies quantities like gain $g$, current per pixel per second $I$, length of time sample in seconds, quantum efficiency QE, and IPC $\alpha$. Reset frames and reset levels can also be set in the config file.

After reading in user specifications and initializing arrays, charge is drawn and accumulated over the total time steps. For the initial time frame, a random realization of charge is drawn from a Poisson distribution with a mean of $QE\times I\times\delta t$, with $\delta t$ being the time between each time step. If the BFE is turned on, a matrix of pixel area defects $W(i,j;t)$ given by Eq. 31 is calculated by convolving a user-specified input kernel $a_{\Delta i\Delta j}$ with the charge distribution over the pixel grid at the given time $t$. Subsequent time frames compound the previous time frames with charge drawn from a Poisson distribution with the mean modified by the pixel area defects, i.e. charge is drawn from a Poisson distribution with a mean of $W\times QE\times I\times\delta t$. If the BFE is turned off, all time frames are accumulated with fully random realizations of charge.

After the charge has accumulated over all time frames, a linear IPC can be applied by convolving the full charge data cube with an IPC kernel. Note that from this stage onward, operations are only performed on the time samples that will be saved (i.e. in the default settings, IPC is applied to the 66 time samples and not the intermediate substeps). Non-linearity $\beta$ can also be applied after the IPC, where $Q(i,j;t)\rightarrow Q(i,j;t)-\beta[Q(i,j;t)]^{2}$.

We create realizations of noise datacubes using nghxrg121212https://github.com/BJRauscher/nghxrg, the HxRG Noise Generator written in Python by Bernard Rauscher (Rauscher, 2015). This software produces white read noise, pedestal drifts, correlated and uncorrelated pink noise, alternating column noise, and picture frame noise and was based on a principle components analysis of the James Webb Space Telescope NIRSpec detector subsystem. Here, we use input parameters tuned to WFIRST configurations.

The final step is to convert the charge into DN by dividing by $g$ and save in an array of unsigned 16-bit integers. The output datacube is saved in fits format with a header containing information about the configuration settings and parameter values used to run the simulation.

The main flat field generation is done as part of the solid-waffle pipeline, which is further described in Section 5; the noise file must be generated separately using the nghxrg package.

4.2 Test bed of simulated flats and darks

We created a set of 10 simulated flat fields and 10 simulated darks to test the characterization and BFE measurement framework presented in this paper. The input parameters were chosen to resemble the real detector data analyzed in Paper II and are summarized in the ‘truth’ column of Table 1. All simulated data cubes are ascending ramps where the signal level in DN increases as time increases and have NAXIS1=NAXIS2=4096 (spatial dimensions) and NAXIS3=66 (readout time frames). The time between each time frame is 2.75 s, $g=2.06$, QE=0.95, and $\beta=0.58$. The IPC kernel is as described in Eq. 9 so that $K_{0,0}=1-4\alpha$, $K_{0,\pm 1}=K_{\pm 1,0}=\alpha$, and all others are zero, with $\alpha=0.0169$. The BFE is turned on, and has a zero-lag component $a_{0,0}=-1.372$ ppm/e. We conservatively set substep=20 to ensure convergence. Table 1 provides symmetrized mean values of $[K^{2}a]_{\Delta i\Delta j}$, the convolution of $K^{2}$ (the auto-convolution of K) and $a_{\Delta i\Delta j}$. The simulated flats each have illumination $I=559$ e/s/pixel, while the simulated darks have illumination $I=0.191$ e/s/pixel, which was chosen so that the resulting slope of the signal vs time matched a typical real dark. Random seeds from 1001-1010 and 2001-2010 were set for the flats and darks, respectively.

Finally, we generated 10 noise data cubes using nghxrg with NAXIS dimensions as specified above, n_out=32 (number of detector outputs), nroh=8 (row overhead in pixels), rd_noise=4 (standard deviation of white read noise in e), pedestal=4 (pedestal drift in e), c_pink=3 (standard deviation of correlated pink noise), u_pink=1 (standard deviation of uncorrelated pink noise), c_ACN=1 (standard deviation of alternating column noise). For simplicity, we did not add any picture frame noise, and we set a bias offset of 19222 e to match a typical real dark. See Rauscher (2015) for further details of the noise recipe. Each of the 10 noise data cubes are combined with a flat and dark each so that each flat/dark pair have a noise realization in common; this simplification should not affect the results of this analysis in any significant way. These simulations were run on Pitzer (Ohio Supercomputer Center, 2018), a supercomputer at the Ohio Supercomputing Center (Ohio Supercomputer Center, 1987).

5 Characterization based on flat fields

We now turn to the practical problem at hand: extracting the calibration parameters ($g$, $\alpha$, $\beta$, $a_{\Delta i,\Delta j}$, etc.) from a suite of flat field and dark exposures. We first provide an overview of our characterization pipeline (known as solid-waffle), and then describe in detail the modules therein. The tools are written in Python 2, with data stored in numpy arrays. Due to the large file size associated with flat fields using multiple up-the-ramp samples (2.2 GB per file for a 66-frame H4RG flat), the full data set is not stored in RAM; instead the fitsio package was used to enable rapid access to small subsets of the data from disk without reading the entire file.131313With standard FITS routines and “usual astronomer writing Python” level of attention to data handling, reading the files can completely swamp the computation time!

Our analysis takes as input $N$ flat fields and $N$ dark images, where $N\geq 2$. The SCA is broken into a grid of $N_{x}\times N_{y}$ “super-pixels,” each of size $\Delta_{x}\times\Delta_{y}$ physical pixels. Statistical properties such as medians, variances, and correlation functions are understood to be computed in each super-pixel. Note that $N_{x}\Delta_{x}=N_{y}\Delta_{y}=4096$ for an H4RG (and 2048 for an H2RG). Super-pixels may be made larger to improve S/N, but this implies more averaging over the SCA so localized features and patterns may be washed out (we will see examples of this in Paper II).

Each super-pixel is processed through “basic” characterization. Following this, it passes through inter-pixel non-linearity (IPNL) determination using the non-overlapping correlation function, and then (optionally) through advanced characterization and other tests. We now describe these steps.

5.1 Basic characterization

The basic characterization step for a super-pixel is a prerequisite to studying all of the more subtle effects in the NIR detectors. It uses four time frames $t_{a}$, $t_{b}$, $t_{c}$, and $t_{d}$, and it does not take into account the diagonal IPC, the brighter-fatter effect, non-linear IPC, or signal-dependent QE.

We first construct the CDS images $S_{ab}(i,j|{\rm F}_{k})$ and $S_{ad}(i,j|{\rm F}_{k})$ within the range of column $i$ and row $j$ in the super-pixel, for each flat Fk. We build a median (over flats $k$) image $f(i,j)={\rm med}_{k=1}^{N}S_{ad}(i,j|{\rm F}_{k})$, and then a pixel mask based on requiring $f(i,j)$ to be within 10% of its median (this time taken over $i,j$). This rejects disconnected or low-response pixels.

Our next step is to perform a reference pixel subtraction. The procedure used here was obtained after some experimentation with DCL data, and is the default in our code, but may require some re-adjustment for other setups. We first find the range of rows $j_{\rm min}...j_{\rm max}=j_{\rm min}+\Delta_{y}-1$ corresponding to the super-pixel, and find the two $4\times\Delta_{y}$ blocks of reference pixels on the left and right sides of the SCA. For each flat exposure ${\rm F}_{k}$, and for each of our two CDS difference images ($S_{ab}$ and $S_{ad}$), we find the median of these $8\Delta_{y}$ pixels, and subtract this from the entire super-pixel. A similar procedure is applied to the dark images Dk. Note that this procedure only adds or subtracts a constant in the super-pixel, and does not correct each individual row.141414Correcting each row would print noise from the reference pixels as additional horizontal correlations. There are row-dependent drifts in the electronics, however we found that these are better eliminated at the correlation function level by either subtracting the correlation function in the darks or by the “baseline subtraction” method described in §5.2.

We next want to compute the raw gain, $\hat{g}_{abad}^{\rm raw}$. To do this, we need to compute the mean signal levels $M_{ab}$ and $M_{ad}$. The current default is to take the reference-corrected image $S_{ab}(i,j|{\rm F}_{k})$, and compute a mean in $k$ followed by a median in $(i,j)$. The variance $V_{ab}$ is obtained by taking each of the $N(N-1)/2$ pairs of flats $(k,\ell)$, with $1\leq k<\ell\leq N$. For each pair, we compute the difference $S_{ab}(i,j|{\rm F}_{k})-S_{ab}(i,j|{\rm F}_{\ell})$, and compute the inter-quartile range (IQR) of the $\Delta_{x}\Delta_{y}$ pixels.151515We use difference images because they are robust against permanent structure in the flat fields, e.g. variations in pixel area or quantum efficiency. The variance is estimated as $({\rm IQR}/1.349)^{2}/2$, as appropriate for a Gaussian (but note that the IQR estimator is robust against outliers, unlike the standard variance estimator), and with a factor of 2 to account for the fact that the flat difference has noise from both flats. The $V_{ab}$ used in the raw gain estimator is the average of the $N(N-1)/2$ estimates obtained from the various flat pairs. These means and variances are then plugged into Eq. (59).

Inter-pixel capacitance is addressed through the flat field auto-correlation method, which we implement as follows. For each of the $N(N-1)/2$ flat pairs, we construct a difference $T(i,j|{\rm F}_{k},{\rm F}_{\ell})=S_{ad}(i,j|{\rm F}_{k})-S_{ad}(i,j|{\rm F}_{\ell})$. We clip the top $100\epsilon$% and bottom $100\epsilon$% of the $T(i,j|{\rm F}_{k},{\rm F}_{\ell})$ map, leaving $100(1-2\epsilon)$% of the pixels unmasked. Then we define a horizontal correlation

[TABLE]

where the average is over pixels where both that pixel $(i,j)$ and its horizontal neighbor $(i+1,j)$ are unmasked. We then compute an averaged horizontal correlation

[TABLE]

Here we have subtracted the correlation from a pair of dark frames (to remove the contribution of correlated read noise), and averaged over the flat pairs. The factor of $\frac{1}{2}$ takes into account the fact that by subtracting two flats, we have doubled the correlation function. Finally, the factor of $f_{\rm corr}$ takes into account the suppression of correlations by the histogram clipping of $T$. It depends on $\epsilon$; for our default choice of $\epsilon=0.01$, we have $f_{\rm corr}=0.7629$. See Appendix A for a derivation of $f_{\rm corr}$. A similar calculation is used to obtain the vertical correlation function $C_{\rm V}$ and the diagonal correlation function $C_{\rm D}$.

Finally, we need a measure of ramp curvature. We construct the difference box

[TABLE]

and perform the usual reference pixel subtraction (based on all $8\Delta_{y}$ “left+right” reference pixels in the same range of rows as the super-pixel). We clip the pixels corresponding to the top and bottom $100\epsilon$% of the histogram of $S_{ab}(i,j|{\rm F}_{k})$ and of $R(i,j|{\rm F}_{k})$, and then compute

[TABLE]

Note that $1+{\tt frac\_dslope}$ is the ratio of the slope of the signal (in DN/frame) in the $cd$ interval relative to the $ab$ interval. For a perfectly linear detector, ${\tt frac\_dslope}$ should be zero. For a non-linear detector, the mean signal is $\langle S_{a}(i,j)\rangle=[It_{a}-\beta_{\rm r}(It_{a})^{2}]/g$, where

[TABLE]

is the ramp curvature (here ${\Sigma_{a}}$ denotes the signal-dependent QE, and enters via Eq. 37). Via straightforward algebra, we can see that the slope difference ratio ${\tt frac\_dslope}$ is expected to be $-\beta I(t_{c}+t_{d}-t_{a}-t_{b})$.

We may now construct an IPC + non-linearity corrected ($\alpha\beta$-corrected) gain $g$, estimated current per pixel $I$, horizontal and vertical IPC $\alpha_{\rm H}$ and $\alpha_{\rm V}$, and ramp curvature $\beta_{\rm r}$ by iteratively solving the system of equations:

[TABLE]

This is 5 equations for 5 unknowns; note that the difference between $\beta$ and $\beta_{\rm r}$ (i.e., the signal-dependent QE term ${\Sigma_{a}}$) has been neglected in the gain and IPC determination. Initializing the system with $g=\hat{g}^{\rm raw}_{abad}$, $\alpha_{\rm H}=\alpha_{\rm V}=\beta=0$, $I=gM_{ad}/t_{ad}$, and solving the above equations in turn for $g$, $\alpha_{\rm H}$, $\alpha_{\rm V}$, $I$, and $\beta$ leads to rapid convergence.

The resulting parameters $g$, $\alpha_{\rm H}$, $\alpha_{\rm V}$, $I$, and $\beta$ contain small residual biases due to the BFE, nonlinear IPC, and signal-dependent QE if these phenomena are present. These will be explored in more detail in Paper II.

5.2 IPNL determination via the non-overlapping correlation function

With the basic parameters in each super-pixel measured, we may now measure the non-overlapping correlation function, $C_{abcd}(\Delta i,\Delta j)$ for $a<b<c<d$. This is almost a direct test for the presence of inter-pixel non-linearities (BFE and NL-IPC), since it contains no contribution from linear IPC, and only small corrections for classical non-linearity ($\beta$) are required. In particular, at zero lag, Eq. (58) can be rearranged to give

[TABLE]

The ramp curvature does not yield an estimate directly for $\beta$, but rather the combination $\beta_{\rm ramp}=\beta-\frac{1}{2}{\Sigma_{a}}$. We also recall that to order $\alpha a$, we have $[K^{2}a^{\prime}]_{0,0}=[K^{2}a]_{0,0}-(1-8\alpha){\Sigma_{a}}$. We can thus write, to ${\cal O}({\alpha a})$:

[TABLE]

Similarly, one may compute the adjacent pixel correlation functions:

[TABLE]

and similarly for the vertical directions.

Equations (71) and (72) show that the non-overlapping correlation function method, as we have implemented it, is sensitive to the $[K^{2}a^{\prime}]_{\Delta i,\Delta j}+[KK^{\prime}]_{\Delta i,\Delta j}$ coefficients. Note that the BFE and NL-IPC appear together, both with $t_{ab}t_{cd}$ time dependence, and the non-overlapping correlation function method provides no way to separate them. This method has only a small correction on the right-hand side due to the ramp curvature $\beta_{\rm r}$, so this method of IPNL determination is not subject to spurious detection due to small errors in the basic parameters ($g$, $I$, and $\alpha$). In most practical situations, we will find that the $4\alpha_{\rm H}\beta_{\rm r}$ correction is smaller than the IPNL, and the $2(1-8\alpha)\beta_{\rm r}$ correction is similar to the IPNL (see Paper II for quantitative details on a WFIRST development detector).

Our pipeline provides results out to a separation of 2 pixels in either the horizontal or vertical directions, i.e., it reports a $5\times 5$ kernel $[K^{2}a^{\prime}+KK^{\prime}]$.

We now turn to the implementation details of $C_{abcd}(\Delta i,\Delta j)$ in the pipeline itself. The correlation function can be determined by the same methods used to compute $C_{\rm H}$ and $C_{\rm V}$. However, we found in initial studies on DCL data that the measurements showed statistically significant deviations depending on which flat was used, which are suspected to be low frequency noise in the data (see horizontal stripes in the dark image and discussion in Paper II). Therefore, the default setting in our pipeline is to filter out the low frequencies via a baseline correction: instead of using the raw correlation function,

[TABLE]

(where the overbar denotes an average and $N_{\rm pair}$ is the number of pixel pairs in the sum), we find a “baseline” contribution:

[TABLE]

where the pair summation runs over $6\leq|\Delta i^{\prime}|\leq 10$, and again $N^{\prime}_{\rm pair}$ is the number of pixel pairs in the sum. That is, the baseline is the correlation function obtained by replacing pixel $(i+\Delta i,j+\Delta j)$ with the average of pixels in the same row but 6–10 pixels left or right (ahead or behind in the readout sequence). Both the leading and trailing regions are used with equal weight, except that (i) the standard 1% outlier rejection is used before taking the covariance, and (ii) the implementation in the code rejects one of these regions if the pixel pair $(i,j)\leftrightarrow(i+\Delta i^{\prime},j+\Delta j)$ would span an output channel boundary. The correction regions are shown schematically in Figure 2. We then define a corrected correlation function:

[TABLE]

5.3 Advanced characterization

While the basic characterization stage is sufficient to provide a pixel mask and the properties (gain, IPC, non-linearity) needed to convert the non-overlapping correlation function to an IPNL measurement, there are several ways it could be improved. The statistical uncertainties in the gain and IPC are significant, especially with small super-pixels. Moreover, if the BFE exists in these detectors (and we will see in Paper II that it does), then it imprints a bias on $g$, $\alpha$, etc., and an iterative process is required to de-bias the final result. The “advanced characterization” tool handles both of these issues.

To motivate our approach to the first issue (noise in the parameters), and understand the improvement in knowledge of gain and IPC that can be achieved, let us first recall the uncertainty in gain and IPC achievable by the “basic” approach. If $t_{ab}\ll t_{ad}$, then in the computation of $\hat{g}_{abad}^{\rm raw}$ from Eq. (59), the uncertainty is dominated by $V_{ad}$. The variance of a Gaussian distribution with $n_{\rm pix}=\Delta_{x}\Delta_{y}$ samples has a fractional uncertainty of $\sqrt{2/n_{\rm pix}}$. Similarly, the correlation coefficient $\rho\sim 2\alpha$ of two adjacent pixels has $2n_{\rm pix}$ samples (counting both vertical and horizontal pairs) and hence an uncertainty of $\sqrt{1/(2n_{\rm pix})}$. With $N-1$ flat pairs, we should thus in principle achieve

[TABLE]

In practice, our pipeline does not do this quite well – the uncertainty in $V_{ab}$ is not negligible, and the use of the IQR carries a factor of 1.64 penalty in error for a Gaussian relative to the “idealized” case.161616See, e.g., DasGupta (2011), §9.5 for a general discussion of this issue. However, if the flat field has $N_{\rm frame}$ samples, and we break it into “sub-flats” of length $\mu$, one might expect that by combining the sub-flats we could achieve an uncertainty that is reduced by a factor of $\sqrt{N_{\rm frame}/\mu}$. Our pipeline does not quite achieve this, but it nevertheless can beat the estimate in Eq. (76). One expects that if $\mu$ is decreased, we should see a reduction in the error (down to the fundamental limit of $\mu=1$). However, the magnitudes of the correlation functions decrease as one decreases $\mu$, and hence we become more sensitive to the subtraction of noise from $C_{\rm H}$ and $C_{\rm V}$. Therefore there is a trade-off in the choice of $\mu$ (and $\mu^{\prime}$ defined below), and we allow the user to set these in the configuration file.

The implementation of these ideas in our pipeline is as follows. First, the user sets the range of frames used (earliest frame $a$ and latest frame $d$), as well as two integers $\mu$ and $\mu^{\prime}$ (with $\mu^{\prime}<\mu$) in the configuration file describing the spacing of time slices used in the gain and IPC determination; typical values would be $\mu^{\prime}=1$ and $\mu=3$. We then compute an averaged correlation function

[TABLE]

where $C_{{\rm H},a+j,a+j+\mu}$ is obtained using the same methodology as in basic characterization using the difference image of frames $a+j$ and $a+j+\mu$. Something similar is performed to compute $\bar{C}_{\rm V}$. Finally, we compute the difference of variances

[TABLE]

where $V_{ef}$ is the variance of the difference of frames $e$ and $f$ as obtained using the same methodology as in basic characterization.

In the advanced characterization stage, the mean information on the ramp is obtained by taking the sequence of differences $M_{a,a+1}$, $M_{a+1,a+2}$, … $M_{d-1,d}$, and performing a linear fit:

[TABLE]

where the sum of the square of residuals is minimized. One then wants to simultaneously solve the equations:

[TABLE]

Here “Err[…]” denotes the contribution to the specified quantity coming from BFE, NL-IPC, and signal-dependent QE (we will consider these shortly; in future versions of the pipeline we may add other effects). Once again, these are 6 equations for 6 unknowns ($g$, $I$, $\alpha_{\rm H}$, $\alpha_{\rm V}$, $\alpha_{\rm D}$, and $\beta_{\rm r}$). A straightforward and effective method is to alternately use the $\Delta\bar{V}$, $c_{1}$, and $c_{0}$ equations to solve algebraically for $I$, $g$, and $\beta_{\rm r}$; and then to use the $\bar{C}_{\rm H}$, $\bar{C}_{\rm V}$, and $\bar{C}_{\rm D}$ equations to solve for $\alpha_{\rm H}$, $\alpha_{\rm V}$, and $\alpha_{\rm D}$.

The advanced characterization pipeline can run in two modes for computing the error terms Err[…]; these are none, bfe, and nlipc. The none mode is the simplest: it sets the error terms to zero. When run on a detector that has, e.g., the BFE, the none mode is subject to similar biases as the “basic” characterization, but can give smaller statistical error.

Given that we will see in Paper II that the BFE is significant for the H4RGs, we included the bfe mode. This computes the error terms Err[…] under the assumption that there is a BFE ($a_{\Delta i,\Delta j}\neq 0$), but with no non-linear IPC ($K^{\prime}_{\Delta i,\Delta j}=0$) or signal-dependent QE (${\Sigma_{a}}=0$). Under these assumptions:

[TABLE]

One must iteratively perform the advanced characterization computation in this section and solve for the $[K^{2}a^{\prime}+KK^{\prime}]$ kernel via the procedure in §5.2 until all parameters are converged.

A similar approach is used for the nlipc mode, where the IPNL kernel is attributed entirely to NL-IPC instead of the BFE. In this case:

[TABLE]

5.4 Characterization of simulated detector data

In Fig. 3, we show the results of applying the aforementioned advanced characterization steps to pairs of simulated flats and darks using the specifications described in § 4. Mean quantities over $N_{good}$ good super-pixels and their statistical errors are provided in Table 1. The latter values are computed as standard deviations on the mean of the $N_{good}$ super-pixels.

Table 1 contains the values of the recovered BFE coefficients obtained after iterative application of the advanced characterization described in this section and the method described in § 5.2 (labeled ‘Method 1’). The time frames used for our fiducial scheme are 3, 11, 13, and 21. solid-waffle solves for $[K^{2}a^{\prime}+KK^{\prime}]$, which reduces to $[K^{2}a^{\prime}]$ since $K^{\prime}=0$ in the simulations. $[K^{2}a^{\prime}]$ values are provided as symmetrical averages for stacks of 3 and 10 simulated flats and compared against the simulation input, where the input $a^{\prime}$ has been convolved with the input $K^{2}$ (auto-convolution of $K$) to get values comparable to what is actually measured in the correlation analysis. In the central value at zero-lag, $[K^{2}a^{\prime}]_{0,0}$, we can see there is a bias of 0.1398 ppm/e for the 10 flat stack relative to the input into the simulation (12.1% bias compared to the input value). We compute the Method 1 BFE coefficients for two alternative time intervals; the first uses time intervals of half the duration of the fiducial scheme and results in a bias of 11.2% in the zero-lag coefficient, while the second uses time intervals of twice the fiducial duration and results in a bias of 20.5%. We note that the changes to $\beta_{\rm ramp}$ in these alternative time setups are much less than a percent.

We have also run the simulation with only BFE and no IPC and no classical non-linearity. In the fiducial 3, 11, 13, 21 time frame analysis setup for 10 simulated flats and darks, we obtain $[K^{2}a^{\prime}]_{0,0}=-1.3225\pm 0.0077$ (stat) ppm/e, which is biased compared to the input value of -1.3720 ppm/e by 3.6%. In this setup, the correct charge per time slice, gain, $\alpha$ and $\beta$ are consistent with the input values (where the latter two are consistent with 0). We suggest the likely source of bias in the BFE coefficients extracted from the simulations is due to exclusion of higher order terms in the interactions among the BFE, IPC, and classical non-linearity, and we will revisit this investigation in future work. Note that such an investigation of higher-order effects has recently been completed for CCDs (Astier et al., 2019).

5.5 Raw gain and equal-interval correlation tests

The above techniques enable us to correct the measured properties (gain, IPC, and non-linearity) for the IPNL – if we know whether to interpret the non-overlapping correlation function as BFE, NL-IPC, or a mixture of the two. Fortunately, the flat field auto-correlations carry enough information to distinguish the sources of IPNL. We cannot do this based on the non-overlapping correlation function, since in that case both BFE and NL-IPC scale as $\propto t_{ab}t_{cd}$, but we can use the scalings of the raw gain $\hat{g}^{\rm raw}_{abad}$ and the adjacent-pixel correlations $C_{adad}(\langle 1,0\rangle)$ as a function of which intervals in the flat field are taken.

5.5.1 Raw gain vs. interval duration

In this case, the key observable is the mean-variance slope, in the form $\hat{g}^{\rm raw}_{abad}$. From Eq. (62), one sees that there should be two time dependences: one that depends on the start time $t_{a}$ and contains only the classical non-linearity $\beta$, and one that depends on the duration pattern ($t_{ab}$ and $t_{ad}$) and depends on both $\beta$ and ${a}_{0,0}$. In this section, we consider the first dependence. We fix $t_{a}$ and fit a linear equation of the form:

[TABLE]

where $C_{0}$ is the intercept and $C_{1}$ is the slope.171717An alternative, which we tried first, is to do a linear fit $\hat{g}^{\rm raw}_{abd}=B_{0}+B_{1}I(t_{ad}+t_{ab})$, and use the slope-to-intercept ratio $B_{1}/B_{0}$. This procedure is not stable because the intercept $B_{0}$ is obtained by extrapolating to $t_{ad}+t_{ab}=0$. There is therefore a strong anti-correlation between the slope and intercept, which results in a noise bias: $B_{1}/B_{0}$ is biased upward by an amount $-{\rm Cov}(B_{0},B_{1})/B_{0}^{2}$. The amount of bias increases as subsets of the data are used. The formulation of Eq. (83) avoids this problem. From Eq. (62), we interpret the slope as

[TABLE]

where the three possibilities on the right are for no IPNL (none), and for the cases where the IPNL is pure BFE (bfe) or pure NL-IPC (nlipc). If there is a measurement of $[K^{2}a^{\prime}+KK^{\prime}]$ from the non-overlapping correlation function, then Eq. (84) can be used to test these hypotheses about its origin.

We compute the raw gain for frame triplets from [1,3,5], [1,3,6],…,[1,5,18] as a function of the signal level accumulated between the first time slice and the time slice $d=5...18$ for the simulated detector data. The top panel of Figure 4 visualizes the results of this test. Each data point is a mean over all super-pixels, with an error bar based on the error on the mean. The dashed line is the bfe interpretation of quantities from Method 1, as given by Eq. 84, and the solid line is the nlipc interpretation. These lines are plotted such that the central values pass through the center of the measurements. The simulated data agree firmly with the bfe slope, as is expected.

In each panel of Figure 4, we also show a systematic error related to the modeling of the non-linearity (“sys nl”). This is based on fitting a 5th order polynomial to the median signal levels in the detector. For both this 5th order curve and the quadratic ($\beta$) model, we computed the expected raw logarithmic gain $\ln g^{\rm raw}_{a,b,d}$ for Poisson statistics, compute the difference, and plot an error bar showing the peak$-$valley range. For the case of these simulated data this systematic is negligible, however we include it in anticipation of the analysis of the real data in Paper II where there may be deviations of the classical non-linearity from the $\beta$ model.

We can also make an estimate of the zero-lag BFE coefficient by re-arranging the left part of Eq. 84 and substituting $\beta_{\rm r}=\beta-\frac{1}{2}{\Sigma_{a}}$:

[TABLE]

Since we did not include non-linear IPC in the simulations, Eq. 85 simplifies to $a_{0,0}+8\alpha a_{<1,0>}$. For the 10 simulated flats, $\hat{a}_{0,0,M2}=-1.3513\pm 0.0144$ (stat) ppm/e. The input value is -1.3341 ppm/e, so these values agree to within 1.2$\sigma$.

5.5.2 Raw gain vs. interval center

A similar test can be carried out by measuring how the raw gain $\hat{g}^{\rm raw}_{abad}$ varies with $t_{a}$ as $t_{ab}$ and $t_{ad}$ are held fixed. We fit:

[TABLE]

In this case, we see that one should have

[TABLE]

Note that the slope $C^{\prime}_{1}$ has no sensitivity to the BFE – the none and bfe cases give identical predictions. It is however sensitive to NL-IPC.

We compute the raw gain for frame triplets from [1,3,5], [2,4,6],$...$,[14,16,18], as a function of the signal level accumulated between the first time slice and the time slice $a=1,...,14$ for the simulated detector data. The middle panel of Figure 4 visualizes the results of this test, showing that the simulated data are again consistent with the bfe slope.

Re-writing Eq. 87 and using the fact that $\alpha^{\prime}=0$, we can also compute $\beta=\frac{1}{2}C_{1}^{\prime}$ and ${\Sigma_{a}}=C_{1}^{\prime}-2\beta_{\rm r}$. $\beta=0.5677\pm 0.0079$ ppm/e, which is very close to the input value of 0.58 ppm/e (within 1.6$\sigma$). Likewise, ${\Sigma_{a}}=-0.0211\pm 0.0158$ is very close to the expected value of 0.

5.5.3 CDS autocorrelation vs. signal

This method uses the equal-interval correlation function in adjacent pixels, Eq. (53). Once the preliminary characterization of the detector has been performed, we may fix the starting time $t_{a}$ and fit the combination $g^{2}C_{abab}(\pm 1,0)/(It_{ab})$ as a function of $t_{ab}$, i.e., we fit

[TABLE]

The slope is given by

[TABLE]

Adding $8\alpha_{\rm H}\beta_{\rm r}$ to the left hand part of Eq. 89 gives

[TABLE]

We measure the IPC via basic characterization of frame triplets from [1,2,3], [1,2,4],…, [1,2,18], and CDS auto-correlations for [frame 3 - frame 1], [frame 4 - frame 1],…, [frame 18 - frame 1]. The bottom panel of Figure 4 visualizes the results of this test on the simulated detector data, which are consistent with the bfe interpretation.

We expect that $[K^{2}a^{\prime}+2KK^{\prime}]_{\langle 1,0\rangle}-\alpha{\Sigma_{a}}$ simplifies to $[K^{2}a^{\prime}]_{\langle 1,0\rangle}$ for the simulated data. This value is $0.1855\pm 0.0032$ ppm/e and can be compared with the value obtained from Method 1 of $[K^{2}a^{\prime}]_{\langle 1,0\rangle}=0.1980\pm 0.0033$ and the input value of 0.2034 ppm/e ($\sim$9% difference between the input and the value obtained with the CDS autocorrelation method).

6 Discussion

In this paper, we present formalism to connect flat field correlations to various detector effects in infrared detector arrays, including non-linear effects such as the BFE and NL-IPC. This formalism is built up through first considering the Poisson statistics in a perfect detector and then including contributions from the IPC kernel, classical non-linearity, BFE, and NL-IPC. In the expression for the combined cross-correlation of two CDS images (sampled at time frames $a,b,c,d$), we consider the leading order interactions, namely $\alpha$, $\alpha^{2}$, $\beta$, $\alpha\beta$, $a$, $\alpha a$, $\alpha^{\prime}$, and $\alpha\alpha^{\prime}$. We discuss two special cases of the combined correlation function: the non-overlapping correlation function ($a<b<c<d$), which has the most sensitivity to the inter-pixel non-linear effects, but cannot by itself distinguish between the BFE and NL-IPC; the equal-interval correlation function ($a=c<b=d$), which is the auto-correlation of a CDS image and is most similar to the flat field statistics available for CCDs. We also discuss features of the raw gain for the case of ($a=c<b<d$), which provides a means of distinguishing between the BFE and NL-IPC interpretations through the different behaviors of these mechanisms as a function of time.

We describe a procedure for characterizing detector arrays and extracting measurements of the IPNL. This involves constructing CDS images, performing a reference pixel subtraction, computing the raw gain, IPC, correlations (in the horizontal, vertical, and diagonal directions), and ramp curvature; we use these to solve for $g$, $\alpha_{\rm H}$, $\alpha_{\rm V}$, $I$, and $\beta$. We show how to use the non-overlapping correlation function to obtain the IPNL and also how to apply an iterative scheme to correct the $g$, $\alpha_{\rm H}$, $\alpha_{\rm V}$, $I$, and $\beta$ for residual biases imprinted by the IPNL.

We validate our methodology on simulated flat fields, which are constructed to imitate characteristics ($g$, $\alpha_{\rm H}$, $\alpha_{\rm V}$, $I$, and $\beta$) of the real detector array tested in Paper II. For this first investigation, we input a BFE kernel (but no NL-IPC). We extract parameters that match the inputs with high accuracy, except for the BFE kernel, for which we obtain a zero-lag component which is biased by 12%. We also show that the raw gain and equal-interval correlation function interpretation tests are successful in distinguishing between the BFE and NL-IPC as the underlying mechanism for the IPNL in the simulations. Given the success in obtaining equivalent inputs and outputs of the other key parameters, namely $\beta$ and $\alpha$, we suggest the 12% bias in the extracted BFE kernel could likely be explained by unaccounted interactions at higher orders that were dropped in the approximations used in this work. The impact of these higher order terms is under investigation and will be addressed in future work.

Acknowledgements

We thank the Detector Characterization Laboratory personnel, Yiting Wen, Bob Hill, and Bernie Rauscher at NASA Goddard Space Flight Center for their efforts enabling the existence and access to the data analyzed in this series of papers, and we thank Chaz Shapiro, Andrés Plazas, and Eric Huff for helpful discussions. We thank Jay Anderson and Arielle Bertrou-Cantou for useful presentations to the Detector Working Group on their analyses of non-linearities in the HST/WFC3-IR and Euclid H2RG detectors. We thank the anonymous referee for helpful suggestions that improved the clarity of this paper. We are also grateful for the use of Ohio Supercomputer Center (1987) for computing the results in this work. AC and CMH acknowledge support from NASA grant 15-WFIRST15-0008. During the preparation of this work, CMH has also been supported by the Simons Foundation and the US Department of Energy. Software: Astropy (Astropy Collaboration et al., 2013, 2018), fitsio (Sheldon, 2019), Matplotlib (Hunter, 2007), NumPy (Oliphant, 2006–), SciPy (Jones et al., 2001–)

Appendix A Clipping correction to the covariance

This appendix considers the correction to the covariance matrix of two jointly Gaussian distributed variables, $X$ and $Y$, when those distributions are clipped. We are interested in the parameter $f_{\rm corr}$ defined by

[TABLE]

We assume that a fraction $\epsilon$ of the data are clipped from both the top and the bottom of the distribution in $X$ and $Y$; if $X$ and $Y$ were independent, this would mean that a fraction $(1-2\epsilon)^{2}$ of the data points survive the clipping, but the fraction that survives may be larger if $X$ and $Y$ are covariant.

The determination of $f_{\rm corr}$ is invariant to linear rescaling of $X$ and $Y$, so without loss of generality, we assume that $X$ and $Y$ both have mean 0 and variance 1. Their “true” covariance is then the correlation coefficient $\rho$. The clipping is equivalent to the restriction of the data at $|X|,|Y|<\xi$, where

[TABLE]

Since the clipped distribution still has $\langle X\rangle=\langle Y\rangle=0$ by symmetry, we are interested in the mean value of $XY$:

[TABLE]

where the denominator is the survival probability of a data point $(X,Y)$, and the probability distribution is

[TABLE]

The covariance can be expanded in a power series in $\rho$; the leading term is

[TABLE]

so that

[TABLE]

The clipped covariances used in this paper to measure IPC are corrected using the leading constant term in $f_{\rm corr}$. The correction factor should converge to 1 as $\epsilon\rightarrow 0$; this is easily verified.

Note that the “correction” is not small: for $\epsilon=0.01$ (i.e. clipping the top 1% and bottom 1% of the distribution) we have $f_{\rm corr}=0.7629$. If one clips more of the distribution, the correction becomes enormous: at $\epsilon=0.025$ we have $f_{\rm corr}=0.5758$. Going the other way, even for $\epsilon=10^{-3}$, the correction is $f_{\rm cor}=0.9587$.