The statistics of electron-hole avalanches

Charge multiplication through avalanche processes is commonly employed in the detection of single photons or charged particles in high-energy physics and beyond. In this report, we provide a detailed discussion of the properties of avalanches driven by two species of charge carriers, e.g. electrons and holes in a semiconductor exposed to an electric field. We derive equations that describe the general case of avalanches developing in inhomogeneous electric fields and give their analytical solutions for constant fields. We discuss consequences for the time resolution achievable with detectors that operate above the breakdown limit, e.g. single-photon avalanche diodes (SPADs) and silicon photomultipliers (SiPMs). Our results also describe avalanches that achieve finite gain and are important for avalanche photodiodes (APDs) and low-gain avalanche detectors (LGADs).


Introduction
Ever since the discovery of electron ionisation processes in gases by Townsend [1], avalanche processes have been used as powerful mechanisms to provide high gain in situations where a small amount of initial charge needs to be detected. In recent decades, advances in the manufacturing of semiconductors have made it possible to exploit avalanche multiplication occurring in the solid state. This has led to the development of different types of semiconductor devices that make use of avalanche charge multiplication in various ways. Detectors based on this principle include avalanche photodiodes (APDs), low-gain avalanche detectors (LGADs), single-photon avalanche detectors (SPADs), and silicon photomultipliers (SiPMs).
Contrary to avalanches developing in gases, where solely electrons drive the avalanche, charge carriers of both polarities (i.e. electrons and holes) can be important in semiconductors. Many properties of electron-hole avalanches are already well understood, including the condition for breakdown [2,3], fluctuations of the avalanche gain [4] and some aspects concerning the average development of the avalanche [5,6]. In this report, we provide a detailed discussion of the properties of electron-hole avalanches that complements the available literature. Figure 1a shows the formulation of the problem. A number of electrons and holes are deposited at x = x 0 in a solid state material of thickness d exposed to a strong electric field E(x). This configuration can be created in a reverse-biased p-n junction with suitable dopant concentrations. Both types of charge carriers then multiply in the (generally position-dependent) electric field, which leads to the formation of an avalanche and the multiplication of the initial charge.

Description of the problem
The total number of charge carriers in the avalanche determines the current induced on the readout electrodes of a particle detector. Detectors such as SPADs and SiPMs operate above the breakdown limit, where the initial charge triggers a diverging avalanche with a certain efficiency. A threshold is applied to the signal, as shown in Figure 1b. The threshold-crossing time fluctuates due to the stochasticity of the avalanche, and these variations determine the time resolution of the device. In the applied electric field, electrons move in the positive x-direction and holes move in the negative x-direction. b) Signal produced by a series of electron-hole avalanches in a thin layer of silicon, obtained from a Monte Carlo simulation. The total number of charges in the avalanche determines the current I(t) induced on the readout electrodes. The time at which the signal crosses the applied threshold is indicated by a vertical dashed line for each avalanche. One avalanche fails to diverge and does not cross the threshold.
We model the development of the avalanche as shown in Figure 2. Both electrons and holes can undergo multiplication reactions and produce additional electron-hole pairs. The probability for an electron to multiply in the small spatial interval dx is taken to be α(x)dx. Analogously, the probability for a hole to multiply is given by β(x)dx. We assume that electrons and holes move with velocities v e (x) Figure 2: Multiplication processes occurring in the material. Electrons (e) move with a velocity ve(x) and can create additional electron-hole pairs as described by the impact ionisation coefficient α(x). Holes (h) move with a velocity v h (x) and can create additional electron-hole pairs as described by the impact ionisation coefficient β(x). The rate with which these processes occur depends on the strength of the local electric field E(x), and thus on the position x.
and v h (x) through the material. These velocities depend on the strength of the local electric field E(x) in the medium, and therefore on the position x. The probabilities for a multiplication to occur in a small time interval dt are then α(x)v e (x)dt and β(x)v h (x)dt. The impact ionisation coefficients α and β encode the details of the underlying scattering processes in the material, and as such are also strongly dependent on the applied electric field E(x). In analogy to avalanches in gases, they are sometimes referred to as "Townsend coefficients". Together with the velocities v e and v h , we treat the Townsend coefficients here as external input parameters and refer to parametrisations for specific materials that are available in the literature.

Outline
This report consists of three main parts. Section 2 discusses the development of an electron-hole avalanche in an infinitely extended semiconductor to which a homogeneous electric field is applied. This allows us to develop the basic characteristics of electron-hole avalanches in isolation. In Section 3, we then return to the practical situation mentioned above and in Figure 1a. For both scenarios, we give a complete treatment of the average development of the avalanche as well as its fluctuations around this average evolution. Taken together, this allows us to quantify the achievable time resolution in the sense of Figure 1b. All of our results are phrased in terms of explicit analytic formulas, or, in situations of greater complexity, efficient numerical procedures.
Our discussions have as their starting point a series of differential equations that encode the dynamics of the avalanche process. These equations are most easily derived from the principles of stochastic mechanics, to which we give a self-contained introduction, along with a list of references, in Appendix A. However, the material presented there is purely supplementary and our results, and all calculations presented in the main body, can be appreciated without consulting this appendix.

Electron-hole avalanche in an infinite semiconductor with a constant electric field
We first consider the idealised scenario of an infinitely extended semiconductor exposed to a constant electric field E, as shown in Figure 3. Then, the Townsend coefficients α and β and the velocities v e and v h become independent of the position x. a) b) Figure 3: a) A one-dimensional semiconductor of infinite spatial extent, exposed to a constant electric field E. Electrons (e) move to the right with a constant velocity ve and holes (h) move to the left with a constant velocity v h . b) Multiplication processes occur in the material at a rate independent of the position x.
We are interested in describing the evolution of the total number of electrons and holes in the medium, N e and N h , as the avalanche progresses. The formulation of the avalanche process given in Section 1.1 introduces two important assumptions that are worthwhile to make explicit. First, we take all charge carriers to be independent of each other, irrespective of the total number of charges present in the avalanche. Second, the multiplication process itself is assumed to be memoryless, i.e. the probability for a charge carrier to undergo a multiplication reaction is independent of its history.
The formation of the avalanche is a stochastic process, i.e. N e and N h at a given time t are understood as random variables. We denote the probability to find an avalanche containing N e electrons and N h at time t, starting from certain initial conditions, as p(N e , N h , t). This probability satisfies the following differential equation, which determines its time evolution, The content of this equation can be appreciated intuitively. There are two ways in which the probability p(N e , N h , t) can be modified in a small time interval dt. First, any charge carrier in an avalanche with N e −1 electrons and N h −1 holes can multiply and add an additional electron-hole pair, thereby increasing p(N e , N h , t). The probability (per unit time) for this to happen is proportional to p(N e − 1, N h − 1, t) and given by the first, positive, term on the right-hand side. Second, any multiplication in an avalanche with N e electrons and N h holes reduces p(N e , N h , t) by an amount proportional to itself. It is given by the second, negative, term. In the physics and statistics literature, Equation 1 is sometimes referred to as a "master equation" for a Markov process. Its structure is a direct consequence of our independence and memorylessness assumptions. We give a formal derivation in Appendix A.1.6.

