text:
# Modeling Charge Transfer Inefficiency in the Chandra Advanced CCD Imaging Spectrometer

## 1 Introduction

This paper is the third in a suite of articles outlining recent efforts to improve the modeling of CCDs used in the Chandra X-ray Observatory (hereafter Chandra). The Advanced CCD Imaging Spectrometer (ACIS) instrument employs bulk front-illuminated (FI) CCDs and back-illuminated (BI) CCDs to give good spectral resolution and good 0.2–10 keV quantum efficiency . These devices, designed and manufactured at MIT’s Lincoln Laboratories , couple with the Chandra mirrors to give users the ability to do spatially-resolved X-ray spectroscopy on spatial scales comparable to ground-based visual astronomy.

In the first paper in this series , Pavlov and Nousek developed a solution to the diffusion equation, including recombination, to describe the process of charge diffusion in field-free regions of CCDs. In the second paper (, hereafter T01), Townsley et al. describe the latest Pennsylvania State University (PSU) CCD simulator, which includes the results from Pavlov and Nousek and many recent innovations in CCD modeling made by our ACIS team colleagues at the Massachusetts Institute of Technology (MIT) Center for Space Research (CSR) , , . T01 also gives introductory material such as the geometry of ACIS devices and the event grading scheme used; readers unfamiliar with the operation of X-ray CCD cameras should see that paper and references therein for details. This paper narrows the focus of T01 to a calibration complication with the ACIS CCDs – charge transfer inefficiency (CTI).

Below we discuss the pronounced CTI effects observed in ACIS flight data. CTI is a phenomenon in which some of the charge resulting from an X-ray photon’s interaction with the CCD is lost as that charge is transferred across the device to the readout nodes. Thus CTI changes the gain across the device. It also causes spatial redistribution of the charge in the pixel neighborhood of an event, as charge traps capture electrons and release them on short timescales during the readout process. This results in charge being redistributed in surrounding pixels within the region recognized as a single event. Charge loss can be substantial (equivalent to the reduction of apparent X-ray energy of over 100 eV), depending on photon energy and the event’s position with respect to the CCD readout nodes.

Due to the unanticipated forward scattering of charged particles (probably mainly $``$100 keV protons , ) by the Chandra mirrors onto the ACIS CCDs, the FI devices suffered degraded performance on-orbit, most pronounced at the top of the devices, near the aimpoint of the imaging array. This radiation damage created charge traps with a variety of time constants that increased the parallel CTI of these devices markedly, resulting in the equivalent of several years of damage in only a few weeks of operation.

As a result of their manufacturing process, the ACIS BI devices have always shown modest serial and parallel CTI effects, but their geometry has protected them from severe radiation damage on-orbit so they show no additional CTI. Although operational changes have halted the rapid degradation of the FI devices, the space environment is likely to increase CTI effects slowly on all chips throughout the life of the mission. Prior to launch we began development of a CTI model for BI devices and a data correction process to remove CTI effects. We have extended this model and corrector to address both BI and FI CTI effects at two operating temperatures.

To ensure that no water ice could collect on Chandra’s instruments, the focal plane temperature of the Observatory was kept at the relatively warm setting of $``$110C from mid-September 1999 to the end of January 2000. Thereafter, it was determined that the Observatory structures had outgassed enough that it was safe to lower the focal plane temperature to $``$120C. At this new colder temperature, the trap population that causes the row-dependent energy resolution in FI devices is partly suppressed , . This step greatly improved the performance of the FI devices by ameliorating the effects of CTI. The BI CTI is very slightly worse at $``$120C than at $``$110C but the improvements in the FI performance completely justify operating at $``$120C.

Given the existing and potential complication CTI causes for ACIS calibration, we have developed a model for CTI and incorporated it into our Monte Carlo CCD simulator. CTI is parameterized as a function of the charge in each pixel, including the effects of de-trapping (charge trailing), shielding within an event (charge in the leading pixels of the 3$`\times `$3 pixel event island protect the rest of the island by filling traps), and non-uniform spatial distribution of traps. This technique partially recovers the degraded energy resolution near the top of the FI chips and it reduces the position dependence of gain and grade distributions. By correcting the grade distributions as well as the event amplitudes, we can improve the response of the FI chips, both in quantum efficiency at high energies and in background rejection.

First we will illustrate the effects of CTI on ACIS performance. Our CTI model is described in Section 3; simulations showing its ability to reproduce CTI are then presented. Next we describe the program that we have developed to partially remove CTI effects from ACIS data and examine the performance of this “CTI corrector.” A preliminary version of the CTI model and corrector was described by Townsley et al. .

## 2 Illustrating CTI with ACIS Data

Due to crystal defects caused in the manufacturing process or by radiation damage, silicon-based detectors such as CCDs can have charge traps, where charge is inefficiently transferred. There is a large body of work on the reasons for and results of charge traps in silicon: see Janesick et al.  for more details on bulk and radiation-induced traps in X-ray astronomy CCDs; see Stetson  for an illustration of CTI effects on CCDs used for visual astronomy. For our purposes, we can think of charge traps simply as potential wells that trap charge on a timescale short compared to the timescale for charge transfer from one pixel to the next. The traps then release charge on some exponential timescale that may be comparable to the charge transfer time (see Gendreau et al.  for a mathematical development of these concepts). Prigozhin et al. have shown that the charge traps in ACIS FI devices consist of multiple trap species characterized by different decay time constants .

Charge release from traps with short time constants results in a phenomenon that we call “charge trailing” – charge in the final detected event appears shifted from its native pixel into the adjacent pixel farther from the readout node. Since events are graded based on the distribution of charge in a $`3\times 3`$ pixel neighborhood around the brightest pixel (see Figure 2 in T01), charge trailing causes events to be recognized with different grades than their original charge distribution would have generated – we call this phenomenon “grade morphing.” Only part of the trapped charge appears trailed into the adjacent pixel(s); the rest of it is lost to the event because the charge in that event was transferred away from the interaction site before the other traps released all of their captured charge.

Our data show that the charge loss and trailing vary with photon energy in a complicated way. This is due to the nature of the traps convolved with the size of the charge clouds. In the BI device, it is also the result of the device geometry. For a BI device, low-energy photons can create large-diameter charge clouds, since many of those photons interact near the surface of the device and suffer substantial charge diffusion as the charge cloud transits the depletion region. High-energy photons create large-diameter initial charge clouds, which also can suffer substantial diffusion, depending on the interaction depth. As charge clouds get bigger, they spread over more pixels, so they can encounter more charge traps as they are transferred across the device.

### 2.1 On-orbit Calibration Data

During radiation belt passage, the Chandra Science Instrument Module is configured so that ACIS is out of the telescope’s field of view. This protects ACIS from further radiation damage and orients the instrument under its calibration source  (called the “External Calibration Source” or ECS here). As a way to monitor the CTI, ECS data are collected and telemetered every orbit. We have combined these ECS data to calibrate our model for CTI and to check the fidelity of our CTI corrector that is based on that model.

The ECS spectrum is well-understood and was described in detail by the MIT/ACIS team in the ACIS Calibration Report . Using their line identifications, Table 1 shows the main spectral features apparent in our composite ECS dataset and notes the ones used in this analysis.

In order to achieve the highest fidelity in our model calibration, we have combined all available data (as of early 2001) from the ACIS ECS at a given focal plane temperature to generate a corrector appropriate for any dataset taken at that temperature. We have checked that the CTI is stable enough over the time periods in question that this approach is sound and that the resulting corrector performs acceptably on all ECS datasets in the sample.

Example spectra for Amplifier 3 (Node D) on the FI device I3 and Amplifier 0 (Node A) on the BI device S3, using a large sample of the appropriate ECS data taken at $``$120C and all rows of the devices, are shown in Figure 1. These spectra reflect the work of the Chandra X-ray Center (CXC) to remove the largest effect of CTI (illustrated in Figures 2 and 4 below), the row-dependent gain. They represent the standard processing applied to all ACIS data and form the fundamental comparitors for our CTI corrector’s performance.

The data were filtered to keep only events with certain “standard” grades. The ACIS event grading scheme can be related to that of the Japanese X-ray satellite ASCA , which is familiar to X-ray astronomers. The standard ACIS grade filtering scheme includes ASCA-like Grades 0, 2, 3, 4, and 6 (“g02346”, ) and was developed (in the absence of CTI) to maximize both sensitivity and spectral resolution (see the “Chandra Proposers’ Observatory Guide,” ). Virtually all standard data products employ this grade filter.

### 2.2 FI CTI

Figure 2 illustrates the effects of CTI on the ACIS I3 chip at focal plane temperatures of $``$110C (left) and $``$120C (right) by showing the row-dependent charge loss and energy resolution for the spectral lines in the ECS. The events are binned by position and energy using the Event Browser software  into an “image” to make the features easier to see.

Note that the charge loss gets more severe at higher event energies (the loci in Figure 2 steepen with line energy, or larger DN values). The data are displayed using the digital value (“Data Number,” DN) of the event amplitude because the standard data processing pipeline applies a correction factor to the event energy in eV (a quantity derived from the event amplitude) as a first-order correction for CTI. Since we wish to illustrate this charge loss effect here, we display the data in their native units. For a rough approximation of a line’s event amplitude in DN, divide its energy in eV by 4. We use I3 as the example FI CCD because it contains the aimpoint for the ACIS imaging array. The other I3 amplifiers and other FI devices show similar charge loss effects.

