## Abstract

In time-correlated single-photon counting (TCSPC), photons that arrive during the detector and timing electronics dead times are missed, causing distortion of the detection time distribution. Conventional wisdom holds that TCSPC should be performed with detections in fewer than 5% of illumination cycles to avoid substantial distortion. This requires attenuation and leads to longer acquisition times if the incident flux is too high. Through the example of ranging with a single-photon lidar system, this work demonstrates that accurately modeling the sequence of detection times as a Markov chain allows for measurements at much higher incident flux without attenuation. Our probabilistic model is validated by the close match between the limiting distribution of the Markov chain and both simulated and experimental data, so long as issues of calibration and afterpulsing are minimal. We propose an algorithm that corrects for the distortion in detection histograms caused by dead times without assumptions on the form of the transient light intensity. Our histogram correction yields substantially improved depth imaging performance, and modest additional improvement is achieved with a parametric model assuming a single depth per pixel. We show results for depth and flux estimation with up to 5 photoelectrons per illumination cycle on average, facilitating an increase in time efficiency of more than two orders of magnitude. The use of identical TCSPC equipment in other fields suggests that our modeling and histogram correction could likewise enable high-flux acquisitions in fluorescence lifetime microscopy or quantum optics applications.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. INTRODUCTION

Time-correlated single-photon counting (TCSPC) time stamps individual photon detections with picosecond resolution, and it has thus been used extensively in photon-starved applications such as quantum communications, lidar depth mapping, fluorescence lifetime imaging (FLIM), and non-line-of-sight imaging [1–7]. Through repeated illumination of a short laser pulse and detection with detectors such as single-photon avalanche diodes (SPADs) connected to fast timing electronics, TCSPC can be used to build up a histogram of photon detection times relative to the illumination time that reveals transient optical phenomena. Single-photon lidar (SPL) takes advantage of the sensitivity of TCSPC to individual photons in order to produce accurate low-light depth and reflectivity images from as little as one photon per laser position [8–10], which is especially critical for long-range imaging [11–13]. However, not all targets in a lidar acquisition will reflect so little light.

Strong ambient illumination or highly reflective targets may cause much more than one photon to be incident on the detector per illumination cycle. In atmospheric monitoring, for instance, backscatter from aerosols over a wide range of altitudes results in detection rates that vary over several orders of magnitude within a single measurement [14]. Unfortunately, a *dead time* follows each photon detection, during which additional incident photons cannot be registered. Especially when the rate of photons incident on the detector approaches or exceeds 1 per illumination, photons arriving during a dead time will be missed, and the detection time histogram will be distorted as a result [15]. Crucially, both SPADs and TCSPC timing electronics have independent dead times, with durations $t_{\text{d}}$ and $t_{\text{e}}$ for the detector and electronics dead times, respectively. A detection will cause both the detector and electronics to be dead, and if $t_{\text{d}} \lt t_{\text{e}}$ it is possible that a photon arriving after the end of detector reset but before the electronics reset causes a SPAD avalanche that is not registered by the electronics but that results in another detector dead time. The key challenge is accurately modeling how the overlapping dead times affect the probabilistic behavior of photon detection, so that dead time distortions can be mitigated.

Existing approaches to handling dead times are inflexible, incompatible with modern TCSPC systems and free-running SPADs, or consider only one source of dead time. The conventional approach is to keep the detection rate at most at 5% of the illumination rate [16], requiring attenuation to ensure that the probability of photons arriving during a dead time is negligible. However, attenuation in order to handle strong signals unnecessarily throws away informative photons that could be detected, increasing the number of illuminations required to sufficiently capture weak signals. Alternatively, the effects of dead time can be avoided by setting $t_{\text{d}} \ge t_{\text{e}}$ and $t_{\text{d}}$ to an integer multiple of the illumination repetition period $t_{\text{r}}$ [17], but $t_{\text{r}}$ and $t_{\text{d}}$ are generally not free parameters, instead being determined by the equipment or the required unambiguous range. Dead-time effects are simple to describe and mitigate for classical TCSPC systems allowing a single detection per illumination cycle, since detection times in different illumination cycles are then statistically independent of each other. Corrections can then be performed from the histogram alone without requiring the actual detection time sequence to be known [18] and can be combined with optimized attenuation to achieve accurate range estimates at flux values significantly higher than the conventional rule of thumb [19]. Desynchronizing the detector start time with respect to the illumination can further prevent classical pileup, especially when the ambient light is the dominant source of photons [20]. Modern TCSPC systems with free-running SPADs can detect multiple photons per illumination cycle, and dead times in one cycle could cause missed photons in a subsequent cycle, so a classical dead-time model only holds if the contribution from ambient light and dark counts is negligible and the signal response has finite duration [21,22], or if prior information is known about the time of the signal response, so that gating can be used [23]. Other dead-time corrections that allow for multiple detections per cycle consider a single source of dead time [20,24,25] and may assume gated operation that interrupts the detection time distribution from reaching a steady state [26–30]. To the best knowledge of the authors, [31] is the only work aiming to perform histogram correction that accounts for both $t_{\text{d}}$ and $t_{\text{e}}$ without the assumption of synchronization. While that work presents useful intuition for describing and correcting for dead time, the probabilistic model is inexact and thus does not perform well at very high photon flux.

In this work, we propose a method to perform ranging with SPL at high flux without bias, which is achieved by exploiting the intricate interplay among $t_{\text{d}}$, $t_{\text{e}}$, and $t_{\text{r}}$ and effectively mitigating the effects of the dead times on photon detection time distributions. We present a more comprehensive model of photon-counting systems and establish the effect of the two sources of dead time on the probability of photon detection when $t_{\text{d}} \lt t_{\text{e}} \le t_{\text{r}} \le 2t_{\text{d}}$; results for other combinations of ($t_{\text{d}}$, $t_{\text{e}}$, $t_{\text{r}}$) with $t_{\text{d}} \lt t_{\text{e}}$ can be derived in a similar fashion. Lidar experiments validate our modeling and demonstrate accurate depth estimation with photoelectron flux up to two orders of magnitude greater than the conventional limit.

## 2. SINGLE-PHOTON-COUNTING OPERATION

