Poisson_CCD: A dedicated simulator for modeling CCDs Free

Charge Coupled Devices (CCDs) have been the workhorse devices for astronomical imaging for some time. George Smith’s Nobel lecture in Ref. 1 gives an excellent summary of the early history. While other detectors are making inroads, CCDs are still the dominant imaging device in astronomical applications. In recent years, thick, fully depleted CCDs with their wide spectral response have been applied to spectroscopic applications and imaging. Although these devices have high quantum efficiency, relatively good linearity, and acceptable dynamic range, they have a number of problematic effects that can impact the precision and accuracy of astronomical data. It is important that these effects are well understood so that they can be accounted for during image processing. To help understand these effects, we have built a dedicated simulator, Poisson_CCD, which solves the electrostatics in the bulk silicon of the CCD and propagates incoming charges down to the collecting wells where they are collected and stored. The simulator has proven very useful for understanding a number of CCD effects, which will be described in this article.

Of course, detailed semiconductor modeling codes already exist, are commercially available, and have been validated against silicon results. What is the purpose of developing yet another simulator? The answer is severalfold. First, commercial semiconductor codes typically use proprietary source codes and are quite expensive, while the code described here is open source and freely available.² Second, using a commercial semiconductor device simulator requires spending quite a bit of time learning to use the code and set up the initial conditions. The code described here sets up the initial conditions for a typical CCD with a few simple configuration parameters. A third reason is that some of the commercial tools have trouble handling grids physically large enough to simulate a meaningful number of pixels. This code circumvents this by using a variable grid spacing, as described in Sec. II A. Also, it is hoped that this code is simple enough that it can be mastered by people who are not semiconductor experts. The target user group is people in the astronomy field who want to answer questions about CCDs without investing a great deal of time learning the details of semiconductor physics.

The code was developed as part of the development effort of the LSST. The LSST, originally known as the Large Synoptic Survey Telescope, is now known as the Simonyi Survey Telescope at the Vera Rubin Observatory and will be used to conduct a 10-year survey of the southern sky known as the Legacy Survey of Space and Time. This instrument is an innovative, large, fast survey facility currently under construction at Cerro Pachon in Chile.³ The digital camera, also currently under construction, consists of approximately 3.2 gigapixels and is the largest digital camera ever constructed. The camera uses fully depleted silicon CCDs that are back illuminated and 100 $μ$m thick in order to optimize quantum efficiency in the near infrared. The imaging area consists of 189 CCDs, with each CCD containing 16 imaging regions laid out in an $8×2$ array. Each imaging region has a pixel array with approximately $500×2000$ 10 $μ$m square pixels, giving 16 megapixels total. Each imaging region also has its own independent amplifier.^4,5 The focal plane contains CCDs from two different vendors, the ITL STA3800C from the University of Arizona Imaging Technology Laboratory⁶ and the E2V CCD250 from Teledyne E2V.⁷ However, although this code was developed and tested against these two CCDs from the LSST project, it has already found more general use on other CCDs,⁸ and the hope is that this will continue.

This article is divided into several sections. In Sec. II, we give an overview of the simulator, describing the basic structure of the simulation volume, how we solve the semiconductor equations and how we treat incoming photons. We also show a number of examples of the outputs available from the simulator. In Sec. III, we review a number of the validation tests that were performed to validate the simulation results against measured data of different kinds, and finally, we conclude.

The simulator performs two basic tasks, as shown in Fig. 1. First, given the charges in the silicon bulk and the boundary conditions determined by potentials applied to the silicon surface, the simulator solves Poisson’s equation numerically to determine the potentials and electric fields in the silicon bulk. The solution to Poisson’s equation is determined using the technique of successive over-relaxation (SOR—see, for example, Ref. 9) and using multi-grid methods (see Sec. II D) to speed convergence. In regions where there are mobile carriers (holes and electrons), quasi-Fermi level methods are used to simultaneously solve for the potentials and free carrier densities in the device. The electric fields are determined by numerically differentiating the electrostatic potential. In general, the simulator only solves for the potentials and free-carrier densities in equilibrium and is not intended to give transient solutions. However, it is possible to repeatedly solve the equations with slight changes in initial conditions in order to give transient results. This technique has been used to generate movies of the CCD charge transport, as discussed in Sec. III F.

