Mitigation of the Brighter-fatter Effect in the LSST Camera

Thick, fully depleted charge-coupled devices are known to exhibit nonlinear behavior at high signal levels due to the dynamic behavior of charges collecting in the potential wells of pixels, called the brighter-fatter effect (BFE). The effect results in distorted images of bright calibration stars, creating a flux-dependent point-spread function that if left unmitigated, could make up a large fraction of the error budget in Stage IV weak-lensing (WL) surveys such as the Legacy Survey of Space and Time (LSST). In this paper, we analyze image measurements of flat fields and artificial stars taken at different illumination levels with the LSST Camera (LSSTCam) at SLAC National Accelerator Laboratory in order to quantify this effect in the LSSTCam before and after a previously introduced correction technique. We observe that the BFE evolves anisotropically as a function of flux due to higher-order BFEs, which violates the fundamental assumption of this correction method. We then introduce a new method based on a physically motivated model to account for these higher-order terms in the correction, and then we test the modified correction on both data sets. We find that the new method corrects the effect in flat fields better than it corrects the effect in artificial stars, which we suggest is the result of sub-pixel physics not included in this correction model. We use these results to define a new metric for the full-well capacity of our sensors and advise image processing strategies to further limit the impact of the effect on LSST WL science pathways.


INTRODUCTION
The Legacy Survey of Space and Time (LSST) will be conducted with the Simonyi Survey Telescope at the Vera C. Rubin Observatory, which is under construction on the summit of Cerro Pachón in Chile.The survey plans to image 20,000 deg 2 of the sky with O(100), 15 second visits over 10 years across six filters (ugrizy) in the optical and near-infrared (NIR) parts of the spectrum.It will map galaxies and optical transients to understand the natures of dark energy and dark matter and their impacts on the formation of structure in the universe (The LSST Dark Energy Science Collaboration et al. 2021).
The instrument for this is a 3.2 gigapixel camera, which contains 201 individual charge-coupled devices (CCDs), with 189 designated specifically for science imaging.The sensors are fully depleted high-resistivity bulk silicon CCDs developed by two separate vendors.One type is made by Imaging Technology Laboratories (ITL), and the other is made by Teledyne e2v (E2V) to similar general architectural specifications.These sensors are arranged by type into 3 × 3 groups called raft-tower modules (RTMs) or rafts, which can each operate as an independent camera.Each CCD sensor is 4 cm × 4 cm made up of sixteen 1 megapixel channels each read out by its own amplifier and readout electronics.Each pixel on the LSST Camera focal plane is 10 µm × 10 µm and has a depth of 100 µm (Lopez et al. 2018).

arXiv:2312.03115v1 [astro-ph.IM] 5 Dec 2023
Much of the design and current development of instrumentation for LSST focuses on reducing the impact of systematic sensor artifacts in order to produce subpercent level precision measurements of cosmological parameters and test currently prevailing thermodynamic models of the universe and theories of dark matter (Albrecht et al. 2006).
To meet these requirements, we need to understand the systematic effects in our sensors at the far limits of the capabilities of our instrumentation (Ivezić et al. 2019).Many of the problematic sensor behaviors are spatially static and can be simply calibrated or modeled as intrinsic constants at a particular location on a sensor and easily removed from raw images.However, this method would not work with locally variable or other signal-dependent effects, and such effects are known to exist (Stubbs 2013;Antilogus et al. 2014;Astier et al. 2019).
Correlations between neighboring pixels have been shown to arise at high signal levels, at which point captured photocharges produce significant transverse electric fields on incoming photocharges (Downing et al. 2006;Holland et al. 2014;Lage et al. 2017).During the integration of an exposure, photo-electrons deflect into neighboring pixels in reaction to quasistatic changes in effective pixel area from the accumulated charges in the potential wells of the pixels, causing the measured light profile of a bright source to differ from that source's intrinsic surface brightness profile.This effect, dubbed the brighter-fatter effect (BFE), broadens intrinsic surface brightness profiles.The magnitude of the BFE depends on the surface brightness profile of the source itself and thus cannot be modeled solely by its location on a sensor.The consequence is that the BFE breaks the critical assumption of experimental imaging analysis that the pixels are independent light collectors that perfectly obey Poisson statistics.
The BFE has been observed in LSSTCam sensors by Antilogus et al. (2014) and Lage et al. (2017), and in detectors used in other astronomical cameras such as Hyper Suprime-Cam (HSC) by Coulton et al. (2018), the Dark Energy Camera used by the Dark Energy Survey (DES) by Gruen et al. (2015), the Wide Field Camera 3 H1RG detector of the Hubble Space Telescope by Plazas et al. (2017), MegaCam by Guyonnet et al. (2015), the Mid-Infrared Instrument (MIRI) on board the James Webb Space Telescope (JWST) by Argyriou et al. (2023), and in the near-infrared (NIR) detectors of the Wide Field Imager of NASA's Nancy Grace Roman Space Telescope (Plazas et al. 2018;Hirata & Choi 2020;Choi & Hirata 2020;Freudenburg et al. 2020;Plazas Malagón et al. 2023).These studies measured a devia-tion from the Poissonian behavior of pixels in flat field images at bright illuminations, which they attributed to the BFE.Several corrections have been proposed by Antilogus et al. (2014), Gruen et al. (2015), and Coulton et al. (2018), the last being a scalar algorithm which is currently used by in the HSC data reduction pipeline and is currently planned to be used in the LSST science pipelines (Bosch et al. 2018;Jurić et al. 2017;Bosch et al. 2019).
In this paper, we will directly measure the BFE in the LSST Camera and our ability to correct it.In §2, we will introduce the theory behind the BFE and its correction by Coulton et al. (2018) and weigh it in light of more recent findings by Astier et al. (2019).In §3, we will describe our laboratory measurements and image processing.In §4.1-4.2, we will explore the dynamic response of our sensors to charges in these data.In §4.3-4.4,we will construct an improved application of the correction by Coulton et al. (2018) and apply it to flat fields and artificial stars.In §5, we will discuss the relative results between these two cases and directly test the underlying assumptions of the correction made by Coulton et al. (2018).Finally in §5.2, we discuss how residual BFEs could influence several of LSST's dark energy science goals and motivate future work on analysis structures to mitigate its impact.

