Mitigation of the Brighter-fatter Effect in the LSST Camera

We parameterize the BFE from flat-field statistics as changes in effective pixel area from distortions of a rectilinear grid, as proposed by Astier et al. (2019). We derive the pixel area distortions from the covariance function evaluated from flat fields. For covariances at a given signal level μ, both expressed in ADU:

$$\begin{eqnarray}&&\begin{array}{l}{C}_{{ij}}(\mu )=\displaystyle \frac{\mu }{g}\left[{\delta }_{i0}{\delta }_{j0}+{a}_{{ij}}\mu g+\displaystyle \frac{2}{3}{[{\boldsymbol{a}}\otimes {\boldsymbol{a}}+{\boldsymbol{a}}{\boldsymbol{b}}]}_{{ij}}{\left(\mu g\right)}^{2}\right.\\ \quad +\,\left.\displaystyle \frac{1}{6}{[2{\boldsymbol{a}}\otimes {\boldsymbol{a}}\otimes {\boldsymbol{a}}+5{\boldsymbol{a}}\otimes {\boldsymbol{a}}{\boldsymbol{b}}]}_{{ij}}{\left(\mu g\right)}^{3}+\cdots \right]+{n}_{{ij}}/{g}^{2},\end{array}\end{eqnarray} \tag{ 2 }$$

where a_ij [1/el] describes the fractional change in pixel area due to changes in the pixel boundaries from accumulated source charges, b_ij [1/el] describes smaller time-dependent processes that weaken or strengthen the BFE as we will explore in a later section, μ is a given mean signal over the image region, g is the sensor gain [el ADU⁻¹], n_ij is a constant term that represents the mean square noise (n₀₀ ∼ 5.7 el², where the non-(0, 0) terms of the noise matrix are contributed by the electronics and overscan), a b is an element-wise matrix multiplication, and ⨂ refers to a discrete convolution operation.

This expansion allows for higher-order terms which are nonlinear BFEs. Astier et al. (2019) found this model to fit their data well, and the contribution from the nonlinear terms was not negligible (with the μ² term accounting for 20% of the loss in variance). Any physically well-motivated correction algorithm would therefore need to take higher-order terms into account.

Coulton et al. (2018) assumed that the deflection field along a single boundary is the gradient of a unitless, 2D scalar kernel K. In analogy with electrostatics, the effect of the kernel on the covariance matrix then takes the form of Poisson's equation with the Laplacian of the kernel (Equation (19) of Coulton et al. 2018):

$$\begin{eqnarray}&&{\widetilde{C}}_{{ij}}=-{\mu }^{2}{{\rm{\nabla }}}^{2}{K}_{{ij}}\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&{\widetilde{C}}_{{ij}}={C}_{{ij}}-\left(\displaystyle \frac{\mu }{g}{\delta }_{i0}{\delta }_{j0}+\displaystyle \frac{{n}_{{ij}}}{{g}^{2}}\right).\end{eqnarray} \tag{ 4 }$$

The true image can thus be recovered from the measured image using Equation (22) of Coulton et al. (2018):

$$\begin{eqnarray}&&\hat{F}({\boldsymbol{x}})=F({\boldsymbol{x}})+\tfrac{1}{2}{\rm{\nabla }}\cdot \left[F({\boldsymbol{x}}){\rm{\nabla }}\left(K\otimes F({\boldsymbol{x}})\right)\right]\end{eqnarray} \tag{ 5 }$$

where $\hat{F}$ is the true image and F is the measured image, and x runs over the pixel space in two-dimensions. The expression includes a factor of 1/2 to average the BFE over the duration of the exposure as the first charge to enter the pixel experiences no BFE and the last charge is assumed to experience twice the average effect.

The correction explicitly assumes that the covariance matrix does not evolve as a function of flux (Coulton et al. 2018), which can be seen by the form of Equation (2). However, Astier et al. (2019) showed their functional form up to O(μ⁴) fits the observed PTC well over the flux range they considered, citing χ²/N_dof = 1.23. The number of observations is the number of flat pairs, and the parameters for just the variance element of the fit are g, n₀₀, and a₀₀, so they use N_dof = N(flat field pairs) − 3. Their good model fit suggests that the PTC shape is contributed by other nonlinear terms included in the covariance model. Whether to include the higher-order BFEs in the kernel, or to somehow weight the correction toward higher signal levels (where the BFE correction will be most important) becomes ambiguous, as there is no special signal level at which to sample ${\widetilde{C}}_{{ij}}$ and construct a kernel. Furthermore, the maximum charge capacity of CCDs are not well-defined in practice or in literature. It then becomes important to measure the signal range over which the Coulton et al. (2018) algorithm can reconstruct the impact of the BFE.

Correcting the BFE in flat-field images allows us to directly test the impact of higher-order effects at different signal levels. And correcting artificial star images allows us to test how well the Coulton et al. (2018) correction can reconstruct shapes at those signal levels. The underlying sub-pixel electric fields and higher-order physics are different in flat images and star images. Comparing the impact of the scalar correction on both flat-field images and star-field images is therefore a direct test of the validity of the underlying assumptions of the scalar correction.

