# Time resolution and efficiency of SPADs and SiPMs for photons and charged
particles

W. Riegler  P. Windischhofer  CERN  University of Oxford

###### Abstract

We give an analytic treatment of the time resolution and efficiency of Single
Photon Avalanche Diodes (SPADs) and Silicon Photomultipliers (SiPMs). We
provide closed-form expressions for structures with uniform electric fields
and efficient numerical prescriptions for arbitrary electric field
configurations. We discuss the sensor performance for single photon detection
and also for charged particle detection.

##  1 Introduction

Single Photon Avalanche Diodes (SPADs) have been used for photon detection and
photon counting since many decades. These semiconductor devices contain a
highly doped p-n junction of $0.5{-}2\,\mu$m thickness, the so called gain
region, that is biased above the breakdown voltage. This means that a single
primary electron or hole entering this region can produce a diverging
avalanche through impact ionization and therefore lead to a detectable signal.
The growth of the diverging avalanche is quenched by the breakdown of the
electric field, which leads to a digital type signal with an amplitude
independent of the number of primary charges. The restoration of the electric
field is governed by a quench resistor or quench circuit external to the
device.

Figure 1: a) p-in-n SPAD with a gain layer only. b) n-in-p SPAD with a gain
layer and a conversion layer.

Two basic SPAD geometries are shown in Fig. 1. For the p-in-n SPAD in Fig. 1a
the high field region is placed close to the surface of the silicon. Photons
with a wavelength of $<500$ nm have an absorption length of less than
$1\,\mu$m, so they will be efficiently absorbed in the gain region and trigger
the breakdown.
For photons with longer wavelengths and therefore larger absorption length,
the geometry of Fig. 1b with a so called ’conversion region’ will be more
efficient. This region consists of a depleted layer of silicon with thickness
ranging from 10-100 $\mu$m that is adjacent to a gain region. The electrons
created in the conversion region will drift to the gain region where they
provoke the breakdown. The sensor can be illuminated from the front and the
back side. These SPADs are nowadays heavily used for LIDAR [1] applications
with wavelengths in the near-infrared above 750 nm.

While originally being a single channel device, the advances in silicon
industry allowed the arrangement of many of these diodes in chessboard
structures. Depending on the readout circuitry these are called SPAD image
sensors or Silicon Photomultipliers (SiPMs). In SPAD image sensors, the SPAD
pixels are read out individually in order to form an image. In SiPMs the
individual pixels are connected in parallel through series resistors into a
single channel. The amplitude of the output signal will then be proportional
to the total number of fired pixels and therefore proportional to the number
photons that have hit the pixel matrix, which represents the function of a
traditional photomultiplier. There are many technological challenges related
to the implementation of such pixel structures, specifically the elimination
of crosstalk between the channels and the minimization of the dark count rate.

Scope and outline
In this report, we discuss the time resolution and efficiency of SPADs and
SiPMs of the two types shown in Fig. 1 and evaluate their performance for the
detection of photons as well as charged particles. Our results are derived
from a series of fundamental equations that describe the movement and the
interactions of the participating charge carriers.
We provide analytic expressions for the simplified structures in Fig. 2 with
constant electric fields. Fig. 2a represents the situation where a photon
interacts directly inside a gain layer of constant electric field $E_{1}$ that
is above the breakdown value. Fig. 2b shows a SPAD where we represent the
conversion region by a silicon layer of thickness $w$ with constant electric
field $E_{0}$ below the breakdown field, which is adjacent to a gain layer
with constant field $E_{1}$.

Figure 2:  Two simplified SPAD structures with constant electric fields
discussed in this report: a) A photon interacts directly in the gain layer and
produces an e-h pair that provokes breakdown. b) A photon produces an e-h pair
in the conversion layer and the electron moves to the gain layer where it
provokes breakdown. A thin gain layer of 0.5-2 $\mu$m is also highly efficient
for charged particle detection.

This report is structured as follows: Section 2 discusses the absorption of
photons in silicon, and the contribution of the conversion layer to the time
resolution. Section 3 then discusses the mechanisms of avalanche breakdown and
the average growth of the avalanche. Section 4 describes the contribution to
the efficiency resulting from the avalanche formation in the gain layer.
Section 5 discusses avalanche fluctuations and their contribution to the time
resolution. This discussion is based to a large degree on the companion paper
[2], which develops the statistics of electron-hole avalanches in great
detail. Section 6 then discusses the performance of SPADs for the detection of
charged particles. We finally drop the assumption of constant electric field
in Section 7 and return to the realistic field configurations of Fig. 1. We
give efficient numerical prescriptions that extend the analytic results
obtained previously.
Details of all calculations are given in several appendices. Although we focus
specifically on devices based on silicon, our results are expected to also
cover the basic geometries for different types of semiconductors.

##  2 Conversion layer

We assume a layer of silicon of thickness $w$ extending from $x=0$ to $x=w$ as
shown in Fig. 3. The probability for a photon to be absorbed between position
$x_{0}$ and $x_{0}+dx_{0}$ is given by
$P(x_{0})dx_{0}=1/l_{a}e^{-x_{0}/l_{a}}dx_{0}$, where $l_{a}$ is the photon
absorption length from Fig. 4a. The efficiency, i.e. the probability for a
photon to convert in the layer, is then given by $p=1-e^{-w/l_{a}}$ and the
numbers are shown in Fig. 4b. Photons of wavelength $<500$ nm are efficiently
absorbed in $<1\,\mu$m of silicon while infrared photons of ${>}750$ nm need
several tens of $\mu$m of silicon to be absorbed efficiently.

Figure 3: The conversion layer of a SPAD or SiPM. A photon is absorbed at
position $x_{0}$ producing an e-h pair. The electron will drift to the gain
layer at $x=w$.

a)  b)

Figure 4:  Absorption length $l_{a}$ for photons of different wavelengths in
silicon [3, 4]. b) Photon absorption efficiency for different values of
silicon thickness.

First we investigate the case where the photon is arriving from the ’left’
side as shown in Fig. 3a. Normalising the conversion probability to the
efficiency, the probability for a photon to be absorbed between position
$x_{0}$ and $x_{0}+dx_{0}$ is

$P(x_{0})dx_{0}=\frac{1}{1-e^{-w/l_{a}}}\frac{1}{l_{a}}e^{-x_{0}/l_{a}}\Theta(w-x_{0})dx_{0}$
(1)

The electron is then moving to the edge of the silicon layer at $x=w$ with a
velocity $v_{e}$, where it arrives at time $t=(w-x_{0})/v_{e}$. The velocity
of electrons and holes in silicon is shown in Fig. 5a and the parametrization
is given in Appendix A.

a)  b)

Figure 5: a) Drift velocity of electrons ($v_{e}$) and holes ($v_{h}$) as a
function of electric field in silicon. The velocity
$v^{*}=2v_{e}v_{h}/(v_{e}+v_{h})$ that is relevant for the avalanche growth in
the gain layer is shown as well. b) Probability for the electron to arrive at
$x=w$ between times $t$ and $t+dt$ for $w=10\,\mu$m, for a photon with
$\l_{a}=1\,\mu$m entering the layer from the left side. The ’no diffusion’
curve refers Eq. 2 and the ’diffusion’ curve refers to Eq. 6.

The arrival time distribution of the electron at $x=w$ is therefore

$\rho_{1}(t)=\int_{0}^{w}P(x_{0})\delta[t-(w-x_{0})/v_{e}]dx_{0}=\frac{w}{l_{a}(e^{w/l_{a}}-1)}\frac{1}{T}e^{\frac{tw}{l_{a}T}}\,\Theta(T-t)$
(2)

where we have expressed the velocity $v_{e}$ by the maximum drift time
$T=w/v_{e}$ of the electrons inside the conversion layer. An example is shown
in Fig. 5b. The variance of the arrival time is then

$\sigma_{t}^{2}=\int_{0}^{T}t^{2}\,P_{2}(t)dt-\left(\int_{0}^{T}t\,P_{2}(t)dt\right)^{2}=T^{2}\left(\frac{l_{a}^{2}}{w^{2}}-\frac{1}{4\sinh[w/2\l_{a}]^{2}}\right)$
(3)

Including the effect of diffusion we use the fact that an electron deposited
at position $x_{0}$ at $t=0$ and moving with an average velocity of $v_{e}$
will be found at position $x$ after a time $t$ with a probability of

$p(x,x_{0},t)dx=\frac{1}{\sqrt{2\pi}\sqrt{2Dt}}\exp\left[-\frac{(x-(x_{0}+v_{e}t))^{2}}{2(2Dt)}\right]dx$
(4)

The standard deviation of the distribution is given by $\sigma(t)=\sqrt{2Dt}$.
The probability for an electron to arrive at $x=w$ between time $t$ and $t+dt$
is then $p(w,x_{0},t)v_{e}dt$. Since the probability of a photon absorption at
$x_{0}$ is given by $P(x_{0})$ from Eq. 1 we find the arrival time
distribution $P_{3}(t)$ as

$\rho_{2}(t)=\int_{0}^{w}P(x_{0})p(w,x_{0},t)v_{e}dx_{0}=$ (5)
$\frac{v_{e}}{2l_{a}}\frac{1}{e^{w/l_{a}}-1}e^{(D+l_{a}v_{e})t/l_{a}^{2}}\left(\mbox{Erf}\left[\left(\frac{\sqrt{Dt}}{l_{a}}+\frac{v_{e}}{2}\sqrt{\frac{t}{D}}\right)\right]-\mbox{Erf}\left[\frac{\sqrt{Dt}}{l_{a}}+\frac{v_{e}}{2}\sqrt{\frac{t}{D}}-\frac{w}{2\sqrt{Dt}}\right]\right)$
(6)

with Erf$(z)=2/\sqrt{\pi}\int_{0}^{z}e^{-t^{2}}dt$. The variance evaluates to

$\sigma_{t}^{2}=T^{2}\left(\frac{l_{a}^{2}}{w^{2}}-\frac{1}{4\sinh[w/2l_{a}]^{2}}\right)+T\frac{2D}{v_{e}^{2}}\left(\frac{1}{1-e^{-w/l_{a}}}-\frac{l_{a}}{w}\right)+\frac{8D^{2}}{v_{e}^{4}}$
(7)

The first term is the one from Eq. 3 due to the varying position of the photon
interaction together with the average drift time, the second term is due to
diffusion. The third term is an artefact of the assumption that a charge
placed at $x=0$ can diffuse into a region of $x<0$, so the term does not
vanish even for $w=0$. In a realistic implementation of a conversion layer,
the region of $x<0$ in Fig. 3 will represent a region where the electric field
drops sharply to low values. The electron might spend a rather long time in
this region before moving back into the high field region or it might even get
lost. We therefore count the cases where the electron moves to this area as
additional small inefficiency and neglect the term. In practice, photons
interacting at the boundary of the conversion layer will cause a long tail in
the time distribution and it becomes a practical question whether to include
the tails in the calculation of the time resolution or count them as
inefficiency. In the limit of large and small values of $l_{a}/w$ the time
resolution approximates to

$\sigma_{t}^{2}=\frac{T^{2}}{12}+\frac{DT}{v_{e}^{2}}\quad l_{a}\gg
w\qquad\sigma_{t}^{2}=T^{2}\frac{l_{a}^{2}}{w^{2}}+\frac{2DT}{v_{e}^{2}}\quad
l_{a}\ll w$ (8)