After solving for the device potentials, the second major task of the simulator comes into play. As incoming photons enter the CCD, they generate hole–electron pairs. The electric field in the CCD separates these charges, and the electrons propagate down to the collecting wells where they are collected, stored, and later counted. The simulator models this carrier transport in a physically realistic way, in order to determine in which pixel a generated carrier ends up. This is very useful for modeling pixel distortions that result from electric fields in the device, either built-in electric fields, such as those due to “tree rings,”¹⁰ or electric fields due to collected charges, such as those that lead to the brighter-fatter (BF) effect.^11–15 Note that in an astronomical CCD, the time between incoming charges (on the order of ms) is typically much longer than the time required for a charge to propagate down to the bottom (on the order of ns). Also, a single charge has little impact on the existing potentials and fields. So it is an excellent approximation to assume that a single charge propagates in a frozen electric field. The simulator is designed so that if multiple charges are being added, the user can choose how often to re-solve Poisson’s equation. A comment here about the nature of the charges is in order. The simulator is designed for N-channel CCDs (with P-type substrate), where the collected charges are electrons, and for the rest of the article, we will make that assumption. There is no physical reason why it will not work for P-channel CCDs (with N-type substrate), but it is not currently designed to support them. Most of the work in adapting it for P-channel CCDs would be nomenclature related, basically interchanging the names of holes and electrons, but adjustment of parameters like mobilities and effective masses would also be needed. This could be done if there is a demand for it.

The simulator is written in C++ and is controlled by a text-based configuration file, which contains all of the information about the silicon volume, pixel sizes, number of pixels, any non-pixel regions, etc. The configuration file also defines the problem being solved. By convention, the configuration file has a .cfg extension, but this is not necessary.  Appendix A lists the configuration parameters. As the simulation progresses, it writes out a number of files. Large files containing information like potentials, charge densities, and electric fields at each grid point are written as high-density HDF5 files, having file extension .hdf5. Several smaller text files, with .dat extension, are written, which contain information on the grids or the number of electrons in each pixel. After the simulation has completed, easily modifiable Python scripts are used to plot out results as desired.

The simulation volume is set up on a fixed three-dimensional rectangular grid, which does not change once the simulation has started. One begins by deciding the number of grid cells in each dimension. Because multi-grid methods are used (see Sec. II D), the number of grid cells in each dimension must be a multiple of 32. Note that as a convention, we refer to the side of the CCD where the circuitry is patterned as the bottom, and the side where the incident light comes in as the top. For the LSST CCDs, which are 100 $μ$m thick and have pixels 10 $μ$m square, a typical resolution is to have 32 simulation grid cells per pixel, so that each grid cell is 0.31 $μ$m. Initially, the grid was defined to be completely uniform in all three dimensions. However, for thick CCDs like those used in the LSST, the potentials and fields change rapidly in the region near the bottom, and only very slowly near the top. When the simulation had enough resolution to be accurate in the rapidly changing region at the bottom, most grid cells were wasted near the top. Of course, adaptive grid methods solve this problem, but also make the code much more complex, which defeated the purpose of having a relatively simple simulator. The solution chosen in this work was to use a non-linear grid in the Z-dimension only. This has proven to give high resolution where needed, without adding significant complexity to the simulation code. Figure 2 shows the scheme. This introduces some additional partial derivatives, as discussed more in Sec. II C, but these values can be pre-calculated and this is much simpler than an adaptive grid scheme.

Because the simulator is not a general purpose semiconductor solver and is intended to model devices with a given structure, certain assumptions are made about the structure of the CCD which simplifies building the device structure. This is one of the reasons why this code is so much simpler than a commercial device simulator. It is assumed that the CCD is a slab of silicon with a given thickness given by the parameter “SensorThickness.” It is assumed that the top surface of the CCD is at a fixed voltage given by the parameter “Vbb.” The bottom surface of the CCD has voltages specified by the various gate potentials. The doping deep in the silicon is assumed constant with a value given by “BackgroundDoping,” although a periodic variation in this doping can be introduced using the “TreeRing” parameters. The doping level is assumed to be modified by the introduction of implants from the bottom side. There are several options for these doping profiles, including a square profile of a given depth or a sum of N Gaussian profiles. The use of one or two Gaussian profiles has been found to accurately reproduce the measured doping profiles on commercial CCDs. For more details on this, see Ref. 16.

Setting up the initial conditions in the periodic pixel array is straightforward and is specified by a relatively small number of parameters that describe the gate voltages and doping levels. Rather than going through these in detail, the reader is referred to the  Appendix A or the “pixel-itl” and “pixel-e2v” examples at Ref. 2.

