Rigorous calibration method for photon-number statistics

Characterization of photon statistics of a light source is one of the most basic tools in quantum optics. Although the outcome from existing methods is believed to be a good approximation when the measured light is sufficiently weak, there is no rigorous quantitative bounds on the degree of the approximation. As a result, they fail to fulfill the demand arising from emerging applications of quantum information such as quantum cryptography. Here, we propose a calibration method to produce rigorous bounds for a photon-number probability distribution by using a conventional Hanbury-Brown-Twiss setup with threshold photon detectors. We present a general framework to treat any number of detectors and non-uniformity of their efficiencies. The bounds are conveniently given as closed-form expressions of the observed coincidence rates and the detector efficiencies. We also show optimality of the bounds for light with a small mean photon number. As an application, we show that our calibration method can be used for the light source in a decoy-state quantum key distribution protocol. It replaces the a priori assumption on the distribution that has been commonly used, and achieves almost the same secure key rate when four detectors are used for the calibration.


I. INTRODUCTION
Autocorrelation measurement of a light source has been known to be a convenient method for investigating the characteristics of the source. It dates back to Hanbury-Brown and Twiss (HBT) who measured [1,2] the correlation of photocurrents from two photodetectors shone by a common light source to determine its second-order intensity correlation function. The meaning of the correlation functions in quantum optics was clarified by Glauber [3,4], who showed that the second-order correlation function G (2) (r 1 , r 2 ; τ ) = I(r 1 , t)I(r 2 , t + τ ) can also be determined from a HBT setup with two photon detectors by measuring the rate of coincidence counts. It was successfully used in the direct observation of the non-classical property of light [5]. In the case of a pulsed light, a proper integration around τ = 0 and a normalization lead to a normalized factorial moment g (2) (0) = n(n − 1) / n 2 of the photon number n in the pulse [6]. Since g (2) (0) = 0 is achieved only by an ideal single photon source (SPS), the measurement of g (2) (0) by the HBT setup has been widely adopted for the characterization of experimentally developed SPSs [7,8]. The extension of the method to higher-order moments g (r) (0) is also straightforward by increasing the number of detectors to r [9,10].
Despite the apparent success for the non-classical light sources, the above characterization method is not particularly suited to the emergent applications in quantum information. The most severe drawback is the fact that the HBT setup with conventional threshold photon detectors provides the accurate value of the factorial moment only in the limit of low detection efficiencies. The value of g (r) (0) from an actual experiment is only approximate, and how it is deviated from the true value is unknown. This is problematic for the applications in-volving the security [11,12], for which rigorous security statements are mandatory. For the computational tasks using photons such as boson sampling [13], reliability of the outcome will eventually be ascribed to that of the constituent components including light sources. Another drawback is that we encounter the moments {g (r) (0)} r much less frequently in the quantum information theory than the photon-number probability distribution {p n } n . Photon-based protocols concern the presence of a photon in a pulse, and hence the relevant quantity is the corresponding probability p 1 . The requirement for the light source used in the decoy-state quantum key distribution (QKD) [14][15][16] is given by a set of inequalities in terms of the photon-number probability distribution {p n } n [17,18]. While the knowledge on {g (r) (0)} r may be converted to bounds on {p n } n , it is not straightforward due to the unboundness of the photon number n. Hence, as a calibration method in the era of quantum information, it is vital that it provides rigorous statements over {p n } n , and preferably it is tight.
In this paper, we propose a calibration method to obtain rigorous bounds on {p n } n≤D−1 using a HBT setup with D threshold photon detectors. It is flexible and versatile since the formula for an arbitrary number of detectors is given in a closed-form expression and the detection efficiencies do not need to be uniform. The explicit formula enables us to cope with various imperfections such as ambiguity in detection efficiencies and statistical uncertainty in the outcomes. Our method is optimal for coherent states and thermal states with a small mean photon numbers. As an application of our calibration methods to QKD, we modify the security analysis of the decoy-state BB84 protocol to accommodate the use of a calibrated source. It serves to bolster the security by lifting the a priori assumption on the source, and we further show that the performance drop is minimal when a arXiv:1710.00457v1 [quant-ph] 2 Oct 2017 1. An implementation of the D = 4 case of our calibration method. The overall detection efficiencies are given by η1 = R1T2η , where T k (k = 1, . . . , 3) and R k (k = 1, . . . , 3) are transmittance and reflectance, respectively. D = 4 setup is used for the calibration.

