Issue 
A&A
Volume 670, February 2023



Article Number  A118  
Number of page(s)  17  
Section  Astronomical instrumentation  
DOI  https://doi.org/10.1051/00046361/202245407  
Published online  17 February 2023 
Correction of the brighterfatter effect on the CCDs of Hyper SuprimeCam
LPNHE, CNRS/IN2P3, Sorbonne Université, Université ParisCité, Laboratoire de Physique Nucléaire et de Hautes Énergies,
75005
Paris, France
email: pierre.astier@in2p3.fr
Received:
8
November
2022
Accepted:
6
December
2022
The brighterfatter effect affects all CCD sensors to various degrees. Deepdepleted thick sensors are seriously affected and the measurement of galaxy shapes for cosmic shear measurements requires an accurate correction of the effect in science images. We describe the whole correction chain we have implemented for the CCDs of the Hyper SuprimeCam imager on the Subaru Telescope. We derive nonlinearity corrections from a new sequence of flat field images, and measure their statistics, namely their twopixel function. We constrain an electrostatic model from flat field statistics that we use to correct science images. We find evidence that some fraction of the observed variance and some covariances is not due to the combination of Poisson statistics and electrostatics – and the cause remains elusive. We then have to ignore some measurements when deriving the electrostatic model. Over a wide range of image qualities and in the 5 bands of the imager, stars in corrected science images exhibit size variations with flux small enough to predict the point spread function for faint objects to an accuracy better than 10^{−3} for the trace of second moments – and even better for the ellipticity and the fourth radial moment. This performance is sufficient for upcoming largescale cosmic shear surveys such as Rubin/LSST.
Key words: instrumentation: detectors / techniques: image processing
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
The brighterfatter (BF) effect refers to a dynamical image distortion that affects CCD sensors. The most spectacular manifestation of the effect is that bright stars appear slightly bigger in size than faint ones, a manifestation that is reflected the very name of the effect. All studies of the effect have attributed it to distortions of the drift electric field sourced by the charges stored in the pixel potential wells during image integration. This modifies the apparent shape of bright objects and the twopoint statistics of uniform exposures, and thick CCDs are more vulnerable to the effect than thinner sensors. Evidence of the effect and the physical explanation can be found in Guyonnet et al. (2015, and references therein, G15 hereafter), together with the relation of the effect with nontrivial flat field statistics. Electrostatic calculations are shown to reproduce the data in Rasmussen et al. (2016) and in Lage et al. (2021) for a specific sensor.
On deepdepleted thick CCDs, bright stars thus generally appear bigger by a few percent than faint stars, compromising the modeling of the image point spread function (PSF) at a level that is not tolerable for largescale cosmic shear measurements (see e.g., Mandelbaum et al. 2018). The correction method proposed in G15 relies on flat field statistics to constrain the correction applied to science images, which mostly consists of correcting the recorded image for dynamically displaced pixel boundaries. The method has been implemented for DECam in Gruen et al. (2015), and for Hyper SuprimeCam (HSC) on the Subaru Telescope in Coulton et al. (2018, C18 hereafter) with some minor differences with respect to G15.
The method proposed in G15 relies on firstorder perturbations both in the modeling of flat field correlations and when correcting science images. In Astier et al. (2019, A19 hereafter), the relation between pixel area alterations and flat field statistics is extended to higher orders, which removes significant biases from the analysis. In the same paper, correcting for nonlinearity of the video chain is shown to play a potentially important role when constraining the BF effect from flat fields. One other evolution since the G15 proposal is that detailed electrostatic calculations have been shown to reproduce the measured flatfield statistics, when the mandatory manufacturing data is available (see, in particular, Lage et al. 2021). However, the CCD vendors do not necessarily release this data, or they do not even have it available to the required level of accuracy. As we show in this paper, there is also some detectable demographic variability among the CCDs of the HSC camera, which are all of a unique type from a single vendor. Gruen et al. (2015) also detected some variability among DECam CCDs. Thus, constraining the image corrections from measurements of the actual sensors is still in order.
In the present paper, we revisit the BF correction for the HSC camera described in C18 and applied to science images in Mandelbaum et al. (2018). We take advantage of a new flatfield sequence that allows us to redetermine both the nonlinearity correction and the twopoint correlation function of flat fields. We apply to these images some potentially important corrections (described in A19) and we fit a variance and covariance model that has been improved since C18. We also propose a different approach for transforming the information extracted from flat fields into the correction of science images. Finally, we test the correction of the images separately over a broad range of image qualities and in the five bands of the camera.
The flow of the paper is as follows. We first detail in Sect. 2 why it is necessary to correct nonlinearities prior to measuring flat field statistics and the nonlinearity measurement itself. In Sect. 3, we describe the measurements of flat field statistics and the fit of the measurements, as well as the variability observed among the sensors. The information extracted from flat field statistics is fundamentally insufficient to correct the science images, thus, we describe in Sect. 4 the electrostatic model we use to derive the correction from the flatfield results and the outcome of different fits we perform. Once we obtain models that allow us to correct science images, we apply those to real data, as described in Sect. 5, and we compare the various outcomes. We face the evidence that the BF correction is inadequate for the y band and we compute a physically motivated reduced correction for this band, which we eventually apply to the science data. In Sect. 6, we compute some PSF modeling quality indicators commonly used in the context of shear estimation and conclude that the quality of our correction exceeds the requirements for a largescale cosmic shear survey such as Rubin/LSST. It also fulfills the less stringent requirements of the PSF fidelity implied by photometry accuracy of highredshift faint supernovae in order to estimate their distances, which is our initial motivation for this work.
2 Nonlinearity correction
2.1 Importance of the nonlinearity correction
Following A19, we assume that the effective area A of a pixel is linearly altered by the charge content of the sensor: (1)
where a_{ĳ} is a characteristic of the sensor and Q_{ĳ} denotes the charge content of the image, which evolves as light integration goes on. The indices i and j refer to distances in pixel units along the serial and parallel directions, respectively. In the same coordinates, δA applies to the pixel located at (i = 0, j = 0). Conventionally, a_{ĳ} is expressed in el^{−1}, Q_{ĳ} in ADU and g is the gain in el ADU^{−1}. We note that, equivalently, we may set g = 1 and express Q in electrons. From parity symmetry, we assume that: (2)
so that we measure the a_{ĳ} on the i, j ⩾ 0 quadrant. If we consider a uniform image, with all Q_{ij} identical, then δA has to be zero, because of translation symmetry. This imposes the following “sum rule”: (3)
where the sum extends to the four quadrants. This sum rule can also be regarded as a consequence of area conservation. Since a_{ij} is the fractional pixel area change for a unit source charge, we refer to these quantities as “area coefficients”. Since samesign charges repel each other, a pixel shrinks as it fills up, and the selfinteraction coefficient a_{00} is negative. For regular CCD operating conditions, all the other coefficients turn out to be positive. The sum rule hence indicates that is much larger (in absolute value) than any other coefficient.
Following Eq. (1), the expression of pixel covariances in uniform exposures reads (see the derivation in A19): (4)
where µ is the average (in ADU) of the uniform exposure. n_{ĳ} refers to noises (expressed in el^{2}), and n_{00} refers to the usual read noise. From this expression, the relation between the variance of uniform exposures and their average (usually called the photon transfer curve or PTC) can be expressed as: (5)
where the first term is the read noise, the second is the Poisson term, and the quadratic term is the first and largest contribution of the BF effect: the slope of the PTC decays with the signal level (because a_{00} < 0). The fact that the PTC of CCDs grows less rapidly than the signal level was initially noted in Downing et al. (2006), who also noted that thanks to positive covariances, grouping data into bigger pixels tends to restore the Poisson behavior. These nontrivial flatfield statistics and the brighterfatter effect per se were unified under the same physics in Antilogus et al. (2014).
The analog readout chain can be affected by nonlinearity and a common distortion is a quadratic contribution: (6)
where µ_{m} and µ_{t} refer to measured and true values, and k is a small quantity describing the distortion. The last term can equivalently refer to either value of µ. For the variances, we have: (7)
Thus, comparing Eqs. (7) and (5), we find that a quadratic nonlinearity, if it is not corrected for, biases the measurement of a_{00}, which is the largest coefficient describing the BF effect. The same calculation applied to covariances shows that there is no effect at the µ^{2} level, hence the measurement of the other a_{ij} coefficients is much less sensitive to a quadratic nonlinearity than a_{00}.
2.2 Measurement of response nonlinearity
The Hyper SuprimeCam instrument is installed at the prime focus of the 8.2m Subaru telescope on the MaunaKea summit in Hawaii. A detailed description of the instrument can be found in Miyazaki et al. (2018) and here we provide only the salient features for our work. The HSC camera contains 104 science sensors, which are deepdepleted Hamamatsu 2k×4k CCDs with four readout channels on each sensor. The sensors are 200 µm thick, the pixel size is 15 µm, which corresponds to 0.185′′, on average, over the focal plane. The saturation of sensors occurs around 120 000 el and the video channels have a typical gain of 3.2 el ADU^{−1}. All channels may a priori exhibit different nonlinearities.
In order to measure the nonlinearity of the camera video chains, we rely on a new sequence acquired on 2019/10/30, of uniform exposures known as “dome flats” obtained at night by illuminating a screen attached to the dome, toward which the telescope is pointed. The signal level is controlled by varying the exposure time. The illumination system is not equipped with a photodiode or any similar device and we use the exposure time as a proxy for the amount of light received by the CCD. Obviously, this is an acceptable proxy only if the illumination system is stable over the whole sequence or if its variations can be modeled. In order to disentangle a drift of the lamp intensity from a genuine nonlinearity, the various exposure times are not ordered timewise in a monotonie way, as shown in Fig. 1.
For a given video channel, we collect from every exposure (labeled i an average signal, µ_{i} (in ADUs), an exposure time, T_{i}, and the (Julian) date at which it was acquired, d_{i}. We minimize the quantity: (8)
where S represents a sevenknot spline meant to describe the lamp variation with time, O is a potential systematic shutter offset, k represents the quadratic nonlinearity, and T_{i} is the required exposure time. The shutter offset is meant to account for a (small) difference between the required and realized exposure times. We do not find any compelling evidence for a global shutter offset, nor a spatial trend in the focal plane, that could arise from such a large blade shutter. Hence, we set O = 0. Since the illumination may evolve along the sequence in slightly different ways across the focal plane, we do not enforce the spline S to be the same for all channels. We eventually compare the outcome of different fits (on the same sensor) in Fig. 2. The fact that these corrections are almost indistinguishable indicates that the nonlinearity and the lamp variations can be robustly separated from this data set.
Figure 3 displays the fit residuals (from Eq. (8))as a function of signal level or of date. Those are globally small. A more flexible model as a function of date would not reduce much these residuals, highly correlated over the channels at the same date. We can attribute this scatter to random fluctuations of the lamp or of the shutter. The fit residuals as a function of signal level do not strongly call for a more flexible or more general nonlinearity model than a quadratic correction: the largest systematic residual is found around 10 000 ADUs and amounts to a few ADUs. However, for some yet unknown reason, we had to discard exposures before d_{i} – 58786 < 0.255 because they exhibit large residuals (tens of ADUs) consistently over channels. One might attribute this to some erratic behavior of the lamp that vanishes after about 0.5 h. While these first exposures are rejected from the nonlinearity analysis, they are used to measure the flat field statistics in what follows.
In Fig. 4, we display the values of quadratic nonlinearity in the HSC focal plane in order to visualize the absence of geometrical trend that could occur from a spatially variable evolution of the illuminating system. The distribution peaks around k = 10^{−7} ADU^{−1}, which corresponds to a correction of ~0.4% at the maximum CCD signal. We go on to show soon that a_{00} averages to −1.3 × 10^{−6} and the channel gains are around 3; if one ignores quadratic nonlinearity (Eq. (7)), the bias affecting a_{00} amounts to ~8%. a_{00} is the best measured area coefficient and essentially drives the scale of the BF correction. We might note that such a bias unavoidably causes a violation of the sum rule and, hence, translating such a set of area coefficients into a correction is somehow arbitrary because sensible image corrections have to derive from areas changes that sum up to zero. We also note that eliminating this trivial cause of sum rule violation allows us to detect more subtle problems in Sect. 4.3.
Fig. 1 CCD signal values as a function of time for a single channel during a sequence of uniform exposures illuminated by a lamp in the dome. The intensities are deliberately not ordered in a monotonie way, so that the drift of the illumination system can be separated from the nonlinearity. The sequence was designed by N. Yasuda (IPMU). 
Fig. 2 Spline correction from the minimization of Eq. (8) for four different channels, normalized to its average. Since this is presumably due to lamp variations, it should be similar for all channels. The maximum difference is in the 10^{−4} range. 
Fig. 3 Fit residuals as a function of date (top) and signal level (bottorn) for all channels of all sensors. There are no obvious needs for a more sophisticated model as a function of date nor signal level. The rms residual is about 3 × 10^{−4}µ. 
Fig. 4 Values of measured quadratic non linearity k (in Eq. (8)) in the focal plane (top) and their distribution (bottom). Most of the channels cluster around 10^{−7}, which corresponds to a correction of ~0.4% around maximum CCD signal. 
Fig. 5 Nearest serial neighbor deferred signals for the four channels of sensor 42. The values are small and vary almost linearly with the signal. The model continuous lines are the cubic polynomial model that we use for correction. 
3 Measurement and fit of flatfield statistics
3.1 Measurements
Measuring the variances and covariances of flat fields allows one to measure the dynamical alteration of pixel areas due to electrostatic distortions; namely, we measure variances and covariances of flat field exposures (which we may call the twopoint function of the exposure) at various signal levels, in order to measure the a_{ĳ} coefficients (Eq. (1)) by fitting Eq. (4) to the measurements. We closely follow the procedure described in A19: we first correct for nonlinearity (see the previous section) and correct for deferred signals.
A deferred signal is by definition collected after (in the serial time chain) the pixel it belongs to. In order to measure those, we closely follow the procedure described in A19, and find that there are no significant deferred signals beyond the nearest serial neighbor. We display a few examples of the nearest neighbor deferred signals in Fig. 5: those are small and almost linear, indicating that they mostly originate from electronics. We tabulate one correction curve per channel, and for correcting the effect, we typically “place back” less than 0.1% of a pixel into the preceding one. Failing to carry out this correction results in the covariance, C_{10}, exhibiting a linear contribution (∝µ), which cannot be absorbed by the model of Eq. (4), hence biasing a_{10}. For our CCDs, at the maximum signal level, the contribution of uncorrected deferred signals to C_{10} would represent about 25% of the electrostatic value.
When computing image statistics, we mask any outlier pixels detected in the images and use the pixel mask applied to science images at the epoch of the sequence. These pixel masks flag pixels that repeatedly depart from their neighbors in uniform exposures and pixels on the sides of the focal plane which are heavily vignetted (or fully blind). We account for the masked pixels in covariance computations in a standard way, described in Appendix A of A19. The bias of variance and covariances due to masking outlier pixels from images themselves is compensated for as described in Sect. 5.2 of A19.
We measure each video channel independently and produce variance and covariance curves using the 66 flat pairs of the sequence. We measure the covariances on the difference of flat pairs acquired at the same signal level to eliminate contributions from illumination nonuniformity or sensor sensitivity variations. These 66 flat pairs were acquired in g band. Contrarily to C18, we carry out the measurements on all HSC science sensors.
Fig. 6 Measurement and fits of the variance (C_{00}) and the two nearest neighbor covariances (C_{10} and C_{01}) for the four channels of sensor 42. The channels have been offset for visual convenience. The decay of C_{00}/µ with µ, originally noted in Downing et al. (2006), is the most obvious effect of the BF effect. The spread of measurements is compatible with shot noise. 
3.2 Fits of the covariance curves
We then fit the measurements using Eq. (4), with weights derived from shot noise. This equation contains terms up to µ^{4}, and including one extra order would change the predictions by a few 10^{−5}, which is about 100 times less that the smallest shot noise. Figure 6 displays the measurements and fits for a few (i, j) pairs and for the four channels of a sensor. One can readily note that the slopes for C_{10} (serial neighbors) and C_{01} (parallel neighbors) differ by a factor of more than 4, which is due to the very different mechanisms that confine the charges into the pixel wells for the serial and the parallel directions: while the boundaries between columns are static and implanted, the boundaries between rows result from the potentials of the parallel clock stripes that are swung during image readout in order to move charges towards the serial register and the output amplifier. The parallel pixel boundaries are then somehow weaker than the serial ones and this is reflected on the pixel covariances of uniform exposures. We perform the measurements in the range i, j < 10 and fit Eq. (4) to all the measurements of a given channel at once.
In order to question whether the fitted model (Eq. (4)) properly describes the data, we plot in Fig. 7 the fit residuals of the variance and the nearest covariances. Since we are studying effects related to sensors, we turned the averages and covariances into gainfree quantities by converting them into electrons, using the gains obtained from the fits in the previous paragraph. An examination of Eq. (4) shows that the coefficients of the µ^{3} terms are entirely determined by combinations of the coefficients of the µ^{2} terms. The residuals in Fig. 7 exhibit some curvature, which is compatible with a mismatch between the coefficients of µ^{2} and µ^{3}. For the residuals of C_{00} (top), the curvature of the residuals represents about 15% of the cubic term of Eq. (4) (thus illustrating the importance of these higher order terms). For C_{01}, the curvature represents about 30% of the cubic term, while they have similar values for the residuals to C_{10}. We will propose later an explanation for the curvature of residuals of C_{00} and C_{01}, but this explanation does not apply to C_{10}.
Fig. 7 Residuals of fits of variance and covariance measurements, all expressed in electrons using the gains from the fits, so that the data from different sensors can be averaged. The red signs are binned averages of the residuals from individual fits. 
Fig. 8 Distribution of area coefficients over the ~400 hundred channels of HSC (blue) and the averages over each sensor (orange). All a_{ij} coefficients are expressed in el^{−1}. We may note that the trends for channels and sensors are similar, which indicates some homogeneity within sensors (except for Σ a, where the scatter is dominated by shot noise). Σ a is the sum for −10 < i, j < 10. 
3.3 Demography of HSC sensors
We may now investigate the demography of channels and sensors regarding the BF effect. In Fig. 8, we display scatter plots of the wellmeasured a coefficients. We can see the correlations, which are similar at the channel and sensor level, indicating that the variations are weaker within sensors. We can note that a_{00} does not correlate to the other coefficients, and that the apparent correlations are in a direction that does not preserve the “sum rule.” This may indicate that they are not due to genuine variations of the BF effect among the sensors. We may also observe that on average the sum of a is positive: (9)
Since all a_{ij} at larger distances are positive, the sum to infinity is even farther from 0. Since < a_{00} >≃ −1.3 × 10^{−6}, the violation of the sum rule is hence at least 5% of a_{00} and we eventually see that it reaches ~10% once we have a model to evaluate the largedistance contributions. This means that the sum of all covariances rises with signal level faster than Poisson, at variance with expectations from statistics (see Eq. (8) in A19). Our measurements are thus affected by some noise, which increases with signal level and contributes as µ^{2}.
Since the variations of wellmeasured a_{ij} is at most 10% across channels, we decide to average the measurements in order to average the shot noise at large distance. Once equipped with this average, we fit an electrostatic model to it that, by construction, satisfies the sum rule. In the next section, we show that the mismatch between the model and the data is indeed localized.
4 Electrostatic fit
4.1 From the area coefficients to science image corrections
We study the BF effect in order to suppress its impacts from the science images. The scheme proposed in G15 consists of evaluating from the science image itself the motions of pixel boundaries with respect to a perfect grid, as well as evaluating via interpolation how much charge was flowing over these pixel boundaries during integration, and then placing this amount of charge back where it belongs.
We can readily note that the measurements we have performed so far constrain the change in the area of a pixel, but do not tell how its shape is altered. If we describe the shape change by different motions of the serial and parallel sides of a pixel, this simplistic shape description already requires two quantities per pixel, while we only have one. Defining a shape change by distinct serial and parallel boundary shifts means in turn that we have to “split” the a_{ij} coefficients (into , , ·, ) and rely on areaconserving symmetries such as . We can also regard a_{ij} as the discrete divergence of the pixel boundary displacement field and we are faced with the illposed problem of determining a vector field from its divergence, in a 2D discrete space.
Various implementations of the correction method proposed in G15 differ in the way they perform this promotion of pixel area change into two directional components. In G15 and in Gruen et al. (2015), some ratios are imposed, which were estimated from the area changes themselves. C18 argue that the motions of pixel boundaries is a 2D vector field that results from the (perturbation) electric field sourced by the charges that represent the image, and, hence, the 2d displacement field is (as is the 3d electric field) curlfree. The “scalar” a_{ĳ} field can then be transformed into a 2D vector curlfree field. While this curlfree assumption seems appealing, we show in Sect. 4.4 that it turns out to be violated by a solution of the Poisson equation. We may guess (and we show below) that the pixel boundaries motions are proportional to the integral over pixel boundaries of the perturbating electric field. These integrals would themselves be curlfree (as is the electric field) if the integration paths were all identical for serial and parallel boundaries. The fundamental anisotropy (e.g., from the ratio of a_{01} and a_{10}) indicates that this is not true at small distances, and questions the curlfree hypothesis applied to the pixel boundaries displacement field. Even without this anisotropy, we might also question whether the curlfree property of the continuous 3D electric field is precisely transferred to a curllike combination of finite differences over the pixel lattice.
Our practical approach is to promote the discrete a_{ĳ} scalar field into a discrete 2D vector field relying on electrostatics; namely, to propose to fit the geometrical parameters of a simple electrostatic model of the perturbating electric field to the data and to evaluate the 2D vector field we need from the model. We note that since the model reports actual pixel areas, the sum rule is by construction satisfied by the outcome of a fit.
4.2 The electrostatic model
We are interested in pixel boundary shifts under the influence of stored charge and here we sketch a firstorder perturbation scheme that allows us to predict these shifts from a few geometrical quantities. The first ingredient is the electric field sourced by the stored charges inside the sensor, which we refer above as the “perturbating” electric field (because it adds to the drift field that defines pixels). We assume that the field sourced by collected charges causes no rearrangement of charge within the bulk of the device. We model this perturbating field as sourced by a charge between two infinite grounded equipotentials figuring the light entrance window (on which the drift voltage is applied) and the parallel clock stripes. This is a textbook problem and using the image charge technique, we can write the corresponding potential as: (10)
where ρ^{2} = x^{2} + y^{2}, t is the sensor thickness, the source charge Q is located at (0,0, z_{q}), and the equipotential planes are at z = 0 (representing the parallel clock stripes) and z = t (the light entrance face of the sensor). We can check that ϕ(ρ, 0) = ϕ(ρ, t) = 0: in these cases, each positive term (first line) of the series has an exact opposite among the negative terms (second line), which is how the image technique works. This expression for the potential converges poorly with n for ρ much larger than t and alternative expressions should then be used (Pumplin 1969). For the fit, we do not need to evaluate the model at ρ > t.
The field that results from the normal operation of the sensor drives the charges into the pixel wells. The drift lines on the pixel boundaries have a peculiar point at which the field is null and the potential has a saddle point. This point is commonly located a few microns away from the clock stripes (at z = 0). These points are located at different heights above serial and parallel boundaries, and this difference contributes to the anisotropy of flat field smalldistance covariances. We display in Fig. 9 a contour crafted to relate the pixel boundary shifts to the source charge using Gauss’s theorem, applied to the total electric field, the drift field plus the perturbation introduced by the charge. This contour is made from four segments (numbers matching Fig. 9). Segment 1 corresponds to the perturbed drift line that separates two pixels and the integral of the transverse electric field is null. Segment 2 represents the short path between the unperturbed zerofield point and the perturbed zerofield point and is very close to a drift line; the integral of the electric field is then a second order quantity, which we ignore. Along segment 3, the integral of the unperturbed field vanishes because it is the unperturbed drift line; we are thus left with the integral of the perturbation field. Over segment 4, the integral of the field is the product of the drift field, E_{d}, times the boundary displacement, d, that we wish to estimate.
So, Gauss’s theorem is finally expressed as: (11)
where z_{0} is the zerofield point altitude, x_{b} and y_{b} are the coordinates of the unperturbed drift line, is the field transverse to the boundary sourced by the charge Q, E_{d} is the drift field (at the top), and the righthand side is the bulk charge contained inside the contour, because the Silicon material can contain residual impurities sourcing a small bulk electric charge. We do not know this righthand side, but we can assume that (to first order) it is proportional to d, our perturbation parameter. So, we eventually obtain: (12)
In other words, the boundary displacement is proportional to the integral of the perturbating (transverse) field over the unperturbed trajectory. This is usually called the Born approximation. In a chargefree material, the normalization is expressed as 1/E_{d}. The expression assumes that charges are produced when light enters the sensor, at z = t. For the reddest bands, we should instead account for the conversion depth of photons in the sensor, as we will do in Sect. 5.3. This expression should be averaged over the pixel side in the direction perpendicular to the Fig. 9.
In order to constrain the rhs of the expression, it is tempting to carry out the measurements at different values of E_{d} by changing the drift voltage applied to the entrance side of the sensor. However, changing the drift actually alters E_{d} but also z_{0} (for both flavors of boundaries), so measuring the sensor under different drift voltages would not help significantly at constraining the model. We may, in principle, predict the proportionality coefficient of Eq. (12) because the (empty CCD) charge density, thickness, applied drift voltage, and drift field are related by basic electrostatics. We nonetheless stick to fitting the global scale because precisely estimating the drift electric field involves at least the clock stripe sizes and their potentials during image integration. In this specific study, we have to fit the overall model scale because we were not able to find the drift and parallel clock potentials applied to the sensors.
We add some flexibility to the model by allowing source charges to be extended and we stick to a very crude model: a uniformly charged rectangle. This refinement, as compared to point sources, only influences the very first neighbors. So, the model eventually has six parameters (see Fig. 10): a global normalization factor, the height of the source charge z_{Q}, the height of the zerofield points over parallel and serial pixel boundaries z_{p} and Z_{s}, and the sides of the (uniform) rectangular source charge a and b.
The calculations of the model also require the thickness of the sensor (200 µm) and the pixel side (15 µm). For the practical implementation of the integrals from Eq. (12), we settle for integrating analytically the terms of the series similar to Eq. (10) for the electric field and emulating the rectangular source by splitting the charge into nine equally spaced point charges (only for n = −1,0,1). We sum the series of the expressions in Eq. (10) up to n = ±11. This seems to be sufficient for up to 10 pixels because increasing to n = 12 only changes outputs at the sixth decimal place. This electrostatic modeling approach was initially presented in Le Breton (2017), which we adapted Figs. 9 and 10 from.
Fig. 9 Contour to which we are going to apply Gauss’s theorem to relate the boundary shift (the length of segment 4) to the electric field sourced by the charge Q. 
Fig. 10 Schematic of the geometry of the electrostatic fit. 
Fig. 11 Data (top), fitted electrostatic model (middle), and difference (bottom). The sign of a_{00} has been flipped in the top plots, but not in the difference. We can note the clear difference of the two nearest neighbors, an anisotropy that is already visible in Fig. 6. We expect a poor fit because the data violates the sum rule, while the fit cannot. 
4.3 Fit of the electrostatic model to the data
We performed a leastsquares fit to the average area coefficient data, using the spread over channels to weight the squares. We first fit the average data obtained at the previous section. The data, fit, and difference are shown in Fig. 11. We can see that the measurements reach a signaltonoise ratio of about 1 at a distance of 8 to 9 pixels from the source. Beyond the small separations, the measurements exhibit a radial symmetry as expected from electrostatics and displayed as well by the model. The residuals do not average to zero and the data is larger than the model. This reflects that the data violates the sum rule but the fit cannot. The sum of the model area coefficients is −7 × 10^{−8} up to i, j < 10 (and tends to 0 with increasing bounds, so the contribution of unmeasured a_{ĳ} is 7 × 10^{−8}). Since the sum of measured a_{ij} is 7 × 10^{−8}, the data violates the sum rule by about 1.4 × 10^{−7}, which is more than 10% of a_{00}· The details of the electrostatic model cannot change the contribution of unmeasured aĳ significantly (at i > 9 or j > 9) and certainly cannot flip its sign. The excess of the data with respect to the model is concentrated on the 3 first serial pixels and the first parallel neighbor.
Some excess of variance and covariance along the serial direction is expected if a video signal experiences rapid gain variation (or “gain noise”): rapid gain changes contribute a variance component that scales as the square of signal level and, hence, artificially increase the value of a_{00} (see Eq. (5)). If gain fluctuations last longer than the time to read a pixel, they also bias covariances along the serial direction. Since the residuals seem to decay along the serial direction, possible gain variations have also to decay rapidly. They cannot be invoked to explain the excess on the first parallel neighbor because while serial pixels are read out microseconds apart, milliseconds separate neighboring lines.
We now attempt a fit that ignores the variance and the two next serial pixels, that is, a_{00}, a_{10} and a_{20}· The fit displayed in Fig. 12 seems acceptable: the residuals are at most 1% of the largest used coefficient a_{01}≃1.6 × 10^{−7}, and an even lower fraction of a_{00}The residuals exhibit some sort of lowlevel “chessboard pattern”, which is not entirely systematic. Assuming it is real, we have no proposition for its source, and we could not invent a small periodic variation of the physical size of pixels that would produce this 2 × 2 pixel pattern of the twopixel correlation functions.
We display the measured and fitted values, and their ratio in Fig. 13. The model reproduces the data well and gives some confidence that the model delivers sensible values for the data we decided not to use. We note that the decay of signal with distance is reproduced very well by the model over more than two orders of magnitude of signal level. In Pumplin (1969), it is shown that the largedistance decay of the perturbating electric field from a point charge in the sensor depends essentially exponentially on the inverse thickness of the sensor. The other geometrical parameters of the model cannot alter significantly the logarithmic slope of this decay at large distances. So, the evaluation of the model at large distances, needed to gauge the compliance of data to the sum rule, is a robust outcome of the model if the slopes of data and model match.
In Table 1, we display the measured and fitted largest area coefficients of the data and the two fitted models. We can see in particular the size of the discrepancies along the serial direction. We can now extract from the fitted model the pixel boundary displacement caused by a unit charge. This allows us to compute for a given science image, the accumulated boundary displacements resulting from all present charges. We will test the quality of the BF effect correction derived from the models on science images in the following section. We now question the curlfree hypothesis, and then discuss the gain noise issue.
Fig. 12 Data (top), fitted electrostatic model (middle), and difference (bottom). This fit ignores the three first serial pixels, marked with crosses in the bottom plot. Compared to Fig. 11, we see that the data excess in a_{01} has disappeared. We note that the color scale of the difference is zoomedin, as compared to Fig. 11. 
Fig. 13 Data and fit (left) as a function of distance, difference normalized to the model (right). The fit ignores the three first serial pixels, which are labeled with crosses in the rhs plot. The dashed line in the lhs plot separates the logarithmic and linear scale regions. 
Values of the largest area coefficients for the data and two fitted models.
4.4 Considering whether the pixel boundary displacement field is curlfree
In C18, the authors assume that the boundary displacement field is curlfree. Our definition of discrete curl is expressed as: (13)
In Fig. 14, we display the discrete 2D curl of the field, the displacement of one boundary (arbitrarily chosen as the parallel one, the serial one is in fact smaller) and their ratio. We used the electrostatic solution displayed in Fig. 12, which (by definition) satisfies the Poisson equation and, hence, has a curlfree 3D electric field. We can readily see that the curlfree assumption may lead to sizable errors (several tens of percent) in the correction.
One reason for the displacement field to be rotational is that drift paths along serial and parallel boundaries do not end at the same distance from the parallel clock stripes. The importance of this difference vanishes with distance because at large distances, the field varies less steeply with z than at short distances and, hence, the fractional contribution to the boundary displacement of the end of the drift path decays with distance. When it comes to science images, and especially when the image quality is good, the shortdistance boundary displacements dominate the effective correction for stars.
4.5 Gain noise
We measure a_{00} = −1.24 × 10^{−6}, while the fit of Fig. 12 indicates a_{00} = −l.37 × 10^{−6}. This means that the model predicts a variance that is smaller than the measured one, and the leading difference is quadratic in signal level. This difference of a_{00} values also contributes to the apparent curvature of the residuals displayed in Fig. 3: it roughly explains the curvature for C_{00} and C_{01} − but not for C_{10}.
While questioning the area coefficients along the serial direction, it is tempting to relate the excess we observe to charge transfer issues; however, we provide a few arguments against this explanation. First, we have measured and corrected deferred signals and we find that the linear dependence displayed in Fig. 5, causes, if uncorrected, a linear contribution to C_{10} (∝ µ); however, the a_{10} excess we observe corresponds to a quadratic contribution to C_{10} (∝ µ^{2}). Second, any mechanism relying on imperfect charge transfer will unavoidably affect C_{00} and C_{10} in comparable and opposite amounts, but the data in Table 1 indicate that the excess of a_{00} is about 20 times larger that the one of a_{10}, and they have the same sign. Third, serial charge transfer inefficiency is usually associated to localized and rare defects that is typically 1 or more frequently 0 along the serial register. The distributions of area coefficients (Fig. 8) do not seem quantized.
So, assuming that this difference of a_{00} values is due to gain variations during the read out, the relative gain variations should be about 3.5 × 10^{−4} rms. Gain variations of the HSC electronics have been studied during the HSC camera fabrication process and are reported in Miyatake et al. (2012). They are evaluated as ~2.4 × 10^{−5}, however, this is only considering the video chain and not the CCD (and perhaps in a context that is different from what the instrument actually faces).
The CCD readout chain of HSC relies on correlated double sampling, an approach that integrates the video signal during logical gates provided by the clocking system. The collected signal is hence vulnerable to fluctuations of the timing of the integration gates. We crosscorrelated the subimages from different channels of the same sensor in order to diagnose synchronous correlations, including variations of the gains (whatever the cause). We did not find any compelling signal for any sensor. We were hoping that crosscorrelations would deliver the size of gain fluctuations, or any other commonmode fluctuation of channel response, thus allowing us to subtract this contributions from flatfield statistics.
In summary, we are not able to provide a convincing cause of the sum rule violation. Because the electrostatic fit indicates an excess of fluctuations along the serial direction, we ignore three area coefficients. This excess is obviously puzzling, albeit compatible with rapid gain fluctuations (with some ad hoc time decay). Fortunately, we do not have to rely in any way on this hypothesis to use the electrostatic model.
Fig. 14 Values of the (discrete) curl, the top boundary displacement and their ratio. One can see that assuming a curlfree displacement field is acceptable only at distances larger than ~3 pixels. 
5 Correction of science images
5.1 Data set and reduction
We process the images of the Cosmos field from the ultradeep part of the Subaru Strategic Program, acquired between 2014 and 2019. The ultradeep part of the survey is geared at obtaining extremely deep images by coadding a large number of observations, as well as detecting highredshift supernovae and measuring their light curves. The sample of images covers the five bands: g, r, i, z and y of the camera, and also covers a broad range of image qualities. We use data acquired with the new i and r filters, which are called r2 and i2 in HSC parlance and we stick to these names. This image sample constitutes a representative playground for testing the quality of the brighterfatter correction on HSC science images.
The image reduction we perform is fairly standard for the brighterfatter correction. We first average the overscan and subtract it from the actual image data. We then correct each image segment for nonlinearity. Then we have to express the image in electrons in order to perform the BF correction, so we multiply each image segment by its corresponding gain (as determined when fitting the covariance data), perform the correction, and divide back each segment by the gain. We note that the BF correction has to deliver an image that has not been corrected by channel gains because our flats contain the gain differences and, hence, they should be applied to an unscaled image. We prefer to use flats that encode the relative gains because their evolution with time encodes a possible evolution of gains. Using approximate gains for the BF correction is clearly a secondorder issue, while using approximate relative gains on a sensor can cause artificial steps in the sky background at the channel boundaries.
The brighterfatter correction consists in computing the boundary shifts by summing all the actions of all image charges on all boundaries. For parallel boundaries, the shifts are expressed as: (14)
and similarly for serial boundaries. The factor of 1/2 accounts for the fact that source charges alter pixel shapes only once they have reached the pixel wells; thus, on average, during half of the integration time. In Eq. (1), the charge on the rhs refers to the charge accumulated so far during the integration and, hence, the timeaveraged boundary shift should be derived from the timeaveraged charge content. The charge to displace from one pixel to its neighbor is computed from the pixel boundary shift and the charge flowing over the same pixel boundary. We compute the latter as the average between the two pixels that share the boundary. We tried a quadratic order interpolator (involving 4 pixels), and did not find a decisive difference. We should note that the vast majority of our images are well sampled (the image quality is larger than 3 pixels) and sharper images could react differently. Ultimately, the correction of the BF effect in the images themselves requires that an accurate estimation of the charge at pixel boundaries can be devised.
Once the BF correction has been applied, we divide back each channel section by its gain and apply the flat field constructed from dome flats accumulated over about one month. For the y band, we have to subtract a fringe pattern constructed from a large set of science images.
On each individual CCD from each exposure, we run SExtractor (Bertin & Arnouts 1996) to detect objects and measure the Gaussian second moments of these objects from an unweighted 2D Gaussian fit to the light distribution. Because we do not integrate the Gaussian over pixels, we use the fast procedure described in Astier et al. (2013) to solve the normal equations. These second moments are similar if not identical to the “SDSS adaptive moments.” We identify stars in the image from the distribution of moments (see Astier et al. 2013). Because we have to accommodate significant variations of the PSF over the field of a single CCD, the cut that selects stars is broad enough to avoid rejecting genuine stars because of the BF effect. We get three moments per star: M_{XX}, M_{YY} and M_{XY}, where X and Y refer to the serial and parallel directions on the sensor. As an indicator of PSF size, we use: (15)
We then need color measurements of our stars, which we obtain by averaging fluxes over images and calibrating instrumental magnitudes over the common footprint of our star catalog and the field D2 from Betoule et al. (2013). This is not a critical step since only small color corrections are involved in what follows and we, hence, do not need colors measured to better than ~0.1 mag.
5.2 Processings and results
In order to assess the variation of star moments with their flux, we report the variation of star sizes with peak flux f_{max}, rather than with the flux itself. The variation with f_{max} is of practical interest, because f_{max} determines if a given star can enter into the PSF modeling. In order to isolate the variation of moments with peak flux we have to account for two other variations: the spatial variation in every exposure (mostly due to optics), and the color dependence of apparent star sizes. We remove the spatial dependence by fitting a secondorder 2D polynomial to the star moments measured on a given CCD in a given exposure. We then interpolate this crude model at the star position and study the residual of the measurement to the model. We model M_{XX} and M_{YY}3 independently. We apply a conservative cut at f_{max} < 35 000 ADU, in order to avoid any effect of sensor saturation. The dependence on f_{max} and color of these residuals are shown in Fig. 15, for a processing in which no BF correction was applied. We see that for stars selected at (M_{XX} + M_{YY})/2 ≃ 3.5 pix^{2}, the total increase of moments with f_{max} is ~0.13 pix^{2}, which is about 3.7% for stars varying from 0 to saturation. This slope depends on the selection of image quality. Figure 15 also indicates that the variation of apparent size with color in the z band are not considerably smaller than the ones with flux. Since the f_{max}color correlation coefficient is about −0.18 (in z band), we should account for colorinduced size variations.
In order to account for both the peak flux and color dependence of the moments, we regress the star moment residuals against both f_{max} and color: (16)
where A stands for X or Y. We then use α (i.e., the BF slope) as a colorindependent indicator of the BF effect intensity (both before and after correction). The linear correction in color can only work within a finite color range, and the selected color indicators and their range are provided in Table 2. For this fit, we ignore the measurements with f_{max}<5000 because the measured moments could be affected by a biased background estimation. Out of about 15 × 10^{6} star measurements, the flux and color cuts retain 2.4 × 10^{6} measurements contributing to the slope measurements that are reported here.
In order to study the quality of the BF correction on real science conditions, we bin stars in moments bins and restrict the range to (M_{XX} + M_{YY})/2 < 10, which roughly corresponds to an image quality of 1.35′′ FWHM. We do not apply a lower cut because the best image qualities are precious for at least cosmic shear. The bin selection operates on the expected moments at the location of the star rather than the measured ones, so that the quantity used to select the bins contents is statistically independent of flux and color.
The measurements of α (accounting for color corrections) for both CCD directions and in IQ_{2} bins are displayed in Fig. 16. We can notice an increase of the slopes with IQ_{2}, as well as the fact that the y band is notably less affected by the BF effect than the other bands. In the y band, the photons convert deeper in the sensor bulk and the charges experience a shorter drift path than for bluer bands, thus reducing the action of the perturbating electric field. For the y band, the starting point of the integral of Eq. (12) is noticeably smaller that the sensor thickness t.
For some time, we ignore this subtlety and correct all bands using a same model for all bands, derived from covariances measured in the g band. We first apply the correction derived from “model 1” (displayed in Fig. 11) and the residual BF slopes are displayed in Fig. 17. Although the BF slopes are significantly reduced with respect to the raw data (Fig. 16), the residual slopes are large at low IQ_{2}, indicating that the correction is underestimated at small distances. This correction leaves about 10% of the BF effect at the best image qualities and about 3% at the upper end.
We then apply the “model 2” correction (Fig. 12), namely, the one that ignores the three first suspicious measurements along the serial direction. The corresponding BF slopes are displayed in Fig. 18. The quality of the correction, in particular at the lowest IQ is improved. Ignoring the y band, we can interpret the figure as a global small overcorrection that leaves BF slopes about 30 times smaller than in the raw data.
Fig. 15 Differences between the star second moments and their expectations from the spatial smoothing, as a function of f_{max} (top) and color (bottom). We have selected stars with 3 < IQ_{2} < 4 pix^{2}, in the z band. The color range was chosen so that the color dependence is roughly linear. One can see that the slope of the top right plot is slightly steeper than the top left plot, reflecting the anisotropy of pixel covariances in flat fields. 
Color indicators and color ranges used for the various HSC bands.
Fig. 16 BF slopes (with color correction) in IQ_{2} bins for the five bands of HSC, without any BF correction, for the serial (top) and parallel (bottom) directions. We can note a roughly linear increase of the slopes with IQ_{2}, and that for the best image qualities, the slope for M_{YY} is larger than for M_{XX}, as expected from the anisotropy of the nearest neighbor a_{ij} (as seen in Fig. 11). 
Fig. 17 BF slopes (with color correction) in IQ_{2} bins for the five bands of HSC, with the BF correction derived from “model 1” (Fig. 11) which uses all the covariance measurements. One can note a residual BF slope essentially independent of image quality, for all bands, except y that deserves a specific treatment (detailed in Sect. 5.3). 
5.3 Brighterfatter correction for the yband
When the energy of photons becomes comparable to the silicon band gap, the absorption crosssection tends to vanish and silicon becomes transparent. The band gap is about E = 1.2 eV corresponding to λ = 1.1 µm. In the y band, the absorption length of photons becomes comparable to the sensor thickness and we can no longer assume that charges produced by converted photons drift all the way from the entrance window (z = t in the coordinates of Sect. 4.2). The flatfield data was acquired in g band, where the mean free path of photons is well below 1 µm, and the approximation of immediate conversion is adequate. This approximation works up to i band, may be questioned for z band and does not seem to apply to y band.
Once we have our electrostatic model, the prediction for boundary motions in y band can be computed by altering the lower bound of the integral of Eq. (12): (17)
where z_{c} is where the photon converts and k is the normalization determined when fitting the model to the area coefficient data (obtained in g band in our case). Since the integrand can vary rapidly with z, , it is unwise to evaluate these integrals at the average conversion depth of photons in the y band. We instead compute the average of d(z_{c}) for a realistic distribution of z_{c}: (18)
where dN/dz_{c} is the distribution of conversion points (and should integrate to unity over the integration domain), and Z_{b} refers to z_{s} or z_{p}, depending on which boundary we are considering. The distribution of conversion depths depends on photon absorption length and some assumption for the object spectrum. The latter may be regarded as inconvenient, but this is a small chromatic dependence of a small flux dependence and, hence, a secondorder effect.
We choose the spectrum of a color 0 (in AB magnitudes) object to compute the conversion depth distribution, noting that z−y = 0 lies inside the observed star color distribution, which has an average of <z−y>≃0.13. For the absorption length as a function of wavelength, we use the expression from Rajkanan et al. (1979) for Silicon at 173 K, which typically predicts 103 µm at λ = 950 urn. The conversion depth distributions for z and y bands is displayed in Fig. 19, and we can note that in y band, the conversions are spread over the whole thickness and the depth distribution is stable when changing the color by about the width of the z − y distribution of stars. In Fig. 20, we compare the electrostatic model 2 integrated over the full thickness with the same model for the y band, using Eq. (17), with the distribution of Fig. 19. The largest relative changes happen at large distances, where the perturbating electric field has a sizable contribution over most of the drift path.
Next, we apply this modified BF correction for the y band to the actual science data. We can see in Fig. 21 that the y band thus behaves similarly to the other bands, which indicates that the reduced BF effect in y likely originates in shorter drift paths in this very red band. This also indicates that for these CCDs, flat field correlations measured in a blue band can be used to predict (via a model) the BF correction for a red band. Alternatively, we could certainly consider measuring flatfield statistics in y band. We have not applied the same treatment to the z band, where the correction is much smaller than for y, and the need for reducing the BF correction is considerably less obvious. We will however apply the BF model for z band in our final image processing.
Fig. 18 BF slopes (with color correction) in IQ_{2} bins for the five bands of HSC, with the BF correction derived from “model 2” (Fig. 12) which uses all the covariance measurements but the three first serial pixels. We face a small over correction, and the BF slopes have been reduced by a factor of about 30 as compared to the raw values of Fig. 16. The y band still requires a specific treatment. 
Fig. 19 Distributions of conversion depth for a zero color object in the HSC CCDs. The y distribution is much flatter and leads to a globally reduced BF effect. Changing the color of the object by about the width of the our z−y distribution does not significantly alter the expected depth distribution in the y band. For the z band, most of the conversions occur at small depths, where the perturbating electric field is small. For the i2 band, the average conversion depth is about 10 µm. 
6 Discussion
6.1 Quality of the BF correction
The practical uses of PSF models crucially depend on the size of the model capable of faithfully representing the actual PSF size for faint objects. Two applications may come to mind: first, the measurement of faint supernovae for cosmological applications where fluxes have to be measured using PSF photometry and supernova fluxes are calibrated against the ones of bright stars; second, the measurement of galaxy shapes where one is interested in the intrinsic shape, namely, “before” it is smeared by the PSF. An inaccurate PSF size results in general in a biased shape. For both of these applications, the gauge is the difference of the PSF model size to the real size of a faint object, relative to the PSF size. The LSST Dark Energy Science Collaboration et al. (2018) set the maximum acceptable PSF size bias (δT/T) to 10^{−3} for the 10yr Rubin/LSST survey for all sources of PSF size bias, where T = M_{XX} + M_{YY} is the trace of the second moment matrix. So, in what follows, the PSF “size” is, in fact, the “area”. Regarding PSF photometry, the relative flux bias caused by PSF size bias reads δf/f = 1/2(δT/T). Shear measurements also require an accurate estimation of the PSF ellipticity. When surveys measure the same object with the same orientation of the sensors, sensorinduced ellipticities are transferred to shape measurements. For the Euclid mission the requirement on ellipticities transferred to objects is 5 × 10^{−5} of the rms (Table 1 of Cropper et al. 2013). In the context of the BF effect, we concentrate on the X/Y ellipticity, namely, (M_{XX} − M_{YY})/T. Cropper et al. (2013) also provide a requirement for the PSF size for Euclid, which reads 5 × 10^{−4} rms, comparable to the Rubin 10^{−3} bound.
In Fig. 22, we see that on corrected images, the relative PSF size difference is below 10^{−3} in all bands at all levels of image quality. In Fig. 23, we see that in the corrected images, the residual PSF ellipticity δ(M_{XX} − M_{YY})/T_{PSF} is mostly below 10^{−4} across the range of image quality and meets the 5 × 10^{−5} rms bound. We thus note that the need to omit three nearby serial pixels when fitting model 2 did not significantly degrade the quality of the BF correction in one direction. In Fig. 24, we display the slope of the radial fourth moment inaccuracies of the PSF, which also induces adverse effects on shear measurements, with a factor on the order of 1, as for the PSF size (see Zhang et al. 2021). The residuals are even smaller than for the PSF size. We note that our residual trends after correction all indicate an overcorrection of the BF effect and all this could thus be improved by adjusting the overall scale of the correction. We also note that these results were obtained in a blind way: we did not change the procedure after having seen them.
Gruen et al. (2015) evaluated their image correction for DECam in ways similar to ours (their Fig. 12): their method overcorrects the size of objects and leaves a negative BF slope, about −1/3 of the uncorrected one. This is certainly too large for a largescale cosmic shear survey. The corrected ellipticities seem significantly better (and probably small enough for a DESlike survey) but they are weakly affected by the BF effect in the test sample (as compared with the top plot of our Fig. 23). Mandelbaum et al. (2018) assessed the correction derived in C18 for HSC images: the (linear) size of corrected stars varies with magnitude by about 0.2% over three magnitudes (an eyeball estimation from their Fig. 6). Then, T_{PSF} would evolve twice as much and this is small enough for the analysis of the first year of the HSC survey. Regarding the ellipticities, trends with regard to the flux are not provided. For both DECam and HSC, the methods overcorrect the BF trend and both assessments ignore chromatic contributions to star size variations. If the color–flux correlation is similar to ours, correcting sizes for color decreases the apparent BF slope, hence degrading the performance of the BF correction in both instances.
We have insisted on the importance of higher order terms in Eq. (4) when analyzing the covariance curves. Those terms result from the integration over time of Eq. (1). Regarding the correction of the images, we did not carry out integration over time, and we just considered that half of the endofexposure image is a good approximation of the source of electrostatic distortions, hence, neglecting second order effects which are important for covariance curves. We anticipate that neglected second order effects in image correction should manifest themselves as structured residuals to the moments versus f_{max} linear fits. We analyzed these residuals in band and image quality bins, summing different bands at the same image quality, and we could not find any hint of departure from linearity. This is fortunate since accounting for next to leading order effects in image correction is more difficult than for the shape of covariance curves.
Fig. 20 Area coefficients for “model 2” (displayed in Fig. 12) based on the same model with shorter drift paths meant to describe the y band (center) and the relative difference at the bottom. The differences tend to increase with distance to the source. 
Fig. 21 BF slopes computed after correcting images with model 2, with the modifications described above for y band. With respect to Fig. 18, only the y data has changed. We may note that all bands now behave in a similar way. 
Fig. 22 Relative change of the PSF size between an average PSF star and a faint object, as a function of T_{PSF}/2. For the PSF size extimator, we use as well the trace of the second moment matrix of the PSF. The top plot refers to the raw data, the bottom one to the data corrected by model 2 with the shortened drift paths in y band. We assume that an average PSF star peaks at 1/3 of the saturation. 
Fig. 23 Change of the PSF ellipticity e_{PSF} = (M_{XX} − M_{YY})/T_{PSF} between an average PSF star and a faint object, as a function of T_{PSF}/2. The top plot refers to the raw data, the bottom one to the data corrected by model 2 with the shortened drift paths in y band. We assume that an average PSF star peaks at 1/3 of the saturation. 
6.2 Possible further developments
While the obtained performance of the BF correction seems sufficient, we may consider potential avenues for improvements. First, in our analysis, the overall normalization of the BF correction model primarily depends on a_{01}, which is the best measured area coefficient actually used in the electrostatic fit (the model 2 fit). This coefficient is slightly biased by an erroneous value for a_{00} through the µ^{3} terms in the covariance fit to a level compatible with the small overcorrection we are facing. It may then seem legitimate to actually tune the overall normalization of the correction in order to bring the average BF slope to 0, at least for the best observing conditions.
Second, the BF correction we implemented is the outcome of an analysis where three area coefficients of the BF effect (the three first serial measurements, including the largest coefficient) had to be ignored because they could not be accommodated by an electrostatic model. Those coefficients had to be derived from the measurements of other area coefficients through an electrostatic modeling. We might then anticipate that for a camera that does not suffer from this noise bias (which is first manifested by a violation of the sum rule), a better correction model could be constructed. In particular, an analysis of test data of the Vera Rubin Observatory camera indicates that the sum rule is satisfied on the integrated instrument.
One appealing – and very different – approach could be to determine the electrostatic model from the science data by fitting the pixel distortion pattern that makes the PSF homothetic. The electrostatic modeling could probably be computed sufficiently rapidly to be inserted into a PSF modeling fitting loop. The concerns about the precision measurement of the twopoint function of flat fields would become pointless.
Fig. 24 Relative change of the radial fourth moment as a function of image quality for all HSC bands, between a faint source and an average PSF star, before (top) and after (bottom) BF correction. We assumed that an average PSF star peaks at 1/3 of the saturation. 
6.3 Computer codes
Our code is split in three different parts: code to measure covariances and fit the covariance curves, code to perform the electrostatic fit, and code to process the images and in particular corrects those for the BF effect. Our public repository^{1} contains the python code for the two first steps. The code that measures covariances did not evolve significantly since it was developed for A19. We publish here the electrostatic modeling code for the first time. Our image correction code is fairly straightforward and would probably run much faster with the convolution in Eq. (14) computed in Fourier space. The parts of the analysis that are not in the repository can be made available upon request.
Acknowledgements
We are indebt to the Subaru Telescope technical staff that very efficiently operates the observatory and its instruments. HSC images are made available on the SMOKA server (https://smoka.nao.ac.jp/index.jsp) that we have used extensively. We perform all our reductions and store our results at the Centre de Calcul de l’IN2P3 (https://cc.in2p3.fr), a computing facility of CNRS. This work benefited from useful discussions with N. Suzuki (LBL) and N. Yasuda (IPMU), and our colleagues from the Paris team. The manuscript eventually benefited from excellent suggestions from our referee.
References
 Antilogus, P., Astier, P., Doherty, P., Guyonnet, A., & Regnault, N. 2014, J. Instrum., 9, C3048 [Google Scholar]
 Astier, P., El Hage, P., Guy, J., et al. 2013, A&A, 557, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Astier, P., Antilogus, P., Juramy, C., et al. 2019, A&A, 629, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Betoule, M., Marriner, J., Regnault, N., et al. 2013, A&A, 552, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Coulton, W. R., Armstrong, R., Smith, K. M., Lupton, R. H., & Spergel, D. N. 2018, Astron. J., 155, 258 [NASA ADS] [CrossRef] [Google Scholar]
 Cropper, M., Hoekstra, H., Kitching, T., et al. 2013, MNRAS, 431, 3103 [Google Scholar]
 Downing, M., Baade, B., Sinclaire, P., Deiries, S., & Christen, F. 2006, Proc. SPIE, 6276, 9 [Google Scholar]
 Gruen, D., Bernstein, G., Jarvis, M., et al. 2015, J. Instrum., 10, C05032 [CrossRef] [Google Scholar]
 Guyonnet, A., Astier, P., Antilogus, P., Regnault, N., & Doherty, P. 2015, A&A, 575, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lage, C., Bradshaw, A., Anthony Tyson, J., & The LSST Dark Energy Science Collaboration 2021, J. Appl. Phys., 130, 164502 [NASA ADS] [CrossRef] [Google Scholar]
 Le Breton, R. 2017, Ph.D. Thesis, Université Pierre et Marie Curie, Paris VI, France [Google Scholar]
 Mandelbaum, R., Miyatake, H., Hamana, T., et al. 2018, PASJ, 70, S25 [Google Scholar]
 Miyatake, H., Aihara, H., Fujimori, H., et al. 2012, Phys. Proc., 37, 1413 [NASA ADS] [CrossRef] [Google Scholar]
 Miyazaki, S., Komiyama, Y., Kawanomoto, S., et al. 2018, PASJ, 70, S1 [NASA ADS] [Google Scholar]
 Pumplin, J. 1969, Am. J. Phys., 37, 737 [NASA ADS] [CrossRef] [Google Scholar]
 Rajkanan, K., Singh, R., & Shewchun, J. 1979, SolidState Electron., 22, 793 [NASA ADS] [CrossRef] [Google Scholar]
 Rasmussen, A., Guyonnet, A., Lage, C., et al. 2016, Proc. SPIE 9915, 99151A [NASA ADS] [CrossRef] [Google Scholar]
 The LSST Dark Energy Science Collaboration (Mandelbaum, R., et al.) 2018, ArXiv eprints [arXiv:1809.01669] [Google Scholar]
 Zhang, T., Mandelbaum, R., & The LSST Dark Energy Science Collaboration 2021, MNRAS, 510, 1978 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 CCD signal values as a function of time for a single channel during a sequence of uniform exposures illuminated by a lamp in the dome. The intensities are deliberately not ordered in a monotonie way, so that the drift of the illumination system can be separated from the nonlinearity. The sequence was designed by N. Yasuda (IPMU). 

In the text 
Fig. 2 Spline correction from the minimization of Eq. (8) for four different channels, normalized to its average. Since this is presumably due to lamp variations, it should be similar for all channels. The maximum difference is in the 10^{−4} range. 

In the text 
Fig. 3 Fit residuals as a function of date (top) and signal level (bottorn) for all channels of all sensors. There are no obvious needs for a more sophisticated model as a function of date nor signal level. The rms residual is about 3 × 10^{−4}µ. 

In the text 
Fig. 4 Values of measured quadratic non linearity k (in Eq. (8)) in the focal plane (top) and their distribution (bottom). Most of the channels cluster around 10^{−7}, which corresponds to a correction of ~0.4% around maximum CCD signal. 

In the text 
Fig. 5 Nearest serial neighbor deferred signals for the four channels of sensor 42. The values are small and vary almost linearly with the signal. The model continuous lines are the cubic polynomial model that we use for correction. 

In the text 
Fig. 6 Measurement and fits of the variance (C_{00}) and the two nearest neighbor covariances (C_{10} and C_{01}) for the four channels of sensor 42. The channels have been offset for visual convenience. The decay of C_{00}/µ with µ, originally noted in Downing et al. (2006), is the most obvious effect of the BF effect. The spread of measurements is compatible with shot noise. 

In the text 
Fig. 7 Residuals of fits of variance and covariance measurements, all expressed in electrons using the gains from the fits, so that the data from different sensors can be averaged. The red signs are binned averages of the residuals from individual fits. 

In the text 
Fig. 8 Distribution of area coefficients over the ~400 hundred channels of HSC (blue) and the averages over each sensor (orange). All a_{ij} coefficients are expressed in el^{−1}. We may note that the trends for channels and sensors are similar, which indicates some homogeneity within sensors (except for Σ a, where the scatter is dominated by shot noise). Σ a is the sum for −10 < i, j < 10. 

In the text 
Fig. 9 Contour to which we are going to apply Gauss’s theorem to relate the boundary shift (the length of segment 4) to the electric field sourced by the charge Q. 

In the text 
Fig. 10 Schematic of the geometry of the electrostatic fit. 

In the text 
Fig. 11 Data (top), fitted electrostatic model (middle), and difference (bottom). The sign of a_{00} has been flipped in the top plots, but not in the difference. We can note the clear difference of the two nearest neighbors, an anisotropy that is already visible in Fig. 6. We expect a poor fit because the data violates the sum rule, while the fit cannot. 

In the text 
Fig. 12 Data (top), fitted electrostatic model (middle), and difference (bottom). This fit ignores the three first serial pixels, marked with crosses in the bottom plot. Compared to Fig. 11, we see that the data excess in a_{01} has disappeared. We note that the color scale of the difference is zoomedin, as compared to Fig. 11. 

In the text 
Fig. 13 Data and fit (left) as a function of distance, difference normalized to the model (right). The fit ignores the three first serial pixels, which are labeled with crosses in the rhs plot. The dashed line in the lhs plot separates the logarithmic and linear scale regions. 

In the text 
Fig. 14 Values of the (discrete) curl, the top boundary displacement and their ratio. One can see that assuming a curlfree displacement field is acceptable only at distances larger than ~3 pixels. 

In the text 
Fig. 15 Differences between the star second moments and their expectations from the spatial smoothing, as a function of f_{max} (top) and color (bottom). We have selected stars with 3 < IQ_{2} < 4 pix^{2}, in the z band. The color range was chosen so that the color dependence is roughly linear. One can see that the slope of the top right plot is slightly steeper than the top left plot, reflecting the anisotropy of pixel covariances in flat fields. 

In the text 
Fig. 16 BF slopes (with color correction) in IQ_{2} bins for the five bands of HSC, without any BF correction, for the serial (top) and parallel (bottom) directions. We can note a roughly linear increase of the slopes with IQ_{2}, and that for the best image qualities, the slope for M_{YY} is larger than for M_{XX}, as expected from the anisotropy of the nearest neighbor a_{ij} (as seen in Fig. 11). 

In the text 
Fig. 17 BF slopes (with color correction) in IQ_{2} bins for the five bands of HSC, with the BF correction derived from “model 1” (Fig. 11) which uses all the covariance measurements. One can note a residual BF slope essentially independent of image quality, for all bands, except y that deserves a specific treatment (detailed in Sect. 5.3). 

In the text 
Fig. 18 BF slopes (with color correction) in IQ_{2} bins for the five bands of HSC, with the BF correction derived from “model 2” (Fig. 12) which uses all the covariance measurements but the three first serial pixels. We face a small over correction, and the BF slopes have been reduced by a factor of about 30 as compared to the raw values of Fig. 16. The y band still requires a specific treatment. 

In the text 
Fig. 19 Distributions of conversion depth for a zero color object in the HSC CCDs. The y distribution is much flatter and leads to a globally reduced BF effect. Changing the color of the object by about the width of the our z−y distribution does not significantly alter the expected depth distribution in the y band. For the z band, most of the conversions occur at small depths, where the perturbating electric field is small. For the i2 band, the average conversion depth is about 10 µm. 

In the text 
Fig. 20 Area coefficients for “model 2” (displayed in Fig. 12) based on the same model with shorter drift paths meant to describe the y band (center) and the relative difference at the bottom. The differences tend to increase with distance to the source. 

In the text 
Fig. 21 BF slopes computed after correcting images with model 2, with the modifications described above for y band. With respect to Fig. 18, only the y data has changed. We may note that all bands now behave in a similar way. 

In the text 
Fig. 22 Relative change of the PSF size between an average PSF star and a faint object, as a function of T_{PSF}/2. For the PSF size extimator, we use as well the trace of the second moment matrix of the PSF. The top plot refers to the raw data, the bottom one to the data corrected by model 2 with the shortened drift paths in y band. We assume that an average PSF star peaks at 1/3 of the saturation. 

In the text 
Fig. 23 Change of the PSF ellipticity e_{PSF} = (M_{XX} − M_{YY})/T_{PSF} between an average PSF star and a faint object, as a function of T_{PSF}/2. The top plot refers to the raw data, the bottom one to the data corrected by model 2 with the shortened drift paths in y band. We assume that an average PSF star peaks at 1/3 of the saturation. 

In the text 
Fig. 24 Relative change of the radial fourth moment as a function of image quality for all HSC bands, between a faint source and an average PSF star, before (top) and after (bottom) BF correction. We assumed that an average PSF star peaks at 1/3 of the saturation. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.