Modelling the Nonlinear Response of Silicon Photomultipliers

TL;DR

This paper introduces a comprehensive statistical model for the nonlinear response of silicon photomultipliers, accounting for various noise sources and recovery effects, validated against experimental data for different SiPMs.

Contribution

It presents a novel, detailed model that accurately predicts SiPM responses under various conditions, improving understanding of nonlinear behaviors.

Findings

Model accurately predicts SiPM response at moderate nonlinearity

Validated with experimental data from different SiPMs and light sources

Provides insights into effects of overvoltage and pulse shape on response

Abstract

A statistical model of the nonlinear response of silicon photomultipliers is presented. It includes losses of both the photodetection efficiency and the gain during pixel recovery periods as well as the effect of correlated and uncorrelated noise. The model provides either the mean output charge of a SiPM for incident light pulses of arbitrary shape or the output current for continuous light. The dependence of the SiPM response on both the overvoltage and the pulse shape is also properly described. The model has been validated for two different silicon photomultipliers using scintillation light pulses from a LYSO crystal as well as continuous light from a LED. Good agreement is found with experimental data at moderate nonlinearity.

Tables

Table 1. TABLE I: Characteristics of the Tested SiPMs.

Parameters	Symbol	1325	1350	Unit
SiPM area	-	$1.3\times 1.3$		mm2
Pixel pitch	-	25	50	$\mu$m
Number of pixels	$N$	2668	667	-
Breakdown voltage	$V_{\mathrm{br}}$	51.80		V
Recom. op. overvoltage	$U_{\mathrm{op}}$	5.00	3.00	V
Dark count rate	$r_{\mathrm{dc}}$			kcps
Avalanche charge	$q$	0.7	1.7	$\times {10}^6$ e
Charact. overvoltage	$U_{\mathrm{ch}}$	2.69	2.68	V
Overvoltage shift	$U_0$	0.66	0.00	V
Maximum PDE	${\epsilon }_{\max }$	0.327	0.597	-
Recovery time	$t_{\mathrm{rec}}$	17	29	ns
Cor. noise probability	-	0.03	0.01	-

Equations

$$\epsilon ={\epsilon }_{\max }\cdot \left[1-\exp \left(-\frac{U-U_0}{U_{\mathrm{ch}}}\right)\right]$$

$$⟨Q⟩=N\cdot q\cdot \left[1-\exp \left(-\frac{\epsilon \cdot n}{N}\right)\right],$$

$$I=\frac{\epsilon \cdot r\cdot q}{1+\epsilon \cdot r\cdot t_{\mathrm{dead}}/N},$$

$$u(t)=U\cdot \left[1-\exp \left(-\frac{t}{t_{\mathrm{rec}}}\right)\right],$$

$${\epsilon }_{\mathrm{rec}}(t)$$

$$q_{\mathrm{rec}}(t)$$

$$P(m)=\frac{(\epsilon \cdot n{)}^m}{m!}\cdot e^{-\epsilon \cdot n},$$

$$F(Q)=\sum _{m=0}^{\infty }P(m)\cdot F_m(Q),$$

$$⟨Q⟩=\sum _{m=0}^{\infty }P(m)\cdot ⟨Q{⟩}_m,$$

$$⟨Q{⟩}_m={\int }_0^{\infty }{F}_m(Q)\cdot Q\cdot \mathrm{d}Q.$$

$$F_2(Q)={\int }_0^{Q}f(Q-q_s)\cdot f(q_s)\cdot \mathrm{d}{q}_s,$$

$$\alpha =\frac{2}{\epsilon }\cdot {\int }_0^{\infty }{\int }_{t_0}^{\infty }p(t)\cdot p(t+t_s)\cdot {\epsilon }_{\mathrm{rec}}(t_s)\cdot \mathrm{d}{t}_s\cdot \mathrm{d}t,$$

$$\begin{gathered} \varphi (q_s)=\frac{2}{\alpha \cdot \epsilon }\cdot {\int }_0^{\infty }{\int }_{t_0}^{\infty }p(t)\cdot p(t+t_s) \\ \cdot {\epsilon }_{\mathrm{rec}}(t_s)\cdot f_s(t_s,q_s)\cdot \mathrm{d}{t}_s\cdot \mathrm{d}t, \end{gathered}$$

$$\gamma =\frac{2}{\epsilon \cdot q}\cdot {\int }_0^{\infty }{\int }_{t_0}^{\infty }p(t)\cdot p(t+t_s)\cdot {\epsilon }_{\mathrm{rec}}(t_s)\cdot q_{\mathrm{rec}}(t_s)\cdot \mathrm{d}{t}_s\cdot \mathrm{d}t.$$

$$\begin{gathered} F_2(Q)={\int }_0^Q\left(1-\frac{1}{N}\right)\cdot f(Q-q_s)\cdot f(q_s)\cdot \mathrm{d}{q}_s \\ +\frac{1-\alpha }{N}\cdot f(Q)+{\int }_0^Q\frac{\alpha }{N}\cdot f(Q-q_s)\cdot \varphi (q_s)\cdot \mathrm{d}{q}_s. \end{gathered}$$

$$⟨Q{⟩}_2=2\cdot q\cdot \left(1-\frac{1}{N}\right)+\frac{(1+\gamma )\cdot q}{N},$$

$$F_m(Q)=$$

$$+\frac{(1-\alpha )\cdot Q}{N\cdot q}\cdot F_{m-1}(Q)$$

$$+{\int }_0^Q\frac{\alpha \cdot (Q-q_s)}{N\cdot q}\cdot F_{m-1}(Q-q_s)\cdot \varphi (q_s)\cdot \mathrm{d}{q}_s,$$

