1. Introduction

Single-photon avalanche detectors (SPADs) are critical for many areas of research in atomic physics, quantum optics, and quantum communication as well as a variety of other fields, including biology at the single-photon level [1–21]. The operating principles of SPADs have been described in detail [22,23]. Briefly, SPADs are photodiodes reverse biased well above their breakdown voltage. In this state, when a single photon impinges on the photodiode and excites an electron to the conduction band, the large bias accelerates the electron to high enough energies to cause impact ionization, resulting in an “avalanche” of current. This avalanche is sensed by a quenching circuit [23], which lowers the bias until the avalanche ends, and then increases the bias back to its normal operating value. Importantly, during the avalanche, an additional photon may arrive and excite a new electron, but the noise on the avalanche is so large that this additional excitation is unresolvable. Thus, SPADs are inherently unable to resolve photon number, unlike other detectors that have at least some photon-number-resolving capability [24–26]. This complicates what might be one of the most obvious uses of devices that can detect single photons: number-state reconstruction. To obtain photon-number information, a common approach is to assume a SPAD produces either zero or one click per experimental attempt, and this “click/no-click” information is the basis of several techniques that use multiple detectors or multiplex a single detector [27–36]. These techniques by their nature increase costs or complexity of implementation.

In this work, we perform photon number-state reconstruction using a single SPAD with no inherent photon-number-resolving capability. In our scheme, light pulses and correlation time scales longer than the detector dead time allow the SPAD to click more than once in a measurement period. This simple scheme, when combined with maximum-likelihood analysis, allows for remarkable performance with minimum cost and complexity, while the requirement that light pulses are longer than the detector dead time is compatible with many modern atomic physics and quantum optics experiments. Our algorithm also accounts for various SPAD imperfections, namely detector efficiency less than unity, dark counts [37], afterpulsing [38–40], and dead times [22].

We first develop a model for the SPAD, discussed in Section 2 (and detailed in Section A of Supplement 1). Our model requires detector characterization, which we detail in Sections B and C of Supplement 1. We use this data to perform the reconstruction via the expectation-maximization-entropy (EME) algorithm [41], reviewed in Section 3. The implementation of this model and the EME algorithm in code have been made open access [42]. To validate the algorithm, we use carefully calibrated coherent pulses, whose number distribution we know a priori. The results are discussed in Section 4. For coherent-state pulses with peak count rates of up to a few Mcounts/s, the total variation distance we achieve is below $10^{-2}$, both for simple square pulses and pulses with nontrivial time profiles. We also explore the algorithm and its limitations at high count rates, approaching one photon per dead time. Even at these high count rates, reasonably good agreement is achieved, depending on detector parameters. We also show that the reconstruction of non-classical light achieves good agreement with an independently measured value of the second-order pulse-averaged autocorrelation function $g^{(2)}(0)$. Thus, we show that a single inexpensive SPAD is capable of performing number-state reconstruction for a wide range of input light statistics.

2. Description and Model of SPADs

Given that the clicks produced by a SPAD do not correspond one-to-one to photons incident on the detector, we adopt the following formalism. Consider a pulse of light with a number distribution $\mathbf {P}$ (equivalent to the diagonal of the Fock-basis density matrix of the light). The ensemble measurement of this light produces a corresponding click number distribution $\mathbf {C}$. We model the detector as a matrix $\mathbb {D}$ relating these two vectors:

The matrix $\mathbb {D}$ in our model is the product of four individual matrices, each of which contains information about a particular imperfection of the SPAD. (1) The effects of the non-unit detector efficiency are captured by a matrix $\mathbb {L}(\eta _0)$, where we define $\eta _0$ as the detector efficiency when the detector is operating normally, i.e., recovered from all previous clicks. (2) The effects of background counts are accounted for by a matrix $\mathbb {B}(p_b)$, where $p_b$ is the probability of a background count in the data collection window. Here background events result from both dark counts [37] and ambient, non-signal photons. (3) Recovery time effects [22,40,43] are due to the quenching circuit turning the bias off and back on after a click, which takes a total time called the “recovery time.” There are two such effects, both accounted for in the matrix $\mathbb {R}$. One is losses due to the detector dead time, when the detector bias is below breakdown so that avalanching does not occur. The other effect is twilight counts, when photons produce clicks while the detector bias is turning back on. We treat the system efficiency as time-dependent after each click (see Fig. S1 of Supplement 1), and account for the time delay typical of twilight counts [43]. Rather than measuring both the time-dependence of the detector efficiency and the delay profile of the twilight counts [43], we assume a model, and measure the effective parameters of that model (see Supplement 1 for details). Thus, we write this matrix as $\mathbb {R}[\gamma (t), D(\tau )]$ to signify that it depends on both the photon profile $\gamma (t)$ and a model for the time-dependent detector efficiency $D(\tau )$. (4) Finally, afterpulsing effects [38,39] are modeled by the matrix $\mathbb {A}(p_a)$, where $p_a$ is the probability of one afterpulse click in the data collection window given the input photon profile. This probability is constructed from the detector’s temporal afterpulsing profile as well as its total probability of afterpulsing $\tilde {p}_a$. The measured values of the detector parameters relevant for all four matrices are summarized in Table 1; the methodology for these measurements is described in Section B of Supplement 1.

Table 1. Measured Parameters for Two SPADs under Test^a

View Table

These matrices are combined into $\mathbb {D}$ by multiplying in a physically motivated order. First, incoming photons are subjected to the loss of $\mathbb {L}$. Then background counts are added using $\mathbb {B}$. Both background clicks and clicks produced by signal photons are subject to recovery time effects encoded in $\mathbb {R}$. Finally, the clicks that remain may afterpulse, modeled by $\mathbb {A}$. Thus, we compute

3. Reconstruction Algorithm

Given the relation $\mathbf {C} = \mathbb {D}\mathbf {P}$, and given the measured click data $\mathbf {C}$ and the matrix $\mathbb {D}$, it seems natural simply to invert $\mathbb {D}$ to find $\mathbf {P}$. However, this direct matrix inversion may not produce a physically reasonable result, i.e., one which sums to unity and has no negative elements. If nothing else, physically unreasonable results may arise due to sampling error. This problem is exacerbated when, for instance, the detection efficiency is low.

Instead, we use the maximum-likelihood principle. The likelihood can be constrained via Lagrange multipliers and can be maximized over the relevant number of variables using any maximization method. Here we use an iterative procedure, a variant of the expectation-maximization (EM) algorithm known as the expectation-maximization-entropy (EME) algorithm [41]. We initialize the algorithm with a uniform distribution $\mathbf {P}^{(0)}_n = 1/(n_{\text {max}}+1)$. Then for all subsequent iterations, the distribution is

4. Results and Discussion

To validate our algorithm, we sent calibrated coherent light pulses into our detectors using the setup shown in Fig. 1. The light ($\lambda =780$ nm) is sent through a set of polarization optics for power control, as well as an acousto-optic modulator (AOM) to shape the pulses and a slow shutter for protecting the SPADs. The light is coupled into a fiber, whose output is passed through a calibrated neutral density (ND) filter with transmittance $T_{\text {ND}} = 0.001412(1)$. The light is then coupled into another fiber, which is inserted into either a SPAD or into a trap detector [45] used for calibration. The trap consists of multiple photodiodes arranged to allow very high efficiency and uniformity, making it well suited as a calibration transfer standard. The trap is read through a high-gain precision transimpedance amplifier (TIA) and an 8-1/2-digit digital multimeter (DMM). SPAD signals are sent to a time tagger module (TTM). Our TTM has a resolution of approximately 164 ps, and we use bins of approximately 1-ns width (6 TTM bins) for all of our reconstructions.

 

