General theory for decoy-state quantum key distribution with arbitrary number of intensities

We develop a general theory for quantum key distribution (QKD) in both the forward error correction and the reverse error correction cases when the QKD system is equipped with phase-randomized coherent light with arbitrary number of decoy intensities. For this purpose, generalizing Wang's expansion, we derive a convex expansion of the phase-randomized coherent state. We also numerically check that the asymptotic key generation rates are almost saturated when the number of decoy intensities is three.


Introduction
The BB84 protocol proposed by Bennett and Brassard [1] has been known as a famous protocol guaranteeing information theoretical security. Its security has been proved by Mayers [2] in the single-photon case. However, implemented channels usually have imperfect quantum state generators, which generate not a perfect single-photon but a mixture of the vacuum state, the single-photon state, and the multi-photon state. In fact, many implemented quantum communication systems are equipped with phaserandomized weak coherent light, which can be regarded as the above type of mixed state [3,4,5,6,7]. In order to guarantee the security in the above imperfect quantum communication systems, Gottesman et al. (GLLP) [8], and Inamori et al. (ILM) [9] have obtained a useful formula for the asymptotic key generation (AKG) rate, which needs the counting rate and the phase error rate of single photon state. If we treat the worst case of these parameters, the obtained AKG rate become worse. In order to estimate these parameters properly, Hwang [10] proposed the decoy method, in which we randomly choose different intensities of the coherent light. In the decoy method, a pulse for final key generation is called the signal pulse while a pulse only for estimation of these parameters is called a decoy pulse. When k intensities are choosed for decoy pulses, we need to prepare k + 1 intensities. Therefore, as was pointed out by Lo et al. [11], for the estimation of the counting rate of single photon state, we have to treat an infinite number of the unknown parameters only with k + 1 equations even if the counting rates of individual number states do not depend on the basis. (Note that the counting rates for multi-photon states depend on the basis as well as the number of photon. However, we consider this special case for simplicity in introduction.) In this case, its estimate can be derived from the restriction for all unknown parameters. Hence, lager number of intensities yields a more precise estimate of the counting rate of single photon state. The phase error rate of single photon state has the same characterization.
In order to reduce the cost of its realization, Wang [12] proposed to choose the minimum decoy intensity to be the vacuum. By using this idea, QKD system can be implemented with k kinds of intensities essentially. In the following, we adopt this assumption. In the case of k = 2, Wang [12] also proposed an expansion of the phaserandomized coherent states ∞ n=0 e −µ µ n n! |n n|, which describe the signal pulse and the decoy pulses. His expansion is different from the conventional expansion by number states in that the signal and decoy states can be described by stochastic mixtures of at most four states. Under his expansion, we obtain three constraint equations with four unknown parameters. He derived an estimate of the counting rate of the single photon state. Ma et al. [13] also derived the same estimate by a different method. Ma et al. [13] and Wang [14] obtained an estimate of the phase error rate of the single photon state, and completed the estimation of AKG rate in the case of k = 2, independently. Ma et al. [13] pointed out that if the second smallest intensity goes to 0, the required parameters can be solved perfectly. So, we obtain the best AKG rate, which is called the universal upper bound ‡. However, as was mentioned by them, it is impossible to control the infinitesimal intensity. Further, even if it is possible, if we choose small intensity, we cannot estimate the required parameters properly because the infinitesimal intensity makes our estimate fragile for statistical fluctuation. Thus, for a precise estimate, we need to fix the minimum difference between two intensities. That is, it is suitable to compare the qualities of our estimates by changing the number of decoy intensities under this constraint. For a good survey, please see Wang [15]. Now, we trun to the improvement of AKG rate. After GLLP-ILM's result, Lo [16] conjectured an improvement of their AKG rate with a rough idea of its proof. Based on his conjecture, Boileau-Batuwantudawe-Laflamme (BBL) [17] discussed a further improvement of the rate by taking into account the effects of dark counts in the detector. They pointed out that the AKG rate in the forward error correction case is different with respect to that in the reverse case. In the following, these rates will be called BBL formulas. Lo's conjecture has been proved by Hayashi [22] and Koashi [18], independently. Hence, it is suitable to apply these formulas to the decoy method.
In this paper, we treat the case of arbitrary number k of decoy intensities with BBL formulas when the intensity can be controlled §. For this purpose, we generalize Wang's expansion of the phase-randomized coherent states, in which k + 1 phase-randomized coherent states are given by convex combinations of k + 2 states, which are called basis states. The idea of this generalization is trivial, however, it is not trivial to derive the concrete forms of k + 2 basis states and to check their positivity. It is also needed to give the AKG rate formula by using the counting rates and the phase error rates of each intensities based on this expansion. For this purpose, we generalize mean value theorem and the concept of "difference". In the conventional mean value theorem, we treat the derivative of a given function f and the difference between outcomes of the function f in two points. In this paper, we extend the concept of "difference" to the case of n + 1 points, and derive its formula based on the n-th derivative. We also derive a formula for the generalized "difference". Using these formulas, we give concrete forms of k + 2 basis states and showed their positivity. We also derive a formula for estimating the counting rate and the phase error rate of the single photon state from the counting rates and the phase error rates of signal and decoy state. Then, the AKG rate formulas with BBL formulas are obtained. Moreover, we numerically compare the AKG rates of the cases of k = 2, 3, 4 and the universal upper bound in a proper model. As our result, the AKG rates of the case of k = 3 attain the universal upper bound within 1% error in casein the forward and reverse error case. That is, further great improvement cannot be expected even if the number k increases more than 3.
The organization of this paper is given as follows. Section 2 gives the AKG rates in both the forward and reverse cases as conjectured by BBL [17], taking into account ‡ Indeed, if we use a good code instead of random coding, there is a possibility to improve this bound. However, if we use random coding it is impossible to improve the universal upper bound. § For the case when the intensity fluctuates, see the papers [19,20] and the back note [26] of the paper [21]. the effect of dark counts. In Section 3, we treat the universal upper bound in a proper channel model. In Section 4, we review the results by Wang [12,14], Ma et al. [13]. In Section 5, we generalize Wang's expansion of the phase-randomized coherent state ∞ n=0 e −µ µ n n! |n n|. In Section 6, we give a lower bound of AKG rate based on the counting rates and phase error rates of respective intensities. In Section 7, we numerically compare the AKG rates in the case of k = 2, 3, 4 and its universal upper bound. In Section 8, the generalized "difference" is introduced, and all theorems are proved.