In the coaxial SPL configuration shown in Fig. 1, both illumination and detection share an optical axis. A laser illumination pulse passes through a polarizing beam splitter and is directed toward one point in a scene by a pair of galvo mirrors. Light reflected back from the scene reaches the SPAD after passing through a band-pass filter at the illumination wavelength and a lens focusing light onto the SPAD’s active area. An optional neutral-density filter may also be inserted in the detector path to reduce the flux incident on the detector. Time differences between photon detections and the most recent laser illumination are recorded by electronics in a TCSPC module and sent to a computer for processing. The following sections describe the basic operation of SPADs and TCSPC timing electronics to elucidate the cause of the dead time in each device.

#### A. SPAD Detectors

SPADs are reverse-biased photodiodes biased above the breakdown voltage with basic operation as follows [2,5,32–34]. When a photon hits a SPAD’s active area, a charge carrier (e.g., a *photoelectron*) may be generated via the photoelectric effect. Due to the reverse bias, generation of one carrier will further cause an avalanche of carriers, resulting in a current that is detectable as a digital signal. In order for the detection circuit to be sensitive to subsequent photon arrivals and to protect the avalanche from damaging the diode, the avalanche is quenched to reduce the bias below the breakdown voltage. In the actively quenched and recharged detectors that we assume in this work, the bias voltage is reset to its initial level after a fixed hold-off time ${t_{\text{ho}}}$, and the duration of the hold-off plus the reset is considered to be the detector dead time $t_{\text{d}}$. Not all photons incident on the detector actually cause an avalanche of electrons. The photon detection efficiency $\eta \in [0,1)$ is the probability that a photon incident on the detector will cause a detection event registered by the TCSPC system, which depends on whether an incident photon generates a carrier and whether that carrier initiates a detectable avalanche [35]. We assume the SPAD bias voltage is sufficiently high that the avalanche triggering probability is approximately unity, such that any photon converted to a photoelectron would cause an avalanche so long as the detector is not dead. Additionally, spurious avalanches occasionally occur that are not due to incident light, such as primary dark counts caused by thermal or tunneling-based carrier generation [36–38] or secondary afterpulses resulting from charge carriers trapped in semiconductor defect sites during an avalanche and then “detrapped” after the SPAD has been reset [39,40]. Our detection model considers the contributions from these sources to be negligible.

#### B. TCSPC Timing Electronics

While the SPAD provides the single-photon sensitivity and fast response, the precise timing is achieved using specialized TCSPC electronics. Classical TCSPC systems rely on an analog time-to-amplitude converter (TAC) and operate based on a stopwatch principle: a capacitor charges in the time between start and stop signals, and the capacitor voltage is read out via a high-resolution analog-to-digital converter (ADC) [2,41]. Importantly, TAC operation limits recording to at most one detection event per cycle. Modern TCSPC systems employ an alternative approach based on time-to-digital converters (TDCs) that compute time differences digitally [42]. TDCs enable a long measurement range, fast digital readout, monolithic integration, and the “multi-stop” capability of recording multiple events within one illumination cycle [43]. Due to speed limitations on digital-clock-based counters, many modern TDC-based TCSPC systems combine a coarse counter with finer timing that interpolates within clock cycles. The highest-resolution TDCs, such as the HydraHarp 400 (PicoQuant) used in this work, have historically required the fine interpolation to be performed by an analog component similar to a TAC [44,45]. A linear voltage ramp or pure sinusoid synchronized with the coarse clock is sampled and processed to achieve 1 ps time bins, while still maintaining the multi-stop and unlimited range capabilities of the coarse counter. Because such precise measurements require the circuitry to have recovered between event recordings, the HydraHarp enforces a dead time of around 80 ns [46]. Recent purely digital TDCs have much lower electronics dead times around 1 ns [47,48], although there is typically a trade-off between the dead time and minimum time bin duration. Still, the dead times of newer detector architectures [3,49,50] have similarly declined, so the relative dead-time durations remain important in determining detection time sequences.

## 3. OCCURRENCE TIME DISTRIBUTIONS

We now consider the timing of photons incident on the SPAD and the effects of the detector and electronics in determining which photons are detected. The steps from photon generation to detection are illustrated in Fig. 2.

#### A. Photoelectron Arrival Time Model

The light intensity incident on a detector [Fig. 2(a)] generates a sequence of photon times described by a Poisson process [Fig. 2(b)] [51,52]. We define those photons successfully converted to charge carriers as photoelectron *arrivals* [Fig. 2(c)], which are likewise described as a Poisson process since $\eta \lt 1$ simply results in a Bernoulli thinning of the arrivals [52]. The Poisson process intensity function $\lambda (t)$ is periodic and inhomogeneous for TCSPC due to the repeated illumination with period $t_{\text{r}}$. It follows from the properties of Poisson processes that the arrival times modulo $t_{\text{r}}$ are independent and identically distributed (i.i.d.) random variables on $[0,t_{\text{r}})$ with probability density function (PDF) proportional to $\lambda (t)$. In general, the intensity is composed of two parts: $\lambda (t) = {\lambda _{\text{s}}}(t) + {\lambda _{\text{b}}}(t)$, where ${\lambda _{\text{s}}}(t)$ is the time-varying intensity of a *signal* process, and ${\lambda _{\text{b}}}(t)$ is the intensity due to *background* (ambient light and dark counts), which is assumed to be a constant ${\lambda _{\text{b}}}$. For SPL with a single surface per laser position, ${\lambda _{\text{s}}}(t)$ is described parametrically in one cycle as the scaled and time-shifted illumination pulse ${\lambda _{\text{s}}}(t) = \alpha \beta \eta s(t - 2z/c)$, where $\alpha \in [0,1]$ is the target reflectivity, combining attenuation effects due to object reflectance, radial falloff, and view angle; $\beta \in [0,\infty)$ is a gain factor due to the optical illumination power; $s(t)$ is the illumination pulse shape normalized to be a valid PDF on $[0,t_{\text{r}})$; $z$ is the depth of the target; and $c$ is the speed of light. For other applications, $s(t)$ would be determined by other relevant parameters such as the fluorescence lifetime in FLIM or additional reflections for multidepth imaging. We refer to the quantities $S: = \int_0^{t_{\text{r}}} {\lambda _{\text{s}}}(t)\text{d}t$, $B: = {\lambda _{\text{b}}}t_{\text{r}}$, and $\Lambda : = \int_0^{t_{\text{r}}} \lambda (t)\text{d}t = S + B$ as the signal flux, background flux, and total flux, respectively, which represent the mean number of photoelectrons generated by each process within one illumination cycle. The flux could be considered equivalent to the detection rate if there were no dead times. The signal-to-background ratio (SBR) is $S/B$.