II. CALIBRATION METHOD FOR THE PHOTON NUMBER DISTRIBUTION OF A PHOTON SOURCE
We consider a generalized HBT setup with D photon detectors, where the light pulses to be calibrated are divided into D parts by beam splitters, and various rates of coincidence detections are recorded. We assume that each detector is a threshold detector with quantum efficiency η (det) i , which is modeled as reporting a presence of nonzero photons after a linear absorber with transmission η (det) i . Figure 1 shows an example for D = 4. The configuration of the beam splitters is not relevant. The only relevant system parameters are the overall efficiency for each detector, η i (i = 1, . . . , D), including the branching efficiencies. We denote their average as η := D −1 i η i . We define the averaged r-tuple coincidence probability c obs,r to represent the observed data. For the r-tuple coincidence, there are D r different combinations of r detectors, and the associated probabilities. The averaged quantity c obs,r is their arithmetic mean. Let us first consider the case where the measured pulse contains exactly n photons. The coincidence detection probability of the first to the r-th detectors is given by . . , r} and |W | is the cardinality of set W . The averaged probability c obs,r for the n-photon input is then given by where we defined ω r,j := D−j r−j / D r and I j := {W ⊂ Z D | |W | = j}. Note that c n,r = 0 if r > n. For the case of a general input pulse with distribution {p n } n , the averaged coincidences should satisfy c obs,r = ∞ n=0 p n c n,r (r = 1, . . . , D).
Our goal is to find rigorous bounds on {p n } n under the constraint of Eq. (2). To describe the formula for the bounds, it is convenient to introduce vectors of order D + 1 as follows. Let c obs := (1, c obs,1 , . . . , c obs,D ), c n := (1, c n,1 , . . . , c n,D ), and c ∞ := (1, 1, . . . , 1). Consider a basis {c j } j∈S specified by the index set S = {0, 1, . . . , D}, and let {d (S) j } j∈S be its reciprocal basis, namely, c j · d , respectively. We postulate that the variations among the efficiencies {η i } i are moderate, and they satisfy i∈W for any W, W ⊂ Z D . Now our main result is stated in the form of the following theorem.
Theorem 1. Under the condition (3), the constraint (2) implies for n = 0, . . . , D − 1 and We first prove the inequalities (4), (7), and (8) involving S. The crux of the proof is to show that c m · d with α W := −ln(1 − i∈W η i ) > 0. Assume that the constants {z j } j are given by z j = ω j · d (S) n , where ω j := (ω 0,j , ω 1,j , . . . , ω D,j ) with ω r,j := 0 for r < j. Let us show that the numbers of zeros of f (m) and its derivative f (m) are D and D − 1, respectively. The prerequisite (3) ensures that there exists {β j } D−1 j=1 such that max W ⊂Ij α W < β j < min W ⊂Ij+1 α W . Let us temporarily assume that all {z j } j are nonzero and have the same sign. Using the notation D γ := e −γm ( d dm )e γm , we have We see that the coefficients of e −α W m for j = 1 and 2 have the same sign. In a similar vein, all the coefficients in D β D−1 · · · D β1 D 0 f (m) have the same sign, implying that this function has no zeros. Since multiplication of e ±γm does not change zeros, Rolle's theorem assures that f (m) should have no more than D zeros. If some of {z j } j were zero or had the opposite sign, we could remove the zeros by a fewer number of operating D γ , contradicting the fact that f (m) has at least D zeros. Hence, we conclude that f (m) has exactly D zeros at m ∈ S \ {n}, and f (m) has exactly D − 1 zeros.
Since Rolle's theorem also implies that the D − 1 zeros of f (m) must be strictly between neighboring zeros of The remaining inequalities (5), (6) and (9) involving the index set S can be proved in a similar way. For Eqs. (5) and (6), the function f (m) has only D −1 known zeros at S \ {n, ∞}, but it is compensated from the fact that z 0 = 0. The detail is given in Appendix A, which concludes the proof of Theorem 1.
Calculation of the bounds p L n and p U n in Eqs. (4)-(9) is straightforward because C and C are triangular matrices. Closed-form expressions for D = 2, 3 and 4 are given in Appendix B in terms of normalized quantities c obs,r := c obs,r /c r,r for convenience, since they are O(1) in the limit of η i → 0. In Fig. 2, we show the bounds obtained from the D = 4 setup when it is applied to an ideal Poissonian source with p n = e −µ µ n /n! along with the true values. As shown in Fig. 2, the bounds are fairly good for µ ≤ 0.5, and gradually become worse as µ gets larger.
When we compare the explicit formulas for different values of D, we notice that sometimes the dominant O(1) terms in η do not change as an increment of D, like Eq. (B2) to Eq. (B8), and Eq. (B10) to Eq. (B18). In general, we can prove that c obs · d (S) n for D = D 0 and c obs · d (S ) n for D = D 0 + 1 coincide in the limit of η → 0 for a given source (see Appendix C). In other words, each bound p U,L n shows a major improvement for every other increment of the number D of detectors used in the set up. The origin of this unexpected behavior may be understood from the threshold or saturation behavior of the detectors, which an adversary may exploit to fool us to believe in a wrong distribution. Given a distribution {p n } n and the corresponding {c obs,r } r , one may change {p n } n only slightly by replacing a O(η D ) potion of the pulses with ones with an extremely large photon number to flood the D detectors. Through this small modification, the adversary can increase c obs,D to any value without changing the rest of {c obs,r } r . Existence of such an attack forces us to trust the observed value of c obs,D only in one direction, which results in improving only half of the bounds. The above argument tells us that, in the derivation of the rigorous bound, it is essential to take the saturation behavior of the detectors into consideration.
The present formula allows us to check the optimality of the obtained bounds. From the first row of the equality for 0 ≤ n ≤ D (or equivalently, if all the S-related bounds are nonnegative), we find that {p n } n defined as p n = c obs · d (S) n for 0 ≤ n ≤ D and p n = 0 for n ≥ D + 1 is a probability distribution, and it fulfills Eq. (2) as is seen from the remaining rows of CC −1 c T obs = c T obs . The inequalities (4), (7) and (8)  for n = n 0 , and p n = 0 otherwise. It saturates the inequalities (5), (6) and (9), while it fulfills Eq. (2) in the limit of n 0 → ∞, implying the optimality of Eqs. (5),(6) and (9).
Thanks to the closed-form expression of the bounds, we may discuss what types of light sources are optimally calibrated. We consider the limit of η → 0. For weak light sources where p n rapidly decreases with n, the most severe condition for the optimality of the S-related bounds is expected to be p L D−1 = c obs · d Equations (12) and (13) form a sufficient condition for the optimality of the S-related bounds. We can obtain the conditions for S -related bounds by replacing D with D−1. The proof of these statements is given in Appendix C. Examples for sources satisfying these conditions are a coherent light source with n < 1 and a thermal light source with n < 1/D. The closed-form expression also allows us to adapt our method easily to the cases where there are ambiguities in the values of η i and c obs,r , simply by calculating the worst-case values and by introducing confidence levels if necessary. It also helps us to estimate how the degrees of such ambiguities will affect the tightness of the bounds. As an example, we consider the D = 4 setup with a uniform efficiency η 1. We assume that the actual distribution {p n } n is similar to that of a weak coherent light source or a single photon source, namely, p n p n+1 for n ≥ 1. It implies p 0 ∼ 1 −c obs,1 and p n ∼c obs,n for n ≥ 1. Since c r,r = r!η r , we have ∆c obs,r /c obs,r = ∆c obs,r /c obs,r − r(∆η/η). By applying it to the dominant terms in Eqs. (B15)-(B22), we obtain ∆p U,L n p n |∆c obs,n | c obs,n + n |∆η| η (n ≥ 1).
These show that the relative errors in c obs,r and η affect the bounds only proportionally, unless p 0 = 0. While the above relations are derived from the explicit form for D = 4, it is not difficult to show that they hold for arbitrary values of D.

