Behavioral Model of Silicon Photo-Multipliers Suitable for Transistor-Level Circuit Simulation

Abstract

Silicon Photomultipliers (SiPMs) are photo-electronic devices able to detect single photons and permit the measurement of weak optical signals. Single-photon detection is accomplished through high-performance read-out front-end electronics whose design needs accurate modeling of the photomultiplier device. In the past, a useful model was developed, but it is limited to the device electrical characteristic and its parameter extraction procedure requires several measurement steps. A new silicon photomultiplier model is proposed in this paper. It exploits the Verilog-a behavioral language and is appropriate to transistor-level circuit simulations. The photon detection of a single cell is modeled using the traditional electrical model. A statistical model is included to describe the silicon photomultiplier noise caused by dark-count or after-pulsing effects. The paper also includes a procedure for the extraction of the model parameters through measurements. The Verilog-a model and the extraction procedure are validated by comparing simulations to experimental results.

1. Introduction

Silicon Photomultipliers (SiPMs) are photo-electronic devices, which are gaining more and more importance, able to detect single photons and to allow the measurement of weak optical signals [1,2]. The device can be used in several applications, such as in 3D image sensors [3,4], in gamma spectroscopy [5,6,7,8], with particular emphasis in positron emission tomography (PET) [9,10,11,12,13,14], and, more recently, in visible-light communication (VLC) [15,16].

The SiPM is made up of N parallel-connected avalanche photodiodes (APDs) operating in Geiger-mode. A quenching or ballast resistor is connected in series to each APD, as shown in Figure 1. To detect optical signals, the SiPM requires a biasing and a read-out circuit, as shown in the simplified diagram in Figure 2. In its idle condition, no current crosses the device, and voltage $V_{\mathrm{bias}}$ sets the SiPM operating point at a voltage which, usually, is 10–20% higher than its breakdown voltage, $V_{\mathrm{br}}$. A photon impinging on the active surface of the i-th APD can induce an avalanche current pulse, $I_{\mathrm{av}}\left(t\right)$, which is amplified to the output node as a voltage pulse, $v_{\mathrm{out}}\left(t\right)=G{I}_{\mathrm{av}}\left(t\right)R_s$. If more photons induce k avalanche currents in k different cells, the SiPM will produce a net current pulse proportional to the superposition of k elementary cells with the (ideal) final result of having the output voltage pulse

$$v_{\mathrm{out}}\left(t\right)={\displaystyle \sum _{i=1}^k}G{I}_{\mathrm{av}}\left(t\right)R_s=Gk{I}_{\mathrm{av}}\left(t\right)R_s$$

(1)

In principle, the resulting photosensor is capable of discerning and counting the occurrence of single photons in a pulse of light. In reality, the ideal response of a SiPM in (1) is altered by several sources of noise, such as dark-count, after-pulsing and optical crosstalk. Furthermore, the counting properties or the noise rejection strongly depends on the front-end electronics used to detect the avalanche process.

In this context, transistor-level simulations that examine the interaction between the SiPM device and the remaining components (biasing and read-out circuits) greatly help the designer. Hence, having a reliable electrical model of the SiPM with the possibility of including its noise characteristic, becomes an important aspect of the design phase.

In this paper, we propose an interesting SiPM model that exploits the Verilog-a behavioral language and is appropriate to transistor-level circuit simulations. The photon detection of a single cell is modeled using the traditional electrical model, but the detection of photons in more than one cell is also covered. A statistical model is included to describe the SiPM noise caused by dark-count or after-pulsing effects (The optical crosstalk was not implemented in our model, as it is negligible with new technology processes. However, our Verilog-a model can be extended to include crosstalk effects easily.). After a brief description of the traditional SiPM modeling, we discuss the basic theory of the proposed model in Section 3 and Section 4. Section 5 introduces a procedure for extracting the parameters of the SiPM model from a proper set of measurements. The Verilog-a implementation is discussed in Section 6. Finally, the Verilog-a model and the extraction procedure are validated in Section 7. The conclusions are discussed in Section 8.

2. Traditional SiPM Device Modeling

The literature reports different electrical models for the APD that constitutes the matrix of a SiPM. One of the first electrical models was introduced in [17] and used SPICE to describe the APD as the parallel connection between a capacitor and the series of a voltage source, a controlled switch and a resistor. The model was improved in [18] where, using the Orcad PSpice, the self-quenching mechanism was emulated by means of three controlled switches. These kind of SPICE-like models have been used fruitfully to design active quenching-and-reset circuits (QRCs) for silicon photodiodes such as those reported in [19,20]. Other APD models used the Verilog-a description language to emulate the dynamic response, the self-quenching mechanism and even the statistical operation of the diode [21,22].

However, especially for large arrays, the simulation time can become prohibitively large if a complex model is used for each single APD of the matrix. Hence, a model resulting from the trade-off between accuracy and simulation speed must be used.