Fig. 1. The experimental setup for reconstruction measurements. Here, 780-nm light is sent through an acousto-optic modulator (AOM) and a shutter, then a system for power control consisting of a half-waveplate (HWP), a polarizing beam splitter (PBS), and a polarizer (Pol) on a motorized rotation mount. The light is then coupled into a fiber, whose output is passed free-space through a calibrated ND filter on a motorized flip mount, then coupled into another fiber, which may be connected to a trap detector [45] or to a SPAD. The trap detector is attached to a high-gain precision transimpedance amplifier (TIA) whose output voltage is read by an externally triggered digital multimeter (DMM), while the SPAD signals are sent to a time tagger module (TTM).

Download Full Size | PDF

To calibrate the average photon number per pulse, we connect the fiber to the trap detector, with no ND filter in the path, and record the voltage. The power incident on the trap is $P_0 = V_{\text {meas}}/(\eta _{\text {trap}}RG)$, where $\eta _{\text {trap}}$ is the efficiency of the trap, $R$ is the detector responsivity, $G$ is the TIA gain, and $V_{\text {meas}}$ is the average voltage reading. The calibrated ND filter is then inserted to reduce the power at the SPAD to $P_{\text {in}} = T_{\text {ND}}P_0$. Finally, while the measurement of $P_0$ is done with continuous wave (cw) light, the final reconstruction measurements are done on pulsed light, so we convert the cw power to a total average energy in a pulse via a factor $f_s$, which accounts for the pulse shape relative to the cw measurement and depends on the time window for the data collection (see Supplement 1). Thus, the expected average number of photons in each coherent pulse is

In addition, for all the reconstructions performed in this work, we report a reconstructed value of the pulse-averaged second-order autocorrelation function $g^{(2)}(\tau = 0)$, estimated as [22]

In Fig. 2, we present the results of our reconstruction algorithm for square coherent pulses, all with roughly the same peak input count rate of approximately 2 Mcounts/s. Results are shown for both SPAD1 and SPAD2. We change the length of the pulse to vary the average photon number, and show the results for three average photon numbers. For each of these pulses, and all coherent pulses used here, we collected $3 \times 10^7$ cycles of data. The pulse shapes as measured by the SPADs are shown in the insets. We used 1-ns binning for the reconstruction process (both here and for all of the reconstructions done in this work) and for the insets of the figure. For all pulses measured here, we find $\Delta < 10^{-2}$, and the calibration differences $\delta \bar {n}$ are all within 2% of the independently measured values, showing the accuracy of the algorithm.

 

Fig. 2. Reconstructed distributions (blue bars) and fitted Poissonians (black dots) for square coherent pulses for (a1)–(c1) SPAD1 and (a2)–(c2) SPAD2. A peak input photon rate of $\approx 2$ Mcounts/s is used for all pulses. For both SPADs and all average photon numbers, the distances $\Delta$ are less than $10^{-2}$ and all calibration differences $\delta \bar {n}$ are less than 2% in magnitude. The recovery time corrections are taken to second order for all results in this figure. The insets show the pulse as measured by the SPADs, with 1-ns bins, which is what we used in the reconstructions. The dip in power at approximately 100 ns after the start of the pulse is likely an artifact of the AOM drive turning on.

Download Full Size | PDF

In addition, in Fig. 3, we validate the time-dependent portions of the formalism by feeding non-square pulses to the SPADs. The pulses are shaped by using an arbitrary waveform generator to control the amplitude of the AOM drive. The shape was chosen to mimic single photons produced by cold-atom ensembles. To control the total energy delivered, we change the power through the first fiber in the setup using the polarizer. In Fig. 3, we show the reconstructions for $\bar {n} \approx 1, 5$ for SPAD1. Again we find $\Delta < 10^{-2}$ and calibration differences within 0.5% in magnitude, validating the time-dependent formalism.

 

Fig. 3. Reconstructed distributions (blue bars) and fitted Poissonians (black dots) of shaped pulses with SPAD1 at two different input photon rates. The peak input rates $r_{\text {max,in}}$ are as indicated. The insets show the pulse as seen by the SPADs, with 1-ns bins. The recovery time corrections are taken to (a) third and (b) sixth order.

Download Full Size | PDF

