Model-independent Characterization of Charge Diffusion in Thick Fully Depleted CCDs

The charge-coupled device (CCD) is one of the most sensitive visible-light detectors available for use in modern astronomical instruments. Because of the increasing cosmological interest in extragalactic sources, considerable effort is devoted to pushing CCD sensitivity further into the near-infrared range. Thick fully depleted silicon CCDs (FDCCDs), originally developed for use in X-ray astrophysics, have been proposed for use in current and upcoming wide-field surveys, e.g., Panoramic Survey Telescope and Rapid Response System (PanSTARRS), Dark Energy Survey, HyperSuprime, and Large Synoptic Survey Telescope (LSST). These FDCCDs offer improved infrared sensitivity, but they experience increased charge diffusion as a result of their thickness. An accurate quantification of charge diffusion is therefore essential when characterizing such devices.

Several techniques for investigating charge diffusion in CCDs have been developed. Optical methods include a virtual knife-edge scan (Karcher et al. 2004) that measures the profile of a finely focused spot, measurement of the line-spread profile using a projected slit image angled with respect to the CCD columns (Robbins 2011), and modulation transfer function measurement using interferometrically projected fringe patterns of variable spatial frequency (Takacs et al. 2010). Other techniques use ionizing radiation, either in the form of X-ray photons from a radioisotope or accelerator sources (Janesick 2001; Kraft et al. 1995) or charged particles such as cosmic-ray muons (Robbins 2011). When a mask can be placed directly on the CCD's sensitive surface, the mesh or virtual pixel scan technique (Tsunemi et al. 1997) gives a detailed and accurate characterization of the point-spread function (PSF).

Rodney and Tonry (2006) introduced a methodology that is very similar to the one used here, relating the histogram of pixel values found in an image exposed to ⁵⁵Fe events to the charge diffusion, and providing a simple statistic that can be used to infer charge diffusion given an average event. This article offers an alternate approach that has two advantages: we need not assume a specific model for how charge carriers diffuse through the silicon, and we can construct two different estimators that bound the true amount of charge diffusion.

In this analysis we work with images resulting from exposure of the CCD to ⁵⁵Fe X-rays. The absorption of a single soft X-ray photon generates charge carriers in a tight cluster (Tsunemi et al. 1999) that diffuse laterally to produce a charge cloud with a Gaussian profile (Janesick 2001; Hopkinson 1987; Groom et al. 2000). The width σ of the cloud grows as carriers drift from the conversion point to the collecting electrodes. A single ⁵⁵Fe image typically contains thousands of these charge clouds.

With a uniform flux of X-rays incident on the CCD's entrance window, photon conversions are distributed uniformly across the pixel array and with an exponential falloff in depth. (The attenuation length of Mn Kα photons is 29 μm in silicon.) The geometry is illustrated in Figure 1. In a back-illuminated device, photon conversions close to the entrance window result in the widest charge clouds, since the charge carriers they produce have the longest transit times to the potential wells on the front side. The profile of conversion depths therefore leads to a cluster width distribution with a sharp cutoff at σ_max, corresponding to photons that convert just at the entrance window. This limiting width is significant because it predicts the amount of diffusion experienced by visible-light photons, which convert within the first few microns of silicon.

Zoom In Zoom Out Reset image size

Two difficulties hamper the analysis of X-ray images to extract σ_max. First, the profiles of X-ray clusters are severely undersampled, given the pixel size and typical operating conditions. Second, it is necessary to examine the ensemble of cluster sizes and develop a robust estimator for only that subset having the largest amount of diffusion. For device characterization purposes, it is most desirable to obtain a strict upper bound on σ_max.

In the following sections we report a new technique for extracting the maximum charge cloud width from an X-ray image. The method is computationally inexpensive, accurate for a wide range of PSF sizes, and avoids the need for any diffusion model. We begin in § 2 with a simple algorithm to select Mn Kα photon conversions within an image. In § 3, we characterize each event with a single number—the "difference ratio." Section 4 considers the aggregate of all difference ratios and proposes two different methods for determining the maximum PSF width. We test the algorithm using simulated data in § 5 and apply it to images from a high-resistivity device in § 6. Conclusions are presented in § 7.

Starting with a background-subtracted image file of an FDCCD exposed to ⁵⁵Fe radiation (Fig. 2), we reconstruct events and select those originating from conversions of 5.988 keV Mn Kα photons, which produce 1590 electron-hole pairs in silicon at 173 K (Lowe and Sareen 2007). In our target FDCCDs the diffusion scale is on the order of the pixel size, so it is sufficient to consider events whose charge is contained in 2 × 2 pixel clusters (in the case of FDCCDs with larger PSFs, a straightforward Gaussian fit to the charge distribution would suffice and this technique would not be necessary). The algorithm proceeds as follows:

Zoom In Zoom Out Reset image size

• 1.  
  Calculate the total charge in every 2 × 2 pixel "box" within the image.
• 2.  
  For every pair of overlapping boxes, immediately reject the one with the smaller total charge, leaving only the boxes that are local maxima.
• 3.  
  Histogram the totals of the remaining boxes (Fig. 3).
• 4.  
  Identify the window centered on the Mn Kα peak position with upper bound at the location of the minimum between the Kα and Kβ peaks, and remove from consideration all boxes with totals outside this window.

Zoom In Zoom Out Reset image size

This method of selection is unbiased in that no event will be represented by multiple boxes: an event may lie within multiple boxes, but only one of those can be a local maximum. The algorithm also tolerates overlapping events, since a box containing multiple events has a total outside the selection window and is ignored altogether.

Suppose that a single X-ray photon is absorbed at the point (x₀,y₀) in an infinite unpixelated plane. Since the diffusion PSF is Gaussian, the flux S is given by

where s is the total signal produced and σ is the width of the charge cloud. Let u and v be the total charge gathered in the positive- x and negative- x half-planes, respectively. Integrating equation (1), we find that

where erf is the error function.¹

For any X-ray event, we define the difference ratio d₁ such that

Note that d₁ is closely related to the first moment of the charge distribution about the coordinate origin.

Alternatively, we may express the difference ratio in terms of x₀ and σ:

which follows from equations (2) and (3). Note that the difference ratio does not depend on s.

In the pixelated case, consider a 2 × 2 box that has been selected by our algorithm in § 2. As before, suppose that B contains an X-ray event located at coordinates (x₀,y₀) relative to its center, and let the PSF width be σ. Using equation (4), we can calculate the ideal difference ratio for this box as d₂(x₀,σ).

If the X-ray charge is fully contained in B, then u and v are exactly equal to the summed fluxes b + d and a + c, respectively. When the PSF is small, this remains a valid approximation, and so u ≈ b + d and v ≈ a + c. Equation (3) can then be applied to this event as well, and we find that d₂(x₀,σ) ≈ d₁(b + d,a + c) for every box B. We now have a practical method to compute the difference ratio of any event without knowledge of x₀ or σ. (We have only considered the difference ratio in the x -axis direction. Clearly, an analogous quantity can be defined in the y -axis direction. Since charge diffusion is isotropic, we do not differentiate between the axes in our analysis, although considering the axes separately can be useful for quantifying anisotropy in nonideal CCDs.)

The primary utility of the difference ratio is that it provides a simple relationship between measured pixel values (a, b, c, and d) and purely theoretical Gaussian profile parameters (x₀ and σ). The value of the difference ratio d for a single cluster does not provide enough information to unambiguously determine the cluster's width σ or centroid location x₀. However, knowing that x₀ is uniformly distributed in and that large σs are proportionately overrepresented in the population, we can estimate σ_max from the distribution of d for the entire ensemble, as shown in the next section.

Once we have computed the x - and y -directed difference ratios for every 2 × 2 box, we plot them in a histogram with area normalized to unity (Fig. 4). At this point, data from multiple images may be combined to reduce statistical errors. The histogram is seen to have a pronounced peak close to the maximum difference ratio that levels off to a plateau as the difference ratio decreases. Two features of this histogram can be used to estimate σ_max —the location of the peak and the height of the plateau.

Zoom In Zoom Out Reset image size

We can determine σ_max from the location of the histogram's peak, which must come from those events with the widest PSF. The variable |x₀| has a uniform distribution between 0 and , but since the error function is approximately linear for small arguments and asymptotically constant for large arguments, a greater density of difference ratios will be obtained when |x₀| is as large as possible. Since σ has a sharp cutoff at σ_max and |x₀| has a sharp cutoff at , a pronounced peak is expected at the region of greatest density. Therefore, the histogram peak is located at . Solving for σ_max, we obtain

where p is the location of the peak and erfinv is the inverse error function.

A second indicator of σ_max is the value of the plateau. For difference ratios near zero, the error function is nearly linear, and . Inverting, we have . Then the number of points with a difference ratio between d_a and d_b is equal to the number of points with |x₀| between and , where σ takes some value less than σ_max. Since |x₀| is uniformly distributed from 0 to , this is simply equal to , where n is the total number of points. Finally, a normalized histogram must have constant value z throughout the region where the error function is linear, and we have . Since the distribution of σ is biased toward σ_max, we have