We assume $D=35$ cm${}^{2}/$s for electrons in silicon. For electric fields in
excess of $5\times 10^{4}$ V/cm the electron velocity is close to saturation
and we have $D/(v_{sat}^{e})^{2}=0.35$ ps. For a conversion layer of
$w=1/10/100\,\mu$m the saturated drift time is $T\approx 10/100/1000$ ps.
For $\l_{a}\gg w$ the probability for the photon conversion position becomes
uniform across $w$, diffusion is negligible and the time resolution is
$\sigma_{t}=T/\sqrt{12}\approx 2.89/28.9/289$ ps for $w=1/10/100\,\mu$m. For
$\l_{a}\ll w$ the photon conversion point is always close to $x=0$, diffusion
will dominate and the time resolution is equal to
$\sigma_{t}=\sqrt{2DT}/v_{e}\approx 1.6/8.37/26.46$ ps for $w=1/10/100\,\mu$m.

If the illumination takes place from the ’right’ side as indicated in Fig. 3b,
the time resolution becomes

$\sigma_{t}^{2}=T^{2}\left(\frac{l_{a}^{2}}{w^{2}}-\frac{1}{4\sinh[w/2l_{a}]^{2}}\right)+T\frac{2D}{v_{e}^{2}}\left(\frac{l_{a}}{w}+\frac{1}{1-e^{w/l_{a}}}\right)+\frac{8D^{2}}{v_{e}^{4}}$
(9)

This expression differs from the previous one only by the diffusion term,
because the distance between the conversion point and $x=w$ is now different.
In the limit of large and small values of $l_{a}$ we have

$\sigma_{t}^{2}=\frac{T^{2}}{12}+\frac{DT}{v_{e}^{2}}\quad l_{a}\gg
w\qquad\sigma_{t}^{2}=T^{2}\frac{l_{a}^{2}}{w^{2}}+\frac{2DT}{v_{e}^{2}}\frac{l_{a}}{w}\quad
l_{a}\ll w$ (10)

For $l_{a}\gg w$ the expression is equal to the one from above, while for
$l_{a}\ll w$ the variance goes to zero because the conversion point is close
to $x=w$. The expression does of course not apply if the absorption length
$l_{a}$ is of the same order or smaller than the gain layer thickness $d$,
because the photons will interact directly in the gain layer.

##  3 Electron-hole avalanches and breakdown, average signal

Figure 6: The primary electrons and holes deposited at $x=x_{0}$ are
multiplying, which results in a diverging avalanche in case the electric field
$E_{1}$ is above the breakdown limit.

An electron drifting inside the conversion layer will move to the gain layer
and trigger an avalanche starting from $x=0$. Alternatively a photon can
convert inside the gain layer and the e-h pair at position $x=x_{0}$ will
trigger the avalanche. To cover both situations, we treat the general case
where $n^{0}_{e}$ electrons and $n^{0}_{h}$ holes are deposited at $x=x_{0}$
at time $t=0$, as shown in Fig. 6.
To derive equations describing the avalanche, we allow for general position-
dependent electric fields $E_{1}(x)$. With the field orientated as shown in
Fig. 6, electrons move to the right and holes move to the left with velocities
$v_{e}(x),v_{h}(x)$. The probability for an electron to create an e-h pair
when travelling a distance $dx$ is $\alpha(x)dx$ while the probability for a
hole to produce an e-h pair over distance $dx$ is $\beta(x)dx$, where
$\alpha(x)$ and $\beta(x)$ are called the impact ionization coefficients. The
values for silicon (Fig. 7) are reported in [7] with the parameters listed in
Appendix A. Since $1/\alpha$ and $1/\beta$ refer to the average distance that
an electron or a hole has to travel in order to produce one additional e-h
pair, we see that only for fields in excess of $2{-}3\times 10^{5}$ V/cm there
is an appreciable probability to provoke an avalanche in a few $\mu$m of
silicon.

a)  b)

Figure 7:  a) Impact ionization coefficient $\alpha$ for electrons and $\beta$
for holes as a function of the electric field. b) Minimum electric field value
provoking breakdown for a gain layer with constant electric field across a
given thickness $d$.

Fig. 8 shows a Monte Carlo (MC) simulation of a few avalanches starting with a
single electron. After some initial fluctuations the avalanche just grows
exponentially. There is also a finite probability that no diverging avalanche
develops.

Figure 8: Monte Carlo simulation for an electron-hole avalanche starting with
a single electron at $x_{0}=0$ for a $1\,\mu$m diode at a field of 4 V/$\mu$m.
After some initial fluctuations the number of charge carriers increases
exponentially. The dominant term from Eq. 30 approximates the average signal
extremely well for times $t>d/v^{*}\approx 11$ ps. When the avalanche is still
small there is a finite probability that no breakdown occurs, as is the case
for the avalanche MC4.

We denote as $n_{e}(x,t)dx$ and $n_{h}(x,t)dx$ the average number of electrons
and holes between position $x$ and $x+dx$ at time $t$ (note that this is
different from the notation used in [2]). These charge densities result in
local average current densities of $j_{e}(x,t)=v_{e}(x)n_{e}(x,t)$ and
$j_{h}(x,t)=-v_{h}(x)n_{h}(x,t)$. By the continuity equation $\partial
j/\partial x+\partial n/\partial t=\sigma$, with $\sigma$ the generation rate,
we therefore have

$\displaystyle\frac{\partial n_{e}(x,t)}{\partial t}+\frac{\partial
v_{e}(x)\,n_{e}(x,t)}{\partial x}$ $\displaystyle=$
$\displaystyle\alpha(x)v_{e}(x)n_{e}(x,t)+\beta(x)v_{h}(x)n_{h}(x,t)$ (11)
$\displaystyle\frac{\partial n_{h}(x,t)}{\partial t}-\frac{\partial
v_{h}(x)\,n_{h}(x,t)}{\partial x}$ $\displaystyle=$
$\displaystyle\alpha(x)v_{e}(x)n_{e}(x,t)+\beta(x)v_{h}(x)n_{h}(x,t)$

The fact that electrons move to the left and holes move to the right gives the
boundary conditions

$n_{e}(0,t)=0\qquad n_{h}(d,t)=0$ (12)

Since Eq. 11 represents a set of linear equations we can use the Ansatz
$n_{e}(x,t)=f(x)e^{St}$ and $n_{h}(x,t)=g(x)e^{St}$ and we find

$\displaystyle Sf(x)+[v_{e}(x)f(x)]^{\prime}$ $\displaystyle=$
$\displaystyle\alpha(x)v_{e}(x)f(x)+\beta(x)v_{h}(x)g(x)$ (13) $\displaystyle
Sg(x)-[v_{h}(x)g(x)]^{\prime}$ $\displaystyle=$
$\displaystyle\alpha(x)v_{e}(x)f(x)+\beta(x)v_{h}(x)g(x)$

with $f(0)=0$ and $g(d)=0$. The multiplication of electrons and holes can lead
to a finite amount of total charge in the avalanche ($S<0$) or it can diverge
and cause breakdown ($S>0$). The boundary between the two regimes is at $S=0$,
so by setting $S=0$ in the above equations and solving them with the given
boundary conditions we find the breakdown condition (Appendix B)

$\int_{0}^{d}\alpha(x)\,\exp\left[-\int_{0}^{x}(\alpha(x^{\prime})-\beta(x^{\prime}))dx^{\prime}\right]dx=1$
(14)

The breakdown condition is independent of $v_{e}(x)$ and $v_{h}(x)$. This
relation is usually derived by evaluating the point at which the gain for a
constant current injected into the gain layer diverges [8]. Evaluating the
breakdown equation for constant $\alpha$ and $\beta$ implies that breakdown
occurs if

$d>\frac{1}{\alpha-\beta}\ln\frac{\alpha}{\beta}$ (15)

The electric field at which breakdown takes place for a gain layer with a
given thickness $d$ is shown in Fig. 7b. For general electric field profiles,
Eqs. 11 can be efficiently solved with numerical methods. The solution can be
given in analytical form for constant values of $\alpha,\beta,v_{e},v_{h}$ and
is derived in [2]. It is represented as an infinite sum of exponential terms
with (generally complex valued) time constants. At least one time constant is
guaranteed to be real-valued. The largest real-valued time constant defines
the long-term behaviour of the avalanche. Above the breakdown limit, this term
determines the rate of exponential growth of the avalanche. Starting with
$n_{e}^{0}$ electrons and $n_{h}^{0}$ holes at position $x_{0}$ at time $t=0$,
it reads

$n_{e}(x,t)=\frac{1}{d}\,a(x)\left[n_{e}^{0}u_{e}(x_{0})+n_{h}^{0}u_{h}(x_{0})\right]\,e^{\gamma
v^{*}t}\qquad
n_{h}(x,t)=\frac{1}{d}\,b(x)\left[n_{e}^{0}u_{e}(x_{0})+n_{h}^{0}u_{h}(x_{0})\right]\,e^{\gamma
v^{*}t}$ (16)

with

$\displaystyle u_{e}(x_{0})$ $\displaystyle=$ $\displaystyle
e^{-a_{1}x_{0}}\sin(k-kx_{0}/d)$ (17) $\displaystyle u_{h}(x_{0})$
$\displaystyle=$ $\displaystyle\,\frac{e^{-a_{1}x_{0}}}{\alpha
d}\left[k\cos\left(k-kx_{0}/d\right)+\lambda_{1}\sin(k-kx_{0}/d)\right]$ (18)
$\displaystyle a(x)$ $\displaystyle=$
$\displaystyle\frac{2v_{h}ke^{a_{1}x}\sin\left(k\frac{x}{d}\right)}{(v_{e}+v_{h})(1+\lambda_{1})\sin
k}$ (19) $\displaystyle b(x)$ $\displaystyle=$
$\displaystyle\frac{2kv_{e}e^{a_{1}x}\left[k\cos\left(k\frac{x}{d}\right)+\lambda_{1}\sin\left(k\frac{x}{d}\right)\right]}{(v_{e}+v_{h})\beta
d(1+\lambda_{1})\sin k}$ (20)

and the constants $v^{*},a_{1},\gamma,k$ are defined by

$\displaystyle v^{*}$ $\displaystyle=$
$\displaystyle\frac{2v_{e}v_{h}}{v_{e}+v_{h}}$ (21) $\displaystyle a_{1}$
$\displaystyle=$ $\displaystyle\frac{\alpha v_{e}-\beta
v_{h}}{v_{e}+v_{h}}+\frac{v_{e}-v_{h}}{v_{e}+v_{h}}\frac{1}{d}\ \,\lambda_{1}$
(22) $\displaystyle\gamma$ $\displaystyle=$
$\displaystyle\frac{\alpha+\beta}{2}+\frac{\lambda_{1}}{d}$ (23)
$\displaystyle k$ $\displaystyle=$ $\displaystyle\sqrt{\alpha\beta
d^{2}-\lambda_{1}^{2}}$ (24)

The parameter $\lambda_{1}$ is the largest real solution of the equation

$\lambda_{1}+\sqrt{\alpha\beta d^{2}-\lambda_{1}^{2}}\cot\sqrt{\alpha\beta
d^{2}-\lambda_{1}^{2}}=0$ (25)

It holds that $-\infty<\lambda_{1}<d\sqrt{\alpha\beta}$. For
$\lambda_{1}<-d\sqrt{\alpha\beta}$ the constant $k$ will become imaginary,
which will still lead to real valued expressions for $n_{e}(x,t)$ and
$n_{h}(x,t)$ with $\sin,\cos$ becoming $\sinh,\cosh$. The functional form of
$a(x)$ and $b(x)$ as well as the equation for $\gamma$ were already derived in
[11]. Fig. 9a shows the functions
$u_{e}(x_{0}),u_{h}(x_{0}),u_{e}(x_{0})+u_{h}(x_{0})$ that determine how the
average growth of the avalanche depends on the position of a primary electron,
hole or e-h pair. They are the mirror images of the functions $a(x),b(x)$ from
Eqs. 19, 20 that determine the distribution of the electrons and holes inside
the gain layer [2].