We used the Bench for Optical Testing as described in Roodman et al. (2018). This test bench envelops the LSSTCam cryostat in a dark box that suppresses ambient light to a level <0.01 electrons per pixel per second and includes a rig to swap different optical projectors for illuminating the focal plane with various light patterns. This includes a spot grid projector, which projects a 49 × 49 grid of stars (approximately 3 cm × 3 cm on the focal plane and approximately 65 pixels of rectilinear spacing between spots) to mimic a star field of various brightnesses, shapes, sizes, and positions (see Figure 2). The entire projector can be moved using remotely controlled motorized stages in all three axes for dithering (XY) across the focal plane and focusing (Z). This projector consists of a Nikon 105 mm f/2.8 Al-s Micro-Nikkor lens, and an HTA Photomask photolithographic mask etched with pinholes. The spot grid pattern is illuminated by an integrating sphere with a 1'' opening, which is fed by a 3 mW fiber-coupled light emitting diode made by QPhotonics with a narrow band output peaked at 680 nm.

Zoom In Zoom Out Reset image size

For this study of the BFE, we selected one top-performing sensor of each type, one ITL (R03-S12) and one E2V (R24-S11), to independently test the BF correction. The R-number refers to the raft row and column location on the focal plane and the S-number refers to the sensor row and column location within the raft. We selected these sensors based on serial Charge Transfer Inefficiency (sCTI) criterion (Snyder & Roodman 2020). Both sensors have measured sCTI below 10⁻⁶ for all signal levels up to 10⁵ el, which is within the requirements for LSST set in O'Connor et al. (2016). We also selected two other sensors, R02-S00 (ITL) and R22-S02 (E2V), which were used only as a secondary checks to compare broader performances between sensor types which we will explain in a later section, and the results we will show and discuss in this paper come from analyses of the two "primary" sensors.

3.3.1. Calibrations

3.3.2. Artificial Stars

3.4.1. Instrument Signature Removal

3.4.2. Photometry

We first assess the dynamic range of the CCDs from our calibration data. We observe a correspondence between rapid changes of the covariance function with different metrics of pixel saturation or "turnoffs" defined in previous literature, which we further correlate with physical mechanisms of charge transport. Typically, the pixel saturation level is defined as the level above which the flat images dramatically lose variance (dC₀₀/d μ < 0) and no longer monotonically increases (the "PTC turnoff" as described by Janesick 2001). Another method defined by Snyder & Roodman (2020) assigns the limit as the flux above which charge can no longer be effectively transferred serially out of the sensor during readout. The serial charge transfer efficiency, sCTE (or rather sCTI in the case of charge transfer in-efficiency described in Rhodes et al. 2010) is calculated using the extended pixel edge response (EPER) method:

$$\begin{eqnarray}&&\mathrm{CTI}(\mu )=\displaystyle \frac{{F}_{\mathrm{overscan}}(\mu )}{{N}_{T}{F}_{\mathrm{lastpixel}}(\mu )}\end{eqnarray} \tag{ 6 }$$

where F_overscan/F_lastpixel is the ratio between the total charge left in the overscan region after an image is transferred out to the last image column and N_T is the number of pixel transfers during that readout, which is calculated for a specific flat image at a given flux level. The critical limit ("sCTI turnoff") is defined as the signal level at which the sCTI crosses above 10⁻⁵. Another method defines the maximum capacity more simply as the brightest recorded pixel ("maximum observed signal"). In our investigation, we discovered another relevant limit, which falls below all of these other limits. It can be identified as a turnoff in the parallel transfer component of CTI ("pCTI turnoff"), as shown in Figure 3, which is calculated analogously to the sCTI in Snyder & Roodman (2020). We define this point of pCTI turnoff by fitting the pCTI(μ) to a flat line above the level of noise (2.5 × 10⁴ el), and below any other features (5.0 × 10⁴ el), and determine where the pCTI deviates by more than 3σ given by the standard error of the fitted data points. We typically measure this limit to be between 0.8 and 1.1 × 10⁵ e⁻, which is below any of our other metrics for CCD full-well capacity.

Zoom In Zoom Out Reset image size

We will refer to the parallel EPER as parallel CTI since it is calculated analogously to serial CTI, however the origin of the turnoff in parallel CTI is still unknown. It is possible that this turnoff stems from the onset of a surface channel effect which persists from pixels in previous exposures that are above a critical signal level.

Each of these turnoff levels represents a different source of pixel correlations that begins to be relevant at different signal levels. And while we find the turnoff levels to vary from sensor to sensor, we find that their order is consistent in all the sensors we tested. The lowest is the pCTI turnoff (0.83–1.13 × 10⁵ el), followed by the PTC turnoff, the sCTI turnoff (1.40–1.75 ×10⁵ el), and the maximum observed signal (>1.50 × 10⁵ el). The different physical effects that determine these levels would bias the fit of our covariance model to the PTC, and our estimation of the strength of the BFE. We therefore limit the ranges of our fits to avoid including these other sources of pixel correlations that could be misattributed to the BFE.

We fit Equation (2) to our measured covariances from the zero signal limit to the pCTI turnoff for each channel in each of our sensors. We fit up to O( a ³) (the "full covariance model") with an 8 × 8 px covariance matrix (i, j = 0, 1, ⋯ 7) as a function of mean flat field signal. A summary of the fitted values for each sensor is shown in Table 1.

Table 1. Fitted Parameters of the Full Covariance Model (Equation (2)) to Calibrated Data (Section 3.3.1)

Cij (μ) Fit Parameters
Parameter	R03-S12 (ITL)	R24-S11 (E2V)	R02-S00 (ITL)	R21-S02 (E2V)
a00	−1.80 × 10−6	−3.11 × 10−6	−1.80 × 10−6	−3.10 × 10−6
a10	1.18 × 10−7	1.63 × 10−7	1.32 × 10−7	1.64 × 10−7
a01	2.54 × 10−7	3.83 × 10−7	2.62 × 10−7	3.78 × 10−7
σ(a00) / ∣a00∣	3.42%	3.44%	7.35%	3.94%
b00	−8.69 × 10−7	−2.95 × 10−7	−8.10 × 10−7	−5.38 × 10−7
b10	−5.36 × 10−6	−3.65 × 10−6	−5.97 × 10−6	−3.01 × 10−6
b01	4.71 × 10−8	2.04 × 10−6	5.42 × 10−7	1.06 × 10−6
σ(b00)/∣b00∣	50.8%	171%	163%	133%
$\bar{g}$	1.70	1.50	1.68	1.50