$$⟨Q{⟩}_m=\frac{N\cdot q}{1-\gamma }\cdot \left[1-{\left(1-\frac{(1-\gamma )\cdot (1+c)}{N}\right)}^m\right],$$

$$⟨Q⟩=\frac{N\cdot q}{1-\gamma }\cdot \left[1-\exp \left(-(1-\gamma )\cdot (1+c)\cdot \frac{\epsilon \cdot n}{N}\right)\right].$$

$$I=\frac{N\cdot q}{2\cdot t_{\mathrm{dead}}}\cdot \left[1-\exp \left(-2\cdot t_{\mathrm{dead}}\cdot (1+c)\cdot \frac{\epsilon \cdot r}{N}\right)\right].$$

$$(1+c)\cdot n=(1+c)\cdot\eta\cdot Y\cdot E=k_{E}\cdot E\$$

$$p(t)=\{\begin{array}{ll}\frac{1}{T}& 0<t<T\\ 0& \mathrm{otherwise}.\end{array}$$

$$\gamma =\frac{2}{T^2\cdot \epsilon \cdot q}\cdot {\int }_0^{T-t_0}{\int }_{t_0}^{T-t}{\epsilon }_{\mathrm{rec}}(t_s)\cdot q_{\mathrm{rec}}(t_s)\cdot \mathrm{d}{t}_s\cdot \mathrm{d}t.$$

$$\gamma \approx \frac{{\epsilon }_{\max }\cdot U}{6\cdot \epsilon \cdot U_{\mathrm{ch}}}\cdot \frac{(T+t_0)\cdot (T-t_0{)}^3}{T^2\cdot t_{\mathrm{rec}}^2}.$$

$$\gamma \approx 1-\frac{2\cdot t_{\mathrm{dead}}}{T},$$

$$\begin{gathered} t_{\mathrm{dead}}=t_0+{\int }_{t_0}^{\infty }\left[1-\frac{{\epsilon }_{\mathrm{rec}}(t_s)\cdot q_{\mathrm{rec}}(t_s)}{\epsilon \cdot q}\right]\cdot \mathrm{d}{t}_s \\ =t_0+t_{\mathrm{rec}}\cdot [\frac{{\epsilon }_{\max }}{\epsilon }\cdot e^{-t_0/t_{\mathrm{rec}}}-\frac{U_{\mathrm{ch}}}{U} \\ +\left(\frac{{\epsilon }_{\max }}{\epsilon }-1\right)\cdot \sum _{n=1}^{\infty }\frac{{\left(U/U_{\mathrm{ch}}\right)}^n}{n\cdot n!}\cdot {\left(e^{-t_0/t_{\mathrm{rec}}}\right)}^n]. \end{gathered}$$

$$p(t)=\frac{1}{{\tau }_2-{\tau }_1}\cdot \left(e^{-t/{\tau }_2}-e^{-t/{\tau }_1}\right)$$

$$\gamma =\frac{{\tau }_2^2\cdot {\gamma }_2-{\tau }_1^2\cdot {\gamma }_1}{{\tau }_2^2-{\tau }_1^2}$$

$${\gamma }_i=\frac{1}{{\tau }_i\cdot \epsilon \cdot q}\cdot {\int }_{t_0}^{\infty }{e}^{-t_s/{\tau }_i}\cdot {\epsilon }_{\mathrm{rec}}(t_s)\cdot q_{\mathrm{rec}}(t_s)\cdot \mathrm{d}{t}_s(i=1,2).$$

$${\gamma }_i\approx \frac{{\epsilon }_{\max }\cdot U}{\epsilon \cdot U_{\mathrm{ch}}}\cdot \frac{{\tau }_i\cdot (2\cdot {\tau }_i+t_0)}{t_{\mathrm{rec}}^2}\cdot e^{-t_0/{\tau }_i}.$$

Full text

I Introduction

The silicon photomultiplier (SiPM) is a high-sensitivity semiconductor photodetector that is increasingly becoming the best choice in many applications owing to its compactness, high gain, insensitivity to magnetic fields, rapid response and exceptional photon-counting resolution. This device consists in an array of Geiger-mode avalanche photodiodes in the $\mu$m scale, hereafter called pixels, which are sensitive to single photons. The pixels are connected in parallel in such a way the SiPM output signal is the sum of the signals from all the pixels. In this work, the SiPM response for an incident light pulse is characterized by the mean output charge, which is given by $\langle Q\rangle=\varepsilon\cdot n\cdot q$ for an ideal SiPM, where $\varepsilon$ is the photodetection efficiency, $n$ is the mean number of photons per pulse impinging on the SiPM surface and $q$ is the mean charge released by a single breakdown avalanche, which characterizes the gain of the device.

The mean avalanche charge $q$ is proportional to the overvoltage $U$, which is defined as the reverse-bias voltage $V_{\rm bias}$ minus the breakdown voltage $V_{\rm br}$ of the SiPM. Both $q$ and $V_{\rm br}$ can be obtained from the output pulse charge spectrum at photon counting levels (see, e.g., [1]). For typical blue sensitive SiPMs based on a p-on-n structure, $\varepsilon$ is described by an exponential function of $U$ [2]. In this work it is assumed

[TABLE]

for $U\geq U_{0}$, otherwise $\varepsilon=0$. Here $\varepsilon_{\rm max}$ is the product of the quantum efficiency and the fill factor (i.e., the ratio of the photosensitive area of a pixel to the entire pixel area), the parameter $U_{\rm ch}$ is a characteristic overvoltage of the SiPM depending on the input wavelength spectrum [2], and $U_{0}$ is a small overvoltage shift that is introduced to better fit experimental data of $\varepsilon$ (see [3, 4] for an interesting studies on the breakdown behaviour of SiPMs). It should be clarified that this empirical expression is not applicable to n-on-p SiPMs [5].

The SiPM response is nonlinear for high photon densities because the number of pixels is finite (typically of the order of 1000, depending on the SiPM) and they take a few tens of ns to recharge after each breakdown avalanche. Besides, SiPMs have uncorrelated noise due to the thermal production of electron-hole pairs in silicon, as well as correlated noise caused by two processes called afterpulsing and crosstalk. Afterpulses are stochastic parasitic avalanches produced in the same pixel where a previous avalanche has been triggered. On the other hand, crosstalk avalanches are produced in nearby pixels either simultaneously with the primary avalanche (prompt crosstalk) or with some delay (delayed crosstalk) [1].

The combination of all these features makes the statistics of avalanche triggering and charge production very complex. As a consequence, an exact statistical description of the response of SiPMs is only possible via Monte Carlo simulations. H.T. van Dam et al. [6] developed a comprehensive statistical model of the SiPM response for exponentially decaying light pulses, accounting for both the correlated noise and the pixel recovery under some approximations. However it requires non-trivial numerical calculations. A simpler statistical model was presented in [7], but it ignores the correlated noise and accounts of the recovery of pixels in an incomplete way. In addition, analytical expressions of the SiPM response in two limit situations are available. First, for light pulses much shorter than the pixel recovery time and in the absence of correlated and uncorrelated noise, the mean output charge $\langle Q\rangle$ is given by the well-known expression (see, e.g., [8])

[TABLE]

where $N$ is the number of pixels. Second, for continuous light and ignoring again the correlated noise, the output current intensity $I$ produced by the SiPM can be approximated by

[TABLE]

where $r$ is the rate of impinging photons and $t_{\rm dead}$ is a certain non-paralyzable dead time accounting for the pixel recovery [9]. Modified versions of the above expressions are being used to describe the response of SiPMs at particular conditions, e.g., [10, 11, 12, 13, 14]. However, a general analytical model is not available yet.

In this paper, a simple model that includes all the above features under some simplifications is described. A preliminary version of this model was first reported in [15], but significant upgrades have been made since then. The model has been validated against experimental data for scintillation light pulses from a LYSO crystal as well as for continuous light from a LED.

II The model

II-A Approach of the problem

In steady state, all the pixels of a SiPM are biased with the supplied overvoltage $U$. Therefore, the probability that an impinging photon fires a pixel is equal to the nominal photodetection efficiency $\varepsilon$ given by (1). If the photon is detected, the mean output charge of the SiPM is equal to $q$ in the absence of correlated noise.

When a breakdown avalanche is triggered in a pixel, the instantaneous overvoltage $u$ of that pixel drops to zero and then grows exponentially as

[TABLE]

where $t$ is the delay time from the avalanche triggering and $t_{\rm rec}$ is the recovery time of the SiPM, which ranges from tens of ns to several $\mu$s depending on the device. During the recovery period, both the probability that a photon triggers a new avalanche in the pixel and the mean charge released by this avalanche (if triggered) depend on $u$. The recovery functions of the photodetection efficiency and the gain can be described by

[TABLE]

for $t>t_{0}$, where the time shift $t_{0}$ is defined by $u(t_{0})=U_{0}$. For $t<t_{0}$, $\varepsilon_{\rm rec}(t)=0$ and $q_{\rm rec}(t)=0$. The time shift $t_{0}$ can generally be assumed to be much lower than $t_{\rm rec}$. However, $t_{0}$ may become significant when $U$ approaches $U_{0}$, making $\varepsilon_{\rm rec}(t)$ to grow more slowly than $q_{\rm rec}(t)$.

If several photons hit the same pixel, the time evolution of avalanche triggering and pixel recovery may be very complicated, especially if the light pulse duration is comparable to $t_{\rm rec}$. Besides, each avalanche can induce correlated noise that also contributes to the total output charge. The probabilities of crosstalk and afterpulsing are proportional to the number of charge carriers in the primary avalanche and also depend on the instantaneous overvoltages in the pixels where these secondary processes may take place. Taking into account that secondary avalanches can in turn induce correlated noise, this results in a complex chain process that depends on the temporal and spatial characteristics of crosstalk and afterpulsing (see [16, 1] for details).

Despite of the complexity of the problem, the whole process of charge production for each impinging photon can be modelled in three steps:

i) The photon has a probability $\varepsilon$ to produce an avalanche seed in the pixel that it hits. This event is independent of whether or not more seeds are produced in the same or other pixels.

