Passive quenching, signal shapes, and space charge effects in SPADs and SiPMs

In this report we study the dynamics of passive quenching in a single-photon avalanche diode. Our discussion is based on a microscopic description of the electron-hole avalanche coupled to the equivalent circuit of the device, consisting of the quench resistor and the junction capacitance. Analytic expressions for the resulting signal shape are derived from this model for simple electric field configurations, and efficient numerical prescriptions are given for realistic device geometries. Space charge effects arising from the avalanche are included using simulations. They are shown to distort the signal shape, but alter neither its basic characteristics nor the underlying quenching mechanism.


Introduction
Single-photon avalanche diodes (SPADs) and silicon photomultipliers (SiPMs) are Geiger-mode avalanche detectors with direct single-particle sensitivity to photons and charged particles [1,2,3]. In these devices, the small amount of charge deposited by the primary particle triggers a self-sustaining electron-hole avalanche which leads to the creation of a detectable electric signal. The avalanche develops in a thin region with very high electric fields, created in a p-n diode which is initially reverse-biased above its breakdown voltage.
The resulting exponential growth of the avalanche must be stopped ("quenched") at a suitable moment so that all charge carriers are removed from the junction and the detector can register the next event.
In passively quenched devices, the diode is decoupled from the bias supply through a series resistance. The charge produced by the avalanche then reduces the electric field in the junction and suppresses the avalanche growth. Quenching occurs when the voltage across the diode drops below the breakdown limit. Once all charges have drifted out of the junction, the original bias voltage is restored.
The resulting voltage and current waveforms form the output signal. The signal properties are ultimately determined by the complicated interplay of impact ionisation, the quenching circuit, and the diode capacitance. The operating conditions are typically chosen such that the achieved signal charge is around 10 4 -10 5 electrons, i.e. very large avalanches must be considered.
Existing analyses of the signal formation in passively quenched SPADs either model the avalanche discharge and quenching purely in terms of an equivalent electrical circuit [4,5,6], or use a coarse effective description of the avalanche development, e.g. in terms of rate equations [7,8]. In this report, we construct a simple analytic quenching model which is directly based on the microscopic dynamics of electron-hole avalanches. It naturally connects to our previous treatment of avalanche statistics given in Ref. [9]. As a result, our model describes the detector signal in terms of the underlying material parameters and the device geometry. It can therefore guide the evaluation of specific detector designs and allows direct comparisons to be made with measurements. This paper is structured as follows. Section 2 briefly reviews the avalanche formation and summarises properties of large electron-hole avalanches used to formulate the quenching model. Analytic predictions for the signal shape and the spatial distributions of charge carriers in the high-field region are derived in Section 3. The impact of the space charge field of the avalanche on the quenching process is then studied in Section 4.