Exact expression for
To solve Equation 1, we consider an avalanche that is initiated by N 0 e electrons and N 0 h holes at t = 0, i.e. p(N e , N h , t = 0) = δ Ne,N 0 e δ N h ,N 0 h . The master equation only couples the probabilities p(N e , N h , t) and p(N e + 1, N h + 1, t), i.e. it relates configurations that differ only by some number of electron-hole pairs. As a first step, it thus makes sense to consider the evolution of the number of additional electron-hole pairs created during the evolution of the avalanche. We label this quantity by the random variable K eh . Then, the probability p(N e , N h , t) can be expressed in terms of K eh as where N e (K eh ) = N 0 e + K eh and N h (K eh ) = N 0 h + K eh are the numbers of electrons and holes present in the avalanche after the creation of K eh additional pairs. Inserting this relation into Equation 1 yields the following evolution equation for p(K eh , t), We now express this equation in terms of the z-transform of p(K eh , t), which is defined as In the z-domain, Equation 3 becomes Since K eh is defined as the number of electron-hole pairs created during the evolution of the avalanche, the initial condition for Equation 4 is p(K eh , t = 0) = δ K eh ,0 . In the z-domain, this becomes P (z, t = 0) = 1. The evolution equation can be solved with the method of characteristics and the solution compatible with the above initial condition is with The parameter A characterises the avalanche, and as we shall see below, many of its properties can be expressed in terms of it. We refer to A as the "avalanche parameter". The quantity ν(t) sets the absolute scale for the (exponential) growth of the avalanche.
To perform the inverse z-transformation of Equation 5, we make use of the relation where Γ(·) is the gamma function. This gives where the approximation holds in the limit of late times, i.e. large avalanches. The avalanche parameter A effectively controls the shape of the distribution, as illustrated in Figure 4. Via Equation 2, this immediately yields the solution for p(N e , N h , t). For many practical applications, the total number N of charge carriers is relevant, N = N e + N h , e.g. for the computation of the signal induced on the readout electrodes in a particle detector.
where K eh (N ) = (N − N 0 )/2 is the number of additional electron-hole pairs that are required to reach N charges from N 0 = N 0 e + N 0 h . Inserting the result obtained for p(K eh , t), we find  Figure 4: Probability p(K eh , t) to find K eh additional electron-hole pairs at time t in an avalanche triggered by a single electron. The probability p(K eh , t) is defined only for integer K eh , but for visual clarity is drawn here as a continuous function. The velocities of electrons and holes are taken to be equal, ve = v h = v. The dashed grey line corresponds to β = 0, i.e. an avalanche in which only electrons multiply. The solid black line shows an electron-hole avalanche with identical α + β, but β = 0, i.e. a non-integer value for the avalanche parameter A.

Moments
The probability p(K eh , t) determines the m-th moment K m eh , which is defined as In particular the first two moments K eh and K 2 eh are important, since they contain information about the average evolution of the avalanche as well as its fluctuations around this average. Through the identities N = N 0 + 2K eh , N e = N 0 e + K eh and N h = N 0 h + K eh , the moments K m eh also give access to the expectations N m e N n h and N m . We just note here the mean and the variance of K eh , which are the most important ones in practice, Summing Equation 8 explicitly for large m can be difficult. An alternative way to obtain the moments N m e N n h directly is to consider them as explicit functions of time. As is shown in Appendix A.1.7, the moment N m e N n h satisfies the following differential equation, The sum runs over the set S = {(l, k)|0 ≤ l ≤ m, 0 ≤ k ≤ n} \ (m, n), i.e. contains all combinations except the one where l = m and k = n simultaneously. An inspection of the right-hand side of Equation 9 shows that it contains only terms N l e N k h where l + k ≤ m + n. Iterating this expression thus gives rise to a closed system of differential equations that can be solved systematically. For the first moments N e and N h , these equations read while for the three second moments N 2 e , N 2 h , and N e N h they are As a convenient way to quantify the fluctuations of the avalanche around its average evolution, we introduce the ratio of the variance of N and the square of its mean, σ(N ) 2 / N 2 . We will refer to this quantity as measuring the "relative fluctuations" of N . Similar measures can be defined for other observables such as N e or N h , and we shall make extensive use of them in what follows. The solutions of Equations 11-13 then reveal the following relation between the avalanche parameter A and the relative fluctuations of the number of charges at late times The solutions also imply that N e and N h are maximally Pearson correlated at all times. We have

Avalanches driven by a single species: the case of integer-valued avalanche parameter
The results given above have important special cases. In case the avalanche parameter A is an integer, the dependence on K eh in Equation 6 simplifies. We get This covers the case of avalanches that are driven by electrons only (β = 0), where we have A = N 0 e and where N e = N 0 e e αvet . This is the well-known Yule-Furry law [7,8] for electron avalanches. Another important special case concerns avalanches that are initiated by a certain number of electron hole pairs, i.e. N 0 e = N 0 h = N 0 eh . In this case, the avalanche parameter becomes A = N 0 eh . The avalanche is then conveniently described in terms of the total number of electron-hole pairs, N eh = N 0 eh + K eh . The statistics of these avalanches is identical to those driven by electrons only, where N eh = N 0 eh e (αve+βv h )t .

Time response function
The current I(t) induced on the readout electrodes in a particle detector can be computed with the Ramo-Shockley theorem [9,10]. If we assume a homogeneous weighting field E w /V w throughout the semiconductor, the signal produced by the avalanche is where e 0 is the elementary charge. We define t th as the time at which the signal crosses an applied threshold of I th . The threshold-crossing time is a random variable whose distribution contains information about the achievable time resolution. The signal in Equation 18 can also be expressed as Making use of this relation, a threshold I th applied to the signal can be converted into an equivalent threshold N th on the total number of charges in the avalanche (which is not generally experimentally accessible). It is thus sufficient to study the statistics of the threshold-crossing time for this case. We define the time response function ρ(N th , t)dt to give the probability that the avalanche crosses a threshold of N th charges in the small time interval [t, t + dt]. It is given by where the additional factor [αv e N e (N th ) + βv h N h (N th )] dt gives the probability that an avalanche with N th charges crosses this threshold in the time dt through the multiplication of an electron or hole. The time response function is properly normalised such that where it is again understood that (N th − N 0 )/2 ∈ N. For large thresholds (and late times), the time response function is approximately given by In this form we can see that a scaling of the threshold N th → N th λ corresponds to a shift of the time argument by t → t − ∆t, with ∆t = log λ αve+βv h . Figure 5 indeed shows that the shape of the time response function changes significantly only for small N th . In the regime of large N th , only a shift in the time argument results.