Setting up the initial conditions in non-periodic regions outside the pixel array is straightforward but more laborious than setting up the pixel arrays. The extents, dopings, applied voltages, and quasi-Fermi levels need to be specified for each region. At present, only rectangular regions are supported. Also, it is assumed that the same doping profiles that are used in the pixel array are used in the surrounding circuitry, so the only options for doping profiles are the channel doping, the channel stop doping, and no doping. Examples of simulations setting up non-periodic regions are the “edge.cfg,” “trans.cfg,” and “io.cfg” files described in Ref. 2, and the results of these simulations are detailed in Sec. III.

In this section, we give a brief description of the methods that are used to solve Poisson’s equation on the grid. We are trying to solve the following equation, where $φ$ is the potential, $ρ$ is the charge density, and $ϵSi$ is the dielectric constant of silicon:

Each of the partial derivatives can be discretized on a rectangular grid with grid spacing h as follows:

Giving for the discretized Poisson’s equation,

This can be turned into an iterative equation, and basically one just iterates until convergence,

which we write in shorthand as follows:

where we define

and we have split $ρ$ into a fixed charge density $ρf$ and a mobile charge density $ρm$⁠. However, $ρm$ is a highly nonlinear function of the potential $φ$⁠, as described below. In quasi-equilibrium, the drift current $JE$ and the diffusion current $JD$ are equal, giving a net current of zero, so we can write (see Sze,¹⁷ for example),

Here, $qe$ is the electron charge, $μn$ is the electron mobility, $n$ is the electron density in electrons per unit volume, and $Dn$ is the electron diffusion coefficient. This gives

We also know from the Einstein relations that the mobility and diffusion coefficient are related by the following equation, where $k$ is Boltzmann’s constant and $T$ is the absolute temperature:

so, integrating both sides,

We take the constant of integration into the exponential and define the quasi-Fermi level $φF$ in terms of the intrinsic carrier density $ni$⁠, giving

So we need to solve the following equation for $φ$⁠, where $φF$ is a constant:

which when discretized is

Because of the strong non-linearity, simply iterating is numerically unstable. The method that works, described in detail by Rafferty et al.¹⁸ and building on the work of Gummel,¹⁹ is to take this last equation as a non-linear equation for $φ(n+1)$ in terms of $φn$ and run Newton’s method “inner loop” to find $φ(n+1)$ at each grid point. Then we iterate to convergence as before. This allows us to simultaneously solve for the potential and the carrier density. Of course, we just described the electron density here, but there is a similar equation for holes, but with opposite signs. Figure 3 shows an example of varying $φF$ on the solution. Note that the quasi-Fermi level is constant in each region containing mobile carriers. In the CCD, each collecting well contains a different number of mobile carriers, so $φF$ is constant in each well but is different from well to well. However, when simulating the device, instead of knowing the value of $φF$⁠, we typically know the number of electrons in each well. So how do we translate from the known number of electrons to the unknown value of $φF$? The code provides two methods, selected by the value of the parameter “ElectronMethod.” With this parameter set to 1, a test simulation is run where the parameter $φF$ (called QFe in the code) is varied through a range. For each value of QFe, we integrate over a pre-determined region to count the number of electrons in the well, building a look-up table of the number of electrons vs QFe. The results for a few values of QFe are shown in Fig. 3. Then the code interpolates to determine the value of QFe that gives the appropriate number of electrons. In practice, one can get close to the desired number of electrons, but the non-linearity causes variations from the desired number, so a second method was developed. When “ElectronMethod” has a value of 2, what is done is to place the correct number of electrons in the well, uniformly distributed in the center of the well. The code then moves the electrons around until the value of QFe is constant in the well. This allows one to get the correct number of electrons in the well without needing to know the value of the quasi-Fermi level.

There is one more complication. As discussed in Sec. II A, a non-linear Z-axis is used to concentrate grid cells in the bottom region where the potentials and charge densities are changing most strongly. In principle, any smooth function can be used for Z-axis mapping. Here we chose an easily differentiable polynomial function of the following form, where z is the linear z coordinate and zp is the non-linear coordinate, which is the actual value used in solving and plotting. Here, $TSi$ is the silicon thickness, and NZExp is a user-defined constant that controls the non-linearity, as described in Fig. 2,

The non-linear z-axis modifies Poisson’s equation from

to

In practice, the additional partial derivatives can be pre-computed, so this is simply accounted for in the discretized equations and has only a minor impact on the speed of iteration.

