Experimental observation of anomalous supralinear response
of single-photon detectors

We performed an analysis of nonlinearity measurement uncertainty $\sigma(\Delta)$ to find its ultimate physical limits and find out whether the data approach the fundamental nonlinearity resolution. The statistical uncertainty was calculated for the case of a balanced experimental setup $R^{\text{det}}_{\mathrm{A}}=R^{\text{det}}_{\mathrm{B}}$, measurement time $T$ = 20 s for each rate, and $N=30$ repeated measurements.

The fundamental limit is represented by the shot noise of the Poisson counting process, where the measured number of events ($R^{\text{det}}T$) exhibits variance equal to the mean value. This implies the standard deviation of the detection rate $\sigma(R^{\text{det}})=\sqrt{R^{\text{det}}/T}$. If we introduce this standard deviation to each detection rate in the nonlinearity formula

through standard error propagation $\sigma^{2}(\Delta)=\sum_{i}\left[\sigma(R^{\text{det}}_{i})\times\partial\Delta/\partial R^{\text{det}}_{i}\right]^{2}$, we obtain the standard deviation

The measured nonlinearity $\overline{\Delta}$ is averaged from $N$ measurements, so the standard deviation of the result is

This limit is compared to the experimental values in Fig. 7. The measured standard errors of the nonlinearity for SPADs meet their lower bound $10^{-4}$ for rates higher than $2\times 10^{5}$ cps (Fig. 7(a-d)). For higher rates, excess noise limits the measurement uncertainty; it nevertheless remains below $\sigma(\overline{\Delta})<10^{-3}$. In this region, the mean nonlinearity value is larger than the uncertainty by at least one order of magnitude. Subfigures (e-l) show the nonlinearity measurement uncertainty of the SNSPD for several different values of bias current. Latching rapidly increases for bias current $I_{\text{bias}}>25$ µA and brings extra uncertainty to the measurement (Fig. 7(k) and (l)).

The shot noise limit is well applicable to Poisson count rates. However, dead-time saturation will further affect the measurement, because the measured count rates become sub-Poissonian as they approach the limit $R^{\text{det}}\to 1/\tau_{\text{NP}}$. Because the standard deviation of accumulated counts is complicated to calculate analytically [38], we are going to assume a sufficiently long integration time so that the number of accumulated counts is large: $\langle n\rangle=R^{\text{det}}T\gg 1$. This is a reasonable assumption for high rates, where saturation occurs. Then, we can make the following derivation.

Upon measuring $n$ detections during a time $T$, we obtain the average time between detections, $\overline{\Delta t}=T/n$. Because each realization of $\Delta t$ is independent, the average time $\overline{\Delta t}$ is normally distributed with a standard deviation $\sigma(\overline{\Delta t})=\sigma(\Delta t)/\sqrt{n}$ according to the central limit theorem. We assume $\Delta t$ to be a dead time $\tau_{\text{NP}}$ plus an exponentially distributed delay with a rate parameter $\lambda$, which makes $\langle\Delta t\rangle=\tau_{\text{NP}}+1/\lambda$ and $\sigma(\Delta t)=1/\lambda$. Because $n\gg 1$, the standard deviation is small relative to the mean; $\sigma(\overline{\Delta t})\ll\langle\Delta t\rangle$. As a result, the mean rate $R^{\text{det}}=1/\overline{\Delta t}$ is also normally distributed and the standard deviation is scaled based on the slope around the mean value,

Due to the large number of detections, on the right side we can simplify $\overline{\Delta t}\approx\langle\Delta t\rangle\approx 1/R^{\text{det}}$, so

The resulting sub-Poissonian bounds are shown in Fig. 8 for two dead time values to illustrate the cases of an actively and passively quenched detectors. The bounds do not make a difference within the data range of Fig. 7.

Here we address the optimal allocation of measurement time to minimize the uncertainty. Given the overall measurement time $T_{O}$ per one sample of $\Delta$, there exists an optimal distribution between measuring $R^{\text{det}}_{\text{A}}$, $R^{\text{det}}_{\text{B}}$, and $R^{\text{det}}_{\text{AB}}$. Assuming a balanced splitting $R^{\text{det}}_{\text{A}}\approx R^{\text{det}}_{\text{B}}$, the minimum uncertainty of $\Delta$ is reached for