Time resolution from time response function
Using the time response function from Section 2.4, we can compute the expected threshold-crossing time t th and its second moment t 2 th can be computed as where the polygamma function ψ k (z) is defined in terms of derivatives of the gamma function as ψ k (z) = d k+1 ln Γ(z)/dz k+1 .  can be obtained. We use this quantity as a measure for the intrinsic time resolution that can be achieved by a device that relies on the avalanche multiplication of some initially present charge. We find The function ψ 1 (z) is monotonically decreasing with z. For large thresholds the second term thus becomes negligible compared to the first, and the time resolution σ th saturates at σ th,∞ = lim N th →∞ σ th (N th ). We get The rapid convergence σ th (N th ) → σ th,∞ , even for moderate thresholds, is shown in Figure 6: the time resolution is essentially identical to its limiting value already for thresholds N th of 30-40. For the special case of an electron avalanche (i.e. β = 0) initiated by N 0 e electrons, we have where the approximation is valid for large N 0 e . This shows that the time resolution improves as σ th ∼ 1/ N 0 e as the primary charge increases. In case the avalanche is initiated by a single electron, N 0 e = 1, we recover the result which was derived in [11] in the context of resistive plate chambers.

Time resolution for large thresholds
For large thresholds, N th → ∞, the time resolution can be computed in a physically more intuitive way that does not require the full time response function ρ(N th , t). At late times, i.e. for large avalanches,  we expect the evolution of the avalanche to follow Note that this is a relation between the two random variables K eh (t) and k ∞ . The latter implements the fluctuations in the absolute size of the avalanche. The exponential scaling is factored out and made manifest. Whenever K eh N 0 , as is the case in the limit we consider here, a similar relation also holds for the total charge content of the avalanche, N .
The random variable describing the threshold-crossing time t th is then given in terms of k ∞ by and its standard deviation σ th,∞ is This shows that the asymptotic time resolution is determined by the fluctuations of the logarithm of the number of charges in the avalanche. The argument as presented above is very general and holds for any multiplication process that leads to (asymptotically) exponential growth. We now make the connection to the specific avalanche process described above in Section 2.1. At late times, the distribution p(K eh , t) in Equation 6 can be written as We see explicitly that the distribution has support only for exponentially growing avalanches of the form of Equation 25. The distribution of the factor k ∞ can be directly read off: it follows a gamma distribution with a shape parameter of A and a rate parameter of unity. As is shown by an explicit calculation in Appendix A.1.8, we have that lim and so we recover the expression for the time resolution already stated in Equation 24.

Electron-hole avalanche in a thin semiconductor with an arbitrary electric field
We now return to the discussion of the development of an electron-hole avalanche in a particle detector of the type described in Section 1. We consider a bounded one-dimensional region 0 ≤ x ≤ d, in which the avalanche develops. This "avalanche region" is depicted in Figure 7a. Practically feasible doping profiles generally lead to a position-dependent electric field E(x) in the semiconductor upon the application of the bias voltage. To take this into account, we allow the Townsend coefficients and velocities to depend on the position, i.e. we assume general profiles α(x), β(x), v e (x) and v h (x). The probability for a charge multiplication reaction to occur thus also depends on the position, as shown in Figure 7b.
The avalanche develops in a semiconductor of thickness d that is exposed to a position-dependent electric field E(x). We take electrons and holes to move with velocities ve(x) and v h (x). Electrons can leave the avalanche region through the boundary at x = d; holes cross the boundary at x = 0. b) The probability for charge multiplication to occur now depends on the position x.
As explained in detail in Appendix A.2, this scenario presents significant technical challenges, which make the computation of explicit distributions such as p(N, t) very difficult. In the following, we will thus focus on the discussion of expectation values of the spatial charge densities n e (x, t) and n h (x, t). These are now random variables indexed by x and t. This will allow us to gain insight into the average spatial distribution of charges in the semiconductor, as well as their correlations across different positions. In all cases, we shall consider the evolution of the avalanche starting from well-defined initial charge densities at t = 0, n 0 e (x) and n 0 h (x).

Average development of the avalanche
We begin with the discussion of the average development of the avalanche. The quantities of interest here are the expectation values n e (x, t) and n h (x, t) , i.e. the average densities of electrons and holes at a specific position x. For brevity, we shall often suppress the explicit position-and time dependence.
As is shown in Appendix A.2.6, these averages satisfy the following system of equations, These equations are the position-dependent analogue of Equation 10: they immediately reduce to the latter for constant Townsend coefficients and velocities upon integration over x and the assumption that the average densities vanish as x → ±∞. They can also be interpreted as continuity equations for the average electron-and hole densities in the avalanche region. The situation in Figure 7a demands the boundary conditions n e (0, t) = 0 and n h (d, t) = 0 as well as the initial conditions n e (x, 0) = n 0 e (x) and n h (x, 0) = n 0 h (x). For general field profiles E(x), Equations 28-29 can be solved numerically. Below, we present their analytical solutions for constant electric fields E, i.e. constant Townsend coefficients and velocities.