#### B. Detection Time PDF Derivation

The presence of detector and electronics dead times necessitates distinction between the arrival and detection processes. Not all photoelectrons that could cause a detector avalanche actually do, since some of them arrive during the SPAD quenching or hold-off period following a previous avalanche [Fig. 2(d)]. Likewise, not all avalanches are registered as detection events by the timing electronics because they occur during the dead time following a previous detection event [Fig. 2(e)]. Consequently, the detection process is not Poisson because the interdetection times are not independent. In the following, we account for the dependence due to dead times and show that the detection time sequence forms a Markov chain, i.e., given the most recent detection time, the PDF of the next detection time is independent of the previous detection times. We note that the detection time sequence being a Markov chain is specific to asynchronous systems operated in free-running mode. Synchronous systems [2,18,19] and those operating with range gates [23,26–30] or multiplexing detector elements [53–57] introduce additional factors that control the photon detection sequence. Still, we emphasize that, for any TCSPC system architecture, modeling the effect of dead times requires considering the entire photon detection sequence, even if parameter estimates are computed only from histograms of detection times.

### 1. Absolute Detection Times

We showed in previous work with free-running SPADs that if $t_{\text{d}} \ge t_{\text{e}}$, the sequence of detection times forms a Markov chain whose limiting PDF describes the empirical detection time histogram [24]. Although it is more complicated to account for both the detector and timing electronics, we derive in Supplement 1, Section 1 that the detection time sequence for $t_{\text{d}} \lt t_{\text{e}} \le t_{\text{r}} \le 2t_{\text{d}}$ is also a Markov chain. Note that the derivation can be extended to other $(t_{\text{d}},t_{\text{e}},t_{\text{r}})$ combinations with $t_{\text{d}} \lt t_{\text{e}}$; we omit the details for simplicity. Let ${\{{T_i}\} _{i \in {\mathbb N}}}$ denote *absolute* detection times measured from the start of the experiment. We show in Supplement 1, Section 1 that ${\{{T_i}\} _{i \in {\mathbb N}}}$ is a Markov chain with transition probability

### 2. Relative Detection Times

Let ${\{{X_i}\} _{i \in {\mathbb N}}}$ denote *relative* times measured with respect to the most recent illumination [e.g., used to form a histogram over $[0,t_{\text{r}})$]. Based on Eq. (1), we further show in Supplement 1, Section 1 that ${\{{X_i}\} _{i \in {\mathbb N}}}$ is also a Markov chain. The expression for the transition probability distribution changes slightly for different combinations of ${x_i}$ and ${x_{i + 1}}$ due to the different domains of the terms in Eq. (1) and how they interact with the modulo effect of the repetition period $t_{\text{r}}$. Specifically, for a particular detection time ${x_i}$ in one cycle, $[0,t_{\text{r}})$ is partitioned into regions related to the first possible cycle during which ${x_{i + 1}}$ could occur in each region; the transition probability also accounts for the possibility of ${x_{i + 1}}$ occurring in those regions during subsequent cycles. Details about the partition are illustrated in Supplement 1, Section 1, and we present an example for the case ${x_i} \in {\cal A} = [0,t_{\text{r}} - t_{\text{e}})$ in Fig. 3 to provide intuition about our derivation. As illustrated in Figs. 3(d)–3(g), we partition the space of ${x_{i + 1}}$ into

Since the Markov chain for the relative detection time sequence is irreducible, recurrent, and aperiodic, there will be a unique limiting distribution [58]. In other words, after enough photon detections, the marginal PDF of the next photon detection time will reach a steady state equal to the limiting distribution of the Markov chain. Thus, computing the limiting distribution gives us a useful approximation to the empirical PDF we get from forming a histogram. Although this computation is difficult to perform in closed form, we can compute an approximation as in [24] by discretizing the state space into $n_{\text{b}}$ bins, forming the Markov chain transition probability matrix ${\bf P}$ from Eq. (4) and computing the discretized limiting distribution ${{\bf f}^{\text{D}}} = {{\bf f}^{\text{D}}}{\bf P}$ as the leading left eigenvector of ${\bf P}$.

#### C. Photoelectron Flux Estimation

The total number of photoelectron arrivals is a Poisson random variable $N_{\text{tot}} \sim \text{Poisson}(n_{\text{r}}\Lambda)$, where $n_{\text{r}}$ is the number of illumination cycles. Similarly, the number of arrivals in a background-only calibration measurement (i.e., with the laser switched off) is ${N_{\text{b}}} \sim \text{Poisson}(n_{\text{r}}B)$. The maximum likelihood (ML) estimator of $\Lambda$ from the arrivals is ${\hat \Lambda ^{\text{AML}}} = N_{\text{tot}}/n_{\text{r}}$; likewise, ${\hat B^{\text{AML}}} = N_{{\text{b}}}/n_{\text{r}}$. However, the Poisson detection process approximation only holds at low flux, and errors from non-Poissonian statistics are incurred if dead-time effects are present. Photoelectron flux estimates that account for dead time typically assume a single source of dead time and constant flux [24,59,60]. Rate estimation at high flux requires not just the detection counts but their absolute time stamps as well. Let ${U_i} = {T_{i + 1}} - {T_i}$ for $i \ge 1$ be the *interdetection times*. Let the subsequence ${\{{U_{{i_k}}}\} _{{{k \ge 1}}}}$ of ${\{{U_i}\} _{i \ge 1}}$ contain all ${U_{{i_k}}}$s such that ${U_{{i_k}}} \gt t_{\text{d}} + t_{\text{e}}$. Isbaner *et al*. suggest in [31] that the “useful” interdetection periods, defined as $R_{k}=\lfloor (U_{i_{k}}-t_{e}-t_{d})/t_{r}\rfloor$, were distributed with an exponential decay that could be estimated with a weighted least-squares fit. Similar to our flux estimator in [24], we rigorously prove in Supplement 1, Section 2 that the ${R_k}$s are i.i.d. geometric random variables with probability mass function $P({R_k} = m) = {e^{- m\Lambda}}({1 - {e^{- \Lambda}}})$ for $m = 0,1,2, \ldots$. The detection rate ML estimate is then given as