It is well known that multi-grid methods’ speed convergence of solutions to Poisson’s equation by getting correct solutions to the long-wavelength modes at a coarser grid where convergence is much more rapid. There is a wealth of literature on the subject, and Briggs²⁰ or Press²¹ give excellent summaries. The basic idea is shown in Fig. 4. In practice, in this code, we have adopted a simpler method. Rather than use “Restriction” to propagate the boundary conditions down to the coarser grid, we simply set up the boundary conditions on all of the sub-grids at the outset of the problem. In addition, we have found that there is little value in running coarser grids than $403$ because at this resolution the problem converges very rapidly. So for a typical problem that has perhaps $3203$ grid cells, we define the finest grid and three subgrids, with the coarsest grid having $403$ grid cells. We then set up the boundary conditions on all of the subgrids, iterate the coarsest grid, “Prolongate” the solution up to the next grid, and continue until we have reached the finest grid.

At this point, it is appropriate to discuss the problem of convergence. Convergence of the SOR algorithm is notoriously slow. Multi-grid methods help a great deal; however, care must still be taken to ensure that the solution has converged adequately for your problem. The parameter “ncycle” controls the number of iterations taken at the finest grid. Each coarser grid increases the number of iterations by a factor of 4. So, for example, a typical problem like one of the “pixel” examples, which has ncycle = 128, the coarsest grid has 1/8 the resolution, and will run $128×43=8192$ SOR cycles at the coarsest grid. Figures 5 and 6 show the convergence of a typical problem. For most problems, a value of ncycle = 64 is adequate. The problems most dependent on accurate convergence have proven to be the pixel distortion simulations like those in Sec. III A. Since we are dealing with very small deviations in the pixel shapes due to the BF effect, it is important to make sure the results have converged.

Compared to modern digital electronic devices, CCDs have relatively large geometries and relatively low electric fields. Because of this, a simple drift-diffusion model for carrier transport was deemed to be sufficient. The good agreement between the simulations and measured data shown in Sec. III A shows that this decison is justified. The basic scheme for modeling carrier transport is shown in Fig. 7. Electrons are assumed to have lattice collisions on a time scale $τ$⁠, which is on the order of picoseconds. At each collision, the electron is assumed to pick up a thermal velocity $Vth$ which is in a random direction. In addition to this thermal velocity, it also has a drift velocity given by $Vdrift=μE$⁠, where the mobility $μ(E,T)$ is calculated as a function of electric field and temperature using the mobility model of Jacobini²². These two velocities are added vectorially and it travels in this direction at this velocity for a time $δt$ until the next collision, and this continues until the electron reaches the bottom. The electron path is logged in the *\_Pts.dat file. If the parameter “LogPixelPaths” is zero, only the initial and final positions are logged. If this parameter is one, the entire path is logged. The thermal velocity has a multiplier (“DiffMultiplier”) which can be used to tune the amount of diffusion. If this is set to zero, diffusion is turned off, with the impact as shown in Fig. 8. A value $DiffMultiplier=2.30$ has been found to accurately reproduce the amount of diffusion seen in $Fe55$ data (see Sec. III B). Since the value of $me∗$ in the code is the bare electron mass, this value is equivalent to an electron effective mass of about 0.19. This is somewhat low, as Green²³ finds a value of 0.27. This discrepancy may be due to the limitations of the simple drift-diffusion model of carrier transport. This value of DiffMultiplier is easily adjusted by the user, however.

The initial electron locations in X and Y can be determined in a number of ways, as determined by the “PixelBoundaryTestType” parameter. These include an equally spaced grid, a Gaussian spot, a random location within a boundary, an $Fe55$ event, or reading in a list of locations. The starting location in Z can be either specified or calculated given a filter band.

Once the simulation has run, the potential, electric fields, and charge carrier densities are available throughout the simulation volume. What one chooses to visualize depends on the problem being studied. Here we have chosen three examples (Figs. 9–11) of the type of data that are available.

To successfully model a CCD, an accurate physical characterization of the CCD is needed. Details such as dopant densities, layer thicknesses, and physical dimensions need to be known, at least approximately. For modeling the CCDs from the two vendors that will be used in the LSST camera, we obtained detailed physical measurements, which have been described in Ref. 16. This has helped to make the simulations described in Secs. III A–III F as physically realistic as possible. Most of the plots in this section can be reproduced using the examples at the Poisson_CCD code site described in Ref. 2.

Regarding the parameters used in the simulations, there has been optimization of the numeric parameters, such as grid size, number of iterations, successive over-relation factor, etc., to ensure convergence and accuracy. However, the physical parameters, such as physical dimensions, voltages, doping levels, oxide thicknesses, etc., have all been directly measured. In particular, these parameters have not been tuned to improve the fit with measured data shown in the following examples.

As has been extensively discussed in the literature,^11–14 as charge builds up in the central region of bright objects, the stored charge repels additional incoming charge and broadens the profile of these objects. The impact of the stored charge on the pixel shapes can be determined by measuring the pixel–pixel covariances on a large number of flat images.^11,15 These covariances are calculated from a large number of flat pairs of varying intensity (see Ref. 15 for example) as