ii) If the seed is produced, it has a probability $\varepsilon_{\rm rec}(t_{s})/\varepsilon$ to develop into an avalanche with mean charge $q_{\rm rec}(t_{s})$, where $t_{s}$ is the time difference between the occurrence of this seed $s$ and the last avalanche triggered in the same pixel. Notice that $\varepsilon_{\rm rec}(t_{s})/\varepsilon=1$ and $q_{\rm rec}(t_{s})=q$ for $t_{s}\gg t_{\rm rec}$, which corresponds to the case that the seed is produced in a pixel with $u=U$.

iii) If the avalanche is triggered, it may induce crosstalk and afterpulsing. For simplicity, these processes are not considered to produce further seeds, but the total charge from the primary avalanche and all the ensuing secondary ones triggered in the same or other pixels is associated to the same seed. The net charge $q_{s}$ from this seed will have a certain probability density function $f_{s}(t_{s},q_{s})$ with mean $q_{\rm rec}(t_{s})\cdot[1+c_{s}(t_{s})]$, where $c_{s}(t_{s})$ stands for the mean value of the relative contribution from secondary avalanches to $q_{s}$. Notice that, besides the explicit dependence on $t_{s}$, $f_{s}(t_{s},q_{s})$ and $c_{s}(t_{s})$ are generally different for each seed because they depend on the occurrence of both primary and secondary avalanches from all other seeds produced during the light pulse. In case that only one seed is produced, the probability density function and mean value of the net seed charge are simpler and will be referred to respectively as $f(q_{s})$ and $q\cdot(1+c)$.