## 4. DETECTION TIME MODEL VALIDATION

#### A. Simulated Data

We first validate that the Markov chain limiting distribution is a close approximation to the empirical distribution of detection times by simulating a sequence of detections for which the parameters are known exactly. For each of the $n_{\text{r}}$ illumination cycles, a sequence of incident photon times is first created by generating Poisson random numbers of signal and background detections, $N_{{\text{s}}} \sim \text{Poisson}(S)$ and ${N_{\text{b}}} \sim \text{Poisson}(B)$, and then drawing those numbers of samples from the respective temporal distributions. Next, the photoelectron arrival times are generated by selecting from the incident times with probability $\eta$; note that we can skip this step by assuming $\eta = 1$. Finally, starting with the first arrival, the sequence of detection times is determined to be those photoelectrons that arrive when neither the detector nor the electronics is dead.

Figure 4 shows histograms of data simulated with $n_{\text{r}} = 5 \times {10^5}$ illuminations, $t_{\text{r}} = 100\;\text{ns} $, $t_{\text{d}} = 50\;\text{ns} $, $t_{\text{e}} = 80\;\text{ns} $, pulse width ${\sigma _p} = 0.5\;\text{ns} $, bin duration ${t_{\text{bin}}} = 50\,\,\text{ps}$, and with the time of flight fixed at $2z/c = 75\;\text{ns} $. Overlaid on the histograms are the PDFs of arrival times from which the photoelectron times were generated; the PDFs predicted as the Markov chain limiting distributions including both $t_{\text{d}}$ and $t_{\text{e}}$; and the predicted PDFs using the model of [24], which only includes $t_{\text{d}}$. The $S$ and $B$ values sum to $\Lambda \gg 5\%$, so the dead times clearly cause distortions in the detection time histogram relative to the arrival time density. These distortions are accurately captured in our model, which predicts both the slight shift of the peak towards earlier detection times as well as the particular shape of the low-probability ripple behavior. Ignoring the effect of $t_{\text{e}}$ prevents the previous model in [24] from fully fitting the data, instead capturing only the dip following the signal pulse in the detection time PDF.

#### B. Experimental Data

We next evaluate the correctness of our model with experimental data from an actual TCSPC system. The SPL system implemented as in Fig. 1 uses a HydraHarp 400 TCSPC module (PicoQuant) with dead time $t_{\text{e}} \approx 80\;\text{ns} $ and a fast-gated SPAD detector module (Micro Photon Devices), which has an adjustable hold-off time between 48 ns and 1 µs. The illumination source, a pulsed diode laser (PicoQuant LDH-series) at 640 nm and with full width at half-maximum (FWHM) pulse duration around 100 ps, was aimed at a Lambertian white target at a fixed distance of around 50 cm. A band-pass filter (Thorlabs FB-640-10) at the operating wavelength was mounted in front of the SPAD to reduce the amount of incident ambient light. A distortion-free pulse shape calibration was acquired using a neutral-density (ND) filter with optical density (OD) 3.0, and high-flux measurements were acquired with no attenuation (OD 0), while the hold-off time ${t_{\text{ho}}}$ was varied between 48 ns and 198 ns. In each case, we assumed that $t_{\text{d}} \approx {t_{\text{ho}}} + 2\;\text{ns} $ and $t_{\text{e}} = 80\,\,\text{ns}$, and thus the dataset included cases of both $t_{\text{d}} \lt t_{\text{e}}$ and $t_{\text{d}} \gt t_{\text{e}}$. A high-flux dataset with the laser turned off was acquired to serve as a background calibration.

The total ${\hat \Lambda ^{\text{HF}}}$ and calibrated background ${\hat B^{\text{HF}}}$ flux values were estimated using the method proposed in Section 3.C. We attempted to perform the same validation as with simulated data using ${\hat S^{\text{HF}}} = {\hat \Lambda ^{\text{HF}}} - {\hat B^{\text{HF}}}$. However, the background calibration ${\hat B^{\text{HF}}}$ did not lead to a good fit between the measured histogram and the predicted effects of dead time for the calibrated low-flux pulse shape. Several factors not included in the acquisition model likely contributed to errors in the calibration:

- • Filter spectral response: The absorptive ND filters (ThorLabs NEK01) used in the experimental measurements have a fairly limited neutral region of 400–650 nm, so background outside of this passband is attenuated by varying amounts, and the relative attenuation of that background is not consistent for ND filters with different OD values.
- • Dark counts: Our modeling assumes all photons not due to signal can be grouped into a homogeneous “background” process, but attenuation of the incident flux only affects background counts from ambient light, with no effect on dark counts. While the dark count rate is often quite low for SPADs (on the order of 100 counts per second for our device), for large attenuation factors applied to the ambient light, extrapolation of the high-flux background intensity from a low-flux measurement is not accurately computed.
- • Afterpulsing: SPADs are typically held off long enough for the probability of an afterpulse detection to be sufficiently small and uncorrelated in time, so they simply appear as an increase in the background as the flux increases [61]. As we note in Supplement 1, Section 3, for shorter hold-off times, the afterpulse time correlation with respect to the most recent avalanche time is non-negligible, especially at high signal flux due to the periodicity of the signal component.

Instead of using poorly calibrated background values, we choose an approximate background value $\tilde B$ that minimizes the Kolmogorov–Smirnov (KS) statistic, a measure of similarity between the empirical and predicted cumulative distribution functions (CDFs) ([62], Chapter 10.2). For a histogram ${\bf h}$ of length $n_{\text b}=\lceil t_{r}/{t_{\text{bin}}}\rceil$ with $n$ total detections, we define the $j$th element of the empirical CDF as