Solution for constant electric field
In this case, Equations 28-29 are equivalent to the following second order differential equation for n e , together with the algebraic relation To find the most general solution of Equation 30, we make the ansatz n e = X(x)T (t). Inserting this into the equation and acting on it with ∂ x ∂ t shows that either X /X or T /T must be constant. The former is incompatible with our boundary conditions. Focusing on the second case, we set T = λ t T , i.e. T (t) = Ce λtt with some constant C. We then get With the following definitions, we can write the solution of Equation 31 as The "eigenvalues" λ for which these solutions satisfy all boundary conditions must satisfy the equation which depends only on the characteristic combination αβd 2 . We will refer to the corresponding solutions f e λ and f h λ as "eigenfunctions". The eigenvalue equation has a finite number of real-valued solutions, as well as an infinite number of complex-valued solutions, which appear in complex-conjugate pairs. The latter can be efficiently obtained by solving Equation 35 numerically, using as initial values. The number of real-valued solution depends on the value of αβd 2 . Their behaviour is shown in Figure 8. For the case αβd 2 = 3, Figure   Inner product and orthogonality of eigenfunctions. We can construct an inner product that makes the eigenfunctions in Equations 33 and 34 orthogonal for different eigenvalues. For arbitrary functions f e (x), As is shown in Appendix B, the eigenfunctions for different eigenvalues are orthogonal (but not orthonormal) w.r.t. this inner product, The normalisation N (λ) is given by General solution of the initial value problem. The most general solution of the initial value problem posed by Equations 28-29 is given by a linear combination of the eigenfunctions, The coefficients C(λ) are determined by the initial condition. Making use of the orthogonality property in Equation 38, we can write them explicitly as If we start with N 0 e electrons and N 0 h holes at the position x = x 0 , i.e. with n 0 A similar formula for the coefficient C(λ 1 ) is mentioned in [12] for the special case α = β and v e = v h . The eigenvalue with the largest real component, labelled as λ 1 , is of particular significance for the behaviour of the avalanche at late times. If γ(λ 1 ) = γ 1 > 0 and C(λ 1 ) = 0, the average grows exponentially on average. The device in which the avalanche forms is then said to operate above breakdown. Below breakdown, i.e. for γ 1 < 0, the total charge created by the avalanche (and therefore the gain) is finite. The conditions under which breakdown occurs are discussed in [2,3].
Induced signal and late-time behaviour. If we assume a constant weighting field throughout the avalanche region, Equation 18 determines the average current induced on the readout electrodes of a particle detector. It is given by The average total numbers of electrons and holes present in the avalanche region at time t, N e (t) and N h (t) , are given by The coefficients N e (λ) and N h (λ) are For a device operated above breakdown, N (λ 1 ) = N e (λ 1 ) + N h (λ 1 ) is related to the normalisation of the fastest-growing component that increases as e γ1v * t and becomes dominant at late times. Figure 10 shows a comparison with results obtained from Monte Carlo (MC) simulations. Already after one transit time (T transit = d/v * ), this approximation describes the average charge content with excellent precision.
Below breakdown, N (λ 1 ) can be used to compute the absolute size of the induced (decaying) signal at late times, as is shown in Figure 11. This is relevant for detectors such as avalanche photodiodes. In this regime, the predominant fraction of the signal is generated within one transit time T transit . For t > T transit , the average signal decays approximately exponentially with a time constant given by

Development of spatial correlations and fluctuations
The spatial development of an electron-hole avalanche proceeds in a very different manner compared to an avalanche that is driven by electrons only. Since electrons and holes move in opposite directions, we can no longer expect the charge densities at different positions to be independent of each other (as is the case for an electron-only avalanche), but should rather anticipate nontrivial correlations. In turn, these will affect the fluctuations of the charge content of the avalanche, and thus of the induced signal.
To study these, we look at the expectation values n e (x)n e (y) , n h (x)n h (y) and n e (x)n h (y) . As is shown in Appendix A.2.6, they satisfy the following system of equations with the boundary conditions n e (0, t)n e (y, t) = n h (d, t)n h (y, t) = 0 and n e (0, t)n h (y, t) = n e (x, t)n h (d, t) = 0 as well as the initial conditions n e (x, 0)n e (y, 0) = n 0 e (x)n 0 e (y), n h (x, 0)n h (y, 0) = n 0 h (x)n 0 h (y) and n e (x, 0)n h (y, 0) = n 0 e (x)n 0 h (y). These equations are the direct position-dependent analogues of the evolution equations for the moments N m e N n h encountered earlier: taking the Townsend coefficients and velocities to be constant and integrating over x and y immediately reproduces  We are interested in computing the covariances cov[n e (x), n e (y)] = n e (x)n e (y) − n e (x) n e (y) , which encode information that is not already contained in the averages n e and n h . The diagonal components (x = y) of cov[n e (x), n e (y)] and cov[n h (x), n h (y)] correspond to the variances of the charge densities, i.e. measure the magnitude of fluctuations at a certain location. Their off-diagonal entries (x = y) contain information about the correlations of the charge densities at different positions. Furthermore, we expect cov[n e (x), n h (y)] to be trivially nonzero if x > y: electron-hole pairs that are created at a common position at a time t c (and are thus positively correlated) move apart such that the electrons are always to the right of the holes for t > t c . We thus make the ansatz where the three functions f ee 1 , f hh 1 and f eh 1 make the expected features manifest. The remaining functions f ee 2 , f hh 2 and f eh 2 encapsulate additional, more indirect, effects. The functions f ee 2 and f hh 2 are symmetric in x and y, while f eh 2 does not have any symmetry property. Since n e and n h are spatial densities, the functions f ee 2 and f hh 2 have dimensions of L −2 , while f ee 1 and f hh 1 have dimensions of L −1 . We have seen above that densities scale as e γ1v * t for large avalanches. Purely on dimensional grounds, we therefore expect f ee 2 and f hh 2 to dominate at late times, while f ee 1 and f hh 1 are important at early times. As an explicit calculation shows, the same is in fact true for f eh 1 and f eh 2 .
As a measure of the fluctuations of the avalanche around its average evolution, we consider again the ratios σ(n e (x)) 2 / n e (x) 2 and σ(n h (x)) 2 / n h (x) 2 . These are the position-dependent versions of the relative fluctuations from Equation 14. Their limiting behaviour for large avalanches, i.e. late times, will play an important role for the discussion of the time resolution in Section 3.3 below. Making use of the scaling laws mentioned above, the relative fluctuations can be expressed entirely in terms of f 2 in this limit, σ(n e (x)) 2 To make progress towards this goal and to find evolution equations for the functions f 1 and f 2 , we insert the above ansatz into the original Equations 44-46. The functions f ee 1 , f hh 1 follow the differential equations with the boundary conditions f ee 1 (0, t) = f hh 1 (d, t) = 0 and the initial conditions f ee 1 (x, 0) = f hh 2 (x, 0) = 0. In the domain x > y, the function f eh 1 is determined by with the boundary condition and the initial condition f eh 1 (x, y, 0) = 0. Finally, the functions f ee 2 , f hh 2 and f eh 2 are governed by the following equations, with the boundary conditions f ee 2 (x, 0, t) = f ee 2 (0, y, t) = 0, f hh 2 (x, d, t) = f hh 2 (d, y, t) = 0, f eh 2 (0, y, t) = f eh 2 (x, d, t) = 0 and the initial conditions f ee 2 (x, y, 0) = f hh 2 (x, y, 0) = f eh 2 (x, y, 0) = 0. This is a system of three coupled differential equations, for which f eh 1 plays the role of an external source. Equations 49-55 are equivalent to the original system in Equations 44-46, but all participating functions and source terms are now continuous. For general field distributions E(x), they must be solved numerically. Below, we discuss the characteristics of their solutions for constant electric fields E.

