High ?delity point-spread function retrieval in the presence of
    electrostatic, hysteretic pixel response
    Andrew Rasmussen
    a
    , Augustin Guyonnet
    b
    , Craig Lage
    c
    , Pierre Antilogus
    b
    , Pierre Astier
    b
    ,
    Peter Doherty
    d
    , Kirk Gilmore
    a
    , Ivan Kotov
    e
    , Robert Lupton
    f
    , Andrei Nomerotski
    e
    , Paul
    O'Connor
    e
    , Christopher Stubbs
    d
    , Anthony Tyson
    c
    , and Christopher Walter
    g
    a
    SLAC National Accelerator Laboratory, Menlo Park, CA, United States
    b
    LPHNE/IN2P3/CNRS, UPMC, France
    c
    University of California at Davis, Davis, CA, United States
    d
    Harvard University, Cambridge, MA, United States
    e
    Brookhaven National Laboratory, Upton, NY, United States
    f
    Princeton University, Princeton, NJ, United States
    g
    Duke University, Durham, NC, United States
    ABSTRACT
    We employ electrostatic conversion drift calculations to match CCD pixel signal covariances observed in ?at ?eld
    exposures acquired using candidate sensor devices for the LSST Camera.
    1
    ,
    2
    We thus constrain pixel geometry
    distortions present at the end of integration, based on signal images recorded. We use available data from several
    operational voltage parameter settings to validate our understanding. Our primary goal is to optimize ?ux point-
    spread function (FPSF) estimation quantitatively, and thereby minimize sensor-induced errors which may limit
    performance in precision astronomy applications. We consider alternative compensation scenarios that will take
    maximum advantage of our understanding of this underlying mechanism in data processing pipelines currently
    under development.
    To quantitatively capture the pixel response in high-contrast/high dynamic range operational extrema, we
    propose herein some straightforward laboratory tests that involve altering the time order of source illumination
    on sensors, within individual test exposures. Hence the word hysteretic in the title of this paper.
    Keywords: CCDs, drift ?elds, charge collection, ?at ?eld statistics, pixel size variation, imaging nonlinearities,
    brighter-fatter e?ect, instrument signature removal
    1. INTRODUCTION
    The literature already includes several instances of how ?at ?eld correlations may be used to correct or compensate
    astronomical data for the so called brighter-fatter (BF) e?ect.
    3
    {
    5
    In a separate work, Niemi et al.
    6
    opted to favor the information available from direct, focused spot measure-
    ments over the indirect information from ?at ?eld correlations, and that the ?ux level dependence to measured
    spot sizes was described as an intrinsic CCD PSF that depends on intensity and wavelength only. A key piece
    of the puzzle that con?icts with this picture is that the BF e?ect nearly vanishes when a single pixel's center is
    illuminated with a spot of sub-pixel diameter:
    7
    the instrument's signature that contributes to systematic errors
    in turn must depend on the incident ?ux distribution at the sensor's entrance window as well as the instanta-
    neous recorded signal distribution as illumination progresses. Consequently, we argue that contributions from
    the instrument cannot be separated from the contributions of the incident ?ux, unless the integration of the
    recorded image is also considered in the process.
    In this contribution, we approach the issue by calculating changes to pixel areas based on families of electro-
    static solutions to Poisson's equation in the drift region.
    8
    A separate e?ort, not discussed in detail here, applies
    a Poisson solver to sensor's semiconductor properties informed by detailed fabrication steps and lithographic
    information, provided by the vendor.
    9
    Send correspondence to A.R. { E-mail: arasmus@slac.stanford.edu; Tel: +1 650 926 2794; Fax: +1 650 926 5566
    Full electronic version of this paper with color ?gures is available on arXiv.org [astro-ph.IM]

    2. DYNAMIC CHANGES TO PIXEL AREAS AND THEIR RELATIONSHIP TO
    MEASURED COVARIANCES
    Direct pixel boundary calculations, if they indeed reproduce available characterization data, are likely to be
    preferable to pixel border shift models because they provide a two-dimensional point-to-pixel partition (they do
    not over- or under-count area elements), and also o?er chromaticity information (pixel boundaries as a function
    of conversion depth) in the dynamic pixel geometric response. Direct boundary calculations properly handle
    built-in nonlinearities in pixel geometries as aggressor signals, due to accumulated conversions, approach full
    well. In the context of the rolled-up model we describe, the recorded signal distribution of an image is used to
    calculate the self-consistent, distorted pixel boundaries in e?ect at the end of the exposure. Better constraints
    on the incident ?ux distribution would naturally result from knowledge of those boundaries.
    In the following, we provide some calculations that extend work discussed in the aforementioned references
    to connect pixel correlations to their theoretical signal level dependence, whether for ?at ?eld or focused spot
    applications. We compute pixel area response to aggressor signal level.
    2.1 Flat ?eld statistical ?uctuations and signal expectation values as a function of lag
    In this discussion, the aggressor is a statistical ?uctuation ? in a recorded ?at ?eld image, that occurs about the
    mean ? and induces changes in neighboring pixels' areas. A neighboring pixel is indicated by its lag from the
    aggressor using two indices, ij, where i and j are the lags along the serial and parallel directions, respectively.
    Because we calculate area variations due to drift during collection and not transfer statistics due to trap pop-
    ulations and channel occupancy, full descriptions of area variations are captured by considering only positive i
    and j. Correspondingly, ij = 00 indicates the aggressor pixel and q
    00
    is the charge signal accumulated there.
    Considering the direct aggressor{victim channel only, a nearby pixel with lag ij has an area at the end of
    integration:
    ?lna
    ij
    (q
    00
    = ?+?j?) ˇ
    d ln a
    ij
    dq
    00
    ?;
    (1)
    and on average this pixel would contain a signal level that is systematically biased by the area distortion. While
    the area distortion is zero at the beginning of the integration (?a
    ij
    = 0 @ t = 0), it should be ?nite by integration
    end. Averaged over all possible trajectories ?a
    ij
    (t), the in?uence of the statistical ?uctuation in pixel ij = 00
    on the bias h?q
    ij
    i is readily isolated:
    h?q
    ij
    (q
    00
    j?)i =
    ?
    2
    (exp(?lna
    ij
    (q
    00
    j?))? 1) ˇ
    1
    2
    ??
    d ln a
    ij
    dq
    00
    =??
    d ln a?
    ij
    dq
    00
    (2)
    where for convenience, we also use the exposure time averaged pixel area a?
    ij
    . The expression for the covariance
    Cov
    ij
    may also be simpli?ed by using the same approximation, and the variance Var ? Cov
    00
    appears:
    ?
    Cov
    ij
    =
    P
    kl
    ?
    kl
    P
    h?q
    k + i;l + j
    i
    kl
    1
    =
    ?
    2
    P
    kl
    ?
    kl
    (exp(?ln
    P
    a
    k + i;l + j
    (q
    kl
    j?)) ? 1)
    kl
    1
    (3)
    ˇ
    ?
    2
    P
    kl
    ?
    2
    kl
    d ln a
    k + i;l + j
    dq
    kl
    P
    kl
    1
    =
    ?
    2
    Cov
    00
    d ln a
    ij
    dq
    00
    :
    (4)
    The correlations are then the covariances divided by the variance, and following through with the same approx-
    imation, they are proportional to the area response d ln a
    ij
    =dq
    00
    and the scaling term ?=2:
    Corr
    ij
    ?
    Cov
    ij
    Cov
    00
    ˇ
    ?
    2
    d ln a
    ij
    dq
    00
    ?
    For consistency in nomenclature, we de?ne the zero lag covariance to be equal to the variance:
    Cov
    00
    =
    P
    kl
    ?
    k
    P
    l
    ?q
    k +0 ;l +0
    kl
    1
    =
    P
    kl
    ?
    2
    kl
    P
    kl
    1
    ? Var:

    Now irrespective of ?nite correlations in the ?at ?elds, Poisson statistics are recovered by re-binning images
    (in the case of data frames), or in the case of ?nite pixel area distortions:
    P
    kl
    ?
    ?
    kl
    P
    ij
    ?
    k + i;l + j
    ?
    P
    kl
    P
    ij
    1
    =
    P
    kl
    ?
    2
    kl
    P
    kl
    1
    +
    P
    kl
    P
    ij 6=00
    ?
    kl
    ?
    k + i;l + j
    P
    kl
    P
    ij 6=00
    1
    = Cov
    00
    +
    X
    ij 6=00
    Cov
    ij
    ? ?:
    Because the expressions for Cov
    ij
    are symmetric under exchange i ! ? i and j ! ? j, the above expression may
    be further simpli?ed to include only the unique quantities Cov
    ij
    for i ? 0;j ? 0:
    ?=
    X
    ij
    Cov
    ij
    = Cov
    00
    + 2
    X
    i ?0 ;j ?0
    (1? ?
    0 i
    ?
    0 j
    )(1+?
    ij
    )Cov
    ij
    ;
    (5)
    where ?
    ij
    is the Kronecker delta. Finally, area is conserved, such that any area lost (or gained) by pixel ij = 00
    is recovered (or ceded) by others:
    X
    ij
    ?a
    ij
    = 0 = (exp(?lna
    00
    (q
    00
    j?))? 1)+
    X
    ij 6=00
    (exp(?lna
    ij
    (q
    00
    j?))? 1):
    (6)
    Equations
    5
    &
    6
    provide a way to express the shape of the mean-variance curve. Solving for Cov
    00
    :
    Cov
    00
    = ?
    0
    @
    1?
    X
    ij 6=00
    Cov
    ij
    1
    A
    =
    ?
    ?
    1+
    P
    ij 6=00
    Covij
    Cov
    00
    ?
    (7)
    ˇ
    ?
    1+
    ?
    2
    P
    ij 6=00
    d ln a
    ij
    dq
    00
    =
    ?
    1?
    ?
    2
    d ln a
    00
    dq
    00
    =?
    X
    1
    n =0
    ?
    ?
    2
    d ln a
    00
    dq
    00
    ?
    n
    :
    (8)
    Evidently, the observed deviation from Poissonian behavior in the mean-variance curve indicates
    d ln a
    00
    dq
    00
    <0
    3
    ,
    8
    and the approximate relation is valid if j
    ?
    2
    d ln a
    00
    q
    00
    j ˝ 1. In general however, Equation
    7
    should remain valid
    independent of this assumption. The general form of this is Equation
    9
    . Routine, accurate gain determination
    may be enabled by ?tting functions of this form to gain-variance measurement pairs that show this curvature.
    10
    2.2 Application to measured covariances
    We apply the equations above to the speci?c regime where statistical ?uctuations ? are much smaller than (and
    tied to) the ?at ?eld ?ux ?. Recall that in the Poisson limit, ?
    2
    ? ?, but this is apparently not correct in
    actual photon transfer curves. High-quality ?at ?eld data sets can be used to generate a pattern of lag (ij)
    speci?c covariances Cov
    ij
    , which are in turn converted into correlations via Corr
    ij
    ? Cov
    ij
    =Cov
    00
    . Generally,
    Cov
    ij 6=00
    scale as ?
    2
    , while Corr
    ij 6=00
    scale as ?. From the modeling side, a statistical ?uctuation translates to
    an aggressor amplitude p? which in turn produces the pattern of area distortions ?a
    ij
    (p?j?) that govern the biases
    in the expression for Cov
    ij
    (Eq.
    2
    ). Taking the rms exposure averaged aggressor to be p? ? z
    chan
    ?q
    e
    , we write
    the following:
    Cov
    ij
    = ???a
    ij
    (p?)
    Cov
    00
    = ?
    2
    =
    ?
    1?
    ?
    ?
    ?a
    00
    (p?)
    (9)
    Corr
    ij
    =
    ?
    ?
    ?a
    ij
    (p?) = (? ? ??a
    00
    (p?)) ?a
    ij
    (p?):
    (10)
    In the above, Corr
    ij
    , ? and ? are measurements and ?a
    ij
    (p?) are compiled from results of the drift calculation.
    A measured correlation pattern Corr
    ij
    , together with estimates for ? and ? for a speci?c ?at ?eld illumination
    were used to ?t an electrostatic drift model for its undetermined parameters. The data were acquired from a

    candidate sensor prototype for the LSST Camera manufactured by e2v: it was described previously in Ref.
    4
    .
    ˜
    2
    was minimized using the Nelder-Mead method with results shown in Figure
    1
    (left). The best-?t parameter
    list is given in Table
    1
    . Four free parameters were jointly estimated in the process. The area distortion model
    appears to provide enough detail to reproduce the measured, anisotropic correlation pattern and fallo? with
    separation. In addition to constraining magnitudes of the periodic barrier dipole moments, the ?t also favors a
    speci?c value for the impurity concentration in the silicon bulk, N
    a
    . The goodness of ?t was acceptable when
    using estimated ? and ?, so the process completed without invoking a gain error parameter in Eq.
    10
    .
    A secondary result of this ?tting procedure is a value for the channel depth. For ? ˘ 65ke
    ?
    , we ?nd
    z
    chan
    = p=? (?q
    e
    ) ˇ 2:37?m. This result is shown below (x
    4.1
    ) to constrain other physical properties of the sensor.
    1
    10
    100
    10
    −4
    10
    −3
    0.01
    lag specific correlation coefficients Corr
    ij
    radial lag (i
    2
    +j
    2
    ) [pix
    2
    ]
    measured_vs_computed_corr_ij_comparision.qdp
    electrostatic computed (surface)
    measured correlations
    electrostatic computed (@ saddle locus)
    0
    5×10
    −4
    10
    −3
    1.5×10
    −3
    0
    5×10
    −4
    10
    −3
    1.5×10
    −3
    parallel coordinate [cm]
    serial coordinate [cm]
    contour_map_surface_deep.qdp
    p = 2.8189 p
    0
    δ
    ln a
    00
    = −0.306626
    δ
    ln a
    00
    = −0.130224
    δ
    ln a
    10
    = 0.0094869
    δ
    ln a
    10
    = 0.0239552
    δ
    ln a
    01
    = 0.032465
    δ
    ln a
    01
    = 0.021333
    δ
    ln a
    11
    = 0.0097589
    δ
    ln a
    11
    = 0.0047034
    Figure 1. Pixel area variations. Left: a comparison of the best-?t electrostatic drift model to measured correlation
    coe?cients for a speci?c ?at ?eld illumination. Measurements of Corr
    ij
    are the dots (black) with error bars, ?lled
    triangles (red) are for shallow conversions, and open triangles (blue) are limiting cases for the conversions occurring near
    the saddle locus along the drift lines. Any observed chromaticity in the ?at ?eld correlations should produce numbers
    that lie somewhere between the shallow and deep limits corresponding to the same lag ij. Best-?t parameters, along with
    quantities that a?ect the interpretation of this model are given in Table
    1
    . Right: a graphical representation of the same
    electrostatic model's pixel boundary distortions, for the nearest neighbors (the aggressor is denoted lag ij = 00). The solid
    (red) lines show the pixel boundaries for cold electrons very close to the backside surface, while the dotted (blue) lines
    show a two-dimensional projection of the saddle locus, where adjacent drift lines diverge to feed the channels belonging to
    adjacent pixels. The dashed (black) lines, together with the plot frame, show positions of the undistorted pixel boundaries,
    i.e., for zero aggressor amplitude. For the purpose of this graphic, the aggressor dipole moment amplitude p was increased
    by a factor of ˘ 500 as compared to the best-?t on the left (while holding all other parameters ?xed).
    3. TUNABLE ELECTROSTATIC INFLUENCES TO PIXEL AREAS
    We have a starting point for more detailed modeling. It is a relatively straightforward task to reproduce ?xed
    pattern features seen in the sensors that would be categorized as cosmetics { but can in fact be traced to pixel
    area distortions. In previous work we have demonstrated success in reproducing observed features seen in ?at
    ?eld illumination via electrostatic modeling: edge rollo?, midline charge redistribution, bright and dark column
    pairs identi?ed as tearing features,
    11
    and bamboo,
    12
    each with self-consistent pixel shifts and elongations that
    accompany the pixel area distortions revealed in the ?at ?elds.
    8
    The important di?erence between this and prior
    work is that the statistical properties of the ?at ?eld illumination, Cov
    ij
    , are also reproduced at the same time.

    Table 1. Parameter list for the best-?t electrostatic drift model (for cold carriers)
    parameter
    value
    units
    comments
    N
    a
    1:11?10
    12
    cm
    ? 3
    acceptor density in depleted Si
    t
    Si
    100
    ?m
    sensor thickness (?xed)
    BSS
    ? 78
    Volts
    backside bias (?xed)
    a
    ˘
    cs
    12:407
    ˘
    0
    b
    channel stop 2-D dipole moment
    ˘
    ck
    2:6425
    ˘
    0
    b
    clock barrier 2-D dipole moment
    p?
    0:0057208
    p
    0
    c
    aggressor dipole moment
    d
    ?
    65230
    e
    ?
    mean signal level in ?at
    ?
    2
    58429
    (e
    ?
    )
    2
    variance in ?at
    d
    a
    constrained by measured X-ray di?usion variation with BSS on a
    similar device
    b
    ˘
    0
    ?10
    ? 6
    q
    e
    c
    p
    0
    ? 10
    5
    ?mq
    e
    d
    exposure averaged, rms aggressor moment is p? ? z
    chan
    ? q
    e
    The best-?t electrostatic parameters given in Table
    1
    can represent ?ducial performance, and small changes
    in their values will consequently a?ect the dynamic, hysteretic response of the sensor. For example, a change in
    clock rail voltage di?erential would induce a proportional change in ˘
    ck
    , a change in backside bias (BSS) would
    induce a change in E~
    BD
    (z) according to Appendix
    B
    , Equation
    13
    , and operation with tearing present (betrayed
    by darker column pairs straddling segment boundaries in ?at ?eld response) would alter ˘
    cs
    11
    and boost the
    BF e?ect due to reduced barriers there. Moreover, parameterization of the aggressor moment p, its dependence
    on signal ? and the geometric response of pixel boundaries (see x
    3
    ) can provide an informed process by which
    recorded images can be used to constrain ?ux distributions incident at the sensor entrance window. A check for
    how Corr
    ij
    vary with BSS between measurement and calculation is given in Figure
    2
    . It demonstrates, via a
    blind test, that counterintuitive dependencies in the measured Corr
    ij
    are reproduced using the simple, far-?eld
    approximate, multipole expansion drift model that explicitly satis?es Poisson's equation. The Dirichlet solution
    approach
    9
    to describe the sensor's photosensitive volume from the clocks and channel stops toward the backside
    surface, may still be required to accurately model or understand other details of sensor operation, but these are
    not de?ned or addressed here. This serves as a proof of concept that the drift model may already be adopted
    for more widespread application in data analysis to reduce sensor systematics.
    3.1 The brighter-fatter template
    Several e?orts could bene?t from this detailed and robust modeling. These include astronomical data reduction
    pipelines and simulation tools (based on image, table, or ray tracing). E?cient implementations require that
    results of the time consuming electrostatic drift calculation be ported to faster simulation or analysis frame-
    works. This was discussed previously (Ref.
    8
    , x5.3) but is brie?y summarized and expanded for completeness in
    Appendix
    D
    . A single calculation result, in the format of a BF template could be applied over broader range of
    aggressor amplitudes, indeed up to the canonical full well depth for these sensors, using the linear perturbation
    model. The work addressed in x
    3.2
    tests this notion.
    3.2 Pixel area dependence on aggressor dipole moment
    By evaluating pixel border distortions for di?erent aggressor strengths p and comparing against a linear scaling
    of the BF template, we essentially test the validity of the following expressions (cf. Appendix
    D
    ), where p
    t
    is the
    aggressor dipole moment for which the template was generated:
    ?~c
    t ( i;j )
    k
    (p) =
    ?
    p
    p
    t
    ?
    ?~c
    t ( i;j )
    k
    (p
    t
    ); 8k
    ?d
    t ( i;j )
    l
    (p) =
    ?
    p
    p
    t
    ?
    ?d
    t ( i;j )
    l
    (p
    t
    ); 8l

    d ln(a
    10
    )/dq
    00
    ij=01
    2.6E−7
    2.8E−7
    3.0E−7
    3.2E−7
    d ln(a
    10
    )/dq
    00
    ij=10
    0.0E−7
    1.0E−7
    2.0E−7
    −1.0E−7
    0
    20
    40
    60
    80
    d ln(a
    ij
    )/dq
    00
    aligned |BSS| (V)
    ij=02
    ij=20
    0.0E−7
    0.5E−7
    1.0E−7
    1.5E−7
    Figure 2. A blind, qualitative comparison between measured correlation coe?cients (left) and computed area distortions
    (right) as the backside bias is varied. Left: this plot was reproduced from Fig. 8 of Guyonnet et al.
    4
    and represents ?at ?eld
    correlations measured for several lags ij 2 f01; 10; 02; 20g for ? ˘ 100 ke
    ?
    . Right: computed per-electron area distortions
    for the same lag selection, for the model summarized in Table
    1
    , ?t for jBSSj = 70V (\aligned" jBSSj = ? 78V
    12
    ). While
    even the Corr
    ij
    pattern shown here (for jBSSj = 70V) isn't entirely consistent with the Corr
    ij
    pattern shown in Figure
    1
    (left) collected for ? ˘ 65 ke
    ?
    , this imperfect comparison shows that counterintuitive results seen in the data are readily
    reproduced in the drift model. These include changes in sign of Corr
    ij
    and their derivatives with respect to jBSSj. The
    ordinate scales on the right hand plots were adjusted to compare directly to the corresponding plots on the left.
    where the instantaneous p=p
    t
    may be as large as twice the ratio between the canonical full well (100ke
    ?
    ) to
    the rms statistical ?uctuation level ? given in Table
    1
    . In the present case, p
    max
    ˘ (2FW=?)p? ˘ 826p?, many
    times the ?uctuation levels sampled by ?at ?eld correlations. Figure
    3
    provides a quantitative comparison of
    the evolution of pixel area variations with aggressor p and shows that the scaled template/linear perturbation
    approach su?ers from signi?cant error. It should be emphasized that the apparent (per lag) nonlinearities
    revealed here are entirely inaccessible to con?rmation via ?at ?eld correlations, unless the operational sensor full
    well is closer to 10
    7
    ke
    ?
    : a natural consequence of ?at ?eld statistics we used to probe the BF e?ect in the ?rst
    place.
    4. POSSIBLE NONLINEARITIES MEASURABLE IN FLAT FIELD COVARIANCES
    Although the nonlinearities suggested in x
    3.2
    are not measurable using ?at ?elds for the reasons described, we
    consider other terms here that could alter parameterization of the electrostatic elements in the drift model, or
    would otherwise a?ect charge accumulation in the receiving channels beyond the cold carrier approximation.
    4.1 Aggressor dipole moment dependence on accumulated conversions
    The connection between charge collected and e?ective aggressor dipole moment in this drift model is the depth
    of the buried channel, z
    chan
    (cf. x
    2.2
    ). A strawman model for the buried channel depth may be constructed
    by superposing two contributors of the axial component of the drift ?eld, one from the bound charge density
    (depleted n-type Si) and one from the in?uence of the free charges in the conductive polysilicon gates, which can
    be treated as an image charge. With z
    0
    equal to z
    chan
    in the limit of an unpopulated channel, we do not allow

    0
    0.1
    0.2
    0.3
    0.4
    pixel_area_evolution.qdp
    −Δa
    00
    (p)
    drift calc
    simple linear
    linear perturbation
    0
    0.02
    0.04
    +Δa
    01
    (p)
    0
    5×10
    −3
    0.01
    0.015
    +Δa
    10
    (p)
    0
    1
    2
    3
    4
    5
    0
    5×10
    −3
    0.01
    0.015
    aggressor dipole moment p [p
    0
    ]
    +Δa
    11
    (p)
    0
    1
    2
    3
    4
    5
    1
    1.2
    1.4
    Δ
    a (drift calc) /
    Δ
    a (linear perturbation)
    aggressor pixel dipole moment p [p
    0
    ]
    pixel_area_ratio_evolution.qdp
    ij=00
    ij=01
    ij=10
    ij=11
    (p
    FW
    )
    0
    5×10
    −4
    10
    −3
    1.5×10
    −3
    2×10
    −3
    0
    5×10
    −4
    10
    −3
    1.5×10
    −3
    2×10
    −3
    parallel transfer coordinate [cm]
    serial transfer coordinate [cm]
    compare_pixel_contour_maps.qdp
    ij=00
    ij=01
    ij=10
    ij=11
    direct calculation (p/p
    0
    = 4.65)
    p/p
    0
    = 5.72E−3 template
    scaled to p/p
    0
    =4.65
    dir:
    δ
    ln a
    01
    = 5.442E−2
    sca:
    δ
    ln a
    01
    = 4.229E−2
    dir:
    δ
    ln a
    11
    = 1.607E−2
    sca:
    δ
    ln a
    11
    = 1.352E−2
    dir:
    δ
    ln a
    00
    = −0.5776
    sca:
    δ
    ln a
    00
    = −0.4955
    dir:
    δ
    ln a
    10
    = 1.551E−2
    dir:
    δ
    ln a
    10
    = 1.061E−2
    Figure 3.
    A study of pixel area distortion dependence on aggressor amplitude p. Upper left: for each of 4 lags,
    ij 2 f00; 01; 10; 11g, measures of the pixel area distortions described in the text are given { the linear perturbation using
    the BF template (dotted line), the direct drift calculation (solid line), and a linear area model where the area distortion
    scales directly with aggressor amplitude p (dashed line). All curves in each plot coincide where p
    t
    = p? = 0:00572 p
    0
    from
    the ?t to ?at ?eld correlations, and are tabulated out to p ˘ 4:7p
    0
    , which corresponds approximately to 100ke
    ?
    in the
    aggressor. In the case of ij = 00 (top tier), a negative sign has been applied to ?a for a better comparison to the other
    curves. Upper right: ratios are plotted to compare the direct drift calculation to the linear perturbation approach using
    the template for each of the lags considered here. While the two approaches mismatch by about 10% for ij = 00 near
    the full well scale, the discrepancies are typically much larger, as much as 50%, for the other lags there. Bottom: pixel
    boundary maps for an aggressor at full well, 100ke
    ?
    (p = 4:65p
    0
    ) to further compare the two methods for these lags.
    The ˘50% discrepancy between the two methods shown for ij = 10 may appear counterintuitive because the area gained
    appears larger than it really is. In fact, ?a
    10
    ˇ ? 0:027?a
    00
    , based on the drift calculation shown here.

    for any spatial extent to the charge cloud and treat only its centroid:
    E~(z
    chan
    ) ? z^ =
    ?
    E~
    B
    (z
    chan
    jz
    0
    )+E~
    image
    (z
    chan
    jz
    0
    )
    ?
    ? z^ = 0
    z^? E~
    B
    (z
    chan
    jz
    0
    ) =
    N
    d
    2?
    0
    ?
    Si
    (z
    chan
    ? z
    0
    )
    z^? E~
    image
    (z
    chan
    jz
    0
    ) =
    ?+?
    4ˇ?
    0
    ?
    Si
    (2z
    chan
    )
    2
    z~ ? z
    chan
    =z
    0
    z~
    3
    ? z~
    2
    = ?
    ?+?
    8ˇN
    d
    z
    3
    0
    :
    The preceding equations describe two real solutions for ?nite z~ (one stable solution for each charge sign) in terms
    of physical properties of the channel, channel depth for zero signal z
    0
    and the signal occupancy ? + ? (mean plus
    aggressor). As the channel becomes populated, the solutions draw closer to one another until they coalesce as an
    in?ection point. Beyond this there is no ?nite, real solution and the channel should empty as quickly as it ?lls
    to expose gate structure layers with conversions. We may de?ne this measure of the full well, FW ˇ
    32 ˇ
    27
    N
    d
    z
    3
    0
    (to potentially constrain physical parameters of the channel using linearity measurements). Solutions for z~ for
    ?+?˝
    32 ˇ
    27
    N
    d
    z
    3
    0
    follow the approximately linear trend:
    z~? 1ˇ?
    ?+?
    8ˇN
    d
    z
    3
    0
    =?
    ?+?
    3:3 ? 10
    6
    ?
    N
    d
    10
    16
    cm
    3
    ?
    ? 1
    ?
    z
    0
    2:36 ?m
    ?
    ?3
    :
    (11)
    If indeed N
    d
    ˘ 10
    16
    cm
    ? 3
    , then we should expect a weakening trend of the coupling between signal and aggressor
    dipole moment of about 3% per 100ke
    ?
    ; and a proportionally weakened coupling if N
    d
    is smaller. This scaling
    of the coupling term z
    chan
    contributes additional, signi?cant detail (depending on N
    d
    ) at the scale shown in
    Figure
    3
    and inclusion of this e?ect is treated in x
    4.3
    .
    4.2 Aggressor driven modi?cation to drift time and di?usion contribution to correlations
    It has been proposed
    4
    ,
    5
    that the presence of collected conversions at the channel can attenuate the axial term
    in the drift ?eld and could help to explain some of the discrepancies seen between an electrostatic treatment
    of pixel borders and observed correlations. Guyonnet
    4
    (x5.3) investigated this further to place upper limits
    (< 4%) on the contribution by di?usion - longer collection times - to the total the BF e?ect for focused spots.
    In the context of the drift model described here, it seems di?cult to accurately isolate longer collection times
    from the e?ect of distorted pixel boundaries that accompany a reduced electric ?eld near the backside window.
    Indeed, an initial survey of temperature dependence of the ?at ?eld correlations produced inconclusive results.
    In any case, the notion we examine here is that a large aggressor could raise collection times feeding into its
    own pixel, e?ectively causing an additional redistribution of charges to neighboring pixels, while the neighboring
    pixels do not reciprocate via the same mechanism (they retain their nominal collection times). If signi?cant, this
    mechanism could be included by adding another term on the right hand side of Equation
    2
    .
    The same drift calculation used to determine pixel boundaries for cold carriers is used to estimate the e?ect
    of the accumulated conversions on subsequent conversions' collection times (and di?usion). Figure
    4
    shows this
    dependence as if the trajectory for cold carriers can be used to compute a drift time and di?usion, for carriers
    with temperature matching the substrate's temperature of T = 173 K. Presence of the barriers and the aggressor
    are clearly seen as cusps in the ˙(~x) ?eld sampled by the crosscuts shown. When averaged over these linear
    traces, positions nominally tied to the central pixel have a modestly increased ˙?
    00
    that is about 1 part in 400
    greater than for positions not tied to the central pixel. We further estimate the net redistributive e?ect when the
    central pixel has a larger di?usion ˙ than its neighbors. With ˙
    nom
    ˘ 0:4 (pixels), we sample and average over
    all possible Gaussian centroids contained within the pixel to compute the expected contribution to that pixel,
    which (for 0:2 < ˙ < 0:6) works out to:
    hC
    00
    (˙)i =
    Z
    1
    0
    dx
    0
    Z
    1
    0
    dy
    0
    Z
    1
    0
    dx
    Z
    1
    0
    dy
    1
    2ˇ˙
    2
    exp
    ?
    ?
    (x? x
    0
    )
    2
    +(y? y
    0
    )
    2
    2
    ?
    ˇ 1:0216347 ? 1:7486435˙ + 0:8989589˙
    2
    :

    4
    4.1
    4.2
    diffusion_assist_bf.qdp
    σ
    along 45
    o
    diagonal
    −40
    −30
    −20
    −10
    0
    10
    20
    30
    40
    4
    4.1
    4.2
    position along segment [μm]
    σ
    along serial transfer axis
    σ
    along parallel transfer axis
    diffusion
    σ
    [
    μ
    m] for fiducial performance
    Figure 4. A calculation to estimate the di?usion assisted contribution to the BF e?ect. Notionally, in addition to the
    pixel boundaries getting distorted, the presence of the exposure averaged aggressor dipole p? can slow down the carrier
    drift toward the channel to increase collection time, boosting the di?usion parameter ˙. For ?ducial sensor performance
    (parameters listed in Table
    1
    but also T = 173 K), the drift calculation was used to sample the launch position dependence
    for di?usion. Left: three linear traces are given { one diagonal and two along the address axes { to show ˙ vs. position.
    In each case, position zero is the location of an exposure averaged aggressor dipole p? that corresponds to a ?nal, recorded
    signal of 100ke
    ?
    . Presence of the electrostatic elements (˘
    cs
    , ˘
    ck
    and p?) in speci?c locations print through to produce
    the modulations shown. They collectively cause slowdowns and de?ections along the ?eld lines. Right: a graphic
    representation of how the linear traces were simulated (sections A-A' , B-B' and C-C' ) where lines of the grid represent
    the barriers that form pixel boundaries. An approximate averaging over these linear traces reveal only a modest (0.25%)
    increase in ˙ for positions that feed the central pixel containing the aggressor over the neighbors.
    The net redistributive e?ect from lag ij = 00 to neighboring pixels would be ?lnq
    00
    ˘ ?˙
    @
    lnhC
    00
    i, or about
    ? 2:2?10
    ? 3
    (for ˙ ˇ 0:4p and ?˙ ˘ ˙=400). If this \missing" signal were recovered and divided evenly between
    the nearest four neighboring pixels on average, the largest in?uence would be seen in lag ij = 10, because the
    exposure averaged area increase there (cf. Figure
    3
    ) is small: ˘ 8?10
    ? 3
    , and would impart a ˘ 7%, di?usion
    assisted excess over the nominal, pixel area distortion-driven Corr
    10
    . The term is far less consequential for lag
    ij = 01 (˘ 2% excess) and for lag ij = 00 (˘ 1% excess). The ?t to the correlations presented above { Table
    1
    and Figure
    1
    { did not include this as a separate term that would tend to dilute the observed anisotropy between
    Corr
    10
    and Corr
    01
    in the model. It is unknown at this point whether its inclusion would improve the quality
    of the ?t, or if tighter constraints on the Corr
    ij
    measurements would motivate its inclusion. In any case, the
    term's scaling and dependence mimics that of the area distortion model described in Equation
    2
    and is currently
    absorbed in the deterministic, detailed pixel boundary model. A separate treatment of this e?ect may be more
    important for devices with smaller backside ?eld strengths, or if these devices were operated with smaller jBSSj.
    We are encouraged by what is supported by Figure
    2
    as a reasonable correspondence between measurement and
    calculation shown for Corr
    ij
    vs. jBSSj { that separate inclusion of this term may be unwarranted.
    4.3 Predictions for nonlinearities in Corr
    ij
    (?)
    In the sections above, several terms were described that a?ect the evolution of the dynamic pixel area distortion
    model that are currently not well constrained by test stand measurements. These generally in?uence the detailed
    relationship between accumulated ?ux and aggressor dipole moment (x
    4.1
    ). In the current case, we have a pattern
    of Corr
    ij
    measured at a single ?at ?eld level ? and variance ?
    2
    , a drift model parameter list that reasonably
    reproduces the correlation pattern when pixel area distortions are mapped via Equation
    10
    , and a set of pixel
    boundary distortions generated for a selection of aggressor levels p. Assuming for the time being that the gain
    was accurately determined and that errors on ? and ? are negligible, we investigate how the computed area
    distortion patterns can be mapped onto an observable set of parameters. By ?rst assuming a value for N
    d
    ,

    expressions for the aggressor p? and z~ are used to determine channel depth in the zero signal limit by solving
    iteratively
    z
    0
    =
    p?
    F
    =?
    F
    1?
    ?
    F
    16 ˇ N
    d
    z
    3
    0
    ;
    where ?
    F
    and ?
    F
    were used in the ?tting procedure (Eq.
    10
    ) used to determine p?
    F
    . Drift calculations for other
    aggressors p?
    k
    produce the pixel area distortions ?a
    ij
    (p?
    k
    ). Self-consistent mean & variance pairs (?
    0
    k
    ;?
    0
    k
    ) are then
    calculated, also iteratively, using the equations:
    ?
    0
    k
    =
    p?
    k
    =z
    0
    ?
    1?
    ?
    0
    k
    16 ˇ N
    d
    z
    3
    0
    ?
    =
    A
    1? B?
    0
    k
    ;
    ?
    0
    k
    =
    ?
    ?
    0
    k
    ?
    2
    0
    B
    B
    @
    1 ?
    ?
    0
    k
    ?
    0
    k
    ?a
    00
    (p?
    k
    )
    ?
    1?
    ?
    0
    k
    16 ˇ N
    d
    z
    3
    0
    ?
    ?
    1?
    ?
    F
    16 ˇ N
    d
    z
    3
    0
    ?
    1
    C
    C
    A
    =
    ?
    ?
    0
    k
    ?
    2
    1?
    ?
    0
    k
    ?
    0
    k
    ?a
    00
    (p?
    k
    )
    1? B?
    0
    k
    1? B?
    F
    !
    :
    This process allows us to predict the detailed shape of the mean-variance curve as well as the mean-correlation
    curves for speci?c lags ij. A similar procedure would be used to allow for (and constrain) a gain error in a
    non-degenerate way. This isn't discussed here, but is straightforward to implement, given an additional set of
    approximate, (?,?) pairs derived from photon transfer curves.
    Figure
    5
    provides a family of curves that predict the signal dependence of the observable quantities ?
    2
    ? Var
    and Corr
    ij
    for di?erent assumptions of N
    d
    , again by assuming that an accurate gain was determined to produce
    ?
    F
    and ?
    F
    . Self-consistent values for z
    0
    for the assumed values of N
    d
    are also given. It turns out that while the
    mean-variance curve is relatively insensitive to the nonlinearities considered in x3.3, the signal level dependence of
    the Corr
    ij
    may the most straightforward indicator for an evolution in the coupling between signal and aggressor.
    5. FALSIFIABLE TESTS OF LINEAR EXTRAPOLATIONS OF THE
    BRIGHTER-FATTER TEMPLATE: HIGH-CONTRAST LABORATORY TESTS
    The situation we described above is that we predict signi?cant departures from linear perturbation models when
    we deal with real, high-contrast/high dynamic range data. Di?culty arises from not being able to sample high-
    contrast/high dynamic range conditions using ?at ?elds alone: the aggressor scale available tops out near the
    square root of the full well. It's surprising, then, that the linear perturbation methods used in astronomical
    pixel data pipelines can correct 90% of these dynamical e?ects: that only 10% of the initial BF e?ect remains
    uncorrected
    13
    after the compensation strategy is applied. This is based only on tracking a single \width"
    parameter for the PSF's intensity dependence, and does not at all capture the platykurtic distortions to the PSF
    pro?le that result from the boundary distortion mechanism. Nonlinear terms due to the variable channel depth
    appear to reduce the BF e?ect by about 6% averaged over the exposure for full well (if N
    d
    = 5 ? 10
    15
    cm
    ? 3
    );
    the direct drift calculations suggest that the linear perturbation underestimates the BF e?ect, anisotropically,
    by 10 ? 20% for high-contrast/high dynamic range exposures reaching the same full well in the pixels receiving
    the highest ?ux. We expect that sensors using smaller electric ?elds or longer drift distances than these LSST
    candidate devices should show correspondingly larger complications.
    We consider some lab measurements that could be performed to test the drift model { and the linear pertur-
    bation template methodology. Because the template is based only on a single aggressor level p^
    F
    =p
    0
    = 0:00572
    (for ?
    F
    = 65ke
    ?
    & ?
    F
    = 242e
    ?
    ), the next-to-leading order terms described in x
    4.3
    are not carried. We describe
    a receiving pixel array with geometric parameters that evolve with time, as the exposure progresses toward full
    well. Ratios of images (long vs. short exposure), incremental di?erence images (subtraction of images with
    adjacent exposure times), etc., are simulated. Pixel areas at the end of exposure are also recorded { to predict
    the e?ect of a ?at ?eld \?ash" at the end of the high-contrast exposure. Di?erences between (high-contrast +
    ?ash) and (high-contrast only) can be used to measure the pixel area ?eld across the array at the end of high-
    contrast exposure. Deviations from these predictions may be interpreted as a superposition of the nonlinearities

    0
    5×10
    4
    10
    5
    1.5×10
    5
    2×10
    5
    −0.2
    −0.1
    0
    0.1
    Cov
    00
    relative deviation from linear
    flat field level [μ
    k
    ]
    area_distortions_to_observables_cov00.qdp
    μ
    F
    (Cov
    00
    k
    )/Cov
    00
    F
    ))*(μ
    F
    k
    ) − 1
    5×10
    4
    10
    5
    1.5×10
    5
    2×10
    5
    −0.05
    0
    relative change in observables (finite N
    d
    ):(large N
    d
    ) − 1
    flat field level [μ
    k
    ]
    area_distortions_to_observables.qdp
    μ
    F
    Cov
    00
    Corr
    ij
    N
    d
    = 1E16 (z
    0
    = 2.39μm)
    N
    d
    = 5E15 (z
    0
    = 2.41μm)
    (N
    d
    =5E21: z
    0
    = 2.37μm)
    Figure 5. Aggressor induced, pixel area distortion calculations (?a
    ij
    (p?
    k
    )) mapped into observable ?at ?eld statistics to
    predict their dependence on ?at ?eld level ?
    k
    , while allowing for an unknown donor density in the channel (N
    d
    ). Left:
    the expected nonlinear term in the variance, ?
    2
    k
    ? Cov
    00
    (?
    k
    ), evaluated for a ?xed channel depth (equivalently, a large
    N
    d
    ). Any existing error in the gain determination should result in photon transfer data plotted o? of the locus shown
    here. The level called out as ?
    F
    is the ?ux level at which the correlation pattern was ?t with an adequate pixel area
    distortion template drawn using the drift calculation. Similar nonlinear terms in calculations of Corr
    ij
    were below the
    percent level with no clear trend, so these were not shown here. Right: corrections to the Cov
    00
    and Corr
    ij
    mappings
    that result from a channel depth that evolves with signal level. Two values for donor density in the channel are shown:
    N
    d
    = 1 ? 10
    16
    and 5 ? 10
    15
    cm
    ? 3
    . Mapping corrections for Cov
    00
    tend to be limited to the 2% level, but corrections to
    Corr
    ij
    are independent of lag ij, and would tend to show small positive o?sets in Corr
    ij
    vs. ?, depending on the sampling
    values that are available.
    described, namely details of the lag-dependent pixel area evolution with aggressor, combined with details of
    the relationship between signal and aggressor, and the evolution of the channel depth. Figures
    6
    and
    7
    show
    predictions for a focused, Gaussian spot and for interferometric two-slit fringe projections, respectively. NB:
    these predictions do not include the detailed, position-dependent drift times for cold carriers (x
    4.2
    ) and only use
    the one-to-one position-to-pixel mapping.
    6. CONCLUSIONS
    By mapping and matching drift calculation results, representing aggressor-victim pixel area distortions, to mea-
    sured ?at ?eld correlations, we whittled down and simpli?ed the BF e?ect to just three electrostatic parameters
    of the drift ?eld that were not already well constrained by X-ray charge cloud di?usion, for this CCD imaging
    sensor. A fourth, self-interaction parameter was also determined, and de?nes the e?ciency by which conversions
    collected in the channel can break the symmetry of the drift ?eld to reposition and distort pixel boundaries
    encountered for subsequent conversions.
    A compact, digital form for the pixel boundary mapping kernel, or Green's function, was used to predict
    results for some high-contrast, high dynamic range illuminations of the sensor that could be tested in the
    laboratory. We expect measurable deviations from these predictions for at least two reasons: (1) that the pixel
    area distortion variation with aggressor amplitude does not increase as the linear, perturbation theory would
    predict, and (2) that the self-interaction parameter coupling should evolve with accumulated conversions as an
    exposure progresses.
    These laboratory measurements, when performed, may provide a basis by which our quantitative understand-
    ing of the BF mechanism can be extended to include the high-contrast, high dynamic range domain needed for
    precision astronomical corrections.

    15
    20
    25
    30
    35
    0
    2×10
    4
    4×10
    4
    cts per pixel
    pixel address [10μm]
    psf_profile_compare.qdp
    first half
    second half (vs. parallel)
    second half (vs. serial)
    15
    20
    25
    30
    35
    0.8
    0.85
    0.9
    0.95
    1
    area [pix
    2
    ]
    pixel parallel address [10μm]
    psf_pixel_area_evolution.qdp
    end of half exposure
    end of full exposure
    Figure 6. A high-contrast, high dynamic range simulation of an idealized PSF using linear perturbation of the template.
    The illumination used is an isotropic Gaussian with 0.7
    00
    FWHM centered on the geometric midpoint of four pixels.
    Integration continues until full well is reached, which corresponds to AB˘15.2 for LSST's r-band in a 15 second exposure.
    Left: a comparison of the accumulated image after the ?rst half of the exposure to the additional accumulation during
    the second half. These are linear traces that pass through the PSF centroid. The BF e?ect is seen by comparing linear
    traces for the second half of the integration against that for the ?rst half. The total number of counts in the traces for
    the second half are typically about 5% lower than for the ?rst half, because counts are also distributed perpendicularly
    to the trace. Right: A comparison of the pixel areas resulting from the PSF integration half-way through and after
    completion. These can be used to estimate structure in the sky background contribution, and errors in the PSF pro?le if
    sky background is subtracted (without using this information). Laboratory data obtained to reproduce these results may
    reveal the di?erences predicted in x
    4.3
    and the limitations of the linear perturbation model used here.
    APPENDIX A. COMPARISONS TO SIMILAR EXPRESSIONS IN THE
    LITERATURE
    Because we aim to generalize and extend what's already in the literature, it should be useful to review here
    expressions of similar quantities that have already been published. Note that in the preceding equations, a
    ij
    indicates pixel areas [e.g., cm
    2
    ], not pixel boundary shift coe?cients perpendicular to boundary axes [e.g.,
    pixel/carrier], as in some of the expressions below.
    A.0.1 Antilogus et al. (2014)
    In their x4.2 treatment of charge responsive pixel boundaries applied to ?at ?eld correlations, the authors don't
    distinguish between instantaneous pixel boundary shifts and those shifts implied by statistics of the recorded
    image { in other words, the exposure averaged boundary shifts. We ?nd a factor of 2 discrepancy between their
    equations and ours if the former interpretation is followed, but perhaps no discrepancy with the latter. Upon
    comparing their equations 4.14 and 4.15 against our approximate expressions (Eqs.
    4
    and
    8
    respectively):
    Cov
    ?
    Q
    0
    i;j
    ; Q
    0
    0 ; 0
    ?
    = 4?V
    X
    X
    a
    X
    i;j
    Cov
    ij
    ˇ
    ?
    2
    Cov
    00
    d ln a
    ij
    dq
    00
    and
    Cov
    ?
    Q
    0
    0 ; 0
    ; Q
    0
    0 ; 0
    ?
    =V +4V?
    X
    X
    a
    X
    0 ; 0
    Cov
    00
    ˇ
    ?
    1?
    ?
    2
    d ln a
    00
    dq
    00
    ˇ?
    ?
    1+
    ?
    2
    d ln a
    00
    dq
    00
    ?
    +O
    ?
    3
    ?
    1
    2
    d ln a
    00
    dq
    00
    ?
    2
    !
    :

    0
    5×10
    4
    10
    5
    @ end of step
    accum. counts
    period1_fringe_integration.qdp
    0.99
    1
    1.01
    rel. 1st step
    counts increment
    0
    20
    40
    60
    80
    100
    0.98
    1
    1.02
    @ end of step
    pixel area
    parallel pixel address [10μm]
    0
    5×10
    4
    10
    5
    @ end of step
    accum. counts
    period2_fringe_integration.qdp
    0.95
    1
    1.05
    rel. 1st step
    counts increment
    0
    20
    40
    60
    80
    100
    0.9 0.95 1 1.05 1.1
    @ end of step
    pixel area
    parallel pixel address [10μm]
    Figure 7. A high-contrast, high dynamic range simulation of fringe projector illumination to accompany laboratory
    measurements. This simulation used linear perturbation of the template and the illumination was chosen to imitate the
    available parameter space. In both cases, the peak:valley ratio is set to 3, and the ?nal maximum counts accumulated is
    near the full well depth of 100 ke
    ?
    =pixel. Postage stamp images to the left show the image accumulated at readout time.
    Upper plots: a fringe period of 28.8 pixels with orientation 82
    ?
    ; Lower plots: a fringe period of 7.6 pixels with orientation
    140
    ?
    . On the right, for each fringe calculation, the plots show the counts accumulated in 11 steps (top tiers), the 10 ratios
    of the incremental counts accumulated divided by counts accumulated in the ?rst step (middle tiers), and the 11 area
    curves computed at the end of each step (bottom tiers). The perturbed area ?elds may be validated by subtracting \fringe
    only" exposures from \fringe+?at" exposures, which may be the most direct way to probe a dynamic pixel distortion
    mechanism. The perturbed area ?eld would be easier to measure, with greater contrast, if the \?at" part of the exposure
    could be applied after the \fringe" part, rather than in a simultaneous exposure.

    We identify the equivalent instantaneous area distortion coe?cients d ln a
    kl
    =dq
    00
    that are a factor of 2 larger
    than the exposure averaged area distortion coe?cients, valid for ?at ?eld applications, at least:
    1
    2
    d ln a
    ij
    dq
    00
    ?
    d ln a?
    ij
    dq
    00
    ˘
    V
    Eq: 4 : 14
    Cov
    00
    4
    X
    X
    a
    X
    i;j
    1
    2
    d ln a
    00
    dq
    00
    ?
    d ln a?
    00
    dq
    00
    ˘
    V
    Eq: 4 : 15
    ?
    4
    X
    X
    a
    X
    0 ; 0
    +O
    ?
    ?
    1
    2
    d ln a
    00
    dq
    00
    ?
    2
    !
    where it appears that the V in their equations 4.14 and 4.15 may have di?erent de?nitions.
    y
    A.0.2 Guyonnet et al. (2015)
    In their x5.2 parameterization of pixel size variations as a function of ?ux, this paper uses largely the same
    notation as in,
    3
    except that boundary shift coe?cients a
    X
    i;j
    are de?ned di?erently by a factor of 4, such that the
    area distortion coe?cients are identi?ed, but the de?nitions for their V in equations 16 and 17 remain to be
    aligned:
    z
    1
    2
    d ln a
    ij
    dq
    00
    ?
    d ln a?
    ij
    dq
    00
    ˘
    V
    Eq: 16
    Cov
    00
    X
    X
    a
    X
    i;j
    1
    2
    d ln a
    00
    dq
    00
    ?
    d ln a?
    00
    dq
    00
    ˘
    V
    Eq: 17
    ?
    X
    X
    a
    X
    0 ; 0
    +O
    ?
    ?
    1
    2
    d ln a
    00
    dq
    00
    ?
    2
    !
    :
    A.0.3 Gruen et al. (2015)
    In their x3.2 (Flat ?eld covariances), the authors cite
    3
    but yet write down slightly di?erent expressions for the
    covariances as a function of lag, and choose a di?erent normalization for the pixel boundary shifts. As before,
    we compare their Equations 3.7 and 3.8 to approximate our approximate forms of our Eqs.
    4
    and
    8
    , respectively:
    Cov(Q
    00
    ;Q
    ij
    ) = 2?
    2
    X
    X = T;B;L;R
    a
    X
    ij
    Cov
    ij
    ˇ
    ?
    2
    Cov
    00
    d ln a
    ij
    dq
    00
    and
    ?Var(Q
    00
    )=Var? ?=? 4?
    2
    ?
    a
    R
    1 ; 0
    +a
    T
    0 ; 1
    ?
    = +2?
    2
    X
    X = T;B;L;R
    a
    X
    00
    Cov
    00
    ? ?ˇ
    1
    2
    ?
    2
    dlna
    00
    dq
    00
    +O
    ?
    3
    ?
    1
    2
    d ln a
    00
    dq
    00
    ?
    2
    !
    :
    Equivalent instantaneous area distortion coe?cients d ln a
    kl
    =dq
    00
    are expressed as:
    1
    2
    d ln a
    ij
    dq
    00
    ?
    d ln a?
    ij
    dq
    00
    ˘
    ?
    Eq: 3 : 7
    Cov
    00
    2
    X
    X
    a
    X
    i;j
    1
    2
    d ln a
    00
    dq
    00
    ?
    d ln a?
    00
    dq
    00
    ˘2
    X
    X
    a
    X
    0 ; 0
    +O
    ?
    ?
    1
    2
    d ln a
    00
    dq
    00
    ?
    2
    !
    y
    If indeed V
    Eq: 4 : 14
    = V
    Eq: 4 : 15
    , then Equations 4.14 and 4.15 can be used together with Eq. 4.4, the sum rule for a
    X
    i;j
    ,
    to recover Poisson statistics, essentially by rebinning. Coe?cients to the a
    X
    i;j
    terms cancel, leaving ? = V . However, we
    believe Equation 4.14 should scale with the recorded variance: V
    Eq: 4 : 14
    6= ?.
    z
    Augustin notes that high-quality ?ts to covariances for (i; j) 6= (0; 0) are achieved if the following relation is assumed:
    2V
    Eq: 16
    = ?+Cov
    ?
    Q
    0
    0 ; 0
    ; Q
    0
    0 ; 0
    ?

    where the authors have assumed Var = ? in the ?at image prior to expressing the change in Var in Eq. 3.8.
    It should be pointed out that all of the above sets of expressions each internally recover Poisson statistics as
    covariances out to large lags ij are summed: The equations of
    3
    ,
    4
    and
    5
    have terms that cancel by subtraction,
    while (our) equations
    3
    and
    7
    are derived from this same principle. Also, the mean{variance relation out to large
    signal levels ? may be calculated recursively if ? ln a
    ij
    (? + ?j?) (Eq.
    1
    ) is computable.
    APPENDIX B. ELECTROSTATIC DRIFT MODEL FOR COLD ELECTRONS
    Here we brie?y summarize the electrostatic drift ?eld calculation, which is described in greater detail else-
    where
    8
    ,
    11
    ,
    12
    and reproduced here for convenience. Figure
    8
    shows the assumed pixel geometry and electrostatic
    elements in the model, but here only depicts a 2?2 pixel region close to the channel. Collected conversions are
    represented by the four \bubble" like structures that hover over positions between pairs of extruded arrows in
    two dimensional symmetric arrangement about the potential wells. The potential wells (\bubbles") do not lie
    in the plane of the front side clock structure because these devices feature a buried channel. The integrating
    and barrier clocks are strip-like equipotentials that extend for long distances along the serial address (i) axis and
    provide boundary conditions that justify utilizing the method of images for a small channel depth z
    chan
    relative
    to other relevant dimensions (distance to positions within the drift region, pixel dimension, and the combined
    width of adjacent integrating clocks). The accumulated conversions constrained to the potential well will appear,
    in the far ?eld approximation, to have an equal and opposite image charge distribution on the opposite side of
    the clock plane, acting together as the perturbative, aggressor dipole ?eld denoted p~
    ij
    .
    Similar arguments are used to describe the far-?eld in?uence of the cannel stop ion implants, which, under
    depleted operation, act as another dipole moment in the z direction with translational invariance along the
    parallel transfer (j) axis. These are denoted ˘~
    cs
    in Figure
    8
    . Finally, adjacent integrating and barrier clocks act
    as dipole moments con?ned to the i-j plane with translational invariance along the serial transfer (i) axis, shown
    as ˘~
    ck
    .
    The predominant component of the drift ?eld is the backdrop ?eld, denoted E~
    BD
    (z), is a one dimensional
    solution of Poisson's equation in the depleted silicon. The in?uence of the periodic and non-periodic contributors
    described above can be added in superposition because they explicitly satisfy Gauss' law everywhere except for in
    those small volumes that contain ?nite bound charge densities ˆ
    b
    not described in a the one-dimensional impurity
    concentration pro?le N(z), and surface charge densities ˙
    f
    arising on semiconductor-conductor interfaces with
    nonzero normal component of the local electric ?eld. The constant of integration for this backdrop ?eld E~
    BD
    (z) is
    chosen such that a zero backside bias BSS implies a zero electric ?eld strength directly inside the backside surface
    of the sensor. In Figure
    8
    then, the backside window is located a large distance directly above the con?guration
    of electrostatic moments shown, i.e. toward where the nine vertical arrows point. Cold carrier pixel boundaries
    undergo shifts in response to changes in the positions and magnitudes of the and charge con?guration moments
    p~
    ij
    , ˘~
    cs;i
    and ˘~
    ck;j
    .
    The equations used for vector integration along the drift ?eld lines are:
    E~
    tot
    (~x) = E~
    BD
    (z) + ?E~ (~x)
    (12)
    E~
    BD
    (z) =
    ?
    1
    ?
    0
    ?
    Si
    Z
    t
    Si
    z
    dzN
    a
    (z)? V
    BSS
    =t
    Si
    ?
    z^
    (13)
    d~l =
    ?
    ?
    ?
    EE~~ ((~~xx))
    ?
    ?
    ?
    ds
    (14)
    ~x
    i +1
    = ~x
    i
    +d~l
    (15)
    t
    coll
    (x~
    0
    ) =
    Z
    ~x
    0
    ~x ? k
    ^=
    z
    ch
    dl
    ?
    e
    (E(z);T) jE(z)j
    :
    (16)
    where the cold electron collection time t
    coll
    (x~
    0
    ) is used to estimate the thermal di?usion at the end of the
    axial drift, using the mobility in the small ?eld limit, ˙
    2
    ˇ 2k
    B
    T=q
    e
    ? ?
    e
    (E = 0;T) ? t
    coll
    (x~
    0
    ). Mobility

    Figure 8. The assumed ?eld geometry in the electrostatic model.
    parameterizations of Jacoboni
    14
    were used, although we believe there is mounting evidence in X-ray illumination
    data for these sensors to suggest that the velocity saturation e?ect for carriers is not as strong as that model
    provides, for the operating conditions in question.
    Details of the electrostatic in?uence by periodic barriers, denoted ˘~
    cs
    and ˘~
    ck
    in Figure
    8
    , are contained in
    the term ?E~ (~x), and are subdominant for positions far from the channel, ~x ? z^ ˛ z
    chan
    , but compete with and
    can ultimately dominate in?uence of the backdrop ?eld E~
    BD
    near the channel.
    x
    The image charge modeling
    strategy used, and also the (in?nite) periodic arrangement of the channel stop and clock barrier potentials are
    explicitly given in Rasmussen
    8
    [xx3.2-3.3] and are not reproduced here.
    APPENDIX C. SHOELACE FORMULAE UTILIZED
    After pixel boundaries are sampled via the drift calculation, they are compiled into lists that comprise polygonal
    representations of the pixels. The following formulae were used to compute geometric parameters for each pixel.
    While direct mapping (e.g., ~x 2 pixel[i; j] vs. ~x 3 pixel[i; j]) is utilized for certain simulation applications via
    e?cient point-in-polygon routines, interpretation of recorded images may be aided with use of ancillary pixel
    information according to the recorded signal distribution in the pixels. Polygon representations of pixels in?u-
    enced by the recorded signal distributions are straightforward to perform if the signal scale speci?c perturbation
    patterns are known. The following shoelace formulae were given in Rasmussen
    8
    [x3.5] that connect pixel area A
    ij
    ,
    pixel astrometric shifts I x
    ij
    & I y
    ij
    , and second moments I xx
    ij
    , I yy
    ij
    & I xy
    ij
    , given an ordered set of N vertices
    (x; y)
    k
    where (x; y)
    N
    ? (x; y)
    0
    . For the pixel boundary calculations represented in this work, we typically worked
    with either 15 or 25 points per side (60 ? N ? 100):
    x
    Any charge con?gurations near the front side potential wells
    9
    that would motivate carrying higher order terms in a
    multipole expansion, are not entertained here. We imagine that such terms would include any ?nite spatial extent in depth
    and width of the channel stop implant, and any spatial extent in 3 dimensions of the accumulated signal carriers collected
    in the potential well. Such higher order terms necessarily would have a shorter range of in?uence (j?Ej ˘ r
    ? s
    ;s ? 4).
    At a level where they might be important in the drift calculations, these terms will also in?uence the shapes of charge
    clouds residing in adjacent wells, complicating the drift calculation. We plan to neglect such terms until there is su?cient
    evidence in the data to suggest their importance.

    A
    ij
    ? +
    1
    2
    N
    X
    ? 1
    k =0
    (x
    k +1
    y
    k
    ? x
    k
    y
    k +1
    );
    (17)
    I xx
    ij
    A
    ij
    ? ?
    1
    12
    N
    X
    ? 1
    k =0
    (y
    k +1
    ? y
    k
    )(x
    2
    k
    +x
    2
    k +1
    )(x
    k
    + x
    k +1
    );
    (18)
    I yy
    ij
    A
    ij
    ? +
    1
    12
    N
    X
    ? 1
    k =0
    (x
    k +1
    ? x
    k
    )(y
    2
    k
    +y
    2
    k +1
    )(y
    k
    + y
    k +1
    );
    (19)
    I xy
    ij
    A
    ij
    ? +
    1
    6
    N
    X
    ? 1
    k =0
    (x
    k +1
    ? x
    k
    )x
    k
    (y
    2
    k
    +y
    2
    k +1
    + y
    k
    y
    k +1
    )
    +
    1
    24
    N
    X
    ? 1
    k =0
    (x
    k +1
    ? x
    k
    )
    2
    (y
    2
    k
    + 3y
    2
    k +1
    + 2y
    k
    y
    k +1
    );
    (20)
    I x
    ij
    A
    ij
    ? ?
    1
    6
    N
    X
    ? 1
    k =0
    (y
    k +1
    ? y
    k
    )(x
    2
    k
    +x
    2
    k +1
    + x
    k
    x
    k +1
    );
    (21)
    I y
    ij
    A
    ij
    ? +
    1
    6
    N
    X
    ? 1
    k =0
    (x
    k +1
    ? x
    k
    )(y
    2
    k
    +y
    2
    k +1
    + y
    k
    y
    k +1
    ):
    (22)
    The sign of Eq.
    17
    corresponds to a speci?c choice of chirality for the polygonal vertex list. The quanti-
    ties above are used to evaluate distortions to pixel area (? ln A
    ij
    ), pixel astrometric shift vectors (e.g., p~
    ij
    ? x^ =
    [I x
    ij
    A
    ij
    ]=A
    ij
    ) and pixel ellipticities (??
    1 ;ij
    = [I xx
    ij
    A
    ij
    ? I yy
    ij
    A
    ij
    ]=[I xx
    ij
    A
    ij
    +I yy
    ij
    A
    ij
    ]; ??
    2 ;ij
    = 2 I xy
    ij
    A
    ij
    =[I xx
    ij
    A
    ij
    +
    I yy
    ij
    A
    ij
    ]). It may be possible for existing pixel data pipelines to be retro?tted to take advantage of such book-
    keeping information when estimating object parameters, particularly for PSF estimation purposes.
    APPENDIX D. A LINEAR PERTURBATION TEMPLATE TO REPRESENT
    DYNAMIC PIXEL RESPONSE
    The proportional pixel boundary shifts laid out by Antilogus et al.
    3
    x4.2 (and subsequently Refs.
    4
    ,
    5
    ), uses
    constructions that linearly accumulate the in?uence of aggressors in the pixel's vicinity. The coupling coe?cients
    are determined using a matrix inversion of constraint equations (containing measured covariances) that utilize
    re?ection symmetries, and a sum rule (e.g., Ref.
    4
    x6.1). Our detailed electrostatic drift calculation may also be
    applied, and we can do so while explicitly guaranteeing the continuity equation and one-to-one mapping between
    a two dimensional continuous position ?eld and pixel address. In other words, the Greens function doesn't su?er
    problems intrinsic to a general arrangement of rectangular pixels that naturally over- and under-claim pixel
    \ownership" of the continuous position ?eld.
    In the same spirit, we apply the Greens function according to the supposition that all de?ections of pixel
    boundaries are perturbations that scale linearly with aggressor amplitude. We refer to application of the linear
    perturbation equations collectively as the BF template. Figure
    9
    illustrates the geometry. The equations used are
    as follows, where ?~c
    t
    k;k +1
    are computed distortion vectors of two adjacent corners for the template aggressor p
    t
    ,
    nominally separated by a single pixel step along the positive m transfer direction e^
    m
    : m 2 f0; 1g, e^
    m
    2 f
    ?
    1
    0
    ?
    ;
    ?
    0
    1
    ?
    g:
    ?~x
    t
    l
    ?e^
    m
    ?
    ?
    ?~c
    t
    k
    +
    ?
    l
    n? 1
    ?
    ?
    s e^
    m
    + ?~c
    t
    k +1
    ? ?~c
    t
    k
    ?
    ?
    ? e^
    m
    ?~x
    t
    l
    ?e^
    ( m +1)mod2
    ? ?~c
    t
    l
    ? e^
    ( m +1) mod 2
    + (?d
    t
    l
    ? ?d
    t
    0
    );
    where ?~x
    t
    l
    , l 2 f0:::n ? 1g are solutions to the electrostatic drift calculation (for template aggressor p
    t
    ) that
    form a locus for the pixel boundary of this border, ?d
    l
    are the boundary distortions perpendicular to the border
    axis in the (m + 1) mod 2 direction as shown. As usual, i and j are the lag indices in the m = 0 and m = 1

    directions respectively, and s is the pixel spacing. The per-lag template quantities ?~c
    t ( i;j )
    k
    and ?d
    t ( ij )
    l
    are then
    compiled by summing over in?uences of aggressors accumulated at the channel according to:
    ?~c
    ( q;r )
    k
    =
    X
    ij
    ?
    p
    q ? i;r ? j
    p
    t
    ?
    ?~c
    t ( i;j )
    k
    ?d
    ( q;r )
    l
    =
    X
    ij
    ?
    p
    q ? i;r ? j
    p
    t
    ?
    ?d
    t ( i;j )
    l
    ~x
    ( q;r )
    l
    ?e^
    m
    ? ((mmod2==0)?q:r)s+
    ?
    ?~c
    ( q;r )
    k
    +
    ?
    l
    n? 1
    ?
    ?
    s e^
    m
    + ?~c
    ( q;r )
    k +1
    ? ?~c
    ( q;r )
    k
    ?
    ?
    ? e^
    m
    ~x
    ( q;r )
    l
    ?e^
    ( m +1)mod2
    ? (((m+1) mod 2 == 0)?q : r) s+?~c
    ( q;r )
    l
    ? e^
    ( m +1) mod 2
    + (?d
    ( q;r )
    l
    ? ?d
    ( q;r )
    0
    ):
    In the above, ?~c
    ( q;r )
    k
    and ?d
    ( q;r )
    l
    are the resulting total perturbations to the pixel corners and boundaries from
    multiple aggressors, respectively; and the boundary pairs ~x
    ( q;r )
    l
    for this border are combined the boundary pairs
    for each of the three other borders to produce a (closed) polygonal description of pixel (q;r). Point-in-pixel
    algorithms, as well as the shoelace formulae of Appendix
    C
    are then readily applied to these distorted pixel
    descriptions.
    The preceding provides a generalized, 2-dimensional application of drift calculation results as linear pertur-
    bations. It is analogous to the simpler, perturbative pixel border shifts described in previous work.
    3
    {
    5
    ACKNOWLEDGMENTS
    This material is based upon work supported in part by the National Science Foundation through Cooperative
    Support Agreement (CSA) Award No. AST-1227061 under Governing Cooperative Agreement 1258333 managed
    by the Association of Universities for Research in Astronomy (AURA), and the Department of Energy under
    Contract No. DEAC02-76SF00515 with the SLAC National Accelerator Laboratory. Additional LSST funding
    comes from private donations, grants to universities, and in-kind support from LSSTC Institutional Members.
    REFERENCES
    [1] Kurita, N., Kahn, S., Stubbs, C. W., Ritz, S., Nordby, M. E., and Riot, V. J., \Large synoptic survey tele-
    scope camera design and construction," in [Advances in Optical and Mechanical Technologies for Telescopes
    and Instrumentation], Navarro, R. and Burge, J. H., eds., Proc. SPIE 9912 (2016). in press.
    [2] Kahn, S., \Final design of the large synoptic survey telescope," in [Ground-Based and Airborne Telescopes
    IV], Hall, H. J., Gilmozzi, R., and Marshall, H. K., eds., Proc. SPIE 9906 (2016). in press.
    [3] Antilogus, P., Astier, P., Doherty, P., Guyonnet, A., and Regnault, N., \The brighter-fatter e?ect and pixel
    correlations in CCD sensors," Journal of Instrumentation 9 (03), C03048 (2014).
    [4] Guyonnet, A., Astier, P., Antilogus, P., Regnault, N., and Doherty, P., \Evidence for self-interaction of
    charge distribution in charge-coupled devices," A&A 575 , A41 (Mar. 2015).
    [5] Gruen, D., Bernstein, G. M., Jarvis, M., Rowe, B., Vikram, V., Plazas, A. A., and Seitz, S., \Character-
    ization and correction of charge-induced pixel shifts in DECam," Journal of Instrumentation 10 , C05032
    (May 2015).
    [6] Niemi, S.-M., Cropper, M., Szafraniec, M., and Kitching, T., \Measuring a charge-coupled device point
    spread function. Euclid visible instrument CCD273-84 PSF performance," Experimental Astronomy 39 ,
    207{231 (June 2015).
    [7] Miyazaki, S. Private communication, discussed at Foreground Physical E?ects on LSST Weak Lensing
    Science: A Workshop on the Impact of the Last Kiloparsec, University of California at Davis, regarding
    measurements made using HSC sensors. (2015).
    [8] Rasmussen, A., \Use of sensor characterization data to tune electrostatic model parameters for LSST sen-
    sors," Journal of Instrumentation 10 , C05028 (May 2015).

    8×10
    −4
    10
    −3
    1.2×10
    −3
    1.4×10
    −3
    1.6×10
    −3
    1.8×10
    −3
    2×10
    −3
    2.2×10
    −3
    9.2×10
    −4
    9.4×10
    −4
    9.6×10
    −4
    9.8×10
    −4
    10
    −3
    parallel transfer coordinate [cm]
    serial transfer coordinate [cm]
    template_pixel_contour_map_pd_+2.8189_detail.qdp
    ij=10
    ij=00
    ij=01
    ij=11
    ij=20
    ij=21
    Figure 9. An illustration to show the element level information content of the pixel distortion template. This is for the
    lower boundary for the pixel with lag ij = 11 (and contains the same information as the upper boundary for the pixel
    with lag ij = 10). The (red) ?lled triangles connected by solid lines form the pixel boundary locus ~x
    t ( i;j )
    l
    and divide the
    pixel areas labeled by their corresponding lag (ij = 00; 10; 20; 01; 11; 21). The (black) solid line connecting ?lled circles
    shows the position of the undistorted pixel boundary (between ij = 10; 11). Open (blue) triangles show positions of the
    distorted pixel corners, with (black) solid lines (?~c
    t (1 ; 1)
    k;k +1
    ) connecting them to the undistorted corners (?lled circles). With
    distorted corner coordinates projected onto the e^
    0
    (serial, horizontal) coordinate and the separation divided equally (n
    samples per side), the e^
    1
    (parallel, vertical) de?ections ?d
    t (1 ; 1)
    l
    of the border are recorded and stored as a template for
    downstream use. With n = 25 as shown here, the template can be stored with modest memory requirements, 32 numbers
    per border per lag, and when compiled for an arbitrary recorded charge distribution, 96 (x; y) pairs per pixel that form
    a closed polygon. Application of the shoelace formulae (Appendix
    C
    ) provides e?cient distillation of this information to
    6 leading geometric terms per pixel. This template was generated for an aggressor amplitude p? = 2:819 p
    0
    : 492 times the
    e?ective exposure averaged, aggressor level ?t described in Table
    1
    . This corresponds to a readout time aggressor signal
    of 119ke
    ?
    , comparable to the full well depth for these sensors.
    [9] Lage, C., \Measurements and simulation of the brighter-fatter e?ect," Talk delivered at Foreground Physical
    E?ects on LSST Weak Lensing Science: A Workshop on the Impact of the Last Kiloparsec, University of
    California at Davis (Dec. 2015).
    [10] O'Connor, P., Antilogus, P., Doherty, P., Haupt, J., Herrmann, S., Hu?er, M., Juramy-Giles, C., Kuczewski,
    J., Russo, S., Stubbs, C., and Berg, R. V., \Integrated system tests of the LSST raft tower modules," in
    [High Energy, Optical, and Infrared Detectors for Astronomy VII], Holland, A. D. and Beletic, J., eds.,
    Proc. SPIE 9915 (2016). in press.
    [11] Rasmussen, A., \Pixel area variations in sensors: a novel framework for predicting pixel ?delity and distor-
    tion in ?at ?eld response," Journal of Instrumentation 9 , C4027 (Apr. 2014).
    [12] Rasmussen, A., Antilogus, P., Astier, P., Claver, C., Doherty, P., Dubois-Felsmann, G., Gilmore, K.,
    Kahn, S., Kotov, I., Lupton, R., O'Connor, P., Nomerotski, A., Ritz, S., and Stubbs, C., \A framework
    for modeling the detailed optical response of thick, multiple segment, large format sensors for precision
    astronomy applications," in [Modeling, Systems Engineering, and Project Management for Astronomy VI ],
    Proc. SPIE 9150 , 915017 (Aug. 2014).
    [13] Lupton, R. Private communication, regarding HSC fork of LSST DM stack (2016).
    [14] Jacoboni, C., Canali, C., Ottaviani, G., and Alberigi Quaranta, A., \A review of some charge transport
    properties of silicon," Solid State Electronics 20 , 77{89 (Feb. 1977).

    Back to top