Properties of large electron-hole avalanches
The development of the electron-hole avalanche is determined by the transport of charge carriers of both polarities through the multiplication region and their interactions with the material. For sufficiently strong electric fields, both electrons and holes can cause impact ionisation and produce additional electronhole pairs. The probability for a charge carrier to undergo a multiplication reaction may be parameterised in terms of the material parameters and, in general, the past history of the carrier. The drifting carriers in the avalanche induce a current I ind on the metallised device contacts, which may be computed with the Ramo-Shockley theorem [10,11]. Impact ionisation is a stochastic process. The avalanche development is thus subject to fluctuations, which are particularly pronounced at early times when the total number of participating charge carriers is still small. For avalanches containing many charge carriers, these fluctuations average out and the further evolution appears deterministic. (The continued evolution of an avalanche containing N charges will generate relative fluctuations with a magnitude proportional to 1/ √ N , as discussed in Section 3.5 of Ref. [9].) It is useful to introduce a current scale I det to separate these two regimes, chosen such that avalanche fluctuations are important whenever I ind I det and a deterministic description of the avalanche is appropriate when I ind I det . The current I det therefore defines the level of fluctuations that are deemed "important" in a particular situation. It implements a purely semantic definition, and expressions for physical observables do not depend explicitly on I det .
Our model of the quenching process is based on the following two general properties of large electronhole avalanches (I ind I det ) developing in a time-independent electric field, illustrated in Fig. 1: (i) The shape of the induced current I ind is identical for each avalanche event. It is directly proportional to the average current I ind , where the average is taken over all avalanche evolutions starting from the same initial conditions, The proportionality factor k is a random variable. The resulting fluctuations around the average evolution determine the time resolution, but are irrelevant for the study of the quenching process. (ii) The average current I ind evolves exponentially as and is parameterised by the dimensionful quantities I 0 and S 1 . The current I 0 depends on the number and the spatial distribution of the primary charges initiating the avalanche. The parameter S 1 represents the asymptotic net growth rate of the avalanche, i.e. it takes into account the generation of charges in the high-field region and their outflow across its boundaries. It depends on the device geometry and the reverse bias voltage V d applied across the junction, i.e. S 1 = S 1 (V d ).
In the following, the breakdown voltage V br is defined such that S 1 (V br ) = 0. This implies that S 1 (V d ) > 0 for V d > V br , i.e. exponentially growing avalanches are possible if the device is operated "above breakdown". For a bias voltage below the breakdown limit, V d < V br , impact ionisation cannot sustain the discharge. The avalanche current then decays exponentially, i.e. S 1 (V d ) < 0.
For a given device and electric field configuration, the parameters S 1 and I 0 are always accessible through Monte Carlo (MC) simulations of the avalanche development (such as the one in Fig. 1); efficient numerical prescriptions for quasi-one-dimensional situations are referenced below. The above properties arise as explicit predictions of the avalanche model studied analytically in a one-dimensional geometry in Ref. [9]. It is briefly reviewed in the following. We use this model to illustrate our results in Sections 3 and 4. However, our description of the quenching process is more general and applies to all scenarios in which the avalanche obeys the properties (i)-(ii) above, including realistic three-dimensional device geometries.
One-dimensional memoryless avalanche model. The diameter D of a typical SPAD pixel is significantly larger than the longitudinal thickness d of the high-field region, typically D/d 10. Close to the centre of a pixel, this results in a quasi-one-dimensional field geometry through which electrons and holes drift in opposite directions with velocities v e and v h , respectively. The drift direction is taken to be along the x-axis (Fig. 2a).
The one-dimensional model of Ref. [9] describes the avalanche development in this geometry. In this model, impact ionisation is parameterised by the impact ionisation coefficients α and β, which represent the probability per unit length for an electron or hole to undergo a multiplication reaction (Fig. 2b). This probability is assumed to depend only on the local electric field strength E(x), i.e. α(x) = α(E(x)) and β(x) = β(E(x)). Any dependence on the history of the charge carriers is neglected, i.e. the avalanche is assumed to be a memoryless stochastic process.
For general field profiles E(x), the growth parameter S 1 may be obtained from the numerical solution of the differential equations describing the evolution of the average charge content of the avalanche, Eqs. 30 and 31 in Ref. [9]. These equations also determine the spatial carrier densities in the junction. As shown in Appendix B of Ref. [3], the breakdown condition S 1 = 0 may be expressed in terms of the impact ionisation coefficients as It is thus equivalent to the breakdown condition of Ref. [12], which evaluates the point at which the gain for a constant current injected into the junction diverges. For constant fields, S 1 may be expressed as the product of the effective carrier drift velocity v * = 2v e v h /(v e + v h ) and the growth factor γ 1 , i.e. S 1 = γ 1 v * . The latter is defined in Eq. 35 in Ref. [9]. The drift velocities for electrons and holes, v e and v h , are themselves dependent on the local electric field. The spatial carrier densities for electrons and holes, n e (x, t) and n h (x, t), are proportional to the functions f e λ1 (x, t) and f h λ1 (x, t), respectively. Their analytic expressions are given in Eqs. 36 and 37 of Ref. [9]. The numerical parameter λ 1 is related to the time constant S 1 and defined in Eq. 38 of Ref. [9].
The proportionality factor k in Eq. 1 is approximately distributed according to a gamma distribution with shape parameter A and scale parameter A −1 (cf. Section 3.6 of Ref. [9]), where Γ(z) is the gamma function, and the avalanche parameter A is defined as for position-independent electric fields. In Eq. 5, N 0 e and N 0 h label the numbers of primary electrons and holes, respectively. For general field profiles E(x), the avalanche parameter may be computed as shown in Eq. 79 of Ref. [9].