The number of seeds $m$ produced by a light pulse has a Poisson distribution

[TABLE]

with $\varepsilon\cdot n$ being the mean number of seeds per light pulse. On the other hand, the total output charge $Q=\sum_{s}q_{s}$ depends on the correlations between the avalanche triggering probabilities and the net charges of seeds. These correlations can be understood as the cause of nonlinearity of the SiPM response. Under this scheme, the probability density function of $Q$ can be expressed as

[TABLE]

where $F_{m}(Q)$ is the probability density function of $Q$ in the case that exactly $m$ seeds are produced. Then, the expectation value of $Q$ is

[TABLE]

where

[TABLE]

In the case that only one seed is produced (i.e., $m=1$), $Q$ is equal to the net charge $q_{s}$ of the seed, hence $F_{1}(Q)=f(q_{s})$ and $\langle Q\rangle_{1}=q\cdot(1+c)$. However, even for this simple case, a rigorous theoretical determination of $c$ would require to describe the development of chains of secondary avalanches, including afterpulses that induce crosstalk, and vice versa. In addition, measurements of all the relevant parameters for the calculation of $c$ are arduous and subjected to large uncertainties (see, e.g., [1]). For this reason, $c$ is taken as a free model parameter that is determined from data in the linear response regime of the SiPM, as will be shown later.

Notice that all the complex dynamics of pixel recovery and the interactions between seeds are enclosed in the functions $F_{m}(Q)$. In the following, $F_{m}(Q)$ and $\langle Q\rangle_{m}$ are derived for an arbitrary number of seeds $m$ under some approximations.

II-B Exact solution for a simple case

Before going to the general situation, let us consider first $m=2$ and no correlated noise. The two seeds are generally produced in two different pixels with $u=U$. In this case, both seeds develop into two independent avalanches with charge distributions $f(q_{1})$ and $f(q_{2})$, because they do not interact with each other in the absence of correlated noise. In addition, the mean net charge is $q$ for each seed because $c=0$. Therefore, the probability density function of the total output charge is given by

[TABLE]

where the integration over $q_{s}$ accounts for all the possible combinations of $q_{1}$ and $q_{2}$ adding up to $Q$. The mean output charge is simply $\langle Q\rangle_{2}=2\cdot q$.

However, in the unlikely event that both seeds are produced in the same pixel at times $t^{\prime}$ and $t^{\prime\prime}$, the first one still develops into an avalanche with charge density function $f(q_{1})$ and mean charge $q$, whereas the second one has a probability $\varepsilon_{\rm rec}(t_{2})/\varepsilon$ to trigger an avalanche, provided that $t_{2}=t^{\prime\prime}-t^{\prime}>t_{0}$. If triggered, the charge distribution of this second avalanche has a probability density function $f_{2}(t_{2},q_{2})$ with mean $q_{\rm rec}(t_{2})$ in the absence of correlated noise.

Now, let us consider that the two seeds are produced in the same pixel randomly with a time distribution $p(t)$, which corresponds to the distribution of the arrival times of photons. The probability that both seeds develop into an avalanche is

[TABLE]

where the variables of integration are $t=\min(t^{\prime},t^{\prime\prime})$ and $t_{s}=\lvert t^{\prime}-t^{\prime\prime}\rvert$, and the factor 2 accounts for both possibilities $t^{\prime}<t^{\prime\prime}$ and $t^{\prime}>t^{\prime\prime}$. Here it is assumed $p(t)>0$ for $t>t_{0}$, otherwise $\alpha=0$.

Similarly, it can be defined a time-averaged charge density function for the second seed $s$ (whichever it is) as

[TABLE]

which is properly normalized so that $\int_{0}^{\infty}\phi(q_{s})\cdot{\rm d}q_{s}=1$. The mean of this distribution can be expressed as $\gamma\cdot q/\alpha$, where the dimensionless parameter $\gamma$ is defined by

[TABLE]

From this result, the probability density function of the sum charge for two seeds uniformly distributed over the SiPM area and with time distribution $p(t)$ in the absence of correlated noise is

[TABLE]

The three terms of the right-hand side of this equation arise from the following contributions: i) the two seeds develop into two avalanches in different pixels, ii) both seeds are produced in the same pixel and only the first one develops into an avalanche, iii) both seeds are produced in the same pixel and develop into two avalanches.

Substituting (15) into (10) results in

[TABLE]

II-C Approximations for the general case

It becomes unwieldy to use the above procedure for arbitrary $m$ and including both correlated and uncorrelated noise. Nevertheless, (15) and (16) can be generalized under some simplifications:

i) For a total output charge $Q$, the mean number of fired pixels is $Q/q$ to first-order approximation, which is equivalent to assume that each avalanche is triggered in a different pixel and has charge $q$. In reality, several avalanches (both primary and secondary ones) may be produced in the same pixel, but this is partly compensated by the fact that the greater the number of avalanches in a pixel, the smaller the mean charge per avalanche as a consequence of the reduction of the pixel gain during recovery periods. Under this approximation, the probability that a given seed $s$ is produced in a busy pixel is $(Q-q_{s})/(N\cdot q)$, where $q_{s}$ is the net charge contribution (if any) from this seed and $Q-q_{s}$ is that from all other seeds. Otherwise, the seed is considered to be produced in a free pixel with $u=U$.

ii) If the seed is produced in a free pixel, it is assumed that the avalanche triggering probability of the seed is unity and that the ensuing secondary avalanches are unaffected by the presence of other seeds. Therefore, the probability density distribution and the expectation value of $q_{s}$ are $f(q_{s})$ and $q\cdot(1+c)$, respectively, where $c>0$ now. Notice that, although some seed interactions are ignored under this approximation, the correlated noise contributes to the nonlinearity via a reduction in the number of free pixels.