Solution for constant electric field
For position-independent Townsend coefficients and velocities, Equations 49 and 50 can be solved directly, and their solutions are and The retarded time t ret and the positionx are defined as For constant Townsend coefficients and velocities, the source term entering into Equations 53 and 54 for f ee 2 and f hh 2 can be expressed in terms of the function It is the symmetric continuation of the function f eh 1 into the region x < y. The system of equations determining the functions f ee 2 , f hh 2 and f eh 2 deserves a bit more thought. First, we introduce the auxiliary function f he 2 (x, y, t) = f eh 2 (y, x, t). Then, this system of equations is of the form ∂ ∂t where F 2 is the vector of functions F 2 = (f ee 2 , f hh 2 , f eh 2 , f he 2 ) T and g is the source term built from f eh 1,sym , g = (βv h f eh 1,sym , αv e f eh 1,sym , 0, 0) T . L x,y is a linear operator that contains only spatial derivatives. If we defineF t 2 (x, y, t) as the solution of the following homogeneous initial value problem, then the solution of the original system in Equation 59 (and therefore the solution for the functions f 2 ) is given by It is easy to see that this is indeed a solution of the inhomogeneous system of equations. Acting on this expression with the operator ∂ t + L x,y yields, by Equation 60 Solution of the homogeneous system. To be able to evaluate the integral in Equation 61, we first need to solve the homogeneous system in Equation 60. Written explicitly, it reads Inserting the ansatz we find that the two functions f e 2 and f h 2 must follow the equations, These are identical to the equations for the expected densities n e and n h studied in Section 3.1. Making use of the results found there, we can express the general solution for f ee 2 , f hh 2 and f eh 2 again as a linear combination of eigenfunctions. They are now indexed by two eigenvalues, λ a and λ b . We define F λa,λ b = (f ee λa,λ b , f hh λa,λ b , f eh λa,λ b , f he λa,λ b ) T as the vector combining the eigenfunctions for the functions contained in F 2 . Its entries are given in terms of the eigenfunctions f e λ and f h λ defined in Equations 33 and 34, The most general solution of the homogeneous system is then a linear combination of these eigenfunctions, where the coefficientsC(λ 1 , λ 2 ) are determined by the initial condition. We derive an explicit expression below for the problem in Equation 60.
Inner product and orthogonality of eigenfunctions. Analogously to before, we again define a convenient inner product. For two arbitrary vectors of functions F = (f ee , f hh , f eh , f he ) T and G = (g ee , g hh , g eh , g he ) T , where each component function is defined on ( Following this definition, the inner product of two eigenfunctions F λa,λ b and F λ a ,λ b for eigenvalues λ a , λ b and λ a , λ b can be expressed in terms of the inner product defined in Equation 37 as, As before, this inner product makes eigenfunctions corresponding to different eigenvalues orthogonal to each other. The normalisation N (λ, t) is the one defined in Equation 39.
Solution of the inhomogeneous system. Making use of the orthogonality of the eigenfunctions, the coef-ficientsC(λ a , λ b , t ) for the initial value problem in Equation 60 are thus given bỹ They combine with the eigenfunctions F λa,λ b to produce the solutioñ

The full solution given by Equation 61 then reads
We see that the solution of the inhomogeneous system has the same character as the solution of the homogeneous system, except that the coefficients multiplying the eigenfunctions now explicitly depend on time.
Late-time behaviour. The behaviour of the solution F 2 as t → ∞ is of particular interest, as it contains information about the asymptotic level of fluctuations of the avalanche. For late times, the fastestgrowing component F λ1,λ1 becomes dominant. The corresponding coefficient C(λ 1 , λ 1 , t) can be read off from Equation 64. It approaches a constant C ∞ (λ 1 , λ 1 ) like This shows that avalanche fluctuations at early times (and their interactions with the boundaries) are important to be able understand the behaviour of the avalanche at late times. Intuitively, this fact can already be appreciated from Figure 1b. The constant C ∞ can be obtained from Equation 63 as In this form, it is well-suited for numerical evaluation, where f eh 1 is taken from Equation 58 and N from Equation 39. This integral also exists analytically, but the explicit expression is too lengthy to be of much practical use.
The late-time solution thus reads, As mentioned below Equation 47, the covariances are given entirely in terms of the functions f 2 in this limit, cov[n e (x), n e (y)] ≈ f ee Keeping in mind that f ee 2 (x, y, t) ≈ C ∞ (λ 1 , λ 1 )f e λ1 (x, t)f e λ1 (y, t) at late times, we see that lim t→∞ cov[n e (x), n e (y)] σ(n e (x))σ(n e (y)) = lim Similarly, we find that and lim t→∞ cov[n e (x), n h (y)] σ(n e (x))σ(n h (y)) = lim The above results have important implications. Equations 67-69 correspond to the Pearson correlations between electron-and hole densities at positions x and y: at late times, the charge distributions inside the avalanche region become maximally correlated across space. This also implies that the absolute numbers of electrons and holes become maximally Pearson-correlated at late times, By the same argument, we see that the relative fluctuations of the charge densities defined in Equation 48 become independent of the position x at late times, and identical to the relative fluctuation of the total charge content, lim t→∞ σ(n e (x)) 2 This demonstrates again that, at late times, the only degree of freedom of the charge densities n e (x) and n h (x) lies in their overall normalisations: once maximal spatial correlation is achieved, the fluctuations become identical for all x. The calculation also shows that the relative fluctuations for electrons and holes in fact become identical at late times, This is in direct correspondence to Equation 14 for the case of an avalanche developing in an infinitely large semiconductor. For a constant electric field, the relative fluctuation is conveniently computed as where C ∞ (λ 1 , λ 1 ) is taken from Equation 65 and C(λ 1 ) from Equation 41. The total charge content of the avalanche, and thus its relative fluctuation, is not generally an experimentally accessible quantity. Of direct relevance for the achievable time resolution is the relative fluctuation of the induced current I(t) at late times, i.e. the ratio lim t→∞ σ(I(t)) 2 / I(t) 2 . Taking the form of the signal from Equation 18 and making use of Equation 70, we find that it again limits to lim t→∞ σ(I(t)) 2 Note that the strategy for the solution of Equations 53-55 and in particular the ansatz in Equation 62 also works in the general case of position-dependent Townsend coefficients and velocities. For late times, the position-dependence of f ee 2 , f hh 2 and f eh 2 thus continues to be given by products of the eigenfunctions of Equations 28-29 (which must now be found numerically). Maximal correlation occurs regardless of the shape of the field profile E(x). Also Equations 70 and 72 continue to hold in the general case and allow to compute the relative fluctuations from the numerical solutions for f ee 2 and f hh 2 .

Time resolution
We focus on the discussion of the time resolution for large thresholds placed on the induced current I(t), which is the case of practical relevance. Our discussion here parallels the treatment given in Section 2.6 for infinitely large semiconductors and constant electric fields.
Before quantifying the time resolution itself, it is worthwhile to point out a conceptual difference to the case studied there. An avalanche developing in a bounded semiconductor will generally diverge only with a finite probability (or "efficiency") < 1. Correspondingly, a fraction 1− of all avalanches will end once all charges have escaped across the boundaries. The (generally unknown) probability distribution p(N e , N h , t) for this case has the structure where 1 − (t) is the fraction of all avalanches that have already ended at time t, and the distribution p S (N e , N h , t) describes the "starters", i.e. those avalanches that have not (yet) ended. Once an avalanche has ended, it can never diverge again, i.e. (t) is monotonically falling with t. Its limiting value = lim t→∞ (t) can be easily computed as shown in [2,3]. For our purposes, the efficiency is merely an external input parameter like the Townsend coefficients α and β. It depends on the applied electric field and on the position x 0 of the charge carrier initiating the avalanche. For the computation of the time resolution, therefore, we must restrict all expectation values to only include those avalanches that actually diverge (and therefore cross the threshold). We use a subscript S to indicate that only the starters are considered. In complete analogy to the steps leading to Equation 26, the time resolution for large thresholds is given by However, no closed system of evolution equations exists for the exact evaluation of σ(log I) S , even in the case of a constant electric field (see the discussion in Appendices A.1.8 and A.2.7). Instead, we proceed approximately. We saw in Section 2.2 that in the simple case of a infinitely extended semiconductor and a constant electric field, the avalanche parameter A is directly related to the relative fluctuation via lim t→∞ σ 2 (N )/ N 2 = A −1 . For this case, Equation 27 then establishes a direct link between the relative fluctuation of the quantity on which a threshold is placed, and the corresponding time resolution. The exact form of this relationship depends on the shape of the underlying probability distribution p(N ).
For the general case considered here, we do not have access to the shape of this distribution, and so the time resolution can not be directly connected to the relative fluctuation of the signal, σ(I) 2 / I 2 , as computed in Section 3.2. However, under the approximation that this relationship is in fact similar to that derived in Section 2.2, the time resolution for large thresholds is given by where the effective avalanche parameter A eff is defined as To evaluate this expression further, we note the following consequence of Equation 73. We have that for any observable O that is homogeneous in N e and N h , i.e. for which N e = N h = 0 =⇒ O(N e , N h ) = 0. This is true in particular for N and N 2 , and so, for late times, N = N S and N 2 = N 2 S . The effective avalanche parameter can therefore be computed in terms of unconditional averages as where the last equality is valid only for constant electric fields. For the general case, σ(N ) 2 / N 2 can be taken from the numerical solution as explained in Equation 70.
Inspecting the expression for the time resolution in Equation 75, we see that its overall scale is set by the inverse of the average asymptotic growth rate of the avalanche, γ 1 v * . It can be computed from Equation 32 and Figure 8 and generally accounts for the predominant part of the material-and fielddependence of the time resolution.
Scaling out this trivial dependency, the dimensionless quantity σ th,∞ γ 1 v * measures the effect of avalanche fluctuations on the time resolution. Figure 12 summarises the results from a series of MC simulations that show the evolution of this coefficient with the Townsend parameter β and the thickness d of the semiconductor. Its dependence on the material properties is usually rather weak. The approximation made in Equation 75, σ th,∞ γ 1 v * ≈ ψ 1 (A eff ), is accurate to within about 10% across a wide range of the ratio α/β. For avalanches that are initiated by an electron, it becomes exact in the limit α β. This case is relevant e.g. for detectors based on silicon. As shown in Figure 13, σ th,∞ γ 1 v * depends only weakly on the position x 0 of the charge that triggers the avalanche. Variations are limited to about 30% as x 0 traverses the medium.
Neglecting the finite thickness of the semiconductor altogether, i.e. taking A eff = A generally reproduces the correct scaling of σ th,∞ γ 1 v * with β. Of course, this approximation does not capture the residual dependence on the thickness d and the initial position x 0 .

Conclusions
We have given a comprehensive description of electron-hole avalanches developing in a thin semiconductor that is exposed to a strong electric field, a structure that is commonly used for the amplification and detection of a small amount of initial charge.
Starting from a series of differential equations, we describe the average evolution of the avalanche as well as its fluctuations around this average. For general inhomogeneous electric fields, these equations can be efficiently solved with numerical methods. For constant fields, we give their analytical solutions, including compact expressions for the average signal induced on the readout electrodes of a particle detector.
We also study the achievable time resolution for devices that operate above the breakdown limit, e.g. SPADs and SiPMs. The time resolution depends on the average growth rate of the avalanche as well as a numerical coefficient of order unity which encodes fluctuations of the avalanche. The latter can be approximated with an accuracy of about 10% using our methods and shows only a relatively weak dependence on the strength of the electric field and the thickness of the semiconductor.

A. Avalanches as stochastic many-body systems
Under the assumptions that all charge carriers are independent of each other, and that multiplication processes are memoryless, the charge avalanche considered in Section 1.1 turns into a Markov process. That is, a full description of the avalanche at a certain time t 0 (the current "state" of the avalanche) contains enough information to allow for the computation of the state of the avalanche at any future time t > t 0 .
A charge avalanche evolves stochastically, and so a complete description of its state is equivalent to a specification of the probability p(A) with which a certain avalanche configuration A occurs. By "avalanche configuration" we mean a sufficiently detailed description of a concrete realisation of an avalanche. In practice, this can take various forms. In the most general (and most complex) case, A contains the spatial positions of all participating charges (or, equivalently, their spatial densities). In many situations, however, merely keeping track of the total number of charges is sufficient. We will consider both of these situations in what follows.
The probability p(A) will generally change with time: we write p(A, t). For a Markov process, the equation describing this time evolution is of the form The two terms on the right-hand side correspond to the two ways in which the probability p(A, t) can change in the small time interval [t, t + dt]. First, an avalanche with a configuration A can evolve into A, thereby increasing p(A, t) by an amount proportional to p(A , t). This is captured by the first (positive) term in Equation 77, where T (A → A) labels the rate at which this transition occurs. The second (negative) term represents avalanches that are evolving away from the configuration A into another configuration A , thereby reducing p(A, t) by an amount proportional to itself.
The set of all transition rates T (A → A ) are collectively referred to as the "transition matrix". The transition matrix elements encode the dynamics of the avalanche process. However, writing them down explicitly can be very complicated, in particular if the description of the avalanche configuration A is very detailed. It is therefore important to set up an efficient bookkeeping mechanism to render the calculations manageable.
In the remainder of this Appendix, we introduce such a formalism, based on the occupation-number representation familiar from quantum mechanics [13,14,15,16]. Note that this does in no way imply that we treat the avalanche as a quantum system. Rather the opposite: the way in which many-body quantum mechanics is conventionally formulated provides a very convenient framework for the description of a general Markov process.
Section A.1 applies these ideas to the case of an electron-hole avalanche evolving in an unbounded semiconductor that is exposed to a constant electric field. Section A.2 then treats the significantly more complicated case of a thin semiconductor and a general position-dependent electric field.

A.1. Electron-hole avalanche in an infinite semiconductor with a constant electric field
For a position-independent electric field E, the Townsend coefficients α and β and the velocities v e and v h become constants. Then, the spatial distribution of the charge carriers in the avalanche is irrelevant and the avalanche configuration is fully specified by the total numbers of electrons (N e ) and holes (N h ) in the medium.

A.1.1. Avalanche configurations and state vectors
As in many-body quantum mechanics, we represent such an avalanche configuration by a vector in a Fock space that we denote as |A = |N e , N h . The "empty" avalanche that does not contain any charges is consequently written as |0, 0 = |0 . These vectors serve as a basis, from which more complicated states can be assembled.
To model the evolution of an avalanche, we will routinely need to transition between basis vectors with different electron and hole contents, e.g. |A = |N e , N h → |A = |N e , N h . To this end, we introduce creation and annihilation operators for electrons (â † e ,â e ) and holes (â † h ,â h ) that act on these states in the way familiar from quantum mechanics. The defining properties of these operators are encoded in their commutation relations (for two operatorsÂ andB, their commutator is [Â,B] :=ÂB −BÂ). They read and The constant appearing on the right-hand side of Equation 79 is chosen by convention and its precise value is not important. The fact that it is nonzero is of great physical relevance: there is one more way of first adding a particle to an arbitrary avalanche and then removing one, than doing it the other way around.
Using the creation and annihilation operators, we can also define number operators for electrons and holes,N e =â † eâe ,N h =â † hâ h . Equations 78 and 79 lead to the following commutation relations for the number operators, All commutators that combine an operator involving electrons with an operator involving holes, e.g.
Making use of the creation operators, we identify the basis states as i.e. as the empty avalanche, to which the correct numbers of electrons and holes have been added. That this assertion is in fact consistent with the definitions ofâ † ,â andN is easy to see. Applying the commutation relations in Equation 80 shows that the states |N e , N h indeed contain a well-defined number of charges, as measured by the number operators, The action of the creation operators on the basis states is trivially given bŷ while the relations in Equation 80 describe the action of the annihilation operators aŝ and in particular,â e |0, 0 = 0,â h |0, 0 = 0. The basis states |N e , N h describe specific configurations A of the avalanche. We represent a generic "state" |ψ(t) of the avalanche as the vector corresponding to the sum of all allowed configurations, weighted by their respective probabilities p(A, t) = p(N e , N h , t), This notion of a general state will turn out to be very useful for the computation of expectation values. Equation 81 can be thought of as the analogue of superposition in quantum mechanics. An avalanche which attains a particular configuration |N e , N h with probability one is thus represented by the corresponding basis vector, i.e. |ψ = |N e , N h . h . This allows us to denote the inner product of two vectors |N e , N e and |N e , N h as N e , N h |N e , N h . With respect to this inner product, the adjoints of the operatorsâ e andâ h areâ † e andâ † h -a fact already anticipated by the notation. We define the vector representing the empty avalanche to be normalised as 0|0 = 1. Then, using the commutation relations in Equations 78 and 79, we can evaluate the inner product between two basis vectors to give

A.1.2. Dual vectors and inner product
This shows that the basis states |N e , N h are orthogonal (but not orthonormal) w.r.t. this inner product. Given a general state |ψ(t) like in Equation 81, we can recover the probability p(N e , N h , t) by taking the inner product of |ψ(t) with the basis state N e , N h |,  (â e +â h ) n n! , which has the important property of being a left-eigenstate of the creation operators. In other words, we have avg|â † e = avg| , avg|â † h = avg| . Its overlap with the empty state is avg|0 = 0|1 + . . . |0 = 1. From this, it follows that also avg|N e , N h = 1.
To compute the expectation in Equation 83, we turn the observable O(N e , N h ) into an operator by inserting the number operators into its arguments: O(N e ,N h ). Acting with this operator on the state |ψ(t) and taking the inner product with avg| computes the expectation O(N e , N h ) , This makes explicit that expectation values are linear functions of the state vector, and thus of the probabilities. This is to be contrasted with quantum mechanics, where the expectation value of an operator O is O = ψ|O|ψ and thus depends quadratically on the amplitudes. An important special case of Equation 84 occurs for O = 1, which shows the normalisation of p(N e , N h , t) as actual probabilities, N h , t).