The parameter $\gamma$ defines the exponential growth of the avalanche and is
shown in Fig. 9b. It has the following properties:

$\displaystyle\gamma$ $\displaystyle=$ $\displaystyle
0\qquad\qquad\qquad\quad\,\,\alpha\beta
d^{2}=\frac{\alpha\beta}{(\alpha-\beta)^{2}}\ln^{2}\frac{\beta}{\alpha}\leq 1$
(26) $\displaystyle\gamma$ $\displaystyle=$
$\displaystyle\frac{\alpha+\beta}{2}-\sqrt{\alpha\beta}\qquad\alpha\beta
d^{2}=1$ (27) $\displaystyle\gamma$ $\displaystyle=$
$\displaystyle\frac{\alpha+\beta}{2}\phantom{+\sqrt{\alpha\beta}}\qquad\alpha\beta
d^{2}=\frac{\pi^{2}}{4}\approx 2.47$ (28) $\displaystyle\gamma_{max}$
$\displaystyle=$
$\displaystyle\frac{\alpha+\beta}{2}+\sqrt{\alpha\beta}\qquad\alpha\beta
d^{2}\rightarrow\infty$ (29)

a)  b)

Figure 9: a) The functions
$u_{e}(x_{0}),u_{h}(x_{0}),u_{e}(x_{0})+u_{h}(x_{0})$ from Eqs. 17, 18 that
determine how the average growth of the avalanche depends on the position of a
primary electron, hole or e-h pair. The values are for a $1\,\mu$m gain layer
at $E_{1}=4.5\times 10^{5}$ V/cm. b) $\gamma$ as a function of electric field
in silicon for different values of the gain layer thickness $d$. At the
breakdown limit we have $\gamma=0$. For higher electric fields $\gamma$
quickly approach $\gamma_{\max}$.

The total number of electrons and holes is then given by
$N=\int_{0}^{d}n(x,t)dx$ and evaluates to

$N_{e}(t)=B_{e}\,\left[n_{e}^{0}u_{e}(x_{0})+n_{h}^{0}u_{h}(x_{0})\right]\,e^{\gamma
v^{*}t}\qquad
N_{h}(t)=B_{h}\,\left[n_{e}^{0}u_{e}(x_{0})+n_{h}^{0}u_{h}(x_{0})\right]\,e^{\gamma
v^{*}t}\qquad$ (30)

where we have

$\displaystyle B_{e}$ $\displaystyle=$
$\displaystyle\frac{2kv_{h}\left[e^{a_{1}d}(a_{1}d-k\cot k)+k\csc
k\right]}{(v_{e}+v_{h})(a_{1}^{2}d^{2}+k^{2})(1+\lambda_{1})}$ (31)
$\displaystyle B_{h}$ $\displaystyle=$
$\displaystyle\frac{2kv_{e}\left[e^{a_{1}d}(k^{2}+a_{1}d\lambda_{1}+k(a_{1}d-\lambda_{1})\cot
k)+k(\lambda_{1}-a_{1}d)\csc k\right]}{(v_{e}+v_{h})\beta
d(a_{1}^{2}d^{2}+k^{2})(1+\lambda_{1})}$ (32)

The total induced current becomes

$\displaystyle I(t)$ $\displaystyle=$ $\displaystyle
e_{0}\frac{E_{w}}{V_{w}}\left[v_{e}N_{e}(t)+v_{h}N_{h}(t)\right]$ (33)
$\displaystyle=$ $\displaystyle
e_{0}\frac{E_{w}}{V_{w}}\left(v_{e}B_{e}+v_{h}B_{h}\right)\left[n_{e}^{0}u_{e}(x_{0})+n_{h}^{0}u_{h}(x_{0})\right]e^{\gamma
v^{*}t}$ (34)

where $e_{0}$ is the electron charge. Here we have assumed a constant
weighting field $E_{w}/V_{w}$ in the region in which the charges are moving.
For the single photon detection using the conversion layer we have only a
single electron at $x_{0}=0$ and therefore $n_{e}^{0}=1,n_{h}^{0}=0$ and the
expression is

$I(t)=e_{0}\frac{E_{w}}{V_{w}}\left[v_{e}B_{e}+v_{h}B_{h}\right]\sin
k\,e^{\gamma v^{*}t}$ (35)

Assuming a velocity $v_{e}\approx v_{h}\approx v_{sat}\approx 0.1\,\mu$m/ps
and a weighting field of $E_{w}/V_{w}=1/d=1/(1\,\mu$m), the current
corresponding to $10^{5}$ charges at $t=44$ ps in Fig. 8 is 1.6 mA.

##  4 Efficiency

In this section we calculate the probabilities $P_{e}(x)$ and $P_{h}(x)$ for a
single electron or a single hole placed at position $x$ in the gain layer to
cause breakdown. We follow [9] to establish the equations for these
quantities. We start by considering a single electron at position $x-dx$ which
moves in positive $x$-direction in the applied electric field. The probability
for it to create a diverging avalanche is $P_{e}(x-dx)$. Between $x-dx$ and
$x$ two things can happen. 1) With a probability of $(1-\alpha dx)$ there is
no multiplication of the electron and then the electron at position $x$
creates a diverging avalanche or 2) the electron is multiplying over the
distance $dx$ and at least one of the two electrons and the hole create
breakdown. This can be written as

$P_{e}(x-dx)=(1-\alpha dx)P_{e}(x)+\alpha
dx\left[1-(1-P_{e}(x))^{2}(1-P_{h}(x)\right]$ (36)

Writing the corresponding equation for $P_{h}(x)$ and expanding for small $dx$
gives

$\displaystyle\frac{d\,P_{e}(x)}{d\,x}$ $\displaystyle=$
$\displaystyle-\alpha(x)[1-P_{e}(x)]\left[P_{e}(x)+P_{h}(x)-P_{e}(x)P_{h}(x)\right]$
(37) $\displaystyle\frac{d\,P_{h}(x)}{d\,x}$ $\displaystyle=$
$\displaystyle\phantom{+}\beta(x)[1-P_{h}(x)]\left[P_{e}(x)+P_{h}(x)-P_{e}(x)P_{h}(x)\right]$

Provided $\alpha(x)$ and $\beta(x)$ are known, these equations can be
integrated with the boundary conditions $P_{e}(d)=0$ and $P_{h}(0)=0$.
Following [10] we define $P(x)=P_{e}(x)+P_{h}(x)-P_{e}(x)P_{h}(x)$ and by
differentiating this expression and using Eqs. 37 we have

$\frac{d\,P(x)}{d\,x}=-(\alpha-\beta)P(x)[1-P(x)]$ (38)

We use the boundary condition $P(0)=P_{e}(0)=p_{0}$, with $p_{0}$ still to be
determined, giving the solution

$P(x)=\frac{p_{0}}{p_{0}+(1-p_{0})\exp\left[\int_{0}^{x}(\alpha(x^{\prime})-\beta(x^{\prime}))dx^{\prime}\right]}$
(39)

Since $P(x)$ can be written as $P(x)=1-[(1-P_{e}(x))(1-P_{h}(x))]$ we see that
$P(x)$ refers to the breakdown efficiency of a single e-h pair, which we will
use later for the efficiency calculation for MIPs. Knowing $P(x)$ and using
$P_{h}(d)=0$ we can integrate Eqs. 37 and have

$P_{e}(x)=1-\exp\left(-\int_{x}^{d}\alpha(x^{\prime})P(x^{\prime})dx^{\prime}\right)\qquad
P_{h}(x)=1-\exp\left(-\int_{0}^{x}\beta(x^{\prime})P(x^{\prime})dx^{\prime}\right)$
(40)

And finally $P_{e}(0)=p_{0}$ gives the equation that allows us to determine
$p_{0}$

$\displaystyle
p_{0}=1-\exp\left(-\int_{0}^{d}\frac{p_{0}\,\alpha(x^{\prime})}{p_{0}+(1-p_{0})\exp\left[\int_{0}^{x^{\prime}}(\alpha(x^{\prime\prime})-\beta(x^{\prime\prime}))dx^{\prime\prime}\right]}dx^{\prime}\right)$
(41)

In general this equation can only be evaluated numerically. This equation also
reveals again the breakdown condition. Close to the threshold of breakdown the
value of $p_{0}$ will be small. For small values of $p_{0}$ the expression
$1-\exp\left[-p_{0}/(p_{0}+(1-p_{0})e^{v}))\right]$ is approximated by
$p_{0}e^{-v}+O(p_{0}^{2})$, so the above relation turns into the breakdown
condition of Eq. 14.
For constant $\alpha$ and $\beta$ the above expressions evaluate to

$\displaystyle P(x)$ $\displaystyle=$
$\displaystyle\frac{p_{0}}{p_{0}+[1-p_{0}]e^{(\alpha-\beta)x}}$ (42)
$\displaystyle P_{e}(x)$ $\displaystyle=$ $\displaystyle
1-e^{-\alpha(d-x)}\left[\frac{(1-p_{0})e^{(\alpha-\beta)d}+p_{0}}{(1-p_{0})e^{(\alpha-\beta)x}+p_{0}}\right]^{\frac{\alpha}{\alpha-\beta}}$
(43) $\displaystyle P_{h}(x)$ $\displaystyle=$ $\displaystyle 1-e^{-\beta
x}\left[(1-p_{0})e^{(\alpha-\beta)x}+p_{0}\right]^{\frac{\beta}{\alpha-\beta}}$
(44)

Eq. 41 that determines $p_{0}$ reads as

$e^{-(\alpha-\beta)d}=\frac{1}{p_{0}}\left[(1-p_{0})^{1-\frac{\beta}{\alpha}}-(1-p_{0})\right]$
(45)

Fig. 10 shows a few examples.

a)  b)

Figure 10: a) Breakdown probability (efficiency) for a single electron and a
single hole deposited at position $x$ inside a gain layer of $d=1\,\mu$m for
two values of the electric field. b) Breakdown probability $p_{0}$
(efficiency) for a single electron placed at $x=0$ for different values of the
gain layer thickness $d$.

##  5 Time resolution

The statistics of electron-hole avalanches and the resulting contribution to
the time resolution are discussed in detail in [2], where the problem is
treated by using the theory of continuous-time Markov processes.

###  5.1 Avalanches in absence of boundaries

For the case of an e-h avalanche in absence of any boundaries and for a
constant electric field, an alternative approach can be used to derive the
avalanche fluctuations [12] that does not require the formalism developed in
[2] and which is closely related to the arguments resulting in Eq. 37. We
define $p_{e}(n,m,\Delta)$ to be the probability to find $n$ electrons and $m$
holes at time $t_{0}+\Delta$ for an avalanche starting with a single electron
at $t_{0}$. For an avalanche starting at $t=0$, there are two ways to reach
this state at the later time $t+dt$. First, the initial electron does not
multiply in the first small time interval $[0,dt]$ (with probability $1-\alpha
v_{e}dt$), but then produces $n$ electrons and $m$ holes during the subsequent
time interval $[dt,t+dt]$. This happens with a probability $p_{e}(n,m,t)$.
Second, the electron already multiplies in the interval $[0,dt]$ (with
probability $\alpha v_{e}dt$) and the resulting two electrons and one hole
multiply into $n$ electrons and $m$ holes during $[dt,t+dt]$. This is written
as

