# On the convergence of the recurrence solution of McIntyre's local and non-local avalanche triggering probability equations for SPAD compact models

Dorian Saint-Pierre, Raphaël Clerc, Remi Helleboid, Denis Rideau

► To cite this version:

Dorian Saint-Pierre, Raphaël Clerc, Remi Helleboid, Denis Rideau. On the convergence of the recurrence solution of McIntyre's local and non-local avalanche triggering probability equations for SPAD compact models. ESSDERC 2022 - IEEE 52nd European Solid-State Device Research Conference (ESSDERC), Sep 2022, Milan, France. pp.277-280, 10.1109/ESSDERC55479.2022.9947121. hal-05416699

HAL Id: hal-05416699

https://hal.science/hal-05416699v1

Submitted on 15 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the convergence of the recurrence solution of McIntyre's local and non-local avalanche triggering probability equations for SPAD compact models

Dorian Saint-Pierre, Raphaël Clerc

Laboratoire Hubert Curien, UMR 5516 F-42023

Université de Lyon, UJM-Saint-Etienne, CNRS, Institut d'Optique

Graduate School, Saint-Etienne, France

raphael.clerc@institutoptique.fr

Rémi Helleboid, Denis Rideau

STMicroelectronics, 850 rue Jean Monnet, BP 16, F-38926

Crolles Cedex, France

denis.rideau@st.com

**Abstract**— McIntyre's local and non-local equations are used to calculate the avalanche triggering probability, useful in the modeling of Dark Count Rates and Probability Detection Efficiency of Single Photon Avalanche Detectors (SPAD). Although non-linear, the local equations have a closed-form solution, while the more rigorous non-local history-dependent equations can be solved only by an iterative approach. However, the convergence of the latest approach is fairly slow, making the modeling of SPAD properties relatively time consuming. An alternative and numerically efficient method is proposed to solve the local equations, avoiding the use of the closed-form solution and its time-consuming root finding step. Several approaches are then discussed to improve the convergence of the non-local history-dependent equations.

**Keywords**— SPAD, avalanche triggering probability, history-dependent equations, McIntyre.

## I. INTRODUCTION

Single Photon Avalanche Diodes (SPAD) are present both in Silicon and III-V technologies and are commonly used to detect time dependent and small optical power signals. They are found in several technologies such as lidars, time of flight 3D imaging, fluorescence lifetime microscopy, optical communications, and particle detection [1] [2].

Biased above the breakdown voltage, electrons and holes pairs created by band-to-band absorption of few photons may trigger an avalanche, leading to a measurable electrical current. It is thus important to determine accurately the avalanche triggering probability versus device design. This question has been addressed at the early stage of the development of this technology in the seventies in the pioneering works of Oldham [3] and McIntyre [4]. In these works, a set of non-linear differential equations was derived [3] and solved [4] in order to calculate the probability $P_p(x)$, i.e. the probability that an electron or a hole generated in x will initiate an avalanche. This model has been successfully used by many different authors to calculate the SPAD figures of merit, such as the Dark Count Rates and Probability Detection Efficiency [5] [6] [7] [8].

As pointed out by Okuto and Crowell [3] [4], and then Hayat [11], the local nature of Oldham and McIntyre formalism ignores the fact that a carrier needs to gain a given threshold energy before being capable of triggering an impact ionization. The typical distance (sometimes referred as dead space [11]) for reaching this threshold energy being of the order of magnitude of ten nanometers (depending of the field and the material), it may not be negligible in highly doped

devices, where the space charge layer is typically hundreds of nanometers long. To address this issue, McIntyre has proposed in 1999 [12] an improved model to account for history-dependent impact ionization in the calculation of avalanche triggering probability, leading to a new set of non-linear integral equations that can be solved by an iterative procedure. This improved approach has also been extensively used in the field of SPAD modeling [13] [14] [15].

Despite their 1D nature, both local and non-local models are relatively time consuming, making the use of these models difficult for compact model applications. As recognized by McIntyre himself "It was found that the convergence is fairly slow close to the breakdown voltage but rapid well above the breakdown voltage" [12]. This work investigates and proposes solutions to accelerate the calculation of both models.