In the regime where $|\Delta|\ll 1$, the approximate ratios are $T_{\text{A}}:T_{\text{B}}:T_{\text{AB}}=0.3:0.3:0.4$. The assumption we made on the way is that the measurement time is long enough so that all the accumulated counts $C_{i}$ ($i=\text{A, B, AB}$) are approximately normally distributed. This holds if $\langle C_{i}\rangle=R^{\text{det}}_{i}T_{i}\gg 1$.

There have been many SPAD response models proposed that model saturation effects and noise [78, 57, 64, 34, 35, 36, 37, 38]. The model used to fit the presented results considers a non-paralyzable dead time $\tau_{\text{NP}}$, dark count rate $R_{0}$, and neglects afterpulsing. The rate of detection events as a function of the incident rate $R$ then follows a standard formula

In Fig. 10, we illustrate the effects of dead time and dark counts, when this model is applied to nonlinearity (1). The left slope represents the background noise characterized by dark counts and the right slope represents dead time saturation. Other theoretical models shown below yield similarly V-shaped nonlinearity $\Delta$ with analogous scaling.

To correct for nonlinearities under this model, one can apply the formula

with $i\in\{\text{A, B, AB}\}$ to the measured detection rates. Fig. 11 shows the nonlinearity values calculated from corrected rates. The same data as presented in the main text are used, and the values of $R_{0}$ and $\tau_{\text{NP}}$ are taken from the least-squares fits (see Fig. 12 and Table 3 below).

If the model (22) holds, then the corrected nonlinearity should be close to zero. However, due to the Poissonian uncertainty (8) and the nonlinear nature of the formula (4), there is a nonzero bias (solid line in Fig. 11). The exact probability distribution of $\Delta^{\text{corr}}$ needs to be computed numerically for $R^{\text{corr}}_{\text{AB}}<10^{2}$ cps, because a normal approximation no longer holds [74].

Like in the main text, each point represents the average from 30 measurements. To test the compatibility of the corrected data and the model (22), we test both the average value and the distribution of the 30 samples for each point. The average nonlinearity is compared to a 95% confidence interval, and the Kolmogorov–Smirnov (KS) p-value of the distribution is denoted by the point color. The p-value is a probability that, assuming the model holds, we would obtain “worse” data than those measured, as parameterized by the KS statistic (see Statistical Methods below). Very low p-values show that statistical errors alone are unlikely to explain the data.

Many points show a disagreement with the theoretical model, with the difference standing out more in comparison to uncorrected data (Fig. 12), especially for $R_{\text{AB}}^{\text{det}}\gtrsim 10^{6}$ cps.

Accordingly, the basic model (22) does not provide a satisfactory explanation of the measured data, so we should justify its use as opposed to other dead-time models, and quantify the effect of afterpulsing.

The effect of dead time in Geiger-mode SPADs is akin to Geiger-Müller (GM) counters [79, 58, 80, 59, 81, 78, 82, 60, 61, 62]. Two basic types of idealized models for dead time have been defined [79, 59, 60]. Namely, it is the paralyzable dead time $\tau_{\text{P}}$ model

and the non-paralyzable dead time $\tau_{\text{NP}}$ model (22). For the non-paralyzable case, each registered detection is followed by dead time, during which no further events are registered. In the paralyzable case, dead time follows every photon absorption, even those that occur within a previous dead time and are not otherwise recorded. This case covers the fact that secondary detections in GM tubes still require quenching, but are not recognized due to low voltage output.

In the case of GM counters, single-parameter dead-time models are just an approximation. Hybrid models were proposed by combining paralyzable and non-paralyzable dead times [63]. There are two variants; the NP-P model

and the P-NP model [61, 83]

We have tested whether these theoretical models can be used empirically to fit the measured nonlinearity data. We found that all hybrid-model fits converge to either the paralyzable or non-paralyzable case. Fig. 12 shows both of these response models applied to the measured nonlinearity of all detectors. A complete list of the measured parameters and best-fit parameters is given in Table 3.