One of the contributions to the distances $\Delta$ of the reconstructions in Figs. 2 and 3 is drifts of up to a few percent in the optical power over the course of the hour-long measurement. The average of a drifting coherent state produces a distribution slightly wider than a Poissonian of the same average number, contributing to the observed deviations and the slight broadening of the reconstructed distributions, as indicated by $g^{(2)}(0)$ values slightly greater than 1. Based on Monte Carlo simulations, we believe these drifts contribute to total variation distances at the level of $1$-$3 \times 10^{-3}$, representing a significant portion of the observed $\Delta$ values.

Further, we test the limits of the algorithm at high input photon rates with $\approx 85$-ns full width at half maximum (FWHM) pulses. The results are shown in Fig. 4 for both SPAD1 and SPAD2 for three input rates. The reconstruction algorithm still performs reasonably well for SPAD1 (Fig. 4(a1)–(c1)), whose afterpulsing probability is small. Note that the highest photon rate here, four photons in 85 ns, is greater than one input photon per recovery time—i.e., the detector is dead for most of the pulse. This makes the near-agreements of the SPAD1 data with the expected distributions in Fig. 4 particularly noteworthy. For SPAD2 (Fig. 4(a2)–(c2)), which has a higher afterpulsing probability, the distributions are less faithful to the expected Poissonians. For both datasets, we observe increasing $\Delta$ between the reconstructed and expected distributions as the input count rate increases; we discuss possible reasons for this in the next section. In addition, we note the counterintuitive result that the largest values of $\delta \bar {n}$ in Fig. 4 occur for the lowest count rates, where corrections for afterpulsing, background counts, and recovery time effects are small. We believe these large $\delta \bar {n}$ values signify an unidentified systematic calibration error of $\sim 5$–$10$% in $\bar {n}_{\text {exp}}$ which affects all short pulses for both SPADs. This positive systematic error, combined with a decrease of $\bar {n}_{\text {fit}}$ away from $\bar {n}_{\text {exp}}$ (leading to a negative contribution to $\delta \bar {n}$), results in the appearance of improving $\delta \bar {n}$ with increasing count rate.

 

Fig. 4. Reconstructed distributions (blue bars) and fitted Poissonians (black dots) for coherent pulses of FWHM $\approx 85$ ns, using (a1)–(c1) SPAD1 and (a2)–(c2) SPAD2. The insets show the pulse as seen by the SPADs, with 1-ns bins. The recovery time corrections are taken to (a) seventh and (b),(c) eighth order. The peak input count rates $r_{\text {max,in}}$ are as indicated in the top center of each figure, corresponding to roughly (a) one photon per seven recovery times, (b) one photon per two recovery times, and (c) one photon per recovery time. The distances $\Delta$ increase with input count rate. We believe this increase is due to interactions between twilight counts and afterpulsing, as discussed in the main text.

Download Full Size | PDF

Finally, we investigated our ability to reconstruct anti-bunched light. Since we do not know the incident photon number distribution a priori, the metrics used above to evaluate our reconstruction of coherent states do not apply. As an alternative, we measure the pulse-averaged value of $g^{(2)}(\tau = 0)$ using a two-detector Hanbury Brown–Twiss (HBT) setup [22,46]. We compare the independently measured $g^{(2)}_{\text {HBT}}$ to the value estimated from the reconstructed distribution using Eq. (7).

To produce anti-bunched light, we used a Rydberg atomic ensemble described previously [46], with size $\approx 20\,\mathrm{\mu} \text {m} \times 8\,\mathrm{\mu} \text {m} \times 8\,\mathrm{\mu} \text {m}$ and Rydberg state $n=112$, and we performed approximately $2 \times 10^7$ data collection cycles. The light we produce is measured with the HBT setup to have a raw $g^{(2)}_{\text {HBT, raw}}(0) = 0.25(1)$ and $g^{(2)}_{\text {HBT}}(0) = 0.21(1)$ after background subtraction [46]. The light pulse shape is shown in Fig. 5(a).

 