The paper is organized as follow. In section 2, the main models and solution algorithms are summarized. Results and discussion are reported in section 3.

## II. NUMERICAL SOLUTION OF MCINTYRE LOCAL AND NON-LOCAL EQUATIONS

### A. Hypothesis and definitions

For the sake of simplicity, the device was modeled as a simple NP junction with constant doping profile throughout this paper. The electrical field and potential were calculated using the total depletion approximation. We believe that a more rigorous solution of the Poisson equation would have only increased the computational burden, without affecting the overall conclusions.

The space charge layer extends from $x = 0$ to $x = W$, and consequently generated electrons are moving from $W$ to 0 (negative $x$ direction) and holes from 0 to $W$ (positive $x$ direction). Following Oldham [16], we define $P_e(x)$ (resp. $P_h(x)$) as the probability that an electron (resp. a hole), generated in $x$ by any mechanism (light, trap, band to band or impact ionization), will initiate an avalanche (i.e. an infinity of successive impact ionization events). The probability that an electron (resp. a hole), generated in $x$, will not initiate an avalanche is referred as $P_{ne}(x)$ (resp. $P_{nh}(x)$), and is simply given as $P_{ne}(x) = 1 - P_e(x)$. The joint probability $P_p(x)$ that an electron or a hole generated in $x$ will initiate an avalanche is thus given by:

$$P_p(x) = 1 - (1 - P_e(x))(1 - P_h(x)) \quad (1)$$
## B. McIntyre local model and its closed-form solution

According to [4] and [16], $P_e(x)$ and $P_h(x)$ are solution of the coupled non-linear differential equations given by:

$$ \frac{dP_e}{dx} = \alpha_e(x)P_p(x)[1 - P_e(x)] \quad (2) $$

$$ \frac{dP_h}{dx} = -\alpha_h(x)P_p(x)[1 - P_h(x)] \quad (3) $$

Where $\alpha_e$ (resp. $\alpha_h$) are the electron (resp. hole) impact ionization coefficient. In this paper, the Reggiani's models for impact ionization for Silicon were used [17], [18].

If neither an electron generated at the entrance of the N region nor a hole at the entrance of the P region could initiate an avalanche, the boundary conditions are thus given by:

$$ P_e(0) = 0 \quad (4) $$

$$ P_h(W) = 0 \quad (5) $$

The conventional approach for solving these equations is detailed in [4] and leads to the following transcendental equation where $P_h(0)$ is the unknown variable:

$$ 1 = (1 - P_h(0)) \exp \left[ \int_{0}^{W} \alpha_h(u) P_p(u, P_h(0)) du \right] \quad (6) $$

This equation can only be solved numerically, which may be computationally heavy. In the framework of the total depletion approximation, we use the Maserjian's procedure [19] to calculate analytically the first integral, but this simplifying method cannot be used in the general case, where $E(x)$ does not have an analytic expression. Moreover, for applied voltage close to the breakdown voltage, such numerical solution requires accuracy as $P_h(0)$ is close to zero, leading to floating point issues.

Once $P_h(0)$ is found, other parameters such as $P_e(W)$ and $P_e(x)$ can be easily deduced from the integration of equations (2) and (3) (see reference [4] for details).

## C. An original iterative method to solve McIntyre local model:

The closed-form solution of McIntyre's local model remains numerically challenging, because of the root search step in Eq. (6). Inspired by the iterative approach of the McIntyre non-local model (see in the following for details), we suggest adopting a similar method to solve numerically the coupled equations (2) et (3). Using boundary conditions (4) and (5), these equations can be rearranged in the following integral form:

$$ P_e(x) = F(P_e, P_h, x) = \int_{0}^{x} \alpha_e(u) P_p(u) [1 - P_e(u)] du \quad (7) $$

$$ P_h(x) = G(P_e, P_h, x) = \int_{x}^{W} \alpha_h(u) P_p(u) [1 - P_h(u)] du \quad (8) $$