As in the main text, we evaluated $\chi^{2}$ to demonstrate that investigated models significantly deviate from the measured nonlinearity. In Table 4, we show the chi-squared per one degree of freedom, $\chi^{2}/\nu$, where $\nu$ is the number of data points minus the number of fitting parameters.

Both the $\chi^{2}$ values and the plots in Fig. 12 show that the paralyzable and hybrid models do not offer a better explanation of the data; nor can they satisfactorily match the trend in the middle of the graphs of SPAD-1–3, where all nonlinearities seem to be systematically lower.

Actively quenched SPADs exhibit afterpulsing and twilight pulsing that affect the mean detection rate [38]. Both effects can be evaluated numerically [84], or – if we neglect the temporal distribution of afterpulses – an approximate rate formula can be used [38],

The new parameters are the mean number of afterpulses per detection $\langle n_{\mathrm{AP}}\rangle$, and the twilight-pulse proportionality constant $\alpha$ that introduces rate dependence. As one would expect, when both of these parameters are zero, the formula (27) is reduced to the basic non-paralyzable model (22).

The key observation here is parameter degeneracy. In the formula (27), the effect of afterpulsing can be substituted by adjusting the dead time and dark count variables in (22). If we put the model parameters in square brackets, the equivalence can be expressed as

Let us now substitute both rate models into nonlinearity (1), assuming balanced splitting $R^{\text{det}}_{\mathrm{A}}=R^{\text{det}}_{\mathrm{B}}$ and the parameter equivalence (28),

One finds that $\Delta_{\text{AP}}\equiv\Delta_{\text{NP}}$, which follows from the equivalence (28) and the scaling of $R$ by a linear factor there. The principle is that both $f_{\text{AP}}^{-1}$ and $f_{\text{NP}}^{-1}$ yield incident rates that differ by the same linear factor, and their ratio does not change upon the multiplication by 1/2. A subsequent application of $f$ therefore gives identical rates for both models. As a result, all fitted models are identical and yield the same $\chi^{2}$ either with or without afterpulsing.

A further consideration would be a more complex afterpulsing model that considers the temporal distribution of afterpulses [38]. This model is based on computationally demanding Monte-Carlo simulations and is therefore unsuitable for least-squares fitting. However, the difference between this full model $\Delta_{\text{full}}$ and the simple model (27) can be evaluated for an example case. In Fig. 13, the difference is shown for SPAD-1 based on the afterpulsing data measured in Ref. 38. Note that we have already established the equivalency of (22) and (27), so we can conclude that the full model introduces a very small correction to the basic model we used to fit the data.

An example case of an actively quenched SPAD, represented by SPAD-1, shows that an afterpulsing correction to our nonlinearity model would yield a relative difference of $\approx 0.1\%$ in $\Delta$. A slightly simplified afterpulsing correction that can be used for data fitting yields no difference from the basic model. For this reason, we can conclude that afterpulses do not play a role in fitting the nonlinearity data of actively quenched SPADs.

The other detectors exhibit a much more complex behavior. Passively quenched SPADs exhibit paralyzable dead time, but also a rising efficiency curve after each detection [85], which makes the model too complex to be parameterized by a few measurable quantities.

SNSPDs show no afterpulsing, but, to our knowledge, their (non-)paralyzability has not been confirmed. Additionally, both SNSPDs and passively quenched SPADs exhibit changing efficiency after each detection, an effect that is both non-negligible and complex. In the case of SNSPDs, the single-photon and two-photon efficiencies also depend on the bias current. For these reasons, no standard rate model has been formulated that would be sufficiently accurate.

Upon exploration of certain empirical models combining non-paralyzable and paralyzable dead time, we found that none of them offer a significant advantage over the basic non-paralyzable model. SPAD-1–3 are non-paralyzable. SPAD-2 and SPAD-3 are better fitted with the paralyzable model, but the difference is very small, as can be seen from Fig. 12. SPAD-P is significantly better fitted by the non-paralyzable model. Consequently, we use a single model (22) for all SPADs in the manuscript.

The most significant supralinearity ($\Delta<0$) among SPADs is exhibited by SPAD-3. Here we use an ad hoc model to fit the data and assess the feasibility of various factors as possible explanations of the supralinearity.

