Photostatistics Reconstruction via Loop Detector Signatures

Photon-number resolving detectors are a fundamental building-block of optical quantum information processing protocols. A loop detector, combined with appropriate statistical processing, can be used to convert a binary on/off photon counter into a photon-number-resolving detector. Here we describe the idea of a signature of photon-counts, which may be used to more robustly reconstruct the photon number distribution of a quantum state. The methodology is applied experimentally in a 9-port loop detector operating at a telecommunications wavelength and compared directly to the approach whereby only the number of photon-counts is used to reconstruct the input distribution. The signature approach is shown to be more robust against calibration errors, exhibit reduced statistical uncertainty, and reduced reliance on a-priori assumptions about the input state.


I. INTRODUCTION
One of the most promising applications of quantum science lies in the field of quantum information processing [1]. Quantum optics provides an extremely convenient test bed for quantum information ideas because quantum effects dominate even at room temperatures and non-classical states are relatively easy to achieve in the laboratory. This has been most impressively illustrated with recent small-scale experimental demonstrations of quantum algorithms [2,3].
Photon-number resolving detectors, as well as the ability to characterise quantum optical states using such detectors, are a key element in the further development of quantum optical information processing systems [4,5]. Due to the low energies carried by photons, contemporary detectors are incapable of resolving photon number with high fidelity. Typically linear PIN photodiodes used as direct detectors of the optical field are limited to mesoscopic photonnumber regimes [6]. Recent results show promise at telecommunications wavelengths using liquid-He cooled diodes [7], bolometric devices [8] and super-conducting nano-structures [9] although such solutions are not yet widespread due to their inherent cost and/or cryogenic cooling complexities.
Consequently, use is commonly made of single photon detector modules (SPDMs), based on avalanche photo-diodes (APDs). A number of techniques have been developed for approximating number-resolving detection using these inherently non-number-resolving detectors [10,11,12,13,14,15,16,17,18,19], all of which are variations on the idea that the optical field is distributed across multiple modes which are measured independently. The experiment is arranged so that the probability of any given mode being populated by more than one photon approaches zero. Thus in the absence of dark counts, afterpulsing or imperfect quantum efficiency, the sum of the detection events across the modes closely approximates the number of photons in the incident state.
In the case of the ideal N -port multiplexer with ideal binary photon detectors, the maximum number of photons able to be resolved is equal to the number of output modes. This limit arises due to excess input photons resulting in > 1 photons in the output modes, which are not distinguishable from the single photon case. In a realistic detector, however, the presence of loss means that n input photons may result in m clicks where n ≥ m. If m ≤ N and n > N then it can be seen that the loss gives rise to a useable increase in the detector's dynamic range. In general though, the accuracy of the reconstruction degrades and the statistical errors increase with increasing amounts of loss. Similarly, fragile non-classical effects are washed out by increased losses, as quantified in [22]. This situation is further complicated by the presence of dark counts and afterpulsing.
There have been many demonstrated solutions to the reconstruction of the photostatistics (photon number distribution) of an input field using binary (on/off) detectors, but they may each be loosely grouped into one of three categories. All realistic photodetectors exhibit loss, which is modelled by preceding an otherwise ideal detector with a beamsplitter of transmission η. The photostatistics of the input state are thus convolved with the Bernoulli (binomial) distribution arising from the stochastic photon selection process performed at the beamsplitter [23].
This relationship forms the basis of the first category of reconstruction techniques -those involving inversion of the Bernoulli distribution. The analytic inversion is cumbersome and requires η ≥ 0.5 [24], an impossible constraint on N-ports with N > 2 which is further exaggerated at telecommunications wavelengths due to the very small quantum efficiencies. Numerical approaches such as maximum likelihood techniques permit more realistic values of η and have been successfully applied in [25,26].
The second category involves reconstruction of the input density matrix, ρ, via conditional probabilities and appears the most common contemporary approach. An N-port detector may be completely described by the corresponding conditional probability matrix P (m|n) which quantifies the probability of observing m clicks, given n incident photons. The experiments of [14,16,27] perform a large number of measurements of an unknown input state, and record the number of clicks obtained from each measurement. By applying either the inverted matrix of P (m|n) values to the observed "number of click" statistics or utilising maximum likelihood techniques, the input photon number distribution may be obtained.
Finally, as observed by [28], an improvement in the accuracy of the reconstruction is possible by noting the coincidences of the m detection events of the N detectors. This is equivalent to recording detection signatures for each measurement and forms the basis of our reconstruction implementation. To date, experimental work in this third category has considered only N = 2 distinct output modes [29] for a number of detection efficiencies η. The purpose of this paper is to present experimental results for N = 9 time-domain modes with a fixed (small) value of η for each mode.

