Absolute efficiency estimation of photon-number-resolving detectors using twin beams

A nonclassical light source is used to demonstrate experimentally the absolute efficiency calibration of a photon-number-resolving detector. The photon-pair detector calibration method developed by Klyshko for single-photon detectors is generalized to take advantage of the higher dynamic range and additional information provided by photon-number-resolving detectors. This enables the use of brighter twin-beam sources including amplified pulse pumped sources, which increases the relevant signal and provides measurement redundancy, making the calibration more robust.


Introduction
Quantum optics enables one to make measurements that are more precise than the fundamental limits of classical optics [1]. Central to this capability are quantum optical detectors, those that are sufficiently sensitive to discern the inherent discreteness of light. These detectors are key to emerging quantum technologies such as quantum imaging and lithography, in which the standard wavelength limit to resolution is surpassed by using quantum states of light and photonnumber sensitivity. However, the majority of quantum optical detectors have a response that sat-urates at only one photon, imposing a significant limitation on the brightness of the optical fields that can be used for such quantum technologies. The result of this binary detector response is a measurement that can discriminate only between zero photons and one or more photons arriving simultaneously; the detector produces an identical response for any number of photons greater than zero. To allow increased brightness of the light sources used in the technologies outlined above (and the concomitant improvement in accuracy that this brings), it is necessary to use detectors that can discern the number of photons incident simultaneously on the detector -a photon-number-resolving detector (PNRD). Indeed, the development of PNRDs is an active area of research including photomultiplier tubes, the extension of avalanche photodiodes (APDs) to higher photon number through time multiplexing [2,3] or spatially-multiplexed arrays [4], transition-edge sensors (TES) [5,6,7], charge-integration photon detectors (CIPS) [8], superconducting single-photon detectors [9,10], visible-light photon counters (VLPCs) [11], and quantum dot field effect transistors [12,13].
Calibration of optical detectors is a difficult problem. The standard approach uses a previously calibrated light source. The drawback of this approach is that errors in the source brightness translate directly into errors in the detector efficiency calibration. The converse is also true, leading to a detector and light-source calibration dilemma. To get around this, brightness calibration is typically based on a fundamental physical process, for example, the luminosity of blackbody radiation of gold at its melting point [14], or heating of a cryogenic bolometer [15]. Such methods are suitable for calibrating bright light sources. In contrast, using this method to calibrate detectors operating at the quantum level (i.e. fields containing only a few photons) requires sources with powers on the order of a femtowatt to avoid saturation. Such sources are impractical.
Nonclassical states of light allow us to circumvent the horns of the dilemma due to their behavior in the presence of loss. Realistic detectors can be modeled by optical loss (i.e. attenuation) followed by unit efficiency detectors. The transformation of classical states (e.g. a coherent state, thermal state, etc.) under loss is parameterized only by source brightness, which scales linearly with the loss. In contrast, nonclassical states change their statistical character upon experiencing loss: correlations in optical phase, intensity, photon number, and electric field transform in a non-trivial way under loss. Based on this fact, Klyshko proposed a way to calibrate detectors based on the statistical character of light rather than its brightness [16]. This approach relies on spontaneous parametric downconversion (SPDC) as a light source in which a photon is simultaneously created in each of two optical modes, usually denoted signal and idler. Since the photons are created in pairs the two output modes are perfectly correlated in photon number. Thus detection of the idler without the simultaneous detection of the signal photon can only be caused by loss in the signal arm. Measuring the detected rate of idler photons R i and photon pairs R c then allows a calibration of the detector efficiency η s and vice versa with i ↔ s. This efficiency estimation was shown experimentally in [17,18,19,20]. The Klyshko scheme is limited by three factors. First, it relies on the implicit assumption that at most one photon pair is emitted at a time, a feature which can be violated if the SPDC is pumped strongly. Indeed more than one pair is often desirable, as in continuous-variable experiments. In this case, the detector efficiency can be obtained by first lowering the pump power and using the Klyshko method [21]. However, the detector is now calibrated outside the regime of its intended use, which could necessitate subsequent assumptions such as the independence of the estimated efficiency from the pump power. Thus, the direct in situ calibration of detectors would be desirable. Second, the Klyshko scheme is primarily designed for single-photon detectors and is not directly translated to PNRDs. And third, this calibration method is highly sensitive to the input state quality and measurement uncertainties. Because the number of measurements taken is exactly the number needed to determine the efficiency, any errors in the measurements will propagate directly into the efficiency estimation.
Here we present a technique for measuring the absolute quantum efficiency of a detector based on the Klyshko method, but explicitly taking into account multiple photon events. This improves the calibration accuracy and allows for calibration of a PNRD. This approach utilizes the PNRD capability to measure the photon-number distribution of an optical mode. Using two PNRDs the joint photon-number statistics between the two electromagnetic field modes, including photon-number correlations and individual photon-number distributions, of the SPDC source can be determined. For each element of the resulting joint photon statistics, one can find a formula giving the detector efficiencies of the two PNRDs. For the zero-and one-click elements, this reproduces the Klyshko result. However, the increased dynamic range of PNRDs allows the use of brighter sources that produce measurable rates in the higher photon-number elements of the joint statistics. We then use optimization techniques to estimate the detector efficiencies from the increased number of measurements. This added redundancy improves the tolerance of our scheme to background light and statistical noise. We experimentally demonstrate this efficiency estimation method with two time-multiplexed PNRDs [3].
We begin by introducing a general treatment of PNRDs and then use this as a basis for describing our generalized Klyshko method. This is followed by the description of an experiment to test the efficiency estimation with PNRDs.