The left-hand panels of Figure 2 are representative of the state of the CTI at $``$110C and apply to data gathered in the first five months of the mission. The improvement gained by operating at $``$120C is apparent in the right-hand panels of Figure 2 which show less line blurring at large row numbers.

The lower panels of Figure 2 show expanded views of the Mn and Ti lines. Note the non-linearity of the charge loss with row number that lends a slight S-shape to the curve: the slope is flatter in the middle third of the CCD than at low and high row numbers. Here the increase in linewidth toward the top of the chip (high row number) is easily seen.

Due to the non-uniform spatial distribution of traps, ACIS CCDs also exhibit column-to-column gain variations that increase in amplitude with energy and row number. These variations contribute substantially to the spectral resolution degradation. Figure 3 illustrates these variations using ECS data from Amplifier 3 of the I3 chip. The data are displayed in energy space now, using the standard processing techniques to remove the gross row-dependent gain variations so that the column-wise variations are more apparent.

### 2.3 BI CTI

Although CTI in ACIS BI devices is caused by a different mechanism than in the FI devices, it causes similar non-uniformities in the data. The two major differences are that BI devices show CTI in the serial as well as the parallel registers and the framestore is also affected. This means that parallel CTI acts over 2048 parallel transfers in BI devices rather than just 1024 parallel transfers as in the FI devices. The amplitude of serial CTI is much larger than parallel CTI in the BI devices, but there are at most 256 serial transfers. Thus both parallel and serial CTI must be modeled in BI devices and the interaction of these two effects results in spatially-complex charge loss and grade morphing behavior.

Figure 4 illustrates CTI in the BI device S3, using only data from Amplifier 0 at $``$120C. BI CTI is only slightly worse at $``$120C than at $``$110C, so only $``$120C data are presented here. The left-hand columns illustrate serial CTI over the full range of the ECS data (top) and for the major spectral lines (bottom). The amplifier readout is at Column 0. The right-hand columns show parallel CTI for the same spectral ranges. Parallel CTI is illustrated using only data from the 64 columns closest to the readout node in order to minimize the effects of serial CTI. As in Figure 2, the data are displayed in their native units (DN) and only standard grades are used.

Again due to the non-uniform spatial distribution of charge traps, BI CCDs show column-to-column gain variations that increase in amplitude with energy and degrade the spectral resolution. These variations are visible in the left-hand panels of Figure 4. Also apparent there are charge loss variations of lower spatial frequency.

### 2.4 Grade Morphing

Since CTI causes charge to be shifted into adjacent pixels, thus causing grade morphing, stringent grade filtering results in non-uniform quantum efficiency across the device. Figure 5 shows the number of detected events in the $``$6 keV Mn K lines in ECS data, using standard grades (g02346, left) and the complementary Grades 1, 5, and 7 (right). The top plots are for the FI device I3 at $``$110C; at $``$120C the FI CTI is reduced and the quantum efficiency is quite uniform across the device for the standard grades. The bottom plots are for the BI device S3 at $``$120C; the behavior at $``$110C is similar. Note the complex interaction between parallel and serial CTI in the BI device: events morph into Grade 7 due to parallel charge trailing (e.g. top left corner of image), but they morph back out of Grade 7 due to serial charge loss (e.g. top of image around column 300). The degree of grade morphing is a complicated function of event energy and the original event grade. The mirroring effect seen in the BI images is due to the alternating direction of serial readout in the four amplifiers.

## 3 The Model

### 3.1 The Primary Model

Our CTI model is phenomenological rather than physical, opting to characterize the effects of CTI in the data rather than to model directly the spatial distribution and time constants of the trap population (see , , and  for examples of such physical models). It consists of a set of amplifier-dependent piecewise-linear equations that relate the quantity of charge lost from a pixel to the original quantity of charge in that pixel. In turn, the amount of charge trailed into the adjacent pixel was determined empirically to be a piecewise-linear function of the charge loss. Trailing is not modeled as an amplifier-dependent relation, rather data from all amplifiers can be combined to derive the fit parameters.

We measure charge loss at a given energy by fitting a line to a plot of the central pixel value (in DN) for only Grade 0 (single-pixel) events as a function of event position on the chip (e.g. row number, called “CHIPY,” for FI devices). All data are fit as an ensemble of points, using sigma-clipping to remove outliers. This is illustrated in Figure 6. We use the Event Browser software  for display and fitting. Here the event data are binned into an image for display only; the actual fitting is done on the events themselves, not on binned data.

We show the Event Browser window here because it illustrates details of the fitting process. The $`+`$’s show the input to the fitting tool: approximate endpoints of the line and a rough guess at the linewidth (marked at CHIPY $``$450). The set of dashed lines at the edges of the spectral line show iterations of the sigma-clipping. The final fit is shown as a solid line.

Charge trailing is measured similarly: for FI devices, we plot the “top center” pixel of the $`3\times 3`$ neighborhood for ACIS Grade 0 and Grade 64 (upward singly-split) events as a function of CHIPY, then perform a sigma-clipped linear fit to the ensemble of events. This is illustrated in Figure 7, where again we show the Event Browser window to illustrate the input to the fitting routine ($`+`$’s), the sigma-clipping (dashed lines), and the final fit (solid line).

Our linear characterization of the data leads to a linear model of the charge lost ($``$) and trailed ($`𝒯`$) for each pixel, i.e.

$$=N_tL(P)$$
(1)
$$𝒯=N_tT(L)$$
(2)
where $`N_t`$ is the number of transfers that the charge packet must undergo to reach the readout node, $`P`$ is the charge in a single pixel, $`L(P)`$ is the charge lost per pixel transfer for an isolated charge packet, and $`T(L)`$ is the charge trailed per pixel transfer for an isolated charge packet. Obviously the physical process of CTI occurs one transfer at a time and the correct model would recalculate the pixel’s charge for each transfer separately (e.g. via a loop in the code), but that is computationally expensive so we do not attempt such a solution.

Accurate modeling of CTI requires accurate knowledge of the shape of $`L(P)`$ and $`T(L)`$, especially for small charge packets (low energies) where we have very little calibration information. We have found that the grades and energies of corrected events depend on the exact shape of $`L(P)`$ and $`T(L)`$ at low energies. For example, a single linear fit to $`L(P)`$ yields a non-zero intercept, implying that there can be substantial charge loss in pixels originally containing very little charge – this is not supported by the data. This single-line model causes the corrector to boost the charge in low-valued pixels, often causing them to exceed the split threshold (13 DN for both BI and FI devices). This disrupts the grade distribution of a set of monoenergetic events and causes the split events in that set to exhibit higher energies than single-pixel events. We have determined empirically that $`L(P)`$ and $`T(L)`$ must pass through the origin to yield realistic grades and energies for monochromatic input photons.

Thus we model $`L(P)`$ and $`T(L)`$ as piecewise linear, continuous curves that are required to pass through the origin (see Figure 8). Three lines are needed for $`L(P)`$ and two for $`T(L)`$. The functions’ values (for both parallel and serial CTI) are stored as lookup tables that are accessed for each pixel of each event. We chose a set of lines as the functional form for these relations because they adequately fit the data, they keep the model simple, and because we can restrict the coefficients to ensure zero charge loss or trailing when there is zero charge present.

Our original model for $``$110C FI data  employed a power-law fit, which worked well for all measurable energies up to $``$6 keV, but underfit the Au L$`\alpha `$ line at 9.7 keV. Application to the first Galactic Center observation made by ACIS showed that this underfitting might be affecting observations of astrophysical iron lines (Y. Maeda, private communication), so the power-law model was abandoned in favor of the piecewise linear model. The new model improves the $``$110C corrector performance at high energies without compromising the low-energy accuracy. For the more modest FI CTI effects seen at $``$120C and for BI devices, the power law fit is adequate  but we choose to use the piecewise linear model for all temperatures and devices for consistency.

For the BI devices, we must model parallel and serial charge loss and trailing independently to match the data. There are charge traps in the imaging array, framestore array, and serial register. Since the FI devices were damaged by radiation that penetrates only a few microns into the device , their (shielded) framestore arrays and serial registers are protected. So, although the CTI is more pronounced on FI devices, it is easier to model than in the BI devices.

The primary FI CTI model consists of a piecewise linear fit to the measurements of parallel charge loss per pixel transfer as a function of line center (in DN units), measured using only the central pixel of Grade 0 events and on a per-amplifier basis. The components cover the ranges 0–100 DN (0–0.4 keV), 100–550 DN (0.4–2.2 keV) and 550–2500 DN (2.2–10 keV) (see the left panel of Figure 8). The mid-energy line is necessary to accomodate a known turnover around these energies (known from calibrations with the soft supernova remnant E0102-72.3). The break is wherever the mid- and high-energy lines meet and is calculated in the code. The low-energy line takes over at 100 DN, a value determined roughly by looking at the energies of corrected split events in the S3 (BI) device. It forces the model through the origin. As mentioned above, this was necessary to avoid anomalous energies for low-energy split events in S3 – in effect, it forces the model to assume zero charge loss for zero-energy events.

An example of the charge trailing model is shown in the right panel of Figure 8. For trailing, only two lines were necessary to fit the data; the low-energy line is forced through the origin as in the loss model. The break between the lines is taken as the loss value where the third line takes over in the loss model (averaged over the four amplifiers). For the example in Figure 8, the two lines in the fit have very similar slopes. Note that the Au L$`\alpha `$ points are omitted from the trailing model. Due to the large charge cloud size and confusion with background events, our method for measuring charge trailing proved unreliable for Au L$`\alpha `$ .