Figure 5 shows the close match between the measured and predicted PDFs for a variety of ${t_{\text{ho}}}$ values. The insets in the plots highlight how even the shapes of lower-probability regions of the measured PDFs are predicted using our Markov chain analysis, including the multiple jumps for $t_{\text{d}} \lt t_{\text{e}}$ and the nearly distortion-free special case when $t_{\text{d}}$ is an integer multiple of $t_{\text{r}}$ as pointed out by [17]. The plots emphasize how the effect of dead times on the detection time distribution is sensitive to many parameters: $t_{\text{d}}$, $t_{\text{e}}$, $t_{\text{r}}$, $S$, and $B$. Additional acquisition data and flux estimates for each hold-off time are found in Supplement 1, Section 3.

## 5. APPROXIMATE HISTOGRAM CORRECTION

The equality that relates the arrival intensity (unknown) and the detection PDF (known up to some quantization and finite sample error) is the Markov chain stationary condition

where $p(y,x;\lambda) = {f_{{X_{i + 1}}|{X_i}}}(x|y)$ from Eq. (4), which is a function of the arrival intensity $\lambda$. Unfortunately, we do not obtain a simplified equation as in [24]. Although it is still possible to derive a gradient method directly from Eq. (9), computing the gradient is tedious. Therefore, we use the fixed-point iteration method for the estimation of $\lambda (x)$. Note from the expression of $p(y,x;\lambda)$ that we can pull out a $\lambda (x)$ factor from the integral; see the ${c_k}$ factors in Eq. (4). Then, by Eq. (9), we obtain the following fixed-point equation of $\lambda$:## 6. RESULTS

We demonstrate the merit of our dead-time compensation approach by forming point clouds combining depth and flux estimates from single-photon lidar. We adjust the experimental setup described in Section 4.B to have hold-off time ${t_{\text{ho}}} = 48\;\text{ns} $, laser power $\approx 4\;\text{mW} $, and illumination period $t_{\text{r}} = 100\;\text{ns} $. Detection times were recorded for $160 \times 240$ pixel raster scans of the indoor mannequin scene shown in Fig. 7(a). Long, low-flux measurements serving as a ground truth proxy were acquired for 0.1 s/pixel with OD 3.0 attenuation. High-flux acquisitions were made with an OD 0.1 ND filter.

Each depth estimate uses a pixelwise log-matched filter as described in [25]. Depth estimates for pixels with no detections are set to zero. Total flux estimates are displayed since the SBR is high and the signal flux is uncertain due to unreliable background estimation. Following the nomenclature of [24], the tested algorithms are as follows:

- • the conventional
**low-flux (LF)**method matches the calibrated pulse shape ${{\bf f}^{\text{A}}}$ to the low-flux histogram ${{\bf h}^{\text{LF}}}$ and uses the Poisson flux estimator ${\hat \Lambda ^{\text{AML}}}$; - • the naïve
**high-flux (HF)**method ignores dead-time distortions by matching ${{\bf f}^{\text{A}}}$ to the high-flux histogram ${{\bf h}^{\text{HF}}}$ and uses ${\hat \Lambda ^{\text{AML}}}$; - • our proposed
**Markov chain histogram correction**(MCHC) applies the proposed flux estimator ${\hat \Lambda ^{\text{DML}}}$ to estimate ${\hat \Lambda ^{\text{HF}}}$ and corrects ${{\bf h}^{\text{HF}}}$ via the method in Section 5 to match to ${{\bf f}^{\text{A}}}$; and - • our proposed
**Markov chain–based PDF prediction**(MCPDF) estimates ${\hat S^{\text{HF}}}$ and ${\hat B^{\text{HF}}}$ to compute the dead-time distorted distribution ${{\bf f}^{\text{D}}}$ to match ${{\bf h}^{\text{HF}}}$.

MCPDF implementation ideally requires precise estimates for $S$ and $B$ at each pixel to predict the dead-time-distorted PDF ${{\bf f}^{\text{D}}}$. However, due to the inaccuracies in background intensity calibration as described in Section 4.B and the fact that the background value changes for each pixel in a coaxial lidar system, obtaining accurate ${\hat B^{\text{HF}}}$ estimates for each pixel is a challenge. Fortunately, as we demonstrated in [25], a rough estimate of $S$ and $B$ may suffice, and the performance gain in terms of depth error is likely to be relatively small after a certain level of accuracy in the $S$ and $B$ estimates. We therefore follow a similar strategy here as proposed in [25], where we compute ${\hat \Lambda ^{\text{HF}}}$ for each pixel assuming a fixed SBR of 9 (i.e., ${\hat S^{\text{HF}}} = 0.9{\hat \Lambda ^{\text{HF}}}$, ${\hat B^{\text{HF}}} = 0.1{\hat \Lambda ^{\text{HF}}}$), and we choose the precomputed PDF ${{\bf f}^{\text{D}}}$ for the nearest of eight quantized levels of ${\hat \Lambda ^{\text{HF}}}$.

The MCHC method is more convenient than MCPDF, in the sense that it does not require background calibration. In fact, our histogram correction algorithm makes no assumptions on the parametric form of the transient information and would thus be broadly applicable to multidepth lidar, fluorescence lifetime estimation, etc. However, our current implementation of MCHC is quite slow and memory demanding, mainly due to the computation and storage of the matrix ${\bf Q}({{\boldsymbol \lambda}^t})$ at each iteration. We leave the study of efficient implementation of MCHC—possibly via general-purpose graphics processing units (GPUs)—for future work.

Full point cloud results in Fig. 7 show depths overlaid with flux estimates. Mean absolute error (MAE) values of the mannequin depth estimates are computed with respect to the ground truth estimate in Fig. 7(b). The results in Figs. 7(c) and 7(d) show typical performance for pixelwise log-matched filter depth estimators as described in [63]. For the low-flux acquisition in Fig. 7(c), few signal photons are detected, so the depth estimator is more sensitive to background photons, which can lead to large errors. The LF approach requires much longer acquisition times [e.g., $n_{\text{r}}{= 10^5}$, Fig. 7(d)] to detect enough signal photons and achieve reasonable accuracy in depth estimation. In the high-flux acquisitions in Figs. 7(e)–7(g), the maximum incident flux estimated at a pixel is around 5 photoelectrons per illumination cycle, which is 100 times greater than the detection rate limit recommended by the 5% rule of thumb. The unattenuated acquisition enables detection of a large number of signal photons in only $n_{\text{r}}{= 10^3}$ illuminations, yielding much smaller errors than the LF method for the same acquisition time, but also demonstrating the advantage of our Markov chain-based modeling for measurements at such high flux. The improvement of MCHC and MCPDF over HF corresponds to mitigating the small shift of the pulse peak due to dead time. A broader pulse would introduce a larger shift and even more distortion. The advantage of our Markov chain modeling is also clearly visible in the flux estimation, as conventional estimates saturate and fail to represent large flux values. A further demonstration of this effect and additional performance quantification can be found in Supplement 1, Section 5.