II. NUMERICAL METHODS
Irrespective of the reconstruction category, the three approaches share the same goal common to all tomography problems -the calculation of the input state from a large data set of marginal distributions. Although the analytic relationship between the measured marginal probability distributions is usually a known function of the input state, this is not the case for the inverse in general. Even if an inverse analytic solution exists it is generally of no use as experimental and measurement uncertainties create an overcomplete data set for which there is no unique solution for ρ. The problem may thus be restated as the determination of the best estimateρ for ρ given the measured data set.
The estimation ofρ is significantly simplified by a priori assumptions of aspects of the statistics of the input state. If it may be assumed that the input state consists of negligible contributions from higher order Fock states, i.e. k|ρ|k ≈ 0 ∀ k ≥ K, then the reconstructed Hilbert space may be restricted to dimension (K +1)×(K +1). Additional constraints such as the expected distribution (e.g. Poissonian) and physicality via Tr[ρ] = 1 andρ kk ≥ 0 ∀ k < K further aid the reconstruction process.
Detection with on/off photodetectors will result in a particular pattern of detector clicks, which we call a signature. The probability, P sig ( d), of the detection signature occurring whereby each of the detectors in vector d trigger, and all others do not is given by [30] where p c (i) represents the probability of an input photon being coupled out of the i th port. Similarly, p loss (j) indicates the probability of photon loss associated with port j. Both of these parameters are a function of the detector architecture and are detailed in [30]. The dark count probability p dc depends upon the particular on/off detector employed. The net probability of detecting m photons, given by summing over all combinations of detection signatures where | d| = m, for | d| is total number of clicks, is thus It is worth noting that P sig ( d) is not necessarily the same for all d where | d| = m. That is, different signatures corresponding to the same measured number of photons needn't have equal probabilities of occurring.
For the three scenarios above, the measured data set is reduced to a vector of observed probabilities for a detection in the relevant basis. In the binomial reconstruction approach, this corresponds to the probability of a clickp click for each of the experimental values of detection efficiency η x . The probability p out (m) of obtaining m output photons depends on the m th and subsequent terms of the input state density matrix ρ. i.e.
Incorporating the independent dark count statistics and assuming a binary on/off detector (typically an SPDM) with a linear response to photon number, Eq. 3 may be integrated over the truncated Hilbert space to give Corresponding expressions can be obtained for the photon counting and signature approaches from Eqs. 1 and 2 respectively, i.e. where As stated earlier, these expressions allow the calculation of the expected statistics for a known density matrix. To reconstruct the unknown state, the statistical distance between the predicted and measured data sets is quantified. A correction is thus determined and applied to the current state prediction to create a new estimate. This process is repeated until either the desired degree of convergence between the datasets is achieved or the incremental changes to the estimated state fall below an arbitrary threshold.
The determination of the optimal means of quantifying the statistical distance between the data sets is an area of ongoing research. Similarly, the calculation of the corresponding amendment to the estimated state is also a nontrivial exercise. A comparison of the popular contemporary methods of maximum likelihood and maximum entropy approaches are discussed in the context of the inversion of photostatistics data in [31].
For our purposes, we quantify the distance between the theoretical and measured data sets via the mean square error between the distributions used as a measure of convergence by [32]. If p andp are the calculated and observed marginal probability distributions in the basis applicable to the reconstruction category, and S is the total data set size, The problem may hence be expressed as a multidimensional numerical minimisation problem, where the goal is the minimisation of ǫ by appropriate selection of the (K + 1) diagonal coefficients of ρ. For this task, the simplex method of Nelder and Mead [33] was chosen for ease of implementation [34] and the requirement for evaluation of ǫ only, and not its derivatives. Although convergence of the algorithm to a global minimum is not guaranteed [35], in practice its convergence was found to be fast and entirely adequate for a wide range of expressions for ǫ.

