1.

Science Motivation

The Earth’s atmosphere is essentially opaque to ultraviolet (UV) radiation at wavelengths shorter than about 300 nm. All sorts of interesting astrophysical phenomena exist to be imaged in the UV, such as emission lines from interstellar gas, redshifted Ly-$\alpha$ from galaxies, and zodiacal light caused by Rayleigh scatter of sunlight. Most of the UV sources mentioned above have very low-UV surface brightness, so observation requires an imaging telescope above the Earth’s atmosphere that integrates for long periods, coupled with an imaging detector with very high efficiency. Typical photon rates are on the order of a photon per $1000\text{pixels}$ per second. As space-based missions are limited in the returned data, compression of data consisting of nondetections is relevant.

Electron-multiplying charge-coupled devices (EMCCDs) are the enabling technology for space-based UV imaging, especially with enhanced UV quantum efficiency. The Teledyne e2v CCD201-20 EMCCD has flight heritage from the Faint Intergalactic Redshifted Emission Balloon.^1–3 Roman coronagraphic instrument has extensively tested and qualified the CCD201-20.⁴ Other space applications using the CCD201-20 are SPARCS,^5–7 currently in production, and the Polarized Zodiacal Light Experiment concept.⁸

2.

What is an EMCCD?

The EMCCD is a CCD modified to achieve high signal-to-noise ratio (SNR) by rendering the read noise effectively zero. Compared to conventional CCDs, EMCCDs have an additional serial register (604 extra charge-coupled “pixels” [see Fig. 1]), where one of the register clocks is replaced by a high-voltage clock (25 to 50 V). The higher voltage causes a multiplication process that stochastically turns one electron into many, resulting in thousands of electrons at the output amplifier. This allows detection of single-photon events by thresholding above amplifier read noise.¹⁰

Fig. 1

Device architecture diagram for the e2v CCD201-20 EMCCD.⁹ The device has two separate output amplifiers, one of which is amplified with the 604 multiplication elements.

In the amplification process, if the charge transfer from one register to the next is done with enough input energy, a signal electron will knock lose another electron. This is the avalanche photodiode (APD) effect. Repeating this process makes an “APD staircase,”¹¹ magnifying the signal by having some extra CCD-like charge transfers explicitly run using overdrive. Since the digitization takes place after magnification, the read noise is proportionally smaller when compared to the original, unamplified signal.

EMCCDs are generally run in “photon counting” (PC) mode, meaning that the detector is read out often enough that the expected signal counts in each pixel are $\ll 1$. After magnification, the final counts are generally either very large or very small, signifying detecting or not detecting signal photons. This is similar to a photomultiplier tube (PMT). The downside, as compared to a PMT or an intensified CCD, is that the dark current and clock-induced charge (CIC) may also add significant noise. The EMCCD data output is then a combination of science signal, dark current, CIC, and readout noise. Dark current decreases with temperature, but sufficient cooling is not always available to make the dark current negligible. CIC is fixed per exposure and is discussed in Sec. 6. Without PC, sufficient EM amplification will make the read noise comparatively small, but the discussion in Sec. 11 shows that in PC mode, the read noise is eliminated with the proper choice of parameters.

3.

Probability Notation

The rest of this paper mainly consists of probability theory related to processing and understanding the signal from these detectors. As the intended target of this paper is engineers and scientists rather than mathematicians, neither mathematical rigor nor mathematically rigorous notation are prioritized. The notation used is as follows:

The various symbols used are, in order of definition:

4.

Poisson Distributions

For independent events arriving at a detector, the Poisson distribution is appropriate to use. For a Poisson-distributed random variable $X$, $E(X)=V(X)=\lambda$, and $X$ is distributed as

Similar to the Gaussian assumption, the SNR goes as the square root of the counts as a result of $E(X)=V(X)$. For no events,

Useful is the cumulative distribution function, which is

Here $Q(n,\lambda )$ is the regularized gamma function; e.g., evaluated in Excel as

5.

Photon Counting