Comparing the proposed methods for this set of experiments and with our current implementation, MCHC is more straightforward to implement than MCPDF as we mentioned above, whereas the parametric signal model used for MCPDF enables a large improvement in speed and a small improvement in accuracy over MCHC. Specifically, to obtain the results shown here, MCHC required 41.6 h with 12 parallel cores on an Intel Xeon Gold 6132 processor and 192 GB of memory, whereas MCPDF required only 8.4 s on a laptop with a 3.3 GHz Intel Core i5 CPU.

## 7. CONCLUSION

This work demonstrates that the conventional wisdom regarding detection rates unduly limits acquisition speeds in TCSPC. While detector and electronics dead times prevent detection time histograms from being directly in proportion to the incident intensity at high flux, we show that a surprising degree of accuracy can be achieved by properly modeling the dead-time effects. Our probabilistic modeling of photon detection in the presence of detector and electronics dead times identified the photon detection time sequence as a Markov chain, which was validated by showing that the limiting distribution matches experimental histograms. We introduced a flux estimator and histogram correction method, which can be used to reconstruct the distribution of photon arrivals from a high-flux acquisition. Our modeling and estimation methods are generally applicable to other fields using similar TCSPC systems. For instance, our histogram correction method could be directly applied to FLIM, and methods of estimating single- or multi-exponential lifetimes could likewise be derived.

To illustrate our theoretical contributions, we demonstrate their application in ranging with single-photon lidar, notably achieving accurate depth imaging results with flux values up to two orders of magnitude greater than the standard 5% rule. Accurate estimation at high flux reduces the number of illuminations needed for ranging, thereby removing barriers to high-speed acquisition and enabling possible deployment of single-photon lidar in autonomous navigation and other real-time applications. Potential directions for future work include reliable background calibration, efficient implementation of our MCHC algorithm, and the incorporation of spatial regularization for even faster depth estimation from fewer photon detections.

## Funding

Charles Stark Draper Laboratory (Draper Fellowship); National Science Foundation (1422034, 1815896); Defense Advanced Research Projects Agency (HR0011-16-C-0030); Google.

## Acknowledgment

This work was completed while Yanting Ma was with Boston University. The authors thank F. N. C. Wong and C. Henley at MIT for their help with the optical setup and BU Research Computing Services for the computing resources.

## Disclosures

The authors declare no conflicts of interest.

See Supplement 1 for supporting content.

## REFERENCES

**1. **J. S. Massa, A. M. Wallace, G. S. Buller, S. J. Fancey, and A. C. Walker, “Laser depth measurement based on time-correlated single-photon counting,” Opt. Lett. **22**, 543–545 (1997). [CrossRef]

**2. **W. Becker, *Advanced Time-Correlated Single Photon Counting Techniques* (Springer, 2005).

**3. **R. H. Hadfield, “Single-photon detectors for optical quantum information applications,” Nat. Photonics **3**, 696–705 (2009). [CrossRef]

**4. **Y. Altmann, S. McLaughlin, M. J. Padgett, V. K. Goyal, A. O. Hero, and D. Faccio, “Quantum-inspired computational imaging,” Science **361**, eaat2298 (2018). [CrossRef]

**5. **C. Bruschini, H. Homulle, I. M. Antolovic, S. Burri, and E. Charbon, “Single-photon avalanche diode imagers in biophotonics: review and outlook,” Light Sci. Appl. **8**, 2047–7538 (2019). [CrossRef]

**6. **D. Faccio, A. Velten, and G. Wetzstein, “Non-line-of-sight imaging,” Nat. Rev. Phys. **2**, 318–327 (2020). [CrossRef]

**7. **J. Rapp, J. Tachella, Y. Altmann, S. McLaughlin, and V. K. Goyal, “Advances in single-photon lidar for autonomous vehicles: working principles, challenges, and recent advances,” IEEE Signal Process. Mag. **37**(4), 62–71 (2020). [CrossRef]

**8. **A. Kirmani, D. Venkatraman, D. Shin, A. Colaço, F. N. C. Wong, J. H. Shapiro, and V. K. Goyal, “First-photon imaging,” Science **343**, 58–61 (2014). [CrossRef]

**9. **D. Shin, A. Kirmani, V. K. Goyal, and J. H. Shapiro, “Photon-efficient computational 3-D and reflectivity imaging with single-photon detectors,” IEEE Trans. Comput. Imaging **1**, 112–125 (2015). [CrossRef]

**10. **Y. Altmann, X. Ren, A. McCarthy, G. S. Buller, and S. McLaughlin, “Lidar waveform-based analysis of depth images constructed using sparse single-photon data,” IEEE Trans. Image Process. **25**, 1935–1946 (2016). [CrossRef]

**11. **A. M. Pawlikowska, A. Halimi, R. A. Lamb, and G. S. Buller, “Single-photon three-dimensional imaging at up to 10 kilometers range,” Opt. Express **25**, 11919–11931 (2017). [CrossRef]

**12. **Z.-P. Li, X. Huang, P.-Y. Jiang, Y. Hong, C. Yu, Y. Cao, J. Zhang, F. Xu, and J.-W. Pan, “Super-resolution single-photon imaging at 8.2 kilometers,” Opt. Express **28**, 4076–4087 (2020). [CrossRef]

**13. **Z.-P. Li, X. Huang, Y. Cao, B. Wang, Y.-H. Li, W. Jin, C. Yu, J. Zhang, Q. Zhang, C.-Z. Peng, F. Xu, and J.-W. Pan, “Single-photon computational 3D imaging at 45 km,” Photon. Res. **8**, 1532–1540 (2020). [CrossRef]