III. 1550 NM LOOP DETECTOR
The TDM loop detector architecture of [15] was selected because of its inherent flexibility in permitting a variable number N of output modes. The architecture is shown in Fig. 1. The experiment allows the various reconstruction methods to be assessed against the Poissonian photostatistics of the coherent output of the diode laser. The laser is biased below threshold, and pulsed briefly into operation by the application of a variable amplitude 20 ns pulse. 95% of the laser output is directed to a NIST-traceable Newport 3227 power meter for calibration/renormalisation purposes. The remaining 5% is then transferred via a series of calibrated fibre attenuators providing ∼ 70 dB of attenuation into the loop detector.
The loop detector comprises the input switch, output coupler and fibre delay loop. The switch was provided by Civcom Inc. and provides a remarkably low insertion loss (1.2 dB) combined with fast (270 ns) switching speed. As the fibre provides 10.044 µs of delay, the switch thus has adequate time to settle in the alternative state (shown by the dotted line of Fig. 1) after admitting the pulse to the detector. The loop consists of 2 km of Corning SMF-28 fibre, exhibiting only 0.8 dB loss at 1550 nm, with each coupler exhibiting 0.5dB of loss. The long delay allows for the recovery of the APD used in the SPDM following a gate pulse, in order to minimise the detrimental effects of afterpulsing. A portion of the beam is tapped off from the circulating field via the fibre coupler at each pass and directed onto an id-Quantique id-200 SPDM operating with a 20 ns gate window and a 100 kHz clock rate. This detector had a quantum efficiency of 10 % and a dark count probability per gate p dc = 9.6 × 10 −4 . N = 14 bins was selected with a coupling coefficient p c = 0.1, yielding a possible truncation error of 23% [30]. In practice, this error is not observed due to the excess loop losses preventing photons making it into the last few bins for the range of input intensities employed. Combined with the η d = 10% quantum efficiency of the id-Quantique id-200 SPDM used, the equivalent detection efficiencies η i for each time bin are shown in Table I. Due to the long delay employed, the loop path length is much greater than the coherence length of the laser. Consequently, the tapped field is uncorrelated with the circulating component within the loop, and any quantum correlations between successive time bins are destroyed. This allows the loop data be considered mathematically classical and photons to be treated as discrete particles.
Accurate timing of the experiment is critical to ensure the APD gate window that defines each time bin is correctly aligned with the pulse circulating within the loop. Control of the timing is achieved via a synchronous finite state machine (FSM) constructed within a field programmable gate array (FPGA). A measurement cycle consisting of the generation of an optical pulse, recording 14 passes of the loop and transferring the data to file for post processing requires 376µs, allowing 2660 distinct measurements per second to be performed. a) 14 15 16 17 18 19 20 21 22 23 24  Compensation for the multitude of electro-optic delays within the experiment is achieved by adjusting the triggering time of the SPDM and the clock frequency of the FSM. Fig. 2 indicates the frequency of detection events in the first (and subsequent) time bins as a function of time delay and clock frequency respectively. The final timing parameters chosen thus allow the detection efficiencies to approach those listed in Table I. Fig. 3 illustrates the resulting exponentially decaying click probabilities obtained for each bin. In practice, it was noted that very few valid detection events were observed in the later time bins. This is as expected due to the vanishingly small effective detection efficiencies. Consequently, the analysis was terminated at N = 9 time bins. To compensate for the induced truncation error, a non-experimentally observable 10 th bin was re-introduced to the expression of Eq. 6 to serve as a "catch all" for photons normally lost to the analysis.
To accurately model the 9 time-bin experimental data, Eq. 6 is thus rewritten and the 10 th element of the coupling and loss vectors are such that p c (10) = 1 − 9 i=1 p c (i) and p loss (10) = 0. As a test of validity, this ensures

