Let us assume a two-dimensional pixelated array with infinitely large area and with a total number of pixels N _pix . Incoming events constitute a homogeneous random sequence following Poisson statistics, both in time and in space, with constant average rate per pixel and per unit time n ₀. Let N ₀(t) be the random variable denoting the total number of incoming events in a pixel element over the time interval from 0 to t. During a time frame of duration Δt, the number of incoming events in a pixel element is therefore N ₀ (Δt), and the cumulative value over a number F _N time frames is N ₀(F _N Δt). Due to the physics of the interaction between the impinging particle and the sensor material, a single incoming physical event can trigger the simultaneous detection in multiple, typically neighbors and pixels. To describe the number of pixels firing due to a single incoming event, we define the random variable event multiplicity Mul . Denoting with E ⋅ the mean value operator, the average effective incoming event rate is n = n 0 E Mul .

In the following sections, quantities are described over the time interval of a frame Δt. The generalization to the sum of a series of F _N time frames is straightforward, given the statistical independence between consecutive frames and the homogeneity of the involved processes. For any random variable X , indeed, it holds E X F N Δ t = F N E X Δ t , (3) V X F N Δ t = F N V X Δ t , (4) where V ⋅ indicates the variance operator.

Let us define the random variable M(Δt) as the total number of recorded counts per pixel in a time frame and the average recorded count rate: m = E M Δ t Δ t . (5) Under the working hypothesis of binary outcome–a pixel can register at most one count per time frame–M(Δt) assumes the nature of a Bernoulli distribution (single-trial binomial distribution) and the probability of a successful outcome, i.e., a pixel registers a count that is at least one effective incoming event hits the pixel: P r M Δ t = 1 = ∑ i = 1 + ∞ P r k = i = 1 − P r k = 0 , where P r k = i is the Poisson probability of having i incoming events in a pixel in a time frame. According to our notation, it takes the following form: P r k = i = n 0 E Mul Δ t i e − n 0 E Mul Δ t i ! . It follows that P r M Δ t = 1 = 1 − e − n 0 E Mul Δ t and therefore, E M Δ t = 1 − e − n 0 E Mul Δ t . (6) The expression for the recorded count rate m in Eq. 5 can be then rewritten as m = 1 − e − n 0 E Mul Δ t Δ t . (7)

The derivation of the counts variance for an individual pixel is straightforward from the binomial distribution and yields the following: V M Δ t = P r M Δ t = 1 1 − P r M Δ t = 1 , = 1 − e − n 0 E Mul Δ t e − n 0 E Mul Δ t . The analysis of the behavior of individual pixels can be extended to the collective behavior of the pixel ensemble by introducing an additional random variable representing the spatially averaged total number of recorded counts: M ̄ Δ t = 1 N pix ∑ i = 1 N pix M i Δ t . Given the linearity of the mean value operator, the mean number of spatially averaged counts is simply E M ̄ Δ t = E M Δ t , while for the variance, the derivation is more involved since pixels may exhibit a correlation. In a previous work [14], we derived an expression for the variance of a correlated ensemble of pixels in counting detectors, based solely on the knowledge of the variance of the single pixel, of the first and second moment of the event multiplicity distribution, and of the counting efficiency η c = m n : V M ̄ Δ t = V M Δ t 1 + E Mul 2 E Mul − 1 η c . Expanding it and using our notation, we obtain V M ̄ Δ t = e − n 0 E Mul Δ t 1 − e − n 0 E Mul Δ t 1 + E Mul 2 E Mul − 1 1 − e − n 0 E Mul Δ t n 0 E Mul Δ t . (8) We would like to emphasize that this expression is completely predictable since it does not contain free parameters to be determined, e.g., by fitting. The only empirical term is the multiplicity distribution that can be retrieved experimentally quite easily ² .

The analysis of the statistical properties of events that occurs randomly, independently, and with uniform probability in a generic mathematical space belongs to the category of Poisson point processes, and the gathering of such events into clusters is known as Poisson clumping or burst. Often, problems in this field necessitate of heuristic approaches to come to an explicit conclusion [15, 16]. During the 1980s, a problematic analogous to ours emerged in the field of chromatography, a chemical analysis technique involving the separation of a mixture into its components, their spatial drift through a system—the velocity depending on the nature of the component—, the reaching of a final stationary state, and the identification of the several possibly overlapping and assumed Poisson-distributed signal peaks. The study of the statistical properties of the recorded peaks in the one-dimensional case led Davis and Giddings to the birth of the statistical model of overlap (SMO) or statistical overlap theory (STO) [17]. Unsurprisingly, this first exact analytical result coincides, space exchanged with time, with the dead time model for the classical paralyzable counting detectors reported in Eq. 1. The STO was then extended (with approximations) to the two-dimensional case by Davis in [18], until the work from Roach [15], developed in the 1940s–1960s during studies on biological and hygienic sciences of phenomena closely related to the overlap of spots in two-dimensional beds, was “re-discovered” in [19]. In the same work, it is also demonstrated that this model provides the best approximation to numerical data, among a series of different other models, and therefore, we took it as reference for our work. Incidentally, it is worth mentioning Davis’s extension of “Roach model” to the generic n-dimensional case [20].