AKG rates with dark counts effect
When the transmitted pulse is given as a mixture of the vacuum, single-photon, and multi-photon, taking into account the effect of dark counts in the detector, BBL [17] conjectured the following AKG rate by use of Lo [16]'s conjecture, which is proved by Koashi [18] and Hayashi [22] with the threshold detector independently. Further, Hayashi [22] allows loss of detector if it is independent of the measurement basis. That is, we assume the same receiver model as Hayashi [22]. Let q 0 and q 1 be the probabilities detecting the the vacuum state and the single-photon state except for dark counts in detector. Let µ be the intensity of the signal pulse, which produces the raw keys. Then, We denote the dark count rate in the detector, the phase error rate concerning the single-photon state counted except for dark counts, the bit error probability of signal pulse, and the counting rate of signal pulse by p D , r 1 , s µ , and p µ , respectively. When we use the phase-randomized coherent light with the intensity µ, e −µ and µe −µ are the probabilities generating the vacuum state and the single-photon state. In this case, Eve has no information concerning Bob's bit of the pulses detected by dark counts and Alice's bit of the pulses whose state is the vacuum state. Hence, the AKG rate with the forward error correction is different from that with the reverse error correction, and they are given as follows. (See BBL [17] or Section IV A of Hayashi [22].) where η(s) is chosen so that 1 − η(s) is the coding rate of the classical error correction code and h(x) is given by In the above formula, we regard the simultaneous event of a dark count and a normal count as a dark count. This is because the collision of both photons disturbs the information of the normal count.
These are the rates per one pulse sent by Alice. Thus, the coefficient 1 2 corresponds to the probability that the basis of Alice coincides with that of Bob. For protocols attaining these rates, see Hayashi [22].
The dark count rate in the detector should be measured before the sending the quantum communication. In the above formulas, the parameters q 0 , q 1 and r 1 are not known a priori. These cannot be estimated if only the pulse with the intensity µ is transmitted. In the decoy method [10], we randomly change the intensity µ and estimate these parameters from the counting rates and error rates of individual intensities. In particular, the parameter q 0 can be estimated by the counting rate p 0 of the vacuum pulses so that we obtain q 0 = p 0 − p D . Once the parameter q 0 is fixed, both the AKG rates I → and I ← are given as monotone increasing functions concerning q 1 (1 − h(r 1 )). Hence, the remaining problem is the estimation of the parameters q 1 and r 1 , which is a more complicated problem.
In order to calculate these parameters, we need simultaneous equations concerning the individual intensities, in which the counting rates except for dark counts in detector and error rates of individual number states are treated as the unknown parameters. For simplicity, we assume that the counting rate except for dark counts in detector does not depends on the basis. Then, we letq n be the counting rate of the state |n n|, andr n be the error rate of the state |n n| with the × basis. Note thatq 1 = q 1 and r 1 = r 1 . Letting p i be the counting rate of the intensity µ i , and s i be the error rate of the intensity µ i with the × basis, we obtain [11] for i = 1, . . . , k and 0 ≤q n ≤ 1, 0 ≤r n ≤ 1.