where $fi,j$ is the difference in flux between the two flats at pixel i,j and $Npix$ is the number of pixels summed over.

We have generated these data on a large number of flat pairs on LSST CCDs, measured on the UC Davis LSST beam simulator,^14,24 and would like to compare these results to simulations. For a dataset of 100 flat pairs, each with $16×106$ pixels with an average signal level of 50 000 electrons, this is approximately $1014$ electrons. Directly simulating these data is out of the question. However, we have found a simple way to run a single simulation that reproduces the measured pixel covariances. In order to do this, we first simulate a situation where one pixel has a fixed amount of charge (typically 100 000 electrons), and all surrounding pixels are empty. After solving for the potential and the resulting electric field, we can track electrons down through the silicon. As the electrons travel down through the silicon under the influence of the electric field, they eventually end up in one of the collecting wells. A binary search is used to find the bifurcation points, where electrons on one side of the bifurcation point end up in one pixel and electrons on the other side of the bifurcation point end up in an adjacent pixel. This binary search is performed with diffusion turned off in order not to introduce a stochastic element into the electron paths. These bifurcation points are identified as the pixel boundaries. This allows us to characterize the distortion in the pixel boundaries, which results from the central pixel charge. A typical result of this process is shown in Fig. 12. The number of vertices used to define the distorted pixel shapes is an adjustable parameter. A minimum of 12 vertices is needed to adequately define the distortion (4 corners plus 2 points per edge). The distorted pixels shown in Fig. 12 and used in Fig. 13 were simulated with a total of 132 vertices (4 corners plus 32 points per edge). The distorted pixel shapes that result are what gives rise to the BF effect, with the central pixel losing area, which is gained by the surrounding pixels.

We find that the area distortions accurately capture the measured pixel–pixel covariances. Figure 13 shows the agreement between the measured pixel–pixel covariances on flat field images and the simulated area distortions, as measured and as simulated on LSST CCDs from both CCD vendors. The agreement is quite good. The asymmetry of the nearest neighbor pixels is correctly modeled, and the simulated values agree with the measurements within the statistical errors. These simulations are run with the “pixel-itl.cfg” and “pixel-e2v.cfg” examples given in Ref. 2. More details on these measurements and simulations are given in Refs. 14 and 25.

CCDs are routinely characterized by exposing the CCD to an $Fe55$ source. (see Ref. 26, for example). The radioactive decay produces x rays with a known energy, with the $Kα$ peak being the strongest peak, typically producing 1620 hole–electron pairs when photoelectrically absorbed in the silicon. These carriers then propagate down to the collecting wells, where they are collected and counted. Because of diffusion, the carriers, which are initially produced in a small volume, spread out and occupy several pixels. To simulate these events, a special module was written. Normally photoelectrons propagate one at a time, without influence from neighboring carriers. But the carriers produced in the $Fe55$ event are produced in a short time, so interactions between the carriers might be important. The code takes the like carrier repulsion and opposite carrier attraction into account, and the parameters “Fe55ElectronMult” and “Fe55HoleMult” can be used to turn off or modify this interaction if desired. Figure 14 shows a typical event, and Fig. 15 shows stacked pixel maps compared between measurements and simulations. The spread of the charge cloud due to diffusion is well modeled. This simulation is run with the “fe55.cfg” example described in Ref. 2.

A large number of simulations have been run to better understand saturation and blooming in the ITL STA3800C device. These simulations have been very helpful to understand the physics of the device because in the simulator one can get information that is simply not accessible to measurement. Figure 16 shows the distinction between the “bloomed full well” condition, where charge begins to bloom above the charge storage regions, and the “surface full well” condition, where the charge blooms along the silicon surface. This effect is discussed in detail in Janesick.²⁶ The simulations in Figs. 17 and 18 reproduce these conditions, and these simulations illustrate the difference between these conditions. Measurements of saturated spots in the surface full well condition, shown in Fig. 19, also show that in the surface full well condition, the charge is lost to traps at the silicon–silicon dioxide interface. Thus, it is apparent that the surface full well condition is to be avoided.

We can go further and quantitatively reproduce measurements of the onset of saturation as a function of parallel low and high voltages, as shown in Fig. 20. By quantifying the barrier height between the storage wells, we show that saturation occurs when the barrier height drops below a certain value. The fit is good except in the strong surface full well condition because the charge loss that occurs is not included in the simulations. It would be possible to modify the simulations to include this effect, but since the surface full well condition is to be avoided, this was deemed to be not worth the effort.