To ease the reading, we outline the basic reasoning at the basis of the Roach model, as described in [19]. The theory starts assuming circular “zones” with radius r ₀, whose centers are randomly distributed on a continuous, two-dimensional space. If the centers of two zones are closer than 2r ₀, they are considered overlapping. A series of overlapping zones is called a spot. The recorded number of spots determined with the following scheme. Selected arbitrarily an initial zone A and its first neighbor B, if the distance between the two centers is greater than 2r ₀, zone A is a singlet spot. Otherwise, the overlapping pair forms either a doublet spot or a higher-order multiplet (e.g., triplet, quartet, quintet, etc.). In this case, a third zone C is individuated such as its center is the closest to either A or B. The shortest distance between C and A or C and B is again compared to 2r ₀. If greater, the pair A and B is a doublet spot. If smaller, A, B, and C form at least a triplet spot. The neighbor-searching procedure is then repeated until the distance of all remaining zones to any of the n zones of the spot is greater than 2r ₀. As a result, a spot consisting of n overlapping zones, or an n-tet, has been isolated. The process is iterated for another arbitrary zone until all zones has been assigned. Figure 1 helps visualizing an example of Roach’s zone selection procedure.

The outcome of the neighbors overlapping test can be modeled with the binomial distribution. For instance, probability p ₁ that the distance between a zone center and the center of the nearest neighboring zone is greater than 2r ₀ (i.e., the first zone is a singlet) corresponds to the probability that no event occurs within a circular area of radius 2r ₀ around the center of the first zone—let us call this coincidence area 4A ₀, where A ₀ is the area of an individual zone. Let us remember that the spatial distribution of the events is assumed to follow Poisson statistics. Introducing now our notation and defining k as the number of events falling within the coincidence area in a time frame Δt, the Poisson probability distribution for the number of events can be written as P r k = i = Q E n 0 Δ t 4 A 0 i e − Q E n 0 Δ t 4 A 0 i ! , where QE is the quantum efficiency, defined as the ratio between the number of detected events and the total number of incoming events assuming no coincidence loss [21], which can be also expressed as Q E = 1 − P r Mul = 0 . (9) The introduction of QE arises from the fact that events that pass completely undetected do not contribute to the overall statistics (this statement is reviewed in Section 2.6). From these premises, it follows that p ₁ can be written as p 1 = P r k = 0 = e − Q E n 0 Δ t 4 A 0 . On the other hand, the probability that the distance between a zone center and the center of the nearest neighboring zone is less than 2r ₀ is the complementary probability 1 − p ₁. The chain of events that brings to an n-tet spot therefore consists of n − 1 sequences in which the nearest neighbor distances are less than 2r ₀ and 1 sequence in which the nearest neighbor distance is greater, breaking the spot connection to the remaining zones. If the inter-zone distances are independent from each other—a condition that, as recognized by Roach, is approximately true and necessary to reduce the two-dimensional overlapping problem to a tractable form—the probability p _n for a zone to be part of an n-tet spot is the product of the individual probabilities: p n = e − Q E n 0 Δ t 4 A 0 1 − e − Q E n 0 Δ t 4 A 0 n − 1 . Since each of the QEn ₀Δt incoming events per unit area (corresponding to a pixel, in our case) per time frame has this probability of forming an n-tet spot, the number of zones contributing to the formation of n-tet spots is QEn ₀Δt p _n . However, because n zones are required to form each spot, we need to weight by 1/n this value in order to obtain the expected number P _n of n-tet spots per unit area per time frame. P n = Q E n 0 Δ t p n / n , (10) = Q E n 0 Δ t e − Q E n 0 Δ t 4 A 0 1 − e − Q E n 0 Δ t 4 A 0 n − 1 / n . The total number of spots P is the algebraic sum of all the n-tet spots, which converges analytically to P = ∑ n = 1 ∞ P n , = Q E n 0 Δ t Q E n 0 Δ t 4 A 0 e − Q E n 0 Δ t 4 A 0 1 − e − Q E n 0 Δ t 4 A 0 .