Figure 3 shows an accurate and widely accepted SPICE-like electrical model for the SiPM array. It was first introduced in [23] and, afterwards, modified or enhanced. Specifically, some parasitic capacitors were added in [24], the simultaneous firing of multiple cells was incorporated in the models in [25,26] and the read-out electronics were included in the model in [27]. These SPICE-like models have been used successfully to design SiPM front-end circuits such as the current-mode front-end in [28] or the fast read-out circuit in [29].

The SPICE-like model is a good trade-off between accuracy and simplicity, but suffers from the following drawbacks:

1.

Several measurement steps are required for extracting the electrical parameters of the model, thus turning out in a cumbersome procedure [24,25,30,31];

2.

Since the photon arrival is emulated by the current generator, $I_{\mathrm{av}}$, the model cannot self-quench the avalanche of the firing diode (specifically, either the designer uses a complex APD model, such as the one in [18], or the APD model requires that the impulse of the avalanche current is calculated a priori);

3.

No statistical phenomena (i.e., dark-count or after-pulsing [17]) are considered in the model.

The proposed model stems from the schematic in Figure 3 and is appropriate to transistor-level circuit simulators. Moreover, differently from a SPICE-like model, it includes:

1.

Multiple concurrent firing cells;

2.

Avalanche self-quenching;

3.

Emulation of statistical phenomena (i.e., dark-count and after-pulsing effects).

The first two characteristics of the proposed model consent to emulate the situation of two photons that generate two avalanches in a very short interval of time. This correctly allows us to model the after-pulsing phenomena where a secondary avalanche triggers before the cell returns to its equilibrium state.

3. SiPM Electrical Modeling

3.1. Equivalent Circuit Transfer Function

The proposed behavioral model is based on the basic core in Figure 4 [32]. The model embodies N independent current generators so that N independent avalanche pulses can be produced concurrently. In the case of a single firing cell, the model reduces to that in Figure 3.

As the circuit is linear, the overall time response can be obtained through the superposition principle from the time response of a single firing cell. Substituting the electrical model in Figure 3 in the read-out circuit in Figure 2 leads to the circuit shown in Figure 5, whose transfer function is

$$H\left(s\right)=\frac{V_{\mathrm{a}}\left(s\right)}{I_{\mathrm{av}}\left(s\right)}=R_s\frac{b_2{s}^2+b_{1}s+1}{a_3{s}^3+a_2{s}^2+a_{1}s+1}$$

(2)

where

$$\begin{array}{cc}\hfill b_1& =R_q(C_d+2C_q)\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill b_2& =R_q^2{C}_q(C_d+C_q)\hfill \end{array}$$

(3b)

$$\begin{array}{cc}\hfill a_1& =2R_q\left(C_d+C_q\right)+R_s\left(C_g+N{C}_d\right)\hfill \end{array}$$

(3c)

$$\begin{array}{cc}\hfill a_2& =R_q\left[N{C}_d{R}_s\left(C_d+2C_q\right)+R_q{\left(C_d+C_q\right)}^2+2C_g{R}_s\left(C_d+C_q\right)\right]\hfill \end{array}$$

(3d)

$$\begin{array}{cc}\hfill a_3& =R_q^2{R}_s\left(C_d+C_q\right)\left[N{C}_d{C}_q+C_g\left(C_d+C_q\right)\right]\hfill \end{array}$$

(3e)

Under a proper approximation of the denominator of (2), a pole-zero compensation reduces the order of the transfer function that turns into

$$H\left(s\right)\approx R_s\frac{\left(R_q{C}_{q}s+1\right)\left(R_q{C}_{F}s+1\right)}{\left(R_s{C}_{T}s+1\right){\left(R_q{C}_{F}s+1\right)}^2}=R_s\frac{R_q{C}_{q}s+1}{\left(R_s{C}_{T}s+1\right)\left(R_q{C}_{F}s+1\right)}$$

(4)

being

$$\begin{array}{cc}\hfill C_T& =C_g+N(C_d\Vert C_q)\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill C_F& =C_d+C_q\hfill \end{array}$$

(5b)

Concerning the transfer function in (4), we note that $C_T$ stands for the total capacitance seen by $R_s$ for very large $R_q$ (i.e., for $R_q\to \infty$). Similarly, $C_F$ stands for the capacitance seen by the resistor $R_q$ of the active APD for a very small $R_s$ (i.e., for $R_s\sim 0$). Lastly, the branch $C_q$–$R_q$ of the active APD, which directly connects the input to the output, gives rise to the zero. For typical values, the transfer function exhibits a fast time constant ($R_s{C}_T$) and a slow one ($R_q{C}_F$).

3.2. Time Response to a Single Firing APD

Referring to Figure 5, a photon impinging on the surface of the SiPM generates a large current, $i_{\mathrm{av}}$, for a small period of time. Assuming a very small firing time, the avalanche current can be modeled as the ideal impulse

$$i_{\mathrm{av}}=Q_{\mathrm{av}}·\delta \left(t\right)$$

(6)