Equations (5) and (6) represent two estimators for σ_max based solely on the shape of the difference-ratio histogram, which we will hereinafter refer to as σ_peak and σ_plat. These estimators make only a few assumptions about the data: that charge diffusion results in Gaussian PSFs that are largely contained within 2 × 2 pixel clusters, that X-ray events are randomly distributed across the surface of the CCD, and that the largest PSF widths are the most common.

To estimate the location of the peak, we use Parzen windows to construct a continuous estimate of the difference-ratio probability density function (Parzen 1962). The weighting function is a Gaussian distribution with standard deviation 0.02, since this is found to produce the best estimate for sets of 5000 events (about 10 images). The peak location is then defined as the mode of the estimated probability density function.

To estimate the height of the plateau, we again use Parzen windows, but some adaptations are made. In order to avoid artificial falloff toward the negative side of the distribution, we construct our estimator on a modified data set comprised of all difference ratios together with the additive inverses of all difference ratios. The weighting function is a Gaussian distribution with standard deviation 0.1 (it can be wider because this part of the histogram is smoother). The plateau height is then given by evaluating the probability density function at 0.

In the following sections we examine the performance of these estimators and the validity of the simplifying assumptions used.

The difference-ratio analysis was applied to simulated images with predetermined σ_max. Values of σ were simulated using a Monte Carlo procedure, which takes into account the absorption probability for 6 keV photons in silicon, the electric field in the CCD, and the transport and diffusion properties of electrons in silicon (Eremin and Li 1995). The simulated distribution of σ values (Fig. 5) was scaled linearly in order to simulate varying amounts of charge diffusion. Values of x₀ and y₀ were chosen uniformly at random between 0 and .

Zoom In Zoom Out Reset image size

Sets of 10 images, each containing 500 X-ray events, were simulated with and without noise for three values of σ_max. The standard deviation of the noise was taken to be of the signal produced by a Kα event (the experimentally determined noise for the particular FDCCD analyzed in § 6). One hundred distinct data sets were created for each set of parameters and analyzed using several algorithms. We determined σ_peak and σ_plat using the difference-ratio algorithms outlined in § 4. The average and standard deviation across the 100 trials for each test are given in Table 1 and Figure 6.

Zoom In Zoom Out Reset image size

The error in σ_peak is caused by two factors: PSF truncation and incorrect determination of peak location. PSF truncation occurs when σ_max is large and charge outside the 2 × 2 analysis box is excluded from the difference-ratio calculation. In the worst case, the difference ratio is calculated as

which is always less than the actual difference ratio of

Second, the histogram peak is sloped asymmetrically, so use of a Gaussian kernel density estimator shifts its apparent location—an effect that is compounded by the addition of Gaussian noise. Due to the behavior of the inverse error function, the algorithm is far more sensitive to errors in peak location for small σ_max, but this error exists for large σ_max as well. The histogram peak is sloped more steeply on the right, so peak location errors tend to shift it left. Therefore, both sources of error increase σ_peak, and σ_peak is always an upper bound on σ_max.

Error in σ_plat stems from the assumption that σ ≈ σ_max. Most values of σ are close to σ_max, but all are at least slightly smaller. When we derive a value for σ_max, we are actually deriving a value for some smaller average σ for the entire spectrum of events. Therefore, we underestimate σ_max whenever we use this method. The exact amount of this underestimate depends on the distribution of σ, but is about 15% for our simulated distribution.

The rms variation between simulation runs is well below the absolute error, indicating that the upper and lower bounds are reproducible. (The lower bound is especially tolerant of noise, because noise does not distort the flat section of the histogram.) Using σ_peak and σ_plat, we obtain a tight upper bound and a robust lower bound on σ_max without assuming any specific distribution of σ.

The methods described in this report have been used to analyze a next-generation FDCCD being developed for the LSST. The FDCCDs under development must balance two conflicting design goals: high quantum efficiency in the near-infrared and small PSF size (Radeka et al. 2009). Since the technology required is still relatively untested, the development of an accurate PSF measurement technique is paramount.