Defining now the random variable E(Δt) as the total number of recorded events per pixel in a time frame, its mean value coincides with p: E E Δ t = P . (11) The recorded event rate e = E E ( Δ t ) / Δ t can therefore be written as e = Q E n 0 2 Δ t 4 A 0 e − Q E n 0 Δ t 4 A 0 1 − e − Q E n 0 Δ t 4 A 0 . (12)

At this point, it is suitable to bring to the attention to the reader that the Roach model was developed in a continuous two-dimensional domain and for the “simple” case of circular-shaped zones ³ , whereas our case involves a discretized domain (pixels) and random event shapes. Nevertheless, results shown in Section 3 suggest that the Roach model still provides an accurate description of the phenomenon, at the price of interpreting A ₀ as a sort of effective quantity to be typically found via curve fitting, on the same line of what is done in [22] for the effective dead time in counting detectors with pulse shapes different from the ideal rectangular one. A rough a priori guess can anyway be attempted, imagining events to be more similar to squares rather than circles—a consequence of the pixelization of the space—with effective area E Mul and side L = E Mul . It is easy to find that this implies a correlation area equal to 4 A 0 = 2 L + 1 2 . (13)

Let us now focus on the variance of the number of recorded events. Due to the complexity of the task in two dimensions, no model could be found neither by the author nor in the consulted literature. However, an analytical solution for the one-dimensional case exists and it deserves some attention. In the context of chromatography, an exact formulation based on statistical arguments was first achieved by Rowe and Davis [23] (see also [24] for a comprehensive summary of one-dimensional models): V P Δ x = λ Δ x e − λ x 0 1 − 2 λ x 0 e − λ x 0 , (14) where P corresponds to the number of recorded events in the one-dimensional space interval Δx, λ corresponds to the incoming event rate per unit space, and x ₀ to the coincidence interval. Few years later, in the context of dead time models for paralyzable counting detectors, Yu and Fessler achieved independently and by analytical methods an equivalent formulation [25]—change Δx with the time interval t and x ₀ with the system dead time τ. In their derivation, the term equivalent to λ e − λ x 0 in Eq. 14 is explicitly identified with the recorded event rate p: p = λ e − λ x 0 , (15) so that Eq. 14 could be rewritten as V P Δ x = Δ x p 1 − 2 x 0 p . (16) By way of analogy, we propose to adapt this expression to the two-dimensional case by identifying p as the recorded event rate e of Eq. 12 integrated over the time frame Δt (in the one-dimensional context, indeed, rate is defined per unit spatial dimension), Δx as the number of pixels N _pix and x ₀ as half the correlation area 4 A 0 2 . This last relation is obtained by equivaleting the expressions for the recorded event rate in one and two dimensions for values of incoming event rates tending to zero, as reported in Appendix 4. Incidentally, it is worth noting that the one-dimensional expression is equivalent to the one for the singlets P ₁ in two dimensions (see Eq. 10 with n = 1), which was, e.g., used to model, alone, the recorded event rate in two dimension in [21]. From the aforementioned proposition, it follows that the variance of the recorded events per pixel and in a time frame can be written as V E Δ t = Δ t e 1 − 4 A 0 Δ t e , (17) = Q E n 0 Δ t 2 4 A 0 e − Q E n 0 Δ t 4 A 0 1 − e − Q E n 0 Δ t 4 A 0 1 − Q E n 0 Δ t 4 A 0 2 e − Q E n 0 Δ t 4 A 0 1 − e − Q E n 0 Δ t 4 A 0 .

The knowledge of the first two momenta of the recorded statistics allows us to compute the corresponding SNR over a time frame, defined as S N R X Δ t = E X Δ t 2 V X Δ t , which, thanks to the properties highlighted in Eqs 3, 4, can be easily generalized to the sum of a series of time frames: S N R X F N Δ t = F N E X Δ t 2 V X Δ t .

More interesting, however, due to its importance and widespread use in the imaging community (see [26] and references therein for an interesting historical overview), is the DQE, defined as D Q E = S N R OUT S N R I N , (18) where SNR_OUT is the signal-to-noise ratio referring to the recorded signal statistics and SNR_IN is the one referring to the incoming signal statistics, which, under the assumption of obeying Poisson statistics, equals to S N R I N F N Δ t = F N n 0 Δ t , (19) which makes the DQE independent on the number of acquired frames.