III. THE DECOY-STATE BB84 PROTOCOL WITH A CALIBRATED SOURCE
As an application of the proposed characterization method, we consider the calibration of the laser light source used in the decoy-state BB84 protocol, which is one of the most frequently demonstrated QKD protocols. Here we focus on the protocol using pulses with three different intensities, termed signal, decoy, and vacuum. Let {p n } n and {p n } n be the photon-number distributions of the signal and the decoy pulses, respectively. We assume that the vacuum pulses contain no photons. The crux of the decoy-state protocol is to estimate the amount of valid signals, namely, the conditional detection probability Y 1 given the sender emits exactly one photon (Y 1 is often called the one-photon yield). Comparison between the observed probability of the detected signal pulses, Q, and that of the detected decoy pulses, Q , establishes a lower bound Y L 1 on Y 1 . Combined with an upper bound e L 1 on the error probability in the single photon emission events, an asymptotic secure key rate is given by [11] where q is the probability of choosing the signal pulse and E is the observed error probability for the signal pulses, with q and E similarly defined for the decoy pulses. The security of the protocol has usually been proved under the a priori Poissonian assumption, namely, p n = exp(−µ) µ n n! and p n = exp(−µ ) µ n n! where µ and µ are the mean photon numbers of the signal and the decoy pulses, respectively. For 1 > µ > µ , the bounds [18] are then given by where the zero-photon yield Y 0 is determined from the observed probability of the detected vacuum pulses, and e 0 = 1/2. Since the Poissonian assumption involves infinite number of conditions, it is impossible to verify it experimentally. The proposals [17,19,20] for the use of light sources other than lasers also rely on an infinite set of conditions. It is thus important to replace such a priori assumptions with experimentally verifiable ones [21,22]. What we seek here is to use experimentally available bounds from Theorem 1. Wang et al. [18] have shown that Eqs. (17) and (18) are still valid by replacing each term by the worst-case values, provided that p L n /p U n ≥ p L 2 /p U 2 ≥ p L 1 /p U 1 holds for all n ≥ 2. We can easily extend it to a bound valid when p L n /p U Eq. (18) is straightforwardly modified to The detailed derivation is given in Appendix D. Figure 3 shows a comparison of the asymptotic secure key rate with the a priori assumption of ideal Poissonian and those with our calibration method. To calculate these curves, we used Eqs. (17) and (18) for curve (a) and Eqs. (19) and (20) for curves (b) and (c) in Fig. 3. For curve (d) with D = 2, Eq. (19) does not hold, and we used a trivial bound shown in Appendix D. We assume that the overall quantum efficiencies of the detectors are uniform and known with an accuracy of 1 percent. From  Fig. 3, we find that the secure key rate improves as D increases, and that for D = 4 it is comparable to that with the a priori assumption of ideal Poissonian.