The (very small) error bars represent only measurement uncertainty (estimated by fitting each line 3 times); we have not attempted to quantify systematics. The high-energy charge loss fit was made using errors in both line centers and charge loss and used all points shown on the plot above 500 DN. The low-energy charge loss fit only includes the Al K$`\alpha `$ and Mn/Fe L-complex points.

A similar approach is taken for the BI devices. First the serial charge loss and trailing model coefficients are determined for each amplifier, also using a 3-component piecewise linear fit for loss and a 2-component piecewise linear fit for trailing. Then the CTI corrector (described below) is run on the data, with coefficients for parallel CTI set to zero. The resulting serial-corrected data are used to measure the coefficients for parallel CTI. Due to non-linearities in the serial CTI, columns near the edges of each amplifier are excluded when computing the linear parts of both the serial and parallel models.

As a check on the functional form of the model at low energies, we measured the low-energy parallel charge loss using the soft supernova remnant E0102-72.3 (“E0102”), a calibration target for Chandra. Figure 9 shows the charge loss model for I3, Amplifier 2, at a focal plane temperature of $``$110C, with the E0102 charge loss values overplotted for the three soft lines measurable with these data.

The E0102 points were obtained from two datasets, Observation ID 49 (source placed at CCD row number $`100`$) and Observation ID 48 (source placed at CCD row number $`500`$). The charge loss per pixel transfer was calculated using only these data, then corrected for the fact that the full-chip charge loss slope is more shallow than that calculated when only using the bottom half of the chip by executing a similar half-chip analysis on ECS data, then determining the correction factor necessary to recover the full-chip slope.

The E0102 data were not used to modify the model. Their agreement with the low-energy model components confirms that the 3-line fit chosen to parameterize the charge loss is reasonable down to $`0.5`$ keV at least. This is adequate for characterizing FI devices but we have no confirmation of our model in the 0.2–0.5 keV range accessible to BI devices.

### 3.2 Shielding Within an Event

This CTI model is “local” – prediction of the amount of charge lost by a particular event does not consider the possibility that another event falling in the same columns but closer to the amplifier readout will have filled some portion of the traps, reducing the charge loss from its nominal value (see  and  for descriptions of such “precursor” models). We do however consider shielding effects that may occur within the nine pixels of a single event. Physical principles suggest that pixels which are not isolated experience loss and trailing that are influenced by their neighbor pixels – what we call the “self-shielding” effect. We hypothesize that the loss mechanism essentially acts only on the portion of a pixel’s charge which exceeds the charge in the previous pixel; similarly the trailing mechanism essentially acts only on the portion of a pixel’s charge which exceeds the charge in the following pixel. Thus a given pixel experiences either loss or trailing, not both – loss if its predecessor has less charge than the pixel in question and trailing if its predecessor has more charge.

The physical process of self-shielding is of course governed by the relationship between the physical volume of silicon occupied by one pixel’s charge packet compared with the physical volume of silicon occupied by its predecessor or successor charge packet. A physically-based model of self-shielding would involve knowing the details of the three-dimensional charge distribution . For our phenomenological model, we instead take a very simple approach to modeling this effect.

For a pixel with charge $`P`$ (“root”) whose predecessor (“lead”) has a smaller charge $`P_{\mathrm{lead}}`$, we compute:

$$_{\mathrm{root}}=N_tL(P)<P$$
(3)
$$_{\mathrm{lead}}=N_tL(P_{\mathrm{lead}})<P_{\mathrm{lead}}$$
(4)
$$_{\mathrm{eff}}=(_{\mathrm{root}}_{\mathrm{lead}})>0$$
(5)
The effective charge loss $`_{\mathrm{eff}}`$ is then subtracted from $`P`$.

For a pixel with charge $`P`$ (“root”) whose predecessor (“lead”) has a larger charge $`P_{\mathrm{lead}}`$, we compute:

$$𝒯_{\mathrm{root}}=N_tT(P)<P$$
(6)
$$𝒯_{\mathrm{lead}}=N_tT(P_{\mathrm{lead}})<P_{\mathrm{lead}}$$
(7)
$$𝒯_{\mathrm{eff}}=(𝒯_{\mathrm{lead}}𝒯_{\mathrm{root}})>0$$
(8)
The effective charge trailed $`𝒯_{\mathrm{eff}}`$ is then added to $`P`$. In this model, events which are almost evenly split between two pixels in the direction of transfer suffer CTI effects almost identical to those experienced by a Grade 0 event with half the energy of the split event.

### 3.3 Deviations from Linearity

Although ACIS FI devices do not exhibit serial CTI, they do show a small charge loss effect which is proportional to the distance between the interaction site and the left or right edge of the CCD (as opposed to the distance to the amplifier readout). This is believed to be caused by the slight drop in clocking voltages necessarily suffered between the middle of the CCD and the supply leads at the chip edges (M. Bautz, private communication). CTI scales the effect so that these gain variations are row- and energy-dependent, with the bottom rows of the chip showing virtually no variation and the top rows varying by up to $``$100 DN at high energies ($``$6 keV).

ACIS BI devices also show low-spatial-frequency deviations from linearity in their charge loss. As shown above, by combining large numbers of ECS observations, column-to-column gain variations can be seen. (Row-to-row variations are not seen in either FI or BI devices.) These effects are also energy-dependent, with high-energy calibration lines in the ECS exhibiting larger-amplitude variations than low-energy lines.

These nonlinearities result in broadened spectral lines for both FI and BI devices even after correction for linear charge loss and trailing. Thus we have developed the concept of a “deviation map” to make moderate gain adjustments to each event that help to regularize the CCD response across the device and narrow the linewidths. Our techniques are reminiscent of flat-fielding techniques used to remove instrumental non-uniformities in visual CCD data.

We model the deviations via a 2-dimensional map ($`D(x_c,y_c,E)`$) with an energy-dependent amplitude:

$$D(x_c,y_c,E)=D_0(x_c,y_c)(1+gE)$$
(9)
where $`x_c`$ and $`y_c`$ are the position coordinates of the event in units of CCD pixels, $`E`$ = event amplitude in DN, $`g`$ = deviation amplitude, and $`D_0`$ is the two-dimensional map of deviations extrapolated to 0 DN. The value of this deviation is calculated for each event based on its $`(x_c,y_c)`$ location and applied as an additive adjustment to the central pixel value. Figure 10 shows images of $`D_0`$ for I3 and S3 at $``$120C. Cuts through the I3 map, given in Figure 11, show that the amplitude of the deviations increases substantially with $`y_c`$ on FI devices, as noted above.

The zero-energy deviation map $`D_0`$ is made for each amplifier in a series of smoothing and averaging steps using constituent maps made in the three strongest lines in the ECS: Al K$`\alpha `$ , Ti K$`\alpha `$ , and Mn K$`\alpha `$ . First the ECS data are CTI-corrected using the linear model described above. The resulting event list is filtered spectrally to keep only events in one of these strong lines. We then make an image of these event amplitudes in DN space, binned by 1 pixel in $`x_c`$ and by 16 pixels in $`y_c`$, with intensity computed as a sigma-clipped mean. This process is repeated for the other two lines, resulting in three images of the deviations covering a substantial range in energy.

For each energy, we then perform a trending analysis of gaussian fits to the spectral line, noting the mean and standard deviation (width) of the spectral line as a function of row number. We perform a linear fit to these trended properties to determine the line energy and width at row zero, where the values are unaffected by CTI. Then we plot the row-zero standard deviations as a function of line energy (in DN) and fit a line, noting its slope (S) and intercept ($`\sigma _0`$). This gives the equation

$$\sigma (E)=\sigma _0+SE=\sigma _0(1+\frac{S}{\sigma _0}E)$$
(10)
which, by analogy to the definition of the deviation map above, yields

$$g=\frac{S}{\sigma _0}.$$
(11)
We obtain an estimate of $`S`$ and $`\sigma _0`$ for each amplifier; these are averaged then divided to obtain $`g`$.

The maps at Al K$`\alpha `$ , Ti K$`\alpha `$ , and Mn K$`\alpha `$ are then smoothed column-wise to reduce the effects of shot noise. Next we generate three estimates of $`D_0`$ using the smoothed maps and average those three estimates to obtain our final $`D_0`$ for that amp. Lastly, the maps for each amplifier are incorporated into an array to make one deviation map for the entire CCD.

This deviation map is used to adjust the amplitude of an event. We compute $`E`$ in Equation (9) above by summing the values of pixels in the $`3\times 3`$ event neighborhood that are above the split threshold (13 DN). The adjustment is then computed from Equation (9), using the event’s $`(x_c,y_c)`$ location to access the correct element of $`D_0`$. The adjustment is made by adding the result to the central pixel of the event. This preserves the event grade; we determined by experimentation that this technique produces the best reconstructed event.

## 4 Simulator Results

The Monte Carlo simulator’s ability to reproduce the effects of CTI is illustrated in the following figures. All data and simulations are for Al K$`\alpha `$ (1.486 keV) with ASCA g02346 grade filtering. In both figures, the left panels show the data and the right give the matching simulations.

The upper panels of Figure 12 show images of the median event amplitude events on an FI device (I3 at $``$120C). Both datasets contain $`10^5`$ events and are binned at $`8\times 8`$ pixels. With this binning, column-to-column variations are not visible, but variations with lower spatial frequency are clear and are reproduced with good fidelity by the simulator.