Passive quenching in the absence of space charge effects
To describe the quenching process, the avalanche evolution summarised in Section 2 is coupled to the electrical equivalent circuit of a single SPAD (Fig. 2c). In this circuit, the p-n junction is represented by its capacitance C d , and the current I ind acts as a parallel current source. (Note that the dynamic resistance R d of the diode is automatically included through the dependence of the carrier drift velocities v e and v h on the electric field and is thus not present as an explicit circuit component.) In the absence of free charge carriers in the junction, I ind = 0 and V d = V supply . The excess voltage V ex is defined as ) which is relevant for the development of the avalanche. To ensure quenching, the decoupling resistor R q must be chosen such that the recharging time constant R q C d is much longer than the time scale 1/S 1 (V supply ) on which the avalanche develops. The initial phase of the quenching process thus remains unchanged as R q → ∞.
In this dynamical situation, the instantaneous growth rate of the avalanche, S 1 (t), becomes itself time-dependent. Provided that V d (t) changes on time scales that are long compared to the internal time constants of the avalanche, the changing growth rate may be approximated as is the function introduced in Section 2 for static situations. This "adiabatic approximation" is shown below to be applicable to practically relevant situations. For R q → ∞ and I ind I det , the combined deterministic system of avalanche and equivalent circuit is then described by the equations More complicated equivalent circuits that include additional parasitic elements or finite values of the quench resistance R q are readily accommodated using the same steps. Eq. 6, which describes the evolution of the avalanche, must then be solved together with the corresponding circuit equations. For the situation discussed here, Eqs. 6-7 may be solved numerically for general S 1 (V d ). This is similar to the rate model constructed in Ref. [7], but in contrast to the calculations performed there automatically includes boundary effects. In many practically relevant situations, V ex /V br 1, i.e. the voltage V d remains close to V br during the full evolution of the system. In this case, the function S 1 (V d ) in Eq. 6 may be replaced by its linear expansion around V br . With the breakdown condition S 1 (V br ) = 0, this gives The parameter K br may be easily obtained through finite differences provided that br . An analytic expression for K br for a memoryless avalanche and position-independent fields is given in Eq. 26 below. If the primary charge deposit occurs at t = 0, Eqs. 6-7 must be solved for t ≥ t det , where t det is the time at which the avalanche current first exceeds I det . The time t det is a random variable whose statistics are further discussed below. Defining V d (t det ) = V det , the relevant initial condition is I ind (t det ) = I det , leading to the following analytic solutions, with the quenching time constant τ q , and the time ∆t, The quantity τ q parameterises the time scale on which V d changes as the avalanche is quenched. For Eqs. 9-10 to be a valid description of the quenching process, i.e. for the adiabatic approximation to be applicable, τ q must be larger than the transit time d/v * . The latter approximates the time scale on which the spatial charge distribution in the junction, and therefore also the instantaneous growth rate S 1 (t), react to changes of V d .
The general expressions in Eqs. 9-10 can be further simplified. For practical devices, the diode capacitance C d typically amounts to a few femtofarads. With the elementary charge e 0 , around C d V ex /e 0 ∼ 10 4 charge carriers need to be produced by the avalanche in the high-field region before any significant discharge of C d occurs, i.e. any avalanche fluctuations remaining during the quenching process amount to a few percent at most. Working to this accuracy, In this limit, Eq. 11 simplifies to V eff ≈ V ex and the quenching time constant is directly given in terms of the excess voltage, and the time ∆t defined in Eq. 13 becomes With this, ∆t/τ q 1 and Eqs. 9-10 simplify to and the adiabatic approximation is valid provided that K br V ex 2v * /d. The current signal I ind attains a peak value of I ind max = C d V ex /τ q at time t max = t det + ∆t. If t = 0 is instead regarded as the position of the peak, the signal shape reads Its evolution is symmetric around this peak with a FWHM of The observable voltage step and The signal charge is thus directly proportional to the excess voltage at which the junction is operated and does not depend on the primary charge deposit which initiated the avalanche. The physical origin of the result in Eq. 21 is easily understood. As long as V d (t) > V br , impact ionisation is self-sustaining and the avalanche grows exponentially. The avalanche, and the induced current I ind , consequently attains its maximum size when V d (t max ) = V br = V supply − V ex . Impact ionisation continues to occur for V d (t) < V br , but is no longer self-sustaining. The current I ind then decays exponentially, reducing V d below V br and increasing the output voltage step ∆V d beyond V ex . In situations where the approximation in Eq. 8 is  appropriate, the collapse of the avalanche exactly mirrors its initial growth; the evolution of V d (t) is then symmetric around V br and so ∆V d = G V ex with G = 2. Note that this conclusion deviates from other commonly-used quenching models, which take the avalanche discharge to stop already at V d (t) ≈ V br and thus lead to G ≈ 1 [4,5,6]. These differences deserve a detailed comparison and discussion, which we carry out at the end of this Section. The timing of the detector signal also carries important information. The average time of maximum current, is experimentally accessible in a setup where the primary charge deposit occurs at a well-defined time, e.g. where charges are deposited by means of a laser pulse or a charged-particle beam. It contains information about the early evolution of the avalanche (through t det ), as well as the quenching process itself (through ∆t). Eqs. 1 and 2 result in the following expression for t det , which, together with Eq. 15, gives for t max The logarithmic expectation log k depends on the magnitude of the avalanche fluctuations; an analytic expression valid for the memoryless avalanche model is derived in Eq. 27 below.
Quenching for one-dimensional memoryless avalanche. We illustrate these results for the case of a simplified SPAD, for which we approximate the field distribution in the multiplication region with a positionindependent, one-dimensional electric field. The field is thus taken to be E(x, V d ) = −x V d /d, wherex is the unit vector along x (cf. Fig. 2a). In this case, the breakdown voltage is related to the breakdown field E br as V br = E br d. The effective drift velocity v * typically saturates in the high-field region, i.e. becomes independent of V d . Using the memoryless avalanche model from Section 2, where α br = α(E br ), β br = β(E br ), α br = dα/dE| E br , β br = dβ/dE| E br , andκ br = d|α br − β br |/2. In the adiabatic approximation, predictions for the charge carrier densities of electrons and holes are proportional to the functions f e λ1(V d ) and f h λ1(V d ) . The linear expansion in Eq. 8 describes γ 1 (V d ) to excellent precision for typical values of V ex . This is visualised in Fig. 3 for silicon, using the parameterisations of the impact ionisation coefficients from Ref. [13].
With the distribution of k in Eq. 4, log k = ψ(A) − log(A) with the digamma function ψ(z), and the expression for t max in Eq. 25 becomes and I 0 may be computed for a given initial charge deposit as shown in the paragraph following Eq. 46 in Ref. [9]. The analytic results derived above are compared to a MC implementation of the memoryless avalanche model, coupled to the quenching circuit of Fig. 2c. The avalanche simulation is performed as described in the Appendix, using the impact ionisation coefficients for silicon from Ref. [13] and saturated drift velocities v e = v h = v * = 0.1 µm/ps [14]. The circuit equations are integrated using the forward Euler method. The induced current I ind is computed assuming the weighting field of a parallel-plate geometry with an electrode spacing of d for the multiplication region. A dedicated conversion layer, common in SPADs optimised for infrared photon detection, is not included. The full simulation is referred to as "memoryless avalanche MC + circuit" (MAMC + circuit). For the comparisons shown here, the thickness d of the high-field region is chosen as d = 0.5 µm and its diameter is D = 10 µm. Assuming a parallel field geometry and using r = 11.7 [15], the junction capacitance is C d ≈ 16 fF. A quench resistor with R q = 200 kΩ is used. The breakdown voltage for this geometry is V br ≈ 20.34 V and with Eq. 26 the model parameter K br evaluates to K br ≈ 0.105 (V ps) −1 . The quenching time constant is τ q ≈ 9.5 ps for V ex = 2 V and τ q ≈ 6.3 ps for V ex = 3 V. The transit time d/v * is 5 ps, i.e. the adiabatic approximation is valid for excess voltages not significantly above 3 V. Fig. 4 shows the induced current I ind for several avalanche events and different excess voltages, obtained with the MAMC + circuit simulation model. Fluctuations in the early development of the avalanche lead to a stochastic displacement of the deterministic current waveform from Eq. 16. The simulated time of maximal current t max is 55.9 ps (40.8 ps) for an excess voltage V ex of 2 V (3 V). These values are well reproduced by Eq. 27, which evaluates to t max ≈ 56.1 ps (40.9 ps) for V ex = 2 V (3 V). Fig. 5 compares the shapes of the approximate solutions in Eqs. 16-17 to the result from the MAMC + circuit simulation model around the time t max , demonstrating very good agreement between the analytic result and the simulation. Small deviations arise as V d crosses the breakdown voltage, resulting in a slight underestimation of ∆V for large excess voltages. This is caused by additional dynamic effects originating from changes in the shapes of the spatial carrier densities of electrons n e (x, t) and holes n h (x, t) that are not included in the adiabatic model. The temporal evolution of these densities is compared in Fig. 6. (For the simulation, the average densities n e and n h are shown to ensure a meaningful comparison with the adiabatic model even at early times where the number of charge carriers is small.) Charge carriers of both polarities are present throughout the entire multiplication region. The highest densities are attained at the corresponding downstream ends, i.e. at x = d for electrons and at x = 0 for holes. The shape of the distributions, and in particular the curvature of n h (x), depends on V d .     Discussion and comparison with other quenching models. The relation ∆V d = G V ex with G = 2 is derived in Eq. 21 as an exact analytic result that holds in situations where the linear expansion in Eq. 8 is appropriate. This conclusion differs from predictions made with simpler, widely-used, quenching models [4,5,6]. In these models, the avalanche is taken to directly transition the device from its initial inactive configuration (V d = V supply and I ind = 0) into a quasi-time-independent, active state where V d = V br and a persistent current I p = V ex /R q flows through the quench resistor. If R q is chosen sufficiently large, the current I p is small enough so that avalanche fluctuations may spontaneously quench the discharge after a brief period of time (the probability per unit time for this to happen is exponentially suppressed in I p [6]). This leads to a voltage step of ∆V d = V ex , i.e. G = 1. The central proposition is that this behaviour can be modelled by the equivalent circuit shown in Fig. 7, where the diode resistance R d is explicitly included and the avalanche is represented by a switch S whose state marks the presence (closed) or absence (open) of an avalanche discharge.
The microscopic analysis conducted here indeed confirms the existence of this active equilibrium state: Eq. 6 implies dI ind /dt = 0 at V d = V br , which reveals I ind = I(R q ) = I p as a valid configuration of the circuit in Fig. 2c. This shows that the simple model in Fig. 7 can adequately describe the possible active (switch closed) and inactive (switch open) equilibrium configurations of the device.
There is, however, no a-priori reason to expect this passive circuit to correctly model the fast transition between these two states. In particular, the avalanche current I ind in Fig. 7 decays exponentially as V d approaches V br , and V d ≥ V br throughout the entire quenching process. This is apparently inconsistent with the definition of V br as the voltage above which each existing carrier produces a diverging number of secondary charges [6]. To the best of our knowledge, the existing literature does not address (or resolve) this contradiction for switch-based models.
The microscopic treatment reveals that it is not the absolute value of the avalanche current (or, equivalently, the total number of charges in the avalanche) that vanishes close to the breakdown point, but rather its growth rate S 1 , thus giving rise to a different transient behaviour and a proportionality factor G = 2. A similar factor G ≈ 2 has also been obtained in Refs. [7,8] through the numerical solution of a system of rate equations describing the avalanche, and in Ref. [16] through a microscopic MC simulation similar to the one presented here. Refs. [7,8] furthermore report measurements of ∆V d and find G > 1.