Photon-number-resolving detectors
Photon-number-resolving detectors (PNRDs) are a class of photodetectors that have a unique response for every input photon-number state within their range. Ideally these responses can be perfectly discriminated. However, the less than perfect efficiency of realistic detectors causes these responses to overlap, and thus does not allow for direct photon-number discrimination. Overlap of detector responses can also arise from the detector electronics (e.g. amplification) or the underlying detector design. Despite this overlap, the linear relationship between the detector response and the input state allows for the reconstruction of the input photon statistics from the measured outcome statistics. This linear relationship is encapsulated by whereρ is the input-state density matrix, P n is the probability for the nth measurement outcome andΠ n is the associated positive operator-value measurement (POVM) operator. Consider the matrices expressed in the photon-number basis. Since PNRDs do not contain an optical phase reference, the off-diagonal elements ofΠ n are zero, meaning the photon-number-resolving detection is insensitive to off-diagonal elements inρ. It is thus useful to write the diagonal elements ofρ, the photon-number statistics, as a vector σ . Similarly, we write the outcome probabilities {P i } as a vector p. In the following, we truncate σ at photon number N − 1, where N is the number of detector outcomes, although this is not strictly necessary. The POVM operators of a general PNRD can be modeled by dividing the detector performance into two components: efficiency and detector design, described by matrices F and L respectively. In the photon number basis, a POVM element can be written as [24] where [F · L (η)] n,m corresponds to the probabilty of detecting n out of m incidence photons. As pointed out in the introduction, detector efficiency η can be modeled by a preceding optical loss l of (1 − η). In the context of a PNRD, loss causes the photon-number statistics to transform according to σ → L (η) σ , where These matrix elements transform the state by lowering photon numbers from j to i, representing the loss of photons through a binomial process with probability l = (1 − η). Although the detector-design component of the model depends on the detailed functioning of the device, a large class of PNRDs -mode-multiplexers -can be treated in the same way. These detectors divide an input optical field mode into many spatial and/or temporal modes and then use single-photon detection on each mode to achieve number resolution. Examples include nanowire superconducting detectors [9], VLPCs [11], intensified charge-coupled devices (CCDs) [22], integrated APD arrays [4,23], and time-multiplexed detectors (TMD). All of these detectors suffer from detector saturation; the one-photon detector response occurs if two or more photons occupy the same mode. This saturation effect is modelled by a detector design matrix F = C, the "convolution" matrix, which has a general but complicated analytic form given in [24] (in the context of the TMD). The form of C depends on relatively few parameters comprising the splitting ratios of the input mode into each of the multiplexed modes and the total number of these modes.
As an example, in a CCD array detector the pixel shape and size defines the detected optical mode. The spatial overlap of the detector mode of each pixel with the incoming optical mode then gives the corresponding splitting ratio. The total number of pixels is the total number of multiplexed modes. As the number of multiplexed modes goes to infinity the C matrix goes to the identity. All PNR detectors can be described by a POVM, which can be reconstructed by detector tomography [24,25]. For PNR detectors that do not rely on mode multiplexing, such as transition-edge superconducting detectors, and single APD detectors, one must factor the POVM elements into an F matrix and an L matrix to apply the calibration procedure described in this paper.
Focusing on the specific case of mode-multiplexed detectors, one can rewrite the Eq. (3) as, where the elements of the ith row of C · L (η) are the diagonals of the ith operator in the POVM set Π i . With PNRDs in two beams, denoted 1 and 2, one can measure not only the individual photonnumber statistics σ 1 and σ 2 , but also the joint photon-number distribution of these two beams. This distribution is written as the joint photon statistics matrix σ , where σ m,n is the probability of simultaneously having m photons in mode 1 and n photons in mode 2. We extend Eq. (5) to relate the probability P m,n , of getting outcome m at detector 1 and outcome n at detector 2, to the joint photon statistics σ , where subscripts indicate the relevant detector and A T is the transpose of A. Joint photon statistics are a measure of photon-number correlations in beams 1 and 2, and are thus sensitive to loss as discussed in the introduction. We use this description of PNRDs to generalize the Klyshko calibration.