Note. We fit up to the pCTI turnoff, which was calibrated per amplifier, and then all parameters are averaged over all amplifiers for each sensor. Sigma parameters are given as the scatter over all amplifiers.

Download table as:  ASCIITypeset image

A good metric of the quality of the fit is the "sum rule" of a_ij , which should be close to zero if the total net pixel area displacement is zero and we are fully capturing the range of the BFE with the chosen size of the a_ij matrix. Figure 4 shows the zero sum rule of the fitted a parameter for different sizes of the matrix in sensor R03-S12 (ITL). We find that beyond 4 px the fractional pixel area change fluctuations are largely dominated by noise, and that modeling out to 8 pixels gives ∑a_ij < 1.04 × 10⁻⁷ el⁻¹, which is 5.8% of ∣a₀₀∣. Without including b in the fit, our sum rule is ∑a_ij < 5.91 × 10⁻⁹ el⁻¹, which is 0.345% of ∣a₀₀∣. The sum rule can deviate from zero as the result of poorly constraining the noise parameter in the fit to Equation (2). The excess sum then results from summing up a noise offset, which dominates the covariance at large lags away from the (0, 0) pixel. In the next section (Section 4.3) we describe how we overcome this issue and enforce the sum rule when ultimately deriving the correction kernel.

Zoom In Zoom Out Reset image size

In addition, the value of a₀₀ varies by 3.42% across amplifiers. Small variability in amplifier-dependent calibration levels (including nonlinearity, overscan correction, etc.) can result in small empirical variation in a₀₀, which should be amplifier-independent. Similar levels are reported for R24-S11 (E2V) in Table 1.

Figure 5 shows the observed a_ij matrix and profile of coefficients averaged over all channels. The a_ij matrix shows a notable anisotropy between the cardinal directions of the pixel coordinate system (a₀₁/a₁₀ ∼ 2.16), as shown in Figure 5, indicating the BFE is stronger in the parallel direction.

Zoom In Zoom Out Reset image size

The other parameter in these higher-order terms, the b_ij matrix, was previously measured to be sub-dominant by Astier et al. (2019), however we find that it is almost the same size as the a₀₀. Only the b_ij values for ${r}_{{ij}}=\sqrt{{i}^{2}+{j}^{2}}\lt 3$ are constrained by the fit to the data. Figure 6 shows the fitted value of the b matrix for R03-S12 (ITL). We measure b₀₀ = (−8.69 ± 1.11) × 10⁻⁷, b₁₀ = (−5.36 ± 0.224) × 10⁻⁶, and b₀₁ = (−4.71 ± 14.5) × 10⁻⁸, all in units of el⁻¹. The element b₀₀ is 48% of the magnitude of a₀₀. The b matrix has the effect of decreasing a_ij if negative and increasing a_ij if positive, corresponding to decreasing or increasing effective pixel areas, and it has the same units as the a_ij matrix. It is not immediately clear how to physically interpret the value of the central b_ij matrix elements or other non-zero elements of b_ij other than as modifiers to the mathematical effect of a_ij that enter in as higher-order terms of ${\widetilde{C}}_{{ij}}(\mu )$ at non-zero signal levels. That is to say b_ij represents nonlinear corrections to a_ij as charge accumulates.

Zoom In Zoom Out Reset image size

Figure 7 shows that ${\widetilde{C}}_{{ij}}(\mu )/{\mu }^{2}$ deviates from a_ij for electron signal levels greater than zero, which must be due to contributions from the nonlinear terms in Equation (2). In correspondence with Astier et al. (2019), we find that terms higher-order than μ¹ are a significant component of the covariance function, accounting for as much as 15%–30% of the loss in pixel variance at 10⁵ el, depending on the sensor tested, and are required for a good fit to the data. Excluding higher-order terms (μ⁴⁺) in the fit contributes a negligible error at least up to the pCTI turnoff level.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The fit to our model leaves an average χ²/N_dof of 13.7 in each amplifier, which is larger (χ²/N_dof = 14.9) when b = 0. This is large since the model fit leaves 0.5% fluctuation in the residual variance (C₀₀) above 5 × 10⁴ ADU (bottom plots of the 00 terms in Figure 7). The fluctuation is characteristic of a residual un-modeled nonlinearity in the Analog Signal Processing Integrated Circuit (described in Juramy et al. 2014; Astier et al. 2019). The nonlinearity takes on the form of a regular periodic function in signal level above 7.5 × 10⁴ el, and it might not have that much of an impact on our full covariance model if the surplus and deficit residuals balance the fit around the true PTC. However, by restricting the range of fluxes fit to those below the pCTI turnoff, we cut the fit range short enough that it might cause the nonlinearity to bias the fit and give us a poor result. For this reason, the quality of the model fit is very sensitive to the flux range. We found that adjusting the maximum signal by ±10³ el (around the pCTI turnoff), we could achieve arbitrarily better χ², but not necessarily a more accurate fit which models the BFE. On the other hand, the non-(0, 0) terms, which are also shown in Figure 7, exhibit good fits, with residuals smaller than the level of read noise and statistical fluctuations. In the next section, when we use the PTC model to construct a BF correction, we will continue to use these PTC fits, which at least generally follow the curves of the PTCs. Regardless, we can still conclude that the higher-order terms are non-negligible and result in nonlinear alterations of ${\widetilde{C}}_{{ij}}/{\mu }^{2}$ as charge distributions grow over time.

