Security of decoy-state quantum key distribution with correlated intensity fluctuations

One of the most prominent techniques to enhance the performance of practical quantum key distribution (QKD) systems with laser sources is the decoy-state method. Current decoy-state QKD setups operate at GHz repetition rates, a regime where memory effects in the modulators and electronics that control them create correlations between the intensities of the emitted pulses. This translates into information leakage about the selected intensities, which cripples a crucial premise of the decoy-state method, thus invalidating the use of standard security analyses. To overcome this problem, a novel security proof that exploits the Cauchy-Schwarz constraint has been introduced recently. Its main drawback is, however, that the achievable key rate is significantly lower than that of the ideal scenario without intensity correlations. Here, we improve this security proof technique by combining it with a fine-grained decoy-state analysis, which can deliver a tight estimation of the relevant parameters that determine the secret key rate. This results in a notable performance enhancement, being now the attainable distance double than that of previous analyses for certain parameter regimes. Also, we show that when the probability density function of the intensity fluctuations, conditioned on the current and previous intensity choices, is known, our approach provides a key rate very similar to the ideal scenario, which highlights the importance of an accurate experimental characterization of the correlations.


I. INTRODUCTION
Quantum key distribution (QKD) offers a way to distribute a secret key over the distance between two communicating parties, Alice and Bob [1][2][3]. When used in conjunction with the one-time-pad encryption scheme [4], QKD allows for information-theoretically secure communications, regardless of the future evolution of classical or quantum computers. This is so because its security is based on the laws of quantum mechanics and does not rely on computational assumptions. In recent years, QKD has progressed very rapidly both in theory and in practice, turning into a flourishing commercial technology that is being deployed in metropolitan and intercity fiber-based networks worldwide [5][6][7][8], including also satellite links [9][10][11][12].
One of the most successful QKD protocols for longdistance transmission is undoubtedly the decoy-state QKD scheme [13][14][15]. It uses phase-randomized weak coherent pulses (PRWCPs) emitted by laser sources to provide a secret key rate that scales linearly with the channel transmission [16]. Since its first theoretical proposal, numerous experimental implementations have been reported in recent years [17][18][19][20][21][22][23][24][25], being the actual distance record over fiber 421 km [26]. Decoy-state QKD has also been demonstrated over satellite links [9,11], and various companies currently implement it in their commercial products. Importantly, the use of decoy-states is also an essential ingredient for other QKD schemes that use laser sources, like e.g., measurement-devide-independent * xsixto@com.uvigo.es † vzapatero@com.uvigo.es (MDI) QKD [27] or twin-field (TF) QKD [28], the latter presently holding a distance record over fiber of 833 km [29]. However, despite these tremendous achievements, there are still certain challenges that need to be overcome for QKD to become a widely used technology. One of these challenges is to increase the secret key rate delivered by current experimental setups, which is severely affected by the limited transmissivity of single photons in optical fibers, as well as by the dead-time of the detectors at the receiver. Possible approaches for this include the use of multiplexing techniques-like e.g. wavelength division multiplexing-that simultaneously transmit several QKD channels over the same fiber [30,31], the adoption of high-dimensional QKD, which can encode many bits of information on a single-photon [32], and the increase of the pulse repetition rate of the sources [33]. Indeed, current decoy-state QKD experimental setups operate at GHz repetition rates [21,[24][25][26], and the trend is to increase their clock frequency even further. However, in such high-speed regime, memory effects in the modulators and electronics that control them create correlations between the generated optical pulses [33,34]. That is, the state of a quantum signal emitted by the source at a certain time instant depends not only on the state preparation settings selected by Alice in that time instant, but also on those selected by her in previous time instants. Importantly, if this effect is not properly taken into account in the security proof of QKD, it might open a security loophole in the form of information leakage [35]. This is so because the settings chosen to encode each quantum signal are also partially leaked through the quantum states of subsequent signals.
So far, most security proofs of QKD have neglected the effect of pulse correlations and assume independent and identically distributed emitted pulses. Therefore, they cannot be used to guarantee the security of high-speed QKD implementations. Only a few recent works partially address this problem. Precisely, the authors of [36,37] studied the case of setting-choice-independent pulse correlations, in which the emitted pulses can be arbitrary correlated between them, but these correlations do not depend on the state preparation setting choices. This scenario may occur, for example, when the temperature of Alice's laser drifts slowly over time due to thermal effects, or when her modulators' power supply fluctuates in time. More recently, various results that address the problem of setting choice-dependent pulse correlations have also been reported. In particular, [35] introduced a security proof that can handle arbitrary correlations that originate from the phase modulator that encodes the bit and basis information of each generated signal. Likewise, a restricted class of nearest-neighbour pulse correlations that arise from the intensity modulator used to prepare decoy states has been studied in [34], where a post-processing technique to treat these correlations is also provided. Also, the authors of [38] introduced a general methodology to treat arbitrary intensity correlations in a decoystate QKD setup. For this, they exploit a fundamental constraint that is a direct consequence of the Cauchy-Schwarz (CS) inequality in Hilbert spaces [35,39]. The case of correlations between the global phases of the coherent states emitted by a laser when is operated under gain-switching conditions has been studied in [40].
A main limitation of the security proof technique presented in [38] is, however, that the delivered secret key rate is significantly lower than that of the ideal scenario without correlations. Only when the intensity correlations are very tiny, the resulting performance can approximate the ideal scenario. In this paper, we improve the general methodology introduced in [38], based on the CS constraint, by combining it with a fine-grained decoystate analysis. In doing so, we achieve a much tighter estimation procedure to determine the relevant parameters that enter the secret key rate formula. This results in a notably enhancement of the achievable secret key rate in the presence of intensity correlations. Indeed, for certain parameter regimes, the maximum attainable distance can be double than that of [38]. In addition, we show that when the probability density function of the intensity fluctuations, conditioned on the current and previous intensity setting choices, is known, our analysis can provide a secret key rate that is very close to that of the ideal scenario. This highlights the importance of properly characterizing intensity correlations in QKD experimental setups. Most importantly, our results can be readily applied to guarantee the security of high-speed decoy-state QKD realizations with current technology, without much penalization on their secret key rate.
The paper is organized as follows. In Sec. II we present the problem of intensity correlations, and explain why the standard decoy-state analysis cannot be directly applied to this scenario. Here, we also emphasize our main con-tributions in relation to previous results. Next, in Sec. III, we provide the main assumptions that we consider in the security analysis. Then, in Sec. IV, we introduce a fine-grained decoy-state parameter estimation method that can be used in the presence of intensity correlations to tightly estimate the relevant parameters that determine the secret key rate. This estimation procedure is then applied to two different scenarios in the following two sections. Precisely, in Sec. V, we consider the general case where Alice and Bob only know the intervals where the intensity fluctuations lie. Coming next, the results for the case where they know the probability density function of the intensity fluctuations are presented in Sec. VI, followed by the conclusions in Sec. VII. To improve the readability of the paper, certain derivations are omitted in the main text and are included in a series of appendices.

II. RELATION WITH PREVIOUS WORK
As already mentioned, intensity correlations refer to the fact that the intensity of an optical pulse generated in a certain round of a QKD protocol depends not only on the intensity setting selected for that round, but also on the intensity settings selected for previous rounds. This effect has been experimentally quantified in [33,34]. It is illustrated in Fig. 1 with a simple example. In general, an eavesdropper (Eve) could exploit these correlations to learn information about the intensity settings selected in previous rounds by measuring the intensity of the pulses emitted in subsequent rounds. This could allow her to make the photon-number detection statistics (i.e., the yields and error rates associated to n-photon pulses) dependent on the intensity setting selected, in so breaking the elementary security premise of the standard decoystate method [13][14][15]. This situation is similar to that of a Trojan-horse attack [41][42][43], in which Eve can learn partial information about the intensities selected to generate the different pulses.
To overcome this problem, there are two main complementary approaches. Precisely, one could implement hardware countermeasures to try to reduce (or even eliminate) the correlations, and one could also include their effect in the security proof of QKD. An example of the first type of approach could be to use various intensity modulators, each of them acting on different nonconsecutive signals, such that their modulation rate effectively decreases to a regime where no correlations are created. An example of the second approach has been recently provided in [34], where the authors studied a restricted type of nearest-neighbour intensity correlations (i.e., this refers to the fact that the actual intensity of each pulse only depends on the intensity settings selected for that time instant and for the previous time instant). In particular, [34] analyzes the case where only certain intensity settings have a significant effect on the actual intensity of a subsequent pulse, and introduces a post-processing technique to guarantee security in this scenario. In addition, a preliminary experimental characterization of the probability density function of the correlations is provided, which seems to indicate that it essentially follows a gaussian-shaped form, though this knowledge is not exploited by the post-processing technique introduced.
More recently, in [38], the authors presented an asymptotic security proof that can deal with intensity correlations of arbitrary range (i.e., not necessarily nearestneighbour pulse correlations) in the decoy-state parameter estimation procedure. Their method is rather general and can be applied when all previous intensity settings can influence the actual intensity of subsequent pulses. The main idea in [38] is to pose a restriction on the maximum bias that Eve can induce between the n-photon yields and error rates associated to different intensity settings, by using the CS constraint [35,39]. Notably, only two parameters, the correlation range ξ and the maximum deviation δ max between the physical intensity and the selected intensity setting, are needed. However, despite the fact that this latter work is quite simple and experimental-friendly, it treats the deviation of each pulse independently of the previous sequence of selected intensities, but instead it takes the maximum possible deviation for the worst-case scenario, thus providing loose bounds. Indeed, the resulting bounds on the secret key rate are significantly lower than those obtained in the absence of correlations. Therefore, the question arises of whether this damaging effect on the key rate is a fundamental feature of the correlations, or rather an artifact of a loose parameter estimation procedure. Illustration of the effect of intensity correlations. The upper subfigure shows a train of optical pulses sent by Alice to the channel with three different intensity values selected at random, which are indicated with the different colours and amplitudes of the pulses. In the absence of intensity correlations, the intensity of each signal is determined only by the intensity setting selected in the corresponding time instant. In the lower subfigure, we show a specific example of nearest-neighbour intensity correlations. Solid lines indicate the pulses actually emitted, while dashed lines indicate the selected intensity settings. Precisely, this example assumes for illustration purposes that when the intensity setting of the previous pulse is lower (higher) than that of the actual pulse, then the actual intensity generated is a bit lower (higher) than that indicated by the setting selected.
In this paper we preserve the essence of the approach introduced in [38] to deal with arbitrary intensity correlations-i.e., we exploit the CS constraint to quantify the maximum bias that Eve may induce between the photon-number detection statistics associated to different intensity settings-but sharpen the decoy-state method with a finer-grained analysis of the yields and errors that keeps explicit track of the record of settings. This, in turn, allows to incorporate finer-tuned CS constraints to the parameter estimation procedure. Putting it all together, our approach enables a noticeable enhancement of the secret key rate when compared to the results in [38].
In addition, we show that, when the probability density function of the intensity fluctuations is known, our fine-grained decoy-state analysis is rather tight, as it can produce a secret key rate that is comparable to that obtainable in the absence of intensity correlations. This is explicitly illustrated by considering a gaussian-shaped probability density function for the correlations, following the preliminary results in [34] (see also [44]). In this scenario, the principal improvement mainly comes from the fact that the knowledge of the probability density function now permits to calculate certain quantitiesthat are necessary to determine the secret key rateprecisely, avoiding the need to use looser bounds that exploit monotonicity arguments.

III. ASSUMPTIONS
For concreteness, below we shall consider a typical polarization encoding decoy-state BB84 protocol with three intensity settings. Nevertheless, our results apply to other encoding schemes, and can be straightforwardly adapted to other decoy-state based QKD protocols with a different number of intensity settings.
In the first place, let us fix the notation that is used throughout the paper. In each round k of the protocol, with k = 1, ..., N , Alice selects an intensity setting a k ∈ A = {µ, ν, ω} with probability p a k , a basis x k ∈ B = {X, Z} with probability q x k and a uniform raw key bit r k ∈ Z 2 = {0, 1}. Without loss of generality, we impose that the intensity settings satisfy µ > ν > ω ≥ 0. Then, she encodes the BB84 state defined by x k and r k in a PRWCP with intensity setting a k , and sends it to Bob through the quantum channel. Importantly, the actual mean photon-number of the pulse might not match the setting a k due to the presence of intensity correlations. To finish with, we assume perfect phase randomization, perfect polarization encoding, and that there are no sidechannels beyond intensity correlations for simplicity.
On the other hand, Bob selects a basis y k ∈ B with probability q y k and performs a measurement described by a positive operator-valued measure (POVM) {M y k ,s k B k } s k ∈{0,1,f } on the incident pulse. Here, B k denotes Bob's k-th incoming pulse, s k stands for Bob's classical outcome and f stands for "no click". As usual, the basis-independent detection efficiency condition is as-sumed, i.e.,M Z,f B k =M X,f B k , and thus we shall simply denote these two operators byM f B k . This assumption could be removed by the use of MDI-QKD [27] or TF-QKD [28].
Let a k = a k , a k−1 , ..., a 1 symbolize the record of intensity settings up to round k, where a j ∈ A for every j, and let α k stand for the actual intensity emitted in round k. We consider that α k is a continuous random variable whose probability density function is fixed by the record of settings a k . From now on, we denote such correlation function as g a k (α k ). Three additional elementary assumptions of our work are listed below. Assumption 1. The presence of correlations does not compromise the poissonian character of the photonnumber statistics of the source conditioned on the value of the actual intensity, α k . Mathematically, this amounts to saying that for any given round k, and for all n k ∈ N, Notably, this assumption is supported by recent highspeed QKD experiments [33,34]. Still, we remark that our analysis could be easily adapted to consider any other photon-number statistics conditioned on the value of α k . Assumption 2. The intensity correlations have a finite range ξ, meaning that the value of the physical intensity of round k, α k , is not affected by those previous settings a k are the relative deviations. Note that, in virtue of Assumption 2, δ + a k and δ − a k only depend on a k and the previous ξ intensity settings.
From these three assumptions, it follows that the photon-number statistics for a given round k and a given record of settings a k are (2) Note that assumption 2 is not explicitly imposed here.

IV. QUANTIFYING THE EFFECT OF INTENSITY CORRELATIONS ON THE DECOY-STATE PARAMETER ESTIMATION PROCEDURE
With the three assumptions stated above and to account for the influence of intensity correlations in the decoy-state analysis, we use the fundamental CS constraint [35,39]. This result poses a natural constraint between the measurement statistics of two non-orthogonal states. Hence, in particular, one can use it to restrict the possible bias that Eve may induce between the yields (error probabilities) associated to different records of settings. The reader is directed to Appendix A for a definition of this statement. The limits we derive with it are shown below.
We define both the yield and the error probability, for any given round k, photon-number n ∈ N and record of settings v 0 , ..., v ξ ∈ A as: (3) In virtue of the CS constraint, for any two distinct intensity settings v 0 and w 0 that could be selected in round k, and for any record of previous settings v 1 , ..., v ξ , in Appendix A it is shown that the associated yields and error probabilities satisfy and where with the function g ± (y, n,v0...v ξ ,r . Crucially, the parameters τ ξ v0w0...v ξ ,n are a lower bound on the squared overlap of the two quantum states with which the two yields are calculated. Explicit expressions for the parameters τ ξ v0w0...v ξ ,n are derived in the Appendix B 2 for the two different scenarios that we consider in Secs. V and VI. To enable the use of linear programming for the decoystate parameter estimation procedure, a linear version of the constraints given by Eq. (4) is needed. As shown in Appendix A, the linearized versions of the CS constraints for the yields and error probabilities can be expressed as and where bothỸ n,v0...v ξ are the reference parameters of the linear approximations, introduced in Appendix C. Finally, the functions G ± are defined as with g ± (y, In this section we now consider a general scenario in which the correlation function g a k (α k ) given by Eq. (2) is unknown. This implies that one cannot compute the photon-number statistics of the emitted pulses explicitly, and thus must impose bounds on them by invoking monotonicity arguments. We consider two cases.
In the first one, the maximum relative deviations δ + (a k ...a k−ξ ) and δ − (a k ...a k−ξ ) depend on the intensity setting choices corresponding to all the previous ξ rounds and the present round, and define an interval that is in general not symmetric with respect to the value of the selected setting. As shown in Appendix B 2, and denoting a N 1 = a 1 . . . a N , the bound for this case reads where the terms a ) and analogously for a (w0)± i . Note that, with this notation, settings a 1 to a k are fixed to a certain value, while settings a k+1 to a i are not. Also, note that the bound given by Eq. (10) does not depend on the photon-number n, i.e. it holds for all n.
Secondly, we consider as well a simplified situation where only a worst-case deviation δ max is known for every possible record of settings. That is, it holds that In this case, the bound simply reads where a ± i = a i (1 ± δ max ). This derivation is also presented in Appendix B 2. Note that Eq. (10) can only be used if one is able to experimentally characterize the maximum relative deviation for every possible combination of settings.
From now on, and in order to keep the discussion simple, we focus on the case of nearest-neighbour intensity correlations. That is, we set ξ = 1. For this purpose, it is convenient to do a slight change of notation, by calling v 0 = a; w 0 = b and v 1 = c as this makes the following section easier to follow. Note that the analysis below can be easily adapted for the case of ξ > 1, and we include simulations for this latter scenario in Sec. V D.

B. Decoy-state method
To introduce the linear programs that perform the parameter estimation, we shall now provide a decoy-state analysis.
Let us start by defining the Z basis gain for a certain pair of intensity settings a and c, and a number of rounds N , as a,c = 1 if in round k, Alice selects an intensity setting a that is preceded by an intensity setting c (in round k − 1), both Alice and Bob select the Z basis, and a click occurs at Bob's side. Thus: Straightforward monotonicity arguments lead to the following bounds on the photon-number statistics p (k) (n|a, c): where a ± c = a(1 ± δ ± (a k =a,a k−1 =c) ). Note that, despite the fact that this notation is similar to that used in Eq. (10), the parameters a ± c are actually not equal to a (γ)± i with γ ∈ A. This is so because here, only the settings a k and a k−1 matter, and they are respectively fixed to a and c. Now, using these intervals in Eq. (13), one obtains for all a, c ∈ A and k = 1, ..., N . Selecting a threshold photon-number for the numerics, n cut , and using the fact that we have that n,a,c .
Notice how replacing Z by X everywhere, one obtains the corresponding analysis for the X basis gains in a certain round k.
We now have to impose similar constraints to the error counts. For that matter, we define the number of X basis error counts with settings a in round k and c in round a,c I {r k =s k } . Then, we have that where we have defined With the same steps as before, it follows that Now, summing over k and dividing by N in both Eq. (17) and Eq. (21), one obtains bounds for the average parameters from the round dependent bounds. Thus, defininḡ we obtain that the final bounds are: C. Linear programs for parameter estimation In this section, we present the linear programs that allow to estimate the relevant single-photon parameters by putting together the decoy-state constraints, already introduced above, and the linearized CS constraints with the notation presented in Eq. (A6).
To begin with, note that the quantities that go into the secret key rate are the average number of single-photon counts associated to the case where Alice selects the signal intensity setting, which we shall denote byZ 1,µ,N , and the average number of phase errors associated to these single-photon counts (which in the case of the BB84 protocol, without state preparation flaws and assuming the asymptotic limit, match the average number of singlephoton error counts in the X basis). Moreover, since below we shall consider the asymptotic secret key regime where p µ ≈ 1, for simplicity we can restrict our estimation to the number of single-photon error counts in the X basis when Alice selects the intensity setting, which we shall denote byĒ 1,µ,N .
Let us considerZ 1,µ,N first. Formally, where for each k, Z is the probability that Alice selects the signal intensity setting µ, both Alice and Bob select the Z basis, and a click occurs at Bob's side. In turn, it is clear that Z (k) 1,µ decomposes as where of course Z Then it trivially follows that: Note that the above equation lower bounds the expected value of the quantity of interestZ 1,µ,N using the yields y 1,a,b,N introduced in the previous section. Thus, the linear program of interest is c + abc,n + m + abc,n y n,a,c,N ≥ y n,b,c,N (a, b, c ∈ A, b = a, n = 0, . . . , n cut ) , c − abc,n + m − abc,n y n,a,c,N ≤ y n,b,c,N (a, b, c ∈ A, b = a, n = 0, . . . , n cut ) , 0 ≤ y n,a,b,N ≤ 1 (a, b ∈ A, n = 0, . . . , n cut ) , with the parameters c ± abc,n and m ± abc,n that arise from the linear CS constraints given in Appendix A. We shall denote the lower bound on Z 1,µ,N obtained with the linear program above byZ L 1,µ,N . Naturally, substituting Z for X everywhere yields the equivalent program for the average number of signal-setting single-photon counts in the X basis.
Proceeding analogously for the average number of signal-setting single-photon error counts in the X basis we have thatĒ where for each round k, E where E (k) 1,µ,h = E (k) 1,µ I {a k−1 =h} . We obtain therefore that: Then an upper boundĒ U 1,µ,N on the quantity Ē 1,µ,N is achieved by using the following linear program with the parameters t ± abc,n and s ± abc,n that arise from the linear CS constraints given in Appendix A.

D. Simulations
As shown in [38], the asympotitc secret key rate is well approximated by for large enough N, as long as the variances of the experimental averages vanish asymptotically (see [38] for a precise meaning of this statement). Here, h(x) denotes the binary entropy function, f EC is the error correction efficiency, the parameterZ µ,N is the gain in the Z basis defined asZ and E tol is the overall error rate observed in the Z basis.
Note that here one cannot simply use the asymptotic secret key rate formula against collective attacks due to the presence of intensity correlations. To be precise, we have that the asymptotic equivalence between the collective and the coherent settings is typically established on the basis of the so-called post-selection technique [45] built on the De Finetti theorem [46,47]. However, the round-exchangeability property required to apply the post-selection technique is generally invalidated by pulse-correlations.
Even though simulations with real data from [34] are presented in Sec. VI, in order to compare our fine-grained analysis with the previous results in [38], we fix the experimental inputs of the linear programsZ a,c,N /q 2 Z p a p c , X a,c,N /q 2 X p a p c andĒ a,c,N /q 2 X p a p c to their expected values according to a typical channel model, which is provided in Appendix C.
For that matter, let η det denote the common detection efficiency of Bob's detectors, and let η ch = 10 −αattL/10 be the transmittance of the quantum channel, where α att (dB/km) is its attenuation coefficient and L (km) is the distance. Also, let p d denote the dark count probability of each of Bob's detectors and let δ A stand for the po-larization misalignment occurring in the channel. This model yields the following results (see Appendix C) and Ē a,c,N for a, c ∈ A and where η = η det η ch . Here, we define the parameter h η,a,c,δA as h η,a,c,δA = e −ηa cos 2 δA − e −ηa sin 2 δA 2 .
We also introduce the parameterĒ a,c,N (Z) , which is equivalent toĒ a,c,N but accounts for the Z basis error clicks. The tolerated bit error rate of the sifted key is set to E tol = Ē µ,N (Z) / Z µ,N . Also, the reference parameters for the linearized CS constraints are fixed by the channel model as well, as indicated in Appendix C.

FIG. 2.
Secret key rate given by Eq.(33) using the analysis provided by Eq. (11). We consider different values for the maximum deviation δmax ∈ {10 −2 , 10 −3 , 10 −4 } between the actual physical intensity α k and the selected intensity setting a k , and two values of the parameter ξ. Precisely, we consider ξ = 1 (i.e. corresponding to the case of nearest-neighbour pulse correlations) and ξ = 2 (corresponding to pulse correlations of range two). The coarse-grained analysis follows the techniques presented in [38].
In particular, in the simulations we take η det = 0.65 and p d = 7.2 × 10 −8 [48]. The attenuation coefficient of the channel is set to α att = 0.2 dB/km, the error correction efficiency to f EC = 1.16 and a channel misalignment of δ A = 0.08 is used. As for the intensities, the weakest intensity setting is set to ω = 10 −4 due to the finite extinction ratio of intensity modulators. For simplicity, due to the large number of constraints included in the linear programs, instead of optimizing both µ and ν to maximize the asymptotic secret key rate K ∞ as a function of the distance L, we select a pair (µ, ν) that roughly maximizes the achievable distance and use that pair for all values of L. Notably, since we are considering the asymptotic regime, K ∞ does not depend on the probabilities of the decoy settings, p ν and p ω , nor on the probability of selecting the X basis, q X , in such a way that setting p µ ≈ 1 and q Z ≈ 1 maximizes K ∞ .
In Fig. 2 we illustrate the effect of the intensity correlations in a rate-distance representation for various values of the maximum relative deviation δ max ∈ {10 −2 , 10 −3 , 10 −4 }, and for ξ = 1 (i.e., for the case of nearest-neighbour pulse correlations) and ξ = 2 (corresponding to pulse correlations of range two). The simulations show that the secret key rate is very sensitive to the deviation δ max . What is more, with realistic experimental data [34], we find that the maximum distance attainable is less than 50 km. Moreover, it is important to remark that, as expected, the fine-grained analysis presented here outperforms the coarse-grained analysis introduced in [38] regardless the value of δ max . Indeed, Fig. 2 shows that when δ max = 10 −2 now the attainable distance is approximately double than that in [38]. Also, note that this figure assumes a single worst-case deviation parameter δ max (following Eq. (11)). The case corresponding to Eq. (10) is illustrated in the next section, when we compare this model with the one where Alice and Bob know the probability density function of the correlations.

A. Characterization
In this section we show that a significant improvement on the secret key rate presented in Fig. 2 can be obtained when Alice and Bob know the probability density function given by Eq. (2). By assuming this, the analysis that we present now is similar to that of the previous section but simpler, as we do not have to bound the photonnumber statistics, but we can solve them numerically.
The main motivation to consider this scenario is that recent decoy-state QKD experiments [34] performed in the high-speed regime seem to indicate that the correlation function is not arbitrary but gaussian-shaped (see also [44]). Triggered by this observation, for concreteness and illustration purposes we shall we assume here a truncated gaussian (TG) distribution for the correlations, which follows from renormalizing a gaussian distribution to a fixed finite interval [λ, Λ]. We note, however, that the analysis presented in this section could be straightforwardly adapted to any other probability density function of the correlations. To be precise, given both the mean and the variance -say, γ and σ 2 respectively-of the parental gaussian distribution together with the truncation interval [λ, Λ], the TG model reads From now on, we shall denote byγ andσ 2 the mean and variance of the TG distribution. We also call γ a k , the average intensity given the sequence a k , and σ a k , the standard deviation given the sequence a k of the parental normal distribution, although below the subscript indicating the record of settings will be omitted unless it can lead to confusion. More details on the truncation model are given in Appendix D.
Importantly, under suitable truncation conditions, this model reproduces the observed gaussian-shaped correlations in a finite support [λ, Λ] (in contrast to the unbounded support of non-truncated gaussian distributions). We note as well that the relevant parameters can be estimated experimentally by monitoring the intensities in long sequences of rounds.
From Eq. (38), the photon-number statistics now read (40) In general, as has been experimentally observed in [33], γ a k = a k for a k ∈ {µ, ν, ω}, i.e. a certain displacement in the mean intensity typically occurs. We account for this shift in the signal and the decoy intensity settings, while we neglect it for the vacuum intensity setting following the observations of [34]. Notably as well, σ 2 a k generally depends on the previous intensity settings, and so far no assumption about the range of the correlations has been made. In this sense, despite the fact that the analysis is valid for an arbitrary range, for simplicity here we only consider nearest-neighbour correlations for illustration purposes, i.e., below we assume ξ = 1.
Following the calculations in Appendix B 2, and combining the analysis above with the TG model here presented, the bound on the squared overlap in the nearestneighbour setting simply reads where we have introduced a threshold photon-number n cut for the numerics.

B. Linear programs for parameter estimation
The techniques for calculating the relevant parameters to evaluate the secret key rate formula using the decoystate constraints in this model are common with those deployed in the model-independent case evaluated in the previous section. The only difference is that certain steps such as the one that leads to Eq. (14) are no longer necessary. Here we solve the photon-number statistics given by Eq. (40) numerically and there is no need to invoke monotonicity arguments to bound them. The resulting linear program for the relevant single-photon yield (error click probability in the X basis) reads:    h n,a,c,N ≥ h n,b,c,N (a, b, c ∈ A, b = a, n = 0, . . . , n cut ) , t − abc,n + s − abc,n h n,a,c,N ≤ h n,b,c,N (a, b, c ∈ A, b = a, n = 0, . . . , n cut ) , 0 ≤ h n,a,b,N ≤ 1 (a, b ∈ A, n = 0, . . . , n cut ) . (43)

C. Simulations
The TG model requires to experimentally determine all of its parameters, namely, the truncation ranges, the mean intensities and the standard deviations of the distribution for every combination of pulses in rounds k and k − 1 (if one assumes, as already mentioned, nearestneighbour intensity correlations), which in turn should be experimentally accessible by monitoring the output of the transmitter in long sequences of rounds.
As discussed in Appendix D, regardless of the value of the mean intensity and the standard deviation, that in principle depend on the present and the previous settings, for illustration purposes we determine the truncation range for these simulations as (λ, Λ) = (γ − tσ,γ + tσ), whereγ is the measured mean intensity of the TG distribution andσ is the measured standard deviation of the TG distribution. Note that ideally, the truncation range could be measured experimentally. We can now make a direct comparison between this model in which we assume that Alice and Bob know the probability density function of the correlations, and the one that uses a maximum relative deviation defined in Eq. (10).
For this, we use Eq. (33) and the channel model presented in Appendix C to evaluate the performance. Moreover, regarding the mean and the standard deviations as a function of the intensity settings in rounds k and k − 1, we take the experimental values reported in [34]. The only exception is the normalized standard deviationσ =σ/γ corresponding to vacuum, where [34] does not report any value and we select for illustration purposes 10 −5 for all possible intensities in the round k − 1. This is so because the fluctuation of the vacuum setting seems to be essentially negligible. In addition, for simplicity, we do not optimize the intensity settings, but fix them to µ = 0.5, ν = 0.2 and ω = 10 −4 , which are the values considered in [34], and we select the parameter t = 4. The relevant parameters are summarized in Table I.
To facilitate a fair comparison between this model and   [34]. The only exception is the normalized standard deviationσ corresponding to vacuum, where [34] does not report any value and we select for illustration purposes 10 −5 . In our simulations we consider that they characterize a TG distribution.
that of Eq. (10), we take as explained in Appendix D, where i, j ∈ A. Here,γ i,j (σ i,j ) refers to the average mean intensity (standard deviation) of round k associated to the record (a k−1 , a k ) = (j, i), whileσ i,j =σ i,j /γ i,j are the normalized standard deviations of the TG model. In this way, we ensure that the physical intensity is bounded in the same interval in both cases. These parameters are obtained from Table I. Fig. 3 illustrates the improvement in the secret key rate when the correlation function is known and corresponds to a TG distribution according to the model that was just introduced. We find that now the results are comparable to those of the ideal case without intensity correlations, which highlights the importance of characterizing the probability density function of the correlations. When comparing the two models studied in this paper, we conclude that knowing the distribution g a k (α k ) allows for a significant improvement of the resulting performance. This is due to the fact that in this scenario one can evaluate the photon-number statistics exactly, which makes the parameter estimation much tighter, improving both the decoy-state constraints and the CS constraints. Finally, Fig. 3 also includes a representation where the standard deviations are twice the values given in Table I, to evaluate the effect that this parameter has on the secret key rate. The fine-grained analysis described in Eq. (10) does not provide key in this latter scenario. Likewise, if we use Eq. (11) in the model-independent case, no key is obtained in neither case.

FIG. 3.
Secret key rate within the TG model assuming the parameters from Table I (solid blue line), and when the normalized standard deviations are twice those illustrated in that table (dash-dotted green line). For comparison, this figure also includes the model-independent case (dashed blue line) where the maximum relative deviation is taken as δ (a k =i,a k−1 =j) = tσi,j for every i, j ∈ A, whereσi,j stands for the normalized standard deviations in the TG model, and we select t = 4. When we consider twice the normalized standard deviations provided in Table I, this latter model provides no key.

VII. CONCLUSIONS
When combining high-speed clock rate QKD transmitters with the decoy-state technique, the presence of intensity correlations in the generated pulses invalidate the central security assumption of this method, namely, that the yields and error rates associated to n-photon pulses are independent of the intensity settings. This problem can be solved by imposing constraints on the intensity-setting dependence of these parameters, which is done here by invoking the so-called Cauchy-Schwarz constraints. Essentially, these constraints arise from the indistinguishability of non-orthogonal quantum states and can be derived on the basis of a minor characterization of the correlations.
In this work, by using a standard decoy-state BB84 protocol with three possible intensity settings, we have evaluated the effect that such correlations have in the secret key rate. We have done so by introducing a finegrained analysis able to handle arbitrary long-range correlations. This approach leads to significantly tighter constraints for the parameter estimation and notably improves the secret key rate achieved by previous works. Also, we have shown that the characterization of the probability density function of the correlations permits to notably improve the resulting performance, which now becomes comparable to that of the ideal scenario without intensity correlations. Putting it all together, the present work provides a solid step forward towards full implementation security in QKD with high performance.

ACKNOWLEDGMENTS
This work was supported by Cisco Systems Inc., the Galician Regional Government (consolidation of Research Units: AtlantTIC), the Spanish Ministry of Economy and Competitiveness (MINECO), the Fondo Europeo de Desarrollo Regional (FEDER) through Grant No. PID2020-118178RB-C21, and MICIN with funding from the European Union NextGenerationEU (PRTR-C17.I1) and the Galician Regional Government with own funding through the "Planes Complementarios de I+D+I con las Comunidades Autónomas" in Quantum Communication.

Appendix A: Cauchy-Schwarz constraint
The CS constraint can be stated as follows: Theorem 1 Let |u and |v be pure states of a certain quantum system. Then, for all positive operatorsÔ ≤ I, where the functions G ± are given in Eq. (6).
This constraint can be used to impose a restriction on the maximum bias that Eve can induce between the nphoton yields and error probabilities associated to different intensity settings. Here we derive a linear version of Eq. (4), enabling the use of linear programming for the decoy-state parameter estimation procedure. Following [38], we have that in virtue of the convexity/concavity of the functions that define the constraints, their first-order expansions around any given reference value provide valid linear bounds as well. Thus, if we focus on the yields, the linearization provides the following bounds: Similarly for the error probabilities, and assuming that we select reference parameters independent of Alice's bit value r k (for a more detailed analysis see [38]) we have: To finish with, note that one can restrict the reference parameters to be round-independent, in such a way that, summing over k and dividing by N in Eqs. (7) and (8), we obtain the average parameters and similarly for the other terms that appear below.
We define the intercepts and slopes as: The linear CS constraints for the round independent case are: and Note that the characterization of the quantum channel is essential, as the tightness of these linear bounds is sub-ject to the adequacy of the selected reference parameters, and so it relies on it. The selected reference values are presented in Appendix C.
Appendix B: Derivation of τ ξ v 0 ,w 0 ,...v ξ ,n In an entanglement based view of the protocol, the global input state describing all the protocol rounds reads: where a N 1 = a 1 . . . a N and equivalently for x N 1 and r N 1 . Also, for every round i, The state |0 E stands for the initial state of Eve's ancillary system. We also define: where C i denotes an inaccessible purifying system with orthonormal basis {|t ni Ci |n i ∈ N} (C i stores the photon-number information of the i-th signal that Alice sends to Bob), B i denotes the system delivered to Bob, |n xi,ri i Bi stands for a Fock state with n i photons encoding the BB84 polarization state defined by (x i , r i ), and the photon-number statistics p ni | ai are defined in Eq. (2).
Let us denote byÛ BE Eve's coherent interaction with systems B 1 , ..., B N and E so thatÛ BE |Ψ represents the global state prior to Bob's measurements. We also refer to Bob's click POVM element in round k aŝ The joint probability p (k) (click, n, v 0 , ..., v ξ , Z) is computed as: A k = {A j |j = k, .., k − ξ},Ã k = {Ã j |j = k} and C k = {C j |j = k}. Note that, for round k, we project on the intensity setting, the basis and the photon-number, and for all previous rounds, we only project on the intensity setting. The unnormalized pure state is defined as: We now re-state Eq. Therefore, in virtue of Bayes rule we obtain: where we have defined Now, from Eq. (3), we recall that the quantity defined above is actually the n-photon yield of round k associated to the record of settings v 0 , ..., v ξ . From the CS constraint, it follows that: for all n ∈ N, w 0 , v 0 , ..., v ξ ∈ A where v 0 = w 0 and k = 1, ..., N . Note that the closer the inner product between the states is to one, the tighter the bounds of the inequality are.
Having reached this stage, we recall that the goal is to set a lower bound on the inner product Ψ where x k = {x j |j = k, } and a k = {a j |j = k, ..., k − ξ}. Then, computing the inner product yields: where we have omitted the subscript a k−ξ−1 in the square root term in parentheses because the photonnumber statistics are independent of this sub-string of settings. Since all the factors previous to round k are equal to 1, and the sums over x k and r N 1 yield In terms of the normalized states, we obtain: That is, this quantity takes the same value given by Eq. (B16) for all n.

Model-independent correlations
Eq. (B16) provides a general formula that can be used as the input of the CS constraint regardless of the correlation function. Here we consider the model-independent scenario where the correlation function is assumed to be unknown. This implies that we cannot evaluate the photon-number statistics on which Eq. (B16) depends. Therefore, we bound these statistics by invoking monotonicity arguments.
Only to keep the notation simple, we consider that the mean physical intensity matches the actual intensity setting, even though according to experiments, there might to be a certain displacement between these two quanti-ties [34]. We emphasize, however, that such shift could be straightforwardly included in our analysis. Let us start by introducing the following shorthand notation: a (τ )± i = a i 1 ± δ ± (ai,...,a k+1 ,a k =τ,a k−1 =v1,...,a i−ξ =v ξ+k−i ) , (B17) with τ ∈ {v 0 , w 0 }.
From the definition of the photon-number statistics and by noticing that e −x x n is strictly decreasing for n = 0 and increasing for n ≥ 1 in x ∈ (0, 1), we have that By combining Eqs. (B16)-(B19), we obtain the bound given by Eq. (10). A simpler bound given by Eq. (11) arises by considering that the relative deviation parameters are all equal. It can be obtained directly from Eq. (10).
It is also worth mentioning that one could use a more exhaustive approach by considering different records of settings for each of the two yields (or error probabilities) to be compared via the CS constraint. However, our numerical simulations suggest that the improvement achieved by doing so is not really significant, while the analysis is much more cumbersome to implement numerically.
If we focus our attention on nearest-neighbour intensity correlations only, and call v 0 = a, w 0 = b and v 1 = c, then from Eqs. (10) and (11)  Following Eq. (B16) and using the same notation as in the previous section for nearest-neighbour intensity correlations, we have that the bound for the different possible combinations of settings when the correlation function is a TG distribution is given by p n k+1 | a k+1 ,a p n k+1 | a k+1 ,b   2 ≡ τ ξ=1 a,b,c,n .
To gain intuition about the form of this distribution and the effect that the truncation range has onσ, in Fig. 5 we plot it for a couple of different intervals [λ, Λ] centered in the mean value. Note that, in this way, the mean is not affected by the truncation.
To compare this model with the model-independent case, we need to relate the truncation interval and δ max , which we recall that fully characterizes the model independent setting. For this matter, we take (λ, Λ) = (γ − tσ,γ + tσ). We state the relation between δ max and the parameter t asγ(1 − δ max ) =γ − tσ. In other words, we demand that the maximum relative deviation matches the extreme of the truncation. Note that, ideally, the interval (λ, Λ) could be measured experimentally.
Regarding the value of t, for the sake of comparison with the model-independent case in the main text, we use t = 4. It is also worth mentioning that modifying the value of t has a negligible effect on the secret key rate K ∞ , as shown in Fig. 6.