Generalizing the Klyshko method
The assumption that only a single photon pair is generated at a time in the Klyshko efficiencyestimation scheme is only valid for very low SPDC pump powers. We can determine in what manner the Klyshko scheme breaks down when more than one pair is produced. Simply assuming an additional contribution of pairs σ 2,2 to the pair rate σ 1,1 changes the efficiency estimate of Eq. (1) toη This overestimates the efficiency for σ 2,2 = 0. Similarly, sensitivity to the detailed photon statistics of the input states acts to degrade the accuracy of the efficiency estimate. Let us consider the effect of background light (e.g. fluorescence from the optical elements and detector dark counts) on the efficiency estimation of Eq. (1). Background photons are uncorrelated between the signal and idler beams, increasing the singles rate, and thus increasing R i . This will make the estimated efficiency of the detectorη s lower than the actual efficiency η s . To generalize the Klyshko scheme to PNRDs and avoid the above limitations we need to determine the true state generated by SPDC photon-pair sources. Ideally, these sources produce a 'two-mode vacuum squeezed state' where λ is proportional to the pump beam energy, and |n s(i) is an n-photon state of the signal (idler) mode. Having only one free parameter, it is tempting to use this state in the generalized efficiency estimation. However, this state can only be generated with careful source design [26]. Instead, photons are typically generated in many spectral and spatial modes in the signal and idler beams [27]. Depending on the number of modes in the beams, the thermal photon-number distribution in Eq. (8) changes continuously to a Poisson distribution [28]. Consequently, it would be incorrect to assume the source produces an ideal two-mode vacuum squeezed state. Still, the number of photons remains perfectly correlated between the two beams. Without access to the number of generated modes we can make only the following assumption about the joint photon statistics of the source where {c i } are arbitrary up to a normalization constant and δ m,n is the Kronecker delta. Inserting the joint photon statistics defined by Eq. (9) into Eq. (6) we arrive at the basis for our generalized Klyshko method. Since C 1 and C 2 of the detectors are known the predicted joint outcome probabilities P are highly constrained having N 2 elements uniquely defined by the N parameters in {c i } and the two efficiencies η 1 and η 2 . Consequently, a measurement of the joint outcome statistics specifies η 1 and η 2 with a large amount of redundancy; the efficiencies are overdetermined. In order to correctly incorporate all measured outcome statistics into the efficiency estimates a numerical optimization approach is used. We minimize the difference G between the measured outcome statistics R and the predicted outcome statistics P, which are determined by {c i } , η 1 , and η 2 This is done by minimizing the Frobenius norm F = Tr (G) 2 1/2 to find the optimal η 1 and η 2 -our estimates of the PNRD efficiencies. Using the Frobenius norm makes this a leastsquares optimization problem over {c i } , η 1 and η 2 , where 0 ≤ η i ≤ 1. However, this method of efficiency estimation is amenable to other optimization techniques such as maximum-entropy or maximum-likelihood estimation. We use the Matlab® function 'lsqnonneg' (with the constraints σ m,m ≥ 0). This efficiency estimation is similar to the optimizations performed in state or process tomography, which are known to be convex problems (i.e. there are no local minima) [29]. However, we are now estimating parameters in both our stateρ and the POVM set Π i . Because this is not guaranteed to be convex [29], we simulated a variety of measured statistics to test for a single minimum. Using C 1 and C 2 for the PNRDs in our experiment (TMDs), with several sets of photon statistics {c i } , and a range of efficiencies η 1 and η 2 , we simulated various outcome statistics. In all cases, the optimization reproduced the correct efficiencies and only a single minimum was observed. Figure 1 shows the result of a typical simulation displaying a global minimum in the optimization residual (the minimum value of F for a given pair of efficiencies η 1 and η 2 ). This is a good indication that efficiency estimation is a convex problem.