In this section, we attempt to linearize the PTCs of our sensors by correcting the BFE out of the flat-field images using a kernel K which includes higher-order BFE components.

To include these higher-order components into K, we construct the kernel (using Equation (3)) by sampling ${\widetilde{C}}_{{ij}}(\mu )$ from our full covariance model at a non-zero signal level, which would therefore include higher-order effects. We will do this through several steps:

• 1.  
  Evaluate ${\widetilde{C}}_{{ij}}(\mu )$ from the full covariance model at several signal levels.
• 2.  
  Use Equation (3), to calculate a kernel for each ${\widetilde{C}}_{{ij}}$.
• 3.  
  Correct the flat pairs in the PTC using each kernel.
• 4.  
  Calculate χ²/N_dof between each corrected PTC and the expected linear behavior, defined by the gain of the uncorrected PTC.
• 5.  
  Determine the kernel that results in the minimum reduced χ² statistic.

By finding the signal level at which the corresponding kernel minimizes the reduced χ² between the resulting corrected PTC and the Poisson expectation, we define a characteristic signal level that can best correct the PTC. We will refer to this signal level as the "characteristic" signal level (μ_*) and the corresponding kernel as the "characteristic" kernel.

We search for μ_* over a range of μ ∈ [10⁴ ADU, pCTI turnoff). We derive corrected PTCs from kernels at six signal levels in this range, in steps of 10⁴ ADU.

We also take a few practical steps to ensure each kernel is accurate. We normalized the kernel at a given signal level to enforce the zero sum rule (similar to charge conservation) by adjusting the value of the central pixel ${\widetilde{C}}_{00}$ so that ${\sum }_{{ij}}{\widetilde{C}}_{{ij}}/{\mu }^{2}=0$, for the measured correlations and the correlation model, integrating distances out to infinity. We risk over-fitting noise if we calculate the kernel boundaries beyond 4 pixels, so we model the correlations beyond 4 px away from the center of the kernel with an empirically measured power law, ${\widetilde{C}}_{{ij}}/{\mu }^{2}\propto {r}^{-\gamma }$ (typically γ ∼ 3.2). However, we observed that inclusion of the correlation model makes only a negligible difference to the kernel and the final correction. The central term of the kernel (K₀₀) is typically decreased by approximately 10% as a result of enforcing this zero-sum rule, and K₀₀ fluctuates across amplifiers by approximately 10% due to variations in calibrations between amplifiers. Enforcing this zero sum rule has a large impact on the correction since the (0, 0) pixel dominates the correction amplitude as we will show in the following sections. In order to avoid discontinuities at amplifier boundaries, we average the kernels per sensor and use the average in the final correction.

We measure the quality of the correction by comparing the corrected PTC to the expected Poisson form. We calculate the reduced χ² as:

$$\begin{eqnarray}&&{\chi }^{2}/{N}_{\mathrm{dof}}=\displaystyle \frac{1}{{N}_{\mathrm{dof}}}\displaystyle \sum _{\mu }\displaystyle \frac{{\left({C}_{00}^{\mathrm{Corrected}}-{C}_{00}^{\mathrm{Poisson}}\right)}^{2}}{{C}_{00}^{\mathrm{Poisson}}}.\end{eqnarray} \tag{ 7 }$$

The Poisson variance is determined using the gain derived from the covariance model fit of the uncorrected PTC. We also use N_dof = N(flat pairs) − 2 because we have N(flat pairs) observations and 2 parameters for gain and noise. We only calculate the residual χ²-statistic from data below 5.0 × 10⁴ el for E2V sensors and below 7.5 × 10⁴ el for ITL sensors to avoid residual NL at higher signals from interfering with our test.

Figure 8 shows the residual reduced χ² as a function of the corresponding signal level. We perform a simple least-squares test by fitting a three-parameter quadratic polynomial to the points across the flux space, and the location of the minimum of this quadratic polynomial determines the characteristic signal level μ_* that that minimizes the reduced χ² statistic for each sensor and best reconstructs the Poisson form of the PTC. The characteristic signal levels are not necessarily exactly half way between zero signal and pixel saturation at the PTC turnoff. The location of the characteristic signal level could be a balance between underestimating the role of higher-order BFEs at low signal and overestimating the role of higher-order BFEs at low signal. We also find that like-sensor types have similar characteristic signal levels within fitting errors between our 4 sensors, and while the reason is not entirely clear, it could be the result of sensor-specific operating conditions or design specifications (e.g., similar doping levels in the substrate material as considered in Rasmussen et al. 2016).

Zoom In Zoom Out Reset image size

Figure 9 shows the ultimate result of applying the empirically determined best BF correction to the flat field images in order to reconstruct the linear PTC. The correction breaks down above the pCTI turnoff for each sensor. Since the kernels for each amplifier are averaged together, the amplifier with the lowest pCTI turnoff determines the maximum signal level that can be corrected. Above 7.5 × 10⁴ el but below the pCTI turnoff, there is still some remaining uncorrected NL from the flat field projector, which prevents us from fully correcting the BFE above that signal level. This did not have an impact on our calculation of μ_* since we only used data points below this level in calculating our χ² statistic specifically to avoid it interfering with our analysis.

Zoom In Zoom Out Reset image size

