Metastable Dark Current in BRITE Nano-Satellite Image Sensors

Dark current in charge-coupled devices (CCDs) is one of the most important sources of impulsive noise present in scientific images. While the dark current originating in the fabrication defects (mainly impurities) is stable and dependent only on temperature, the one present in the proton-irradiated sensors shows a range of metastable states which makes calibration of images almost impossible. In this paper, we show an extended analysis of such metastabilities present in Kodak KAI 11000M CCD sensors employed in the BRITE (BRIghtest Target Explorer) astrophysical mission over 7 years of in-orbit work. Our collection of dark current characteristics has an unprecedented time span, large temperature range and high number of investigated pixels. A special methodology based on the Gaussian mixture model was proposed for identification and characterization of the metastable states in the dark current. We identified several interesting properties of the metastability and found an experimental rule for the dark current in tristable defects. The results shed a new light on the dark current problems, its modeling and the mitigation in an image sensor working in space.


Introduction
Charge-coupled devices (CCDs) are the most widespread image sensors used in scientific applications where precise measurements of light is critical. These sensors have been employed for years since its invention by Smith and Boyle in 1969. The light-gathering technology was further adopted in CMOS (complementary metal-oxide-semiconductor) image sensors which are currently being intensively developed reaching the capabilities of CCDs. These two types of detectors are currently utilized in all imaging applications, both consumer and scientific ones. The dark current (DC) is one of the parameters of both CMOS and CCD devices [1]. This is the bias charge which accumulates in each pixel even without the photons falling on the sensor. There are two main sources of DC: diffusion of charge from the substrate and point-like defects located in the depletion region including impurities or atomic vacancies. While the first one shows a similar behavior across whole sensor, the second exposes strong nonuniformities between pixels [2]. A quantization of DC histogram could be observed in first CCD devices, as the number of point defects followed Poisson distribution [3,4]. In most of the currently fabricated image sensors, including CCDs and CMOS sensors, the diffusion DC is dominant as the impurities are successfully eliminated during the fabrication process.
Dark-frame subtraction is a typical procedure utilized in the compensation of DC. One can understand "dark-frame" to mean that the image is taken with the shutter closed such that no light falls on the matrix. The dark frame should be taken with the same exposure time and the same gain of the output amplifier. This method is almost always utilized in astronomy and in all applications

BRITE CCD Sensors
In this paper, we investigated the behavior of the DC in five proton-irradiated CCDs of the BRITE (BRIghtest Target Explorer) constellation of satellites operated in space environment for almost 7 years. The BRITE mission consists of the five nanosatellites, 20 × 20 × 20 cm each, developed by Canada, Austria and Poland and launched in late 2013 [15]. Their names and ownership are: Austria-BRITE Austria (BAb), UniBRITE (UBr), Canada-BRITE Toronto (BTr), Poland-BRITE Lem (BLb) and BRITE Heweliusz (BHr). Small "r" or "b" in the abbreviation indicates the type of utilized optical filter: red or blue. The aim of the mission is the observation of the brightest stars in the sky for long (>3 months per star) uninterrupted periods [16]. These long sets of data result in the unprecedented photometric accuracy which allows an international team of astronomers to investigate the tiny fluctuations (inaccessible for the ground-based telescopes) of stellar intensities. For more information on the BRITE mission and techniques used for data processing, interested readers can refer to the following works [15,17,18] as well as to the large collection of astronomical papers based on the data gathered by the BRITE constellation.
The image sensors installed in BRITE satellites are full-frame, interline, monochrome KODAK KAI11000M sensors. The full information about their parameters are given in Table 1. These particular sensors were used due to their size, since larger observed portion of sky allows to measure simultaneously more bright stars increasing the mission's yield. With the employed lens system, the total field of view equals 30 • × 20 • . The sensor is not thermally stabilized, which is unfortunate due to the inability to perform classic dark frames subtraction. On the other hand, this is fortunate because the observed strong swings of operating temperature allow for performing all thermal analysis as presented in our study. The image sensors are not fully read out. Only small (e.g., 30 × 30 pixels) windows are extracted around the target stars to reduce the amount of data to be transferred to the ground stations. We show several examples of such windows in Figure 1. Note the presence of a central stellar profile and noise cause by the DC background. Two types of offsets are visible: single bright pixels which are the result of the thermal charge accumulated during an exposure and the brighter columns related to the generation of the DC during vertical transfer of charge in readout cycle. While the former shows the instabilities investigated in this work, the latter does not depend on the exposure time and can be easily extracted in the photometric pipeline by the iterative calculations of the median intensities in window's columns (see Section 3.2 in [17]). All of the BRITE data collected by the satellites are organized in so-called setups, which are related to the configuration of downloaded windows and observed objects. Some of the setups are shorter, since the configuration had to be modified. Another is as long as 3 months. Currently (2020,) we have obtained over 273 setups over 7 years of mission, which gives us an enormous amount of data for performing temporal analysis of the DC.
Let us recall the most important conclusions drawn from our previous research [14]. We observed the linear increase of the mean DC with time over the whole CCD matrix. The growth was significantly lower for the satellites which utilized either 10 mm Borotron or 2 mm tungsten shields (BHr or BTr, respectively). There was some initial high jump of the DC just after the launch, which means that a large dose of radiation was absorbed by matrices while elevating the telescopes to the orbit. Finally, we identified two thermal coefficients of the DC in sensors. First, the found activation energy, ∆E, was 0.68 eV, which is a value frequently found in image sensors irradiated in CCDs. Second, the Meyer-Neldel rule [19], which is still not fully understood, was confirmed in the sensors, and its constant, E MN , was estimated at 24.8 meV.