Universal upper bounds on concrete noise model
In a real optical communication system, when the intensity is µ i the dark counts rate in the optical fiber or the detector is p 0 , which is equal to the counting rate of the vacuum pulse, the counting rate p i and the error rate s i are given by where the parameter α is the counting rate of the single photon except for dark counts in the the detector or interfusion on optical fiber, and the parameter s is the error rate of single photon for normal detection. In an implemented channel, the parameter is given by the function of the length of transmittance L as where θ is the efficiency of the detector, a 0 is loss coefficient in the detector, and a 1 is loss coefficient in the communication channel. When the detection probability of the n-photon state is 1 − (1 − α) n + p 0 , the relation (6) holds. When the error probability of the n-photon state is (7) holds. Hence, the quantities q 1 = α + p 0 − p D and r 1 = (3) and (4). Thus, we obtain where p( µ) = (p(µ 1 ), . . . , p(µ k )), s( µ) = (s(µ 1 ), . . . , s(µ k )). Therefore, applying (1) and (2), we obtain the following universal upper bounds of the AKG rates in the forward and reverse cases: respectively. When p 0 = p D = 0, both rates are equal to which can be approximated by when α is sufficiently small. As is shown in Appendix C, the optimum intensity is characterized as If we choose intensities of decoy pulses suitably, the estimated parameters q 1 and r 1 are close to α and s. Hence, it is suitable to choose the intensity of the phase-randomized coherent light producing the raw keys among the interval [0, 1].

Wang's expansion
We review the previous results concerning the estimation of the parameters q 1 and r 1 .
In this section, we consider only the case of p D = 0 because these results treated only this case. In order to avoid infinite number of unknown parameters in (3), Wang [12] proposed the following expansion of the states ∞ n=0 e −µ i µ n n! |n n| in the case of k = 2: where Based on this expansion, we define the parameters q 2 , q 3 , r 2 , and r 3 as for j = 2, 3. Then, we have the following equations instead of (3) and (5): Wang [12] gave a lower bound of q 2 in the following way. First, he calculated the maximum value q 2 U of q 2 with the conditions (22), (23), and q 1 , q 3 ≥ 0: Using (22), he essentially obtained a lower bound q 2,min of q 1 : Wang [14] also obtained an upper bound b 1,max of On the other hand, Ma et al. [13] treated the case where only two kinds of intensities are used for estimation of q 1 and r 1 . Their results can be translated into our case of k = 3 by putting the counting rate p 0 of the vacuum pulse into the parameter q 0 . Assume that µ 1 + µ 2 < µ 3 , µ 1 + µ 2 < 1, and µ 1 < µ 2 < µ 3 . Then, they obtained an lower bound of q 1 : and an upper bound of q 1 r 1 : Taking the limit µ 1 → 0, they obtained the the lower bound q 2,min of q 1 in the case of k = 2. Also, similarly, they independently obtain the upper bound b 1,max of q 1 r 1 . Consider how tight the bounds q 2,min and b 1,max are. As will be mentioned in Theorem 2, the minimum value of q 1 under the conditions (22), (23) and 1 ≥ q 2 , q 3 ≥ 0 is calculated as Similar to (22) and (23), instead of (4), the relations hold. Then, the maximum value b 1,2 max of q 1 r 1 under the conditions (27), (28), and 1 ≥ q 2 r 2 , q 3 r 3 ≥ 0 is calculated as From (27), (28), max , i.e., the bound b 1,2 max is better than b 1,3 U . Wang [14] proposed that the vacuum pulse and two kinds of intensities are used for estimating the parameters, and another intensity is used for signal pulse. In this method, it is possible to use the counting rate and the error rate of the signal pulse. Hence, it can be expected to improve the AKG rate by taking into account the the counting rate and the error rate of the signal pulse. That is, it is needed to discuss the case of k = 3. In this case, Ma et al.'s bound q 1,3 L can be applied. There is a possibility to improve existing bounds by extending Wang's expansion (14) and (15) to the case of k = 3. Further, as was pointed out by Ma et al. [13], we can expect that the AKG rates are close to the upper bounds I → (µ; p 0 , p D , α, s) and I ← (µ; p 0 , p D , α, s) if the number k of intensities is sufficiently large. In the following section, we concentrate to derive better estimates of q 1 and r 1 q 1 .  (14) and (15): Theorem 1 Assume that µ 1 < . . . < µ k . Define the basis state ρ i (i = 2, . . . , k + 1) as