BRIGHTER-FATTER CORRECTION
Several studies have shown that the pixel covariances in flat field images can be used to measure the BFE in the sensors of other optical and near-IR survey cameras (Antilogus et al. 2014;Gruen et al. 2015;Guyonnet et al. 2015;Coulton et al. 2018;Astier, Pierre & Regnault, Nicolas 2023;Hirata & Choi 2020;Freudenburg et al. 2020;Plazas et al. 2017;Plazas et al. 2018).The covariance can be numerically formalized between an arbitrary central pixel, x = (0, 0), assumed to be far from the edge of the detector, and a neighboring pixel x ′ = (i, j) in the difference of two flat field exposures, F 1 and F 2 , with the same nominal signal level and µ = avg(F 1 , F 2 ): The directionality on each sensor is defined in regard to the CCDs' front-side architecture, where the vertical, û = (0, 1), direction is along the "parallel" readout direction, and the horizontal, û = (1, 0), direction is along the "serial" readout direction.Figure 1 shows an example of a densely sampled measurement of the variance function C 00 (µ) from a single sensor (ITL) on the LSST focal plane.If all the pixels are independent from each other, one can expect a simple linear relationship from the Poisson statistics as where δ ij is Kronecker's delta, g (el/ADU) is the gain conversion between physical electrons and recorded counts (analog-digital units or ADU) by the analog-to-digital signal converter (ADC), and n ij [e 2 ] is the noise at a pixel (i, j).The measurement of C 00 departs from the linear relationship at high signal levels as excesses due to Poisson noise are suppressed due to the BFE since charges would be deflected out of pixels with excess accumulated charge and into neighboring pixels with less accumulated charge, resulting in a correlation between pixels at higher flat field signal levels.We measure a 30% loss in the expected variance function (photon transfer curve or PTC) near pixel saturation (Figure 1), which is consistent with the size of the effect observed by Astier et al. (2019), Coulton et al. (2018), Antilogus et al. (2014), Gruen et al. (2015), Guyonnet et al. (2015), Coulton et al. (2018), andAstier, Pierre &Regnault, Nicolas (2023).
If charge is conserved, the sum of all the covariances away from this central pixel recovers all the pixel-level Poissonian sensor behavior lost due to the BFE (first proposed by Downing et al. 2006).The central pixel and the immediate neighbors should contribute the most to recovering the noise budget and integrating further away from the central pixel should add vanishingly smaller contributions.Figure 1 also shows that accounting for the covariances 8px away from (0, 0) in both directions on the focal plane allows us to reconstruct the Poissonian behavior in flat field exposures within statistical fluctuations up to fluxes above which the variance begins to drop off due to pixel saturation (around 1.35 × 10 5 el), for the two sensors that we tested.

Modeling Pixel-Area Changes from Flat Field Statistics
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: where a ij [1/e − ] describes the fractional change in pixel area due to changes in the pixel boundaries from accumulated source charges, b ij [1/e − ] describes smaller time-dependent processes that weaken or strengthen the We assume that the covariance matrix has parity symmetry about the central pixel, and we sum the covariance matrix from |i, j| < 2 and |i, j| < 8, and we show that summing out to 8 pixels fully reconstructs the Poissonian behavior.We used 342 total exposure pairs with ×1.025 log-spacing in signal space, which were corrected with the basic instrument signature removal defined in Appendix A. We group the data into 100 bins from 0 to 1.75 × 10 5 el.This was calculated for one sensor made by ITL on the LSST focal plane, which had an average gain of 1.702 el / ADU across all amplifiers.
BFE as we will explore in a later section, µ is a given mean signal over the image region, g is the sensor gain [e − /ADU], n ij is a constant term that represents the mean square noise (n 00 ∼ 5.7e 2 , where the non-(0, 0) terms of the noise matrix are contributed by the electronics and overscan), ab is a direct matrix multiplication, and ⊗ refers to a discrete convolution operation.This expansion allows for higher-order terms which are non-linear BFEs.Astier et al. (2019) found this model to fit their data well, and the contribution from the non-linear terms was not negligible (with the µ 2 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.

Deriving a Scalar Correction with Higher-Order Effects
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): where The true image can thus be recovered from the measured image using equation 22 of Coulton et al. (2018): where 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(µ 4 ) fits the observed PTC well over the flux range they considered, citing χ 2 /N dof = 1.2 for the fit to the variance element over all sensor channels in LSSTCam sensors, where they use N dof = N(flat field pairs)−3.This suggests that the PTC shape is contributed by other non-linear 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 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 and linearizing the PTC allows us to directly test the impact of higherorder, flux-dependent effects.And correcting artificial star images allows us to test how well the Coulton et al. (2018) correction can reconstruct shapes at all signal levels.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.