IV. SUMMARY AND DISCUSSION
We have presented a calibration method for photonnumber distribution of a light source and shown the ex- plicit formula for rigorous bounds on probabilities for small photon numbers. We believe that our calibration method makes a significant contribution to the quantum optics toolbox, and is especially useful in quantum information, such as in applications involving security and in computational tasks whose outcome cannot be verified efficiently.
For the measurement of the photon-number distribution of a light source, other methods are known such as the one using homodyne tomography [23] and the one using a photon-number-resolving detector [24][25][26]. Compared to a threshold detector, these devices can singly respond to two or more photons. This leads to a more compact setup for the calibration, and will be useful for rough estimates. On the other hand, these devices tend to require many system parameters to model them, which may make it difficult to produce a rigorous bound. They cannot avoid saturation behavior either, namely, a photon-number-resolving detector can resolve photons only up to a certain photon number, and homodyne tomography usually introduces a threshold manually to avoid artifacts. On these grounds, we believe that the HBT setup has an advantage of being made up of components, each of which is represented by a simple model. Of course, this inevitably raises a question of how closely an actual detector behaves to the threshold detector model. An obvious deviation is the dark counting, which may be treated as an estimation problem of genuine coincidence rates from the actually observed coincidence rates including the contribution of the dark counts. It is also important to verify the validity of the threshold detector model in future experimental researches.

ACKNOWLEDGMENTS
This work was funded in part by ImPACT Program of Council for Science, Technology and Innovation (Cabinet Office, Government of Japan), Photon Frontier Network Program (Ministry of Education, Culture, Sports, Science and Technology), and CREST (Japan Science and Technology Agency).
The formula for the uniform case of η = η 1 = η 2 = η 3 = η 4 is simply given by setting ξ i,j = 0 for all i, j.