iii) If the seed is produced in a busy pixel, the probability that an additional avalanche is triggered is assumed to be equal to $\alpha$, as given by (12), that is, only interactions between two seeds are considered. The distribution of the net charge of the last seed (provided that it develops into an avalanche) is still referred to as $\phi(q_{s})$, which is defined by (13) but with the caveat that it should include the charge contribution from correlated noise. Assuming again that secondary avalanches are unaffected by the presence of other seeds, the mean net charge of the seed is $\gamma\cdot q\cdot(1+c)/\alpha$ in this case.

iv) The function $p(t)$ needed to calculate $\alpha$ and $\gamma$ from (12) and (14) is no longer the distribution of the arrival times of photons, because the actual time distribution of avalanches depends on both the temporal characteristics of the correlated noise and the interactions between seeds. Instead, $p(t)$ is approximated by the normalized output pulse signal, since it does include these effects. Note, however, that it may be necessary to perform deconvolution if the output pulse is of similar length to the single-photon response function. In appendices A and B, useful expressions for $\gamma$ are derived for typical pulse shapes, namely rectangular and double exponential pulses.

v) The effect of uncorrelated noise in the SiPM response is approximately included following the procedure described in appendix C. However, the corrections to be made do not affect the derivation and the interpretation of the main features of the model.

From the above considerations it follows this recurrence equation for $m>1$:

[TABLE]

where the three terms of the right-hand side have a similar meaning to those of (15). This equation provides the probability that the last seed of the sequence is produced in either a free or busy pixel and its contribution to the total output charge. Substituting (II-C) into (10) leads to

[TABLE]

II-D Mean response to light pulses

The above result in combination with (9) leads to the following simple expression for the mean output charge of a SiPM for light pulses with a mean number $n$ of photons:

[TABLE]

In Fig. 1, the properties of this expression are illustrated in a plot of the normalized mean charge $\langle Q\rangle/(N\cdot q)$ versus $(\varepsilon\cdot n)/N$ for different values of $c$ and $\gamma$.

For $(\varepsilon\cdot n)/N\ll 1$, the normalized mean charge increases linearly with $(\varepsilon\cdot n)/N$ with a slope $1+c$. In most practical cases, the input light intensity is however characterized by a quantity $X$ proportional to $n$ (e.g., the energy deposited in a calorimeter). Therefore, the factor $1+c$ can be included in a coefficient $k_{X}$ such that $\langle Q\rangle\approx\varepsilon\cdot q\cdot k_{X}\cdot X$ for $k_{X}\cdot X\ll 1$. If $\varepsilon$ and $q$ are known, $k_{X}$ can be determined from data in the linear response region of the SiPM. It should be clarified that, although the integration time window is large enough to contain the full output pulse, the measured value of $\langle Q\rangle$ may not include small contributions from afterpulses with long delay (see, e.g., [1, 17]).

When $(\varepsilon\cdot n)/N\rightarrow\infty$, the normalized mean charge approaches $1/(1-\gamma)$, where $\gamma$ is calculated from (14) and quantifies the average charge losses due to interactions between seeds produced in the same pixel. For a pulse much longer than the recovery time, $\gamma$ approaches one. This makes that (19) is approximately linear even for $(\varepsilon\cdot n)/N>1$, since each pixel may be fired and recovered several times during the pulse. In the opposite case of a very short pulse (i.e., the deconvolution of the output pulse has a length much shorter than $t_{\rm rec}$), one gets $\gamma\approx 0$ and the saturation level is unity. If $c=0$ can be assumed too, then (19) reduces to (2), which corresponds to the limit situation for an instantaneous input light pulse and no correlated noise.

The model can also be adapted to describe the SiPM response to continuous light by assuming a rectangular pulse of length $T\gg t_{\rm rec}$ (see appendix A) and making the substitutions $r=n/T$ and $I=\langle Q\rangle/T$, where $r$ is the incident photon rate and $I$ is the output current. This results in

[TABLE]

where $t_{\rm dead}$ is defined by (26) and it can be interpreted as an effective dead time. When $U$ is large enough so that $t_{0}\ll t_{\rm rec}$ and $\varepsilon\approx\varepsilon_{\rm max}$, then $t_{\rm dead}$ approaches $t_{\rm rec}$. Notice that, in the absence of correlated noise (i.e., $c=0$), the expansion of (20) to first order in $(\varepsilon\cdot r\cdot t_{\rm dead})/N$ is equivalent to (3) for a non-paralyzable detector with dead time $t_{\rm dead}$.

III Experimental verification

The model was validated against experimental data for two SiPMs of the S13360 series from Hamamatsu, namely the type numbers S13360-1325CS and the S13360-1350CS [18]. These detectors were chosen to have the same area but different pixel size, in such a way that nonlinearity is stronger in the 1350 SiPM (pixel size of 50 $\mu$m and $N=667$) than in the 1325 SiPM (pixel size of 25 $\mu$m and $N=2668$). The characteristics of these SiPMs are summarized in table I. The breakdown voltage $V_{\rm br}$, the dark count rate $r_{\rm dc}$ and the mean avalanche charge $q$ are those given in the datasheet of the devices for the recommended operating voltage $U_{\rm op}$ and at a reference temperature of $25^{\circ}$C. The parameters $U_{\rm ch}$, $U_{0}$ and $\varepsilon_{\rm max}$ were obtained by fitting (1) to the data of the photodetection efficiency in the range 2 V$<U<$ 9 V for 450 nm photons and at $25^{\circ}$C also available in the datasheet. Fit residuals were lower than 3%. The values of the recovery time $t_{\rm rec}$ and the overall probability of correlated noise at $U_{\rm op}$ were taken from [1].