DATA ACQUISITION AND IMAGE PROCESSING
Our measurements of the BFE are derived from datasets taken during the fifth electro-optical testing period for the LSST Camera (informally referred to as Run 5), which took place on the integrated camera focal plane at the SLAC National Accelerator Laboratory in December 2021.In this section we describe the laboratory setup, data acquisition, and analysis methodology.This includes the configurations of the sensor testbed, acquisitions of calibration datasets and artificial star data, image processing pipeline, and photometry.

Laboratory Setup
We used the Bench for Optical Testing (BOT) as described in Newbry et al. (2018).This test bench envelops the LSST Camera 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 3cm × 3cm 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.

Sensors Tested
For this study of the BFE, we selected one topperforming 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 for all channels up to 10 5 e − , well below 10 −6 , which is within the needed 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.

Calibrations
We remove the static sensor effects from our lab exposures by subtracting calibrated images, which include 20 bias images to remove fixed-pattern noise, 20 dark images to eliminate thermal currents, and 342 densely sampled pairs of flat images from low flux (50 e − /px) to high flux (1.75 × 10 5 e − /px) with a log-scale increment of ×1.025 in the SDSS i -band regime (700-800 nm) to measure the covariance function, deferred charge, and linearity.For the flat pair acquisitions, we also have National Institute of Standards and Technology (NIST)calibrated photo-diode measurements to provide an independent accurate measurement of the throughput of the light in the flat-field acquisition.
We note that the SDSS i -band (700 nm-800 nm), in which we measure our covariance function, covers a region of the electromagnetic spectrum outside the wavelength of the artificial star data (680 nm).Chromatic dependence has not yet been rigorously studied in LSST sensors, however, based on preliminary studies by Astier, Pierre & Regnault, Nicolas (2023) in HSC, we believe any potential chromatic dependence of the BFE is negligible at these wavelengths that are comparatively close together.Further study of the BFE in the LSST Camera will investigate these hypothesized effects.

Artificial Stars
We projected artificial stars onto four sensors using the spot grid projector.An example exposure is shown in Figure 2. We project these spots onto the sensors during an exposure and fix the integrated light level by adjusting the integration time.We took a sequence of images of artificial stars at increasing intensities by varying the integration time in 15 steps from 5 s to 200 s to cover the full range in flux of our flat field calibrations (the full dynamic range of our CCDs), with 40 images at each exposure level.We split these image acquisitions in four quadrants around each sensor, changing the projector position 4 times across the whole imaging sequence.This allows us to measure many spots at a variety of positions and signal levels.We also ignore the first image after the projector changes position to avoid images with distorted stars that result from shaking the projector.

Instrument Signature Removal (ISR)
We processed each raw image using a complete sequence of the standard LSST Science Pipelines (w.2023.29 releases) and our calibrated Run 5 testing data products (Jurić et al. 2017;Bosch et al. 2019).The standard Instrument Signature Removal (ISR) removes spatially characteristic systematic effects in the sensors, which includes bias subtraction with medianper-row overscan subtraction, dark current subtraction, non-linearity correction, and bad-pixel masking, and serial deferred charge.A detailed list of configurations used in the software are given in Appendix A.
We correct analog-to-digital converter (ADC) nonlinearity (Downing et al. 2006) using a calibration correction informed by the densely sampled flat pairs compared with an independent measure of flux by the photodiode measurements.Additionally, we corrected the charge transfer inefficiency in the horizontal (serial) direction across the image using the same method as described in Snyder & Roodman (2019) of Astier et al. (2019).We also applied bad-pixel masks and interpolated over the locations of bright and dark pixels and any cosmic-ray artifacts, though these pixels were excluded from calculations of the flat-field covariances.We intentionally turned off the saturated pixel repair to avoid any non-physical distortions introduced by the repair.It is at this stage that we have the option to correct the BFE as well.The detailed set of options we chose for ISR is provided in Appendix A.
We then apply the BFE correction itself.As described and recommended in the implementation of the correction (Coulton et al. 2018), we recursively applied the kernel to each image until the total added charge in each step falls below a threshold of 10 electrons (the "convergence condition"), which typically takes no more than 2 or 3 iterations.

