Skewness-based characterization of silicon photomultipliers

Characterization of SiPMs is an objective of high importance in almost any research, development, and application of these unique photon detectors. Two decades of characterization method developments resulted in a comprehensive and elegant methodology based on a precisely resolved number of fired SiPM cells or the SiPM spectrum. Spectrum fitting procedures are very sensitive to a proper selection of many initial guess values including parameters of minor importance and happen to be unreliable. Moreover, conventional methods are useless for degraded or unresolved SiPM spectra. This challenge was anticipated for SiPM-based detectors of high-luminosity particle calorimeters where a severe radiation degradation of the SiPMs would co-exist with their acceptable performance. To address it, one can measure a mean and variance of SiPM charge responses relying on a priori known Excess Noise Factor of the SiPM but its stability during SiPM degradation is uncertain. This study proposes and evaluates a new approach for the SiPM characterization based on the first three statistical moments (mean, variance, and skewness) of its charge response including skewness. It assumes the Generalized Poisson distribution of the number of fired cells and provides simple closed-form analytical expressions for the parameter estimation. It allows determining a gain, a number of detected photons, and a probability of correlated events directly from raw data. The skewness-based characterization is anticipated to be especially useful for mass testing and continuous monitoring of SiPMs in large-scale experiments due to its simplicity and robustness. Initial evaluations of the method show promising results.


Introduction
One of the main directions in the development of photodetectors has been dedicated to improving the low light level a e-mail: vinogradovsl@lebedev.ru (corresponding author) sensitivity down to single photons. An internal multiplication of charge carriers in the photodetector has been recognized as a straightforward measure to achieve this goal by amplifying a single photon response above electronic noise.
Research and development of vacuum PMTs and solidstate APDs revealed a gain dispersion with a corresponding Excess Noise Factor (ENF) of the gain as an inherent drawback of these devices due to the stochastic nature of the charge multiplication process [1][2][3]. Theoretical studies of photon detection statistics of these photomultipliers resulted in a basic understanding of how their performance, first of all, in terms of signal-to-noise ratio (SNR) and amplitude resolution, depends on the main parameters -quantum efficiency, gain, excess noise factor, dark and electronic noise [4][5][6].
On the other hand, a development of the photomultiplier characterization methods became supported by analytical models of the photon detection statistics, mostly related to the first and second statistical moments, i.e. mean and variance of the photoresponse.
Particularly, such an approach has been proposed and initially evaluated for PMTs [7], and then utilized as the "ENF method" for continuous monitoring and calibration of the gain of the MAGIC Cherenkov telescope PMTs [8]. The method is based on the mean and variance of measured charge in a response to a short LED pulse with known ENF values for each PMT. In such a case the PMT ENF is precisely measured in a laboratory by the mean and variance of a single electron response (SER) charge.
A problem to develop a high gain APD with a good SNR is caused by the strong dependence of the ENF on the mean gain of avalanche multiplication processes in semiconductors. It has been resolved in Solid State or Silicon Photomultipliers [9,10]. Indeed, even for the earlier SiPM generations, the ENF of the gain was as low as 1.01 [11] and 1.003 [12]. Nevertheless, the performance of the SiPM happens to be degraded by stochastic processes of crosstalk and afterpulsing due to its multipixel architecture and high gain. Thus, instead of the gain dispersion in the linear-mode APDs, the crosstalk and afterpulsing in SiPMs happen to play a similar role producing specific ENF [13][14][15].
Characterization of the SiPM parameters has been in the focus of many studies for two decades and became a wellestablished and elegant methodology mostly based on a precisely resolved number of fired SiPM cells in the so-called photo-electron and dark-electron spectrum (see, for example, open materials of the International Conference on the Advancement of Silicon Photomultipliers (ICASiPM) [16] and the comprehensive overview by R. Klanner [17]). At the same time, it is still subjected to intense development and discussions, especially concerning the characterization of the radiation-degraded SiPMs [18,19]. Irradiation gradually increases the dark count rate of the SiPMs resulting in degradation of their performance and making unresolvable their single photo-electron spectra.
A challenge of characterization and monitoring of SiPMs with blurred and even unresolved single photo-electron spectra is well recognized in ongoing projects on novel particle calorimeters, e.g. CERN CMS high granularity calorimeter Phase 2 upgrade for the high luminosity LHC, where, as anticipated, a severe radiation degradation of SiPMs with a GHz-level dark count rate (DCR) would co-exist with their acceptable photoresponse performance. The most important SiPM parameters to be measured and monitored for the calorimetry -a mean avalanche multiplication gain g, a mean number of primary fired SiPM cells or detected photons (also so-called photoelectrons) η, and probability p to observe correlated events (optical crosstalk and afterpulsing) in the detection time gate.
The resolved spectrum allows to determine the mean gain g by the neighboring peak distances, the detected photons η by the area of the zero-event peak assuming Poisson statistics of the photons, and the probability p by deviation of the spectrum from Poisson statistics for the non-zero-event peaks. Typically, the afterpulsing contribution in the gated charge spectra is negligible with respect to the crosstalk, and we can attribute p to the probability of crosstalk. To extract the main parameters of the SiPM with the unresolved spectrum, just a few approaches of rather limited performance have been considered till now. First, one can determine g and η by the ENF method [8] with a priori known ENF. The ENF of modern SiPMs with very low excess noise of an avalanche multiplication is approximately equal to ENF of the correlated events with a major contribution from the crosstalk. Assuming the crosstalk ENF stability during radiation-induced degradation of the SiPM, we can rely on measured ENF values before irradiation (while the SiPM spectrum is resolved). The ENF method was considered as a possible way to monitor the SiPM gain during long-term operations of the CERN CMS high granularity calorimeter [20,21]. Then it has been successfully evaluated as the regular monitoring procedure [22,23].
The second approach is based on the Generalized Poisson (GP) distribution model of the SiPM spectra [24]. In general, it allows fitting all parameters at once, however, while the guess values for the fit could be reasonably determined from resolved spectra or just known a priori. GP model shows good agreement with various experiments, for example, burnin tests of 3000 devices at the CERN CMS [25] and SiPM studies in a framework of the Cherenkov Telescope Array collaboration [26,27].
Moreover, the GP model has been advanced by the University of Hamburg group to account for complex contributions from fully and partially integrated dark single electron response (SER) pulses of a single-exponential shape, delayed crosstalk, and afterpulsing [28]. The advanced GP model describes the non-Gaussian shape of the spectrum peaks and demonstrates an excellent agreement with experiments even for unresolved SiPM spectra of about 50 mean photoelectrons.
This study presents a new approach to the characterization of SiPMs based on the method of moments including the third central statistical moment or, equally, the skewness of the measured charge. This approach relies on the assumptions on the Generalized Poisson distribution of the SiPM charge response and negligible ENF of the avalanche gain.
This approach provides three main parameters of the SiPM -the gain, the mean number of primary events, and the probability of correlated events in a charge integration time gate (presumable prompt optical crosstalk) -in a very simple and easily automated procedure without any need to fit the SiPM spectra.
Moreover, in general, it seems to be applicable to the unresolved spectra using just raw data. However, there is a severe limitation on the DCR noise contribution with respect to the mean number of detected photons, namely, the mean number of dark events in a time gate should not exceed 0.1η to keep the correctness of the method according to zero-DCR simulations as presented in Sect. 4.2.