Thermal Dependencies of Dark Current in CCDs
The DC in CCD pixels follows the well-known thermal formula [20]: where I d is the amount of dark charge per pixel per time interval, G is called a prefactor, ∆E is the activation energy of thermal generation, k is the Boltzmann constant and T is the temperature in Kelvins. To allow for finding activation energy, the relation is usually presented as the logarithm of dark charge against 1/kT, which is called the Arrhenius plot. An example of such plot is given in Figure 2a. This corresponds to the following transformation of (1): By log(), we understand the logarithm is base e; this convention is kept throughout the whole paper. By fitting a linear regression to such a dependency, one can obtain an estimation of the activation energy ∆E which can be identified with various types of defects. Our previous work [14] shows that the activation energy was found at 0.68 eV, which was close to the known energy of phosphorus-vacancy pair in the silicon, 0.7 eV. Importantly, the observed activation energies in BRITE CCDs ranged from 0.4 to 0.8 eV which may result from modulation of DC production by the electric field (i.e., defects may be present in different locations/depths within pixel structure). It was observed that the prefactor G obeys the Meyer-Neldel rule (MNR) [19]. This is an experimental dependency that is still not fully understood but is already found in CCDs [21]. The G factor follows the equation: where G and E MN are positive MNR constants (prefactor and Meyer-Neldel energy, respectively). The first parameter may be significantly different from pixel to pixel. The second should be rather constant; this was indeed observed in our matrices [14]. The obtained value of E MN equaled 24.8 meV.
Eventually, the DC equation can be expressed as: In Figure 2, we selected three plots from our previous study [14], which should be presented here for better readability of the current paper. The following plots show: (a) an exemplary logarithmic dependency of DC in a random pixel showing that it follows the equation Equation (2); (b) the histogram of calculated activation energies with a peak around 0.68 eV; (c) the plot confirming NMR rule present in our sensors with the linear fit assuming E MN = 24.8 ± 0.1 meV.