A.1.4. Time evolution
The remaining aspect is to implement the time evolution of the stochastic system, i.e. to describe the time dependence of |ψ(t) , or equivalently, of p(N e , N h , t). To this end, we introduce the time translation operatorÛ (dt) = 1 + dtĤ which acts on the state |ψ(t) asÛ (dt) |ψ(t) = |ψ(t + dt) . The operatorĤ encodes the dynamics of the system, and for this reason we refer to it as the Hamiltonian. From the definition of the time translation operator, we obtain the evolution equation for the state |ψ as For the special case O = 1, we obtain the relation 0 = avg|Ĥ|ψ(t) .
SinceĤ |ψ(t) = 0, we have the condition avg|Ĥ = 0 that any valid Hamiltonian must satisfy. Therefore, Equation 87 can also be written in the form We shall find that this evolution equation is extremely useful. It can be used to extract important information about the system even in complicated situations where Equation 85 itself can no longer be solved directly.

A.1.5. Hamiltonian for electron-hole avalanche
In order to derive the form of the Hamiltonian for an electron-hole avalanche, it is useful to expand on the connection between Equation 85 and the general form of the time evolution for a Markov process in Equation 77.
To this end, we write the time evolution operatorÛ (dt) in terms of the basis states |N e , N h , The operator matrix elements p(N e , N h → N e , N h ; dt) in this basis correspond to the probabilities that a transition |N e , N h → |N e , N h happens within a small time interval dt. This is seen by acting with Û (dt) on a general state vector |ψ(t) , where have shifted the notation N e → N e , N h → N h in the last line. This shows that the operatorÛ (dt) as expanded in Equation 89 indeed effects an infinitesimal time translation t → t+dt when applied to any state. The probabilities p(N e , N h → N e , N h ; dt) are related to the transition matrix elements as follows, We make use of this correspondence to write down the Hamiltonian for the electron-hole avalanche. It is particularly convenient to work with the expansion in Equation 89. Within a small time interval dt, the avalanche configuration |N e , N h can evolve into the configuration |N e + 1, N h + 1 (with a probability of order O(dt)) or remain itself (with a probability of order unity). Transitions to other states are possible within dt, but with probabilities of order at least O(dt 2 ). A transition occurs if any one of the N e electrons creates an electron-hole pair, or if any one of the N h holes creates an electron-hole pair. Since all of these processes occur independently, the corresponding probability is Similarly, the probability for the avalanche not to leave the state |N e , N h within dt is We can thus express the time evolution operatorÛ (dt) aŝ Writing We can now read off the HamiltonianĤ, which iŝ We have derived the Hamiltonian by starting from the transition probabilities p(N e , N h → N e , N h ; dt). One could have also started from a "microscopic" point of view, i.e. by considering the interactions that each charge can undergo. We shall not expand on this picture here, but merely hint at the related diagrammatic representation of Equation 91. We represent the creation and annihilation operators as lines leaving or entering any given diagram, as shown in Figure 14a. The representation of the four terms in the Hamiltonian, expressed in terms ofâ † andâ are shown in Figure 14b: they correspond exactly to all possible interactions available to a charge carrier. The probability per unit time for an interaction to happen appears as a numerical coefficient in the Hamiltonian and plays the role of a coupling constant.
Putting these terms together immediately leads to Equation Ne N e ≈ 1. This suggests that we might find an observable which is not generally identical to logN , but which limits to this quantity for large avalanches, and which does obey a closed evolution equation.
To this end, we consider the observableŝ with the polygamma functions ψ k . They indeed limit to log N for large avalanches, lim Ne,N h →∞f where C 0 is an irrelevant constant that will drop out. This means that The commutators off 1 andf 2 with the Hamiltonian are trivial, and so are their time evolution equations, Together with Equation 94, this shows that The asymptotic time resolution thus becomes as stated in the main body.
A.2. Electron-hole avalanche in a thin semiconductor with an arbitrary electric field We now proceed to including the spatial development of the avalanche. This is important in case the avalanche develops in a bounded region where the electric field E(x) depends on the position x. This leads to position-dependent Townsend coefficients α(x) and β(x) and velocities v e (x) and v h (x).