A detailed experimental characterization of both SiPMs was previously reported in [1]. After correcting for the amplification factor used in those measurements, the values of $q$ obtained from the pulse charge spectra at photon counting levels are consistent within 10% with the ones given in the datasheet. The correlated noise is dominated by afterpulsing and delayed crosstalk in both SiPMs, which have optical barriers between pixels that drastically reduce prompt crosstalk. The sum of the probabilities of the three components of correlated noise at $U_{\rm op}$ is shown in the table as a reference value, but the charge contribution of the correlated is effectively characterized by the parameter $c$ in this model.

The SiPMs were biased using a Hamamatsu C12332 driver circuit that incorporates a compensation system for the temperature dependence of the SiPM gain. The output signal was registered with a digital oscilloscope Tektronix TDS5032B with math functions, including pulse integration and histogramming. Two types of light sources were used. First, the SiPM was coupled to a scintillation crystal irradiated by $\gamma$ rays. Second, the SiPM was illuminated with continuous light from a LED.

III-A Scintillation light pulses

A LYSO(Ce) scintillation crystal of $3\times 3\times 20$ mm3 with BaSO4 reflector was used. The decay time and light yield of the LYSO(Ce) material are 42 ns and 29 photons/keV, respectively, according to the manufacturer’s data report [19]. Silicone grease was used for the optical coupling to the SiPM. Different radioactive sources (22Na, 60Co, 137Cs and 226Ra) were used, providing a range of $\gamma$-ray energies from 300 to 2100 keV. An integration time window of 400 ns was set on the oscilloscope, which was large enough to contain the full output pulse (see Fig. 2), but still small so that contributions from uncorrelated noise can be neglected ($r_{\rm dc}\cdot T<0.1$).

Output charge spectra for the different radioactive sources were obtained at several overvoltage values for both SiPMs. The optimal values of the output charge resolution for the 137Cs photopeak at 662 keV were 13%(FWHM) for the 1325 SiPM at $U=13$ V and 7%(FWHM) for the 1350 SiPM at $U=5$ V, although the energy resolution was somewhat worse owing to the nonlinear response of the SiPMs, as discussed in [15]. In particular, the resolution degraded significantly at the highest energies for the 1350 SiPM. Nevertheless, seven photopeaks were clearly resolved in the complex spectrum of 226Ra for all $U$ values, except for the photopeaks at 295 keV and 351 keV, which were not well distinguished from the Compton continuum in some cases.

In Fig. 3, the resolved photopeak positions are represented against the $\gamma$-ray energy $E$. The normalized mean output charge $\langle Q\rangle/(N\cdot q)$ is shown to ease the comparison between both SiPMs. The statistical uncertainties of the photopeak positions are smaller than 5%, whereas the systematic uncertainty of $\langle Q\rangle/(N\cdot q)$ was estimated to be 10% due to the uncertainty of $q$. Nonlinearity is more apparent in the 1350 SiPM, because it has less pixels. Similar measurements of the nonlinearity of Hamamatsu SiPMs coupled to LYSO crystals were previously reported in [13, 14].

The number of impinging photons per scintillation pulse is unknown, but it can be assumed to be proportional to $E$. Therefore, the substitution

[TABLE]

was made in (19), where $Y$ is the light yield of the scintillation and $\eta$ is the light collection efficiency, which should be the same for both SiPMs because they have equal area. The coefficient $k_{E}$ is the only fit parameter and it is basically determined by the data in the linear region in which $\langle Q\rangle\approx\varepsilon\cdot q\cdot k_{E}\cdot E$.

For the calculation of the parameter $\gamma$, $p(t)$ was obtained by averaging many output pulses at each $U$ value and normalizing the integral to unity. The recorded pulses were checked to have essentially the same shape regardless of their amplitude. As long as the applied overvoltage was not much higher than the recommended value $U_{\rm op}$, $p(t)$ was well fitted by a double exponential function (27) with $\tau_{1}=12$ ns and $\tau_{2}=45$ ns for the 1325 SiPM, and with $\tau_{1}=15$ ns and $\tau_{2}=60$ ns for the 1350 SiPM. However, the pulse length was found to increase rapidly with overvoltage when it exceeded a certain value, namely 20 V for the 1325 SiPM and 10 V for the 1350 SiPM, as illustrated for the latter SiPM in Fig. 2. When the overvoltage was increased even further, the signal became very noisy and it was not possible to measure the output pulse charge. This pulse widening coincides with a sudden increase in $\langle Q\rangle/(N\cdot q)$ at high overvoltage (see Fig. 3), which is attributed to the fact that the correlated noise is high enough to produce a self-sustained chain of secondary avalanches. This effect was accounted for by calculating $\gamma$ from (28) using the fitted values of $\tau_{1}$ and $\tau_{2}$ at each overvoltage. Results are represented by filled circles connected by solid lines (in blue for the 1350 SiPM and in red for the 1325 SiPM) in the upper plot of Fig. 4. Instead of $\gamma$, the saturation level $1/(1-\gamma)$ is shown. Notice that $1/(1-\gamma)$ increases smoothly with overvoltage until the correlated noise starts to become significant.

Solid lines in Fig. 3 represent the results from (19) when using these calculated values of $\gamma$ and fitting $k_{E}$ to data. The model describes correctly the SiPM response and its dependence on overvoltage as long as nonlinearity is not very strong. For the 1350 SiPM, the predicted saturation level is about 4% lower than the measured one, but discrepancies are within the experimental systematic uncertainties.

The fitted values of $k_{E}$ are represented by filled circles connected by solid lines in the lower plot of Fig. 4. Statistical uncertainties are smaller than the data point size. At low overvoltage where the correlated noise is very small, $k_{E}$ is approximately 3.5 keV*-1* for both SiPMs. Taking $Y=29$ photons/keV and $c=0$ in (21) results in a light collection efficiency of $\eta=12\%$, which is reasonable considering that the SiPM area covers only 18% of the crystal surface.