Image Processing
To obtain the DC of a pixel, several steps of image processing have to be performed. They are necessary to compensate for the bias introduced by the DC generated during the charge transfer (column offsets visible also in Figure 1). Moreover, the pixels included in the stellar profiles have to be masked, since their intensities are obviously elevated above the DC level. These steps have been already explained in details in [14] in Section 3.1; thus, the interested reader is kindly referred to the original work. In brief, the undertaken steps include (1) iterative calculation of median intensities in columns, (2) image denoising using the 3 × 3 median filter, and (3) thresholding and dilatation of the detected stellar profile. The example outcomes of the consecutive image processing steps are given in Figure 3. The following images show: original window, estimated column offsets, window after the compensation, median-filtered image and, finally, the pixels detected as covered by the stellar profile which have to be excluded from further analysis. The result after performing image processing is a set of measurements of pixels DC including: (1) generated charge, (2) temperature of a CCD, (3) exposure time, (4) position of a pixel in matrix and (5) name of a satellite. Due to the large amount of gathered image data, in each window, these information are stored only for the evident defective pixels whose intensities are larger than 9σ RON , where σ RON stands for the readout noise (the value of 9 was used to be sure that the subsequent analysis includes only pixels significantly damaged by protons). Additionally, the average DC is calculated using all pixels in a window without signals from the star. This quantity was investigated in our previous study [14], showing a linear increase for all five satellites over the first 4 years of the mission.
We found that that most of the identified defective pixels expose instabilities in Arrhenius plots. Instead of an expected single line along which the points should group, we observed a range of lines as presented in examples in Figure 4.
One can see that the plots are nearly parallel, which means that the activation energy remains constant, while the G prefactor changes between some temporarily stable states. Unfortunately, the estimation of activation energy by linear regression (even a robust version) is unreliable in such pixels; therefore, we had to develop new routines for fitting ∆E. We used the equation Equation (4) as the starting point and the already estimated MNR constant E MN = 23.8 ± 0.1 meV. We found that the DC usually remains stable for relatively long times, changing its state (i.e., its G parameter) only from time to time. The temporal stability of DC allows to identify the activation energy by finding the value of ∆E which minimizes the following absolute median derivative of logG as follows: where: i is a sample number, Ω is the range of samples and ∆E is the estimated value of the activation energy. The process of finding the best ∆E is presented in Figure 5 wherein the plot on the right side show the dependency of optimized factor O in function of assumed activation energy with a clear minimum around 0.74 ± 0.02 eV. The DC measurements are presented in the upper left plot while the same plot after applying normalization using Equation (5) is given below. It should be noted that the initial measurements are done in varying temperatures, while after normalization (Equation (5)), the temperature-related effects are eliminated and the metastable states of log G with the transitions appear evident.  We provide some examples of mestability with more than two states in Figure 6. Some of the detected faulty pixels in data obtained in the recent years of the mission show the metastability with a larger number of states which may be linked with the presence of several 2-or 3-stable defects in a pixel.

Identification of Metastable States
The goal of this study is to find and investigate the locally stable values of log G. We use the histograms of log G values calculated after normalization with Equation (5); then, we apply the Gaussian mixture model to fit a collection of Gaussian curves to the data of each setup. The fitting is performed 10 times for each defective pixel with a limit of 1000 iterations, random starting points and a rising number of Gaussians components from 2 to 12. The model quality is evaluated with the Akaike information criterion (AIC) [22] for each solution. The modeling with unknown number of distributions in a mixture is usually problematic because the models with higher numbers of Gaussians usually fit better to the data. Fortunately, the AIC coefficient allows for finding a trade-off between the information loss and the simplicity of a model. After visual investigations, we admit that this kind of criterion works well for our data. In Figure 7, we show several outcomes of the modeling of log G distributions. The calculation of models for the whole data set required almost three weeks of uninterrupted work of an 8-core I7 CPU.
After the GMM (Gaussian mixture model) modeling, we obtained a set of information about each pixel, from which the most important are: the number of metastable states, the mean values of log G for each state and the relative proportion of components which represents the division of time spent in the individual states. These quantities were evaluated and their results are presented in the next section.