In addition to simulations of the pixel arrays, simulations can be done of the peripheral circuitry, as shown in Secs. III D–III F. At the edges of the pixel array, lateral electric fields from the surrounding circuitry can introduce pixel boundary shifts, which lead to measurable astrometric shifts. This effect has been characterized on the UC Davis LSST Optical Simulator.^24,27 We have then built a simulation of a portion of the pixel array that extends to the chip edge to compare the measurements to the simulations. In addition to the pixel array, it is necessary to build into the simulation the appropriate “fixed voltage regions” at the edge of the chip. Figure 21 shows the setup of this simulation. The bending of the equipotential lines near the edge of the pixel array betrays the presence of a lateral electric field thats deflects incoming electrons. Figure 22 shows these simulated paths (with diffusion turned off), and the edge deflection is apparent. In addition to this deflection, there is a second effect that influences the astrometric shift, which is that as the measured spot begins to “fall off” the edge of the array, there is a resulting shift in the opposite direction. These two competing effects lead to the astrometric shifts as seen in Fig. 23. The shift has been characterized for a number of different measurement conditions. The simulation, while not perfect, captures all of the trends correctly. These simulations are run with the “edge.cfg” example in Ref. 2.

This section shows simulations that were done to model the output transistor of the ITL STA3800C. Figure 24 shows the simulated and measured $Id−Vg$ characteristics. We obtained a relatively good fit of the transistor turn-on. This simulation is run with the “trans.cfg” example in Ref. 2.

Some other tests have been run which have given results that are qualitatively reasonable, but which have not been compared with quantitative measurements, and three of these are reviewed here. The first of these are known as “tree rings.” As is well known, periodic dopant variations introduced during the growth of the silicon boule can lead to measurable variations in flat fields and introducing astrometric pixel shifts (see, for example, Ref. 10). This effect has been successfully simulated, as shown in Fig. 25. This simulation is run with the “treering.cfg” example in Ref. 2.

A second interesting simulation is of the backside substrate connection $VBB$⁠. It is often questioned how this bias voltage, which is only connected to the CCD frontside, is conducted to the backside when the CCD is fully depleted. The answer is that there is an undepleted region near the chip edge that serves this purpose. It is interesting to simulate this connection, and the result is shown in Fig. 26.

While the simulator solves for the condition of the CCD in equilibrium and does not do transient simulations, it is possible to simulate transient effects by repeatedly solving for the state of the CCD with small incremental changes to the boundary conditions. These can then be stitched together to form a movie of the results. An example of the parallel charge transport is shown in Fig. 27. Several movies constructed in this way are available in Ref. 2.

We have presented a software package optimized for simulating astronomical CCDs. The code has been optimized to be as physically realistic as possible and accurately simulates many aspects of CCD behavior. The package is open source and freely available² and has been validated against a number of different types of CCD measurements. It has proven useful for analyzing sensor effects in the CCDs being used to construct the LSST focal plane, and insights gained from running these simulations are being incorporated into the software stack to be used for instrument signature removal in the LSST images. It is hoped that it will find additional uses in the future.

Kirk Gilmore’s help in setting up the hardware and software for the CCDs from both vendors has also been invaluable. Perry Gee has provided key support in software and networking. Merlin Fisher-Levine has patiently helped with our understanding and implementation of the LSST data management software stack. David Kirkby contributed code for including different filters. Financial support from the Department of Energy (DOE) under Grant No. DE-SC0009999 and Heising-Simons Foundation under Grant No. 2015-106 are gratefully acknowledged. This paper has undergone internal review in the LSST Dark Energy Science Collaboration. The internal reviewers were Andrew Rasmussen, Douglas Tucker, and Dan Weatherill, and their inputs are much appreciated.

The data that support the findings of this study are either available within the article or in the “measurements” directory at the code repository, Ref. 2.

The authors have no conflicts of interest to disclose.

C.L. developed the Poisson_CCD code, ran most of the studies described here, and drafted the paper. A.B. carried out the study of astrometric shifts at array edges and provided test and debug support. J.A.T. supervised the project.