Experimental setup
To experimentally demonstrate and test our efficiency estimation method, time-multiplexed detectors (one possible realization of a PNRD) were employed. In a TMD, the input optical state is contained in a pulsed wavepacket mode. The pulse is split into two spatial and several temporal modes by a network of fiber beam splitters and then registered using two APDs [3]. APDs produce largely the same response for one or more incident photons. The TMD overcomes this binary response by making it likely that photons in the input pulse separate into distinct modes and are thus individually registered by the APDs [3]. The TMD is a well-developed technology, which makes this an ideal detector to test our approach to detector efficiency estimation. The convolution matrix C for this detection scheme is calculated from a classical model of the detector using the fiber splitting ratios [3], and is also reconstructed using detector tomography [ 24,25]. Loss effects in TMDs have also been thoroughly investigated [30]. In our experiments, the TMDs have four time bins in each of two spatial modes, giving resolution of up to eight photons, with a possible input pulse repetition rate of up to 1 MHz. Field-Programmable-Gate-Array (FPGA) electronics are used to time gate the APD signals with a window of 4 ns, which significantly cuts background rates. The joint count statistics R are accumulated by the electronics and transferred to a computer for data analysis. The experimental setup is centered on a nearly-two-mode SPDC source [26] as depicted in Fig. 2. The twin beam state produced by SPDC in a potassium dihydrogen phosphate (KDP) nonlinear crystal consists of two collinear beams with orthogonal polarizations. The pulsed pump (415 nm central wavelength) driving the SPDC is a frequency-doubled amplified Ti:Sapphire laser operating with a 250 kHz repetition rate. A pick-off beam is sent to a fast photodiode (PD) that is used to trigger the detection electronics. Dichroic mirrors (DM) and a red-pass color glass filter (RF) are used to separate the blue pump from the near-infrared (830 nm central wavelength) SPDC light. A polarizing beam splitter (PBS) is used to separate and direct the two co-propagating downconversion beams into separate single-mode fibers (SMFs) connected to the time-multiplexed detectors (TMD1 and TMD2). The joint statistics R of the two TMDs are recorded for a range of pump powers between 1 and 55 mW in order to estimate the two TMD efficiencies at each power. To examine the spectral response of the detector efficiency, one could tune the wavelength of the SPDC source by adjusting the pump wavelength and the crystal orientation.