Space charge effects
For sufficiently large avalanches, the space charge field E sc (x, t) caused by the charge carriers in the avalanche may become comparable to the field E ext (x, V d (t)) produced as a result of the externally applied voltage V d . The electric field profile E that is relevant for the development of the avalanche is now the sum of both components, i.e. E = E ext + E sc . A full description of the effects of space charge (SC) on the quenching process requires a three-dimensional simulation of the avalanche together with the computation of a self-consistent solution for the space charge field at every time step. In particular, E sc depends on the spread of the avalanche in the direction transverse to the externally impressed field E ext .
The situation is studied here in a simplified way as illustrated in Fig. 8. The evolution of the avalanche continues to be simulated by the MAMC + circuit simulation model, with impact ionisation resulting from the drift of charge carriers along the x-axis. A parameterised model is used to describe the development of the avalanche in the transverse yz-plane, defining a radially symmetric charge profile ρ(r, x, t). Different physical effects may be included in this parameterisation, e.g. charge carrier diffusion as well as drift due to nonvanishing radial electric field components. The prescribed charge distribution ρ is then used to compute the space charge field E sc along the x-axis, which, by symmetry, only has a nonvanishing x- x. This (one-dimensional) field distribution is used to determine the impact ionisation coefficients α(x, t) and β(x, t). The finite transverse extent of the avalanche is thus taken into account when computing its space charge field, but it is neglected when evaluating the backaction of this field on the further evolution of the avalanche. This is analogous to the "1.5dimensional" model of space charge effects in resistive plate chambers (RPCs) discussed in Ref. [17]. As before, the field E ext is assumed to be position-independent, i.e. E ext = −x V d /d. The resulting simulation model is referred to as "MAMC + circuit + SC".
For the results presented here, only transverse diffusion is considered for the construction of ρ(r, x, t). This simplified model is likely to overestimate the true effect of space charge on the development of the avalanche. First, the space charge field E sc itself is overestimated: the x-component of the field Figure 8: For the computation of the space charge field, the high-field region is taken to be delimited by two conducting plates of infinite transverse extent, located at x = 0 and x = d. Charge carriers are constrained to move along the x-axis. This charge distribution is extended into a three-dimensional charge density ρ(r, x, t), illustrated here in the plane x = x 0 . The field Esc(r, x, t) is evaluated on the x-axis as the electric field generated by ρ(r, x, t). The external field Eext(t) also contributes to the total electric field E(x, t) relevant for the development of the avalanche.
originating from the chosen radially symmetric charge density ρ(r, x, t) attains a maximum at r = 0 (cf. Eqs. 28 and 29 below). All drifting charge carriers are thus subject to the maximal field, while regions of lower field strengths at r > 0 are not probed by the avalanche. Second, the transverse size of the avalanche is underestimated by only considering diffusion, but ignoring carrier drift. This leads to a narrower distribution of charges and also enhances the computed space charge field.
Transverse spread of the avalanche and space charge field. Charge carrier diffusion is modelled by prescribing a normal distribution for the transverse distribution of charge carriers, while the longitudinal distribution continues to be given by the densities n e (x, t) and n h (x, t). The densities ρ e and ρ h that are relevant for the computation of the space charge field are thus The transverse extent is governed by the transverse diffusion coefficients D T for electrons and holes, i.e. σ e (t) = 2D T,e t and σ h (t) = 2D T,h t. Values of D T,e = 15 cm 2 /s [18] and D T,h = 11 cm 2 /s [19] are used in the following. To the best of our knowledge, no measurements of the transverse diffusion coefficients are available in the literature for the very strong longitudinal fields |E ext | > 20 V/µm that occur above breakdown. The selected values correspond to a field of 5 V/µm in silicon. The transverse diffusion coefficients depend only relatively mildly on the longitudinal electric field, and it is expected that these values give at least an approximate description of diffusion in the relevant regime. This is sufficient for the qualitative assessment of space charge effects attempted here. The electrostatic Green's function for the parallel-plate geometry of Fig. 8 derived in Refs. [20,21] is used to compute the space charge field E sc .
Quenching dynamics with space charge. Fig. 9 shows the evolution of the carrier densities, the electric field, and the induced signal. The same simulation parameters and material properties as in Section 3 0.
0.25 0.5 0  are used. At early times, space charge effects are of subleading importance, and the resulting charge carrier densities are described by the adiabatic model in good approximation. The separation of charge caused by E ext leads to a space charge field which reduces the electric field in the centre of the gain layer and enhances it close to its boundaries. (Space charge cannot alter the voltage V d measured across the junction contacts, i.e. d 0 ds · E sc = 0.) This results in increased impact ionisation rates and carrier densities near the edges of the multiplication region. For the chosen material parameters and device geometry, this effect overcompensates the suppressed avalanche growth in regions with reduced field. The signal charge and the voltage step ∆V d increase by about 50% compared to the situation without space charge effects, but the signal shape is not significantly altered. This is very different from the space charge effects in RPCs studied in detail in Ref. [17]. Also in the RPC geometry the field integral throughout the sensitive volume remains unchanged so there will be regions with reduced field and regions with increased field. Since in a gas detector there is however only electron multiplication and no ion multiplication, the electrons end up only in these low field regions and the gain is suppressed by many orders of magnitude.
Measurements of the voltage step ∆V d in SPADs can thus serve as a direct probe of the importance of space charge effects for the avalanche development.