Prototype FDCCDs have been obtained and tested (O'Connor et al. 2008). Device 106-06 is a 100 μm thick backside-illuminated FDCCD with 13.5 μm pixels. The electron transport model given by Eremin and Li (1995), which takes velocity saturation and temperature dependence into account, was used to predict PSF size. Separately, sets of 50 fully depleted X-ray images were acquired at bias voltages ranging from -30 V to -50 V in 2 V increments. These images were analyzed using the Q -ratio method to give a PSF estimate and analyzed using the difference-ratio technique to give upper and lower bounds (Table 2 and Fig. 7).

Zoom In Zoom Out Reset image size

In previous work (Rodney and Tonry 2006), an alternative measure was proposed for determining σ_max from ⁵⁵Fe data. Those authors considered the quantity

for all 3 × 3 pixel boxes identified as containing Kα hits. Thus, 0 < Q < 1, and Q = 0 for an event whose charge is completely contained in a single pixel. The aggregate Q ratio for an image is formed by superimposing all Kα hits. To derive σ_max from Q, they performed simulations using a simple two-parameter diffusion model, with (x,y,z) conversion locations drawn from the appropriate distributions. From multiple simulations the relation between σ_max and Q can be tabulated. Values of σ_max found in this way are denoted as σ_Q. This computation is also carried out using the diffusion model of Eremin and Li (Fig. 8), and the results are denoted by σ_Q^'.

Zoom In Zoom Out Reset image size

The Q -ratio method provides an easily computed relative measure of diffusion that can be used for comparing devices of different thickness, doping, and substrate bias, as shown in Rodney and Tonry (2006). Its absolute calibration scale, however, relies on a simplified diffusion model (Rodney and Tonry's eq. [3]), which assumes that the diffusion width depends on the conversion point's height z above the front-side gates as , with a correction applied for conversions within some height above the gates for which σ = 0. There are some conditions for which the functional form of σ(z) will depart from the assumed model. An example is shown in Figure 8. The graph compares the σ(z) dependence from Rodney and Tonry (2006) with a model that includes the electric field z dependence and the nonlinear velocity-field characteristic as

where D is the transverse diffusion coefficient, v is the electron velocity from Eremin and Li (1995), and E(z) = V/d - e·N_A(d - 2z)/2. The device modeled is just above full depletion, having thickness 100 μm, resistivity 3 kΩ·cm, barrier-field height 13.5 μm, and V_sub = -35 V. The surface diffusion σ_max is 5.4 μm, and both models are normalized to this value at z = 0. Note that the full model has smaller σ(z) near z = 0, which will result in a higher estimate of σ_max than the simplified model. In general, this discrepancy is largest for thick devices biased at low overdepletion voltages; for well-overdepleted sensors the agreement between the two models is much closer. The barrier-field focusing effect introduced near the gate has little effect for this device, as the absorption probability in this region is only 1.8%. Nonnormal X-ray incidence will increase the fraction of events converting at shallower depths in Si, which will increase the aggregate Q ratio, since those clusters have the widest diffusion.

The present model simply seeks to identify the maximum-size clusters, which can be recognized by the distribution of first moments as described previously, regardless of the details of the σ(z) and P_abs(z) functions. However, it provides only upper and lower bounds on σ_max. While these bounds are shown to be robust, they leave a range of uncertainty that can be large for low-S/N images or for undersampled or truncated PSFs. A combination of difference ratio, Q ratio, and optical techniques may be advantageous when the most accurate results are required.

We have developed an algorithm to quickly and accurately estimate the contribution of charge diffusion to the PSF of thick fully depleted CCDs. The method uses a statistical analysis of the shapes of charge clusters in soft X-ray images to determine upper and lower bounds on the maximum PSF width, thereby bracketing the PSF expected for visible light. Simulations have been used to estimate the accuracy of the method over a range of PSF widths and with varying amounts of noise. All stages of the algorithm are computationally efficient—the analysis of a 100 Mpixel image takes less than 1 minute on a single-core desktop computer. Finally, the algorithm does not rely on a physical model of carrier transport and can be successfully applied to CCDs with widely varying charge diffusion characteristics. An analysis of an LSST prototype device shows good consistency with model predictions and with alternate methods of analysis. These properties make the difference ratio a useful tool for studying charge diffusion within CCDs.

Thanks go to Veljko Radeka, the Instrumentation Division Chair, and the Brookhaven National Laboratory Office of Educational Programs for allowing this research to proceed. The Large Synoptic Survey Telescope (LSST) is a public-private partnership. Funding for design and development activity comes from the National Science Foundation, private donations, grants to universities, and in-kind support at Department of Energy laboratories and other LSST Corporation Institutional Members. Support of the W. M. Keck Foundation for sensor development is gratefully acknowledged. This work is supported by in part the National Science Foundation under Scientific Program Order 9 (AST-0551161) and Scientific Program Order 1 (AST-0244680) through Cooperative Agreement AST-0132798. Portions of this work are supported by the Department of Energy under contract DE-AC02-76SF00515 with the Stanford Linear Accelerator Center and contract DE-AC52-07NA27344 with Lawrence Livermore National Laboratory.