When modeling CTI for the event simulator, we add noise to the pixel adjustments to model the stochastic nature of traps. We find that a simple Poisson noise model on the number of electrons lost is not sufficient for FI devices, so we introduce an energy-dependent tuning parameter that is multiplied by the event’s row number to account for the additional line broadening with increasing distance from the readout node that is seen in the ECS data. The effects of this are apparent in the lower panels of Figure 12. Row-dependent non-linearities in the CTI give the data a slight S shape in this plot; this feature is also seen in the simulation.

The data here show more scatter than the simulation; this is due in part to the background in ECS data and to Al K$`\alpha `$ events that occurred in the framestore region of the CCD. These framestore events do not suffer CTI and appear as a faint horizontal band in the lower left panel of Figure 12.

The upper panels of Figure 13 show the same Al K$`\alpha `$ map as Figure 12, now for the BI device S3. As the figure clearly shows, the four amplifiers in the S3 device have slightly different responses. This is reproduced well by the simulation, as are the CTI-induced spatial variations. The column-wise variations in the Al K$`\alpha `$ event amplitudes are shown in the lower panels of Figure 13. Again the simulator does a reasonable job of reproducing these variations. Both the imaging and the framestore sections of BI CCDs suffer CTI, so no horizontal band of CTI-free events is seen in the lower left panel of Figure 13.

## 5 The Corrector

As mentioned above, we have used our Monte Carlo modeling of CTI effects to build a CTI corrector. This program attempts to remove those effects from individual ACIS events, taking a forward-modeling approach. For each observed event $`O`$ we first hypothesize the nine pixel values of a corresponding clean event $`C`$, i.e. the event that would have been obtained if no CTI effects were present. This clean event is passed through the CTI model, producing a model event $`M`$ which estimates what we should have observed if our clean event hypothesis is correct. The hypothesis (the nine pixel values of $`C`$) is then adjusted using the differences between the observed event and the model event in a simple way:

$$C_{new}=C_{old}+(OM).$$
(12)
This iterative process continues until the model and observed events differ by less than 0.1 DN (for each of the 9 pixels). When convergence is obtained, C is returned as the CTI-corrected event.

Tests show that when a limit of 15 iterations is imposed within the corrector a trivial number of events fail to meet the 0.1 DN convergence criterion; this limit is imposed as a safety in the code and the number of iterations needed for convergence (typically 3–5) is recorded for every simulated event. The pixel values and thus the event amplitudes produced by the corrector are real numbers.

Because the corrector operates on event islands (3$`\times `$3 pixel arrays) rather than on event amplitudes, it has the potential to improve the grades in an event list, leading to a different grade distribution. Thus the CCD quantum efficiency appropriate for grade-filtered corrected data will differ from that appropriate for similarly grade-filtered uncorrected data. A technique for determining this QE adjustment is presented in Section 7.

The CTI corrector consists of a set of IDL  programs and device-dependent parameter files that instantiate the model for each amplifier on each CCD, for a given epoch of the Chandra mission. Obviously observations must be configured to telemeter the $`3\times 3`$ pixel event island in order for the corrector to modify that island. The corrector operates on the “Level 1” event list, one of the standard data products, and requires no additional information regarding the observation. The parameter tuning presumes that the CCDs were operated in their standard full-frame mode, with a frame readout time of 3.24 seconds. Since the trap time constants interact with this readout time to affect the average number of filled traps, observations made in any of the “sub-array” modes (including Continuous Clocking Mode) will exhibit different CTI and cannot be corrected using this full-frame tuning.

The code is available and may be of interest to other groups concerned with the results of radiation damage on X-ray and visual CCD detector systems and to Chandra/ACIS users. We have recently completed tuning our CTI corrector for all chips in the ACIS-I array and chips S2 and S3 in the ACIS-S array, for focal plane temperatures of $``$110C and $``$120C. We provide the parameter files used by the corrector code on our webpage, so any Chandra/ACIS user with data on one or more of the allowed CCDs may experiment with the PSU CTI corrector.

Response matrices have been generated for these device/temperature combinations as well, using the CTI model as part of the Monte Carlo CCD simulator to generate simulated events at a large number of monochromatic energies. The simulator applies CTI to the events, then the corrector is used to remove the non-random effects of CTI. The resulting events are combined into a (now polychromatic) instrumental spectrum, which represents a single row of the response matrix. An example was shown in T01. These response matrices are also available to Chandra/ACIS users.

## 6 Corrector Results

### 6.1 Spectral Resolution

The most fundamental result of CTI correction is the improvement in spectral resolution. This is illustrated in Figure 14 for an FI device at two ECS energies. This figure also illustrates the major improvement in FI spectral resolution made by reducing the ACIS focal plane temperature from $``$110C to $``$120C. The CTI corrector improves the resolution at both focal plane temperatures and all energies for the uniform illumination of the ECS. Celestial sources, which tend to be more pointlike and thus suffer less from the variations across the CCD than ECS data, will show less dramatic improvement.

Figure 15 shows the ECS spectral resolution improvement with CTI correction on a BI device, now as a function of energy rather than position, as the BI devices do not exhibit position-dependent resolution. Only $``$110C data are shown here; $``$120C results are similar. Again the spectral resolution is improved at all ECS energies. Below 1.5 keV, though, the device’s intrinsic resolving ability dominates the CTI effects. Here, CTI correction will not markedly improve the spectral resolution. As for FI devices, the spectral resolution improvement for more pointlike celestial sources will be less dramatic.

Comparing Figures 14 and 15 shows that the best ACIS spectral resolution is achieved by placing targets on the lower half of ACIS FI devices. In these regions, the FI resolution surpasses the BI resolution. Without CTI correction, only the lowest quarter of the FI devices surpassed the BI resolution. This improvement is particularly important for targets with extended emission or point sources distributed across the field, such as star formation regions  and the Chandra Deep Field observations of the X-ray background .

### 6.2 Spectral and Spatial Comparisons

To give more detail on the effects that the CTI corrector has on ECS data, figures identical to those in Section 2 are presented below.

Figure 16 shows CTI-corrected ECS spectra. For ease of comparison, in this one case we reproduce Figure 1 in the top two panels and show the CTI-corrected spectra below. The spectral resolution is improved for both FI (left) and BI (right) spectra.

The performance of the corrector in removing the observed spatial variation of event amplitude in FI ECS data is shown in Figure 17. The corrector has largely removed the gain changes across the device; compare to Figure 2. Since the spectral broadening with row number is primarily the result of a random process involving one of the trap species , it cannot be completely suppressed in our reconstructed events.

Figure 18 illustrates the remaining non-uniformities after the CTI corrector, with its deviation map, is used to remove the column-dependent gain variations caused by CTI. Comparing this to Figure 3, it is clear that accounting for column-to-column gain variations improves the quality of the data.

Now turning to a BI device, Figure 19 shows the ECS data from S3, Amplifier 0, at $``$120C, after CTI correction. The spatial effects of CTI are much more subtle on BI devices than on FI devices; the most obvious improvement compared to Figure 4 is in the lower left panel, where the column-to-column gain variations have been largely suppressed by the CTI corrector.

Finally, Figure 20 shows the grade-specific 6 keV images from Figure 5 after CTI correction. The position dependence of undesirable grades (1, 5, and 7) is much reduced.

### 6.3 Event Grades

As mentioned above, grade morphing due to CTI is a substantial concern for FI observations early in the mission, due to the higher focal plane temperature of $``$110C. This problem was largely solved by lowering the focal plane temperature, but it is worth noting here the effects at $``$110C for archival researchers.

Figure 21 shows the effect of grade-morphing on the $``$110C FI detector quantum efficiency, for the Mn K$`\alpha `$ and Mn K$`\beta `$ lines combined. At this energy, many events are lost forever because they morphed into grades that were not telemetered. For low-energy events (e.g. Al K) even at $``$110C, the charge is more spatially concentrated, leading to more single-pixel events. CTI causes these to morph into grades that are still captured in the standard grade filtering, so the QE at the top of the device is not affected by CTI.

For this high-energy dataset, keeping all grades, there are 12% fewer events at the top of the chip compared to the bottom due to grade morphing and subsequent on-orbit grade rejection. For standard grade filtering with standard processing, this loss is augmented to 20% at the top of the chip. This is because grade morphing caused events to migrate out of acceptable grades. By applying the standard grade filter after applying the CTI corrector, though, there are only 13% fewer events at the top of the chip compared to the bottom. Since this is similar to the curve for all grades, we can conclude that we have recovered virtually all useful events available on the ground.

For the BI devices, grade morphing is an issue even at $``$120C. The success of the CTI corrector in regularizing the grade distribution for the top half of a BI device (where grade morphing is most pronounced), using on-orbit ECS data on the S3 chip at $``$120C, is shown in Table 2, where original and corrected ASCA-like grade distributions are listed for the Al, Ti, and Mn calibration lines.

Some events are not recovered into their original grades because the corrector is unable to recover charge that has been eroded by CTI down below the split threshold. At that level, the charge is indistinguishable from noise and the corrector is purposely not allowed to include such pixels in its reconstruction.

A good test of grade recovery is to compare the number of Grade 2 events to the sum of Grade 3 and Grade 4 events – for devices with no CTI, these quantities should be roughly the same. Applying this diagnostic to Table 2 shows that the corrector is quite useful in recovering appropriate grades.

Grade-dependent BI spectra are shown in Figure 22: the left-hand panel gives the results for standard processing and the right-hand panel those for the CTI corrector. Between 1 and 7 keV, it appears that including ASCA grades 1, 5, and 7 might improve the signal-to-noise ratio in CTI-corrected data, as the ECS spectral lines dominate over the background spectrum. In general, though, CTI correction has restored the original intent of ASCA g02346 grade filtering, to suppress the particle background and improve signal-to-noise across the full spatial and spectral range of the device.