Photometry
We measured the light profile and shape statistics of each artificial star from each exposure, deriving the source centroid, flux, and shape using the Galsim (Rowe et al. 2014) re-implementations of the Gaussianweighted adaptive moments algorithm (HSM) defined by Hirata & Seljak (2003) and tested independently by Mandelbaum et al. (2005).The HSM photometrics were derived for a square stamp around each star's centroid, with a side length equal to the spot grid spacing, which should entirely capture the star and any diffraction rings without presuming anything about the stars' shapes (see the bottom of Figure 2).The metrics include the diagonal elements of the quadrupole image moment tensor (here we denote them I xx , I yy , and I xy ) where x is the serial direction along the orientation of the gates and y is the parallel direction along the orientation of the channel stops.These statistics are useful because they are directly related to metrics of cosmic shear distortions.

Improved Benchmark Measurements of Full-Well Capacity
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 00 /dµ < 0) and no longer monotonically increases (the "PTC turnoff" as described by Janesick 2001).Another method defined by Snyder & Roodman (2019) 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 inefficiency described in Rhodes et al. 2010) is calculated using the extended pixel edge response (EPER) method: 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 −5 .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 (2019).We define this point of pCTI turnoff by fitting the pCTI(µ) to a flat line above the level of noise (2.5×10 4 e − ), and below any other features (5.0 × 10 4 e − ), 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 − 1.1 × 10 5 e − , which is below any of our other metrics for CCD full-well capacity.
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 5 e − ), followed by the PTC turnoff, the sCTI turnoff (1.40 − 1.75 × 10 5 e − ), and the maximum observed signal (> 1.50 × 10 5 e − ).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.

Fitting the Covariance Function
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 3 ) (the "full covariance model") with an 8 × 8px covariance matrix (i, j = 0, 1, • • • 7) as a function of mean flat field signal.
Figure 4 shows the zero-sum rule of a 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, where a ij < 1.04 × 10 −7 el −1 , which is 5.8% of |a 00 |.Without including b in the fit, our sum rule is a ij < 5.91 × 10 −9 el −1 , which is only 0.345% of |a 00 | and fully captures the zero-sum rule.The value of a 00 varies by 3.42% across amplifiers simply due to small differences in amplifier operation.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., and the resulting fit parameters for both sensors are shown in Table 1.The a ij matrix shows a notable anisotropy between the cardinal directions of the pixel coordinate system (a 01 /a 10 ∼ 2.16), as shown in Figure 5, indicating the BFE is stronger in the parallel direction.This anisotropy changes as a function of signal also due to the higher-order terms in equation 2.
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 roughly the same order of magnitude as the a ij matrix, and is a sensitive variable in the full model fit.Figure 6 shows Figure 7 shows that for some i, j, C ij (µ)/µ 2 deviates from a ij for electron signal levels greater than zero, which must be due to contributions from the non-linear terms in equation 2. In correspondence with Astier et al. (2019), we find that terms higher order than µ 1 are a significant component of the covariance function, accounting for as much as 15-30% of the loss in pixel variance at 10 5 el, depending on the sensor tested, and are required for a good fit to the data.Excluding higher-order terms (µ 4+ ) in the fit contributes a negligible error at least up to the pCTI turnoff level.
The fit to our model leaves an average χ 2 /N dof of 13.7 in each amplifier, which is larger (χ 2 /N dof = 14.9) when b = 0, which is unusually large as the model fit leaves 0.5% fluctuation in the residual variance (C 00 ) above 5 × 10 4 ADU (bottom plots of Figure 7).The fluctuation has the same shape as a residual un-modeled nonlinearity in the Analogue Signal Processing Integrated  Circuit (ASPIC) (described in Astier et al. 2019;Juramy et al. 2014).The non-linearity takes on the form of a regular periodic function in signal level above 7.5 × 10 4 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 non-linearity 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 3 el (around the pCTI turnoff), we could achieve arbitrarily better χ 2 , 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.By visual inspection of the PTC fit, it appears that the model accurately follows the PTC curve.In the next section, when we use the PTC model, we will discuss how we mitigate the impact of the poor NL correction.
Regardless, we find that the higher-order terms are non-negligible and result in non-linear alterations of C ij /µ 2 as charge distributions grow over time.