Conclusions
In a passively quenched SPAD, the diverging avalanche discharges the detector capacitance, which in turn lowers the electric field in the multiplication region and reduces the probability for impact ionisation. The dynamics of this process is here described quantitatively based on a microscopic model of the avalanche development, applicable as long as space charge effects are negligible. The main conclusions of this treatment are: • The duration of the avalanche discharge scales with the quenching time constant τ q (cf. Eq. 14).
It is determined in terms of the excess voltage V ex and the parameter K br (cf. Eq. 8). The latter may be extracted starting from the field profile E(x) around the breakdown point. It is known analytically for position-independent fields (cf. Eq. 26).
• The avalanche attains its maximum size as the voltage V d measured across the diode crosses the breakdown voltage. The total voltage step ∆V d caused by the discharge is therefore ∆V d = 2V ex . The total signal charge is Q = 2V ex C d , with the junction capacitance C d (cf. Eqs. 21 and 22).
• Analytic expressions for the evolution of V d (t) and the current I ind (t) induced by the avalanche are available in Eqs. 16 and 17. Analytic formulae for the spatial charge carrier densities also exist.
The space charge field produced by the charge carriers in the avalanche enhances the electric field (and therefore impact ionisation) close to the boundaries of the multiplication region and reduces it in its interior. Charge carriers are present throughout the entire junction, and the avalanche growth may, in principle, get enhanced or suppressed by the altered field configuration. A simplified MC simulation of a representative silicon device geometry shows a modest enhancement of the signal charge as a result of space charge. This is in stark contrast to gas detectors such as RPCs, where space charge considerably reduces the gas gain.
Our results in Eqs. 21 and 22 deviate from the prevalent description of passive quenching, in which the avalanche discharge is modelled as a bistable switch. We encourage further experimental work to investigate and clarify this discrepancy.