Consider an ideal photon-counting detector with infinite gain. If any photons at all are detected, then the detector returns 1, otherwise 0. Let the random variable $X_W$ be the number of electrons in the well of a pixel. Using Eq. (2), the expected value (and second and all higher moments) of this process is then

In practice, a single exposure is not useful. Rather, $\lambda$ is estimated by taking a sequence of $N$ exposures, each in PC mode. Define the $\epsilon$ as the probability of detection

The number of exposures with a signal detection are then given by the binomial distribution

Here $c$ is the total number of counts for the $N$ exposures. If $\epsilon$ is ½, then this is the distribution of the number of heads after $N$ flips of a fair coin. What we want to know is $\epsilon$, which is on the wrong side of the equation. Using Bayes’ theorem,

We need a prior for $P(\epsilon )$ and select the most commonly used prior for the binomial distribution, the beta distribution, given by

The denominator is the beta function $B(\alpha ,\beta )=\mathrm{\Gamma }(\alpha )\mathrm{\Gamma }(\beta )/\mathrm{\Gamma }(\alpha +\beta )$, and the distribution is valid over $\{\alpha >0,\beta >0\}$. A common choice is the uniform prior $\{\alpha =1,\beta =1\}$. For a reader unfamiliar with Bayesian statistics, the prior can be used to bias the result using prior knowledge (thus the name). However, the actual choice of prior makes little difference with sufficient data. More useful is to reverse this last statement; if the result of an observation depends significantly on the choice of prior, then the data quantity is insufficiently robust. Using the above prior,

We did not need to explicitly calculate the Bayesian denominator $P({\mathrm{PC}}_N=c)$. The denominator is determined by noticing that the binomial distribution combined with a beta distribution is another beta distribution, so the new denominator is the normalization factor from the new beta distribution.

The expression for the moments of $\lambda$ is found using the probability distribution for $\epsilon$ to find the powers of $\lambda$ written as a function of $\epsilon$

A change of variable gives, to within a sign, the formula for the geometric moments of the beta distribution

So, the mean and variance of $\lambda$ are (within a sign) the geometric mean and geometric variance of a beta distribution, respectively. The expected value can be written in terms of the digamma function $\psi$

The polygama function ${\psi }_n$ is a transcendental function available in language math libraries or evaluated with a series approximation.

6.

Optimal Exposure Time

$\lambda$, the expected number of counts in an exposure, depends on the exposure time, while the relevant quantities for an observation are the rate of photons per unit time and the total exposure time available. The obvious choice is to have a very large number of exposures and make the exposure time arbitrarily small. However, in practice, an EMCCD has a noise source called CIC that is added to every exposure,^3,12 where an extra charge is added during clocking of the charge-coupled pixels. CIC has been measured as low as $7\times {10}^{-4}{\mathrm{e}}^{-}/\text{pixel}/\text{readout}$,¹³ while achieving is $5\times {10}^{-3}{\mathrm{e}}^{-}/\text{pixel}/\text{readout}$ is relatively easy.

CIC is not constant across the detector. For a $1\mathrm{k}\times 1\mathrm{k}$-pixel frame transfer device, the number of pixel-to-pixel transfers can vary from $1\mathrm{k}$ to $3\mathrm{k}$ (frame transfer plus row number plus column number) before arriving at the readout and amplification chain (see Fig. 1).

For this section, we define the following variables:

The variable $\eta$ is the fundamental astronomical parameter being measured, but includes dark current. Our expected counts per exposure $\lambda$ is then

The expected photon-counting event detection rate is

To avoid the difficult polygama function, we will find the optimal number of exposures $N$ by maximizing the SNR with respect to $\epsilon$ rather than $\lambda$. This should be acceptable, as CIC is only critical for small $\lambda$. Choose the number of exposures $N$ to be $N\gg \alpha$ and $N\gg \beta$ (or choose $\alpha =\beta =1$). Since $\epsilon$ is beta distributed, we know the mean and variance, and thus the SNR

The SNR evaluated at the expected $\overline{\epsilon }$ is then