Average Dark Current
First, we would like to update on the current status of the progress of DC in BRITE CCDs relative to the results presented 2 years ago. The average DC was calculated per pixel per second at room temperature using the approach presented in [14], i.e., the Arrhenius plot for the average DC in a given setup was fitted with a robust linear fit to eventually find its value at 25 • C. The results are given in Figure 8. The time scale was adjusted so that zero appears at the launch of each satellite. The Arrhenius plots for the average DC do not show any RTS-like grouping as in Figure 4 since the pixels are analyzed together and the effects cancel out.
We can observe that the linear increase of DC generation rate is still observed as it was in [14]. We found that the values of DC growth were calculated the wrong way in the previous study-they were not normalized by the total number of pixels but by the number of defective pixels. This leads to significantly lower values presented in the current work. The observed linear growth implies that there is no visible mechanism of defects annealing in the operating room temperatures for this type of defects. This fact shows that our defects are different from those of the Hubble Space Telescope [23], where most hot pixels are annealed at 30 • C. They are probably of the type presented in [5] where the authors report the annealing process at 100 • C. Unfortunately, we are currently unable to rise the temperature to this levels due to the safety reasons. The annealing may be considered in the future, if the number of defects disabled valuable astronomical observations. The linear growth of the DC in CCDs working in space is also reported by authors in [24]. The two shielded satellites (BHr and BTr), as well as the one unshielded (UBr), show a similar increase of DC with time. However, the three corresponding plots differ greatly. The BHr shielded with 10 mm Borotron behaves stably over time and the plot shows low scatter. This means that the DC shows similar statistics over the whole matrix, so the positions of rasters do not have a great impact on the calculated mean DC. This is much different in the BTr satellite, where strong scatter indicates that, although the mean DC shows similar growth, the calculations depend heavily on the selected set of rasters. Interestingly, we found that the observed accuracy of the astronomical photometry in BTr is closer to that obtained in the unshielded satellites. This is because the unshielded satellites and BTr have a similar number of defective pixels, while in BTr, only the DC generation is lower. The last of the three satellites showing slower growth is UBr-the most mysterious one. For the first 4 years of mission, the points were arranged along the line showing a growth similar to the one seen in the remaining unshielded satellites. However, the points from the 3 most recent years show the average DC much smaller than predicted; they are even at the same level. This is the only one satellite which shows this kind of saturation. There is no significant difference in temperature of operation among satellites which would result in a bigger chance for annealing, so the source of this phenomenon remains inexplicable. Interestingly, the three satellites exposing slower generation growth include the red filter. They are installed in the closest proximity of CCD matrix, so if the red filters and the blue ones differ in the absorption properties, they may have some additional impact on the observed degradation process.

Statistics of Metastable Pixels
We evaluated the number of blinking pixels within a set of defective pixels, i.e., the ones whose DC exceeds 9× standard deviation of readout noise. From this set, we extracted pixels for which the GMM returned at least two metastable states and analyzed their relative number for each satellite. The results of this analysis are presented in Figure 9. The plot on the left side shows the distribution of blinking pixels with a given minimal number of states. We can see that most of identified pixels, around 80%, are indeed blinking with at least two metastable states. This is true for both shielded and unshilded satellites. However, the BHr shielded with 10 mm Borotron is significantly different from the rest. Although the percentage of all unstable pixels (i.e., with at least two states found) is similar, the contribution of pixels with a higher number of states is visibly lower (e.g., the pixels with more than five states are 20% of the total, while in the other satellites, they constitute more than 30%).
The plots in the middle and on the right side of Figure 9 show the evolution of the total number of defective and blinking pixels, respectively. BTr follows the growth of the unshielded satellites, although its total DC is significantly lower as presented in Figure 8 and in Table 2. As expected, the percentage of damaged pixels increases with time-slower for BHr and faster for the other satellites. Currently, we reach a total of 10-15% of pixels affected by higher DC in BAb, UBr, BLb and BTr, while in BHr, this value is closer to 5-10%. Importantly, the total percentage of blinking pixels (at least two metastable states) in the defective set remains constant for all years of mission; this can be seen in the right plot of Figure 9. This value oscillates around 80%. The spread of points in the two plots is natural and appears due to the limited number (15)(16)(17)(18)(19)(20) of rasters (usually 28 × 48 pixels around the brightest stars) downloaded from each image, which, in combination with a low percentage of detected pixels, introduces the observed level of uncertainty.