$\displaystyle p_{e}(n,m,t+dt)$ $\displaystyle=$ $\displaystyle(1-\alpha
v_{e}dt)p_{e}(n,m,t)$ (46) $\displaystyle+$ $\displaystyle\alpha
v_{e}dt\sum_{i=1}^{n}\sum_{j=1}^{i}\sum_{r=1}^{m}\sum_{s=1}^{r}p_{e}(n-i-j,m-r-s,t)p_{e}(i,r,t)p_{h}(j,s,t)$
(47)

where $p_{h}(n,m,\Delta)$ is the probability that an avalanche starting with a
single hole at time $t_{0}$ produces $n$ electrons and $m$ holes at
$t_{0}+\Delta$. Writing the corresponding equation for $p_{h}(n,m,t)$ and
expanding for small $dt$ results in the equations defining $p_{e}(n,m,t)$ and
$p_{h}(n,m,t)$. The equations have a structure similar to Eq. 37. Their
solution is given in Appendix C and equal to the one derived in [2]. Having
$n_{e}^{0}$ electrons and $n_{h}^{0}$ holes at $t=0$, the probability to have
$n$ additionally created e-h pairs at time $t$ is

$p(n,t)=\frac{\Gamma(A+n)}{\Gamma(A)\Gamma(1+n)}\left(\frac{1}{\nu(t)}\right)^{A}\left(1-\frac{1}{\nu(t)}\right)^{n}\qquad\sum_{n=0}^{\infty}p(n,t)=1$
(48) $A=\frac{n_{e}^{0}\alpha v_{e}+n_{h}^{0}\beta v_{h}}{\alpha v_{e}+\beta
v_{h}}\qquad\nu(t)=e^{\lambda_{t}t}\qquad\lambda_{t}=\alpha v_{e}+\beta v_{h}$
(49)

The average number of e-h pairs and the variance are

$\overline{n}(t)=[\nu(t)-1]A\qquad\sigma_{n}^{2}(t)=\nu(t)[\nu(t)-1]A$ (50)

In the continuous approximation for $n$, which is valid for large $n$ (and
therefore for large avalanches i.e. late times) we have

$p(n,t)=\frac{n^{A-1}}{\Gamma(A)}\,\left(\frac{A}{\overline{n}(t)}\right)^{A}e^{-nA/\overline{n}(t)}\qquad\int_{0}^{\infty}p(n,t)dn=1$
(51)

with

$\overline{n}(t)=Ae^{\lambda_{t}t}\qquad\sigma_{n}(t)=\frac{1}{\sqrt{A}}\overline{n}(t)$
(52)

a)  b)

Figure 11: a) A threshold is applied to the signal. The fluctuations of the
avalanche size result in a fluctuation of the threshold crossing time, which
defines the time resolution. b) Time response function $\rho(t)$ for different
values of the parameter $A$.

In order to measure the ’signal time’ we apply a threshold to the signal,
which is proportional to the total number of charge carriers. Fig. 11a shows
how the avalanche fluctuations lead to fluctuations of the threshold crossing
time, which determine the time resolution. The probability that the signal
crosses the threshold of $n$ e-h pairs between time $t$ and $t+dt$, the so-
called time response function, is given by

$\rho(n,t)dt=\lambda_{t}\,\frac{\Gamma(1+A+n)}{\Gamma(A)\Gamma(1+n)}\left(\frac{1}{\nu(t)}\right)^{A}\left(1-\frac{1}{\nu(t)}\right)^{n}dt\qquad\int_{0}^{\infty}\rho(n,t)dt=1$
(53)

The average threshold crossing time and its variance are given by

$\overline{t}=\frac{1}{\lambda_{t}}\left[\psi_{0}(n+A+1)-\psi_{0}(A)\right]\qquad\sigma=\frac{1}{\lambda_{t}}\sqrt{\psi_{1}(A)-\psi_{1}(n+A+1)}$
(54)

where $\psi_{0}(x)=d\ln\Gamma(z)/dz$ is the digamma function and
$\psi_{1}(x)=d^{2}\ln\Gamma(z)/dz^{2}$ is the trigamma function. For large
numbers of $n$ the above time response function approximates to

$\rho(n,t)=\frac{\lambda_{t}}{\Gamma(A)}\exp\left[A\ln
n-A\lambda_{t}t-ne^{-\lambda_{t}t}\right]$ (55)

and the average threshold crossing time and the time resolution approximate to

$\overline{t}=\frac{1}{\lambda_{t}}\left[\log
n-\psi_{0}(A)\right]\qquad\sigma=\frac{1}{\lambda_{t}}\sqrt{\psi_{1}(A)}$ (56)

so the time resolution becomes independent of the threshold. This can be
understood when looking at Fig. 11a and it is a well-established fact for
detectors like Resistive Plate Chambers [13], where avalanche fluctuations
dominate the signal characteristics [14, 15]: scaling the threshold by a
constant $c_{1}$ will just shift the time response function by $\Delta t=(\ln
c_{1})/\lambda_{t}$ without altering its shape. If we are not interested in
the absolute time of the threshold crossing but just the time variations we
can arbitrarily set $n=1$ and have the time response function

$\rho(t)=\frac{\lambda_{t}}{\Gamma(A)}\exp\left[-A\lambda_{t}t-e^{-\lambda_{t}t}\right]\qquad\int_{-\infty}^{\infty}\rho(t)dt=1$
(57)

An example of the time response function for different values of the parameter
$A$ is shown in Fig. 11b. The function $\sqrt{\psi_{1}(x)}$ approaches
$\sqrt{1/x+1/x^{2}}$ for small and large values $x$ and it is within 10% in
the entire range of $x>0$, as can be seen in Fig. 12a. The values of
$\sqrt{\psi_{1}(A)}$ for a primary electron, a primary hole and a primary e-h
pair are shown in Fig. 12b.
We see that when starting with a primary hole, the time resolution is
significantly worse when compared to a primary electron or a primary e-h pair.
An electron travels on average a distance of $1/\alpha$ and a hole travels on
average a distance of $1/\beta$ before creating an additional e-h pair. Since
$\alpha>\beta$ in silicon the holes will travel significantly longer than the
electrons before triggering the breakdown and therefore the time fluctuations
are larger. The geometries of SPADs and SiPMs produced from silicon are
therefore built such that the electrons trigger the avalanche in the gain
region. Charge carrier mobilities are quite different in other semiconductors
and the arrangement of doping layers can therefore differ.

a)  b)

Figure 12: a) The trigamma function $\psi_{1}(x)$. b) Values of
$\sqrt{\psi_{1}(A)}$ for different initial conditions. Starting with a single
e-h pair we have $A=1$ and $\sqrt{\psi_{1}(1)}=\pi/\sqrt{6}\approx 1.28$.

###  5.2 Avalanches in a gain layer of finite thickness

Our central result in [2] is the conclusion that the finite thickness of the
gain layer will to first order not affect the avalanche fluctuations but will
only affect the average growth of the avalanche. This approximation works best
if the primary charge carrier has the larger impact ionization coefficient,
i.e. if the primary charge in the gain layer is either an electron or an e-h
pair for silicon. Then, the probability to have a total number of $n$
electrons and $m$ holes at time $t$ in the gain layer, starting with
$n_{e}^{0}$ electrons and $n_{h}^{0}$ holes at $x=x_{0}$ at time $t=0$, is

$p(n,t)\approx[1-\varepsilon(x_{0})]\delta(n)+\varepsilon(x_{0})\frac{n^{A-1}}{\Gamma(A)}\,\left(\frac{A}{N_{e}(t)}\right)^{A}e^{-nA/N_{e}(t)}\qquad
m=\frac{B_{h}}{B_{e}}n\qquad A=\frac{n_{e}^{0}\alpha v_{e}+n_{h}^{0}\beta
v_{h}}{\alpha v_{e}+\beta v_{h}}\qquad$ (58)

with $N_{e}(t)$ from Eq. 30 and $B_{e},B_{h}$ are from Eqs. 31, 32. The number
of electrons $n$ and the number of holes $m$ are considered continuous
variables and taken to be fully correlated in this approximation (as shown in
[2], this is strictly true only at late times). The efficiency
$\varepsilon(x_{0})$ that the $n_{e}^{0}$ electrons and $n_{h}^{0}$ holes at
$x{=}x_{0}$ trigger a diverging avalanche is

$\varepsilon(x_{0})=1-[1-P_{e}(x_{0})]^{n_{e}^{0}}[1-P_{h}(x_{0})]^{n_{h}^{0}}$
(59)

with $P_{e}(x_{0})$ and $P_{h}(x_{0})$ from Eqs. 43, 44. The corresponding
time response function in reference to Eq. 57 is then

$\rho(t)\approx\frac{\gamma v^{*}\,h^{A}}{\Gamma(A)}\exp\left[-A\gamma
v^{*}t-he^{-\gamma v^{*}t}\right]$ (60)

with

$h=\frac{A\,\varepsilon(x_{0})}{B_{e}\left[n_{e}^{0}u_{e}(x_{0})+n_{h}^{0}u_{h}(x_{0})\right]}$
(61)

Here we have divided $N_{e}(t)$ by the efficiency in order to account for the
fact that only diverging avalanches are crossing the threshold. The
corresponding time resolution due to the avalanche fluctuations, when keeping
$n_{e}^{0},n_{h}^{0},x_{0}$ constant, is

$\sigma=\frac{\sqrt{\psi_{1}(A)}}{\gamma v^{*}}$ (62)

where $\gamma$ and $v^{*}$ are from Eq. 21 and Eq. 23. The value of $1/\gamma
v^{*}$ for different values of the thickness $d$ of the gain layer is shown in
Fig. 13a.

a)  b)

Figure 13: a) The factor $1/\gamma v^{*}$ for different values of the gain
layer thickness. b) Contributions to the time resolution for a photon
interacting in the gain layer. The value of $\pi/\sqrt{6}$ indicated by the
horizontal line is due to avalanche fluctuations and the other lines are due
to fluctuations of the position of the photon conversion assuming a uniform
distribution. The two components have to be added in squares, so we see that
for thick gain layers and high fields the position dependence starts to
dominate over the avalanche fluctuations.

For a photon interacting in the conversion layer (Fig. 2b), the electron
arriving at the gain layer will therefore start an avalanche from $x_{0}=0$
and the time resolution contribution of the gain layer will be

$\sigma\approx\frac{\sqrt{\psi_{1}(A)}}{\gamma v^{*}}\qquad A=\frac{\alpha
v_{e}}{\alpha v_{e}+\beta v_{h}}$ (63)

As seen in Fig. 12 we have $\sqrt{\psi_{1}(A)}\approx 1.5$ in the entire
electric field range, so e.g. for a $1\,\mu$m SPAD at $4\times 10^{5}\,$V/cm
the contribution to the time resolution will be around 6 ps. For higher fields
a time resolution of 1 ps should theoretically be achievable.
In case the photon interacts inside the gain layer and produces an e-h pair,
the varying position of the photon interaction will also contribute to the
time resolution. Having one e-h pair at position $x_{0}$ as shown in Fig. 2a
we have $n_{e}^{0}=n_{h}^{0}=1$ and therefore $A=1$,
$\varepsilon(x_{0})=P(x_{0})$ from Eq. 42 and the time response function is

$\rho(t,x_{0})=\gamma v^{*}\,h(x_{0})\exp\left[-\gamma
v^{*}t-h(x_{0})e^{-\gamma v^{*}t}\right]\qquad
h(x_{0})=\frac{P(x_{0})}{B_{e}\left[u_{e}(x_{0})+u_{h}(x_{0})\right]}\qquad$
(64)

The probability of conversion at position $x_{0}$ is given by

$p_{1}(x_{0})dx_{0}=\frac{1}{1-e^{-d/l_{a}}}\frac{1}{l_{a}}e^{-x_{0}/l_{a}}dx_{0}$
(65)