Setting the derivative with respect to $N$ to zero gives the transcendental equation

Over the typical region of ${10}^{-3}<{\lambda }_{\mathrm{CIC}}<5\times {10}^{-3}$, ${\lambda }_{\mathrm{opt}}$ is well approximated by the fit

So, for $\mathrm{CIC}\sim 0.001$, the optimal number of well counts per exposure is ${\lambda }_{\mathrm{opt}}\sim 0.045$ counts for best SNR. Given that CIC for current state-of-the-art, EMCCD controllers is in the range of 0.001 to 0.005, $\lambda$ for optimal SNR is restricted to (0.045 to 0.10). Since CIC varies across the detector, ${\lambda }_{\mathrm{opt}}$ can never be more than a compromise. Other considerations might lead to using shorter exposure times than for optimal SNR. For example, bright stars in the field can cause “blooming,” where the wells of the amplification staircase elements overfill and spill charge into adjacent elements. The exposure time may also be limited by the spacecraft stability for applications where spatial resolution is important.

7.

Comparison to Standard Mode

For bright sources, the detector frame rate may not be high enough to achieve the optimal count rate $\lambda$ per exposure. In this case, observing with standard mode (SM) might provide better SNR. The SNR in SM is

Here $\hat{\mu }$ is the readout noise for standard mode, which may be different than the readout noise $\mu$ for PC mode if the EMCCD device has separate readout amplifiers (see Fig. 1). Given the total number of expected events $\eta T$ for an observation and a known ${\lambda }_{\mathrm{CIC}}$, the number of PC-mode exposures is

The expected number of PC counts is

Assuming $N$ is large enough that the prior is irrelevant and replacing $c$ with the expected counts, the PC-mode SNR is

Equations (21) and (24) can be compared to choose which mode is best. A worked example is shown in Fig. 2. In practice, the standard-mode exposure time $T$ is limited by cosmic rays, whereas for PC mode the exposure time $T/N$ is limited by readout speed.

Fig. 2

Assuming ${\lambda }_{\mathrm{CIC}}=0.002{\mathrm{e}}^{-}$ and $\mu ={\hat{\mu }}^{ }=6{\mathrm{e}}^{-}$ for this plot, PC mode always has better SNR than standard mode. The staircase effect for lower counts in PC mode is real and due to rounding up to the nearest $N$ integer number of exposures.

8.

Amplification

The amplification process of EMCCD amplification is similar to the amplification process in an APD.^11,14–18 The amplification of a single electron from Matsuo 1984¹¹ (hereafter “Matsuo”) is given by a recurrence relation for each stage $m$ by

$Q$ is the probability of an electron generating a secondary. The distribution of counts for each stage $P_m(n)$, where $n$ is the number of counts, depends on the distribution of counts at the previous stage. What happens if the initial distribution of counts is not restricted to a single count, or $P_0(n)\ne {\delta }_{1,n}$? Then the recurrence relation spectacularly fails. However, the most important issue with this recurrence relation is that it physically makes no sense; the relation connects the opposite tails of the distribution to generate the next stage. So, although the relation may hold true, it is more of a mathematical curiosity than anything else.

The APD staircase effect can be described as a process where, at each step, there is a chance when transferring the charge that an electron generates a second electron. Now make two assumptions. First, the electrons are not aware of each other, equivalent to the process at each step being linear in the number of electrons generated on average. Second, assume that an original electron can never generate two or more additional electrons, equivalent to assuming that the energy of an electron, once spent on generating a second electron, is gone. With these two assumptions, we can write down a recurrence relation.

Consider, at some step $m$ in the APD staircase, the distribution of electrons $P_m(n)$. How was this distribution generated if the previous step had $k$ electrons? Clearly, $n-k$ electrons each generated a single electron (first assumption), while $k-(n-k)=2k-n$ electrons did not. Since the electrons are independent (second assumption), the applicable distribution is the binomial distribution

This expression looks unwieldy; however, the right-hand side is nothing more than a fixed matrix multiplication to get from one stage of the APD to the next. The matrix is