where

$$Q_{\mathrm{av}}=v_e(0-)\left(C_d+C_q\right)=v_e(0-)C_F$$

(7)

is the charge issued by the active APD of the matrix and $v_e(0-)$ represents the excess voltage (i.e., the voltage across $C_d$ at $t=0-$) [17].

Defining $V_{\mathrm{br}}$ as the APD breakdown voltage, if no cells are active and the circuit is in steady-state, the excess voltage is $v_e(0-)=V_{\mathrm{bias}}-V_{\mathrm{br}}=V_E$. In this case, the output response in the s-domain is

$$V_{\mathrm{a}}\left(s\right)=H\left(s\right)Q_{\mathrm{av}}=\frac{R_s{V}_E{C}_F\left(R_q{C}_{q}s+1\right)}{\left(R_s{C}_{T}s+1\right)\left(R_q{C}_{F}s+1\right)}$$

(8)

whose equivalent time-domain response results

$$v_{\mathrm{a}}\left(t\right)=A_1{\mathrm{e}}^{-\frac{t}{{\tau }_1}}+A_2{\mathrm{e}}^{-\frac{t}{{\tau }_2}}$$

(9)

being

$$\begin{array}{cc}\hfill {\tau }_1& =R_q{C}_F\hfill \end{array}$$

(10a)

$$\begin{array}{cc}\hfill {\tau }_2& =R_s{C}_T\hfill \end{array}$$

(10b)

$$\begin{array}{cc}\hfill A_1& =\frac{R_s{V}_E{C}_F\left({\tau }_1-{\tau }_z\right)}{{\tau }_1\left({\tau }_1-{\tau }_2\right)}=\frac{R_s{C}_d{V}_E}{R_q{C}_F-R_s{C}_T}\hfill \end{array}$$

(10c)

$$\begin{array}{cc}\hfill A_2& =\frac{R_s{V}_E{C}_F\left({\tau }_z-{\tau }_2\right)}{{\tau }_2\left({\tau }_1-{\tau }_2\right)}=\frac{\left(R_q{C}_q-R_s{C}_T\right)\frac{C_F}{C_T}{V}_E}{R_q{C}_F-R_s{C}_T}\hfill \end{array}$$

(10d)

and ${\tau }_z=R_q{C}_q$. At $t=0$, voltage $v_a\left(t\right)$ reaches its maximum value that set the “1-photon” level. It is

$$v_a(0+)=\frac{C_q}{C_T}{v}_e(0-)=\frac{C_q}{C_T}{V}_E=1 \mathrm{ph}$$

(11)

Finally, as far as the firing cell excess voltage, $v_e\left(t\right)$, is concerned, it behaves as the voltage of a charging capacitor characterized by the time constant ${\tau }_1=R_q{C}_F$ with $v_e(0+)=0$ and $v_e(\infty )=V_E$, that is

$$v_e\left(t\right)=V_E\left(1-{\mathrm{e}}^{-\frac{t}{{\tau }_1}}\right)$$

(12)

which defines the recharge phase of the firing cell.

It is worth noting that voltages $v_a\left(t\right)$ and $v_e\left(t\right)$ behave differently. Specifically, since $v_a\left(t\right)$ is dominated by the fast time constant, ${\tau }_2$, it returns to its final value very quickly. Conversely, $v_e\left(t\right)$ is dominated by the slow time constant, ${\tau }_1$, and it returns more slowly to its final value.

4. SiPM Statistical Modeling

The proposed statistical model is based on the theory discussed in [33]. Primary dark-count and after-pulsing phenomena are described in a manner that can be easily implemented in the Verilog-a behavioral language, making the model appropriate to transistor-level circuit simulations.

4.1. Dark-Count

Thermal or field-assisted mechanisms (i.e., band-to-band tunneling, direct and phonon-assisted tunneling, etc.) randomly generate carriers in the APD depletion layer. The number of carriers that are produced in one second and that, under favorable conditions, may produce an undesired avalanche defines the primary dark-count (DC) rate.

The primary dark-count process can be modeled by a Poissonian distribution with average value, $\mathrm{DCR}$. Therefore, two consecutive dark-count events are temporally spaced by the time interval, $\Delta t_{\mathrm{dc}}$, whose exponential distribution is

$$P_{\mathrm{DC}}\left(t\right)=A_{\mathrm{DC}} {\mathrm{e}}^{-\frac{t}{{\tau }_{\mathrm{dc}}}}$$

(13)

being ${\tau }_{\mathrm{dc}}=1/\mathrm{DCR}$ the mean time between two dark-count events. The number of observed DC events can be evaluated from (13) as

$$N_{\mathrm{DC}}={\int }_0^{\infty }{A}_{\mathrm{DC}}{\mathrm{e}}^{-\frac{t}{{\tau }_{\mathrm{dc}}}} \mathrm{d}t=A_{\mathrm{DC}}{\tau }_{\mathrm{dc}}$$

(14)