and the time response function including the fluctuation of the conversion
point is therefore defined by

$\overline{\rho}(t)=\int_{0}^{d}p_{1}(x_{0})\rho(t,x_{0})dx_{0}$ (66)

The time resolution $\sigma^{2}=\int t^{2}\overline{\rho}(t)dt-(\int
t\overline{\rho}(t)dt)^{2}$ has then two components, a contribution
$\sigma_{av}$ from avalanche fluctuations and a contribution $\sigma_{pos}$
from the varying position of the primary e-h pair in the gain layer

$\displaystyle(\gamma v^{*})^{2}\sigma_{av}^{2}$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}}{6}$ (67) $\displaystyle(\gamma
v^{*})^{2}\sigma_{pos}^{2}$ $\displaystyle=$
$\displaystyle\int_{0}^{d}p_{1}(x_{0})\left[\ln
h(x_{0})\right]^{2}dx_{0}-\left(\int_{0}^{d}p_{1}(x_{0})\ln
h(x_{0})dx_{0}\right)^{2}$ (68)

In the case of a very small photon absorption length i.e. $l_{a}\ll d$ the
primary e-h pair will always be created at the very edge of the sensor which
is equal to the situation of $x_{0}=0$ and the time resolution is given by
$\sigma_{av}$. In the other extreme of $l_{a}\gg d$ there will be a uniform
distribution for the position of the photon interaction in the gain layer and
the contribution to the time resolution is given in Fig. 13b. Whether the
avalanche fluctuations or the position fluctuations dominate depends on the
gain layer thickness and the electric field. We can conclude that the time
resolution is well approximated by $\sigma_{1}\approx c_{0}/\gamma v^{*}$,
where the main dependence is given by the variation of $\gamma$ with sensor
thickness and electric field (Fig. 13a) while $c_{0}\approx 1{-}3$ and $v^{*}$
is saturated at $\approx 0.1\,\mu$m/ps.

The term $\sigma_{pos}$ can also be derived directly from the average
avalanche growth. Neglecting avalanche fluctuations the primary e-h pair will
simply trigger an average avalanche according to Eq. 30

$N(t)=B_{e}\frac{u_{e}(x_{0})+u_{h}(x_{0})}{P(x_{0})}e^{\gamma v^{*}t}$ (69)

where we have divided by the efficiency to account for the avalanches that do
not cross the threshold. Applying a threshold $N_{thr}$ to this signal gives a
threshold crossing time of

$t(x_{0})=\frac{1}{\gamma v^{*}}\left[\ln N_{thr}-\ln h(x_{0})\right]$ (70)

and the variance $\overline{t^{2}}-\overline{t}^{2}$ is the one from Eq. 68.

##  6 Charged particle detection with SPADs

Figure 14: A charged particle is leaving clusters of e-h pairs in the gain
layer.

SPADs can also be used for detection of charged particles. Charged particles
interacting with silicon produce clusters of e-h pairs along their track, with
an average distance of $\lambda\approx 0.21\,\mu$m for MIPs. A SPAD sensor of
1 $\mu$m or 2 $\mu$m thickness will therefore be highly efficient to charged
particles and there is no need for a conversion layer. The probability
$p_{clu}(n)$ for having $n>0$ e-h pairs in a single cluster is approximately
given by a $1/n^{2}$ distribution, but there are significant deviations at
small numbers of $n$ in silicon. Fig. 15a shows the cluster size distribution
for silicon as calculated with HEED [16].

a)  b)

Figure 15: a) Cluster size distribution i.e. the probability for a single
cluster to contain $n$ electron-hole pairs as calculated with HEED [16]. b)
Probability $p(n,d)$ to find $n$ e-h pairs inside a silicon layer of thickness
$d$. The solid lines refer to an $1/n^{2}$ distribution in both plots.

Assuming a sensor thickness $d$ and an average distance between clusters of
$\lambda$, the average number of clusters in the sensor is $n_{0}=d/\lambda$.
Assuming the cluster size distribution $p_{clu}(n)$, the probability $p(n,d)$
to find $n$ e-h pairs in the sensor can then be calculated by using the
Z-transform as [15]

$P_{clu}(z)=\sum_{n=1}^{\infty}\frac{p_{clu}(n)}{z^{n}}\qquad
G(z)=\frac{e^{d/\lambda P_{clu}(z)}-1}{e^{d/\lambda}-1}\qquad
p(n,d)=\frac{1}{n!}\left[\frac{d^{n}G(1/z)}{dz^{n}}\right]_{z=0}$ (71)

It is shown in Fig. 15b.

###  6.1 Efficiency

First we want to calculate the efficiency for a charged particle to cause
breakdown in the gain layer. Since $P(x)$ from Eq. 39 is the probability for a
single e-h pair at position $x$ to trigger a diverging avalanche, the
probability $Q(x)$ that this single e-h pair does not cause breakdown is

$Q(x)=1-P(x)=\frac{1}{1+\frac{p_{0}}{1-p_{0}}\exp\left[-\int_{0}^{x}(\alpha(x^{\prime})-\beta(x^{\prime}))dx^{\prime}\right]}$
(72)

The probability $dq$ that there is no breakdown caused by the charged particle
traversing the interval $[x,x+\Delta x]$ is given by a) the probability
$1-\Delta x/\lambda$ that there is no interaction in $\Delta x$ and b) the
probability that there is an interaction but it does not lead to breakdown,
which we can write as

$dq=\left(1-\frac{\Delta x}{\lambda}\right)+\frac{\Delta
x}{\lambda}\sum_{n=1}^{\infty}p_{clu}(n)Q^{n}(x)\\\ $ (73)

From this, the probability $q$ that the charged particle does not create a
diverging avalanche in any of the slices $\Delta x$ is derived in Appendix D
and the efficiency $p=1-q$ is given by

$p=1-\exp\left[-\frac{1}{\lambda}\left(d-\sum_{n=1}^{\infty}p_{clu}(n)\int_{0}^{d}Q(x)^{n}dx\right)\right]$
(74)

If $Q(x)=0$ i.e. if an e-h pair deposited at $x$ will definitely cause
breakdown, we have $p=1-e^{-d/\lambda}$, which is the correct probability that
there is at least one interaction within $d$. The evaluation for constant
$\alpha,\beta$ is given in Appendix E and shown in Fig. 16 together with the
efficiency for a single electron starting at $x_{0}=0$ in the gain layer.
Above the breakdown field the efficiency rises steeply to the maximum level
$1-e^{-d/\lambda}$. For a gain layer thickness of $d>1\,\mu$m the efficiency
is larger that 99 % above the breakdown limit.

Figure 16: Efficiency for a MIP to provoke breakdown. Beyond the breakdown
field, the efficiency rises sharply from 0 to the level of $1-e^{-d/\lambda}$.
The thin lines correspond to the efficiency for a single electron at
$x_{0}=0$.

###  6.2 Time resolution

We first consider the simpler case where the MIP produces a single cluster of
$m$ e-h pairs at position $x_{0}$. Then, we have $A=m$,
$n_{e}^{0}=n_{h}^{0}=m$ and the time response function is

$\rho(t,m,x_{0})\approx\frac{\gamma v^{*}}{(m-1)!}\
h(x_{0})^{m}\exp\left[-m\gamma v^{*}t-h(x_{0})e^{-\gamma v^{*}t}\right]$ (75)
$h(x_{0})=\frac{1-[1-P(x_{0})]^{m}}{B_{e}[u_{e}(x_{0})+u_{h}(x_{0})]}\approx\frac{1}{B_{e}[u_{e}(x_{0})+u_{h}(x_{0})]}$
(76)

In case the cluster size varies according to $p_{clu}(m)$ and the probability
to have the cluster at position $x_{0}$ varies according to $p_{1}(x_{0})$,
the time response function becomes

$\overline{\rho}(t)=\sum_{m=1}^{\infty}\int_{0}^{d}p_{clu}(m)p_{1}(x_{0})\rho(t,m,x_{0})dx_{0}$
(77)

and the related time resolution has three contributions

$\displaystyle(\gamma v^{*})^{2}\sigma^{2}$ $\displaystyle=$
$\displaystyle\sum_{m=1}^{\infty}p_{clu}(m)\psi_{1}(m)$ (78) $\displaystyle+$
$\displaystyle\sum_{m=1}^{\infty}p_{clu}(m)\psi_{0}(m)^{2}-\left(\sum_{m=1}^{\infty}p_{clu}(m)\psi_{0}(m)\right)^{2}$
(79) $\displaystyle+$ $\displaystyle\int_{0}^{d}p_{1}(x_{0})[\ln
h(x_{0})]^{2}dx_{0}-\left(\int_{0}^{d}p_{1}(x_{0})\ln
h(x_{0})dx_{0}\right)^{2}$ (80)

The first term represents the average of the avalanche fluctuations, where
$\psi_{1}(m)$ is decreasing from $\psi_{1}(1)=\pi^{2}/6$ to zero with $\approx
1/m$ dependence. For the cluster size distribution in silicon from Fig. 15a
this first term evaluates to $\approx 0.36$.
The second term represents the fact that an avalanche starting with $m$ e-h
pairs will on average grow as $me^{\gamma v^{*}t}$. The term evaluates to
$\approx 1.39$ for the cluster size distribution in silicon, significantly
larger than the contribution from the avalanche fluctuations.
The third term represents the dependence on the position of the primary
cluster and for a uniform probability this term evaluates to the values
already shown in Fig. 13b.
In conclusion we therefore observe that assuming a single e-h cluster at a
random position in the gain layer, the avalanche fluctuations are negligible
and only the average growth of the avalanche as well as the position
dependence play a role. In a regime where the contribution from the position
dependence is negligible ($<1$) only the fluctuation of the total charge is
important, which results in a universal dependence of the time resolution on
the thickness of the gain layer

$(\gamma^{2}v^{*})^{2}\sigma^{2}=\sum_{m=1}^{\infty}p(m,d)\psi_{0}(m)^{2}-\left(\sum_{m=1}^{\infty}p(m,d)\psi_{0}(m)\right)^{2}$
(81)

Here, $p(m,d)$ is the probability that the passing MIP produces $m$ e-h pairs
in the gain layer of thickness $d$. The resulting time resolution is shown in
Fig. 17a. Since $\gamma v^{*}\sigma$ is close to unity for typical dimensions
of the gain layer, the time resolution of a SPAD for a MIP is essentially
defined only by $\gamma$ and $v^{*}$ and the values are the ones given in Fig.
13a. Since the cluster size distribution $p_{clu}(m)$ has a long tail towards
large values of $m$, the same is true for the time response function. The
standard deviation of the threshold crossing time is then generally not
identical to the parameter $\sigma$ extracted from a Gaussian fit to the
distribution of the threshold crossing time. Both measures are compared in
Figure 17a.

a)  b)

Figure 17:  a) Time resolution for a MIP when the fluctuations of the cluster
position can be neglected, for different values of the gain layer thickness
$d$. The solid lines refer to $p_{clu}(n)$ for silicon from Fig. 15a, while
the dashed lines assume $p_{clu}(n)\sim 1/n^{2}$. The time resolution is
quantified by the standard deviation of the threshold crossing time (“r.m.s.”)
as well as by the parameter $\sigma$ extracted from a Gaussian fit. b) Time
resolution for a MIP for a gain layer of $0.5,1,2\,\mu$m thickness, taking the
fluctuation of the cluster positions into account. The dashed lines refer to
the numbers from a).