Estimated efficiencies
For each pump power we determine the optimum detector efficiencies that are consistent with the measured statistics at both PNRDs, shown in Fig. 3. Three different regimes are observed. At low powers (up to 6 mW) the estimated efficiency increases with power. Between powers of 6 and 40 mW the estimates appear constant. At 40 mW there is a sudden jump in the estimated efficiency (approximately twice the previous value); above this power the estimates remain constant. Continued investigation revealed that the second-harmonic generation (SHG) process qualitatively changed its behavior at 40 mW: the increased pump power induced unwanted higher-order nonlinear effects, resulting in the generation of additional frequency components other than the second harmonic and a change in the spatial mode structure. This changed both the transmission of the short wave pass filter (DM and BF) and the efficiency of the fiber coupling into our detectors. Since this behavior is outside the scope of our investigation we omit this data from further discussion.
The TMD efficiency should be independent of the average photon number of the state and thus the pump power. This is true in the second region of Fig. 3 but not the first. By reconstructing the joint photon-number distribution of the input state (σ m,n ) using the estimated efficiencies, one can gain insight into the estimation accuracy of the detector efficiency. This serves as a partial check for our assumption that the number of photons in the two beams is equal. In Fig. 4, we show the reconstructed joint photon statistics for two pump powers in regimes 1 and 2. In the second regime, the photon-number distribution is largely diagonal: only 10% of the incident photons arrive without a partner in the other beam. In contrast, the state in the low power regime has significant off-diagonal components, with 43% of the photons arriving alone. This suggests that at low powers the reference state is corrupted by background photons, possi- bly fluorescence from optics in the pump beam path, pump photons leaking through our filters, or scattered pump photons penetrating our fiber coatings. Contributions from dark counts are expected to be negligible, since the specified dark count rates of our detectors are significantly lower than these other effects. In the next section we attempt to remove this background in order to better estimate the efficiency.