A.2.1. State vectors
In this case, an arbitrary configuration A of the avalanche is labelled by the spatial densities of electrons and holes, n e (x) and n h (x). The number of electrons in a small interval [x, x + dx] is then n e (x)dx and the number of holes in the same region is n h (x)dx. The basis vectors spanning the space of all avalanche configurations are then |A = |n e (x), n h (x) .
Again, we introduce creation and annihilation operators to transition between different configurations. These operators now add or remove a charge at position x, i.e. they pick up a position index. As before, their action is defined by the commutation relations Using the creation operators, we can again build up the basis vectors starting from the empty configuration |0 , We also have the number density operatorŝ Using relations 95 and 96, we see that the number operators obey the following commutation relations Their action on the basis vectors is thus given bŷ n e (y) |n e (x), n h (x) = n e (y) |n e (x), n h (x) ,n h (y) |n e (x), n h (x) = n h (y) |n e (x), n h (x) , as expected. A general state |ψ(t) is again built as a linear combination of basis vectors, weighted by the probability to find the corresponding avalanche configuration, p(A, t). Since A is labelled by the functions n e (x) and n h (x), this probability is now formally a functional, p(A, t) = p[n e (x), n h (x), t], where we use square brackets to denote the functional dependency. A generic state |ψ(t) is then given by the functional integral |ψ(t) = Dn e (x)Dn h (x) p[n e (x), n h (x), t] |n e (x), n h (x) .