Next, we consider the general case where the MIP produces a variable number of
clusters, all of which fluctuate in size. Since avalanche fluctuations are
negligible we perform this calculation by following the discussion around Eq.
70. We divide the gain layer into $N+1$ slices of thickness $\Delta x=d/(N+1)$
and assume that a charged particle leaves $m_{n}$ primary e-h pairs in the
$n^{th}$ slice. The efficiency for this case is very close to unity, and the
additional small dependence on the cluster distribution is neglected here.
Then, the total average number of charges in the avalanche becomes

$N_{tot}(t)=\sum_{n=0}^{N}m_{n}h\left(n\Delta x\right)e^{\gamma v^{*}t}$ (82)

with $h(x)$ taken from Eq. 76. Applying a threshold $N_{thr}$ to this signal
and shifting it by a constant offset of $\ln N_{thr}$ we find a threshold
crossing time of

$t(m_{0},m_{1},...,m_{N})=-\frac{1}{\gamma
v^{*}}\ln\left[\sum_{n=0}^{N}m_{n}h\left(n\Delta x\right)\right]$ (83)

The probability $p(m,\Delta x)$ to find $m$ e-h pairs in a slice $\Delta x$ is
given by

$p(m,\Delta x)=\left(1-\frac{\Delta x}{\lambda}\right)\delta_{m0}+\frac{\Delta
x}{\lambda}p_{clu}(m)$ (84)

so the average threshold crossing time and the second moment are

$\displaystyle\overline{t}$ $\displaystyle=$
$\displaystyle-\lim_{N\rightarrow\infty}\frac{1}{\gamma
v^{*}}\sum_{m_{0}=0}^{\infty}\sum_{m_{1}=0}^{\infty}...\sum_{m_{N}=0}^{\infty}p(m_{0},\Delta
x)p(m_{1},\Delta x)...p(m_{N},\Delta
x)\ln\left[\sum_{n=0}^{N}m_{n}h\left(n\Delta x\right)\right]$ (85)
$\displaystyle\overline{t^{2}}$ $\displaystyle=$
$\displaystyle\lim_{N\rightarrow\infty}\frac{1}{(\gamma
v^{*})^{2}}\sum_{m_{0}=0}^{\infty}\sum_{m_{1}=0}^{\infty}...\sum_{m_{N}=0}^{\infty}p(m_{0},\Delta
x)p(m_{1},\Delta x)...p(m_{N},\Delta
x)\ln^{2}\left[\sum_{n=0}^{N}m_{n}h\left(n\Delta x\right)\right]$ (86)

These relations are evaluated in Appendix F, giving a contribution to the time
resolution of

$(\gamma
v^{*})^{2}\sigma^{2}=\left[\int_{0}^{\infty}w_{0}(y)\ln^{2}ydy-\left(\int_{0}^{\infty}w_{0}(y)\ln
ydy\right)^{2}\right]$ (87)

with

$W_{0}(s)=\frac{\exp\left[\frac{d}{\lambda}\,\frac{1}{d}\int_{0}^{d}P_{clu}\left(s\frac{h(x)}{h(0)}\right)dx\right]-1}{e^{d/\lambda}-1}\qquad\qquad
w_{0}(y)=L^{-1}[W_{0}(s)]$ (88)

Here, $P_{clu}(s)$ is the Laplace transform of the cluster size distribution
and the operator $L^{-1}$ denotes the inverse Laplace transform. The
evaluation is shown in Fig. 17b. The dashed lines show the time resolution
from Fig. 17a where the contribution from position fluctuations is neglected.
We see that for fields around the breakdown limit the effect from the
positions variations is small and it increases with increasing field. Overall,
the time resolution for MIPs stays within $\sigma_{t}=(0.8{-}1.5)/\gamma
v^{*}$ for the parameters investigated.

##  7 Realistic field configuration

In this section we finally discuss a realistic field configuration of a SPAD
and we apply the insights from all previous sections to assess its
performance. Fig. 18a shows an example for the electric field in a SPAD
created by a highly doped p-n junction. The specific functional form is
defined in Appendix A, Eq. 95. With the impact ionization and drift velocity
parameters from Appendix A we obtain $\alpha(x)$, $\beta(x)$, $v_{e}(x)$ and
$v_{h}(x)$ as explicit functions throughout the sensor. The impact ionization
coefficients $\alpha(x)$ and $\beta(x)$ are shown in Fig. 18b. For the
purposes of our discussion here, we define $x_{1}=0.4\,\mu$m and
$x_{2}=1.9\,\mu$m as the boundaries of the gain layer, which thus has a
thickness of $d=1.5\,\mu$m.

a)  b)

Figure 18: a) Electric field in a realistic SPAD or SiPM. b) Impact ionization
coefficients inside the sensor with parameters from Appendix A.

#### Efficiency

For this geometry, the integral in Eq. 14 evaluates to 1.39 which is larger
than unity and therefore guarantees that breakdown can take place. To find the
efficiency of the sensor we solve Eqs. 37 numerically, using as boundary
conditions $P_{h}(x_{1})=0$ and $P_{e}(x_{2})=0$. The solution is shown in
Fig. 19a together with the corresponding efficiencies obtained from a MC
simulation of the avalanche development. The efficiency for a MIP passing this
sensor can be calculated with Eq. 74 using $Q(x)=1-P_{eh}(x)$. With the
cluster size distribution from Fig. 15 and $\lambda=0.21\,\mu$m, the
efficiency evaluates to $p=1-8.5\times 10^{-4}$. A SPAD of this kind is a
highly efficient detector for a MIP.

a)  b)

Figure 19: a) Breakdown probability for a primary electron ($P_{e}$), primary
hole ($P_{h}$) and primary e-h pair $(P_{eh})$ placed at position $x_{0}$ in
the gain layer. b) Time at which the average total number of charges in the
avalanche (proportional to the average signal) crosses a threshold of $10^{4}$
charges. The avalanche is initiated by a primary electron, a primary hole or a
primary e-h pair placed at position $x_{0}$ in the gain layer. For both plots,
the numerical solutions give rise to the lines, the markers correspond to the
values from MC simulations.

#### Average signal and contribution to time resolution

To study the average growth of the avalanche, we solve Eqs. 11 with
$n_{e}(x_{1})=0$ and $n_{h}(x_{2})=0$ as boundary conditions. This yields the
average charge densities in the gain layer, $n_{e}(x,t)$ and $n_{h}(x,t)$. The
average total charge present in the gain layer can be obtained through a
numerical integration of these densities. This quantity is proportional to the
average signal produced by the avalanche. As shown in Section 3, it grows
exponentially as $e^{St}$. The time constant $S$ can be directly extracted
from the numerical solution and evaluates to $S=0.48$ ps.

In case the position $x_{0}$ of the primary charge fluctuates, it generates a
contribution to the time resolution according to Eq. 80. The magnitude of this
effect can be estimated from Fig. 19b, which shows the time at which the
average signal crosses the applied threshold. If the position $x_{0}$ of the
initial charge is uniformly distributed, the resulting contribution to the
time resolution is 2.7/2.1/2.4 ps for an initial electron, an initial hole,
and an initial e-h pair.

#### Avalanche fluctuations and contribution to time resolution

According to the discussion in Section 6.2, we expect the contribution of the
time resolution from fluctuations in the avalanche development,
$\sigma_{\mathrm{av}}$, to be of the order of $1/S\approx 2$ ps. A more
precise estimate of the time resolution takes into account the primary charge
initiating the avalanche. Following the discussion leading to Eq. 62, we
approximate $\sigma_{\mathrm{av}}\approx\sqrt{\psi_{1}(A)}/S$. This formula
neglects effects due to the finite size of the gain region and was originally
derived for constant impact ionization coefficients, which enter into the
computation of the parameter $A$. For position-dependent electric fields, the
largest values of $\alpha$ and $\beta$ in the gain layer are relevant for the
formation of the avalanche fluctuations and can be used to compute $A$. As the
comparison with results from MC in Fig. 20 shows, this estimates the time
resolution for the case of an initial electron to within at most 20%. As
expected, the approximation becomes worse if the initial charge includes
holes, which have a low impact ionization coefficient. In this case,
corrections due to the finite size of the gain layer become important. These
can also be computed numerically, as shown in [2], but the calculations are
more involved.

Figure 20: Time resolution $\sigma_{\mathrm{av}}$ due to avalanche
fluctuations, for an avalanche initiated by a primary electron, a primary hole
or a primary e-h pair. The primary charge is placed at position $x_{0}$ in the
gain layer. The approximation from Eq. 62 neglects the $x_{0}$-dependence and
gives rise to the horizontal lines. The markers correspond to the values
obtained from MC simulations.

##  8 Conclusions

We have performed a detailed study of the time resolution and efficiency of
SPADs and SiPMs for the detection of photons and charged particles. Our
discussions start from a series of differential equations, which cover the
conversion and the drift of charges in the conversion layer as well as the
formation of the avalanche in the gain layer. For arbitrary electric field
profiles, the equations for the average avalanche development as well as the
breakdown efficiency can be easily solved using numeric solvers. The
calculation of the avalanche fluctuations and their impact on the time
resolution is more involved in this case and a detailed discussion is given in
[2]. We have provided analytic solutions for the case of constant electric
fields.

For the detection of single photons, the contribution of the conversion layer
of thickness $w$ to the time resolution for constant electric field is
$\sigma=w/(v_{e}\sqrt{12})$. This is valid in case the photon absorption
length is much larger than the conversion layer thickness, which simply
corresponds to a uniform distribution of the photon conversion point inside
the layer.

The contribution of the gain layer to the time resolution has the general form

$\sigma=\frac{c_{0}}{\gamma v^{*}}$ (89)

with $c_{0}=0.8{-}2.5$ for silicon and a gain layer thickness of
$0.5{-}2\mu$m. This relation holds for single photon detection and MIP
detection. It also extends to realistic non-uniform electric fields. Both
contributions from avalanche fluctuations as well as the variation of the
photon or MIP conversion point in the gain layer are captured. The constant
$\gamma v^{*}$ determines the average growth of the avalanche according to
$N(t)\propto e^{\gamma v^{*}t}$, with $v^{*}\approx 0.1\mu$m/ps and $\gamma$
saturating at $\gamma_{max}\approx\alpha+\beta$ at high fields. It should be
possible in practice to limit this contribution to the level of a few
picoseconds at high fields.

The efficiency of a SPAD or SiPM for photons has many contributions, including
the photon conversion probability, the geometry and fill factor of the sensor
as well as the breakdown probability in the gain layer. The contribution from
the breakdown probability can be easily calculated by numerically solving the
related equations. For SPADs with a conversion layer or for photons with
absorption length $<1\,\mu$m that are absorbed close to the edge of the gain
layer this efficiency quickly approaches values close to 100% when biasing the
sensor beyond the breakdown field.

SPADs or SiPMs with a gain layer of $1{-}2\,\mu$m thickness should be highly
efficient for MIP detection. A dedicated conversion layer is not necessary.

This report discussed ’one-dimensional’ sensors. For realistic implementations
of SPAD pixel sensors and SiPMs the boundaries of the pixels together with all
the elements for limitation of optical crosstalk make up complex three
dimensional electric fields. To study these sensors, the 3D field map together
with a full MC simulation with programs like Garfield++ [17] has to be used
and our results can serve as benchmarks for these simulations.

##  9 Appendix A

###  9.1 Velocity of electrons and holes in silicon

The velocity of electrons and holes in silicon is shown is parametrized by

$v_{e}(E)=\frac{\mu_{e}\,E}{\left[1+\left(\frac{\mu_{e}\,E}{v^{e}_{sat}}\right)^{\beta_{e}}\right]^{1/\beta_{e}}}\qquad\qquad
v_{h}(E)=\frac{\mu_{h}\,E}{\left[1+\left(\frac{\mu_{h}\,E}{v^{h}_{sat}}\right)^{\beta_{h}}\right]^{1/\beta_{h}}}$
(90)