Compensation for background light
We investigate two methods of dealing with background light. In the first, a parameterized background contribution is incorporated into the efficiency estimation procedure. The second attempts to measure the outcome distributions due to the background alone and then subtract these from the efficiency data. As pointed out in Section 3, the efficiency estimation is greatly overdetermined. One expects that a modeled background could be added to our input state model, Eq. (9), without jeopardizing the convergence of the optimization. This background would be entirely uncorrelated, possibly making its contribution to the joint outcome statistics easily distinguishable from the SPDC contribution. We model the photon-number statistics of the background in each beam by a Poisson photon-number distribution d(n) = α n exp(−α)/n! [31], which is fixed by a single parameter -the average background photon number, α = n B . Consequently, only two additional parameters (one for each beam) enter into the efficiency estimation, keeping it overdetermined. Unfortunately, we theoretically found that the loss in one beam transforms the outcome statistics in a manner similar to background in the other beam. Thus, the problem is no longer convex; there is a set of equally optimal points {( n B , η)}. To show this we compare two different two-mode number-correlated states σ A(B) , similar to the state in Eq. (9), that undergo the addition of uncorrelated background light and loss respectively. The addition of background to beam 1 of state σ A is given by the convolution of the Poisson distribution with input state, σ A (defined by arbitrary c A ). The elements of the background-added joint photon-number statisticsσ A are given by the sum of all possible ways to add photons from the Poisson background, d, to the first beam of the initial state, σ A , that add up to a particular number of photons in the final state,σ Ã This can be written in terms of a matrix product where the elements of D ( n B ), D m,n = d(n − m) are the probabilities of having an additional n − m photons from the Poissonian background. To compare this background-added case with the situation in which there is no background, but there is loss, we assume a loss l = (1 − η) in beam 2 of the second state σ B (defined by arbitrary c B ). We then attempt to show that for some σ A and σ B (each with perfectly correlated photon statistics) the two resulting states are equal, that is Elimination of c B and c A from the resulting equations is facilitated using a computer algebra program, and we find that there is a range of l and n B ≤ 1 that solve Eq. (13), indicating that it is not possible to simultaneously fit for background and efficiency. This emphasizes the fact that one cannot distinguish between loss and Poisson background for the number-correlated states. Figure  the loss curve converges to line of slope 0.068 suggesting that even PNRDs with an infinite photon number range would not allow one to distinguish loss and background. This suggests that another approach to background light should be used. Another approach to address the background contribution to the efficiency estimation attempts to subtract an independently measured background contribution from the click statistics. For each pump power the pump polarization is rotated by 90 degrees, extinguishing the SPDC and allowing the measurement of the joint outcome statistics due to background alone. Generally, the statistics of two concurrent but independent processes is the convolution of the statistics of the two processes. However, the measured outcome probabilities P M of both light sources combined is not a simple convolution of the background outcome probabilities P B and the twin-beam outcome probabilities P S This is due to the fact that there is an effective interaction between the background and twin-beam signal detection events due to the strong detection nonlinearity of the APDs (i.e. the detectors saturate at one photon). However this saturation effect can be eliminated by first applying the inverse of the C matrices to the measured statistics which gives the estimated photon-number statistics prior to the mode multiplexing. These nonmode-multiplexed statistics can be assumed independent so that, We then use the convolution theorem to find where F indicates the Fourier transform, and the matrix division is element by element. Using P S we estimate the efficiency as previously by the method in Section 3.
To test the accuracy of this background subtraction method we reconstruct the joint statistics at the same powers as in Fig. 3. We find that in both cases the off-diagonal components are significantly reduced. In the low power regime, now only 16% of incident photons are not part of a pair, and at higher power this becomes just 4%. With the background subtracted, the estimated efficiencies are plotted in Fig. 7 and are now in better agreement with the expected constant detector efficiency through the first two regions. However, the estimates still drop off as the power goes very low. Although we measure the background, we do not do so in situ; we need to change the apparatus by rotating a half wave plate (HWP). Thus, we are not guaranteed that this is equal to the background present during the efficiency estimation. Furthermore, errors in background measurements are more significant at low powers as background then forms a larger component of the outcome statistics.
Also plotted in Fig. 7 is the Klyshko efficiency. In contrast to the increased dynamic range of our method the standard Klyshko efficiency increases with pump power, evidence that higher photon numbers in the input beams distort the estimated efficiency.
The average efficiency across the second region was found to be 9.4% ± 0.4% for detector 1 and 8.0% ± 0.4% for detector 2, where the errors are the standard deviations. These relatively  low efficiencies are only partly due to the quantum efficiency of the avalanche photodiodes themselves, which is specified to be 60% ± 5% at our wavelength. Bulk crystal SPDC sources ordinarily emit into many spatial modes, which makes coupling into a single-mode fiber difficult and inefficient. Typical coupling efficiencies are less than 30% [32]. As a final check of the reconstructed photon statistics, we plotted the average reconstructed photon number as a function of pump power in Fig. 8. The relationship is linear, as expected when taking into account the higher dynamic range of a TMD in comparison to standard APDs [33].

Conclusion
Despite being one of the seven base SI physical quantities, the working standards for luminous intensity have a relative accuracy of only 0.5% [34], compared to 10 −12 and better for the second [35]. Relying on a light beam of a known intensity, the efficiency calibration of detectors is similarly limited. Quantum states of light give us the opportunity to bypass the working standards and calibrate detectors directly. Indeed, the photon has been suggested as an alternative to the candela as the definition of luminous intensity [36]. We use twin-beam states and their perfect photon-number correlations to measure the efficiency of a photon-number-resolving detector. This type of state can be produced at a wide range of wavelengths from many different types of source including nonlinear crystals, optical fibers, periodically poled waveguides, and atomic gases. This method has the advantage that it only assumes perfect photon-number correlation and does not assume the state has a specific photon-number distribution, nor even that it is pure. Despite these seemingly detrimental assumptions, the efficiency estimation presented has a large amount of redundancy leading to a relatively small absolute error of 0.4%. We show that this measurement is independent of the average photon number of the state, unlike the Klyshko method, making it more widely applicable. In particular, it is ideal for characterizing photon-number resolving detectors, and for use with bright reference states. PNR detectors are undergoing rapid development via a number of competing technologies. These detectors will play an important role in precision optical measurements and optical quantum information protocols where photon-number resolution is necessary for large algorithms. Future development of direct efficiency calibration should focus on the issue of background, which can corrupt the state and thus the calibration.