States Analysis
In this section, we would like to concentrate on the statistics of states over detected metastable pixels. The results of time analysis of these pixels presented in the previous section indicates that the population of such pixels keeps the same internal proportions over the lifetime of the mission (only the total number of such pixels increases). This allows us to analyze the results of GMM together, regardless of the exact time a pixel was measured. In our set, each single pixel has K identified states with the estimated respective prefactors G i (i = 1, ..., K) and activation energy ∆E. Since the values of G show a range of two magnitudes, we decided to operate in logarithmic scale, i.e., we use log G instead. Summarizing, we collected a great number of measurements and unprecedented set of various observations of the metastablities in the DC. The set contains over 200,000 pixels with identified RTS behavior in function of both time and temperature.
First, we start with the statistics of log G against fitted activation energy. The 3D histograms of dependency between the identified activation energy and the values of the prefactors are presented in Figure 10. In this analysis, we include the pixels with at least two metastable states which constitute the majority (∼80%) of detected faulty pixels. The cumulative histograms of the data projected on both the x and y axes are plotted at the respective sides of each 3D histogram. This technique of data presentation was adopted in all consecutive plots in this part of the paper. We can see that most of the data points are concentrated around the activation energy of 0.7-0.75 eV with a peak at approximately 0.72 eV and log G in between 6-7. The shielded satellites have smaller log G closer to 6. In these two satellites, we can see that the smaller values of prefactor are present, even as small as log G = 3. In contrast, in BAb and UBr, the values lower than 4 are not present, while in BLb, some appear only occasionally. This observation indicates that utilization of shielding reduces the intensity of DC due to the lower values of prefactor log G. For bistable defects, we also analyzed the dependency between the prefactors of higher and lower states. The resulting histograms are presented in Figure 11. It can be seen that the differences log G 1 − log G 2 show a Poisson-like distribution with a long tail. The most frequent values of this difference equal approximately 0.1; however, values two or three times larger are also present. Interestingly, the values are virtually never higher than 0.5. The shielded satellites BTr and BHr have the average value of the difference lower than 0.1, while the others have this value slightly elevated. Still, the number of defects in shielded satellites is still not big enough to draw any evident conclusion (the histograms are rather noisy). In BAb and UBr, we can see that the difference increases as the prefactor log G rises. This is, however, not observed in the remaining satellites. Figure 11. Dependencies between the prefactors in bistable defects.
Next, we checked the proportion of states in which bistable defects persist; the corresponding histograms are presented in Figure 12. We can see that there is no dependency between the value of higher state (log G 1 ) and the percentage of time a pixel spends in it. However, we observed a small but evident predominance of high state, indicating an average proportion of 60%/40% between the two states. We could not see any differences between the satellites, except for maybe a slightly flatter distribution of calculated proportions in BTr.
After checking bistable defects we concentrated on the pixels with identified 3 metastable states. These two groups-bistable and tristable-are our groups of interest since they have to be related to the presence of a single defect in a pixel. Their number is also relatively big in our data to allow for drawing reliable conclusions. The pixels having a larger number of states are likely to be the effect of multiple 2-or 3-stable defects, and their number is smaller, thus they were not investigated further. For tristable defects, we checked how the prefactors for the three states behave. After the GMM modeling, we sorted the prefactors in the descending order: log G 1 , log G 2 and log G 3 and checked the dependency between the differences log G 1 − log G 2 and log G 2 − log G 3 . Although the corresponding histograms seen in Figure 13 are more noisy due to the number of detected tristable pixels, we found that the points are mostly arranged along the line x = y, which indicates an interesting experimental rule given by the following formula: This rule was not observed so far by any researcher; thus, it needs more investigation and physical interpretation in the future. We can only speculate that the mechanism that governs the transitions from state A to B is the same as for the transition from B to C. Our rule indicates that the relative decrease of the DC is the same for both transitions. Unfortunately, the metastability itself is still poorly understood, and several competitive explanations exist with no clear winner. The Meyer-Neldel rule is also an experimental rule which has no evident physical interpretation.