**14. **R. A. Barton-Grimley, R. A. Stillwell, and J. P. Thayer, “High resolution photon time-tagging lidar for atmospheric point cloud generation,” Opt. Express **26**, 26030–26044 (2018). [CrossRef]

**15. **M. Patting, M. Wahl, P. Kapusta, and R. Erdmann, “Dead-time effects in TCSPC data analysis,” Proc. SPIE **6583**, 658307 (2007). [CrossRef]

**16. **D. V. O’Connor and D. Phillips, *Time-Correlated Single Photon Counting* (Academic, 1984).

**17. **A. Cominelli, G. Acconcia, P. Peronio, M. Ghioni, and I. Rech, “High-speed and low-distortion solution for time-correlated single photon counting measurements: a theoretical analysis,” Rev. Sci. Instrum. **88**, 123701 (2017). [CrossRef]

**18. **P. B. Coates, “The correction for photon ‘pile-up’ in the measurement of radiative lifetimes,” J. Phys. E **1**, 878–879 (1968). [CrossRef]

**19. **A. Gupta, A. Ingle, A. Velten, and M. Gupta, “Photon-flooded single-photon 3D cameras,” in *Proc. IEEE Conf. Comput. Vis. Pattern Recog.* (2019), pp. 6770–6779.

**20. **A. Gupta, A. Ingle, and M. Gupta, “Asynchronous single-photon 3D imaging,” in *Proc. IEEE Int. Conf. Comput. Vis.* (2019), pp. 7908–7917.

**21. **F. Heide, S. Diamond, D. B. Lindell, and G. Wetzstein, “Sub-picosecond photon-efficient 3D imaging using single-photon sensors,” Sci. Rep. **8**, 17726 (2018). [CrossRef]

**22. **R. A. Barton-Grimley, J. P. Thayer, and M. Hayman, “Nonlinear target count rate estimation in single-photon lidar due to first photon bias,” Opt. Lett. **44**, 1249–1252 (2019). [CrossRef]

**23. **A. K. Pediredla, A. C. Sankaranarayanan, M. Buttafava, A. Tosi, and A. Veeraraghavan, “Signal processing based pile-up compensation for gated single-photon avalanche diodes,” arXiv:1806.07437 (2018).

**24. **J. Rapp, Y. Ma, R. M. A. Dawson, and V. K. Goyal, “Dead time compensation for high-flux ranging,” IEEE Trans. Signal Process. **67**, 3471–3486 (2019). [CrossRef]

**25. **J. Rapp, Y. Ma, R. M. A. Dawson, and V. K. Goyal, “Dead time compensation for high-flux depth imaging,” in *Proc. IEEE Int. Conf. Acoust., Speech, and Signal Process.* (IEEE, 2019), pp. 7805–7809.

**26. **S. Johnson, P. Gatt, and T. Nichols, “Analysis of Geiger-mode APD laser radars,” Proc. SPIE **5086**, 486384 (2003). [CrossRef]

**27. **P. Gatt, S. Johnson, and T. Nichols, “Dead-time effects on Geiger-mode APD performance,” Proc. SPIE **6550**, 65500I (2007). [CrossRef]

**28. **P. Gatt, S. Johnson, and T. Nichols, “Geiger-mode avalanche photodiode ladar receiver performance characteristics and detection statistics,” Appl. Opt. **48**, 3261–3276 (2009). [CrossRef]

**29. **P. Zhao, Y. Zhang, W. P. Qian, and Y. Xuan, “Investigation of Geiger-mode detector in multi-hit model for laser ranging,” Sci. China Technol. Sci. **58**, 943–950 (2015). [CrossRef]

**30. **Z. Li, J. Lai, C. Wang, W. Yan, and Z. Li, “Influence of dead-time on detection efficiency and range performance of photon-counting laser radar that uses a Geiger-mode avalanche photodiode,” Appl. Opt. **56**, 6680–6687 (2017). [CrossRef]

**31. **S. Isbaner, N. Karedla, D. Ruhlandt, S. C. Stein, A. Chizhik, I. Gregor, and J. Enderlein, “Dead-time correction of fluorescence lifetime measurements and fluorescence lifetime imaging,” Opt. Express **24**, 9429–9445 (2016). [CrossRef]

**32. **M. D. Eisaman, J. Fan, A. Migdall, and S. V. Polyakov, “Invited review article: single-photon sources and detectors,” Rev. Sci. Instrum. **82**, 071101 (2011). [CrossRef]

**33. **M. J. Stevens, “Photon statistics, measurements, and measurements tools,” in *Single-Photon Generation and Detection*, A. Migdall, S. V. Polyakov, J. Fan, and J. C. Bienfang, eds. (Academic, 2013), chap. 2, pp. 25–68.

**34. **C. J. Chunnilall, I. P. Degiovanni, S. Kück, I. Müller, and A. G. Sinclair, “Metrology of single-photon sources and detectors: a review,” Opt. Eng. **53**, 081910 (2014). [CrossRef]

**35. **M. Ghioni, S. Cova, F. Zappa, and C. Samori, “Compact active quenching circuit for fast photon counting with avalanche photodiodes,” Rev. Sci. Instrum. **67**, 3440–3448 (1996). [CrossRef]

**36. **A. C. Giudice, M. Ghioni, S. Cova, and F. Zappa, “A process and deep level evaluation tool: afterpulsing in avalanche junctions,” in *Eur. Solid-State Device Res. Conf.* (IEEE, 2003), pp. 347–350.

**37. **M. A. Itzler, R. Ben-Michael, C. F. Hsu, K. Slomkowski, A. Tosi, S. Cova, F. Zappa, and R. Ispasoiu, “Single photon avalanche diodes (SPADs) for 1.5 µm photon counting applications,” J. Mod. Opt. **54**, 283–304 (2007). [CrossRef]

**38. **Y. Xu, P. Xiang, and X. Xie, “Comprehensive understanding of dark count mechanisms of single-photon avalanche diodes fabricated in deep sub-micron CMOS technologies,” Solid-State Electron. **129**, 168–174 (2017). [CrossRef]