In its modern acceptance ⁴ , DQE is conceived in the two-dimensional domain of the spatial frequencies by applying Fourier transform to the output signal and noise. Use of frequency analysis requires the system to exhibit the properties of shift-invariance, linearity and wide-sense stationary statistics. In general, rarely, a system satisfies all the requirements in a rigorous manner. For example, a pixelated detector is not strictly shift-invariant, unless shifts are by an integer number of pixels. In addition, in a counting system, the noise is stationary only in conditions of uniform illumination as it depends on the incoming signal itself. Therefore, care must be taken when interpreting the results as they might be approximations of the true properties of the system [27, 28]. In our specific case, the requirement of linearity is clearly not satisfied. Since a straightforward extension of the theory is not well-established (see [29] for a nice compendium on possible generalizations of the concept of DQE to non-linear systems), we limit our analysis to the degenerate zero-frequency case denoted DQE₀ ⁵ .

The DQE₀ of the recorded counts as a function of the incoming event rate can therefore be retrieved using Eqs 6, 8, 18, 19 in obtaining D Q E 0 M ̄ = 1 − e − n 0 E Mul Δ t n 0 Δ t e − n 0 E Mul Δ t 1 + E Mul 2 E Mul − 1 1 − e − n 0 E Mul Δ t n 0 E Mul Δ t , while one of the recorded events using Eqs 11, 17–19 obtaining D Q E 0 E = Q E 2 n 0 Δ t 4 A 0 e − Q E n 0 Δ t 4 A 0 1 − e − Q E n 0 Δ t 4 A 0 1 − Q E n 0 Δ t 4 A 0 2 e − Q E n 0 Δ t 4 A 0 1 − e − Q E n 0 Δ t 4 A 0 .

Additionally, derived quantities of particular interest, which we would like to mention, are the area occupancy (AO) and the coincidence loss fraction (CLF). The AO is defined as the ratio between the average number of counting pixels per frame and the total number of pixels, coinciding with the average number of recorded counts in a time frame. A O = N pix E M Δ t N pix , = E M Δ t . Information on the AO can have implications on the detector readout mechanism design and optimization, as well as on the choice of suitable data compression algorithms. The coincidence loss fraction is defined as the ratio between the number of “lost” events, i.e., not recorded, and the total number of incoming events. C L F = Q E n 0 − e Q E n 0 , (20) = 1 − Q E n 0 Δ t 4 A 0 e − Q E n 0 Δ t 4 A 0 1 − e − Q E n 0 Δ t 4 A 0 and reflects the counting efficiency of the system. In Eq. 20, the physical incoming event rate n ₀ in the denominator is scaled by QE as QEn ₀ is the actual ideal event rate recorded in the absence of coincidence loss.

i. The study of the counting statistics of pixels featuring binary counts presented here complements a picture already including at least the classical paralyzable and non-paralyzable counting modes [25] and a particular case of the non-paralyzable mode [14], where the paralysis is avoided, thanks to a circuital stratagem called “instant retrigger.”

To validate the analytical models, we used results of numerical simulations carried out with an improved version of the Monte Carlo tool used in [13, 30, 31]. The first step of the workflow was the creation of a statistically relevant pool of electron tracks in the semiconductor sensor. A total of 40 million tracks were generated using FLUKA ⁷ —a Monte Carlo particle transport and interaction suite [32, 33]—storing, for each of them, the three-dimensional spatial coordinates with an accuracy of 1 µm and the amount of energy released therein.

We then processed each individual electron track with a custom-developed numerical code mimicking the physics of the charge collection and signal formation at the pixelated electrode. The generated charge distribution of each track segment was thus propagated through the remaining sensor thickness to the pixels, and a Gaussian blurring was added to reproduce the effect of thermal diffusion, with a total width depending on an initial intrinsic contribution and, under the assumption of the constant electric field, to a contribution depending on the total travel length. The charge collected by each pixel was converted into energy and a counting threshold was applied, if the energy is higher, the pixel counts 1; otherwise, it counts 0. A random fluctuation representing the electronic noise, normally distributed and assumed uncorrelated among the pixels, was added to the signal. The response of both the sensor and the readout electronics has been assumed uniform in space. At this point, it was possible to extract the statistical distribution of the event multiplicity. Then, for each value of a series of incoming electron rates, a set of 2000 independent images (frames) was generated. The total number of impinging electrons for each individual frame was chosen randomly according to the suitable Poisson statistics and then they were uniformly and randomly distributed across the sensor surface, which covered an area of 2048 × 2048 pixels. A pixel cluster recognition algorithm was then applied to each frame to isolate single events. In order to be considered isolated, two clusters need to be separated by at least one empty pixel.

We chose a case study realistic for CryoEM applications consisting of a 25-µm-thick MAPS detector, covered with 5 µm of Al accounting for the metal layers and featuring square pixels of size 10 µm. The counting threshold energy was 1 keV, the electronic noise was 200 eV rms, and the thermal diffusion between 1 and 3 µm rms. The frame rate was assumed 1 kfps. The energy of the impinging electrons was 300 keV, with values of incoming rates in the range 1–200 el/s/pix.