### 6.4 Gain linearity

Once the CTI corrector is applied to adjust the amplitude and grade of each ACIS event, the DN-to-eV conversion is recomputed using the ECS lines. A linear fit is made to obtain the gain (slope) and offset (y-intercept) for each amplifier. We can then compute the difference between the measured ECS line energies and the values predicted by this linear model; these residuals are a measure of the non-linearity of the amplifier after correcting for CTI in the data. Figure 23 shows the ECS line energy residuals for all amps on I3 (left) and S3 (right) at $``$120C.

These residuals show that a linear gain model is reasonable for CTI-corrected ECS data on FI devices at all relevant energies; deviations are seldom larger than $`\pm 10`$ eV. For the BI device S3, the deviations are larger at energies below 2 keV. The ECS line complex at $``$680 eV is faint and more information is needed to search for gain non-linearities in the important BI spectral range 0.2-0.5 keV. Unfortunately virtually no calibration information at such low energies is available at this time.

Figures 24 and 25 show the position dependence of gain variations remaining after CTI correction. Three median energy maps are given, for the three main ECS lines. These greyscale maps are displayed such that black is $`+1\sigma `$ from the mean energy and white is $`1\sigma `$ from the mean energy. The spectral lines are shown in the right panels.

For the FI device I3 at $``$120C, Figure 24 shows that there is some remaining structure following amplifier boundaries. There is no clear correlation between energies, however; regions that are light at one energy are neutral or dark at another. Thus the algorithm, although imperfect, appears to be removing gain variations in a reasonable way.

Figure 25 shows that the same is true for the BI device S3. Non-random gain deviations exist, but they are not correlated across energies. It might be possible to generate a second-stage deviation map to try to remove these residual features, but our attempts using this dataset (roughly one year of ECS observations) have been unsuccessful.

## 7 Quantum Efficiency Adjustments

As mentioned above, grade adjustments by the CTI corrector combined with subsequent grade filtering result in CCD quantum efficiency improvement. This change in QE must be incorporated into the ACIS calibration products so that source fluxes can be estimated properly. Since the CTI-induced grade morphing changes with position on the CCD, the most direct way to account for these QE changes is via the Quantum Efficiency Uniformity (QEU) files, one component of the Auxiliary Response File (ARF) used in spectral fitting. The ARF incorporates the vignetting function of the telescope, the positionally-averaged QE of the CCD, and position-dependent deviations from that average QE via the QEU files.

Two approaches to adjusting the QEU files for CTI-corrected data are necessary: one for FI devices that were effectively CTI-free before launch, and a second for BI devices that have always suffered measurable CTI. These are described below.

### 7.1 Estimating QE Loss with Respect to an Undamaged Device

This approach is used for FI devices. We estimate how the QE of undamaged devices changes when the damage is applied and the resulting events are corrected for CTI (for example by the simulation experiments shown in Figure 26). Such a loss is of course a function of both energy and position on the CCD. Then the QEU files that were constructed for undamaged devices are adjusted downward to account for unrecoverable CTI losses. For each energy $`E_i`$,

$$QE(E_i,X_j,Y_k)=QE_0(E_i)QEU(E_i,X_j,Y_k)qe\mathrm{\_}loss(E_i,X_j,Y_k)$$

In this formulation, the $`QE_0`$ & $`QEU`$ terms are the mean QE and QE uniformity products produced by standard processing for undamaged devices, and the $`qe\mathrm{\_}loss`$ term expresses the QE lost by the damage plus correction process. For ACIS FI devices the $`qe\mathrm{\_}loss`$ term is not actually a function of X (CCD column number). The pre-launch FI QE is very spatially uniform below $``$8 keV, consistent with unity to within measurement errors . Thus for ACIS FI devices the $`QEU`$ term above is set to unity.

### 7.2 Estimating QE Recovery with Respect to a Damaged Device

This approach is used for BI devices. We estimate the QE recovery of the CTI corrector via simulation/correction experiments, as shown in Figure 27. This recovery is also a function of both energy and position on the CCD. Then QEU products that were constructed for damaged devices are adjusted upward to account for the increased QE after CTI correction. So for each energy $`E_i`$,

$$QE(E_i,X_j,Y_k)=QE_0(E_i)QEU(E_i,X_j,Y_k)qe\mathrm{\_}recovery(E_i,X_j,Y_k)$$

In this formulation, the $`QE_0`$ and $`QEU`$ terms are the mean QE and QE uniformity products generated by standard processing for damaged devices, and the $`qe\mathrm{\_}recovery`$ term expresses the QE recovered by the correction process.

## 8 Model Performance on Celestial Data

There is an uncertainty introduced by using on-orbit calibration data to tune this or any other model: we must assume that the CTI measured by the on-orbit calibration source is representative of that present in celestial data. This assumption could be incorrect, even for celestial data obtained at nearly the same time as the calibration data, due to differences in the count rates (photons plus particle background) between these two configurations. It is important to remember, however, that most ACIS events are due to particle interactions and are not even telemetered, so the event list supplied to the user is not representative of all the potentially trap-filling interactions that occurred. Since most of these particle events are from cosmic rays that enter from all angles and are not absorbed by the Observatory structure, it is reasonable to assume that calibration and celestial data have comparable particle fluxes, thus a comparable trap population.

We have used the instrumental Au L$`\alpha `$ (9.7 keV) and Ni K$`\alpha `$ (7.5 keV) lines as rough diagnostics of the corrector’s performance on celestial data. These spectral features are due to fluorescence from hardware components. The lines are faint, but in long exposures it is possible to measure their gain changes with position on the chip just like we measure lines in the calibration data. The tests show that the corrector works fairly well even at these high energies, reducing the charge loss per pixel transfer by a factor of $``$10.

Another implication of tuning the corrector to the ECS data is that we have only minimal calibration information below the Al K$`\alpha `$ line at 1.486 keV and no information below the Fe-Mn L complex at $``$680 eV. We have extrapolated to lower energies using fits of sensible functional forms, but it is likely that our assumptions break down for small amounts of charge. Understanding the behavior of the CCDs between 0.2 and 1.5 keV is important for obtaining accurate hydrogen column densities, comparing ACIS results to findings from other X-ray missions, and inferring valid astrophysics for the large number of soft sources Chandra is examining.

We can use celestial data to test and improve our calibration of the CTI at low energies. The best available target is the supernova remnant E0102-72.3 (“E0102”) in the Large Magellanic Cloud, used as a soft calibration source for Chandra. We used this source in Figure 9 above to confirm that the piecewise-linear CTI model was reasonable at low energies.

Figure 28 shows the difference between measured and true line energies from ECS data and the E0102 soft spectral lines. Energies for the E0102 lines were provided by R. Edgar (private communication) and can be inferred from Rasmussen et al. . This shows that the model can recover reasonable line energies for both FI and BI devices down to $``$0.5 keV. The I3 plot includes E0102 points obtained from combining Observation ID’s 136, 140, 439, 440, 1314, 1315, 1316, and 1317. The S3 plot includes E0102 points obtained from combining Observation ID’s 141, 1308, 1311, 1530, 1531, and 1702.

Since E0102 has been observed at several off-axis positions on I3, we can compare its standard and CTI-corrected spectra as a function of the position-dependent spectral resolution. Figure 29 shows these spectra, with the spectral resolution increasing with off-axis angle (the I-array aimpoint is at the top of the I3 chip, where CTI-induced spectral resolution degradation is at its worst). All spectra used ASCA g02346 grade filtering, a binsize of 10 eV, and were smoothed with a Gaussian of width 1 bin.

Near the aimpoint, the CTI corrector changes the low-energy gain substantially and improves the spectral resolution. At 4 arcmin off-axis (about half-way up the I3 chip) the spectral lines are noticeably narrower. The CTI corrector applies a more modest gain change and still improves the spectral resolution compared to standard processing. At 7 arcmin off-axis (only 100 rows from the bottom of the chip, where CTI is minimal), standard and CTI-corrected spectra are very similar in gain and spectral resolution. The line energies between the three observations are comparable in CTI-corrected data.

## 9 Summary

We have developed a phenomenological Monte Carlo model of CTI in both front- and back-illuminated ACIS CCDs and used it to devise a CTI corrector for ACIS observations of celestial X-ray sources. The PSU CTI corrector has been tuned for the ACIS CCDs most often used for imaging: I0–I3, S2, and S3. This tuning uses a composite dataset of observations of the ACIS External Calibration Source. Both major focal plane temperatures ($``$110C, September 1999 – 29 January 2000, and $``$120C, 29 January 2000 – present) are calibrated.

The CCD simulator and CTI corrector IDL source code are provided to the community<sup>1</sup><sup>1</sup>1http://www.astro.psu.edu/users/townsley/cti/. Also available are the parameter files that instantiate the simulator and corrector tuning for each device, the modified QEU files, and a suite of response matrices for use with CTI-corrected data.

The corrector works on a Level 1 event list from ACIS observations using the full-frame readout mode. The $`3\times 3`$ or $`5\times 5`$ pixel event islands must be included in order for the corrector to modify the event grade. The user is then required to perform the other filtering steps (grade filtering, applying Good Time Intervals, etc.) to recover the CTI-corrected Level 2 event list.

Refinements of the model since the original report of this work (in early 2000)  have led to several improvements in CTI-corrected data:

* Modeling the charge loss per pixel transfer as three lines improves the fidelity of corrected event amplitudes at high energies without compromising the low-energy results.
* Removing column-to-column gain variations (to the degree possible with this model) improves the energy resolution on both FI and BI CCDs.
* Including charge trailing and shielding within a single event makes it possible to compensate for CTI-induced grade morphing at all event locations.
* Using a more extensive set of calibration observations of the soft SNR E0102-72.3, we have improved the model tuning in the spectral range 0.5–1.5 keV.

Since this corrector does not account for precursor events, we cannot totally eliminate the row-dependent energy resolution of the FI devices. Thus position-dependent response matrices are still necessary for accurate FI spectral fitting. By regularizing the gain and grade distributions across these CCDs, though, we have reduced the number of matrices needed for typical analysis tasks to fewer than 10. For BI devices, the improvement is even better; since S3 does not exhibit position-dependent energy resolution, a single S3 response matrix is usually adequate for CTI-corrected data. This is a great simplification for users, especially those working on extended targets such as supernova remnants or clusters of galaxies.

Lowering the ACIS focal plane temperature to $``$120C achieved dramatic reduction in the effects of CTI on FI devices. This reduced the FI spectral resolution degradation and its row-dependence, moderated the gain changes across the device, and reduced the charge trailing, resulting in much less high-energy quantum efficiency loss at the top of FI devices. Since a large fraction of ACIS data were obtained at $``$120C, we have concentrated on those data in this paper. It is important to remember, however, that the early days of the mission included observations of some very exciting targets and those data are now publicly available. So for the sake of archival research, we have invested substantial effort in tuning the CTI corrector at $``$110C.

The CTI corrector has been used for a variety of Chandra investigations, including:

* Resolving the X-ray background emission in the Chandra Deep Field North, e.g. , ,
* Star formation regions, e.g. ,
* Supernova remnants, e.g.
* The starburst galaxy M82
* A search for low-luminosity active galactic nuclei .

All parameter files, source code, and instructions comprising the PSU CTI corrector and the PSU CCD simulator are available on the Web<sup>2</sup><sup>2</sup>2http://www.astro.psu.edu/users/townsley/cti/. We encourage the use of these tools and welcome feedback.

Other software mitigation techniques for CTI are under development at the Chandra X-ray Center and at MIT, with support by Chandra Project Science at Marshall Space Flight Center and the rest of the Chandra team. The MIT/ACIS group also is working to devise hardware mitigation techniques that may further reduce the effects of the radiation damage. Prigozhin et al. ,  are developing a more physical model of the charge traps that is likely to replace this phenomenological model. We are confident that these efforts will eventually lead to a more robust CTI correction scheme for ACIS data.

Financial support for this effort was provided by NASA contract NAS8-38252 to G. P. G., the ACIS Principal Investigator. This work made use of the NASA Astrophysics Data System. We appreciate the time and gracious response of the anonymous referee. The External Calibration Source data were provided by Allyn Tennant ($``$110C) and the CXC ($``$120C). We are very grateful for their efforts to organize and compile this large body of information. The extensive data analysis necessary to tune the CTI corrector and simulator for each CCD and temperature would not have been practical without the help of many members of the PSU/ACIS team. We thank Xinyu Dai, Sarah Gallagher, Konstantin Getman, Varsha Gupta, Ann Hornschemeier, Bulent Kiziltan, Karen Lewis, Sangwook Park, Divas Sanwal, Michael Sipior, and Yohko Tsuboi for their hard labor and George Chartas and Yoshitomo Maeda for their ideas. We thank our MIT/ACIS colleagues, especially Catherine Grant, Mark Bautz, and Gregory Prigozhin, for many helpful discussions and suggestions regarding CCD modeling and for their lucid explanations of the device physics of radiation damage and charge traps. We thank the CXC for providing a forum for discussion of CTI and other ACIS calibration issues and we especially thank Paul Plucinsky for his help with the analysis of E0102-72.3 data, for testing our data products, and for reminding us, now many years ago, that CTI was a force to reckon with.


format: text

text:
# Specific heat of the Coulomb glass

## 1 Introduction

The Coulomb glass, a prominent model of disordered systems of localised particles with long-range interaction, has been extensively studied for almost three decades. It plays an important role as a semi-classical model for dilute impurity bands when quantum interference as well as thermal excitation to the conduction band or from the valence band can be neglected. Various aspects of this model have been investigated since the pioneering works by Pollak (1970), Srinivasan (1970), Efros and Shklovskii (1975), and Kurosawa and Sugimoto (1975); for related reviews see Shklovskii and Efros (1984), Efros and Shklovskii (1985), Pollak and Ortuño (1985), and Ortuño et al. (2001). Nevertheless, this area of research has not lost its attraction, as is shown, for example, by recent studies of aging effects in electric conduction (Vaknin et al. 2000), of the possibility of a phase transition (Grannan and Yu 1993, Vojta and Schreiber 1994, Díaz-Sánchez et al. 2000b), and by combined investigations of dc conductivity and tunnel conductance (Sandow et al. 2001).

An important and characteristic feature of the Coulomb glass is the many-valley structure of its configuration space, not only concerning single-particle hops, but also with respect to some classes of more complex excitations. In this respect it resembles the Ising spin glass and other disordered systems, see e.g. Mézard et al. (1987). This configuration-space structure is related to the existence of a huge number of so-called pseudo-ground states, as was observed in the first numerical Coulomb glass simulations (Baranovskii et al. 1976). The spectrum of these pseudo-ground states has recently been investigated in detail by Kogan (1998). Moreover, the related importance of many-particle excitations is directly illustrated by Fig. 3 in Schreiber and Tenelsen (1993). In particular, the huge number of local minima has been shown to cause a glass like relaxation behaviour (Schreiber et al. 1996, Pérez-Garrido et al. 1999) and non-ergodic effects (Díaz-Sánchez et al. 2000a).

However, under certain conditions, also “equilibrium” observables can reflect the many-valley character of the configuration space. We focus here on the specific heat $`c`$ (heat capacity per site), in particular on its temperature dependence $`c(T)`$ for various dimensions $`d`$. It has been studied numerically by several groups (Baranovskii et al. 1982, Davies et al. 1982, Davies et al. 1984, Schreiber and Tenelsen 1993, Tenelsen and Schreiber 1994, Möbius and Pollak 1996, Möbius and Thomas 1997), mainly by means of the Metropolis procedure (Metropolis et al. 1953). However, this approach is inappropriate in the low-temperature region. One of the reasons for its failure is a direct consequence of the many-valley structure of the configuration space: the simulations usually get stuck in some local minimum. Consequently, other methods have to be used. Thus, almost complete sets of system states were obtained for the direct calculation of the low-$`T`$ behaviour of $`c(T)`$ by means of different algorithms in Schreiber and Tenelsen (1993), Tenelsen and Schreiber (1994), Möbius and Pollak (1996), Möbius and Thomas (1997), compare also Mochena and Pollak (1991) and Pérez-Garrido et al. (1997). Unfortunately, these specific heat studies have so far been restricted to the three-dimensional lattice model.

To obtain an analytical description of $`c(T)`$, Baranovskii et al. (1980) considered the excitation from the ground state by independent single-particle hops. In the case of randomly positioned sites, the adequacy of this approach was supported by numerical experiments (Baranovskii et al. 1982). However, for the three-dimensional lattice case, Möbius and Pollak (1996) showed (though only for one value of the disorder strength) that this description considerably underestimates $`c(T)`$.

The aim of our present work is to study the temperature dependence of the specific heat $`c(T)`$ in the general case, considering in particular also that $`T`$ region which is not accessible to usual Metropolis simulations. We have simulated one- to three-dimensional samples with both random and lattice structures, and, in the latter case, we have varied the strength of the disorder. To illuminate the puzzling situation described in the previous paragraph, we have sought for the conditions under which $`c(T)`$ can be understood in terms of independent single-particle hops. For that purposes, the computer code which was developed in Möbius and Thomas (1997) for the three-dimensional lattice case has been rewritten and improved to handle arbitrary structures.

## 2 Model

All situations which are to be studied can be described by the Hamiltonian

$$H=\underset{\alpha }{}\epsilon _\alpha n_\alpha +\frac{1}{2}\underset{\alpha \beta }{}\frac{(n_\alpha K)(n_\beta K)}{\left|\stackrel{}{x}_\alpha \stackrel{}{x}_\beta \right|}$$
(1)
where the first sum represents an on-site disorder contribution and the second sum describes the Coulomb interaction, see e.g. Shklovskii and Efros (1984). The sums are over the $`N`$ sites of the sample. These sites positioned at $`\stackrel{}{x}_\alpha `$ are arranged either on a lattice (lattice model), or randomly (random-position model), the more realistic case. The occupation number $`n_\alpha \{1,0\}`$ indicates whether or not site $`\alpha `$ is occupied by one of the mobile charges – note that, in a donor band, $`n_\alpha =1`$ corresponds to a missing electron. In the lattice model, the $`\epsilon _\alpha `$ describe only the on-site disorder potential; they are uniformly distributed between $`B/2`$ and $`B/2`$ where the parameter $`B`$ characterises the strength of the disorder. In this case, for electroneutrality, background charges $`K`$ are attached to each site where $`K=n_\alpha `$. However, for an n-type impurity band emulated by means of the random-position model, the $`\epsilon _\alpha `$ denote the potentials of the charged acceptors and we have $`K=0`$. Elementary charge, spatial density of the $`\stackrel{}{x}_\alpha `$, and dielectric constant are all taken to be 1.