Correction of the BFE in Flat Fields
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 C ij (µ) from our full covariance model at a non-zero signal level, which would therefore include higher-order effects.We will do this by: 1. Calculating C ij (µ) from the full covariance model at several signal levels.
2. Use equation 3, to calculate a kernel for each C ij .
3. Correct the PTC using each kernel.
4. Calculate the χ 2 /N dof between the corrected PTCs and the expected linear behavior defined by the gain of the uncorrected PTC.
5. Determine the kernel that best corrects the PTC.
This defines a characteristic signal level that can best correct the PTC (minimize the χ 2 ).We will refer to this signal level as the "ideal" signal level (µ * ) and the corresponding kernel as the "ideal" kernel.Residual χ 2 values from a scan over the scaleconversion factor parameter space for all four sensors.We calculate the χ 2 for each amplifier from the residuals below the pCTI turnoff between the corrected PTC and the expected linear behavior derived from the uncorrected PTC gain for each amplifier.We then sum all the χ 2 values for each amplifier.A simple 3-parameter quadratic has been fit to the points across the parameter space to determine the factor that minimizes the χ 2 statistic.This shows that we have a characteristic signal level that reconstructs the PTC, and this level is similar in sensors of the same type.
We search a range of µ * ∈ [ 10 4 ADU, pCTI turnoff ) in steps of 10 4 ADU.This is the range within which we expect the χ 2 curve as a function of signal level to follow a simple quadratic as it will not be significantly impacted by any systematic other than the BFE.We also only calculate the residual χ 2 -statistic from data below 5.0 × 10 4 el for E2V sensors and below 7.5 × 10 4 el for ITL sensors to avoid residual NL at higher signals from interfering with our test.We calculate the χ 2statistic from the variance component, exactly as Astier et al. (2019), which we described in §2.2.
We also take a few practical steps to ensure an accurate kernel.We normalized the kernel to enforce the zero sum rule (similar to charge conservation) by adjusting the central value of the central point (0, 0) of the kernel so that ij a ij = 0, for the measured correlations and the correlation model integrated out to infinity.We risk over-fitting noise if we calculate the kernel boundaries beyond 4 pixels, so we model the correlations beyond 4px away from the center of the kernel with an empirically measured power law, C ij /µ 2 ∝ r −γ (typically γ ∼ 3.2).However, we've observed that inclusion of the correlation model makes only a negligible difference to the kernel and the final correction.This central term (K 00 ) is typically decreased by approximately 10% as a result of this zero-sum rule, and K 00 fluctuates across channels/amplifiers by approximately 10% just due to differences in performance 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.
Figure 8 shows the result of a simple least-squares test which determines the best flux level to reconstruct the PTC.Most of the ideal signal levels are around half the range between zero signal and the pCTI turnoff, which could be a balance between underestimating the role of higher-order BFEs at low signal and potential deviations in the underlying assumptions of the Coulton et al. ( 2018) correction scheme at higher signal levels.We also find that like-sensor types have similar ideal signal levels within fitting errors between our 4 sensors, and while the reason is not entirely clear, it could be the result of a physical parameter such as similar operating conditions or an artifact of the manufacturing (e.g.similar doping levels in the substrate material as considered in Rasmussen et al. 2016).
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 4 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 χ 2 statistic specifically to avoid it interfering with our analysis.
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 00 decreased by 94.9% and the residual anisotropy was corrected closer to unity (a 01 /a 10 = 0.735).For the E2V sensor, a 00 decreased by 97.1% and the residual anisotropy was also corrected closer to unity (a 01 /a 10 = 0.942).This shows that the BFE correction is properly modeling the strength and anisotropy of the BFE in flat-fields below the pCTI turnoff.

Correction of the BFE in Artificial Stars
We also measure the BFE on artificial stars and apply the same ideal correction derived from flat fields.The BFE in stars could be a stronger effect than in the ear- lier 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.Firstly, there is a wide illumination background in the artificial star images that peaks around 0.5 e − /px/s, which is evident by eye in the top image of Figure 2.This is likely caused by excess transmission of light through the photo-lithographic mask.The background light level is much smaller than the Poisson noise of each star and is unlikely to significantly affect our photometry.However the shape-fitting algorithm is known to be sensitive to background light levels (Hirata & Seljak 2003;Refregier et al. 2012;Okura & Futamase 2018), especially in the low signal-to-noise (SNR) regime, which will matter for our stars with peak signals below 2.5 × 10 4 e − , so we model and subtract the background in each exposure using a 65px box size and 5px median filter, masking out the positions of the artificial stars.
Secondly, 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 cancelled out even if the diffraction rings of different spots have non-negligible overlap.
Thirdly, 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 15s of exposure, which corresponds to a relatively low peak signal level of 2.5 × 10 4 e − .The distortions are correlated across the grid with the major axes are 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 §5.1.2.
In Figure 4.4, 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 15-second 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 signal-to-noise ratio due to the extended background.We therefore measure the sizegrowth 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.
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 > 2.5 × 10 4 ) and below the pCTI turnoff (which is our linear size-growth regime) with the function: 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 µν,15s 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.
For the ITL sensor, the size of our artificial stars grow by 1.75% from their intrinsic size (T P SF = I xx + I yy ) 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-19% stronger in the parallel direction than in the serial direction for these stars.
Figure 4.4 also show the result of applying the correction derived using the ideal kernel.We correct T = (I xx + I yy ) by 89.9% in ITL and 94.1% in E2V, which is worse than the correction on a 00 in flat fields (95-97%), and the correction preserves 67.2-67.8% of the anisotropy in the BFE.The fact that the correction of I xy is worse than the other components is likely an artifact of it being a small in magnitude.Table 2 breaks down the detailed fit parameters between the two sensors we studied.
We also show the different turnoff points of pixel fullwell discussed in §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.