The parameters from [5] are $\mu_{e}=1417$ cm2/Vs, $\mu_{h}=471$ cm2/Vs,
$\beta_{e}=1.109$, $\beta_{h}=1.213$ and $v^{e}_{sat}=1.07\times 10^{7}$ cm/s
and $v^{h}_{sat}=0.837\times 10^{7}$ cm/s at 300 K .

###  9.2 Impact ionization coefficients for electrons and holes in silicon

The impact ionization parameters $\alpha$ and $\beta$ as reported in [7] are
given by

$\alpha(E)=\alpha_{\infty}e^{-a/E}\qquad\beta(E)=\beta_{\infty}e^{-b/E}$ (91)

with

$\displaystyle\alpha_{\infty}=7.030\times 10^{5}\,cm^{-1}$ $\displaystyle
a=1.231\times 10^{6}\,V/cm\qquad 1.75\times 10^{5}\leq E\leq 6.0\times
10^{5}\,V/cm$ (92) $\displaystyle\beta_{\infty}=1.582\times 10^{6}\,cm^{-1}$
$\displaystyle b=2.036\times 10^{6}\,V/cm\qquad 1.75\times 10^{5}\leq E\leq
4.0\times 10^{5}\,V/cm$ (93) $\displaystyle\beta_{\infty}=6.710\times
10^{5}\,cm^{-1}$ $\displaystyle b=1.693\times 10^{6}\,V/cm\qquad 4.00\times
10^{5}\leq E\leq 6.0\times 10^{5}\,V/cm$ (94)

###  9.3 Functional form of electric field

For a realistic SPAD we assume the following electric field:

$E(x)=E_{0}\,\exp\left[1-(x-\mu)/\sigma-e^{-(x-\mu)/\sigma}\right]\qquad
E_{0}=5\times
10^{5}\,\mbox{V/cm}\quad\mu=1\,\mu\mbox{m}\quad\sigma=0.5\,\mu\mbox{m}$ (95)

##  10 Appendix B

For $S=0$ Eq. 13 reads as

$\displaystyle[v_{e}(x)f(x)]^{\prime}$ $\displaystyle=$
$\displaystyle\alpha(x)v_{e}(x)f(x)+\beta(x)v_{h}(x)g(x)$ (96)
$\displaystyle-[v_{h}(x)g(x)]^{\prime}$ $\displaystyle=$
$\displaystyle\alpha(x)v_{e}(x)f(x)+\beta(x)v_{h}(x)g(x)$ (97)

with boundary conditions $f(0)=0$ and $g(d)=0$. Subtracting the two equations
gives

$[v_{e}(x)f(x)]^{\prime}+[v_{h}(x)g(x)]^{\prime}=0\qquad\rightarrow\qquad
v_{e}(x)f(x)=-v_{h}(x)g(x)+c_{1}$ (98)

Inserting this expression into Eq. 96 we have

$[v_{h}(x)g(x)]^{\prime}-v_{h}(x)g(x)[\alpha(x)-\beta(x)]=-c_{1}\alpha(x)$
(99)

with the general solution

$v_{h}(x)g(x)=\frac{c_{2}-c_{1}\int_{0}^{x}\alpha(x^{\prime})\exp[-\int_{0}^{x^{\prime}}(\alpha(x^{\prime\prime})-\beta(x^{\prime\prime}))dx^{\prime\prime}]dx^{\prime}}{\exp[-\int_{0}^{x}(\alpha(x^{\prime})-\beta(x^{\prime}))dx^{\prime}]}$
(100)

The condition $f(0)=0$ implies $c_{1}=c_{2}$ and $g(d)=0$ then implies

$\int_{0}^{d}\alpha(x)\exp\left[-\int_{0}^{x}(\alpha(x^{\prime})-\beta(x^{\prime}))dx^{\prime}\right]dx=1$
(101)

This is the general breakdown condition which is independent of the electron
and hole velocities. Expressing $v_{h}g$ from Eq. 98 and inserting this
expression into Eq. 97 we have

$v_{h}(x)g(x)=-v_{e}(x)f(x)+c_{1}\qquad[v_{e}(x)f(x)]^{\prime}-v_{e}(x)f(x)[\alpha(x)-\beta(x)]=c_{1}\beta(x)$
(102)

and therefore

$v_{e}(x)f(x)=\frac{c_{1}\int_{0}^{x}\beta(x^{\prime})\exp[-\int_{0}^{x^{\prime}}(\alpha(x^{\prime\prime})-\beta(x^{\prime\prime}))dx^{\prime\prime}]dx^{\prime}+C}{\exp[-\int_{0}^{x}(\alpha(x^{\prime})-\beta(x^{\prime}))dx^{\prime}]}$
(103)

The condition $f(0)=0$ gives $C=0$ and $g(d)=0$ gives

$\int_{0}^{d}\beta(x)\exp\left[-\int_{0}^{x}(\alpha(x^{\prime})-\beta(x^{\prime}))dx^{\prime}\right]dx=\exp\left[-\int_{0}^{d}(\alpha(x^{\prime})-\beta(x^{\prime}))dx^{\prime}\right]$
(104)

which is equal to

$\int_{0}^{d}\beta(x)\exp\left[\int_{x}^{d}(\alpha(x^{\prime})-\beta(x^{\prime}))dx^{\prime}\right]dx=1$
(105)

As shown in [8] Eq. 101 and 105 are identical because

$\int_{0}^{d}(\alpha(x)-\beta(x))\exp\left[-\int_{0}^{x}(\alpha(x^{\prime})-\beta(x^{\prime}))dx^{\prime}\right]dx=1-\exp\left[-\int_{0}^{d}(\alpha(x^{\prime})-\beta(x^{\prime}))dx^{\prime}\right]$
(106)

##  11 Appendix C

The relation that determines the number of electrons $n$ and the number of
holes $m$ at a given time when starting with a single electron is defined in
Eq. 46 as

$\displaystyle p_{e}(n,m,t+dt)$ $\displaystyle=$ $\displaystyle(1-\alpha
v_{e}dt)p_{e}(n,m,t)$ $\displaystyle+$ $\displaystyle\alpha
v_{e}dt\sum_{i=1}^{n}\sum_{j=1}^{i}\sum_{r=1}^{m}\sum_{s=1}^{r}p_{e}(n-i-j,m-r-s,t)p_{e}(i,r,t)p_{h}(j,s,t)$

Establishing the corresponding equation for $p_{h}(m,n,t)$ and expanding for
small $dx$ we have

$\displaystyle\frac{1}{v_{e}}\frac{dp_{e}(n,m,t)}{dt}$ $\displaystyle=$
$\displaystyle-\alpha p_{e}(n,m,t)$ $\displaystyle+$
$\displaystyle\alpha\sum_{i=1}^{n}\sum_{j=1}^{i}\sum_{r=1}^{m}\sum_{s=1}^{r}p_{e}(n-i-j,m-r-s,t)p_{e}(i,r,t)p_{h}(j,s,t)$
$\displaystyle\frac{1}{v_{h}}\frac{dp_{h}(n,m,t)}{dt}$ $\displaystyle=$
$\displaystyle-\beta p_{h}(n,m,t)$ $\displaystyle+$
$\displaystyle\beta\sum_{i=1}^{n}\sum_{j=1}^{i}\sum_{r=1}^{m}\sum_{s=1}^{r}p_{h}(n-i-j,m-r-s,t)p_{h}(i,r,t)p_{e}(j,s,t)$

The Z-transform of these equations is

$\frac{1}{v_{e}}\frac{\partial P_{e}(z_{1},z_{2},t)}{\partial t}=-\alpha
P_{e}(z_{1},z_{2},t)[1-P_{e}(z_{1},z_{2},t)P_{h}(z_{1},z_{2},t)]$ (110)
$\frac{1}{v_{h}}\frac{\partial P_{h}(z_{1},z_{2},t)}{\partial t}=-\beta
P_{h}(z_{1},z_{2},t)[1-P_{h}(z_{1},z_{2},t)P_{e}(z_{1},z_{2},t)]$ (111)

The equations have a structure similar to 37 and we therefore form the
function

$P(z_{1},z_{2},t)=1-P_{e}(z_{1},z_{2},t)P_{h}(z_{1},z_{2},t)$ (112)

Differentiating this equation and using the above relations gives

$\frac{\partial P}{\partial t}=(\alpha v_{e}+\beta v_{h})(1-P)P$ (113)

with the solution

$P(t)=\frac{e^{(\alpha v_{e}+\beta v_{h})t}}{e^{(\alpha v_{e}+\beta
v_{h})t}+c_{1}}$ (114)

The initial conditions that there is one electron at $t=0$ for $p_{e}$ and one
hole for $p_{h}$ reads as $p_{e}(n,m,t=0)=\delta_{n,1}\delta_{m,0}$ and
$p_{h}(n,m,t=0)=\delta_{n,0}\delta_{m,1}$ and we have therefore

$P_{e}(z_{1},z_{2},t=0)=\frac{1}{z_{1}}\qquad
P_{h}(z_{1},z_{2},t=0)=\frac{1}{z_{2}}\qquad\rightarrow\qquad
P(z_{1},z_{2},t=0)=1-\frac{1}{z_{1}z_{2}}$ (115)

and

$P(z_{1},z_{2},t)=\frac{e^{(\alpha v_{e}+\beta
v_{h})t}(z_{1}z_{2}-1)}{1+e^{(\alpha v_{e}+\beta v_{h})t}(z_{1}z_{2}-1)}$
(116)

We can now write Eqs. 110 and 111 as

$\frac{d\ln P_{e}}{dt}=-\alpha v_{e}P\qquad\frac{d\ln P_{h}}{dt}=-\beta
v_{h}P$ (117)

Integrating the equations with the above initial conditions we finally have

$\displaystyle P_{e}(z_{1},z_{2},t)$ $\displaystyle=$
$\displaystyle\frac{1}{z_{1}}\left[\frac{z_{1}z_{2}}{1+e^{(\alpha v_{e}+\beta
v_{h})t}(z_{1}z_{2}-1)}\right]^{\frac{\alpha v_{e}}{\alpha v_{e}+\beta
v_{h}}}$ (118) $\displaystyle P_{h}(z_{1},z_{2},t)$ $\displaystyle=$
$\displaystyle\frac{1}{z_{2}}\left[\frac{z_{1}z_{2}}{1+e^{(\alpha v_{e}+\beta
v_{h})t}(z_{1}z_{2}-1)}\right]^{\frac{\beta v_{h}}{\alpha v_{e}+\beta v_{h}}}$
(119)

In case the avalanche is initiated by $n_{e}^{0}$ electrons and $n_{h}^{0}$
holes at time $t=0$, we are interested in the probability $p_{0}(n,m,t)$ to
find $n$ electrons and $m$ holes at time $t$. In terms of $p_{e}(n,m,t)$ and
$p_{h}(n,m,t)$, it is expressed as an iterated convolution in analogy to the
right-hand sides of Eqs. 11 and 11. In the $z$-domain, $P_{0}(z_{1},z_{2},t)$
reads

$\displaystyle P_{0}(z_{1},z_{2},t)$ $\displaystyle=$ $\displaystyle
P_{e}(z_{1},z_{2},t)^{n_{e}^{0}}P_{h}(z_{1},z_{2},t)^{n_{h}^{0}}$ (120)
$\displaystyle=$
$\displaystyle\frac{1}{z_{1}^{n_{e}^{0}}}\frac{1}{z_{2}^{n_{h}^{0}}}\left[\frac{z_{1}z_{2}}{1+e^{(\alpha
v_{e}+\beta v_{h})t}(z_{1}z_{2}-1)}\right]^{\frac{n_{e}^{0}\alpha
v_{e}+n_{h}^{0}\beta v_{h}}{\alpha v_{e}+\beta v_{h}}}$ (121)