In order to reduce finite-size effects and to exclude the influence of surfaces, we have modified the interaction term in our simulations by imposing periodic boundary conditions and using the minimum image convention (Metropolis et al. 1953, Davies et al. 1982). Thus our sample is an elementary cell in a (super)lattice, and the interaction energy of two sites is assumed to be determined by the shortest distance within this repeated representation. Moreover, we take into account particle exchange with a reservoir. That means we consider the grand canonical potential $`h=H\mu _\alpha n_\alpha `$, where the chemical potential $`\mu `$, if it is not already fixed by symmetry reasons, is determined by relaxing the closed system before starting the Metropolis part of the simulation. Considering the grand canonical ensemble rather than the canonical has several advantages. Correlation times are reduced by emulating in this manner also long-range hops. Moreover, in particular at high $`T`$, the influence of finite-size effects is diminished. However, one must remember that any $`T`$-dependence of $`\mu `$ may influence $`c(T)`$, as it does in metals. We have taken care to exclude such effects in the tests of our code.

## 3 Algorithm

In order to treat correlations completely, we have studied the specific heat of the Coulomb glass by numerical simulations. The standard procedure for this is the Metropolis algorithm (Metropolis et al. 1953) which consists of substituting time averaging for ensemble averaging. For related surveys see e.g. Binder (1984) or Binder (1997), and references therein. However, in this approach, at low $`T`$ and close to phase transitions, large difficulties arise from diverging correlation times. In the low-$`T`$ case, these difficulties are reduced to a large extent by the hybrid algorithm developed in Möbius and Thomas (1997). This bridges the gap between Metropolis simulation and analytical statistical mechanics. In the present work, we use an improved version of this procedure.

The difficulties mentioned above are tackled here in two ways following Möbius and Thomas (1997). First, in order to handle the problem of the great portion of unsuccessful attempts at low $`T`$, we do not perform random trials of state modifications at discrete time values as in the original Metropolis approach. Rather, on a continuous time scale, we focus on the state modifications which are really performed, though each such state investigation is clearly more elaborate than the trials in the standard procedure: To obtain the mean dwell time at the current state, the transition probabilities to all of its neighbours in the configuration space are calculated; the next state (in time) is picked randomly according to the ratios of the various transition probabilities, see Möbius and Thomas (1997). This approach resembles the n-fold way algorithm (Bortz et al. 1975). It becomes increasingly advantageous in comparison to the standard Metropolis procedure as $`T`$ is lowered. Moreover, problems which can arise from time discretization (Adam et al. 1999) are avoided in this way.

Second, our hybrid algorithm takes advantage of the fact that thermodynamic averages are independent of the dynamics considered provided only that detailed balance is guaranteed. Thus, in order to reduce the problems arising from the huge number of local minima, it decomposes the configuration space into low- and high-energy parts, and utilises very fast auxiliary transitions additionally introduced between the states in the low-energy part. Therefore, the sample simulated thermalises very rapidly with respect to the low-energy part of the configuration space as soon as its state enters this region: Due to the fast auxiliary transitions, the dwell time in this low-energy part as a whole is large in comparison to the dwell time in one of the states contained in it. Thus all quantities can be assumed to have the corresponding thermal averages of the low-energy subset of the configuration space while the sample stays in this configuration space region. Moreover, the escape probability from this subset does not depend on the entrance state. For details we have to refer to Möbius and Thomas (1997). For related approaches see Frantz et al. (1990) and Krauth and Pluchery (1994).

In the present work, we have modified these ideas in three ways. We have used an improved algorithm for obtaining an almost complete set of low-energy system states. In Möbius and Thomas (1997), this subset was constructed by first sampling local minima, and then stepwise improving this subset by adding all the low-energy states found during the simulation, see also Schreiber and Tenelsen (1993). Here, we have extended the subset construction as introduced in Díaz-Sánchez et al. (2000a). That means, we first sample local minima; then we improve this backbone of the low-energy subset by thermal cycling (Möbius et al. 1997) (one could also compare this optimisation procedure with compaction of a granular medium, such as ground coffee, by shaking). Finally, we search systematically the surroundings of all low-energy system states that have been found so far for further states (Möbius and Pollak 1996, Pérez-Garrido et al. 1997). Also in our new code, we add the states visited during the Metropolis simulation to the low-energy states set if their energy is lower than the highest energy value in this subset. That this happens very seldom in practice, confirms the efficiency of the three-stage low-energy subset construction procedure, described above.

As has often been stressed by Mike Pollak, see e.g. the review Pollak and Ortuño (1985), cluster excitations are expected to play an essential role in the low-$`T`$ dynamics of the Coulomb glass. In a direct way, one can include only small, compact clusters into the Metropolis dynamics. Unfortunately, due to the long-range character of the interaction, the usual cluster algorithms as proposed in Swendsen and Wang (1987) seem not to be applicable. Therefore, we use another approach here: Comparing repeatedly small numbers of states from the low-energy subset (mainly picked randomly), we search for groups of sites which change occupation only simultaneously, see Kawashima and Suzuki (1992). These clusters are tabulated and included into the simulation dynamics.

Finally, in order to prevent the low-$`T`$ improvements, described above, from slowing down too much the computer code in the high-$`T`$ case, we have introduced a $`T`$ dependent neighbourhood pruning: The number of the transitions included in the Metropolis dynamics is appropriately reduced with increasing temperature.

## 4 Numerical results

A series of tests has been performed to check and optimise the new code. The quality of the sets of low-energy states was confirmed for various degrees of disorder by series of runs in which we treated the same sample repeatedly. Moreover, we simulated lattice samples without disorder and reproduced the ordered ground states. In these tests, when we considered half-filled three-dimensional samples of $`(2m)^3`$ sites ($`m`$ integer), we observed the NaCl structure to be the ground state only for odd $`m`$. The six-fold degenerate layered structure has slightly lower energies for arbitrary even $`m`$. At present, we cannot answer the question whether this strange feature of the long-range Coulomb interaction is only an artifact of the boundary conditions, or whether it has some physical meaning. However, it is remarkable in this context that layered structures have been observed in quench simulations of one-component plasmas (Ogata and Ichimaru 1989). Additionally, we have considered lattice samples without disorder where only a fraction $`2^p`$ of the sites is occupied ($`p`$ integer), studied analytically in Hunt and Pollak (1986). For $`d=3`$, $`p=2`$, and $`m=4`$, we obtained ground states with layered structure.

We turn now to the discussion of the $`c(T)`$ data obtained in the simulations. We have compared our results on the random-position model with those obtained previously by Baranovskii et al. (1982): we observed very good agreement except for slightly larger peak values in our case. (In the present study, the smallest considered $`T`$ value is lower than in Baranovskii et al. (1982) by a factor 16.) Furthermore, also the comparison with lattice data from Möbius and Thomas (1997) showed perfect agreement within the statistical accuracy.

Finite-size effects can considerably complicate Coulomb glass simulations due to the long-range character of the interaction. Such problems are well-known from the study of the single-particle density of states where the sample size is related to a cut-off energy of the interaction spectrum and thus to a reliability limit in the energy dependence of this density of states (Baranovskii et al. 1979, Möbius et al. 1992). In consequence, the thermal smearing of the Coulomb gap is connected with a correlation length diverging as $`T0`$ (Hunt 1990). However, those correlations and corresponding processes which are essential for the formation of the Coulomb gap need not also be the most important for the specific heat. Other state modifications rather than single-particle hops can be far more important for this quantity. This assertion is suggested by the numerical data from Möbius and Pollak (1996) showing that, for $`d=3`$, $`c(T)`$ data seem to be reliable down to $`T`$ values far smaller than the cut-off energy of the interaction spectrum.

To illuminate this point, we have plotted $`c(T)`$ for various sample sizes in Fig. 1, considering the lattice model with fixed disorder strength with $`d=1`$, 2, and 3. For $`d=1`$, the curves for $`N=12`$, and 400 sites agree perfectly (within the statistical accuracy of our data) above as well as below the cut-off of the interaction energies, $`2/N`$, of the smaller sample. On the other hand, considerable finite-size effects are present for $`N=3`$ over the range $`T<1`$. Thus the specific heat is mainly determined by short-range single- or multi-particle hops here. There are large $`c(T)`$ fluctuations in the smaller samples at the lowest $`T`$ values (Tenelsen and Schreiber 1994, Möbius and Pollak 1996). These fluctuations seem to be “localised” in the big samples.

In the case $`d=2`$, significant finite-size effects are present for $`N=9`$ and $`N=25`$. However, the deviation between the data for samples of 100 and 900 sites is very small, just at the accuracy limit. For $`d=3`$, finite-size effects are clearly larger than in the previous cases. However, the results for samples of 512 and 1000 sites agree well with each other.

The convergence properties shown in Fig. 1 yield a rough estimate of the correlation radius: $`l_{\mathrm{corr}}=3\mathrm{}5`$. This length seems to increase “slowly” with decreasing $`T`$.

Fig. 1 shows that the finite-size effects, which arise from the Coulomb interaction, rise with increasing $`d`$. To illuminate the interplay of interaction and dimension directly, Fig. 2 compares $`c(T)`$ data for various $`d`$ from lattice simulations with fixed disorder strength with the interaction-free case. In the latter situation, analytical calculations yield that $`cT`$ for $`T0`$ and that $`cT^2`$ if $`T\mathrm{}`$, in agreement with the numerical data. Switching on the Coulomb interaction causes $`c`$ to decrease at low $`T`$ and to increase at high $`T`$, the more the larger $`d`$ is. In the low-$`T`$ region, $`c(T)`$ decreases faster than linear as $`T0`$. Though $`c(T)`$ cannot be approximated by a simple power law in the low-$`T`$ region studied, our data raises the question of whether $`c(T)`$ might tend to a linear dependence again if $`T`$ is sufficiently lowered.