We fit the full covariance model (up to the pCTI turnoff) to the BF-corrected PTC for both sensors. In the ITL sensor, the average a₀₀ decreased by 94.9% and the residual anisotropy was corrected closer to unity (a₀₁/a₁₀ = 0.735). For the E2V sensor, a₀₀ decreased by 97.1% and the residual anisotropy was also corrected closer to unity (a₀₁/a₁₀ = 0.942). This shows that the BFE correction captures most of the strength and anisotropy of the BFE in flat-fields below the pCTI turnoff.

We rely on the Coulton et al. (2018) algorithm that assumes the kernel is constant with respect to the flux. Our result finds this assumption is not true. Using a single kernel to correct multiple signal levels is intrinsically flawed because a single scalar kernel cannot capture both the BFE at low signal and the BFE at high signal. However, we cannot use multiple kernels to correct images of point-spread function (PSF) stars with more complicated light profiles. All previous literature that has investigated a correction of the BFE derived their correction from flat fields measurements. It is easily calibrated from flat fields and is far less computationally expensive than solving Poisson's equation for each image. If we use flat images to derive a correction for the BFE, it is important to find the kernel that can best reconstruct the Poisson noise of flat images and understand how well the BFE measured in flat fields can approximate the effect in star images.

In this section, we apply the characteristic correction derived from flat fields to correct images of artificial stars. The BFE in stars could be a stronger effect than in the earlier case of flat fields since the contrasts (pixel-to-pixel differences) are larger and therefore the gradients of the charge distribution and the divergence of the deflection fields are also larger.

To start, we found several problematic systematic effects in the images which needed to be addressed. First, there is non-uniform, disk-shaped background in the artificial star images that peaks around 0.5 el px⁻¹ s⁻¹, which is caused by illumination through the photolithographic mask and evident by eye in the top image of Figure 2. The background light level is much smaller than the Poisson noise of each star, and it is unlikely to significantly affect our photometry of bright stars. However the shape-fitting algorithm is known to be sensitive to background light levels, especially in the low signal-to-noise ratio (S/N) regime (Hirata & Seljak 2003; Refregier et al. 2012; Okura & Futamase 2018). We attempted to model and subtract the background in each exposure using a 65 px box size and 5 px median filter while masking out the footprints of the artificial stars. We attempted to remove the background in order to accurately measure the shape of the BFE at low signals, but there appeared to be a residual background that we failed to model by this procedure in the dimmest images, resulting in the increase of the second moment of stars in the low S/N region. We performed this background subtraction even though we ultimately ignored these low signal points in measuring the slope of the BFE in stars. The analyses with and without the background subtraction are not measurably different.

Second, we observe Airy diffraction rings on each of the stars due to the pinholes of the spot projector mask, which can be seen in the bottom image of Figure 2. However, as we will discuss further, we will study the differences between spot sizes at any signal level with respect to the same star's intrinsic shape at low signal and the impact of the detailed structure will be canceled out even if the diffraction rings of different spots have non-negligible overlap.

Third, the focal-plane-to-projector orientation and optical aberration from the lens caused minor distortions of the artificial stars, especially away from the center of the grid, and they cannot be perfectly parameterized by isotropic 2D Gaussians. We observed that the central and brightest stars of the grid hardly vary in intrinsic shape between exposures at the same or different projector positions. To avoid any other systematic effect related to the intrinsic alignment of the spots themselves, we therefore select only the top 5% brightest stars in total flux in each image. These stars happen to correspond to the same pinholes in the spot projector between in each image. In order to avoid including stars that land on the edges of the sensor, which skews their shapes, we also clipped stars beyond 3σ of the mean at each exposure time. This cut only removed O(10) data points in a few of the exposures. We ultimately selected approximately 80 spots in each exposure for our analysis. With 40 images at each exposure level, we use approximately 3200 spots/exposure level.

Figure 10 shows the approximate intrinsic shapes of the artificial stars used in this analysis, which we derive from their shapes at 20 s of exposure, which corresponds to a relatively low peak signal level of 2.5 × 10⁴ el. The distortions are correlated across the grid with the major axes consistently 7% larger than the minor axes, and the major axes of all the spots are aligned at a −24 deg angle relative to the serial readout direction (horizontal and to the left in Figure 10). We do not believe this affects our result, as we will discuss further in Section 5.1.2.

Zoom In Zoom Out Reset image size

In Figure 11, we show the growth of all three independent components of the second-moments matrix (I_xx , I_yy , I_xy ) separately for each sensor type. We chose the mean star shapes derived from the 20 s exposures as a fiducial point from which to show the growth of the size that result from the BFE since the lowest flux stars suffers from low S/N due to the extended residual background light. We therefore measure the size-growth of each star, corresponding to an individual hole in the projector mask, relative to its intrinsic shape. We also use the peak signal of each spot instead of the integrated signal as a measure of the signal level so that we can relate our spot photometry to different measurements of the pixel full-well capacity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To quantify the magnitude of the BFE, we determine the growth of the second moment with respect to this peak signal level. We fit the data points above the influence of our background light (${f}_{\max }\gt 2.5\times {10}^{4}$) and below the pCTI turnoff (which is our linear size-growth regime) with the function:

$$\begin{eqnarray}&&{I}_{\mu \nu }-{I}_{\mu \nu ,20\,{\rm{s}}}={\alpha }_{\mu \nu }{f}_{\max }+\,{c}_{\mu \nu },\end{eqnarray} \tag{ 8 }$$

where the indices run over the x and y , and α_{μ ν} determines the relative strength of the BFE, and c_{μ ν} just captures the difference between our fiducial point at I_{μ ν,20 s} and the true intrinsic shape of our spots in the limit of zero signal, which is important for fitting but unimportant for deriving the strength of the BFE.