The inverse Z-Transform of this expression also gives access to $p(n,t)$ shown
in Eq. 48, which is defined as the probability to find $n$ e-h pairs that are
created in addition to the initial $n_{e}^{0}$ electrons and $n_{h}^{0}$
holes.

$p(n,t)=\frac{\Gamma(A+n)}{\Gamma(A)\Gamma(1+n)}\left(\frac{1}{\nu(t)}\right)^{A}\left(1-\frac{1}{\nu(t)}\right)^{n}\qquad\sum_{n=0}^{\infty}p(n,t)=1$
(122) $A=\frac{n_{e}^{0}\alpha v_{e}+n_{h}^{0}\beta v_{h}}{\alpha v_{e}+\beta
v_{h}}\qquad\nu(t)=e^{(\alpha v_{e}+\beta v_{h})t}$ (123)

##  12 Appendix D

We divide the sensor into $N+1$ slices of thickness $\Delta x=d/(N+1)$. The
probability that there is no breakdown caused by the particle traversing the
slice $[x,x+\Delta x]$ is given by the probability $1-\Delta x/\lambda$ that
there is no interaction in $\Delta x$ and the probability that there is an
interaction but it does not lead to breakdown.

$\displaystyle dq$ $\displaystyle=$ $\displaystyle\left(1-\frac{\Delta
x}{\lambda}\right)+\frac{\Delta
x}{\lambda}\sum_{n=1}^{\infty}p_{clu}(n)Q^{n}(x)$ (124) $\displaystyle=$
$\displaystyle 1-\frac{\Delta
x}{\lambda}\left(1-\sum_{n=1}^{\infty}p_{clu}(n)Q^{n}(x)\right)$ (125)
$\displaystyle:=$ $\displaystyle 1-\frac{\Delta x}{\lambda}f(x)$ (126)

The probability $q$ that there is no breakdown in any of the slices of $\Delta
x$ throughout the sensor is then

$\displaystyle q$ $\displaystyle=$ $\displaystyle\left[1-\frac{\Delta
x}{\lambda}f(0)\right]\left[1-\frac{\Delta x}{\lambda}f(\Delta
x)\right]\left[1-\frac{\Delta x}{\lambda}f(2\Delta
x)\right]...\left[1-\frac{\Delta x}{\lambda}f(N\Delta x)\right]$ (127)
$\displaystyle\ln q$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{N}\ln\left[1-\frac{d}{N\lambda}f(nd/N)\right]$ (128)
$\displaystyle\approx$ $\displaystyle\sum_{n=0}^{N}-\frac{d}{N\lambda}f(nd/N)$
(129) $\displaystyle\approx$
$\displaystyle-\frac{1}{\lambda}\int_{0}^{\infty}f(x)dx$ (130)
$\displaystyle=$
$\displaystyle-\frac{1}{\lambda}\left(d-\sum_{n=1}^{\infty}p_{clu}(n)\int_{0}^{d}Q(x)^{n}dx\right)$
(131)

and the probability $p=1-q$ that the sensor is efficient is therefore

$p=1-\exp\left[-\frac{1}{\lambda}\left(d-\sum_{n=1}^{\infty}p_{clu}(n)\int_{0}^{d}Q(x)^{n}dx\right)\right]$
(132)

##  13 Appendix E

$Q(x)=1-P(x)=\frac{1}{1+\frac{p_{0}}{1-p_{0}}\exp\left[-(\alpha-\beta)x\right]}$
(133)
$\int_{0}^{d}Q(x)^{n}dx=\frac{1}{\alpha-\beta}\left[H\left(n,1+\frac{p_{0}\exp[-(\alpha-\beta)d]}{1-p_{0}}\right)-H\left(n,1+\frac{p_{0}}{1-p_{0}}\right)\right]$
(134)
$H(n,y)=\int\frac{dy}{y^{n}(1-y)}=\ln\frac{y}{1-y}-\sum_{m=1}^{n-1}\frac{1}{my^{m}}$
(135)

##  14 Appendix F

We assume the probability $p_{clu}(n)$ to be continuous in $n$ (which we can
imagine by expressing it as a sum of delta functions centered at integer
values of $n$) and write the expression in Eq. 84 as

$p(n,\Delta x)=\left(1-\frac{\Delta x}{\lambda}\right)\delta(n)+\frac{\Delta
x}{\lambda}p_{clu}(n)$ (136)

We can then replace the sums in Eq. 85 by integrals and have

$\overline{t}(N)=-\frac{1}{\gamma v^{*}}\int dm_{0}\int dm_{1}...\int
dm_{N}p(m_{0},\Delta x)p(m_{1},\Delta x)...p(m_{N},\Delta
x)\ln\left[\sum_{n=0}^{N}m_{n}h\left(n\Delta x\right)\right]$ (137)

We change variables according to

$m=\frac{1}{h_{0}}\sum_{n=0}^{N}m_{n}h_{n}\quad\rightarrow\quad
m_{0}=m-\frac{1}{h_{0}}\sum_{n=1}^{N}m_{n}h_{n}$ (138)

where we have written $h_{n}=h(n\Delta x)$, which gives

$\displaystyle\overline{t}(N)$ $\displaystyle=$ $\displaystyle-\frac{1}{\gamma
v^{*}}\int dm\left[\int dm_{1}...\int
dm_{N}p\left(m-\frac{1}{h_{0}}\sum_{n=1}^{N}m_{n}h_{n},\Delta
x\right)p(m_{1},\Delta x)...p(m_{N},\Delta x)\right]\ln(h_{0}m)$ (139)
$\displaystyle=$ $\displaystyle-\frac{1}{\gamma v^{*}}\int
w(m)\ln(h_{0}m)\,dm$

with

$w(m)=\int dm_{1}\int dm_{2}...\int
dm_{N}p\left(m-\frac{1}{h_{0}}\sum_{n=1}^{N}m_{n}h_{n},\Delta
x\right)p(m_{1},\Delta x)...p(m_{N},\Delta x)$ (140)

The Laplace transform of this expression is

$W(s)=P(s,\Delta x)P\left(\frac{h_{1}}{h_{0}}s,\Delta
x\right)P\left(\frac{h_{2}}{h_{0}}s,\Delta
x\right)...P\left(\frac{h_{N}}{h_{0}}s,\Delta
x\right)=\exp\left[\sum_{n=0}^{N}\ln P\left(\frac{h_{n}}{h_{0}}s,\Delta
x\right)\right]$ (141)

With $P(s,\Delta x)$ being the Laplace transform of Eq. 136

$P(s,\Delta x)=1+\frac{\Delta x}{\lambda}(P_{clu}(s)-1)$ (142)

we have

$\displaystyle W(s)$ $\displaystyle=$
$\displaystyle\exp\left(\sum_{n=1}^{N}\ln\left[1+\frac{\Delta
x}{\lambda}(P_{clu}\left(s\frac{h_{n}}{h_{0}}\right)-1)\right]\right)$ (143)
$\displaystyle\approx$ $\displaystyle\exp\left(\sum_{n=1}^{N}\frac{\Delta
x}{\lambda}\left[P_{clu}\left(s\frac{h_{n}}{h_{0}}\right)-1\right]\right)$
(144) $\displaystyle=$ $\displaystyle
e^{-d/\lambda}\exp\left[\frac{1}{\lambda}\int_{0}^{d}P_{clu}\left(s\frac{h(x)}{h(0)}\right)dx\right]$
(145)

Normalizing to the probability that there is at least one interaction inside
the gain layer we have

$W_{0}(s)=\frac{\exp\left[\frac{d}{\lambda}\,\frac{1}{d}\int_{0}^{d}P_{clu}\left(s\frac{h(x)}{h(0)}\right)dx\right]-1}{e^{d/\lambda}-1}$
(146)

The first and second moment of the threshold crossing time are therefore

$w_{0}(y)=L^{-1}[W_{0}(s)]\qquad\overline{t}=-\frac{1}{\gamma v^{*}}\int
w_{0}(y)\ln[h(0)y]dy\qquad\overline{t^{2}}=\frac{1}{(\gamma v^{*})^{2}}\int
w_{0}(y)\ln^{2}[h(0)y]dy$ (147)

Comparing the expression to Eq. 71 we see that the effect of the position
dependence $h(x)$ is equivalent to a change of the cluster size distribution
according to

$\overline{P}_{clu}(s)=\frac{1}{d}\int_{0}^{d}P_{clu}\left(s\frac{h(x)}{h(0)}\right)dx$
(148)

##  15 Bibliography

## References

  * [1] K. T. Son, C. C. Lee, Multiple-Target Laser Range finding Receiver Using a Silicon Photomultiplier Array, IEEE Trans. Instrum. Meas. 59, 3005-3011, (2010)
  * [2] P. Windischhofer, W. Riegler, The statistics of electron-hole avalanches, arXiv:2012.11285, (2020)
  * [3] M. A. Green and M. Keevers, Optical properties of intrinsic silicon at 300 K, Progress in Photovoltaics, Vol. 3, 189-192, (1995)
  * [4] M. A. Green, Self-consistent optical parameters of intrinsic silicon at 300 K including temperature coefficients, Solar Energy Materials & Solar Cells 92, 1305-1310, (2008)
  * [5] C. Canali et al., Electron and hole drift velocity measurements in silicon and their empirical relation to electric field and temperature, IEEE Trans. Electron Dev. 22, 1045, (1975)
  * [6] Synopsis, Inc., Sentaurus Device User Guide Version D-2010.03
  * [7] R. Van Oversträten and H. de Man, Measurement of the ionization rates in diffused silicon p-n junctions, Solid-St. Electron. 13, 583, (1969)
  * [8] G. E. Stillman, C.M. Wolfe, Avalanche Photodiodes, Semiconductors and Semimetals, Chapter 5, Vol. 12, 291-393, (1977)
  * [9] W. Oldham, R. Samuelson, P. Antognetti, Triggering Phenomena in Avalanche Diodes, IEEE, Vol. 19, No. 9, 1056-1060, (1972)
  * [10] R. J. Mcintyre, On the Avalanche Initiation Probability of Avalanche Diodes Above the Breakdown Voltage, IEEE Trans. Electron. Dev., Vol 20, No. 7, 637-641, (1973)
  * [11] L. H. Holway, Electron-Hole Avalanches with Constant Ionization Coefficients, IEEE Trans. Electron. Dev., Vol. 26, No. 6, 991-993, (1979)
  * [12] W. Legler, Die Statistik der Elektronenlawinen in elektronegativen Gasen, bei hohen Feldstärken und bei grosser Gasverstärkung, Z. Naturforschg. 16a, 253-261, (1961)
  * [13] R. Santonico, R. Cardarelli, Development of Resistive Plate Counters, NIMA 187, 377, (1981)
  * [14] A. Mangiarotti, A. Gobbi, On the physical origin of tails in the time response of spark counters, NIMA 482, 192-215, (2002)
  * [15] W. Riegler, Time response functions and avalanche fluctuations in resistive plate chambers, NIMA 602, 377-390, (2009)
  * [16] I. Smirnov, HEED, Program to compute energy loss of fast particles in gases, Version 1.01, CERN
  * [17] H. Schindler, R. Veenhof, et al., Garfield++ simulation of tracking detectors, https://garfieldpp.web.cern.ch/garfieldpp/