Convex expansion of mixed state ∞
Then, ρ i is positive semi-definite, and Note that the coefficient µ 2 i n−2 t=1 (µ i − µ t )Ω n is positive. The quantities Ω 2 and Ω 3 coincide with those defined in (18) and (19). Hence, we can check that the above expansion of the case of k = 2 reproduces Wang's expansion (14) and (15). Further, Ω 4 and Ω 5 are calculated as and ).

General case
In the previous subsection, we treat the case where p D = 0 and the counting rates and the error rates of the × basis are equal to those of the + basis. Since the average state of the state ρ j with the × basis is different from that the + basis, we have to treat the parameter q j+k of the + basis as a parameter different from the parameter q j of the × basis in the general setting. In this subsection, we extend the discussion of the above subsection to the general case. Letq n × andr n × be the counting rate except for dark counts in detector and error rate of the number state |n n| with the × basis, andq n + andr n + be those with the + basis. Instead of (31) and (32), the parameters q 2 , . . . , q 2k+1 , and r 2 , . . . , r k+1 are introduced by Note that q j is the rate of counting except for dark counts. Instead of (33) -(35), we have the relations for i = 0, . . . , 2k, and for i = 1, 2, . . . , k. In addition, r i belongs to [0, 1] for i = 1, . . . , k, and q i belongs to [0, 1 − p D ] for i = 0, . . . , 2k + 1 because the dark counts occur with the probability p D . Therefore, we can estimate of the ranges of q 1 and r 1 from these conditions. Here, we defined the matrix (P j i ) i=0,...,2k,j=0,...,2k+1 defined by where the k-dimensional vectors Y and Z and the k × k matrix X are given the above. For examples, this matrix in the case of k = 1, 2, 3 is given as

Asymptotic key generation (AKG) rates
In this section, based on the expansion (30), we derive a lower bound of I( p, s). Since the matrix P j i has no inverse matrix, it is impossible to derive q 1 and r 1 from the conditions (36) and (37) uniquely. Then, in order to evaluate I( p, s), we introduce the quantities q 1,k min and b 1,k max by Thus, when µ k is the signal intensity, In the following, we calculate b 1,k max and q 1,k min as follows. Theorem 2 Define the quantities q j,min , q k+j,min , and b j,max by for j ≤ k. Then, the relations q j,min = min q=(q 0 ,...,q 2k+1 ) T hold. Therefore, q 1,k min and b 1,k max are calculated as q 1,k min = max{q 1,min , . . . , q 2k,min } (45) In order to calculate q j,min and b j,max in the channel model (6) and (7), we define the quantity ǫ j α (µ 1 , . . . , µ j ): This quantity can be characterized by the following theorem.

Comparison of AKG rates
In order to compare estimates of q 1 and q 1 r 1 , we assume that p i = p i+k . As is mentioned in Section 4, q 1,2 min and b 1,2 max give the best estimates of q 1 and q 1 r 1 among known estimates in the case of k = 2.
Next, taking into account the dark count effect, we consider the effect of increase of the number k of intensities with the forward and reverse error correction. In these comparsions, it is assumed that all detections with the vacuum pulse are caused by the dark count effect in the detector, i.e., p D = p 0 .  (4) and (5) in this graph.
In order to discuss the forward error correction case with p D = p 0 , we replace the definitions of q 2,min , b 1,min , q 1,3 min , b 1,3 min , q 1,4 min , and b 1,4 min in the above table because these depend on the value p D . The transmission rates of the five cases (2)  In order to discuss the reverse error correction case with p D = p 0 , we replace the definition of I(µ, q 1 , b 1 ) as I(µ, q 1 , b 1 ) := 1 2 (µe −µ q 1 (1 − h( b 1 q 1 )) + p D − p(µ)h(s(µ))) in the above table. The transmission rates of the five cases (2), (3.2), (3.3), (4), and (5) with the reverse error correction are given in the Fig 3. The acheivable transmission length in the forward case is (2)230.7km, (3.2)231.3km, (3.3)232.5km, (4)233.2km, (5)233.3km. That is, incereasing the number k from 2 to 3, yields increasing the acheivable transmission length increases with 1.8 km, while incereasing the number k from 3 to the infinity yields increasing it only with 0.8 km. These comparisions indicate that the AKG rate is almost saturated in the case of k = 3. Further, these graphs (Figs 1, 2, and 3) show that our AKG rate max 0.3<µ 3 I(µ 3 , q 1,3 min (0.1, 0.2, µ 3 ), b 1,3 max (0.1, 0.2, µ 3 )) is better than Wang's proposal in the case of k = 3.  The optimal signal intensity with the reverse error correction is calculated as