4.2. Carrier Release and after-Pulsing

When an avalanche occurs, some deep energy levels may entrap a few carriers, which are then released after a characteristic lifetime. If the deep level releases the carrier when the APD is prepared to detect another photon, a new avalanche process may be activated, thus causing an afterpulse correlated with the previous one. Obviously, this after-pulsing (AP) effect alters the detector performance and is strongly undesired [34].

Different deep levels, each with its characteristic lifetime, can determine the specific behavior of an APD. Generally, four deep levels are sufficient to describe the afterpulsing mechanism thoroughly [35]. Nevertheless, a carrier released from a deep level with a lifetime smaller than ${\tau }_1=R_q(C_d+C_q)$ is not able to trigger an avalanche (in fact, the APD has not yet finished the quenching process) and, therefore, only the slowest levels play a significant role.

The AP process is modeled as follows. We assume that, after a primary avalanche pulse, a carrier may be trapped into a deep level with a constant probability, $P_{\mathrm{trap}}$. Hence, with respect to $N_{\mathrm{DC}}$, the number of trapped carriers is

$$N_{\mathrm{trap}}=N_{\mathrm{DC}}{P}_{\mathrm{trap}}$$

(15)

Once trapped, the carrier is randomly released. Specifically, the interval of time between the primary pulse and the carrier release event, $\Delta t_{\mathrm{cr}}$, follows the exponential probability density function (pdf)

$$p_{\mathrm{CR}}\left(t\right)=A_{\mathrm{CR}} {\mathrm{e}}^{-\frac{t}{{\tau }_{\mathrm{cr}}}}$$

(16)

where ${\tau }_{\mathrm{cr}}=1/A_{\mathrm{CR}}$ is the trap lifetime. The triggering probability, $P_{\mathrm{trig}}$, establishes the probability of a carrier released from a deep level to generate an avalanche pulse. The triggering probability strongly depends on the excess voltage, $v_e$, and, although a complete model is reported in [36], a commonly used expression for this parameter is

$$P_{\mathrm{trig}}=1-{\mathrm{e}}^{-\frac{v_e}{{\eta }_T{V}_{\mathrm{br}}}}\hspace{1em}\mathrm{for}{v}_e\ge 0$$

(17)

where ${\eta }_T$ is an empirical parameter [37,38]. In our model, we further simplify this expression assuming $v_e$ small with respect to ${\eta }_T{V}_{\mathrm{br}}$, so that the triggering probability reduces to

$$P_{\mathrm{trig}}\approx \frac{v_e}{{\eta }_T{V}_{\mathrm{br}}}\hspace{1em}\mathrm{for}{v}_e\ge 0$$

(18)

During the recharge phase, the excess voltage is not constant and is defined by (12), so that the triggering probability changes with time. Considering this situation, the distribution related to the after-pulsing process depends on the number of traps, on the probability to release a carrier and on the probability to trigger an avalanche, that is

$$P_{\mathrm{AP}}\left(t\right)=N_{\mathrm{trap}}·p_{\mathrm{CR}}\left(t\right)·P_{\mathrm{trig}}\left(t\right)=\frac{N_{\mathrm{DC}}{P}_{\mathrm{trap}}{A}_{\mathrm{CR}}{V}_E}{{\eta }_T{V}_{\mathrm{br}}}\left(1-{\mathrm{e}}^{-\frac{t}{{\tau }_1}}\right){\mathrm{e}}^{-\frac{t}{{\tau }_{\mathrm{cr}}}}$$

(19)

Assume that, at $t=0$, a primary pulse occurs followed, during the recharge phase, by a subsequent afterpulse at $t=t^{*}$. In this situation, (12) defines the excess voltage, $v_e$, and, using (11), we obtain for the maximum output voltage

$$v_a\left(t^{*}\right)=\frac{C_q}{C_T}{v}_e\left(t^{*}\right)=\frac{C_q}{C_T}{V}_E\left(1-{\mathrm{e}}^{-\frac{t^{*}}{{\tau }_1}}\right)$$

(20)

a value that is lower than the 1-photon level.

The counting circuit will detect a photon if the output voltage exceeds the threshold of the read-out circuit, $V_{\mathrm{th}}=\alpha (C_q/C_T)V_E$, which is defined as $\alpha$ times the 1-photon level, with $\alpha <1$. However, in the case of an after-pulsing event, if the output voltage remains lower than the threshold, the read-out circuit is unable to detect the avalanche. This happens for any after-pulsing event that occurs before time ${\tau }_{\mathrm{th}}$ so that $v_a\left({\tau }_{\mathrm{th}}\right)<V_{\mathrm{th}}$. It is easy to show that

$${\tau }_{\mathrm{th}}={\tau }_{1}ln\left(\frac{1}{1-\alpha }\right)$$

(21)

Since no after-pulsing event is detected for $t<{\tau }_{\mathrm{th}}$, the overall dark-count noise distribution is evaluated by summing (13) to (19), that is