Such a set of equations can be solved by iteration. The first direct procedure (referred as “basic iteration”) consists in starting from a guess value $P_{e,0}(x)$, $P_{h,0}(x)$ and deducing solutions by several iterative steps according to the recurrence relations:

$$ P_{e,n+1}(x) = F(P_{e,n}, P_{h,n}, x) \quad (9) $$

$$ P_{h,n+1}(x) = G(P_{e,n}, P_{h,n}, x) \quad (10) $$

However, this brute force method is generally diverging for the local model. To avoid this problem, two other approaches have been also investigated. The first one consists in re-injecting only a small fraction $\beta$ of the solution (over-relaxation method) in the calculation of the n+1 step:

$$ P_{e,n+1}(x) = \beta F(P_{e,n}, P_{h,n}, x) + (1 - \beta) P_{e,n}(u) \quad (11) $$

$$ P_{h,n+1}(x) = \beta G(P_{e,n}, P_{h,n}, x) + (1 - \beta) P_{h,n}(u) \quad (12) $$

The second alternative method follows the principles of the well-known Newton-Raphson procedure. A first order power expansion is used:

$$ P_{e,n+1} = F(P_{e,n}, P_{h,n}) + \frac{\partial F}{\partial P_e}(P_{e,n+1} - P_{e,n}) + \frac{\partial F}{\partial P_h}(P_{h,n+1} - P_{h,n}) \quad (13) $$

$$ P_{h,n+1} = G(P_{e,n}, P_{h,n}) + \frac{\partial G}{\partial P_e}(P_{e,n+1} - P_{e,n}) + \frac{\partial G}{\partial P_h}(P_{h,n+1} - P_{h,n}) \quad (14) $$

Interestingly, each term of the Jacobian matrix can be analytically derived, leading to:

$$ A = \frac{\partial F}{\partial P_e} = \int_{0}^{x} \left[ \alpha_e \left( 1 - 2P_{p,n}(u) \right) \right] du \quad (15) $$

$$ B = \frac{\partial F}{\partial P_h} = \int_{0}^{x} \left[ \alpha_e \left( 1 - P_{e,n}(u) \right)^2 \right] du \quad (16) $$

$$ C = \frac{\partial G}{\partial P_e} = \int_{x}^{w} \left( \alpha_h(u) (1 - P_h(u))^2 \right) du \quad (17) $$

$$ D = \frac{\partial G}{\partial P_h} = \int_{x}^{w} \left( \alpha_h(u) (1 - 2P_{p,n}(u)) \right) du \quad (18) $$

Equations (13) and (14) can be re-arranged by inverting the Jacobian matrix, leading to:

$$ P_{e,n+1} = \frac{(D + 1)F_n - BG_n + (A(D + 1) - BC)P_{e,n} + BP_{h,n}}{A + D + AD - BC + 1} \quad (19) $$

$$ P_{h,n+1} = \frac{(A + 1)G_n - C F_n + C P_{e,n} + ((A + 1)D - BC)P_{h,n}}{A + D + AD - BC + 1} \quad (20) $$

Results obtained using these three approaches to solve local McIntyre equations will be discussed in the “results and discussions” section.

## D. Iterative solution of McIntyre non-local model:

McIntyre proposed a model to account for history dependent ionization probability in the calculations of avalanche triggering probability [12]. In this approach, the probability that an electron (resp. a hole), generated in x, will not initiate an avalanche is referred as $P_{ne}(x)$ (resp. $P_{nh}(x)$) is given by:

$$ P_{ne}(x) = F(P_{ne}, P_{nh}, x) \quad (21) $$

$$ = P_{se}(x|0) + \int_{0}^{x} p_{e}(x|u) P_{ne}(u)^2 P_{nh}(u) du \quad (21) $$

$$ P_{nh}(x) = G(P_{ne}, P_{nh}, x) \quad (21) $$