We will refer to the size of stars using the canonical metric T = I_xx + I_yy . In addition, we will define the anisotropy of the BFE by the metric α_yy − α_xx . Note that this quantity measures the anisotropy of the BFE, not the anisotropy of the stars, which is typically defined using the pseudovector quantity of ellipticity (e₁, e₂), where e₁ = (I_xx − I_yy )/(I_xx + I_yy ) and e₂ = 2I_xy /(I_xx + I_yy ).

For the ITL sensor, the size of our artificial stars grow by 1.75% from their intrinsic size near the pCTI turnoff, and they grow by 3.25% for the E2V sensor. However, this size growth is asymmetric, just as we observed in the case of flat fields. For the sensors we measured, the BFE (α_{μ ν} ) is consistently 16% stronger in the parallel direction than in the serial direction for these stars.

Figure 11 also show the result of applying the correction derived using the characteristic kernel. We correct T = (I_xx + I_yy ) by 89.2% in ITL and 94.8% in E2V, which is worse than the correction on a₀₀ in flat fields (958%–97%), and we correct 31.48%–37.8% of the anisotropy in the BFE. We measure the I_xy component to be poorly corrected only because it is a small component of the star shape and measurement of the BF strength along this component is also small. Table 2 breaks down the detailed fit parameters between the two sensors. We also show the different turnoff points of pixel full-well discussed in Section 4.1 to show that the BFE has impacts over the full dynamic ranges of our sensors. Interestingly, we observe that the PTC turnoff has little effect on the shapes of the corrected spots in our study compared to the physical effects at the sCTI turnoff and the pCTI turnoff, which results in a small deviation from the linear relation as shape of the spots increases due to deferred charge above this level.

The charges are not lost when they are deflected, and a good physically motivated correction should conserve charge. While we implement the zero sum rule on a_ij and a zero-sum condition on the whole of K by adjusting the (0, 0) term (the way Coulton et al. 2018, defines the algorithm, and as it was previously implemented in the LSST science pipelines), the correction only conserves charge in the continuous limit ($\left\langle \delta F\right\rangle =0$). It does not conserve charge on small scales due to the application of the kernel in Equation (5). Figure 12 shows the difference in the measured charge of stars after applying the correction. We used 35 px pixel aperture for photometry, which is a large enough aperture compared to the PSF size (6 px FWHM) to capture the total footprint of each star without needing to assume a particular shape for the profile. The brightest unsaturated spots (below the pCTI turnoff) have a small surplus δ < 0.02% of total integrated charge compared to the same star in the uncorrected image. Charge non-conservation was more apparent in high-contrast images such as our artificial stars than in flat-fields, which would suggest that the error is the result of a local deviation of the correction from Gauss's Law as high-contrast images have larger deflection field divergences in the application of the correction. Improvements which enforce Gauss's law locally (on each pixel boundary) are are currently in development and undergoing a detailed study.

Zoom In Zoom Out Reset image size

In this section, we examine the assumptions that allow ${\widetilde{C}}_{{ij}}(\mu )$ to take the form μ²∇² K. The Coulton et al. (2018) approach assumes (1) that there are no higher-order terms that result in nonlinear BFE, and (2) the deflection field produced by accumulated charge in a pixel is conservative (zero curl).

5.1.1. Nonlinear BFE Components

Coulton et al. (2018) defines the kernel (Equation (3)) such that the gradient along one direction (∇_i K) represents the transverse displacement field along that direction. The impact of the BFE on flat-field covariances is the expected Poissonian behavior modified by the change in pixel area so that ${\widetilde{C}}_{{ij}}/{\mu }^{2}$ models an effective fractional area displacement matrix, and in the limit of zero signal ${\widetilde{C}}_{{ij}}/{\mu }^{2}\to {a}_{{ij}}$ (the full covariance model parameter of Astier et al. 2019).

The implication of a fixed kernel K in the Coulton et al. (2018) model is that ${\widetilde{C}}_{{ij}}/{\mu }^{2}$ does not change with signal level. Figure 7 shows the limited validity of the assumption that ${\widetilde{C}}_{{ij}}/{\mu }^{2}$ is independent of signal level. We show in the figure that ${\widetilde{C}}_{{ij}}/{\mu }^{2}$ agrees with the a_ij matrix by Astier et al. (2019) in Equation (2) only in the limit of zero signal, and that near the pCTI turnoff it deviates by 20%–30%, depending on the sensor.

The deviation in the PTC is well-modeled by the higher-order (nonlinear) terms of the full covariance model (Equation (2)). The main contributing higher-order term of the full covariance model has two components ( a ⨂ a + a b ) with two parameters ( a and b ). Coulton et al. (2018) do not distinguish between a and b or include their higher-order components.

The full covariance model does not a priori assume that these components have any specific physical mechanism. The a component is usually used as the principal BFE component and the b component strengthens (positive element) or weakens (negative element) the pixel boundary shifts nonlinearly with growing charge accumulation. The b parameter is mostly negligible at low flux as the fit quality our our calibration data to Equation (2) is slightly worse if we fix b to zero than if we leave it as a free parameter (per-amplifier average reduced weighted χ²/N_dof = 13.7 with b and χ²/N_dof = 14.9 without b ). While both fits had some residual NL, the primary difference between the fits occurs at high signal levels (near pCTI turnoff).