Dark Current Modeling
Using the relations observed in our data, we can derive the way of modeling DC in CCD sensors of the BRITE nanosatellites. This will allow us to better imitate the impulsive noise in our CCDs to evaluate future processing pipelines. This model can be also utilized by the authors of existing and new denoising techniques dedicated for removal of impulsive noise. With the presented steps, one can simulate any scenerio of DC present in an orbital CCD in a given temperature and at a defined time spent in the orbit. Due to the high complexity of our data, we made some simplifications, like fitting with 2D normal distribution, for better accessibility of the model for interested researchers.
First, we start with the definition of the percentage of defective pixels, N d . This value can be extrapolated using the initial growth of such pixels presented in the middle plot of Figure 9. Currently, the highest N d are at 15%. The positions of pixels are uniformly random across the whole CCD matrix. Next, we have to distribute the number of states between the defective pixels. We can use a generation of random variables utilizing the inverse transform method. The cumulative distribution function (CDF) has been already presented in the left plot of Figure 9. We found that the CDF for all satellites can be successfully fitted with the conditional exponential relation given by the following formula: The parameter σ for BLb, UBr, BTr and BAb is 0.2, while for the BHr, it is 0.15. Finally, we have to decide how to distribute the prefactors log G if a pixel has a given number of metastable states:

1.
For pixels having a single state, we use the distribution of both activation energies and prefactors as presented in Figure 10. It can be fitted with 2D normal distribution centered at log G = 5.59 (σ log G = 1.43) and with ∆E = 0.72 (σ ∆E = 0.15).

2.
For pixels having two states, the highest one should be generated using the distribution given above, while the second, lower state should be calculated using the relation obtained experimentally from our data (based on the data in the collective histogram of (log G 1 − log G 2 ) presented in Figure 13): where P (λ) denotes a random number generated with the Poisson distribution and with the expected value of λ.

3.
For pixels with 3 metastable states, the first two states should be generated using the method given for the bistable pixels, while the third state should be calculated using the already discussed rule Equation (7). 4.
The pixels with a higher number of states are relatively rare, and we propose to simplify the procedure by taking random states from the distribution of single-state pixels. In reality, these pixels may appear due to the combination of 2-and 3-stable defects with some probable overlapping of states. However, this is still unclear, and we did not find direct evidences for this statement in our data. A deeper research on such pixels is planned in the nearest future.

5.
Since there is no evidence of higher probability of observation of a given state, its selection can be made with a uniform probability.

Conclusions
In the paper, we analyzed the behavior of pixels affected by the DC in BRITE nanosatellite charge-coupled devices. The five CCD sensors have been working in orbit for 7 years developing 5-15% of defective pixels. The most of these pixels expose multistability, which means that the thermal generation of charge shows sudden transitions between 2, 3 or more metastable states, depending on a pixel. The time span of the mission, together with the number of 200,000 analyzed pixels, allowed us to draw several interesting conclusions about the mutlistability. The most important conclusions from our study are listed below: 1.
The average DC grew linearly for all 7 years of the mission, reaching currently 5-15% of matrices affected by DC, depending on the shielding strategy.

2.
Most of the pixels showing elevated DC also exposed the multistability. Around 80% of all defective pixels are unstable, and this percentage has been observed in all satellites for the whole mission's lifetime.

3.
The shielding with 10 mm Borotron reduced the number of defects leading to both, a lower number of faulty pixels and a lower number of pixels exposing complex multistability (>2 metastable states).

4.
The shielding with 2 mm Tungsten reduced the average DC within the matrix on BTr satellite. However, the number of affected pixels was similar to the one observed in unshielded sensors. We suspect that this shield reduced the energy of protons but did not absorb them totally.

5.
From the two points given above, a clear conclusion for the future missions of small satellites can be drawn: keep place for the shielding with light materials, as they seem to be the best solution for avoiding problems with both blinking pixels and high levels of dark current. 6.
In bistable defects, both states are observed nearly with the same probability, with only a slight dominance of higher state. 7.
The difference between the prefactors log G 1 and log G 2 in bistable defects is usually between 0.1-0.2. 8.
The defects showing three metastable states can be characterized by a rule for the prefactors: G 1 /G 2 = G 2 /G 3 . This relation should be explored in the future research to find the physical explanation.
We also derived a simple model for generating impulsive noise in images affected by the DC, imitating the scenario of BRITE CCDs. The model allows for simulating the behavior of pixels in a defined temperature and as a function of degradation degree of a matrix. Since it follows the distribution of noise in real-world sensors, it can be utilized for evaluation of techniques dedicated for image filtration in space applications.