$$P_N\left(t\right)=A_{\mathrm{DC}}{\mathrm{e}}^{-\frac{t}{{\tau }_{\mathrm{dc}}}}+u\left(t-{\tau }_{\mathrm{th}}\right)A_{\mathrm{AP}}\left(1-{\mathrm{e}}^{-\frac{t}{{\tau }_1}}\right){\mathrm{e}}^{-\frac{t}{{\tau }_{\mathrm{cr}}}}$$

(22)

where $u\left(t\right)$ is the Heaviside step function and

$$A_{\mathrm{AP}}=\frac{N_{\mathrm{DC}}{P}_{\mathrm{trap}}{A}_{\mathrm{CR}}{V}_E}{{\eta }_T{V}_{\mathrm{br}}}=\frac{A_{\mathrm{DC}}{\tau }_{\mathrm{dc}}{P}_{\mathrm{trap}}{V}_E}{{\eta }_T{V}_{\mathrm{br}}{\tau }_{\mathrm{cr}}}$$

(23)

5. Parameter Extraction

5.1. Electrical Model

The extraction procedure starts with the evaluation of the quenching resistor, $R_q$. When the SiPM in Figure 1 is biased in the forward region, the static I-V characteristic is

$$V=R_q\frac{I}{N}+n{V}_{T}ln\left(\frac{I}{N{I}_s}\right)$$

(24)

where n is the junction emission coefficient and with the usual meaning for the remaining parameters. The evaluation of $R_q$ is easily accomplished through a fitting procedure.

For the second step of the extraction procedure, we make use of the SiPM connected as in Figure 2. In this procedure, we measure and record the curves of the output voltage, $v_{\mathrm{out}}\left(t\right)$, obtained when the avalanche pulses are generated by shots of single APDs. To ensure that only one APD is generating the avalanche, the SiPM is maintained under dark condition and only dark-count events are recorded. A digital oscilloscope is sufficient to do the job. Then, we find the parameters ${\tau }_1$, ${\tau }_2$, $A_1$ and $A_2$ of the time-domain response in (9) by a fitting procedure with the curves of the measured voltage, $v_{\mathrm{a}}=v_{\mathrm{out}}/G$, being G the known gain of the pre-amplifier in Figure 2.

Next, using (10), we find ${\tau }_z$

$${\tau }_z=\frac{{\tau }_1{\tau }_2\left(1+\frac{A_1}{A_2}\right)}{{\tau }_2+{\tau }_1\frac{A_1}{A_2}}$$

(25)

and capacitors $C_q$, $C_d$ and $C_g$, that result in

$$\begin{array}{cc}\hfill C_q& =\frac{{\tau }_z}{R_q}\hfill \end{array}$$

(26a)

$$\begin{array}{cc}\hfill C_d& =C_F-C_q=\frac{{\tau }_1}{R_q}-C_q\hfill \end{array}$$

(26b)

$$\begin{array}{cc}\hfill C_g& =C_T-N(C_d\Vert C_q)=\frac{{\tau }_2}{R_s}-N(C_d\Vert C_q)\hfill \end{array}$$

(26c)

In conclusion, the excess voltage is

$$V_E=\frac{R_q{C}_F-R_s{C}_T}{R_s{C}_d}{A}_1$$

(27)

and the breakdown voltage results in

$$V_{\mathrm{br}}=V_{\mathrm{bias}}-V_E$$

(28)

5.2. Statistical Model

The SiPM statistical model relies on relationship (22). From the simulator point of view, we can characterize the SiPM noise by specifying $P_{\mathrm{trap}}$, ${\eta }_T$, $V_{\mathrm{br}}$, ${\tau }_{\mathrm{dc}}$, ${\tau }_1$ and ${\tau }_{\mathrm{cr}}$. It is worth noting that other parameters (i.e., $V_E$ and ${\tau }_{\mathrm{th}}$) are also responsible for the overall noise distribution of the detector, even if they do not model the noise properties of the intrinsic device.

The procedure for extracting the parameters of the statistical model requires the evaluation of $V_{\mathrm{br}}$ and ${\tau }_1$ as described in Section 5.1. Then, choosing the excess voltage, $V_E$, and the threshold $\alpha$ so that ${\tau }_{\mathrm{th}}$ is set as in (21), we record the time intervals between two consecutive shots and evaluate the corresponding distribution. A proper fitting procedure allows us to find the remaining parameters in (22), in particular, $A_{\mathrm{DC}}$, ${\tau }_{dc}$, $A_{\mathrm{AP}}$ and ${\tau }_{\mathrm{cr}}$. The steps to measure and record the time intervals and to obtain the corresponding distribution is reported in [1] and is not included in this text, so we assume the distribution as known.

Finally, from (23), we can associate the remaining unknown parameters, $P_{\mathrm{trap}}$ and ${\eta }_T$, as