The non-(0, 0) components of b are also negative except for the nearest parallel component, These elements are also asymmetric, with the best-fit measurement of b₁₀/b₀₁ = 2.82, which could result from a growing asymmetry in pixel boundary shifts with charge accumulation. As charge accumulates, if the stored charge cloud distorts anisotropically, further charges would preferentially be deflected in a particular direction, producing a nonlinear charge feedback effect.

It has also been proposed in Magnier et al. (2018) and Holland et al. (2014) that different spatial distributions of space charges in the silicon could increase the effect of diffusion by altering the drift time of the incoming photocharges. Holland et al. (2014) believed this effect to be small for detectors operating with large applied voltages. In addition, other works have described the measured covariances well by only considering alterations in charge trajectories without considering any impacts on diffusion (Lage 2019; Astier & Regnault 2023). Typically, diffusion does not result in deviations from Poisson variance since diffusion does not change the fact that electrons collect independently into pixel bins. However, if some source charge lowers the drift field and produces and increases the diffusion scale, then the electrons no longer diffuse into neighboring pixels independently and could therefore produce a deviation from Poisson variances in the PTC. Extra diffusion due to collected source charge would therefore affect the PTC in a way that typical diffusion does not. If this effect was present in our sensors, this increase in diffusion would be indistinguishable from the BFE itself. It is likely that any BF-induced diffusion effects would get picked up and accounted for by the coefficients in our covariance model.

Figure 13 summarizes the correction of artificial stars in R03-S12 (ITL) for different kernels evaluated from different signal levels of the full covariance model, which encode the influence of varying amounts of the higher-order BFEs. The correction based on using ${\widetilde{C}}_{{ij}}(0)\sim {a}_{{ij}}$ over-corrects the artificial star images, and the correction based on ${\widetilde{C}}_{{ij}}({\mu }_{\mathrm{pCTI}})$ under-corrects the images. Some kind of mechanism to alter the BFE from the pixel boundary shift at non-zero signal level is needed. An adjustment of the strength is needed by constructing a kernel from the measured covariance model at some higher signal level. Scanning the whole covariance model to find the characteristic signal level that reconstructs the expected flat-field behavior (as we showed in Section 4.3) is the best method that we tried to balance the assumptions of the kernel approach and the higher-order BF contributions.

Zoom In Zoom Out Reset image size

The characteristic kernel however, still does not fully correct the BFE in flat fields or artificial stars (Figures 9 and 11). The correction also leaves differing residuals in different sensors (Table 2). It could be some error in enforcing the zero-sum rule, some other unmodeled source of measured pixel correlations (such as NL), or we could be missing some component in the physical model of the BFE in the correction. Astier et al. (2019) pointed out that all of the correction methods reported in literature only correct 90% of the initial effect in the images (Gruen et al. 2015; Guyonnet et al. 2015) and that all of the proposed correction techniques assume no contribution from higher-order terms. The flux-scanning might be applicable to these other correction techniques.

5.1.2. Non-zero Curl Components

The correction in Coulton et al. (2018) is entirely based on a scalar kernel. Therefore, the calculation of displaced charge in a given pixel, in application of Equation (5), only depends on the divergence of the kernel at discrete pixel centers. While the scalar kernel can capture anisotropy in boundary displacements along the two directions, it cannot fully capture the realistic complexity of distortions to all four boundaries on a pixel. Astier & Regnault (2023) have investigated modeling all four boundary shifts in a pixel from 3D electrostatic solutions in HSC sensors and measured a non-zero and non-negligible, discrete curl-component of the displaced boundaries sourced by the intrinsic anisotropy in the BFE. The curl component was also measured from simulations in Rasmussen et al. (2016). However, no vector-like curl-component of the deflected pixel boundaries field is included or correctable in the Coulton et al. (2018) model. In fact, they explicitly assume that this component is zero.

The validity of the zero-curl assumption in the LSSTCam can be partially constrained by comparing the correction of flat field covariances and artificial stars. The large pixel-to-pixel contrasts in star images requires accurate modeling of the small-scale physics to converge to a good correction. On the other hand, the small pixel-to-pixel contrasts in flat fields does not require the same accuracy at small scales to achieve a good correction. In the application of the kernel to artificial stars, the correction reaches the total convergence condition before it can properly redistribute charge in the most adjacent pixels, whereas in flat fields, with small pixel contrasts, the correction in the adjacent pixels makes up a larger total fraction of the charge redistribution at each iteration of the kernel. The correction in flat fields therefore redistributes more charge in nearby pixels, in the same number of iterations, before meeting the convergence condition than it does in artificial star images.

The correction might be correcting the anisotropy of the BFE better in flat fields than in stars because the kernel models the local physics better where unmodeled components are small. The BFE in the parallel direction is consistently 215%–234% larger than it is in the serial direction when measured from the a₀₁/a₁₀ in flat fields. And the parallel direction is consistently 16% stronger when measured by α_yy /α_xx in the size-growth of artificial stars. In Figure 13, only about 38% of the anisotropy in stars gets corrected by the algorithm. On the other hand, the anisotropy in flat field covariances is corrected by 60%–65%. The Coulton et al. (2018) algorithm could be correcting the anisotropy better in the case of flat fields in part because it better models the local physics in flat fields, where the unmodeled curl component would likely be smaller.

If we take an extreme case and purposefully apply the transposed kernel to the raw images—that means purposefully applying the correction with the anisotropic component of the kernel in the wrong direction (Figure 13)—the final anisotropy is roughly unchanged with the uncorrected case. This suggests that in artificial star images, most of the correction is dominated by the central pixels of the kernel, however we can see from charge conservation in the a matrix (Figure 4) that most of the charge redistribution that occurs is due to the pixels that are not immediately adjacent to the central pixel. This could be further evidence that the kernel is not accurately capturing the physics of the BFE at small distances.