The influence of disorder strength is studied in Figs. 3 and 4, showing $`\mathrm{log}_{10}c`$ and $`c`$, respectively, in dependence on $`\mathrm{log}_{10}T`$. The observation that, for high $`T`$, increasing $`B`$ reflects as some shift of $`c(T)`$ towards higher $`T`$, which was reported for $`d=3`$ in Möbius and Thomas (1997), is made here also for $`d=2`$ and $`d=1`$. However, at low $`T`$, there is a weaker, but more complicated dependence on the disorder strength. In all cases, the random-position model yields broader $`c(T)`$ curves than the lattice model. At low $`T`$, these data are close to the results obtained for low disorder strength in the lattice model.

A particular feature of $`c(T)`$ is only seen clearly when $`c`$ is plotted on a linear scale as in Fig. 4: In the lattice model, for $`B=0.5`$ and $`d=2`$ as well as $`d=3`$, the $`c(T)`$ curve is narrower and considerably higher than for the other $`B`$ values. The question, whether or not this behaviour might be a precursor of a phase transition for still lower $`B`$ (Grannan and Yu 1993, Vojta and Schreiber 1994, Díaz-Sánchez et al. 2000b) cannot be answered yet. Corresponding simulations are extremely CPU-time consuming.

Finally, we address the question of which processes determine $`c(T)`$ as $`T0`$. To clarify which contributions to the specific heat yield single-particle hops (SPH) and multi-particle hops, respectively, we have calculated the value $`c_{\mathrm{SPH}}`$, which one obtains when considering SPHs to be the only relevant, independent excitations of the ground state. The $`T`$-dependence of the quotient of $`c_{\mathrm{SPH}}`$ and of the “exact” $`c(T)`$ data is represented in Fig. 5 for various $`d`$ and $`B`$, for the lattice as well as for the random-position model. The SPH approximation obviously begins to fail when $`T`$ is sufficiently increased for a significant number of sites to belong simultaneously to more than one of the excited SPHs. Therefore, Fig. 5 focuses on that region, where the corresponding overlap correction does not exceed $`0.1c`$; quotient values for higher $`T`$, included in certain cases only, show that the breakdown of the above condition is correlated with a qualitative change of the quotient’s $`T`$ dependence. The maximum applicability $`T`$ of the SPH approximation decreases with increasing sample size. Therefore, for $`d=1`$ and 2, medium sample sizes are considered here. Nevertheless, quotient values could only be obtained for a $`T`$ region which cannot be studied by standard Metropolis simulations, and, for the random-position model in the case $`d=1`$, the applicability region of the SPH approximation includes only one data point.

Fig. 5 shows that the SPH approximation, when applicable, clearly underestimates $`c`$ for all but one case; the random-position model for $`d=1`$ is the only exception. In more detail, in the two- and three-dimensional lattice models, multi-particle contributions are particularly important for weak disorder. In the random-position model, the relevance of multi-particle excitations increases with the dimension.

Thus, Fig. 5 implies the conclusion that, in agreement with Mike Pollak’s view on the Coulomb glass model (Pollak and Ortuño 1985), multi-particle processes yield an essential contribution to the specific heat in the low-temperature limit.

Acknowledgments

Part of this study was performed during a stay of AM at the California State University, Bakersfield. He thanks the US National Science Foundation for financing this visit under grant DMR-9803686, and the colleagues in Bakersfield for their hospitality. Moreover, we are much obliged to A. Hunt for critical discussions, and to U. Kleinekathöfer, S. Köstlmeier, and M. Schreiber for the possibility to perform part of the simulations at workstations of the Technische Universität Chemnitz, and for various related help.

References

Adam, E., Billard, L., and Lançon, F., 1999, Phys. Rev. E, 59, 1 212.

Baranovskii, S.D., Efros, A.L., Gelmont, B.L., and Shklovskii, B.I., 1976, J. Phys. C: Solid State Phys., 9, 2021.

Baranovskii, S.D., Efros, A.L., Gelmont, B.L., and Shklovskii, B.I., 1979, J. Phys. C: Solid State Phys., 12, 1023.

Baranovskii, S.D., Shklovskii, B.I., and Efros, A.L., 1980, Zh. Eksp. Teor. Fiz., 78, 395 \[Transl.: 1980, Sov. Phys. JETP, 51, 199\].

Baranovskii, S.D., Uzakov, A.A., and Efros, A.L., 1982, Zh. Eksp. Teor. Fiz., 83, 756 \[Transl.: 1982, Sov. Phys. JETP, 56, 422\].

Binder, K., editor, 1984, Applications of the Monte Carlo Method in Statistical Physics, in Topics in Current Physics, Vol. 36, (Berlin: Springer).

Binder, K., 1997, Rep. Prog. Phys., 60, 487.

Bortz, A.B., Kalos, M.H., and Lebowitz, J.L., 1975, J. Comp. Phys., 17, 10.

Davies, J.H., Lee, P.A., and Rice, T.M., 1982, Phys. Rev. Lett., 49, 758.

Davies, J.H., Lee, P.A., and Rice, T.M., 1984, Phys. Rev. B, 29, 4260.

Díaz-Sánchez, A., Möbius, A., Ortuño, M., Neklioudov, A., and Schreiber, M., 2000a, Phys. Rev. B, 62, 8 030.

Díaz-Sánchez, A., Ortuño, M., Pérez-Garrido, A., and Cuevas, E., 2000b, phys. stat. sol., (b), 218, 11.

Efros, A.L., and Shklovskii, B.I., 1975, J. Phys. C, 8, L49.

Efros, A.L., and Shklovskii, B.I., 1985, Electron-electron interactions in disordered systems, edited by Efros, A.L., and Pollak, M., (Amsterdam - Oxford - New York - Tokyo: North Holland), p. 409.

Frantz, D.D., Freeman, D.L., and Doll, J.D., 1990, J. Chem. Phys., 93, 2769.

Grannan, E.R., and Yu, C.C., 1993, Phys. Rev. Lett., 71, 3 335.

Hunt, A., and Pollak, M., 1986, Phil. Mag. B, 53, 353.

Hunt, A., 1990, Phil. Mag. Lett., 62, 371.

Kawashima, N., and Suzuki, M., 1992, J. Phys. A: Math. Gen., 25, 1055.

Kogan, S., 1998, Phys. Rev. B, 57, 9 736.

Krauth, W., and Pluchery, O., 1994, J. Phys. A: Math. Gen., 27, L715.

Kurosowa, T., and Sugimoto, H., 1975, Progr. Theor. Phys. Suppl., 57, 217.

Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H., and Teller, E., 1953, J. Chem. Phys., 21, 1087.

Mézard, M., Parisi, G., and Virasoro, M.A., editors, 1987, Spin glass theory and beyond, (Singapore: World Scientific).

Möbius, A., Richter, M., and Drittler, B., 1992, Phys. Rev. B, 45, 11 568.

Möbius, A., and Pollak, M., 1996, Phys. Rev. B, 53, 16 197.

Möbius, A., and Thomas, P., 1997, Phys. Rev. B, 55, 7 460.

Möbius, A., Neklioudov, A., Díaz-Sánchez, A., Hoffmann, K.H., Fachat, A., and Schreiber, M., 1997, Phys. Rev. Lett., 79, 4 297.

Mochena, M., and Pollak, M., 1991, Phys. Rev. Lett., 67, 109.

Ogata, S., and Ichimaru, S., 1989, Phys. Rev. Lett., 62, 2 293.

Ortuño, M., Talamantes, J., Cuevas, E., and Díaz-Sánchez, A., 2001, Phil. Mag. B, 82, in press.

Pérez-Garrido, A., Ortuño, M., Cuevas, E., Ruiz, J., and Pollak, M., 1997, Phys. Rev. B, 55, R8630.

Pérez-Garrido, A., Ortuño, M., Díaz-Sánchez, A., and Cuevas, E., 1999, Phys. Rev. B, 59, 5 328.

Pollak, M., 1970, Discuss. Faraday Soc., 50, 13.

Pollak, M., and Ortuño, M., 1985, Electron-electron interactions in disordered systems, edited by Efros, A.L., and Pollak, M., (Amsterdam - Oxford - New York - Tokyo: North Holland), p. 287.

Rujan, P., 1988, Zeitschr. f. Physik B (Condensed Matter), 73, 391.

Sandow, B., Gloos, K., Rentzsch, R., Ionov, A.N., and Schirmacher, W., 2001, Phys. Rev. Lett., 86, 1 845.

Schreiber, M., and Tenelsen, K., 1993, Europhysics Lett., 21, 697.

Schreiber, M., Tenelsen, K., and Vojta, T., 1996, J. Lumin., 66, 521.

Shklovskii, B.I., and Efros, A.L., 1984, Electronic properties of doped semiconductors, (Berlin - Heidelberg \- New York - Tokyo: Springer).

Srinivasan, G., 1970, Phys. Rev. B, 4, 2581.

Swendsen, R.H., and Wang, J.S., 1987, Phys. Rev. Lett., 58, 86.

Tenelsen, K., and Schreiber, M., 1994, Phys. Rev. B, 49, 12 662.

Vaknin, A., Ovadyahu, Z., and Pollak, M., 2000, Phys. Rev. Lett., 84, 3 402.

Vojta, T., and Schreiber, M., 1994, Phys. Rev. Lett., 73, 2 933.