Appendix
To simulate the development of the avalanche for equal carrier drift velocities, v e = v h = v * , the avalanche region x ∈ [0, d] is partitioned into N equidistant bins of width ∆x = d/N . The units are chosen such that charges drift a distance of ∆x during the simulation time step ∆t = ∆x/v * . The bin centre of bin j is located at coordinate x j and the simulation time steps are indexed as t i = i∆t.
The electric field distribution E(x, t) = −E x (x, t)x and the resulting impact ionisation coefficients α(x, t) and β(x, t) are discretised into E j (t i ) = E x (x j , t i ), α j (t i ) = α(x j , t i ), and β j (t i ) = β(x j , t i ). The charge carrier densities n e (x, t) and n h (x, t) are represented by the number of electrons N j e (t i ) and holes N j h (t i ) in each bin. These carriers are assumed to be uniformly distributed within the bin.
In the absence of impact ionisation, the N j e holes present in bin j at time t i drift into bin j + 1 at t i+1 , and the N j h holes move into bin j − 1. Impact ionisation is simulated as illustrated in Fig. 10. During the time interval ∆t, the N j e electrons (N j h holes) in bin j produce N j e→eh (N j h→eh ) additional electron-hole pairs. These numbers are drawn from a Poisson distribution with rate parameters α j ∆x and β j ∆x, respectively, N j e→eh (t i ) ∼ Po N j e (t i ); α j (t i )∆x , N j h→eh (t i ) ∼ Po N j h (t i ); β j (t i )∆x . Impact ionisation occurs uniformly along ∆x. The N j e→eh electrons created in bin j at t i drift into bin j + 1 at t i+1 , while the N j e→eh holes spread across multiple bins: N j→j−1 e→eh holes reach bin j − 1, N j→j e→eh remain in bin j, and N j→j+1 e→eh enter bin j + 1. The same consideration holds for the N j h→eh electrons created from hole-initiated impact ionisation events, leading to N j→j−1 h→eh , N j→j h→eh , and N j→j+1 h→eh electrons in bins j − 1, j, and j + 1, respectively. These numbers are drawn from a multinomial distribution with event probabilities of p j→j−1 = 0.25, p j→j = 0.5, and p j→j+1 = 0.25, respectively, The updated charge carrier distributions are then computed as N j e (t i+1 ) = N j−1 e (t i ) + N j−1 e→eh (t i ) + N j+1→j h→eh (t i ) + N j→j h→eh (t i ) + N j−1→j h→eh (t i ), N j h (t i+1 ) = N j+1 h (t i ) + N j+1 h→eh (t i ) + N j+1→j e→eh (t i ) + N j→j e→eh (t i ) + N j−1→j e→eh (t i ).
Assuming a uniform weighting field E w /V w = 1/d, the current I ind (t i ) induced by the drifting charges is computed as where e 0 is the elementary charge.
In the continuum limit where α · ∆x → 0 and β · ∆x → 0, this simulation model converges to the avalanche model described in Section 2. For the results presented in the main body, the avalanche region is divided into N = 500 bins.