$$P_{\mathrm{trap}}=\frac{A_{\mathrm{AP}}{\tau }_{\mathrm{cr}}}{A_{\mathrm{DC}}{\tau }_{\mathrm{dc}}}\times \frac{{\eta }_T{V}_{\mathrm{br}}}{V_E}$$

(29)

Observe that $V_E/\left({\eta }_T{V}_{\mathrm{br}}\right)$ stands for the triggering probability when the excess voltage is at its nominal value, $V_E$.

6. Verilog-a Model

We used the description language Verilog-a to implement the proposed SiPM model. Verilog-a is the analog extension of the well-known Verilog HDL (Hardware Description Language) [39] and is a standard tool in industrial design environments for integrated circuits. More specifically, all the most widespread simulation tools (ELDO, HSPICE, SPECTRE, etc.) can perform mixed behavioral/transistor-level simulations that makes the proposed model appropriate to several simulation and design contexts.

The main core of the Verilog-a code is reported in Figure 6. The SiPM device is modeled as a four-terminal module, a, k, phN and phT, where the last two terminals are used to emulate the arrival of multiple photons. The number of APD cells that compose the SiPM matrix is set by the macro N in line 2.

The section “Parameters and variable declaration” defines and the initializes several model parameters. The model is defined by electrical (Rq, Cd, Cq, Cg, Vbr, delta) and statistical parameters (Ptrap, etaT, tau_dc, tau_cr). Parameter delta is used in the @(initial_step) section (lines 12–16) to define the avalanche current as

$$I_{\mathrm{av}}=\frac{Q_{\mathrm{av}}}{\delta }=\frac{(C_d+C_q)V_E}{\delta }$$

(30)

Note, however, that this parameter is used just to evaluate the order of magnitude for $I_{\mathrm{av}}$, as the actual turn-off depends on the circuit evolution.

The “Analog description code” performs the electrical and the statistical model described in Section 3 and Section 4. Here, the @(initial_step) section (lines 12–16) initializes all the required variables.

Referring to Figure 4, the codes inside the two generate loops (lines 18–19 and 27–39) are replicated N times and are written for the generic j-th cell. Hence, the electrical behavior of each cell is described in lines 37–38 by

$$\begin{array}{cc}\hfill I_{k,c_j}\left(t\right)& =\frac{V_{k,c_j}\left(t\right)}{R_q}+C_q\frac{\mathrm{d}{V}_{k,c_j}\left(t\right)}{\mathrm{d}t}\hfill \end{array}$$

(31a)

$$\begin{array}{cc}\hfill I_{c_j,a}\left(t\right)& =a_{v_j}{I}_{av}+C_d\frac{\mathrm{d}{V}_{c_j,a}\left(t\right)}{\mathrm{d}t}\hfill \end{array}$$

(31b)

while the contribution of the grid capacitor $C_g$ is set in line 41, outside the generate loop, as

$$I_{k,a}\left(t\right)=C_g\frac{\mathrm{d}{V}_{k,a}\left(t\right)}{\mathrm{d}t}$$

(32)

In the code, the j-th APD cell is turned on/off by setting/resetting the flag $a_{v_j}$. More specifically, the j-th APD cell is quenched as modeled in the “Turn-off” section, that is when the excess voltage, $v_{e_j}$ crosses the zero in the falling edge (lines 30–31). The avalanche in an APD cell can be turned on because of a photon arrival, a dark-count or an after-pulsing event. These three events are modeled in the respective sections of code and are described below.

6.1. Photon Arrival

The “Photon arrival” section of the Verilog-a code in Figure 6 is expanded in Figure 7. The code is executed when the voltage at terminal phT crosses $1 \mathrm{V}$ in the positive edge (line 2). The number of photons, and the corresponding number of cells to be triggered, is set by the voltage at terminal phN and, after truncation, is stored in the variable Nph (line 4). If the number of arriving photons is higher than the number of APD cells of the SiPM matrix, the device is in saturation and all the cells are prepared to be turned-on by setting all the N temporary avalanche flags, $a_{v{t}_i}$, to ‘1’ (lines 5–6). Conversely, if the number of arriving photons is lower than the number of cells, Nph cells are chosen randomly and are prepared to be turned on by setting the respective temporary avalanche flags $a_{v{t}_i}$ to ‘1’ (lines 8–20). In lines 22–31, for each cell to be triggered, after checking that no avalanche is taking place and that the APD is biased above breakdown, the avalanche is started. Meanwhile, with a probability established by the trapping probability, $P_{\mathrm{trap}}$, a new carrier release time is stored in the variable $t_{{\mathrm{cr}}_i}$ for the i-th cell. The carrier release time is chosen from an exponential distribution with mean ${\tau }_{\mathrm{cr}}$ and shall be used for executing the “After-pulsing” code, as described in Section 6.3.

6.2. Dark-Count