For the 1325 SiPM, $k_{E}\approx 9$ keV*-1* at $U=21$ V, which translates into $c\approx 1.5$, that is, the charge contribution of correlated noise is as high as 60% of the output charge, suggesting that chains of secondary avalanches are produced efficiently at such a high overvoltage. The overall probability of correlated noise is 0.5 at $U=21$ V for this SiPM according to data from [1]. Therefore, if no more than one secondary avalanche could be induced by each primary avalanche, then $c$ should be at most 0.5 at this overvoltage.

For the 1350 SiPM, the fitted values of $k_{E}$ do not show the expected monotonous increase with overvoltage. This may be explained by the fact that $\varepsilon$ was obtained by fitting the approximated expression (1) to data for 450 nm photons, not for the specific emission spectrum of the LYSO crystal.

Tests were also made where $\gamma$ was calculated using the deconvolution of the output pulse signal from the single-photon response function, which can be approximated by an exponential function with decay time of 9 ns for both SiPMs. The main effect of the deconvolution is to make the rising edge of the output pulse faster (i.e., $\tau_{1}$ decreases down to a few nanoseconds). As a result, the saturation level is reduced by around 9% and 6% for the 1325 SiPM and 1350 SiPM, respectively, whereas $k_{E}$ is basically unchanged.

For consistency, the model was also fitted to the data taking both $\gamma$ and $k_{E}$ as a free parameters. The best fit is represented by dotted lines in Fig. 3, and the fitted values of $1/(1-\gamma)$ and $k_{E}$ are shown as open circles connected by dotted lines in Fig. 4. For the 1325 SiPM, this fit gives almost identical results to the previous one, except for the fact that the uncertainties of the values of $1/(1-\gamma)$ are 20% or larger, since the nonlinearity is small in this energy range. For the 1350 SiPM, the fitted values of $1/(1-\gamma)$ are slightly higher than the calculated ones, resulting in a better agreement of the model with experimental data.

III-B Continuous light

To characterize the SiPM response for continuous light, a green LED (Kingbright L-53GD) was placed at a few centimeters from the SiPM. The mean output current was measured as a function of the LED forward current $I_{\rm LED}$ at different overvoltage values for both SiPMs. Measurements were also performed placing the LED farther away from the SiPM to ensure that it works in the linear region, checking that the light output of the LED was proportional to $I_{\rm LED}$. The uncorrelated noise was found to be negligible in all the measurements.

Results are shown in Fig. 5, where the normalized output current $(I\cdot t_{\rm rec})/(N\cdot q)$ is shown to ease the comparison between both SiPMs. The correlated systematic uncertainty was estimated to be 12% due to uncertainties of $q$ and $t_{\rm rec}$. Notice that the response curves are of a very similar shape to those shown in Fig. 3 for scintillation light pulses.

The substitution $(1+c)\cdot r=k_{\rm LED}\cdot I_{\rm LED}$ was made in (20), where $k_{\rm LED}$ is a fit parameter that has a similar meaning as the parameter $k_{E}$ used for scintillation pulses. The saturation level of $(I\cdot t_{\rm rec})/(N\cdot q)$ is $t_{\rm rec}/(2\cdot t_{\rm dead})$, where the parameter $t_{\rm dead}$ was calculated from (26). The best fit of (20) to the data is represented by solid lines in Fig. 5. The model describes approximately the experimental data, although it somewhat underestimates the saturation level for the 1350 SiPM.

The calculated values of $t_{\rm rec}/(2\cdot t_{\rm dead})$ and the fitted values of $k_{\rm LED}$ for both SiPMs are represented by filled circles connected by solid lines in Fig. 6. As mentioned above, $t_{\rm dead}$ approaches $t_{\rm rec}$ when $U$ is large, so that the upper limit of $t_{\rm rec}/(2\cdot t_{\rm dead})$ is $1/2$. Notice that $k_{\rm LED}$ depends on overvoltage in a very similar way to $k_{E}$ for both SiPMs, which supports that chains of secondary avalanches have an important role in the SiPM response.

The model was also adjusted to data by fitting both $t_{\rm dead}$ and $k_{\rm LED}$. The resulting fit is shown by dotted lines in Fig. 5 and the fitted values of $t_{\rm rec}/(2\cdot t_{\rm dead})$ and $k_{\rm LED}$ are represented by open circles connected by dotted lines in Fig. 6. As in the case of the scintillation data, $t_{\rm dead}$ is only poorly determined for the 1325 SiPM because its response is approximately linear in this range. For the 1350 SiPM, the fitted values of $t_{\rm rec}/(2\cdot t_{\rm dead})$ for $U>1$ V are systematically larger than the calculated ones by about 20%. Nevertheless, both the fitted and the calculated values show the same overvoltage dependence, except for $U=1$ V where the SiPM response is nearly linear.

IV Conclusions

A model of the nonlinear response of SiPMs that includes losses of both the avalanche triggering efficiency and the gain of pixels during recovery periods as well as the effects of correlated and uncorrelated noise has been developed. The simple expression (19) has been derived, which describes the mean output charge of a SiPM for light pulses of arbitrary shape as long as nonlinear effects and correlated noise are moderate. In addition, the modified expression (20) has been obtained for continuous light. The model reduces to the well-known expressions (2) and (3) for instantaneous light pulses and for continuous light, respectively, when both correlated and uncorrelated noise are ignored.