Generalization of the concept "difference"
In this section, in order to prove theorems given in above sections, we generalize the concept of "difference". In the conventional mean value theorem, we focus on the difference f (x 2 )−f (x 1 ) between two points x 1 and x 2 for a given real-valued function f . When we treat n+1 points x 1 , . . . , x n+1 , we introduce the quantity ∆ n+1 as a generalization of the difference. This generalized difference is characterized by n-th derivative by the following generalization of the mean value theorem: Lemma 1 Any n-differentiable function f on R and any n + 1 points x 1 < . . . < x n+1 satisfy the equation: where ∆ n is the n-th simplex, and p n is the standard uniform probability measure on ∆ n . Hence, there exists a point z ∈ [x 1 , x n+1 ] such that Its proof will be given in Appendix A. Therefore, the quantity ∆ n+1 f (x 1 , . . . , x n+1 ) can be regarded as a generalization of difference. This quantity satisfies the following recurrence formula.
Its proof will be given in Appendix B.
Using this formula, we can prove the following by induction: In particular, when 0 ≤ k ≤ n − 2, which can be also checked by Lemma 1. When k = −1, we obtain For example, when , the relation (55) yields that Since the j-th derivative of x n is 0 for n ≤ j − 1, the relation (56) holds for Using these formulas, we obtain interesting characterization of the following two k × k matrixes A k and B k : That is, where A i k is the i-th column vector of the matrix A k and B k:i is the i-th row vector of the matrix B k .
These matrixes play an important role in the proof of Theorems 1 and 2, and are the inverse matrixes of each other, i.e., which is equivalent with The equation (57) is trivial in the case of l ′ ≥ l. When l ′ < l, Applying Lemma 1 to the case of n = l − l ′ + 1 and f (x) = 1, we can show that LHS of (59) is equal to 0. Now, we modify the matrix A k as follows: That is, The inverse matrix C −1 k is characterized as follows: The equation (60) can be checked as follows. Using (55), we obtain Thus, we obtain (60).

Conclusion and further improvement
We have discussed the AKG rates with phase-randomized coherent light by the decoy method, in which the number k of possible intensities is arbitrary. For this purpose, by generalizing Wang's expansion, we have derived a convex expansion of the phaserandomized coherent state, which allows us to parameterize Eve's operation using 3k +3 parameters even in the general case. Thanks to this parameterization, lower bound of AKG rate has been obtained with k decoy intensities. Also, assuming that the noise in the quantum channel is described by a natural model, we have derived upper bounds independent of the number k of decoy intensities as the universal upper bound. It has been numerically demonstrated that the AKG rate is close to the universal upper bound in the case of k = 3, i.e., the AKG rates cannot be improved so much even if we prepare a number of intensities k larger than 3. Further, we have discussed the intensity maximizing our upper bound of the AKG rate. It has been shown that this optimum intensity is always less than 1. We have also characterized the relationship between the AKG rate and intensities that do not generate raw keys. These results indicate how to choose intensities in an implemented quantum key distribution system, in which the decoy method is applied. Unfortunately, this paper does not treat the security with the finite-length code. However, we will treat this issue in the papers [24,21], in which our expansion (30) and the matrix P k play an essential role.
For a further improvement of AKG rate, we can assume that some error happens in the generator or the detector. In this case, there is a relation r 1 ′ = (1−p S )r 1 +p S (1−r 1 ) between the observed error rate r 1 ′ in the × basis and the error rate r 1 ′ out side of the generator and the detector in the × basis, where p S is the probability that the error in the × basis occurs at generation or detection in the single-photon state. That is, it is suitable to substitute . Now, we introduce new parameters b = 1 − a 1 or 1 − a n+1 and b i = a i b (i = 2, . . . , n). Then, .