The “Dark-count” section of the Verilog-a code in Figure 6 is expanded in Figure 8. The code is executed when the absolute simulation time equals the value assigned to the variable $t_{\mathrm{dc}}$ (line 2). First, the cell to be turned on is chosen randomly (line 4). Then, if no avalanche is taking place and the i-th APD is biased above breakdown, the avalanche process is started setting $a_{v_i}=1$. Meanwhile, with a probability set by the trapping probability, $P_{\mathrm{trap}}$, a new carrier release time is stored in the variable $t_{{\mathrm{cr}}_i}$ using an exponential distribution with mean ${\tau }_{\mathrm{cr}}$ (lines 5–12). Finally, a new the simulation time for the next dark-count event is stored in the variable $t_{\mathrm{dc}}$ using an exponential distribution with mean ${\tau }_{\mathrm{dc}}$ (line 14). Lines 17–19 are required to schedule the next dark-count event correctly.

6.3. After-Pulsing

Referring to the main code in Figure 6, the “After-pulsing” section is replicated N times thanks to the generate loop in lines 18–39. The corresponding Verilog-a code is reported in Figure 9 and refers to the j-th cell. The code for the j-th cell is executed when the absolute simulation time equals the value assigned to the variable $t_{{\mathrm{cr}}_j}$ (line 2). If no avalanche is taking place and the j-th APD is biased above breakdown, the triggering probability, $P_{\mathrm{trig}}$, is evaluated and used to start the avalanche and schedule a new carrier release event at time $t_{{\mathrm{cr}}_j}$ with a probability set by the trapping probability, $P_{\mathrm{trap}}$ (lines 6–12). Lines 17–19 are required to schedule the next carrier release event correctly.

7. Model Validation

The proposed Verilog-a model was verified using two 10 × 10 SiPM prototypes provided by STMicroelectronics. We shall refer to them as SiPM-A and SiPM-B. They were used for testing the electrical and the statistical model, respectively.

The photo of the SiPM-A device is shown in Figure 10. The value of the quenching resistor was determined with the procedure in Section 5.1. The static I-V characteristic of the SiPM was fitted with (24), as depicted in Figure 11, thus obtaining $R_q=292.6 \mathrm{k}\mathsf{\Omega }$. The other electrical parameters were determined with the detector assembled as in Figure 2 with $V_{\mathrm{bias}}=31 \mathrm{V}$ and $R_s=25 \mathsf{\Omega }$. The pre-amplifier exhibits a 46-dB gain and a 4-GHz high-frequency bandwidth. Voltage $v_{\mathrm{out}}$ was measured and recorded so to obtain a set of 500 curves. Following the procedure in Section 5.1, the measured data $v_{\mathrm{a}}=v_{\mathrm{out}}/G$ were fitted with the double-exponential function in (9) and the parameters in

Table 1. SiPM-A fitted parameters from experimental data.

Parameter	Value
${\tau }_1$	$49.96 \mathrm{ns}$
${\tau }_2$	$2.89 \mathrm{ns}$
$A_1$	$208.9 \mathsf{\mu }\mathrm{V}$
$A_2$	$423.9 \mathsf{\mu }\mathrm{V}$

were determined. Finally, the electrical parameters reported in

Table 2. Extracted parameters of the electrical model for the SiPM-A device.

Parameter	Value
$R_q$	$292.6 \mathrm{k}\mathsf{\Omega }$
$C_q$	$26.78 \mathrm{fF}$
$C_d$	$143.96 \mathrm{fF}$
$C_g$	$113.37 \mathrm{pF}$
$V_{\mathrm{br}}$	$28.27 \mathrm{V}$
$Q_{\mathrm{av}}$	$466.5 \mathrm{fC}$

were obtained.

The SiPM-A model was implemented in Verilog-a using the code described in Section 6. The avalanche current, $I_{\mathrm{av}}$, was set as in (30) using $\delta =10 \mathrm{ps}$. The SiPM-A device was simulated, and the time response was compared with a second set of measurement data (500 curves), as shown in Figure 12. Apart from the gradual starting slope of the measured data, caused by the finite bandwidth of the oscilloscope, it is apparent that the simulated response matches very well with the measured results, thus validating both the model and the extraction procedure.

As far as the SiPM-B is concerned, following the same procedure used for the SiPM-A device, we set $V_{\mathrm{bias}}=31.5 \mathrm{V}$, $R_s=25 \mathsf{\Omega }$ and extracted the electrical parameters summarized in

Table 3. Extracted parameters of the electrical model for the SiPM-B device.

Parameter	Value
$R_q$	$1062 \mathrm{k}\mathsf{\Omega }$
$C_q$	$171.43 \mathrm{fF}$
$C_d$	$34.286 \mathrm{fF}$
$C_g$	$338 \mathrm{pF}$
$V_{\mathrm{br}}$	$29.5 \mathrm{V}$
$Q_{\mathrm{av}}$	$411.43 \mathrm{fC}$

.

Then, maintaining the same excess voltage, $V_E=2 \mathrm{V}$, and setting the threshold of the read-out circuit to 0.5-photon level ($\alpha =1/2$), we recorded the time intervals between two consecutive shots and evaluated the corresponding distribution, reported in Figure 13 (measured data).