Two principal parameters are introduced in the model. First, the parameter $\gamma$ (or $t_{\rm dead}$ in the case of continuous light) gives the saturation level of the SiPM response and it is calculated by a numerical integration knowing the output pulse shape $p(t)$, the pixel recovery time $t_{\rm rec}$ and the photodetection efficiency $\varepsilon$. Simplified expressions for $\gamma$ have also been derived for both rectangular pulses and double exponential pulses. Second, the parameter $c$ represents the relative charge contribution of the correlated noise to the output charge. This parameter is taken as a free parameter that is determined by fitting the model to data in the linear response region of the SiPM. Taking both parameters as a free parameters is nearly equivalent to some previous fitting models used to describe the nonlinear response of SiPM [12, 13, 14], but the present model provides an adequate interpretation of these parameters, including their dependence on both the overvoltage and the output pulse shape.

The model was validated against experimental data of two different SiPMs at moderate nonlinearity using scintillation light pulses from a LYSO crystal. The output pulses were found to be significantly widened when the correlated noise became important, which also caused an increase in the saturation limit of the SiPM response. This effect was accounted for by using the output pulse shape measured at each supplied voltage for the calculation of the parameter $\gamma$. The model describes the data within the experimental uncertainties.

The model also describes adequately the SiPM response for continuous light, although the saturation level was underestimated by about 20% for one of the SiPMs in this case. More tests are necessary to understand this discrepancy. However, the model describes well the data when both $c$ and $t_{\rm dead}$ are taken as free parameters in the fit. In addition, the overvoltage dependence of the fit parameters is consistent with that obtained for scintillation light pulses.

An adequate parameterization of the chains of secondary avalanches will be useful to provide a theoretical estimation of both the parameter $c$ and the distortion of the output pulse shape. This could be especially relevant for applications where the output pulse height (or a similar parameter) is used instead of the integrated pulse charge.

Appendix A Rectangular pulses

A normalized rectangular pulse can be expressed as

[TABLE]

Substituting this into (14) results in

[TABLE]

For very short pulses such that $T\ll t_{\rm rec}$ and assuming also $t_{0}\ll t_{\rm rec}$, the first-order Taylor expansions of $\varepsilon_{\rm rec}(t_{s})$ and $q_{\rm rec}(t_{s})$ about $t_{s}=0$ can be used in (23), resulting in

[TABLE]

On the other hand, for very long pulses, i.e., $t_{\rm rec}\ll T$, the expansion of (23) to first order in $t_{\rm rec}/T$ and $t_{0}/T$ gives

[TABLE]

where $t_{\rm dead}$ is defined by

[TABLE]

Appendix B Double exponential pulses

In many cases, the output pulse can be fitted by a double exponential function

[TABLE]

with $t\geq 0$ and $\tau_{1}<\tau_{2}$. The rise time and the fall time of the pulse are approximately given by $\tau_{1}$ and $\tau_{2}$, respectively, when $\tau_{1}\ll\tau_{2}$.

Substituting (27) into (14) results in

[TABLE]

where

[TABLE]

Here it has been assumed that the integration time window $T$ is large enough so that $p(t)\approx 0$ for $t>T$.

For $\tau_{i}\ll t_{\rm rec}$, (29) reduces to

[TABLE]

On the other hand, for $t_{\rm rec}\ll\tau_{i}$, (29) reduces to

[TABLE]

where it has been used that $1-[\varepsilon_{\rm rec}(t_{s})\cdot q_{\rm rec}(t_{s})]/(\varepsilon\cdot q)$ goes to zero much faster than $\exp\left(-t_{s}/\tau_{i}\right)$. Therefore, for a very long pulse such that $t_{\rm rec}\ll\tau_{1}<\tau_{2}$, one gets

[TABLE]

The particular case of $\tau_{1}=0$ in (27) corresponds to a single exponential function with decay time $\tau=\tau_{2}$, for which $\gamma$ reduces to $\gamma_{2}$, which is given by (29).

Appendix C Corrections for uncorrelated noise

The uncorrelated noise can be readily accounted for by supposing a low rate $r_{\rm bg}$ of background photons in addition to those from the light source. This rate $r_{\rm bg}$ can be characterized by measuring the background output charge $\langle Q\rangle_{\rm bg}$ in a given integration time window $T$. Assuming linearity, it is obtained that

[TABLE]

where $\varepsilon\cdot r_{\rm bg}$ can be regarded as a continuous rate of background seeds, each one contributing with a mean net charge $q\cdot(1+c)$. It should be clarified that the dark count rate $r_{\rm dc}$, usually given in datasheet specifications, is expected to be slightly larger than $\varepsilon\cdot r_{\rm bg}$, because $r_{\rm dc}$ includes afterpulses surpassing the detection threshold.

The effect of the uncorrelated noise is twofold. In the first place, the effective number of photons $n_{\rm eff}=n+r_{\rm bg}\cdot T$ should be substituted for $n$ in (19). In the second place, the effective time distribution $p_{\rm eff}(t)$ including uncorrelated noise is

[TABLE]

where $\chi=(r_{\rm bg}\cdot T)/n_{\rm eff}$ so that $\int_{0}^{T}p_{\rm eff}(t)\cdot{\rm d}t=1$. Assuming $\chi\ll 1$, the effective parameter $\gamma_{\rm eff}$ is calculated from

[TABLE]

The first term of the right-hand part of this equation includes only the contribution of seeds produced by the light pulse, where $\gamma$ is calculated from (14) assuming that $T$ is large enough to contain the full pulse, and the second term arises from interactions of seeds produced by the light pulse with background seeds. To first order, interactions between background seeds can be neglected.

Corrections for uncorrelated noise are only expected to be significant for long light pulses requiring $T\gtrsim 1$ $\mu$s, because $r_{\rm dc}$ is typically lower than 1000 kcps [18].

Acknowledgment

This work was supported by Spanish MINECO under the contract FPA2017-82729-C6-3-R.