IV. CHARACTERISATION
In addition to the bin efficiencies listed in Table I, our detector can be completely characterised by its conditional probability matrix. The first few terms are listed in Table II and plotted graphically in Fig. 4.  Immediately apparent is the departure from the ideal P (m|n) = δ m,n , arising as a direct consequence of the loop losses and the poor InGaAs APD efficiency. The quantum efficiency of the id-200 can be doubled at the cost of an order of magnitude increase in dark counts, i.e. at η d = 10%, p dc = 9.66 × 10 −4 , while at η d = 25%, p dc = 2.75 × 10 −2 . However, dark counts adversely affect fidelity to a greater extent than reduced quantum efficiency [30]. Consequently, the lower of the id-200 detection efficiencies were used.
The remnant dark count "noise floor" is visible at the exponential tail of Fig. 3 count probabilities. Fig. 5 provides further insight into the limitations of the photodetection process, illustrating the effect of disconnecting the loop at point X as shown in Fig. 1 for η d = 10% and η d = 25% respectively. The non-zero detection probabilities in bins ≥ 2 are the result of dark counts and afterpulsing. Afterpulsing describes the mechanism by which charges trapped within defects of the APD semiconductor give rise to spontaneous avalanche events during the next gate pulse, even in the absence of a detected photon. Unfortunately, since avalanches induced by a valid photon detection event are indistinguishable from APD defect/intrinsic dark count avalanche events, afterpulsing can artificially inflate the effective dark noise floor of many subsequent detector time bins. For fixed APD temperature and gating period, the afterpulsing probability scales non-linearly with η d [36]. This effect is clearly shown in Fig. 5(b), where photons are detected only during time bin 1 and afterpulsing is responsible for the exponential decay to the dark noise background, which is substantially greater than the η d = 10% case of Fig.  5(a). In order to quantify the effects of afterpulsing, the asymptotic background value dark noise was assumed to be independent of detection events and was subtracted away from the open loop data. Fig. 6 indicates the result of comparing the probability of a click in each bin to the one before for the data of Fig. 5(b). A significant 60% of detection events in subsequent time bins are thus observed to result solely from afterpulsing effects. The approximately constant envelope of probabilities suggests the simple "regenerative afterpulsing" model described in [37] adequately models our experiment. Assuming p dc ≪ 1, the detected p click data in bin i may thus be corrected to first order by subtracting a portion of the counts p a detected in a prior time bin, i.e.
Note that this correction aids reconstruction via the binomial method only, as integration of the comprehensive afterpulsing model of [37] into the expressions of Eqs. 8 and 5 is non-trivial.   7 shows the corrections of Eq. IV applied to a typical p click data set. Note that bin 2 consistently experiences the greatest correction as a consequence of afterpulsing following the significant number of detection events in the first time bin. With η d = 10%, p a was measured at a modest ∼ 3%. Given our inability to presently account for afterpulsing in the signature based reconstruction analysis, η d = 10% was chosen for all subsequent experimental data runs.