Fig. 5. Reconstruction results for anti-bunched light using SPAD2. (a) Light pulse as seen by the SPAD, with 1-ns bins. The vertical black dashed lines indicate the 200-ns window used for reconstruction. (b) Reconstructed distribution $\rho _n^{\text {(recon)}}$. The recovery time corrections are taken to second order. (c) Ratio of the reconstructed state to a coherent state of the same average photon number, $\bar {n} = 0.036$, whose number distribution we denote $\rho _n^{\text {(coh)}}$. The $n=2$ and $n=3$ components of the reconstructed distribution are much smaller than those expected of a coherent state (gray dashed line), indicating anti-bunching. The one-standard-deviation uncertainties shown are computed through Monte Carlo sampling, as detailed in Supplement 1.

Download Full Size | PDF

The reconstructed distribution is shown in Fig. 5(b). To perform a reconstruction whose value of the correlation function $g^{(2)}_{\text {recon}}$ can be compared with that measured by the HBT setup, we reduce the reconstruction window to the one shown in vertical dashed lines in Fig. 5(a), which contains most of the photon pulse but avoids the tail. This tail is primarily coherent leakage into the SPADs, which is removed by the background subtraction in the calculation of $g^{(2)}_{\text {HBT}}$ but not by the background correction of our reconstruction algorithm. The resulting reconstructed distribution is compared with a coherent state of the same average photon number in Fig. 5(c). In the latter, the $n=2$ and $n=3$ components are well below what is expected of a coherent state, indicating the anti-bunched nature of the light. The reconstructed distribution is calculated to have $g^{(2)}_{\text {recon}}(0) = 0.2{2(2)}$, which is in excellent agreement with the independently measured value of 0.21(1). Thus, we demonstrate our algorithm’s ability to accurately reconstruct anti-bunched light.

5. Limits of the Algorithm

While the results presented here indicate a successful reconstruction process, the algorithm involves a number of assumptions and limitations. We describe the most important of these here, with a detailed discussion in Section A.V. of Supplement 1.

First, our algorithm is effective when both the light pulse and its correlations have time scales greater than a few detector recovery times. We have demonstrated the algorithm’s effectiveness for coherent pulses of lengths of at least a few detector recovery times, and for a pulse of non-classical light with pulse-averaged $g^{(2)}(0) < 1$. If the pulse or its correlations have time scales less than several detector recovery times, we expect the algorithm to fail. A clear example is that of broadband thermal light, whose bunching time scale is typically much less than 1 ns; in this case, photons “bunched” with previous photons will not register additional clicks, as the detector will be dead, and this will prevent us from correctly reconstructing the light statistics. Some light sources produce bunching with longer time scales [47], which we believe will be able to be reconstructed by our algorithm.

Second, our model of the detector weakens at high count rates. We believe the most important example of this breakdown, at least for our datasets, is that the model of independence between afterpulsing and twilight counts loses validity at high count rates. Although these two effects have independent origins, they “interact” by occurring on a single detector. This means that if the detector is dead, afterpulses from previous clicks may be blocked; it also means that in the few nanoseconds immediately after the end of a recovery time, the detector can only register one click even though both an afterpulse and a twilight count are theoretically possible. Both of these effects lead to an overcompensation of afterpulsing by our model, which should lead both to a decreasing $\delta \bar {n}$ with increasing count rate and a narrower reconstructed distribution. Both of these effects are observed in Fig. 4. This and other effects, such as non-Markovianity in afterpulsing, are discussed in detail in Subsection A.V. of Supplement 1. Thus, our algorithm, while still working well even approaching one photon per dead time, performs best for count rates that are small compared with the inverse of the detector dead time.

6. Conclusion

In conclusion, we have presented a method for performing number-state reconstruction with a single SPAD. The method works when the pulses of interest span at least a few multiples of the detector dead time, a regime that fits many of today’s quantum optics experiments. Further, the algorithm performs reasonably well even with input rates up to approximately one photon per recovery time. At these very high count rates, we have explored some of the ways in which the algorithm loses accuracy, which may be remedied in future work. This method could be extended to two or more detectors, yielding more information than standard multi-detector techniques and allowing reconstruction over a broader array of input light sources, such as sources with correlations shorter than the detector dead time.