A.2.2. Translations
The movement of charge carriers through the semiconductor is central to the formation of the avalanche in position space. We introduce dedicated operatorsT e andT h that implement the effects of spatial translations on the basis vectors |n e (x), n h (x) separately for electrons and holes.
Under a general translation T ∆ , the position of a charge carrier changes as The translation is characterised by the function ∆(x) that describes the magnitude of the (generally position-dependent) shift. The translation operators thus depend on ∆(x) in a functional way. We writê Under a translation of the type of Equation 97, a generic density n(x) transforms as T ∆ n(y)dy = n(x)dx =⇒ T ∆ n(y) = n(x) where we have used the shorthand notation T ∆ n(x) to represent the transformed density. For constant shifts, ∆(x) = ∆, the Jacobian factor 1 + d∆ dx −1 vanishes and the density gets translated as y = x + ∆.
The transformation in Equation 98 is homogeneous, which means that the empty state is invariant under translations,T e |0 = |0 andT h |0 = |0 , as expected. Equivalently, we can describe the translation operators by giving their commutation relations with the number density operatorsn e (x) andn h (x). Using the definitions in Equations 99-100, we find In these expressions, the action of a derivative on a generic operatorÔ(x) is to be understood as which is inspired by the fact that, for smooth expectations,

A.2.3. Global observables
Starting from the number density operatorsn e (x) andn h (x), we can also construct operators that are sensitive to the total number of charges in a finite region R, Their commutation relations with the creation and annihilation operators are given by [N e (R),â e (x)] = −â e (x) · 1 R (x), where 1 R is the indicator function of R, i.e. 1 R (x) = 1 if x ∈ R and 0 otherwise.
In particular, the important property avg|n e (x), n h (x) = 1 holds as before.
To compute the expectation value O[n e (x), n h (x)] , we again replace the arguments with the number density operators, apply the resulting operator to the state |ψ(t) and take the inner product with avg|.

A.2.5. Time evolution operator for electron-hole avalanche
We again introduce the time evolution operatorÛ (dt) = 1 + dtĤ that produces an infinitesimal time evolution of the state vector |ψ(t) . The time evolution of the state |ψ(t) is still given by Equation 85. It now takes the form of a functional differential equation. The time evolution equation for the expectation value O[n e (x), n h (x)] in Equation 88 generally becomes a partial differential equation.
In the present case,Û (dt) is easier to write down than the HamiltonianĤ itself. Since [Ô, 1] = 0, we see that dt Ô ,Ĥ = Ô ,Û (dt) for an arbitrary observableÔ, and so both operators essentially contain the same amount of information. We start from Equation 90 for the position-independent situation and generalise this expression to the general case considered here. Since charges at different coordinates multiply independently, it is easy to account for the position-dependence itself: we need to replaceN e →n e (x),N h →n h (x), α → α(x), β → β(x), v e → v e (x), v h → v h (x) and integrate over x. In the language of Figure 14, this introduces a separate interaction vertex for each position x. The translation operators introduced in Section A.2.2 can be used to implement the movement of electrons and hole in opposite directions, and so the final expression forÛ (dt) readŝ

A.2.6. Evolution equations for moments
To derive evolution equations for the first moments n e (x) and n h (x) , we require the commutation relations (dropping contributions of order O(dt 2 ) or higher) Extracting the term of order dt and computing the average avg| · |ψ immediately leads to Equations 28-29. For the second moments n e (x)n e (y) , n e (x)n h (y) and n h (x)n h (y) , the necessary commutators read (again dropping terms of order O(dt 2 ) or higher) For efficiencies ≈ 1, the left-hand side becomes proportional to σ(log N ) S at late times. Compared to Equation 94, the right-hand side does not vanish, but implements the effects induced by the boundary. It essentially depends on the degree of correlation between the charges currently present in the avalanche region, and the relative rate at which charges disappear across the boundaries. Note that this equation is not closed, i.e. the right-hand side contains variables whose evolution is not determined by the equation itself.

B. Orthogonality of eigenfunctions
Making use of the definition of the inner product in Equation 37, we find for λ = λ I(λ, λ ) = f λ e (·, t) f λ h (·, t) C. Analytic expressions for f ee 1 and f hh