V. ANALYSIS OF RESULTS
In order to evaluate the proposed signature-based reconstruction method, a number of coherent state input tests were performed. 1.5 × 10 5 fixed amplitude pulses were acquired by the loop detector for each of a range of calibrated input intensities. The data was then processed to form two probability distributions -p click (i) for time bins (and thus effective η) i = 1 . . . 10 andp sig ( d) N =9 for all 2 N = 512 signatures.
The photostatistics of the input states were assumed to be Poissonian. For each input amplitude, the MSE given by Eq. 7 was calculated from the experimental data set and a range of predicted data sets generated by Eq. 4 and Eq. 6 for the binomial and signature schemes respectively. To create the predicted data sets, a spectrum of estimated density matricies were generated with Poisson-distributed photon numbers wheren was swept over the interval [0, ∼ 15] photons/pulse. Each of the evaluated input states were renomalised for a unity trace to account for the truncated Hilbert space, with the diagonal coefficients thus given by By means of illustration, Fig. 8 shows the predicted magnitudes of the first 30p sig probabilities for an input coherent state corresponding ton = 8 photons/pulse. These probability values were then compared to the measured experimental probabilities via Eq. 7 to determine the associated MSE. The indicated binary signature vectors correspond to the base-10 signature index, stored in bin-reversed order, i.e. time bin 1 is the right most bit [39]. For example, the signature vector corresponding to a detection in the first and third time-bins is d = {0, 0, 0, 0, 0, 0, 1, 0, 1}. This may be represented as a base-2 (binary) number, i.e. 000000101 2 , which has an equivalent base-10 (decimal) signature index 5 10 . Signature index 6 10 thus corresponds to m = 2 detection events -one in each of the second and third time  bins and so on. As expected from the P (m|n) data of Section IV, the probability of multiphoton detection events decreases rapidly with increasing m, hence the tiny value shown for index 23 10 = 000010111 2 , with m = 4. Fig. 9 illustrates the calculated values for ǫ against mean photon-numbern for the experimental binomialp click and p sig signature probabilities for an input state ofn = 6.5 ± 0.8 photons/pulse. A third trace illustrating the results of considering thep click data corrected for afterpulsing as per Eq. IV is also shown. Note that the latter trace exhibits a very much smaller minimum value of ǫ. This suggests a much closer agreement between the modelled and detected data, and the hence the validity of the applied correction. The minimum value of ǫ for a given reconstruction scheme indicates the value ofn that results in optimum agreement between the predicted and experimentally observed detection probabilities. The offsets between the minima of the three curves result from the inherent sensitivity of the binomial approach to the absolute values of η used in Eq. 4. While approaches to self-calibration have been proposed to mitigate these effects [6,38], the signature reconstruction method appears relatively insensitive to variations of individual values of η. For example, 5% variations of the switch loss term t s shifted the signature minimum by a corresponding amount, yet induced > 10% shifts in the binomialp click minimum. Note also that the signature ǫ value is a more sharply peaked function, potentially improving the accuracy of photon number estimation.
The significance of the departures of the estimated values ofn from the actual input state is examined in Fig. 10. binomial and signature approaches of Fig. 9. For each case, the input distribution is shown for comparison. The error bars were determined from the calibrated experimental uncertainty of the input amplitude. The severe consequences of the calibration-induced errors of the binomial scheme are starkly visible, with the estimated distribution significantly different from the input state. Conversely, the estimated signature statistics closely agree with the actual input statistics.   11 indicates the inferred photon numbers corresponding to the minima of the MSE ǫ error term for the two approaches over a range of input amplitudes. The ideal 1:1 relationship is also shown. The signature approach consistently returns values within the error bounds of loop detector losses and laser intensity calibration. The linearly increasing error corresponds to a systematic calibration error of ±15% or 0.6dB, readily accounted for in the multiple insertion loss uncertainties arising from fibre-fibre interconnections.
Finally, we investigate the performance of the signature reconstruction approach when no a priori assumptions are made regarding the form of the input photostatistics. In this case, we seek to find the density matrix ρ that minimises the MSE error ǫ by searching (K + 1)-dimension parameter space for an optimum set of diagonal coefficients. While we make use of the same coherent state experimental data used in the analysis above, the diagonal photon-number terms of the reconstructed density matrix are considered independent from one another and not bound to assume Poissonian form. To this end, we employed the numerical multidimensional simplex minimisation algorithm described in Section II to determine an optimal estimate of ρ that minimises ǫ. The expression for MSE shown in Eq. 7 was modified to incorporate unity trace constraints in a manner that makes no other assumption about the physicality of the estimated state. The convergence properties of this expression in this application have not been comprehensively investigated. Similarly, Eq. 11 makes no attempt at optimisation of the likelihood or entropy of the reconstructed state [31]. This is identified as an area for further investigation. Nevertheless, the results of reconstruction of the input state photostatistics with a)n i = 4.6 ± 0.8 and b) 6.5 ± 0.8 photons/pulse. are quite impressive, as shown in Fig. 12.
To permit validation of the reconstructed distributions, the mean photon numbers of the reconstructed density matrices were determined vian = K n=0 nρ nn . Poissonian distributions with the corresponding mean photon numbers are shown beside the reconstructed distributions of Fig. 12 for comparative purposes. The agreement between the coefficients is remarkably close, confirming the practical suitability of the proposed signature approach for the reconstruction of photostatistics.

VI. CONCLUSION
We have the described the use of photon count signatures from an N-port photon counting device for the purposes of reconstructing the photon number probability distribution of an input state, and demonstrated this technique experimentally on a 9-port system. We have shown that exploiting signature probabilities offers superior performance to approaches reconstructing photon number probabilities from the measured binomial statistics. Our method offers advantages in reduced statistical uncertainty, relative insensitivity to calibration errors and reduced reliance on a-priori assumptions about the input state.