Once the matrix $B_{nk}$ is formed with a size corresponding to the largest electron counts of interest, the recurrence relation becomes

The final amplified distribution is then

Computing $P_m(n)$ is only a few lines of MATLAB code,

% Create the transform matrix for a single stage of a staircase APD

Q = gain^(1/stages) - 1;

for n = (0:nmax);

 for k = (0:nmax);

  if((k <= n) && (n-k <= k))

   B(k+1,n+1) = exp(gammaln(k+1) - gammaln(n-k+1) - gammaln(2*k-n+1) + (2*k-n)*log(1-Q) + (n-k)*log(Q));

  end % k

 end % n

% Transform initial distribution to the final distribution

P = Po*BB^stages; % Po and P are row vectors here

This recurrence relation gives identical results to Matsuo for an initial distribution of a single electron $P_0(n)={\delta }_{1,n}$, while also, unlike Matsuo, working for any initial distribution.

9.

Comparison to Basden

Basden 2003¹⁹ (hereafter “Basden”) gives the formula for the distribution of final amplified counts $X_A$ as [Basden, Eq. (1)]

$G$ here is the gain of the entire APD staircase, so the gain at stage $m$ would be $G_m=(1+Q{)}^m$.

Basden’s derivation is based on assuming the first assumption above (independent electrons) and a few approximations. We write $\tilde{P}$ because the formula is not a probability distribution; the sum over all values is not equal to unity. The formula has the functional form of an Erlang distribution, but the Erlang distribution is a continuous distribution in $a$. The Erlang distribution is the distribution of the sum of $w$ independent identically distributed (IID) exponentially distributed random variables. This is a hint that the Basden equation could be derived as a sum of IID geometrically distributed random variables. In fact, a plot of the Matsuo distribution for a single initial electron does look very similar to an exponential distribution.

The recurrence relation of Eq. (25) deviates from that of Eq. (26), mostly for larger counts (Fig. 3). Basden explicitly states that their approximation should not be used in this region.

Fig. 3

The Basden approximation is close to the binomial model for moderate counts but fits poorly for the less likely counts. For this plot, the gain is 100, with 50 stages. “Matsuo” is Eq. (25), “Basden” is Eq. (30), and “Binomial” is Eq. (26).

In the limit as the number of stages $m$ becomes large and the amplification $Q$ becomes small for fixed $G$, a closed-form solution in terms of $G$ based on Eq. (26) should exist. At a minimum, the recurrent relation can be written in terms of the eigenvectors and eigenvalues

Since the matrix $B_{nk}$ is triangular, the diagonal elements are the eigenvalues. Taking the limit as $m$ goes to infinity for fixed $G$ gives the eigenvalues as

In the limit, the eigenvectors $V$ are independent of $G$.

10.

Combining Signal with Read Noise

In practice, we add the amplified well counts to the read noise $X_R$. Read noise is also Poisson distributed, so

The expected value of the read noise $\mu$ is trivially measured by taking a zero-exposure-time frame or alternatively running the analog-to-digital converter without advancing the CCD charge. Since the read noise and amplified counts are independent, the distribution of the sum read out by the detector is

Inserting in the expressions from Eqs. (1), (29), and (33)

11.

PC Cutoff Setting

The algorithm used for PC is that $c$ is considered a single count if higher than some threshold $\tau$. The expected value (and second and all higher moments) of this process is then

Since the argument of the summation is a probability distribution over $c$, this is equal to

Substituting in Eq. (35) and swapping the summations gives

The sum over $m$ is again the regularized gamma function from Eq. (3)

The optimal choice for the cutoff parameter $\tau$ is such that the read noise does not bias the PC process. This is equivalent to setting $\tau$ such that false positives equal false negatives. With this choice, lost counts from low-amplification signal counts balances added counts from high noise. The PC method is then an “ideal” photon counter, as in Sec. 5. One can prove that setting false positives to false negatives is equivalent to the read noise having no effect, and so the expected value of $\epsilon$ with noise is equal to $\epsilon$ without noise, or