Method of moments for SiPM charge response
Total SiPM charge response Q on to a light pulse is formed as a sum of contributions from photo-response events Q photo , dark events Q dark , and a noise charge of front-end electronics Q el collected during the charge integration time gate t gate . All these charges represent random variables.
To characterize photodetection of the SiPM, typical measurements are carried out with light turned on to get Q on and with the light turned off to get Q of f . Throughout the paper, all charges are given in units of electrons, not in Coulombs. Let us turn to statistical properties of the charges: Where N is a random number of the SiPM avalanche events or fired pixels, G i is a random gain value of each SER event arrived in a random time t i , and h(t) is a filter function in the theory of filtered point processes which represents the SER shape normalized to a unit area.
Let us distinguish the variables under our discussion on an example of the gain. We designate this random variable by a capital letter G and its expectation (mean value, non-random parameter) by a small letter g so that E[G] = g. However, we only observe realizations of the random variable Gparticular values G i , and the set of G i allows us to estimate g by an estimatorĝ using some estimation procedure, e.g., averaging.
We assume, as it used to be in most experiments on SiPM charge spectra, that the integration time t gate is set to be longer than the duration of the light pulse and SER decay time, and the time gate position is properly adjusted to fully integrate photoresponse pulses. In such a case the charge Q photo is equal to a sum of N gains G i . The random number N is composed by a random number of primary events (or primaries for short) N prim and chains of correlated events K : Where K j is a random number of correlated events per single primary event j. In the Eq. (4) we assume that all correlated events originated from the photoresponse primaries and contributing to Q photo are fully integrated. Prompt optical crosstalk presumable corresponds to this case, and, throughout the paper, we equally use the term crosstalk for the correlated events with such properties.
Assuming Poisson distribution of the primaries N prim ∼ Poisson(η), where η is a mean, statistical moments of Q photo can be expressed as Where μ 1 is the first raw moment of the charge, μ 2 is the second central moment, σ is the standard deviation, g = E[G] is a mean of random variable G, κ = E[K ] is a mean number of crosstalk events per single primary, ENF of gain F G and ENF of crosstalk F K are defined as for any random variable X : Equations (5)- (7) have been used as grounds for the ENF method for PMTs (see Sect. 1, [8]) with an assumption of negligible correlated effects (i.e. κ = 0 and F K = 1) as well as for the two-moment approach for SiPMs [20]- [25], [29] with an assumption of negligible gain fluctuations (i.e. F G = 1). Indeed, for modern SiPMs σ G ∼ 0.03g and F G ∼ 1.001). Therefore, the gain G is assumed to be a non-random variable: Next assumption: the crosstalk is a Poisson branching process with a total number of events of Borel distribution: 1 + K ∼ Borel(λ) where λ is a mean number of events in a single generation of the branching process initiated by a single preceding event. It relates to the probability p to get at least one crosstalk per single primary as follows: Therefore, the Poisson primaries with the Borel crosstalk events result in GP distribution of the total number of fired cells: N ∼ G P(η, λ). And specific expressions for the mean and ENF of the crosstalk of the GP model (see Table 2 in [24]) provide the next simplifications: Therefore, the parameters λ, p, κ, and F K could equally be used and easily substituted. Let us note the expression (10) is often used to define "Excess Charge Factor" of the correlated effects misinterpreting its sole origin in an extra number of events while the charge (namely, gain) of the events could be even lower than the full SER gain due to incomplete recovery of the SiPM pixels in short-term afterpulsing events.
Finally, the skewness method expressions for the population moments of the GP random variable Q photo = N · G (where the gain G = E[G] = g) are as follows [30]: Where μ 3 is a third central moments and γ is a skewness of the Q photo (t gate ) expressed by definition as Let us discuss relations between Q photo and the measured charges Q on and Q of f as far as the natural way to determine the Q photo moments is a subtraction: The Eq. (15) is valid if one of the options (A or B or C) is fulfilled: (A) Q on and Q of f are independent random variables; (B) Q el = 0 (practically negligible) and Q dark is a GP random variable; (C) Q of f = 0 (practically negligible).
The case A presumes independence of results of two independent experiments with and without light, and it is a property of any stochastic system with a linear response (i.e. output Z on input X and Y : Z (X + Y ) = Z (X ) + Z (Y )). The SiPM has some limited linear range, and non-linearity and saturation of the response destroy the independence, say, both Q on and Q of f approach to the same saturation level defined by a number of SiPM pixels or occupancy of the pixels as discussed in [15][16][17]. Moreover, the GP model is applicable only for the linear SiPM response because in general, the number of fired pixels obeys binomial distribution.
The case B could happen due to (typically) low RMS of the electronic noise with respect to the signal (μ 2 (Q el ) μ 2 (Q dark ) and μ 2 (Q photo ) and (typically) normal (Gaussian) distribution of the random variable Q el of zero skewness. However, Q dark in case B should be formed by fully integrated SER charges (Poisson number of primaries) with a negligible amount of incompletely integrated SER with arrival times before t = 0 and t = t gate . It is possible in a specific case of low DCR and large t gate , and of course, with the linear SiPM response.
The case C is a trivial one; it seems to be a possible option in low-temperature SiPM experiments.
Throughout the paper, we assume the validity of the Eq. (15) with the general case A, thus, we get rid of concerns about probability distributions of the electronic noise Q el and the dark charge Q dark with some possible portion of incomplete integrated SERs. So, the Q of f moments are assumed to be completely subtracted from the Q on moments; they do not affect the skewness-based method and, therefore, will not be considered in detail.
The method of moments, introduced by P. Chebyshev and developed by K. Pearson, appears to be rather simple especially if the above moment equations (11)-(13) could be solved for the parameters η, λ, and g: where θ = (η, λ, g). Indeed, in our case the parameters can be expressed in a closed form: In the two-moment approach, the ENF has to be known a priori. Equation (16) eliminates this requirement because F K (Eq. 9) follows directly from the mean, sigma, and skewness of the measured data. Moreover, considering that the term or, equally, μγ σ is equal to 1 for a pure Poisson distribution, it can be interpreted as an "Excess Skewness Factor". Indeed, for the pure Poisson random variables λ = 0, g = 1 and μ 1 = μ 2 = μ 3 = η as can be seen in the Eqs. (11)-(13), therefore μ 1 μ 3 Finally, since the sample moments of the measured charge m i are estimators of the population moments μ i , the estimators of the distribution parameters of Q photo areθ = (η,λ,ĝ) defined by the Eqs. (16)- (18) where μ i are replaced by m i .

Statistical errors of estimators
The main characteristics of statistical errors of the parameter estimatorθ are its bias B(θ), variance V ar(θ), and mean squared error M SE(θ): Here S E(θ) is the standard error of the parameter estimator; B rel (θ) and S E rel (θ) are the corresponding relative errors. The Eqs. (19)- (20) should be expressed by the GP distribution parameters. It can be done in two steps. The first step is to approximate the statistical moments of a function of random variablesθ = f (m) using Taylor expansion series. Let us limit the expansion up to the first significant term which provides the most valuable and sufficient information in typical practical cases (V ar(θ) -first order, B(θ) -second order): Second step is related to the so-called "moments of moments" problem: calculation of the population moments of random variables -sample moments m i , namely Cov(m j , m k ) for all combinations of j, k = 1, 2, 3. The results of our interest have already been derived using the power sums technique and could be explored in a framework of Wolfram Mathematica software as well as in the book "Mathematical Statistics with Mathematica" [31].
Calculating derivatives of the function f (m) and applying known expressions for Cov(m j , m k ) from [31], we can derive the error expressions (24) and (25) in an analytical form. For example, Cov(m 1 , m 2 ) = μ 3 n , where n is a sample size. However, the final expressions are rather long, and their convenient approximation is presented in Appendix A.

Monte Carlo simulations of statistical errors
As expected, the accuracy of the skewness method should be rather sensitive to the sufficiency of statistics required to extract information related to the skewness, which is mainly associated with the upper tail of the charge distribution.
Monte Carlo (MC) simulation is a convenient tool to analyze statistical behavior and characteristics of random variables. It allows checking the statistics sufficiency for large sample sizes substituting long-lasting experiments affected by various instabilities. In this study, the MC simulations have been performed using custom GP and Borel random number generators based on an inversion method algorithm as discussed in ch. 16 "Computer Generation of Lagrangian Variables" [30].
For each given sample size n of the random variable Q photo (t gate ), we have repeated our MC simulations (MC experiments) nn = 100 times and extracted 100 values of estimated parameters (η, λ, g) to calculate their biases and standard errors within a 10% uncertainty on our check. Regretfully, this limitation has been applied as the MC simulations are rather time-consuming.
Finally, we crosschecked the MC simulations and the analytical results on V ar(θ) and B(θ). A good agreement between the MC and analytical results as presented in Sect. 4.1 allows us to rely on the analytical modeling in the whole parameter space.

Experimental technique
The block diagram of some typical setups for the SiPM characterization experiments including charge response measurements with short laser pulses as well as in the dark conditions is shown in Fig. 1.
In this particular case, a PicoQuant PDL 800-B/LDH-P-C-440 laser with a wavelength of 442 nm and a pulse duration of 60 ps was used as a light pulse source. Uniform exposure of the entire photosensitive area of the sample was provided by a diffuse glass plate. The SiPM output HF signals were acquired using a Mini-Circuits ZX60-4016ES+ preamplifier with a gain of 10 and a bandwidth of 4 GHz, and a Tektronix TDS7104b digital oscilloscope with a sampling rate of 10 GHz. A Keithley SMU-236 precision source-meter was used to supply voltage and to monitor the SiPM current. This setup has been developed by the author and his colleagues at the Lebedev Physical Institute for SiPM studies.
Typical SiPM response acquisition by the digital oscilloscope is shown in Fig. 2. The Hamamatsu MPPC S10362-11-050C was used for the initial evaluation of the skewnessbased method. This MPPC series is well known for the SiPM community, and its specification is described later in the corresponding experiments.

Processing of experimental data
The characterization of SiPMs is mainly based on the SiPM spectrum with resolved events. Determination of the mean number of detected photons by the probability of the zeropeak events Pr(0) assuming Poisson statistics of the photons has been well established for PMTs (e.g. [1]). The PMT gain is used to be determined by the single-event peak position in single photoelectron spectra. The same approach has been Fig. 2 The acquisition of the SiPM response on a short few-photon laser pulse using Hamamatsu MPPC S10362-11-050C applied to the first resolvable SiPM spectra to determine the mean number of primaries η and the gain g [10,32]. Probability of the crosstalk p were determined by a dependence of DCR on a discriminator threshold (e.g. [33,34]): The Eq. (26) defines a probability to get one or more correlated events per single primary described by a complementary cumulative distribution function (CCDF) of the prompt events. However, this standard approach does not account for occasional coincidences of two and more simultaneous independent dark events as proposed in [35], but for typical experimental conditions with a sub-MHz DCR level, such a correction is negligible.
Extraction of the probability of crosstalk from the light SiPM spectrum of photon detection by deviation of the single-event peak Pr(1) from Poisson statistics has also been demonstrated (e.g. [14]) using the expression: Here, Pr(0) and Pr(1) are the probabilities of observing zero and one event in a single detection calculated by the relative areas of the corresponding SiPM spectrum peaks.
Therefore, let us define the methods pointed out above as the standard ones, where the Eqs. (26) and (27) are used in the light and dark measurements correspondingly, so that the skewness estimates under study can be referenced to them. More details and comprehensive considerations of characterization methods as well as advanced techniques could be found in [17] and materials from the International Conference on the Advancement of Silicon Photomultipliers [16].

Simulation results
First, the absolute standard errors S E(θ) and the biases B(θ) of the parameter estimates, as well as their relative forms S E rel (θ) and B rel (θ) are calculated for η = 10 and η = 1 in their dependence on the sample size n. Let us note that the S E(g) and B(g) are linear functions of the gain because G is assumed to be a non-random variable which value is equal to the mean g: G = g = const (see Eqs. 11-13). The relative biases of all parameters (Figs. 3, 4, 5) are a few times lower than the corresponding relative standard errors (Figs. 6, 7, 8) for the sample size of 100 becoming about an order of magnitude lower for the sample size of 1K and more. This fast decay corresponds to the fact that B(θ) ∼ 1/n and V ar(θ) ∼ 1/n (see Eqs. 25, 24 of the same structure with respect to Cov(m j , m k )) and, therefore, S E(θ) ∼ 1/ √ n. Second, the standard errors have been analyzed in their dependence on a mean number of primaries η. This parameter is selected because it can be varied by the intensity of the light pulse independently of gain and crosstalk. On the contrary, gain or crosstalk can only be varied by changing overvoltage, which affects all parameters simultaneously, complicating the analysis of the model. The question is to find a value or a range of optimal light pulse intensities with the minimal errors of the skewness method.
Results of the modeling are presented in Figs. 9, 10, 11. Some dispersion of MC data points over analytical ones cor- responds to a 10% uncertainty due to the rather limited number of MC simulation experiments nn = 100.
Overall, the optimal light intensities correspond to a range from about one to four mean detected photons, and could be extended to ten primaries and more without a dramatic increase of the errors.

Monte Carlo simulations with dark noise
The third analysis has been focused on the capability of the skewness method to characterize SiPM parameters if the SiPM spectra are degraded and become unresolvable due to increased DCR. Let us note that the dark noise has not been accounted for in the previous simulations.
The dark charge data has been generated by integration of SER pulses h(t − t i ) of bi-exponential temporal shape (rise time of 1 ns, fall time of 20 ns) in 100 ns time gate. The arrival times t i have been uniformly distributed in the time gate as well as in the earlier time interval of 500 ns to form a contribution to the collected charge from the SER tails of the preceding events. The uniform distribution of t i is equivalent to the Poisson distribution of a number of primary dark events in any time interval, and the dark-initiated prompt crosstalk events have been simulated by Borel distribution of the total events originated from the single primary with the same arrival time.
Thus, the resulting Poisson sum of Borel events has a Generalized Poisson distribution. Any event generates a random charge of Gaussian distribution with a mean g and a sigma of 0.03 g (F G ∼ 1.001). Electronic noise has been accounted for as a contribution of additive white Gaussian noise with a zero mean and a sigma of 0.1 g. The dark and electronic noise charges are presumed to be independent of the photo-event charges (case A of validity of Eq. 15).
The light charge data has been formed as a sum of the dark charges described above and completely integrated photoresponse charges with corresponding Generalized Poisson distribution of the events. Delayed crosstalk and afterpulsing effects have not been accounted for in the simulations. The integrated dark (η = 0) and light (η = 3) charge datasets with crosstalk λ = 0.1 are presented by the histograms in Figs. 12 and 13. The values of DC R = 10 MHz, 30 MHz, and 100 MHz are chosen to be equal to 1, 3 and 10 of the mean dark primaries at the time gate of 100 ns, and the corresponding spectra could be specified as resolvable, hardly resolvable and unresolvable. The histograms are shown just to visualize the spectra, while the parameter extraction has been based on the raw data with the subtraction of the dark moments from the light moments accordingly Eq. (15).
The relative biases of the parameter estimates from the MC simulations (Figs. 14, 15, 16) are randomly distributed (due to a limited number of MC experiments nn = 100) around the analytical model values without the dark noise  (Figs. 17, 18, 19) which are more sensitive indicators of the deviations, so we see that they start to increase for DC R > 3 MHz (0.3 dark primaries).
To clarify the possibility of decreasing the statistical errors at high DCR, the MC experiments with 30 detected photons have been simulated. The high level of η shifts the starting point of the increased bias deviations to about 100 MHz (Figs. 14, 15, 16), but, as predicted by the analytical model and shown in Figs. 9, 10, 11, it increases the low DCR level of the deviations.
However, large-signal MC experiments for η = 30 provide lower statistical errors than those for η = 3 if DC R > 30 MHz. It means, the optimal conditions and parameters of the skewness method should be adapted to the DCR level, say, roughly, to keep DC R · t gate ≤ 0.1η. An attempt has been made to introduce some DCR-related adjustment to the zero-DCR analytical model considering V ar(θ) (Eq. 25) as follows: In particular, the simple result can be derived for V ar(m 1 ): where the covariance is equal to zero, since the light and dark measurements are independent. Thus, the corrected expression for V ar(m 1 ) includes doubled dark noise contribution: V ar(m 1 ) = V ar(m 1 ( photo)) + 2V ar(m 1 (dark)) Fig. 20 The spectrum of the MPPC S10362-11-050C HPK and its GP and pure Poisson model fits. The crosstalk probability determined by the standard method is 8.3%, by the GP fit -9.0% Where a correction coefficient F DC R coincides with the ENF of DCR introduced in [15]: It adjusts only one contribution in the Eq. (25) from V ar(m 1 ), and full consideration of DCR-related contributions to statistical errors appears to be beyond the scope of this introductory article.
Therefore, the zero-DCR expressions for V ar(θ) have been approximated by F DC R ·V ar(θ), and the corresponding relative standard errors are presented in Figs. 17, 18, 19 by the dotted lines. This simplified correction explains the slow initial increase of the MC plots up to DC R · t gate ≤ 0.3η but appears to be insufficient for the higher DC R values. Probably, the higher-order moments with a superlinear dependence on η (see Eq. A.2, Appendix A) and their variances and covariances are responsible for that.
The questions of how the skewness method is affected by DCR and how much its applicability can be extended are very important to be addressed in future studies on this subject.

Experimental results
Systematic studies of the SiPM spectra on the probability distribution of the SiPM responses were related to the first model of the correlated events, the geometric chain model [14]. The experimental method of the studies has been discussed above; the setup is outlined in Fig. 1.
The SiPM under study was Hamamatsu MPPC S10362-11-050C with an area of 1 mm 2 , a pixel size of 50 µm. The MPPC response is shown in Fig. 2. Later, these results have also been used for the evaluation of the GP model [24], as shown in Fig. 20. The GP model fit function has been constructed by convolution of the discrete G P and continuous Where x is an X-axis charge value, M peak is a number of the spectrum peaks taken into account (formally, M peak should be infinite; in practice, it is set a few times larger than a number of the resolvable peaks), μ 0 and σ 0 are a mean and standard deviation of the zero-event peak (pedestal) correspondingly. The probability of crosstalk for these data is 8.3% using the standard method (Eq. 27) and 9.0% by fitting the GP model (Eq. 32). However, both the geometric and the GP distributions approach each other at low crosstalk p and, finally, to the Poisson one at p = λ = 0.
Participating in the CERN CMS HCAL upgrade to develop and evaluate calibration and monitoring procedures for the HEP17 endcap (HE) module (see [36]), the MEPhI SiPM Lab team collected and analyzed about 200 SiPM spectra. HEP17 SiPMs are Hamamatsu MPPCs dedicated to CMS HE Phase 1 upgrade S10943-4732 of 15 µm pixel size, diameter 2.8 mm. The readout of the SiPMs was performed using the QIE charge integration and encoding ADC cards (version 11) developed by Fermilab (see [36,37]). The readout data from the QIE11 is continuously transmitted to the backend electronics for data processing. A calibration unit with LED inside was used for the custom light illumination of the SiPMs.
The vast majority (94.5%) of the spectra have been properly fitted by Eq. (32), including intentionally blurred ones by constant background illumination to mimic the high DCR anticipated by the CMS community for radiation-degraded SiPMs [22,25]. A typical example of the spectrum of 30K events at room temperature and its processing results are shown in Fig. 21.
A dedicated experiment on an initial evaluation of the skewness-based method has been carried out in collaboration with the KETEK team for the SiPM PM3350-WB of 3 × Histograms of the dark and light measurements are shown in Figs. 22 and 23 for visual perception. GP fitting has not been applied to them due to the non-Gaussian form of the peaks, as well as because the focus of the study is a direct calculation of the mean, standard deviation, and skewness of the raw data by the skewness method. The non-Gaussian form of Fig. 24 The crosstalk probability of the SiPM PM3350-WB measured by standard and skewness-based methods. The trendline is defined by a second-order polynomial function Fig. 25 The gain of the SiPM PM3350-WB KETEK measured by standard and skewness-based methods. The trendline is defined by a linear function the peaks, namely, pronounced left-side tails of the zero and all following peaks, is expected to be caused by the AC coupling of the amplifier. An initial positive portion of the bipolar photoresponse contributes to the peak bodies, whereas some portion of the DCR SER tails of negative polarity arrived before the integration gate, contributing to the left-side peak tails. Moreover, the KETEK PM3350-WB SiPMs are characterized by a long-tail SER with a slow recovery time of about 100 ns which complements this explanation. The pronounced left-side tails of the first-and higher-order peaks could also be caused by a low-gain SER of partially charged (incompletely recovered) cells; however, such an effect does not explain the similar tail of the zero-order peak.
The main experimental results of the study are shown in Figs  We also estimated a mean number of dark primaries and corresponding DC R (Fig. 27) in two different ways: -by processing waveforms, extracting interarrival times of dark events, and applying CCDF method as the standard one (see details in [38,40]) to find a mean interarrival time of Poissonian events ∼ 1/DC R; -by calculating a mean number of dark primaries in t gate = 17 ns from raw data using the skewness method.
Let us note the good agreement between both the skewnessbased parameter estimations and the standard ones in the overvoltage range 1.4-4.4 V, especially for the gain, primaries, and DCR. For the primaries (Fig. 26) a trendline has been calculated using a single-parameter model of an avalanche triggering probability proposed by Otte [39], where the so-called Otte parameter of the fit is found to be 0.049.
These results are summarized in Table 1 where the experimental relative errors were calculated by mean square errors M SE of the scatter of the corresponding (standard and skewness) points about trendlines. The trendlines were fitted to the points measured by the standard method to obtain a minimum of the M SE (so-called MMSE method).
The theoretical relative errors of the skewness method were calculated using the expressions of the full analytical model (Eqs. 24, 25) applying the correction factor F DC R . Let us note the rather small influence of this correction on the errors B and S E ∼ √ F DC R ≤ 1.08 (F DC R ≈ 1.17 for the maximal observed DC R ≈ 4.3 MHz at the overvoltage of 6.4 V).
The experimental errors of the skewness method are in good agreement with their theoretical (model) values for all parameters except for a dark crosstalk value of 14% (obviously, due to one outlier point at 3.4 V -see Fig. 24). The skewness errors are also in a good agreement with the standard method errors which are the minimal ones by definition of the MMSE method. And the standard method errors have reasonable values. Indeed, for example, if one directly counts 50K Poisson events, then S E rel (η = 0.8) = 1/ √ n · η = 0.5%. In case of the SiPM characterization by Pr(0) of the zero peak, some uncertainty in the separation of the zero and the first peaks makes S E rel (η) higher up to 1%. The same level of S E rel (η) ≈ 1% is reported in the comprehensive SiPM study [39]. However, the main concern is the correctness of the skewness method in the high-overvoltage range 4.4-6.4 V as summarized in Table 2. The experimental errors considerably exceed their theoretical predictions for all parameters except DC R.
First, let us consider these deviations using the analytical model results. The model (see approximated Eqs. A.3, 1) properly predicts the bias signs: B(λ) is negative, B(ĝ) is positive, B(η) is slightly negative (or close to zero at η ∼ 1) as observed for corresponding parameters (Figs. 24,25,26). But the values of the relative biases at an overvoltage of 6.4 V ( p ≈ 27%, g ≈ 0.03, η ≈ 0.8) should be less than 0.1% for all parameters while the observed relative biases are much higher (up to ∼ 30% for the crosstalk).
The model predicts that the standard errors are much larger than the biases discussed earlier, and the last columns of Tables 1 and 2 are determined by contributions S E rel at the level of a few percent. However, the observed deviations in the high-overvoltage range of 4.4-6.4 V cannot be reasonably explained by the analytical model "as is".
Second, let us consider MC simulations with dark noise to check whether the observed DC R could be a source of incorrectness. As discussed above (Sect. 4.2), if DC R · t gate ≤ 0.1η then the influence of DC R on skewness S E is rather low and can be approximated by the correction factor F DC R . At an overvoltage of 6.4 V with DC R ≈ 4.3 MHz and η ≈ 0.8, the mean number of dark primaries in t gate = 17 ns is about 0.07 < 0.1 η. Therefore, the correction of the standard errors by √ F DC R ≈ 1.08 seems to be a valid estimation of the dark events contribution. The skewness biases for DC R · t gate ≤ 0.1η are negligible and not affected by DC R (see Figs. 14, 15, 16). Thus, the observed deviations of the skewness results from the standard ones at the high overvoltages, especially large biassness of the crosstalk values, appeared to be a few times larger than the MC simulation results and analytical predictions for the low DCR range.
Overall, the reasons for these discrepancies are not clearly understood yet.

Discussion
This study was aimed, first and foremost, at the development and analysis of the theoretical grounds of the skewness-based SiPM characterization method based on the GP model. The GP distribution of the SiPM response is confirmed by many independent experiments with various SiPMs. However, fitting the GP model with the data requires many parameters of negligible or minor importance (e.g. baseline offset μ 0 , electronic noise RMS σ 0 , ENF of gain F G ). As soon as the SiPM charge response corresponds to the GP model, the charge moments should also correspond to those of the GP distribution as well. Therefore, the method should work in such cases.
However, there are other known models of the SiPM response. First, the advanced GP model [28] has been discussed earlier. Second, the 'nearest neighbors' model explicitly accounts for the ignition of pixels adjacent to the primary one by the crosstalk process [41]. Third, the Erlang distribution model has been proposed as just a better fit function for some experiments of the FACT project [42]. Some of the published results show minor differences between the models and the GP distribution (e.g. [41]). Probably the same situation would be reflected in the moments and, finally, in the parameters estimated by the skewness method. Therefore, it makes sense to consider a possible extension of the skewness method to other distribution models.
As anticipated, the skewness method could have many advantages for large-scale experiments in automation, mass testing, and continuous monitoring of SiPMs: -utilizing the simplest processing of the raw data; -predicting optimal measurement conditions by simple approximations of statistical errors (Appendix A); -providing guess values of the key parameters for reliable fitting of the spectrum; -providing some immunity to dark and electronic noises up to the level of about 0.1-0.3 noise events in the integration time gate per single photo-primary event.
Remark: the limiting noise level is simulated for dark events and the same is expected also for an equivalent noise charge in SER units. Unfortunately, an initial expectation that the method could be applied for the characterization of radiation-degraded SiPM with unresolvable spectra as a result of high DCR cannot be confirmed. However, a spectrum barely resolvable (see Fig. 13) with parameters η = 3, λ = 0.1, n = 10K , DC R = 30 MHz (equal to 3 dark primaries in the time gate) has been characterized with relative bias B rel (λ) ∼ 0.1, relative standard error S E rel (λ) ∼ 0.8 (that is, absolute S E(λ) ∼ 0.08 at λ = 0.1) and much lower errors for the gain and primaries (see Figs. 14,15,16,17,18,19).
Moreover, characterization of the definitely unresolvable spectrum of 30 primaries with DCR of 30 MHz has been demonstrated by MC simulations (a histogram is not presented, but all results are also shown in Figs. 14,15,16,17,18,19). In this case, the number of primaries (signal) is ten times higher, the statistical errors are predictably larger according to the model, but an immunity to DCR is also ten times higher (shift of the S E rel rise from ∼ 5 MHz for η = 3 to ∼ 50 MHz for η = 30).
Therefore, more studies are anticipated on the feasibility of the method and its possible optimization for high DCR.
We also note the ability of the method to characterize dark SiPM parameters by Q dark from measurements of Q of f with subtraction of the pedestal Q el . Remarkably, that the SiPM parameters obtained by dark measurements are rather close to the light ones with comparable errors, especially for the gain (Fig. 25, Table 1). And the DC R values determined by the skewness method agree well with the standard values (Fig. 27).

Conclusion
A new method of the skewness-based characterization of the main SiPM parameters -a mean number of primary events, gain, and probability of correlated events -has been proposed, theoretically studied, and experimentally evaluated. The method provides the extraction of these parameters from the raw data of the SiPM response using the simplest processing algorithm: calculating the mean, standard deviation, and skewness of the data set. In this study, the method has been implemented using the assumption of the Generalized Poisson distribution of the SiPM response, which became a well-recognized model confirmed by many experiments for a variety of SiPMs.
The developed analytical model of the method errors predicts an acceptable accuracy and precision of the SiPM characterization starting from the sample size (number of charge measurements) of 100K events in a typical range of the SiPM parameters: the mean number of primaries ∼ 1-10, crosstalk probability up to 30%, and for arbitrary gain. That is, it provides the absolute standard error of the probability of crosstalk of less than 1% as well as the relative standard errors of the gain and the mean number of primaries of less than 0.03. The standard errors are scaled as an inverse square root of the sample size, thus the precision could be improved at a cost of larger statistics. The biases of the parameter estimates are much lower than the standard errors starting from 100 measured events because they are inversely proportional to the sample size. The skewness statistical errors could be improved ten times for n = 1M.
Monte Carlo simulations with DCR ranging from 1 to 100 MHz revealed an influence of the dark noise on the skewness method errors: -is negligible if a mean number of dark events at the time gate per detected photon is less than 0.1 (DC R < 0.1η/t gate ), and the errors correspond to the analytical model; -is low and predictable if DC R < 0.3η/t gate , and the errors increase with respect to their zero-DCR analytical values in correspondence with the correction factor F DC R ∼ 1.6 times; -is high and unpredictable by the model if DC R > 0.3η/t gate , and the standard errors increase by about an order of magnitude.
The initial experimental validation of the skewness method allows us to conclude that in the overvoltage range 1.4-6.4 V it agrees with the standard ones, namely: for the gain with the standard error of 1.8%, for the crosstalk probability with the standard error of 5.4% (light measurements, not dark ones), for the mean number of detected photons with the standard error of 1.5%, and for DC R with the standard error of 2.5%. The experimental errors of the skewness method are also within the limits of the analytical predictions.
In the overvoltage range 4.4-6.4 V the skewness method disagrees with the standard ones for all parameters except DC R. The origins of the disagreements are not yet clearly understood. availability of these data, which were used by the author under a nondisclosure agreement, and so are not publicly available.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: Appendix section
We consider the bias B(θ) and variance V ar(θ) of the parameter estimatorsθ = (η,λ,ĝ) given by the expressions 24 and 25.
The partial derivatives of the function f i (m 1 , m 2 , m 3 ) i = 1, 2, 3 can be obtained from the corresponding expressions (16)- (18). It is a straightforward calculation, so let us skip its consideration in detail.
Covariance matrix Cov(m j , m k ) j, k = 1, 2, 3 of the sample moments is defined by known expressions (see, for example, [31]). They are reproduced below in a simplified form assuming a large sample size (n 1): The covariance matrix A.1 comprises the population moments up to the sixth order. They can be found in the section "Generalized Poisson distribution" of the book [30] and can be retrieved by the Wolfram "Mathematica" software. Moments of orders 1-3 are given in the theoretical section by the Eqs. (11)(12)(13); the higher-order moments are shown below: (1 − λ) 10 + 1 + 52λ + 328λ 2 + 444λ 3 + 120λ 4 (1 − λ) 11 g 6 (A.2) Calculating the sums in the Eqs. (24) and (25) using the covariance matrix A.1 with arguments from Eqs. (12)- (13) and (A.2) we derive the final expressions for the bias B(θ) and variance V ar(θ). Regretfully, being rather long, they could hardly be conveniently displayed. However, the first order Taylor expansion for λ around λ = 0 allows to get rather compact approximations of the errors of reasonable accuracy if λ 1 (λ < 0.1 is a typical case for modern SiPMs): Let us note some features of the results: -The biases and variances of the estimatesη andλ do not depend on the gain; -B(ĝ) and V ar(ĝ) are linear and quadratic functions of the gain correspondingly; therefore, the relative errors do not depend on the gain; -B(λ) is a strictly negative function, moreover, the higher η the lower bias (see Fig. 3); -B(η) is close to zero at η = 1 and λ 1 (see Fig. 5); The approximations are found to show relative deviations from the exact calculation results of about 10% or less if the primaries and crosstalk parameters are within the ranges 1 < η < 10 and 0 < λ < 0.1 correspondingly. Higher crosstalk values inevitably make accuracy worse, for example, up to about 50% for λ ∼ 0.2. However, the exact analytical results should be used in this case.