The figures in this articles were generated with v1.0 of the Poisson_CCD code provided in Ref. 2. There are a total of 14 examples included with the code. Each example is in a separate directory in the data directory and has a configuration file of the form *.cfg. The parameters in the *.cfg files are commented to explain (hopefully) the purpose of each parameter, and a detailed listing of all configuration parameters is in the  Appendix A. Python plotting routines are included with instructions below on how to run the plotting routines and the expected output. The plot outputs are placed in the data/*/plots files, so you can see the expected plots without having to run the code. If you edit the .cfg files, it is likely that you will need to customize the Python plotting routines as well. The estimated run times given here are what was seen on a dual-core 1.4 GHz Intel Core i5 laptop computer.

• Example 1: data/smallpixel/smallpixel.cfg

  1. Purpose: A single pixel and surroundings. The central pixel contains 100 000 electrons. No electron tracking or pixel boundary plotting is done. The subgrids are saved so one can look at convergence. This is useful for getting things set up rapidly

  2. Syntax: src/Poisson data/smallpixel/smallpixel.cfg

  3. Expected run time: <1 min.

  4. Plot Syntax: python pysrc/Poisson_Small.py data/smallpixel/smallpixel.cfg 0

  5. Plot Syntax: python pysrc/Poisson_Convergence.py data/smallpixel/smallpixel.cfg 0

• Example 2: data/pixel0/pixel.cfg

  1. Purpose: A $9×9$ grid of pixels at low resolution (ScaleFactor = 1). The central pixel contains 100 000 electrons. No electron tracking or pixel boundary plotting is done.

  2. Syntax: src/Poisson data/pixel0/pixel.cfg

  3. Expected run time: $≈$ 1 min.

  4. Plot Syntax: python pysrc/Poisson_Plots.py data/pixel0/pixel.cfg 0

  5. Plot Syntax: python pysrc/ChargePlots.py data/pixel0/pixel.cfg 0 2

• Example 3: data/pixel-itl/pixel.cfg

  1. Purpose: A $9×9$ grid of pixels at higher resolution (ScaleFactor = 2). The central pixel contains 100 000 electrons. The parameters are set up for the ITL STA3800C CCD. This should give physically meaningful results and is what was used in the published papers. After solving Poisson’s equation, electron tracking is done to determine the pixel distortions.

  2. Syntax: src/Poisson data/pixel-itl/pixel.cfg

  3. Expected run time: $≈$ 30 min.

  4. Plot Syntax: python pysrc/Poisson_Plots.py data/pixel-itl/pixel.cfg 0

  5. Plot Syntax: python pysrc/ChargePlots.py data/pixel-itl/pixel.cfg 0 2

  6. Plot Syntax: python pysrc/VertexPlot.py data/pixel-itl/pixel.cfg 0 2

  7. Plot Syntax: python pysrc/Area_Covariance_Plot.py data/pixel-itl/pixel.cfg 0

• Example 4: data/pixel-e2v/pixel.cfg

  1. Purpose: A $9×9$ grid of pixels at higher resolution (ScaleFactor = 2). The central pixel contains 100 000 electrons. The parameters are set up for the E2V CCD250 CCD. This should give physically meaningful results and is what was used in the published papers. After solving Poisson’s equation, electron tracking is done to determine the pixel distortions.

  2. Syntax: src/Poisson data/pixel-e2v/pixel.cfg

  3. Expected run time: $≈$ 30 min.

  4. Plot Syntax: python pysrc/Poisson_Plots.py data/pixel-e2v/pixel.cfg 0

  5. Plot Syntax: python pysrc/ChargePlots.py data/pixel-e2v/pixel.cfg 0 2

  6. Plot Syntax: python pysrc/VertexPlot.py data/pixel-e2v/pixel.cfg 0 2

  7. Plot Syntax: python pysrc/Area_Covariance_Plot.py data/pixel-e2v/pixel.cfg 0

• Example 5: data/pixel1/pixel.cfg

  1. Purpose: A $9×9$ grid of pixels at higher resolution (ScaleFactor = 2). The central pixel contains 100 000 electrons. The parameters are set up for the ITL STA3800C CCD. This is included as an example of ElectronMethod = 1, which iterates to find the value of the electron Quasi-Fermi level. After solving Poisson’s equation, electron tracking is done to determine the pixel distortions.

  2. Syntax: src/Poisson data/pixel1/pixel.cfg

  3. Expected run time: $≈$ 45 min.

  4. Plot Syntax: python pysrc/Poisson_Plots.py data/pixel1/pixel.cfg 0

  5. Plot Syntax: python pysrc/ChargePlots.py data/pixel1/pixel.cfg 0 2

  6. Plot Syntax: python pysrc/VertexPlot.py data/pixel1/pixel.cfg 0 2

  7. Plot Syntax: python pysrc/Area_Covariance_Plot.py data/pixel1/pixel.cfg 0

• Example 6: data/smallcapf/smallcap.cfg

  1. Purpose: A simple field oxide capacitor for testing convergence.

  2. Syntax: src/Poisson data/smallcapf/smallcap.cfg

  3. Expected run time: $≈$ 1 min.

  4. Plot Syntax: python pysrc/Poisson_CapConvergence.py data/smallcapf/smallcap.cfg 0

• Example 7: data/smallcapg/smallcap.cfg

  1. Purpose: A simple gate oxide capacitor for testing convergence.

  2. Syntax: src/Poisson data/smallcapg/smallcap.cfg

  3. Expected run time: $≈$ 1 min.

  4. Plot Syntax: python pysrc/Poisson_CapConvergence.py data/smallcapg/smallcap.cfg 0

• Example 8: data/edgerun/edge.cfg

  1. Purpose: A simulation of the pixel distortion at the top and bottom edges of the ITL STA3800C CCD.

  2. Syntax: src/Poisson data/edgerun/edge.cfg

  3. Expected run time: $≈$ 15 min.

  4. Plot Syntax: python pysrc/Poisson_Edge.py data/edgerun/edge.cfg 0

  5. Plot Syntax: python pysrc/Pixel_Shift.py data/edgerun/edge.cfg 0

• Example 9: data/iorun/io.cfg

  1. Purpose: A simulation of the output transistor region of the ITL STA3800C CCD, including the last few serial stages.

  2. Syntax: src/Poisson data/iorun/io.cfg

  3. Expected run time: $≈$ 5 min.

  4. Plot Syntax: python pysrc/Poisson_IO.py data/iorun/io.cfg 0

  5. Plot Syntax: python pysrc/ChargePlots_IO.py data/iorun/io.cfg 0

• Example 10: data/treering/treering.cfg

  1. Purpose: A simulation of a $9×9$ pixel region with a sinusoidal treering dopant variation introduced, so one can see the resulting pixel distortions.

  2. Syntax: src/Poisson data/treering/treering.cfg

  3. Expected run time: $≈$ 20 min.

  4. Plot Syntax: python pysrc/VertexPlot.py data/treering/treering.cfg 0 4

• Example 11: data/fe55/pixel.cfg

  1. Purpose: A simulation of a $5×5$ pixel region with a single Fe55 event above the center pixel. To determine the distribution of pixel counts, one needs to run many of these simulations with a random center point.

  2. Syntax: src/Poisson data/fe55/pixel.cfg

  3. Expected run time: $≈$ 5 min.

  4. Plot Syntax: python pysrc/Poisson_Fe55.py data/fe55/pixel.cfg 0 40

• Example 12: data/transrun/trans.cfg

  1. Purpose: A simulation of the output transistor of the ITL STA3800C CCD. In order to plot the I–V characteristic, this needs to be re-run with different gate voltage values.

  2. Syntax: src/Poisson data/transrun/trans.cfg

  3. Expected run time: $≈$ 10 min.

  4. Plot Syntax: python pysrc/Poisson_Trans.py data/transrun/trans.cfg 0

  5. Plot Syntax: python pysrc/ChargePlots_Trans.py data/transrun/trans.cfg 0

  6. Plot Syntax: python pysrc/MOSFET_Calculated_IV.py data/transrun/trans.cfg 0

• Example 13: data/satrun/sat.cfg

  1. Purpose: A simulation of the saturation level, as determined by the barrier lowering. Again, this needs to be run multiple times with varying parallel voltages.

  2. Syntax: src/Poisson data/satrun/sat.cfg

  3. Expected run time: $≈$ 20 min.

  4. Plot Syntax: python pysrc/Poisson_Sat.py data/satrun/sat.cfg 0

  5. Plot Syntax: python pysrc/ChargePlots.py data/satrun/sat.cfg 0 9

  6. Plot Syntax: python pysrc/Barrier.py data/satrun/sat.cfg 0

  7. Plot Syntax: python pysrc/Plot_SatLevel.py data/satrun/sat.cfg 0

• Example 14: data/bfrun/bf.cfg

  1. Purpose: This builds up a Gaussian spot in the center of a $9×9$ pixel region by sequentially solving Poisson’s equation, adding electrons, and repeating. As written, this is done 80 times, building up the central pixel charge to about 80 000 electrons.

  2. Syntax: src/Poisson data/bfrun/bf.cfg

  3. Expected run time: $≈$ 8 h.

  4. Plot Syntax: python Poisson_Plots.py data/bfrun/bf.cfg 80

  5. Plot Syntax: python ChargePlots.py data/bfrun/bf.cfg 80 9

  6. Plot Syntax: python Plot_BF_Spots.py data/bfrun/bf.cfg 0 NumSpots 80 (Assumes this was run NumSpots times with random center locations)