Testing the photon-number statistics of a quantum key distribution light source

A commonly held tenet is that lasers well above threshold emit photons in a coherent state, which follow a Poissonian statistics when measured in photon number. This feature is often exploited to build quantum-based random number generators or to derive the secure key rate of quantum key distribution systems. Hence the photon number distribution of the light source can directly impact the randomness and the security distilled from such devices. Here, we propose a method based on measuring correlation functions to experimentally characterise a light source's photon statistics and use it in the estimation of a quantum key distribution system's key rate. This promises to be a useful tool for the certification of quantum-related technologies.


I. INTRODUCTION
Quantum random number generators (QRNGs) [1][2][3] and quantum key distribution (QKD) [4][5][6] are the first quantum-related technology to leap out of the lab and reach the maturity necessary for the market.
The goal of a QRNG is to generate unpredictable numbers based on the laws of quantum physics.The typical example is a single, indivisible, photon impinging on a beam splitter [1,7,8], which ideally provides a uniformly distributed random bit, 0 or 1, depending on the output port it emerges from.By measuring the photons with a pair of photodetectors, it is possible to extract random strings that can be suitably post-processed and employed for cryptographic applications as well as for lotteries, gambling and scientific simulations.Other common implementations are based on events that are expected to follow a Poissonian distribution, for example the arrival time of photons emitted by a coherent light source [2,[9][10][11][12].In other cases, as for QRNGs based on the phase noise of a laser [13][14][15][16][17][18], the Poissonian nature of the source can evidence the good functioning of the randomness-generating mechanism.
The aim of QKD, on the other hand, is to generate shared randomness between two distant parties, traditionally called Alice and Bob, through the exchange of quantum signals.The most common light source in QKD systems is a laser followed by an attenuator.The attenuation level depends on the particular protocol implemented, but in general it is set so that the intensity of the emitted pulses approaches the single-photon regime, which is captured by the following condition: In Eq. (1), µ is the mean photon number of each emitted pulse and p n is the probability to emit n photons.The calibration of the mean photon number is essential to guarantee the security of QKD systems.If µ is too large, * marco.lucamarini@crl.toshiba.co.uk the secret information is redundantly encoded in n > 1 photons, allowing an eavesdropper (Eve) to access the secret information.Therefore µ has to be carefully set and it is typically monitored in real time in existing QKD systems [19].However, this is not the only security requirement and the statistics of the source, represented by p n in Eq. ( 1), also plays a crucial role.To gain some insight into this problem, consider Fig. 1, where three different types of photon number distributions are depicted, all displaying an average photon number µ = 1.Fig. 1(a) is the Poissonian case, with p n = e −µ µ n /n!, which typically stems from a laser operated well above threshold [20].The resulting distribution is composed of statistically independent events and the number of photons in each pulse varies in a random way.In the super-Poissonian case, on the other hand, represented in Fig. 1(b), the light is composed of photons bunched in the same optical pulse.The typical example is a thermal source, which features a photon number statistics equal to p n = µ n /(1 + µ) 1+n .Finally, we show in Fig. 1(c) the ideal case of a true single-photon source, which features p 1 = 1, p n =1 = 0 and displays deterministic emission of exactly one photon per pulse.
Security-wise, the case (c) is the most secure, as there never is more than one photon carrying the secret information.On the contrary, the case (b) is the least secure, because the photon bunching effect favours the redun-arXiv:1711.00440v1[quant-ph] 1 Nov 2017 dant encoding of the information.This shows that even if the mean photon number µ is well characterised, the photon number distribution probability p n can still determine the insecurity of a QKD system.It is therefore important to find ways to characterise the photon statistics p n of a quantum-related light source.
For that, one method would be to use photon number resolving detectors or, equivalently, space-multiplexed or time-multiplexed threshold detectors [21].However, in this case, a precise calibration of the detection efficiency would be required, which is far from trivial.A second, more practical, method is the one proposed here.We perform a generalised Hanbury Brown-Twiss (HBT) experiment [22,23] to estimate the normalised correlation functions g (m) (τ ) [24,25] of the light source, with m = {2, 3, 4}.These functions, evaluated at time delay τ = 0, are closely related to the photon statistics of the light source.For example, it is well-known that if p n follows a Poissonian statistics, the light source has all the g (m) (0) equal to unity [20], irrespective of the order m.Although experimentally it is not possible to measure the correlation functions to all orders, we will show that measuring them up to the fourth order is sufficient to determine tight bounds on the secure key rate of a QKD system.To apply this method, we will work in the low detection efficiency approximation.This allows us to treat our threshold detectors as linear and makes our estimate independent of the exact efficiency of the detectors used in the experiment.The details of this approach are given in the next Sec.II.In Sec.III we will describe the experimental implementation, whereas in Sec.IV we will exploit our experimental results to estimate the secure key rate of a QKD system [26].

II. NORMALISED CORRELATION FUNCTIONS
In the present work, we are interested in the characterisation of the photon statistics p n of a light source that is suitable for quantum-related technologies like a QRNG or a QKD system.We consider a scenario where the characterisation is worked out using detectors that are only loosely calibrated.In particular, we assume that the tester has guarantees that the efficiency of his detectors is positive but, at the same time, smaller than a certain threshold η 0 .This allows us to treat threshold detectors as if they were linear and makes our result independent of the exact value of the detectors' efficiency.This assumption is perfectly reasonable if the source is tested in a trusted environment, where there is no eavesdropping ongoing during the test.It is then reasonable to assume that the tester has a certain control over his detectors and can arbitrarily decrease their efficiency below η 0 using a suitable attenuator.
Under this assumption, we can write the "normalized correlation functions" [24,25], or "degrees of coherence" [20], g (2) (0), g (3) (0), g (4) (0) of a stationary light source as: where â † and â are the boson construction and destruction operators of a single quantum of light, respectively, and the bra-ket notation indicates the average operation over the states emitted by the source.
When the source is pulsed, which is often the case in QKD systems, the above relations can be written in a discrete form [27,28]: g (4) [0] = â † where k denotes the integer value representing the pulse number.In the following, we will use a compressed notation g m to indicate the quantity g (m) [0], dropping the zero time delay whenever it is unnecessary to specify it.The normalized photon correlation functions in Eqs. ( 5)-( 7) represent the quantities measured in our experiment.Under the low detection efficiency approximation they are attractive from an experimental point of view as they can be measured using threshold single photon detectors.In Sec.III we will describe the experiment performed to estimate these functions.
To connect the discrete correlation functions g m with the photon statistics p n , we consider the density matrix of the state emitted by the source, expressed as a sum of photon number states: Upon rewriting Eqs. ( 5)-( 7) in terms of the density matrix in Eq. ( 8), and using the compressed notation previously introduced, we obtain: where µ has been defined in Eq. ( 1).From the above equations, it is apparent that the value of µ is necessary to find out the n-photon probability p n .For example, it is intuitive that for large µ the p n will be centred around larger values of n.However, the exact determination of µ is not necessary for estimating the correlation functions g m .Therefore we can treat µ as a parameter of the theory, known to the user Alice who prepares the states for the QKD protocol.It is worth noting that µ is usually determined quite easily and precisely in existing QKD systems by means of a calibrated power meter.

III. EXPERIMENTAL SETUP AND RESULTS
A schematic of the experimental arrangement used to characterise the g m functions is shown in Fig. 2. In order to connect our method with a prominent application, the light source to be characterised has been taken directly from a GHz-clocked QKD prototype [19].
It consists of a DFB laser diode linked to a 1 GHz pulse generator and driven by a DC current source.A small DC current is supplied to the laser diode from the DC current source, which causes a minor amount of spontaneous emission.The pulse generator produces a square wave with approximately 2V amplitude which modulates the DFB laser diode on top of the DC current.Under normal operating conditions, the modulation drives the laser diode in a gain switched mode, i.e., above and below lasing threshold.This causes a train of phaserandomised optical pulses to be produced at a frequency of 1GHz [18].The resulting optical pulses are then attenuated to the single photon level using an electronically controlled attenuator (ECA), so that the condition expressed by Eq. ( 1) is satisfied in our experiments.
To analyse the correlation functions, a four channel HBT interferometer was constructed.Photons from the QKD laser source under test are transmitted into a 50:50 beamsplitter whereupon they encounter a second 50:50 beamsplitter.Four self-differencing avalanche photodiodes (SD-APD) [29] synchronised to the QKD laser source detect the transmitted photons.Care was taken in the experiment to guarantee a detection efficiency well below 1% in each detector.The subsequent electrical output from each SD-APD is then discriminated before being sent to a correlator card where time stamps are assigned to the photon arrivals.The time stamps from the correlator were analysed by a custom built program.The program constructed two, three and four-fold photon coincidences from the photon arrival times.From this data, the resulting correlation functions up to the fourth order were evaluated.

A. Dependence on laser drive current
We start our analysis by examining the effect of changing the laser drive current on the normalised second order correlation function g 2 .This allows us to gain some insight on how the photon statistics might change for a laser diode under different operating conditions.Fig. 3(a) shows the DC laser current vs optical output power dependence for the QKD laser source.The 1GHz square wave modulation is applied at all times.Note the optical power is measured before the ECA using a standard optical power meter.As expected, the optical output rises sharply after the lasing threshold, around 4.5 mA, has been exceeded.Fig. 3(b) (points) shows the resulting normalised second order correlation function minus one (g 2 − 1) as measured in the single photon regime (i.e. after the ECA).Below lasing threshold, g 2 − 1 ≈ 1, indicating photon bunching, as might be expected for a thermal-like source.Around threshold, g 2 − 1 starts to reduce.As the DC current is further increased, the fluctuations begin to die out and g 2 − 1 tends to 0, indicating the laser emission is approaching Poissonian-like photon statistics, for which g 2 = 1.
To gain more insight into the functional dependence of g 2 − 1 on laser intensity, we fit the experimental data in Fig. 3(b) with a non-linear oscillator model [? ], using a single free parameter, the laser threshold intensity I 0 .This model has been successful in describing g 2 around lasing threshold for gas lasers [? ] provided the lasers under study can be approximated as single mode.The fit gives I 0 = 12.67 µW corresponding to a DC current of 4.49 mA (solid curve in Fig. 3(b)).The fit appears reasonable for the experimental points around I = I 0 and up to I = 6.Beyond this laser intensity, the model overestimates g 2 − 1.Alternatively, we can extract a rough indication of lasing threshold I 0 from Fig. 3(a The above analysis suggests the underlying model of the non-linear oscillator qualitatively describes the experimental data, if not quantitatively.Most likely the ideal model lies between the two curves in Fig. 3(b).This indicates the excess g 2 − 1 for the experimental points with the three highest intensities is due to additional fluctuations on top of the noise predicted by the model.Note, as we show below, we do not rely on any physical model to construct bounds on the secure bit rate, see Sec.IV.
For use in QKD, it is important to choose an appropriate DC current level.Too low gives large photon bunching and therefore a high incidence of multi-photon events, which greatly facilitates photon number splitting attacks by an eavesdropper.Too high a DC current yields Poissonian behaviour but increases the laser CW optical background.A high CW optical background increases the phase correlation between successive optical pulses and destroys the phase randomisation of the QKD laser source [18].
We choose a compromise such that under normal operating conditions, i.e. during QKD, the DC current is set to approximately 6.5 mA which corresponds to g 2 − 1 of the order of 10 −3 , as shown in Fig. 3(b).
Next we characterise higher order correlations.By way of comparison higher order correlations are measured for two cases.Firstly, below lasing threshold at a DC current of about 4.0 mA.Secondly, above lasing threshold, where the DC current is 6.5 mA corresponding to normal operating conditions.

Variation of gm for laser below threshold
We first examine the variation of the higher order correlation functions when the laser is operated below threshold.A DC current of 4.0 mA is selected.functions g 2 , g 3 and g 4 and the numerical values are listed in the following equations: The source displays a high level of bunching, approaching a thermal distribution, hence we term it "quasi-thermal".The peak at time zero grows with the order m, as shown with red bars in Fig. 4(b).Also plotted are the predictions of the non-linear oscillator model for the two laser threshold intensities shown in Fig. 3(b).Both scenarios predict slightly stronger correlations for all orders compared to what was experimentally measured.This behaviour can be intuitively understood by assuming the source is mainly bunched but also possesses a smaller Poissonian component since the measured g 2 = 1.6985 is less than the value g 2 = 2 expected for a purely thermal source.The bunching contribution naturally gives stronger correlation values g m as m is increased but at the same time the smaller Poissonian component contributes to reducing these correlation values.

Variation of gm for laser above threshold
We now turn our attention to the variation of the higher order correlation functions when the laser is operated above threshold.A DC current of 6.5 mA is selected for this purpose.Fig. 5(a) shows the measured normalised correlation functions g 2 , g 3 and g 4 and the numerical values are: This time there is no clear peak at time zero, which indicates the laser is qualitatively Poissonian.However, plotting the g m as a function of m reveals an ever so slight increase in correlation despite the large error bars, Again we plot the predictions of the non-linear oscillator model for the two laser threshold intensities shown in Fig. 3(b).For g 2 and g 3 the measured data is between the two predictions.For g 4 both predictions slightly underestimate g 4 .The reason for this is unclear but the error bar for the g 4 measured point is almost an order of magnitude greater than the error bars for g 2 and g 3 .This stems from the four-fold coincidence count rate being considerably lower than three or two-fold coincidence count rates.Furthermore, g 4 − 1 is quite small, < 0.1.This can be compared with g 4 − 1 in Fig. 4(a), which is around a factor of hundred larger.
Note for the measured correlation functions we have also considered corrections from afterpulsing, dark counts and dead time.However these corrections are well within the size of one error bar, so have consequently been safely neglected in the current analysis.

IV. WORST-CASE BOUNDS FROM CORRELATION FUNCTIONS AND SECURE KEY RATE
Let us now consider how the experimental results from the previous section can be used to determine the secure key rate of a QKD system.For that, we rewrite the key rate equation for the efficient [30,31] decoy-state [32][33][34] BB84 protocol [4] in the finite-size scenario [26,35,36] in the following way: With minor modifications, Eq. ( 18) can be easily adapted to other BB84-like protocols.To write it, we considered a protocol where Alice prepares phase-randomised weak coherent states in the two bases Z and X, randomly chosen with probabilities p Z and p X = 1 − p Z , respectively, and in the three intensity classes u (signal), v (decoy) and w (vacuum), randomly chosen with probabilities p w , p v and p u = 1−p v −p w , respectively.However, the users extract the secure key bits only from the basis Z and from the class u.The quantity R is the fraction of secure bits per signal distilled by the system; when multiplied by the system's clock rate, it provides the amount of secure bits in the time unit; e 1,X is the upper bound to the singlephoton QBER in the X basis and y 1,Z is the lower bound to the 1-photon yield in the Z basis.This last quantity is multiplied by the lower bound to the 1-photon probability, p 1 , thus providing an overall lower bound to the 1-photon gain of the protocol.Finally, f is the inefficiency of error correction, Q Z and E Z are the gain and the QBER, respectively, measured in the Z basis and ∆ is a parameter related to the finite-size effect.
The crucial difference between the rate in Eq. ( 18) and previous QKD key rates resides in the term p 1 .Usually, a Poissonian distribution is assumed for decoy-state QKD, entailing that p 1 = µe −µ [32][33][34].However, by applying our method, we can now drop this assumption and replace it with the lower bound p 1 that is measured in the experimental characterisation of the light source.Similar bounds can be obtained for the other probabilities p 0 , p ≥2 , as well as for the other parameters y 0 , y 1 , e 1 usually present in decoy-state QKD.The details on the constrained numerical optimization providing such bounds are given in the Appendix.
As an example, in Fig. 6 we provide the bounds for p 0 , p 1 and p 2 drawn from the experimental values reported in Eqs. ( 15)-( 17) and we compare them with the corresponding probabilities from the Poissonian distributions, for various experimental conditions.On the topleft diagram, we consider the case where Alice prepares a mean photon number µ = 0.42, a common value for the signal state in decoy-state QKD, and uses only the g 2 function from Eq. ( 15) as a constraint.In this case, the comparison with the ideal Poissonian probabilities normally used in decoy-state QKD shows that the bounds are quite loose.Therefore, we expect a dramatic reduction of the key rate obtained by assuming that the distribution is Poissonian.However, if more constraints are added, the bounds become tighter and tighter, as can be noted from the top-right and bottom-left diagrams in Fig. 6.This should entail a close-to-ideal key rate in these cases.Finally, when the mean photon number is FIG. 7. Secure key rate of the efficient decoy-state BB84 protocol in the finite-size scenario, for the ideal case of a Poissonian source (black line) and for the real case where the correlation functions up to the fourth order have been measured.The values of the correlation functions have all been drawn from the experiments.The blue, green and wine lines correspond to the source described in Sec.III A 2, Eqs. ( 15)-( 17), whereas the red line is for the source in Sec.III A 1, Eqs. ( 12)- (14).The bounds for this latter source, which displays a quasi-thermal distribution, are given in the inset, together with the Poisson probabilities for comparative purposes.The mean photon number for the signal and the decoy states in the protocol have been set to u = 0.42 and v = 0.02, respectively, and the security parameter is < 10 −10 .
decreased to µ = 0.1, a value commonly used for decoy states in QKD, and all the constraints are applied, the match with the Poissonian probabilities is nearly perfect (see bottom-right diagram in Fig. 6), implying a negligible reduction in the key rate.
In the estimation leading to Fig. 6, we have considered that the experimental errors follow a Gaussian distribution and we have taken γ = 7 standard deviations to define the related confidence interval.This corresponds to a probability = 2.56 × 10 −12 that the true experimental value falls outside the confidence interval.In turn, once all the confidence intervals related to the constrained optimisation are taken into account [37], this probability generates a global security parameter for the QKD protocol in the finite-size scenario [38,39] well below 10 −10 .
As a prominent application of our method, we use the bounds just derived to estimate the secure key rate of a QKD system according to Eq. (18).The results are shown in Fig. 7.It can be seen that when the g 2 function is the only constraint in the problem (wine-coloured line), there is a considerable decrease in key rate and maximum distance achieved by protocol.On the contrary, when also the g 3 and g 4 are taken into account (green and blue lines), the key rate is very close to the ideal one.Intuitively, this can be explained by the fact that each additional constraint improves the bounds on the p n .When only g 2 is considered, we can only reliably bound the first term in the photon statistics, p 0 , which does not enter the key rate equation (18), or it plays only a minor role in other versions of the BB84 protocol [36].However, when both g 2 and g 3 are considered, we can bound p 0 and p 1 , the latter being the most important term in the secure key rate.This explains the large gap in the plotted key rates between the orange line (g 2 only) and the red line (g 2 and g 3 together).
In the same figure, we also plot in red colour the key rate corresponding to the quasi-thermal source described in Sec.III A 1, which is obtained by plugging Eqs. ( 12)-( 14) into the corresponding optimisation problem.The upper and lower bounds are shown in the inset of Fig. 7.It can be noticed that the lower bound for the zero-photon fraction of this distribution (leftmost yellow bar) is higher than the corresponding Poisson probability (brown bar), clearly showing that the distribution is not Poissonian.Nevertheless, it is still possible to extract from it a key rate that remains positive up to a distance of about 40 km in optical fibre.Larger distances might be achieved upon optimisation of the system's parameters.This positive key rate is remarkable given the highly bunched nature of the source and represents an experimental evidence of a positive key rate drawn from a quasi-thermal source.
On the other hand, this result should be treated as a caveat for the correct operation of a laser source.If the laser is unawarely operated close-to-threshold, the users will not realise that the correct key rate is the one corresponding to a quasi-thermal source (red line in Fig. 7) and will keep using the privacy amplification corresponding to a Poissonian source.This entails that less privacy amplification than necessary will be performed, compromising the security of the system.The calibration of the QKD system's light source and the detailed description of its operative parameters represent a desirable best practice to prevent these kind of security risks.

V. CONCLUSION
In quantum-related technology, claims about randomness and security crucially depend on a set of assumptions that have to be carefully met in the implementation.In the present work, we introduced and experimentally realised a test for one of the most common assumptions in QRNG devices and QKD systems, i.e., the Poissonian nature of the light source.
The method is based on a generalised HBT experiment [22] with single-photon detectors in the lowefficiency approximation.If the detection efficiency is known to be small, the exact efficiency of each detector simplifies in the corresponding expressions of the normalised correlation functions, Eqs. ( 5)- (7).This lets us bypass an accurate calibration of the detectors, which is often a cumbersome process.
On the other hand, our method makes it possible to test a light source upon which the tester has some limited prior knowledge.For example, we showed that the source in this work qualitatively follows a well-established laser model.However, the technique does not make use of this model and makes no stringent assumptions on the physics underlying the photon source under test.This identifies a practical way to undertake calibration tests, which are normally performed in the protected environment of a laboratory, where the assumptions on the potential presence of a malevolent agent are mild.
The results of the experimental tests, given in Eqs. ( 12)-( 14) for a laser operated close to threshold and in Eqs. ( 15)-( 17) for a laser above threshold, show that photon statistics of the light source crucially depend on the driving current.When the laser is operated above threshold, the photon distribution is very close to Poissonian, which is especially true when the mean photon number is low.In this case, the normalised correlation functions up to the fourth order are all close to 1.This implies that the bounds to the Poissonian distribution, depicted in Fig. 6, are quite tight.Nevertheless, when plugged in a typical key rate equation for a QKD system, these bounds generate a key rate that can deviate from the ideal one.This is shown in Fig. 7, which suggests that measuring at least the g 2 and g 3 correlation functions is essential to avoid a major reduction of the secure key rate.This is feasible and practical, as we have shown, thus promoting the present method as a useful tool to guarantee the security of QKD systems.Close to threshold, the distribution of the emitted becomes quasi-thermal and displays a highly bunched behaviour.However, even without optimising the system's parameters, we could obtain a positive key rate over a few tens of kilometres of optical fibre.bottom-right diagram in Fig. 6, where the mean photon number is 0.1, i.e., 5 times bigger than the value considered here).
The constraints in Eqs.(B1)-(B3) define an optimisation problem that is nonlinear in the unknowns.This obstacle can be overcome by replacing the unknowns p(µ) n with their bounds p (µ) n and p (µ) n , drawn from the optimisation problem in Eq. (A1) [41].Moreover, the sums contain an infinite number of terms, which is not suitable for a practical scenario.This can be easily circumvented using the same technique demonstrated for the problem in Sec.A and in the literature [42], i.e., by cutting the sums to a finite value N cut and then bounding the residual.As a result, we obtain the following constraints: where we have indicated the upper and lower bounds to the experimental yields for the mean photon number µ as Y (µ) and Y (µ) , respectively, and the upper bound for the residual sum as Γ µ .

FIG. 1 .
FIG. 1. Schematic of photon number distributions in a train of optical pulses with average photon number µ = 1.(a) Poissonian source; (b) super-Poissonian source; (c) single-photon source.Each coloured dot represents a single photon.Vertical lines identify the optical pulses in the same time slot.

FIG. 2 .
FIG.2.Schematic of the experimental arrangement for measuring normalized correlation functions up to the fourth order.ECA: electronically controlled attenuator, SD-APD: selfdifferencing avalanche photodiode.

FIG. 3 .
FIG. 3. (a): Laser DC current vs optical output power from DFB laser diode.Above threshold, around 4.5 mA, there is a steep rise in output power.(b): g2 − 1 as a function of laser optical intensity I, normalised to threshold optical intensity I0 with I0 = 12.67µW.Points: Experimental data.Solid curve: fit of the non-linear oscillator model using I0 as a free parameter.Dashed curve: non-linear oscillator model with I0 = 5.38µW.

Fig. 4 (FIG. 5 .
FIG. 5. (a): Measured normalised correlation functions g2, g3 and g4 for a DC current of 6.5 mA.(b): n-th order normalised correlation functions as a function of the order n.Red bars: experimental data; blue bars: prediction of non-linear oscillator model when I0 = 5.38 µW; wine bars: prediction of non-linear oscillator model when I0 = 12.67 µW.

Fig. 5 (
Fig. 5(b), red bars.On this occasion the source appears to have a minor bunched-like component, hence much weaker correlation values g m are observed as the order m is increased, compared to the previous case.Again we plot the predictions of the non-linear oscillator model for the two laser threshold intensities shown in Fig.3(b).For g 2 and g 3 the measured data is between the two predictions.For g 4 both predictions slightly underestimate g 4 .The reason for this is unclear but the error bar for the g 4 measured point is almost an order of magnitude greater than the error bars for g 2 and g 3 .This stems from the four-fold coincidence count rate being considerably lower than three or two-fold coincidence count rates.Furthermore, g 4 − 1 is quite small, < 0.1.This can be compared with g 4 − 1 in Fig.4(a), which is around a factor of hundred larger.Note for the measured correlation functions we have also considered corrections from afterpulsing, dark counts and dead time.However these corrections are well within the size of one error bar, so have consequently been safely neglected in the current analysis.