We cannot conclude for certain the presence of a non-zero or nonnegligible curl component of the BFE distortions in LSSTCam from these data since there are not enough degrees of freedom to derive all four boundary shifts on a pixel from flat field statistics alone, even accounting for shared boundaries. It would require further modeling as is done in Astier & Regnault (2023), which includes modeling the 3D electrostatics in a pixel with each boundary shift as a free parameter, in order to fully constrain the curl component as a potential systematic.

We would like to note here that the measurement of anisotropy in the BFE in the sensors could be affected by the intrinsic shape of the stars and their alignment on the pixel grid if the BFE is sensitive to the gradient of light between pixels. The stars are intrinsically shaped and oriented on the pixel grid as shown in Figure 10 (based on their shapes at low signal level). Although there is little variation among the spots shapes in our sample, we attempt to measure how the BFE changes with respect to this intrinsic star shape in Figure 14. We measure the slope of the BFE individually for each star rather than the ensemble, as a function of intrinsic star shape, and we find that there is a small (but slightly beyond estimated error) decrease in the strength of the BFE in more extended sources. Since our stars are intrinsically larger in the y-direction than the x-direction, we would expect that α_xx in our sensors is slightly smaller than what we report and α_yy in our sensors is slightly larger than we report, therefore the overall anisotropy in our sensors is likely larger than what we report. However, given that there is only a small variation in intrinsic spot shapes, the error bars are on the same order of magnitude as our measurements in Figure 10, and it is impractical to extrapolate our measurements for perfectly round sources.

Zoom In Zoom Out Reset image size

The impact of BF residuals on future LSST science remains to be quantified, and steps can be taken to reduce the sensitivity of science to BFEs. BF can impact LSST weak lensing (WL) science by directly biasing survey targets, or indirectly by improper PSF estimation (Liaudat et al. 2023). BF distortion effects on the bright stars in our study are of O(1%), which is larger than the O(0.1%–1%) distribution of shape distortions caused by weak gravitational lensing (WL) in galaxy clusters, which is an observable sensitive to dark matter gravitational potentials along the line of sight and the dark energy equation of state (Huterer 2010; The LSST Dark Energy Science Collaboration et al. 2021).

The LSST Dark Energy Science Collaboration et al. (2021) science document requires us to reconstruct the true PSF size to a level of δ T_PSF/T_PSF < 10⁻³ in the co-added PSF in Y10, which in our lab-simulated stars is only satisfied below 2.5 × 10⁴ el in the ITL sensor and 2.3 × 10⁴ el in the E2V sensor, which covers only about 30% of the dynamic range of these sensors. The current LSST science requirements do not specify an error budget for the I_xy component.

The artificial stars in this study were nearly 2× larger than the expected LSST PSF size (07 or 3.5 px FWHM), and we would expect the BFE and the residual BFE to be worse for smaller PSFs due to the stronger charge gradients in the pixel plane. We also anticipate that bluer bands will experience an augmented BFE due to a shorter photon conversion depth within the sensor, thereby resulting in a strengthened BFE from the increased drift length (observed by Astier & Regnault 2023, in HSC and ongoing work to measure the effect in LSSTCam sensors). In addition, higher-order radial moment errors in the PSF are known to also impact WL observables, which is also a subject of future investigation (Zhang et al. 2021; Zhang et al. 2022).

Luckily, the BF correction benefits with the incorporation of higher-order BFEs into current correction schemes, and informed star selection from improved sensor characterization.

The correction of LSST PSF-like stars is good enough for LSST WL science requirements below a certain peak signal level, which will have to be chosen with regard to pre-defined accuracy requirements when selecting a PSF star sample. We will therefore need to apply the methods in this paper to all LSSTCam sensors to determine sensor-specific calibrated thresholds that will inform PSF star selection. Improvements to the correction will allow for higher thresholds and the inclusion of more stars with higher statistics and result more accurate and precise PSF estimation.

There are several steps that could be added and tested to the LSST Science Pipelines that could improve the correction:

• 1.  
  In the scalar model defined by Coulton et al. (2018), change the form of Equation (5) ($\hat{F}=F+\tfrac{1}{2}{\rm{\nabla }}\cdot V$, where V = F∇(K ⨂ F)) to enforce Gauss's law on small scales such that the charge lost over the area of one pixel is exactly matched by the sum of the charges lost over each boundary (X) of the pixel (such that ∫∇ · Vdx² = Σ_X V_X ). This will ensure charge conservation and improve the local modeling of charge transport due to the BFE (L. Miller 2024, Lance.Miller@physics.ox.ac.uk, in private communication).
• 2.  
  Include a curl component of the deflection field. This could be done by using a electrostatically derived vector-based model fit to the pixel–pixel correlations in flat fields, with each pixel boundary shift being a free parameter. Such a method is proposed in Astier & Regnault (2023) and tested for HSC, but it should be tested in LSSTCam.
• 3.  
  Improve the NL correction at higher signal levels in LSSTCam. The best-fit model parameters of the PTC are somewhat sensitive to the residual NL for these sensors since we fit up to the pCTI turnoff, which is above the signal level we start to observe significant signal NL. In addition, the PTC scanning method we describe is sensitive to the signal range on which we calculate the χ² for the same reason. While we mitigate this by limiting the signal range, the method we propose here does require the adequate removal of NL and other sources of pixel correlations, or it at least it requires plenty of flat-pairs at low enough signal levels that one can still accurately characterize the BFE without the impact of these other effects at high signal levels.