$$ = P_{sh}(x|W) + \int_{x}^{W} p_{h}(x|u) P_{nh}(u)^2 P_{ne}(u) du \quad (21) $$
$p_e(x'|x)$ is the probability of first ionization by an electron and is given by:

$$p_e(x'|x) = \alpha(x'|x) P_{se}(x'|x) \quad (22)$$

Where $P_{se}(x'|x)$ is the survival probability, i.e., the probability that the electron will travel from $x'$ to $x$ without undergoing an ionizing collision:

$$P_{se}(x'|x) = \exp\left(-\int_{x}^{x'} \alpha(x'|u)du\right) \quad (23)$$

Several approaches have been proposed to deduce the non-local ionizing coefficient $\alpha(x'|x)$ from the local one, based on the dead space model [20] or the introduction of effective fields [12],[21]. Here we used the effective field definition proposed by C. Nichetti et al. [21], relying on the formalism of energy balance equations. The mean free path parameters $\lambda_e$ for electrons and $\lambda_h$ for holes have been taken constant and equal to 20 nm for simplicity. It has been numerically checked that local and non-local models give identical solutions when $\lambda_e$ and $\lambda_h$ tends to zero.

As in the previous section for the local model, three approaches (basic iteration, over-relaxation, and Newton Raphson) have been implemented to achieve the convergence of the iterative procedure. Regarding the Newton Raphson method, the same equations (19) (20) can be used, replacing the Jacobian matrix coefficients A, B, C, D by the appropriate expressions.

### III. RESULTS AND DISCUSSION

In this section, the convergence of both the local and non-local iterative models are examined. The error between two consecutive steps is calculated with the r.m.s. error.

Fig. 1 shows the convergence of the three methods implemented for the iterative local equations for two different guess value. A moderately doped abrupt PN junction ($N_d = 7. 10^{17} cm^{-3}$, $N_a = 3 10^{17} cm^{-3}$) was considered. The first guess value is given by:

$$P_{e,0}(x) = \frac{x}{W} \quad \text{and} \quad P_{h,0}(x) = 1 - \frac{x}{W}$$

The second guess value is the exact solution obtained by direct solution of the local McIntyre equations (see section II.B).

<figure><img src="image_2.png" /><figcaption>Fig. 1. Error between successive step versus number of iterations solving the local McIntyre equations, using a linear approximation as a guess value (Guess value 1, left figure) or the closed form solution (Guess value 2, right figure). Doping levels are $N_d = 7. 10^{17} cm^{-3}$, $N_a = 3 10^{17} cm^{-3}$, the excess voltage is 5 V, $\lambda_e = \lambda_h = 20$ nm.</figcaption></figure>

The overall error at the initial step is lower in the second case (exact guess value) than in the first (linear approximation). The basic scheme is diverging in the first case, while the over-relaxation scheme ($\beta = 0.2$) is slowly

converging. The Newton Raphson scheme however shows a fast and smooth error decrease versus number of iterations, reaching an error of $10^{-5}$ after only 15 iterations in the first case, and 8 iterations in the second. Similar results were obtained for all the doping and excess voltage considered, showing that the Newton-Raphson implementation is an extremely efficient way of solving the McIntyre local equations.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Error between successive step versus number of iterations solving the non-local McIntyre equations, using a linear approximation as a guess value (Guess value 1, left figure) or the closed form solution (Guess value 2, right figure).</figcaption></figure>

Conclusions significantly differ when non-local equations are considered (see Fig. 2). In this case, the basic iteration procedure is the most efficient scheme independently of the guess value. The same guess values than in the previous case were used, except that the local direct solution is no longer the exact solution of the non-local equations but only a rough approximation. Contrary to the local case, the change between each step of the non-local iterations is very small, preventing divergence. In this situation, the over-relaxation and Newton-Raphson only slow down the convergence. A “hybrid” method was also investigated, which consists in applying a Newton-Raphson step every 3 steps of basic iteration.

To compare, the correct way of counting a step is to count 4 (one Newton-Raphson step and 3 basic iteration steps) for one hybrid step: this counting method has been used for “Hybrid 2” in Fig. 2, while an incorrect counting (which consists in counting as one step for one Newton-Raphson step and 3 basic iterations steps) would artificially give the impression that the Hybrid method is more efficient (“Hybrid 1” in Fig. 2). Using the proper way of counting iterations (“Hybrid 2”), the hybrid method does not show any convergence improvement compared to the basic method.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Number of iterations required to reach an error lower than 5. 10-5 between two successive steps using the non-local model versus excess voltage : using the linear solution as a guess value (GV: Linear solution), or using the solution at the previous voltage as guess value, starting from the maximum voltage (GV: previous voltage 1) or starting from the minimum voltage, (GV : previous voltage 2).</figcaption></figure>

As in most application, avalanche triggering probabilities are calculated for several applied voltages, we also investigate
the possibility to use the result obtained at the previous applied voltage as guess value for the next one, instead of the solution of the linear equations, as usually done. It turns out that this approach is more efficient, as seen in Fig. 3. Moreover, in this case, it is recommended to start the simulation by the highest applied voltage, and not the lowest, as low voltage simulations, whatever the guest value, are much more time consuming.

## IV. CONCLUSIONS

Several original findings about the numerical solution of the local and non-local McIntyre avalanche triggering probability in SPAD device have been discussed in this paper. First, it appears that the local equations can be very efficiently solved by a Newton-Raphson algorithm. The Jacobian matrix was calculated and inverted analytically, accelerating calculations. We believe that this approach is easier to implement and more efficient than the closed form solution derived by McIntyre.

Regarding the non-local equations, contrary to the local case, the basic implementation of the iterative solutions was found the most efficient approach compared to over-relaxation, Newton-Raphson or even hybrid method. Moreover, this implementation can be accelerated using the previous excess voltage as guess value (and not the local solution as it is usually done). In this case, it is recommended to start with the maximum voltage (where the local approximation is a good guess value) and not with the minimum voltage.

Finally, the local equations have been found to be a good approximation of the non-local one at high excess voltage. Further investigations are required to confirm if the local approximation is sufficiently accurate for all practical case of application in real device design.

## ACKNOWLEDGMENT

The authors would like to thank A. Pilotto and P. Palestri from university of Udine (Italy) for helpful discussions. This work has been supported by the IPCEI Nano 2022 program and by the French ANR through the project GeSPAD under Grant ANR-20-CE24-0004.

## REFERENCES

[1] M. Ghioni, A. Gulinatti, I. Rech, F. Zappa, and S. Cova, "Progress in silicon single-photon avalanche diodes," IEEE J. Sel. Top. Quantum Electron., vol. 13, no. 4, pp. 852-862, 2007, doi: 10.1109/JSTQE.2007.902088.

[2] G. S. Buller and R. J. Collins, "Single-photon generation and detection," Meas. Sci. Technol., vol. 21, no. 1, 2010, doi: 10.1088/0957-0233/21/1/012002.

[3] W. G. Oldham, R. R. Samuelson, and P. Antognetti, "Triggering Phenomena in Avalanche Diodes," IEEE Trans. Electron Devices, vol. 19, no. 9, pp. 1056-1060, 1972, doi: 10.1109/T-ED.1972.17544.

[4] R. J. McIntyre, "On the Avalanche Initiation Probability of Avalanche Diodes Above the Breakdown Voltage," IEEE Trans. Electron Devices, vol. 20, no. 7, pp. 637-641, 1973, doi: 10.1109/T-ED.1973.17715.

[5] T. C. De Albuquerque et al., "An analytical solution for McIntyre's model of avalanche triggering probability for SPAD compact modeling and performance exploration," Semicond. Sci. Technol., vol. 36, no. 8, 2021, doi: 10.1088/1361-6641/ac00d0.

[6] M. M. Vignetti, F. Calmon, P. Lesieur, and A. Savoy-Navarro, "Simulation study of a novel 3D SPAD pixel in an advanced FD-SOI technology," Solid. State. Electron., vol. 128, pp. 163-171, 2017, doi: 10.1016/j.sse.2016.10.014.

[7] A. Gulinatti, I. Rech, M. Assanelli, M. Ghioni, and S. Cova, "A physically based model for evaluating the photon detection efficiency and the temporal response of SPAD detectors," J. Mod. Opt., vol. 58, no. 3-4, pp. 210-224, 2011, doi: 10.1080/09500340.2010.536590.

[8] A. Panglosse, P. Martin-Gonthier, O. Marcelot, C. Virmontois, O. Saint-Pé, and P. Magnan, "Modeling, simulation methods and characterization of photon detection probability in CMOS-SPAD," Sensors, vol. 21, no. 17, pp. 1-17, 2021, doi: 10.3390/s21175860.

[9] Y. Okuto and C. R. Crowell, "Energy-conservation considerations in the characterization of impact ionization in semiconductors," Phys. Rev. B, vol. 6, no. 8, pp. 3076-3081, 1972, doi: 10.1103/PhysRevB.6.3076.

[10] Y. Okuto and C. R. Crowell, "Threshold energy effect on avalanche breakdown voltage in semiconductor junctions," Solid State Electron., vol. 18, no. 2, pp. 161-168, 1975, doi: 10.1016/0038-1101(75)90099-4.

[11] M. M. Hayat, W. L. Sargeant, and B. E. A. Saleh, "Effect of Dead Space on Gain and Noise in Si and GaAs Avalanche Photodiodes," IEEE J. Quantum Electron., vol. 28, no. 5, pp. 1360-1365, 1992, doi: 10.1109/3.135278.

[12] R. J. McIntyre, "A new look at impact ionization-part I: A theory of gain, noise, breakdown probability, and frequency response," IEEE Trans. Electron Devices, vol. 46, no. 8, pp. 1623-1631, 1999, doi: 10.1109/16.777150.

[13] G. J. Rees and J. P. R. David, "Nonlocal impact ionization and avalanche multiplication," J. Phys. D. Appl. Phys., vol. 43, no. 24, 2010, doi: 10.1088/0022-3727/43/24/243001.

[14] J. S. Ng, C. H. Tan, and J. P. R. David, "A comparison of avalanche breakdown probabilities in semiconductor materials," J. Mod. Opt., vol. 51-9, no. 10, pp. 1315-1321, 2004, doi: 10.1080/09500340408235274.

[15] M. M. Hayat et al., "Breakdown probabilities for thin heterostructure avalanche photodiodes," IEEE J. Quantum Electron., vol. 39, no. 1, pp. 179-185, 2003, doi: 10.1109/JQE.2002.806217.

[16] W. G. Oldham, P. Antognetti, and R. R. Samuelson, "New method for breakdown voltage determination in p-n junctions," Appl. Phys. Lett., vol. 19, no. 11, pp. 466-467, 1971, doi: 10.1063/1.1653774.

[17] S. Reggiani, E. Gnani, A. Gnudi, M. Rudan, S. Member, and G. Baccarani, "Low-Field Electron Mobility Model for Ultrathin-Body SOI and Double-Gate MOSFETs With Extremely Small Silicon Thicknesses," vol. 54, no. 9, pp. 2204-2212, 2007.

[18] M. Rudan, R. Katilius, S. Reggiani, E. Gnani, and G. Baccarani, "Impact-ionization coefficient in silicon at high fields - A parametric approach," J. Comput. Electron., vol. 7, no. 3, pp. 151-154, 2008, doi: 10.1007/s10825-008-0184-8.

[19] J. Maserjian, "Determination of avalanche breakdown in pn junctions [4]," J. Appl. Phys., vol. 30, no. 10, pp. 1613-1614, 1959, doi: 10.1063/1.1735012.

[20] M. M. Hayat and B. E. A. Saleh, "Statistical Properties of the Impulse Response Function of Double-Carrier Multiplication Avalanche Photodiodes Including the Effect of Dead Space," J. Light. Technol., vol. 10, no. 10, pp. 1415-1425, 1992, doi: 10.1109/50.166785.

[21] C. Nichetti et al., "An improved nonlocal history-dependent model for gain and noise in avalanche photodiodes based on energy balance equation," IEEE Trans. Electron Devices, vol. 65, no. 5, pp. 1823-1829, 2018, doi: 10.1109/TED.2018.281