The data can be fit with a modified version of formula (22),

where the term $a_{0}$ corresponds to a constant dark count rate, $a_{1}$ corresponds to a constant detection efficiency, and $\Phi$ is the incident photon flux. The last term (with $b$ = 1.004 and $R_{1}$ = $10^{6}$ cps) is empirical and can be attributed to either of these parameters as a hypothetical rate-dependent term.

Fig. 14 depicts the hypothesis (A), and the corresponding dependence of dark counts (B) and efficiency (C) that would be required to match the hypothesis. The dark count rate dependence is arbitrary up to a linear factor $a\Phi$, but we can calculate the absolute minimum dark count rate directly from the data (green points in Fig. 14B).

Dependence on the dark count rate was ruled out by an independent pulsed measurement, where the count rate can be inferred from detections registered outside the pulse window. The dark count rate was found to be below $30$ cps for a signal of $10^{5}$ cps mean detection rate, whereas the nonlinearity data would require at least $R_{0}(10^{5}~{}\text{cps})$ > 1000 cps. This means that rate-dependent dark count rate is not the cause behind the supralinearity of SPAD-3.

If we attribute the supralinearity to rate-dependent efficiency, we are left with relative changes depicted in Fig. 14C. While efficiency decreasing with rate has recently been reported in Ref. 65 for InGaAs SPADs subjected to pulsed signals, here the efficiency would need to increase. For a SPAD, the efficiency increases with bias voltage, which should be fixed in its steady state by a DC voltage supply. However, we need to rule out any temporal dependence of the bias voltage, which is addressed in the next section.

A possible explanation for the supra-linear response of a SPAD would be time-dependent efficiency after each detection. When an actively quenched detector resets, the overall response may fluctuate [57]. In order for a detector to exhibit supra-linearity, the dominant effect would have to include decreasing efficiency over time after each detection. The slope of this decrease needs to be significant on time scales close to inverse mean count rates where supra-linearity occurs ($\sim$10–100 µs). This could, in principle, be caused by a settling bias voltage, which affects detection efficiency.

In an attempt to roughly estimate the necessary magnitude of such an effect, we modeled the efficiency with an exponential. The results are shown in Fig. 15. To achieve supralinearity roughly similar to the data observed on SPAD-3, the relative decrease in efficiency would have to be $2.5\%$ over more than 0.1 ms. However, as the quenching electronics operate on a GHz time scale, such a slow and significant effect is extremely unlikely.

This effect can be ruled out based on interarrival histograms of detections under continuous illumination. For constant efficiency, interarrival times follow an exponential distribution, save for dead time (here 57 ns) and afterpulsing effects that take place on time scales <$10$ µs. For time-dependent efficiency, the interarrival time $\Delta t$ is distributed with the probability density

where $R$ is the incident rate and $\eta$ is a relative term describing the changes in detection efficiency. Interarrival histograms then sample this distribution.

Let us assume the hypothesis shown in Fig. 15 that $\eta(\Delta t)=1+0.025\exp(-15\Delta t\times\mathrm{ms}^{-1})$. Then, for $\Delta t\gtrsim 0.2$ ms, $\eta(\Delta t)\approx 1$ and the distribution will scale as an exponential—$p(\Delta t)\propto\exp(-R\Delta t)$—which is the same scaling as for constant efficiency. On shorter time scale, $\Delta t\lesssim 0.2$ ms, the distribution will be more complex.

Data in Fig. 16 show the measured interarrival histograms relative to the two models—constant and settling efficiency. Because the difference is small, the histogram values $p_{\text{meas}}(\Delta t)$ are divided by the expected values $p(\Delta t)$ to visually distinguish the two cases more clearly. A horizontal trend means that the scaling corresponds to the expectation up to a constant proportionality factor. It can be seen that the data follow a negative exponential for constant efficiency, whereas the comparison to settling efficiency does not result in a horizontal data trend for $\Delta t<0.2$ ms. This provides good evidence that the efficiency settling effect is not the reason for supralinearity of SPAD-3.

Here we summarize the factors that we rule out as causes of supralinearity, as demonstrated on SPAD-3.