**39. **S. Cova, A. Lacaita, and G. Ripamonti, “Trapping phenomena in avalanche photodiodes on nanosecond scale,” IEEE Electron Device Lett. **12**, 685–687 (1991). [CrossRef]

**40. **M. A. Itzler, X. Jiang, and M. Entwistle, “Power law temporal dependence of InGaAs/InP SPAD afterpulsing,” J. Mod. Opt. **59**, 1472–1480 (2012). [CrossRef]

**41. **L. M. Bollinger and G. E. Thomas, “Measurement of the time dependence of scintillation intensity by a delayed-coincidence method,” Rev. Sci. Instrum. **32**, 1044–1050 (1961). [CrossRef]

**42. **G. W. Roberts and M. Ali-Bakhshian, “A brief introduction to time-to-digital and digital-to-time converters,” IEEE Trans. Circuits Syst. II, Exp. Briefs **57**, 153–157 (2010). [CrossRef]

**43. **J. Kalisz, “Review of methods for time interval measurements with picosecond resolution,” Metrologia **41**, 17–32 (2004). [CrossRef]

**44. **M. Wahl, H. J. Rahn, T. Röhlicke, G. Kell, D. Nettels, F. Hillger, B. Schuler, and R. Erdmann, “Scalable time-correlated photon counting system with multiple independent input channels,” Rev. Sci. Instrum. **79**, 123113 (2008). [CrossRef]

**45. **M. Wahl, “Modern TCSPC electronics: Principles and acquisition modes,” in *Advanced Photon Counting: Applications, Methods, Instrumentation*, P. Kapusta, M. Wahl, and R. Erdmann, eds. (Springer, 2015), pp. 1–21.

**46. **M. Wahl, personal communication (2019).

**47. **H. Fedder, S. Oesterwind, M. Wick, F. Olbrich, P. Michler, T. Veigel, M. Berroth, and M. Schlagmüller, “Characterization of electro-optical devices with low jitter single photon detectors—towards an optical sampling oscilloscope beyond 100 GHz,” in *Proc. Eur. Conf. on Opt. Commun.* (IEEE, 2018), pp. 1–3.

**48. **M. Wahl, T. Röhlicke, S. Kulisch, S. Rohilla, B. Krämer, and A. C. Hocke, “Photon arrival time tagging with many channels, sub-nanosecond deadtime, very high throughput, and fiber optic remote synchronization,” Rev. Sci. Instrum. **91**, 013108 (2020). [CrossRef]

**49. **X. Michalet, A. Cheng, J. Antelman, M. Suyama, K. Arisaka, and S. Weiss, “Hybrid photodetector for single-molecule spectroscopy and microscopy,” Proc. SPIE **6862**, 68620F (2008). [CrossRef]

**50. **J. Huang, W. Zhang, L. You, C. Zhang, C. Lv, Y. Wang, X. Liu, H. Li, and Z. Wang, “High speed superconducting nanowire single-photon detector with nine interleaved nanowires,” Supercond. Sci. Technol. **31**, 074001 (2018). [CrossRef]

**51. **B. Saleh, *Photoelectron Statistics: With Applications to Spectroscopy and Optical Communications*, Vol. 6, of Springer Series in Optical Sciences (Springer-Verlang, 1978).

**52. **D. L. Snyder and M. I. Miller, *Random Point Processes in Time and Space* (Springer, 2012).

**53. **J. Arlt, D. Tyndall, B. R. Rae, D. D.-U. Li, J. A. Richardson, and R. K. Henderson, “A study of pile-up in integrated time-correlated single photon counting systems,” Rev. Sci. Instrum. **84**, 103105 (2013). [CrossRef]

**54. **C. Zhang, S. Lindner, I. M. Antolovic, J. M. Pavia, M. Wolf, and E. Charbon, “A 30-frames/s, 252 × 144 SPAD flash LiDAR with 1728 dual-clock 48.8-ps TDCs, and pixel-wise integrated histogramming,” IEEE J. Solid-State Circuits **54**, 1137–1151 (2019). [CrossRef]

**55. **S. W. Hutchings, N. Johnston, I. Gyongy, T. Al Abbas, N. A. Dutton, M. Tyler, S. Chan, J. Leach, and R. K. Henderson, “A reconfigurable 3-D-stacked SPAD imager with in-pixel histogramming for flash LIDAR or high-speed time-of-flight imaging,” IEEE J. Solid-State Circuits **54**, 2947–2956 (2019). [CrossRef]

**56. **A. R. Ximenes, P. Padmanabhan, M.-J. Lee, Y. Yamashita, D.-N. Yaung, and E. Charbon, “A modular, direct time-of-flight depth sensor in 45/65-nm 3-D-stacked CMOS technology,” IEEE J. Solid-State Circuits **54**, 3203–3214 (2019). [CrossRef]

**57. **P. Padmanabhan, C. Zhang, and E. Charbon, “Modeling and analysis of a direct time-of-flight sensor architecture for LiDAR applications,” Sensors **19**, 5464 (2019). [CrossRef]

**58. **S. P. Meyn and R. L. Tweedie, *Markov Chains and Stochastic Stability* (Springer, 2012).

**59. **I. M. Antolovic, S. Burri, C. Bruschini, R. Hoebe, and E. Charbon, “Nonuniformity analysis of a 65-kpixel CMOS SPAD imager,” IEEE Trans. Electron Devices **63**, 57–64 (2016). [CrossRef]

**60. **A. Ingle, A. Velten, and M. Gupta, “High flux passive imaging with single-photon sensors,” in *Proc. IEEE Conf. Comput. Vis. Pattern Recog.* (2019), pp. 6753–6762.

**61. **G. Humer, M. Peev, C. Schaeff, S. Ramelow, M. Stipčević, and R. Ursin, “A simple and robust method for estimating after pulsing in single photon detectors,” J. Lightwave Technol. **33**, 3098–3107 (2015). [CrossRef]

**62. **A. M. F. Mood, F. A. Graybill, and D. C. Boes, *Introduction to the Theory of Statistics*, 3rd ed. (McGraw-Hill, 1974).

**63. **J. Rapp and V. K. Goyal, “A few photons among many: unmixing signal and noise for photon-efficient active imaging,” IEEE Trans. Comput. Imaging **3**, 445–459 (2017). [CrossRef]