Charge Conservation
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 (⟨δF ⟩ = 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 35px pixel aperture for photometry, which is a large enough aperture compared to the PSF size (6px 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.

Testing the Assumptions of Scalar BFE Theory
In this section, we examine the assumptions that allow C ij (µ) to take the form µ 2 ∇ 2 K.The Coulton et al. (2018) approach assumes (1) that there are no higherorder terms that result in non-linear BFE, and (2) the deflection field produced by accumulated charge in a pixel is conservative (zero curl).5.1.1.Non-Linear BFE Components Coulton et al. (2018) defines the kernel (equation 3) such that the gradient along one direction (∇ i K) rep- resents 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 C ij / µ 2 models an effective fractional area displacement matrix, and in the limit of zero signal C ij / µ 2 → 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 C ij / µ 2 does not change with signal level.Figure 7 shows the limited validity of the assumption that C ij / µ 2 is independent of signal level.We show in the figure that C ij / µ 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 (non-linear) terms of equation 2. The main contributing higher-order term of the full covariance model has two components (a⊗a+ab) with two parameters (a and b).The a component is the principal BFE component and the b component strengthens (positive element) or weakens (negative element) the pixel boundary shifts 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 χ 2 /N dof = 13.7 with b and χ 2 /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).However, Coulton et al. (2018) do not distinguish between a and b or include their higher order components.
The physical meaning of b is not as clear as a.The negative b 00 element could have the same impact as a "space-charge effect" as the drift field decreases due to the accumulated charge and therefore increases the drift time of the electron and the time that the BFE has to act on those incoming charges as they drift into the pixel.The effective area of the central pixel effectively decreases due to the space-charge effect and would therefore produce negative best-fit values of b 00 .
The non-(0,0) components of b are also negative except for the parallel terms, which could also be the result of the space-charge effect interfered with by residual pCTI.These elements are also asymmetric, with the best-fit measurement of b 10 /b 01 = 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 as a type of charge feedback effect.
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 higherorder BFEs.The correction based on using C ij (0) ∼ a ij over-corrects the artificial star images, and the correction based on Cij (µ 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 ideal signal level that reconstructs the expected flatfield behavior (as we showed in §4.3) is the best method to balance the assumptions of the kernel approach and the higher-order BF contributions.
The ideal 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 still leave on the order of 10% of the initial effect in the images (Guyonnet et al. 2015;Gruen et al. 2015) and that all of the proposed correction techniques assume no contribution from higher-order terms.

Non-Zero Curl Components
The correction algorithm in Coulton et al. (2018) is entirely based on scalar products and therefore no vectorlike curl-component of the deflection field is included or correctable in this model.The correction in Equation 5only contains a divergence component, however on small scales there could be a non-negligible curl-component.A non-zero curl component could be generated from the intrinsic anisotropy that exists in the sensors.
The validity of the zero-curl assumption is directly tested by comparing the correction of anisotropy in flat field covariances and artificial stars.Flat field pixel correlations have a strong dependence on the divergence component of the BFE and a comparatively small dependence on the curl component if any exists due to the spatial uniformity of the charge distribution.Pixel correlations in artificial star images by contrast could depend strongly on a curl component of the BFE.The BFE in the parallel direction is consistently 215-234% larger than it is in the serial direction when measured from the a 01 /a 10 in flat fields.And the parallel direction is consistently 16-19% stronger when measured by the parallel second moment of the artificial stars.
In (Figure 13), only about 32% 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 corrects anisotropy better in the case of flat fields simply because it better models the impact of local physics on flat field statistics, where the curl component is smaller.
This could also impact the total correction level between flat fields and stars.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 that capture the anisotropy, 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.The correction in flat fields therefore redistributes more charge in the same number of iterations before meeting the convergence condition than it does in artificial star images.If we take an extreme case and purposefully apply the transposed kernel to the raw images-that means 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 is further evidence that the curl component is not negligible.
An improvement to the algorithm might then consider a vector approach, which deals with each pixel boundary.However, 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 (a N i,j = a S i,j+1 and a E i,j = a W i+1,j ), so this method would require further modeling with each boundary left as a free parameter.Work is ongoing in the development and testing of other approaches like this (such as the one proposed by Astier, Pierre & Regnault, Nicolas 2023).
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.

Future Work
The impact of BF residuals on future LSST science remains to be quantified, and steps can be taken to reduce The artificial stars in this study were nearly 2× larger than the expected LSST PSF size (0.7" or 3.5px 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, Pierre & Regnault, Nicolas 2023, in HSC and ongoiong 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.

Proposed Improvements
The correction of LSST PSF-like stars is good enough for LSST weak lensing science requirements below a certain signal level, and this helps to inform the selection of PSF stars to estimate the PSF.We will therefore need to apply the methods in this paper to all LSSTCam sensors since different sensor types and sensors can have different characteristics overall.
In addition, 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 ( F = F + 1 2 ∇•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 ∇ • V dx 2 = Σ X V X ).This will ensure charge conservation and improve the local modeling of charge transport due to the BFE (Lance Miller, Lance.Miller@physics.ox.ac.uk, in private communication).If the lack of charge conservation originates from the small scale deviations from physics, then this implementation should fully conserve charge.
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, Pierre & Regnault, Nicolas (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 χ 2 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.

CONCLUSIONS
In this paper we measured the BFE in two types of sensors in the LSSTCam focal plane from flat-field and artificial stars in the lab.We then quantified and qualified the scalar BF correction proposed by Coulton et al. (2018) with improvements to account for sources of higher-order pixel correlations.
First, we measured dense flat pair data to characterize the PTC, and we found that the BFE gets stronger with signal due to higher-order, non-linear BFE components in a model by Astier et al. (2019), which includes flux-dependent drift fields and feedback mechanisms on accumulating charges.The BF correction proposed in Coulton et al. (2018), which uses a scalar kernel derived from the pixel covariance matrix, relies on the assumption of linearity.The merit of the correction therefore depends on the signal level from which we calculate the covariance matrix and derive the kernel, which is a priori ambiguous.
Secondly, we introduced a new procedure to derive the kernel in a way that includes the contribution of nonlinear BFE terms in the correction.We determined the kernel that best reconstructs the linear PTC.We then used this kernel to correct ithe BFE n both sensor types by 95-97%, but was only able to correct up to a flux we identify as the parallel deferred charge turnoff at 10 5 el, where we can no longer efficiently transfer charge along an image column during readout.The flux range of pixels that can be used for weak lensing science is severely limited by this cutoff.
Thirdly, we corrected the BFE in artificially produced stars by 96 − 99% in the serial direction and by 85 − 90% in the parallel direction.The BFE was observed to be anisotropic (16%) relative to the CCD coordinate system, and the correction preserved 67% of this anisotropy.
We eliminated other systematics that could have interfered with these results, including non-circular artificial stars and charge-conservation, to identify the source of these effects as intrinsic BFEs in the sensor itself and limitations of the assumptions of the Coulton et al. (2018) correction.Our data suggests that the higherorder BFEs are not negligible and that the zero-curl assumption of the Coulton et al. (2018) algorithm is not fully valid.
Our findings also motivate a detailed study on more realistic PSF stars and how measurement errors from BF could ultimately impact cosmology and other science goals.Ultimately, it is important to characterize the sensitivity of cosmological parameters to observables biased with BFEs.Even with state-of-the-art correction techniques, the residual effects could represent a significant component of the systematics error budget for cosmological analyses of LSST observations.APPENDIX

A. INSTRUMENT SIGNATURE REMOVAL CONFIGURATIONS
In Table 3 we present the configurations we used with the LSST Science Pipelines for ISR of the spot images (Jurić et al. 2017;Bosch et al. 2019).For spot identification and characterization, we use an LSST Science Pipeline wrapper called the Mixed Calibration Optics Analysis Test Library (MixCOATL).Our data processing also utilized code (with default configurations) for matching and labeling spots to specific holes in the lithographic mask of the optical projector (GridFitTask).This code is described in Esteves et al. (2023) and integrated into the LSST Science Pipelines.For each task, including the calibration and ISR tasks, we used the default configuration parameters.
After this sequence of ISR, we modeled and subtracted the background light produced by the spot grid projector.Given that the background light produces photocharges that get collected in the potential wells of our sensor, it contributes to the BFE, and should be taken into account when applying the scalar kernel and correcting an image.However, this component should be removed before shape-fitting.We performed this background subtraction separately for both the corrected and uncorrected images before calculating and cataloging the final shape statistics.We performed this using a 65 × 65px median filter to match the artificial star spacing, masking over the fitted footprints of the artificial stars.

Figure 1 .
Figure1.The variance and integrated covariance matrix vs. flux for a series of flat field images on an LSST camera sensor, averaging the covariance matrix for all amplifiers.We assume that the covariance matrix has parity symmetry about the central pixel, and we sum the covariance matrix from |i, j| < 2 and |i, j| < 8, and we show that summing out to 8 pixels fully reconstructs the Poissonian behavior.We used 342 total exposure pairs with ×1.025 log-spacing in signal space, which were corrected with the basic instrument signature removal defined in Appendix A. We group the data into 100 bins from 0 to 1.75 × 10 5 el.This was calculated for one sensor made by ITL on the LSST focal plane, which had an average gain of 1.702 el / ADU across all amplifiers.

Figure 2 .
Figure 2. Top: an example of a spot image (20 s exposure with 680 nm LED) that has been processed by standard instrument signature removal discussed in §3.4.1.It shows the uneven (circularly symmetric) illumination pattern from the projector on a single sensor (ITL R03-S12).Bottom: an image subsection showing the Airy function systematic due to the long wavelength regime (solid red), and the corresponding 65 px × 65 px stamp (dashed-red) which sets the bounds used to derive calculate each star's shape.

Figure 3 .
Figure 3. Measured parallel transfer CTI as functions of flux for the two sensors, with the range of turnoffs among the 16 channels highlighted.The anomalous amplifier in sensor R03-S12 (AMP06/C15) was ignored for all of our analyses due to abnormal pCTI behavior.

Figure 4 .
Figure4.The charge conservation condition of the aij model for the cases where bij is allowed to vary or held fixed at 0. We calculate the sum over the a matrix.This was calculated for the ITL sensor R03-S12.

Figure 5 .
Figure 5.The pixel fractional area change as a matrix derived from the full covariance model (equation 2).The top panel shows the asymmetry in direction (a01/a10 ∼ 2.16), and the bottom panel shows the decay of the matrix into noise past 3 − 4 pixels.We also show the serial and parallel pixels in red and blue to show the difference in the model in both directions on the sensor.The central term is a00 = −1.80 × 10 −6 .Errorbars are errors on the mean across all 16 amplifiers.This was calculated for an ITL sensor on the LSST focal plane (R03-S12).

Figure 6 .
Figure 6.The measured bij charge-displacement component derived from the full covariance model (equation 2).We also show the serial and parallel pixels in red and blue to show the difference in the model in both directions on the sensor.It shows the relative increase (positive) or decrease (negative) shift in the accumulated charges' centroids.This was calculated for an ITL sensor on the LSST focal plane (R03-S12).

Figure 7 .
Figure 7.The higher-order term of the covariance for the ITL sensor (R03-S12).Top: the analytically calculated first-order component of the Astier et al. (2019) covariance model, showing the increase in the fractional pixel area change with signal level and the corresponding fitted and scaled parameters for both sensor types.Bottom: the fractional residuals between the measured covariances and the corresponding covariance model.The vertical lines mark the various measures of full-well conditions.

Figure 8 .
Figure8.Residual χ 2 values from a scan over the scaleconversion factor parameter space for all four sensors.We calculate the χ 2 for each amplifier from the residuals below the pCTI turnoff between the corrected PTC and the expected linear behavior derived from the uncorrected PTC gain for each amplifier.We then sum all the χ 2 values for each amplifier.A simple 3-parameter quadratic has been fit to the points across the parameter space to determine the factor that minimizes the χ 2 statistic.This shows that we have a characteristic signal level that reconstructs the PTC, and this level is similar in sensors of the same type.

Figure 9 .
Figure9.Photon transfer curves reconstructed from corrected flat field images using the corresponding ideal kernels for the ITL sensor (R03-S12) on the left and the E2V sensor (R24-S11) on the right.On the bottom, we include the residuals from the expected linear behavior.Although we use a detector averaged kernel, we compute this for each of the 16 amplifiers.The average gain, included as a parameter in the post-correction full covariance model fit, increased by 0.35% in the ITL sensor and 16% in the E2V sensor.

Figure 10 .
Figure 10.Distributions of the intrinsic shapes of the artificial stars used in these analyses (after selection cuts).Each distribution is centered on its mean and plotted with the same bin size.We derived these values from a collection of low-brightness images at the same exposure time (15 s), and these distributions contain about 3200 spots.The distributions indicate a small average ellipticity (a/b < 0.07) and correlated orientation angle ( θ = −24.1 deg) due to the relative alignment of the spot projector with the mask and focal plane, indicating that the major axes are almost parallel with the serial gate structure.Angles are given by tan(2θ) = 2Ixy/(Ixx−Iyy) and measured relative to the serial transfer direction.

Figure 11 .
Figure 11.The second moments of the artificial stars taken from images vs. their peak pixel value for an ITL sensor (R03-S12).Left column is the uncorrected data and the right column is corrected data.The rows correspond to Ixx, Iyy, and Ixy from top to bottom.Each blue point is a single spot from a single exposure.Error bars are standard errors of the distribution about the mean at each exposure time in moment and flux, and a line has been fitted to the 15-50 s points, which corresponds to 2.5 − 9.0 × 10 4 el (just below the pCTI turnoff).The vertical shaded regions represent the range of full-well capacities of the 16 amplifiers on the sensor using the four different methods.The blue shaded region is the boundary to model the shape to within 1 part-per-thousand.

Figure 12 .
Figure 12.The quantity f = r<35px F / r<35px F is the ratio of integrated flux within 35px of the centroid of each spot, after/before the BFE correction is applied on each of the spots in the ITL sensor (R03-S12) data.We only plot spots with peak signals below the pCTI turnoff.Error bars are standard errors of the distribution at each exposure time.

Figure 13 .
Figure13.The strength of the BFE measured from artificial stars (fit to equation 7) before/after applying scalar corrections derived from different configurations of Cij for sensor R03-S12 (ITL).We show the fitted slopes for the uncorrected artificial stars as well as for the artificial stars corrected using the covariance model sampled at zero signal, sampled at the pCTI turnoff (minimum value of all amplifiers), the derived ideal signal level from our full PTC scan, and for the last case, we also show the correction with the kernel transposed (purposefully applied in the wrong orientation), which shows how sensitive the correction is to the fixed-plane BF anisotropy.The shaded area represents the total allowed residual BFE (αxx + αyy) to meet the LSST science requirements (The LSST Dark Energy Science Collaboration et al. 2021) on TP SF = Ixx + Iyy to the pCTI saturation limit (10 5 el) for stars of size ⟨TP SF ⟩ = 20 px 2 , under the assumption that our artificially produced stars are real, measured PSF stars.

Figure 14 .
Figure14.The decrease in the influence of the BFE in more extended sources, which would have smaller flux gradients in any direction of the pixel plane.We show the slope of the BFE by component, derived from the 81 individual spots used in our analysis, where each point corresponds to a single hole in the spot mask, relative to their intrinsic shapes.We derived these slopes from the same calculation in equation 7. The errors shown are the fitting errors of the slope.These were calculated for the ITL sensor (R03-S12).
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 P SF /T P SF < 10 −3 in the co-added PSF in Y10, which in our lab-simulated stars is only satisfied below 2.5 × 10 4 el in the ITL sensor and 2.3 × 10 4 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.

Table 1 .
Fitted parameters of the full covariance model (equation 2) to calibrated data ( §3.3.1).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.