Let $\hat{\tau }$ be the optimal $\tau$ that makes Eq. (39) satisfied. Then

Over the range of typical values for $\lambda$, $\mu$, and $G$, the optimal $\hat{\tau }$ is well fit to within 5% by

$\hat{\tau }$ only has a weak dependence on $\lambda$ and $G$. The residuals of the fit are due to the $\mu$ dependence not being sufficiently well-fit with a power law. The residuals are shown in Fig. 4.

Fig. 4

The residuals of the Eq. (42) fit to the cutoff $\tau$ in Eq. (41) are within 2%. The scatter is due to variations in $\lambda$ over 0.04 to 0.1 and $G$ over 10 to 1000.

12.

Compression of Photon-Counting Mode Data

For a typical photon-counting exposure, the expected image data are a list of binary values, the vast majority of which are zeros. Since CIC may be a significant contributor to the signal and is unstructured random noise, compression algorithms based on entropy (structure in the image) are likely not optimal. On the other hand, an algorithm that depends on most of the data being zero might be.

If the events in each pixel are independent, then the statistics are Poisson distributed with a “rate” or expected number of pixels between events, of $\lambda$. For a detector where the noise is CIC-dominated and CIC is $\sim 0.001{\mathrm{e}}^{-}/\mathrm{pix}/\text{frame}$, the number of pixels between each photon detection is about $1/(0.001)=1000$.

The simplest compression method for sparse binary data is to count zeros between events. This nominally gives a compression ratio of approximately the number of pixels divided by the number of bits needed to encode the number, or $\sim (1/\lambda ){/\log }_2(1/\lambda )$. A compression ratio roughly on the order of 100:1 should be possible for the above CIC alone.

This algorithm is a variation of run-length encoding (RLE), which is to count the number of bits before a change. The heritage of the RLE algorithm stretches back to at least the early 1960s.²⁰ For mostly empty data, the appropriate standard compression algorithm is the run-length limited (RLL) variation, which has been common for encoding hard disks and optical storage since the 1980s, and only counts zeros. We use $(0,b)$ RLL.

The encoding algorithm for $\mathrm{RLL}(0,b)$ is as follows:

For the data generated by the photon-counting process, we can compute the expected compression ratio. Equation (1) gives the probability for an event (1) and a nonevent (0) for an image pixel. Assuming independent events, the probability of $n$ zeros followed by a nonzero is

The expected number of inputs bits consumed by a $b$-bit encoding is then the probability of all zeros (consuming one encoding value), or some run of zeros followed by a nonzero

Simplifying,

Using the derivative of the formula for the sum of a geometric series, this simplifies to

The final compression ratio is then $E(\text{bits consumed})/b$. The boundaries in $\lambda$ for the ranges of optimal compression are found by looking for the $\lambda$ where compression $(\lambda ,b)$ = compression $(\lambda ,b+1)$. The results are shown in Table 1.

Table 1

Optimal code length for given counts per exposure λ.

13.

Conclusions

The following steps specify the operating mode (PC or SM), expected SNR, optimal exposure time, and data volume. For a space-based astronomical observatory, knowing these parameters is critical to mission design. The steps to set up an EMCCD for an observation are

Fig. 5

An encoding with 5-bit codes has the highest compression ratio for $0.053<\lambda <0.095$, which encompasses most of the range of $\lambda$ for optimal SNR from Sec. 6. The most likely lossless compression assuming optimal SNR is thus from 2:1 to 3:1. The optimal compression ratio for $\lambda <0.01$ is well fit by the approximation $0.264\times {\lambda }^{-0.825}$.

For a given observation of time $T$, the number of exposures $N$ and exposure time $T/N$ must be chosen to ensure that the required SNR is met for all image pixels, taking into account variations in CIC, dark current, and signal brightness across the detector. The optimal photon-counting cutoff $\hat{\tau }$ can be chosen per pixel, either in postprocessing or adaptively as exposures are taken. Using these equations, even bright (but unsaturated) stars will have correctly computed expected values.