Using the data in Table 3, we computed ${\tau }_1=218.4 \mathrm{ns}$ and ${\tau }_{\mathrm{th}}=151.5 \mathrm{ns}$ from (10a) and (21), respectively. Then, we fit the noise distribution in (22) to the measured data and extracted the statistical parameters shown in

Table 4. Extracted parameters of the statistical model for the SiPM-B device.

Parameter	Value
$A_{\mathrm{DC}}$	$1617.8$
$A_{\mathrm{AP}}$	$638.4$
${\tau }_{\mathrm{dc}}$	$2658 \mathrm{ns}$
${\tau }_{\mathrm{cr}}$	$187.8 \mathrm{ns}$
${\tau }_{\mathrm{th}}$	$151.5 \mathrm{ns}$
${\tau }_1$	$218.5 \mathrm{ns}$

. The noise distribution obtained from the fitting procedure is plotted in Figure 13 (measurement fit) and matches very well the measured data with a goodness of fit set by ${\chi }^2/\mathrm{NDF}=1.00$.

We used the Verilog-a code described in Section 6 to implement the simulation model for the SiPM-B device. In this model, we imposed $P_{\mathrm{trig}}=0.5$, when $V_E=2 \mathrm{V}$, that is, when the SiPM is biased at its nominal excess voltage. To do so, we simply set ${\eta }_T=0.13559$ from (18). We also obtained $P_{\mathrm{trap}}=0.05575$ from (29).

The statistical simulation model was validated as follows. We assembled the circuit as in Figure 2 ($G=46 \mathrm{dB}$) and ran a transient simulation of 180 ms recording the time intervals between two consecutive shots that exceeded the 0.5-photon level. In the process, no photon arrival was emulated, so that the evolution remained characterized by the noise of the SiPM only (i.e., dark-count and after-pulsing events). A portion of the transient simulation is reported in Figure 14, where the output voltage of the pre-amplifier in Figure 2 is considered. The 1-photon level is about 40 mV so that we recorded all the signals that exceeded the 20-mV threshold.

Figure 15 shows the different behaviors of two after-pulsing events both caused by antecedent dark-count pulses. In the case depicted in Figure 15a, the AP event occurs after ${\tau }_{\mathrm{th}}$ seconds with respect to the antecedent DC pulse. Hence, as predicted from (20), its maximum output voltage, although lower than the 1-photon level, is still able to overcome the 20-mV threshold and is regularly detected. Conversely, in Figure 15b, with respect to the DC event, the AP pulse occurs before ${\tau }_{\mathrm{th}}$ seconds. The SiPM generates a maximum voltage lower than the 20-mV threshold, and the event is not recorded.

The overall simulation produced 68475 pulses that were organized in the discrete distribution reported in Figure 16 (simulation data). Then, we fitted the noise distribution in (22) to the simulation data and obtained the values in

Table 5. SiPM-B fitted parameters from simulation data.

Parameter	Value
$A_{\mathrm{DC}}$	$1311.1$
$A_{\mathrm{AP}}$	$502.6$
${\tau }_{\mathrm{dc}}$	$2664 \mathrm{ns}$
${\tau }_{\mathrm{cr}}$	$186.6 \mathrm{ns}$
${\tau }_{\mathrm{th}}$	$151.1 \mathrm{ns}$
${\tau }_1$	$218.5 \mathrm{ns}$

, with a goodness of fit set by ${\chi }^2/\mathrm{NDF}=1.00$.

Using the values in Table 5, in Figure 16 we plotted the noise distribution obtained from the fitting of the simulation results (simulation fit). In the same figure, we also plotted the noise distribution obtained using the values in Table 4, that is, from the fitting of the measured data (measured fit). This latter figure required a proper scaling of $A_{\mathrm{DC}}$ and $A_{\mathrm{AP}}$ to have a fair comparison. The percentage error between the two noise distributions (measurement and simulation fit) maintains below 0.47%. The excellent matching between the simulation data and both the measurement and the simulation fits provided a complete validation of the statistical model as well as of its Verilog-a implementation and of the extraction procedure.

8. Conclusions

In this paper, we proposed and developed an interesting SiPM model appropriate to transistor-level simulation. The model exploits the description language Verilog-a, the analog extension of the well-known Verilog HDL (Hardware Description Language), a standard tool in industrial design environments for integrated circuits. Since all the most widespread simulation tools (ELDO, HSPICE, SPECTRE, etc.) can perform mixed behavioral/transistor-level simulations, the proposed model is suitable to several simulation and industrial design contexts. The model makes use of the traditional electrical model to describe the shot of one or more cells of the matrix. A statistical model is also included that allows to emulate the noise behavior in terms of dark-count and after-pulsing phenomena. The theoretical background of the model is introduced. Then, the Verilog-a implementation of the model and the procedure for extracting all the necessary parameters are also discussed. Lastly, both the extraction procedure and the model are validated from experimental data.