Quantum key distribution overcoming practical correlated intensity fluctuations

Abstract

Intensity correlations between neighbouring pulses open a prevalent yet often overlooked security loophole in decoy-state quantum key distribution (QKD). As a solution, we present and experimentally demonstrate an intensity-correlation-tolerant QKD protocol that mitigates the negative effect that this phenomenon has on the secret key rate according to existing security analyses. Compared to previous approaches, our method significantly enhances the robustness against correlations, notably improving both the maximum transmission distances and the achievable secret key rates across different scenarios. By relaxing constraints on correlation parameters, our protocol enables practical devices to counter intensity correlations. We experimentally demonstrate this first practical solution that directly overcomes this security vulnerability, establish the feasibility and efficacy of our proposal, taking a major step towards loophole-free and high-performance QKD.

Introduction

Quantum key distribution¹ (QKD) stands at the forefront of secure communication protocols², as it enables two remote users, Alice and Bob, to share secret keys whose security is guaranteed by the principles of quantum mechanics^{3, 4, 5–6}. However, QKD security proofs typically extract secret bits from the raw data originating from single-photons, and on-demand high-quality single-photon sources at telecom wavelengths are not available yet. A popular solution to solve this pressing issue is the decoy-state method^{7, 8–9}, which provides the same secret key rate scaling as single-photon sources by means of using laser sources emitting phase-randomized weak coherent pulses (PRWCPs). Indeed, this technique is a standard tool in current QKD implementations^{10, 11, 12–13}.

An important breakthrough direction for advancing QKD is to increase its secret key rate, for which fast operating QKD systems are being developed^{10, 11, 12, 13, 14, 15, 16–17}. However, the implementation security problem¹⁸ brought by high-clock-rate decoy-state QKD systems cannot be ignored. Due to memory effects in the devices and the electronics that control them, high-clock-rate decoy-state QKD systems face a troublesome implementation security problem: intensity correlations^{16,19, 20, 21, 22, 23, 24–25}. This means that the intensity setting of any given round may influence the actual intensity emitted in subsequent rounds, resulting in a partial distinguishability of the intensity settings. This breaks a core assumption of the decoy-state method, posing an underestimated threat to the security of QKD^{24, 25–26}. To address this problem, various security analyses have been proposed^21,24,25. However, practical devices^{27, 28–29} often struggle to meet the criteria set by these analyses (such as the magnitude of the correlations) leading to low or even vanishing secret key rates.

Here, we solve this crucial limitation by proposing an approach that we name intensity-correlation-tolerant QKD, which is capable of mitigating the intensity correlations problem in QKD devices. By adding a local monitor, our protocol enables common devices to achieve notably higher secret key rates and longer transmission distances than previous solutions in the presence of this type of correlations. Importantly, we experimentally demonstrate the feasibility and effectiveness of our approach. This advancement is a significant step towards loophole-free and high-performance QKD.

Results

Assumptions

To characterize intensity correlations and fluctuations, we assume a general model from ref. ²⁴. Precisely, in any given round k, the correlations do not compromise the Poissonian character of the photon-number statistics of the source conditioned on the value of the actual intensity, α_k. Nevertheless, the latter does not match the selected intensity setting, say a_k, but it is also influenced by the previous ξ settings, a_k−ξa_k−ξ+1…a_k−1, for a certain correlation range ξ. That is, we model α_k as a random variable, in such a way that every possible record {k, ξ}: = a_k−ξa_k−ξ+1,…,a_k determines a conditional probability distribution for α_k, with expectation ${\overline{\alpha }}_{\left\{k,\xi \right\}}$. Here, we describe this randomness using the relative deviation δ_{{k, ξ}} —such that ${\alpha }_k={\overline{\alpha }}_{\left\{k,\xi \right\}}\left(1+{\delta }_{\left\{k,\xi \right\}}\right)$— and assume that the probability density function of δ_{{k, ξ}}, g_{{k, ξ}}, is only nonzero in a certain record-dependent interval $\left[{\delta }_{\left\{k,\xi \right\}}^{-},{\delta }_{\left\{k,\xi \right\}}^{+}\right]$. For later convenience, we define ${\delta }_{\max ,a_k}^{\mathrm{rand}}={\max }_{\left\{k,\xi \right\}}\left\{\mid {\delta }_{\left\{k,\xi \right\}}^{\pm }\mid \right\}$ and ${\delta }_{\max ,a_k}^{\mathrm{corr}}={\max }_{\left\{k,\xi \right\}}\left\{\mid 1-{\overline{\alpha }}_{\left\{k,\xi \right\}}/a_k\mid \right\}$, which denote the magnitude of sequence-dependent random fluctuations and the size of correlation of the average intensity, respectively.

We make the following assumptions on the intensity correlations and fluctuations.

Assumption 1: The intensity correlations do not compromise the Poissonian character of the photon-number statistics of the source^16,21,22,27. That is to say, given the actual intensity α_k prepared in the k-th round, the conditional photon-number statistics satisfy

1

$p_n{\mid }_{{\alpha }_k}=\frac{{\mathrm{e}}^{-{\alpha }_k}{\left({\alpha }_k\right)}^n}{n!}.$

Assumption 2: The intensity correlations have a finite range ξ, meaning that the intensity setting a_k−i of the (k − i)-th round does not influence α_k if i > ξ.

Assumption 3: Let us introduce the shorthand notation {i, j} to describe the record of settings a_i−ja_i−j+1. . . a_i. Each record {k, ξ} determines a conditional probability distribution for α_k, with expectation ${\overline{\alpha }}_{\left\{k,\xi \right\}}$. For convenience, we describe this randomness using the relative deviation δ_{{k, ξ}}, such that ${\alpha }_k={\overline{\alpha }}_{\left\{k,\xi \right\}}\left(1+{\delta }_{\left\{k,\xi \right\}}\right)$. The probability density function of δ_{{k, ξ}}, which we denote by g_{{k, ξ}}, is only nonzero in the record-dependent interval $\left[{\delta }_{\left\{k,\xi \right\}}^{-},{\delta }_{\left\{k,\xi \right\}}^{+}\right]$.

We remark that, by definition,

2

${\int }_{{\delta }_{\left\{k,\xi \right\}}^{-}}^{{\delta }_{\left\{k,\xi \right\}}^{+}}{g}_{\left\{k,\xi \right\}}\left(\delta \right)\delta \mathrm{d}\delta =0.$Also, for a given record {k, ξ}, the conditional photon-number statistics of round k satisfy

3

$p_n{\mid }_{\left\{k,\xi \right\}}={\int }_{{\delta }_{\left\{k,\xi \right\}}^{-}}^{{\delta }_{\left\{k,\xi \right\}}^{+}}{g}_{\left\{k,\xi \right\}}\left(\delta \right)\exp \left(-{\overline{\alpha }}_{\left\{k,\xi \right\}}\left(1+\delta \right)\right)\frac{{\left[{\overline{\alpha }}_{\left\{k,\xi \right\}}\left(1+\delta \right)\right]}^n}{n!}\mathrm{d}\delta .$Finally, for simplicity, we introduce ${\delta }_{\max ,a_k}^{\mathrm{rand}}={\max }_{\left\{k,\xi \right\}}\left\{\mid {\delta }_{\left\{k,\xi \right\}}^{\pm }\mid \right\}$ and ${\delta }_{\max ,a_k}^{\mathrm{corr}}={\max }_{\left\{k,\xi \right\}}\left\{\mid 1-{\overline{\alpha }}_{\left\{k,\xi \right\}}/a_k\mid \right\}$.

Enhanced decoy-state method

Importantly, in the presence of intensity correlations, the n-photon yield and error rate associated to different intensity settings can be distinct. To address this issue, we introduce an enhanced decoy-state method. Crucially, we consider finer-grained decoy-state constraints by grouping the measurement counts at Bob’s receiver according to the record of settings rather than using the last setting alone as in standard analyses. Specifically, our protocol involves imposing constraints in two key aspects. The first is a photon number constraint. Similar to the standard decoy-state method, we truncate the photon number at n_cut to establish upper and lower bounds. However, due to security concerns arising from correlations, we need to classify the detection rates and error rates according to the sequence {k, ξ}. Despite this, the constraints alone are insufficient for parameter estimation, which leads us to introduce a second constraint. Specifically, we employ a mathematical tool called the Cauchy-Schwarz (CS) constraint^24,25 to set limits on the bias in detection statistics associated with different intensities.

As originally observed in ref. ²⁴, in the presence of intensity correlations, the yields and error probabilities in any given round may differ for different records of settings. As a consequence, the decoy-state method^{7, 8–9} alone does not enable a tight parameter estimation, and additional constraints are required. In accordance with ref. ²⁴, we address this problem by using the so-called Cauchy-Schwarz (CS) constraints, a tool previously exploited in refs. ^25,26 too. Specifically, the CS constraint restricts the bias between the measurement statistics of two different quantum states when subject to the same measurement. The CS constraint can be stated as follows.

Theorem 1

Let $\mid a\left(\right.⟩$ and $\mid b\left(\right.⟩$ be two pure states of an arbitrary Hilbert space $\mathcal{H}$. For any operator $\widehat{O}$ on $\mathcal{H}$ such that $0\le \widehat{O}\le 1$,

4

$\begin{array}{c}\hfill G_{-}\left(\mathrm{Tr}\left[\widehat{O}\mid a\left(\right.⟩\left.⟨\right)a\mid \right],{∣⟨a\mid b⟩∣}^2\right)\le \mathrm{Tr}\left[\widehat{O}\mid b\left(\right.⟩\left.⟨\right)b\mid \right]\\ \hfill \le G_{+}\left(\mathrm{Tr}\left[\widehat{O}\mid a\left(\right.⟩\left.⟨\right)a\mid \right],{∣⟨a\mid b⟩∣}^2\right)\end{array},$where

5

$\begin{array}{c}{G}_{-}=\left\{\begin{array}{cc}{g}_{-}\left(x,y\right)\hfill & x>1-y\hfill \\ 0\hfill & x\le 1-y\hfill \end{array}\right),\hfill \\ G_{+}=\left\{\begin{array}{cc}1\hfill & x\ge y\hfill \\ g_{+}\left(x,y\right)\hfill & x<y\hfill \end{array}\right),\hfill \end{array}$with $g_{\pm }\left(x,y\right)=x+\left(1-2x\right)\left(1-y\right)\pm 2\sqrt{x\left(1-x\right)y\left(1-y\right)}$.

Proof: See the Supplementary Materials of ref. ²⁶.

Essentially, Eq. (4) allows to set quantitative bounds on the detection statistics arising from different records of settings. However, because of their non-linearity, these bounds cannot be directly plugged into decoy-state linear programmes. Notwithstanding, in virtue of the convexity of the G_±(x, y) functions, suitable linearizations of the CS constraints follow. In particular, for any reference value c ∈ [0, 1], we have

6

$\begin{array}{c}{G}_{-}\left(c,{∣⟨a\mid b⟩∣}^2\right)+G_{-}^{\prime }\left(c,{∣⟨a\mid b⟩∣}^2\right)\left(\mathrm{Tr}\left[\widehat{O}\mid a\left(\right.⟩\left.⟨\right)a\mid \right]-c\right)\hfill \\ \le \mathrm{Tr}\left[\widehat{O}\mid b\left(\right.⟩\left.⟨\right)b\mid \right]\hfill \\ \le G_{+}\left(c,{∣⟨a\mid b⟩∣}^2\right)+G_{+}^{\prime }\left(c,{∣⟨a\mid b⟩∣}^2\right)\left(\mathrm{Tr}\left[\widehat{O}\mid a\left(\right.⟩\left.⟨\right)a\mid \right]-c\right)\hfill \end{array},$where

7

$\begin{array}{c}{G}_{-}^{\prime }=\left\{\begin{array}{cc}{g}_{-}^{\prime }\left(x,y\right)\hfill & x>1-y\hfill \\ 0\hfill & x\le 1-y\hfill \end{array}\right),\hfill \\ G_{+}^{\prime }=\left\{\begin{array}{cc}0\hfill & x\ge y\hfill \\ g_{+}^{\prime }\left(x,y\right)\hfill & x<y\hfill \end{array}\right),\hfill \end{array}$with $g_{\pm }^{\prime }\left(x,y\right)=-1+2y\pm \left(1-2x\right)\sqrt{y\left(1-y\right)/x\left(1-x\right)}$. Obviously, although the value of c can be set arbitrarily, it has a direct impact on the tightness of the constraints. In this regard, our preferred choices for these reference parameters are discussed in section ‘Reference values for the linear CS constraints’.

One of the key ideas for addressing intensity correlations is that, in the context of general attacks, setting CS constraints on the measurement statistics of a fixed round k requires to compute the inner product between specific quantum states across all N protocol rounds^24,25 (although respectively conditioned on the two records of settings whose statistics are to be compared). Define this quantity as the intensity correlation parameter ${\tau }_{a_k,a_k^{\prime },n}^{\xi }$, which satisfies

8

$\begin{array}{c}\hfill {\tau }_{a_k,a_k^{\prime },n}^{\xi }:=\left\{\sum _{a_{k+1}\in A}\dots \sum _{a_{k+\xi }\in A}\left(\prod _{i=k+1}^{k+\xi }{p}_{a_i}\sum _{m=0}^{\infty }\right)\right)\\ \hfill {\left(\left(\sqrt{p_m^{\mathrm{L}}{\mid }_{a_{i-\xi }\dots a_k\dots a_i}{p}_m^{\mathrm{L}}{\mid }_{a_{i-\xi }\dots a_k^{\prime }\dots a_i}}\right)\right\}}^2\end{array},$for any given round k, any pair of distinct settings a_k and $a_{k^{\prime }}$, and any photon number n, where $p_n^{\mathrm{L}}{\mid }_{\left\{k,\xi \right\}}$ is the lower bond when we estimate p_n∣_{{k, ξ}}. The detail proof of the intensity correlation parameter ${\tau }_{a_k,a_k^{\prime },n}^{\xi }$ can be seen in section ‘Intensity correlation parameter’. The intensity correlation parameters are calculated in section ‘Calculation of correlation parameters’.

In the first place, the standard decoy-state constraints can be stated in the form

9

$\frac{Z_{\left\{k,\xi \right\}}}{q_Z^2{\prod }_{i=k-\xi }^k{p}_{a_i}}=\sum _{n=0}^{\infty }{p}_n{\mid }_{\left\{k,\xi \right\}}{y}_{n,\left\{k,\xi \right\}},$where the left-hand side is the probability of a click occurring in round k conditioned on a Z-basis match and the record of settings being {k, ξ}, and y_{n,{k, ξ}} denotes the corresponding n-photon yield. As is customary in decoy-state analyses, one can select a threshold photon-number n_cut in order to split Eq. (9) into two complementary bounds,

10

$\begin{array}{c}\frac{Z_{\left\{k,\xi \right\}}}{q_Z^2{\prod }_{i=k-\xi }^k{p}_{a_i}}\ge \sum _{n=0}^{n_{\mathrm{cut}}}{p}_n^{\mathrm{L}}{\mid }_{\left\{k,\xi \right\}}{y}_{n,\left\{k,\xi \right\}}\left(a_{k-\xi },\dots ,a_k\in A\right),\hfill \\ \frac{Z_{\left\{k,\xi \right\}}}{q_Z^2{\prod }_{i=k-\xi }^k{p}_{a_i}}\le \sum _{n=0}^{n_{\mathrm{cut}}}{p}_n^{\mathrm{U}}{\mid }_{\left\{k,\xi \right\}}{y}_{n,\left\{k,\xi \right\}}+1\hfill \\ -\sum _{n=0}^{n_{\mathrm{cut}}}{p}_n^{\mathrm{L}}{\mid }_{\left\{k,\xi \right\}}\left(a_{k-\xi },\dots ,a_k\in A\right).\hfill \end{array}$Analogously, for the error statistics (say, of the X basis) we have

11

$\begin{array}{c}\frac{E_{\left\{k,\xi \right\}}}{q_X^2{\prod }_{i=k-\xi }^k{p}_{a_i}}\ge \sum _{n=0}^{n_{\mathrm{cut}}}{p}_n^{\mathrm{L}}{\mid }_{\left\{k,\xi \right\}}{h}_{n,\left\{k,\xi \right\}}\left(a_{k-\xi },\dots ,a_k\in A\right),\hfill \\ \frac{E_{\left\{k,\xi \right\}}}{q_X^2{\prod }_{i=k-\xi }^k{p}_{a_i}}\le \sum _{n=0}^{n_{\mathrm{cut}}}{p}_n^{\mathrm{U}}{\mid }_{\left\{k,\xi \right\}}{h}_{n,\left\{k,\xi \right\}}+1\hfill \\ -\sum _{n=0}^{n_{\mathrm{cut}}}{p}_n^{\mathrm{L}}{\mid }_{\left\{k,\xi \right\}}\left(a_{k-\xi },\dots ,a_k\in A\right),\hfill \end{array}$where $E_{\left\{k,\xi \right\}}/\left(q_X^2{\prod }_{i=k-\xi }^k{p}_{a_i}\right)$ is the probability of a click and a bit error occurring in round k conditioned on a X-basis match and the record of settings being {k, ξ}, and h_{n,{k, ξ}} denotes the corresponding n-photon error yield.

On the other hand, the linearized CS constraints that arise from Eq. (6), when applied to the yields, can be written as^24,25

12

$\begin{array}{c}{c}_{\left\{k-1,\xi -1\right\},a_k,a_k^{\prime },n}^{+}+m_{\left\{k-1,\xi -1\right\},a_k,a_k^{\prime },n}^{+}{y}_{n,\left\{k-1,\xi -1\right\}{a}_k}\ge y_{n,\left\{k-1,\xi -1\right\}{a}_k^{\prime }},\hfill \\ c_{\left\{k-1,\xi -1\right\},a_k,a_k^{\prime },n}^{-}+m_{\left\{k-1,\xi -1\right\},a_k,a_k^{\prime },n}^{-}{y}_{n,\left\{k-1,\xi -1\right\}{a}_k}\ge y_{n,\left\{k-1,\xi -1\right\}{a}_k^{\prime }},\hfill \end{array}$when $a_{k-\xi },\dots ,a_k,a_k^{\prime }\in A;a_k\ne a_k;n=0,\dots ,n_{\mathrm{cut}}$, with

13

$\begin{array}{c}{c}_{\left\{k-1,\xi -1\right\},a_k,a_k^{\prime },n}^{\pm }=G_{\pm }\left({\tilde{y}}_{n,\left\{k-1,\xi -1\right\}{a}_k},{\tau }_{a_k,a_k^{\prime },n}^{\xi }\right)\hfill \\ -G_{\pm }^{\prime }\left({\tilde{y}}_{n,\left\{k-1,\xi -1\right\}{a}_k},{\tau }_{a_k,a_k^{\prime },n}^{\xi }\right){\tilde{y}}_{n,\left\{k-1,\xi -1\right\}{a}_k},\hfill \\ m_{\left\{k-1,\xi -1\right\},a_k,a_k^{\prime },n}^{\pm }=G_{\pm }^{\prime }\left({\tilde{y}}_{n,\left\{k-1,\xi -1\right\}{a}_k},{\tau }_{a_k,a_k^{\prime },n}^{\xi }\right),\hfill \end{array}$where the reference values $\left\{{\tilde{y}}_{n,\left\{k,\xi \right\}}\right\}$ of the linearization are provided in section ‘Reference values for the linear CS constraints’.

In a similar fashion, the resulting linearized CS constraints for the error statistics read

14

$\begin{array}{c}{t}_{\left\{k-1,\xi -1\right\},a_k,a_k^{\prime },n}^{+}+s_{\left\{k-1,\xi -1\right\},a_k,a_k^{\prime },n}^{+}{h}_{n,\left\{k-1,\xi -1\right\}{a}_k}\ge h_{n,\left\{k-1,\xi -1\right\}{a}_k^{\prime }},\hfill \\ t_{\left\{k-1,\xi -1\right\},a_k,a_k^{\prime },n}^{-}+s_{\left\{k-1,\xi -1\right\},a_k,a_k^{\prime },n}^{-}{h}_{n,\left\{k-1,\xi -1\right\}{a}_k}\ge h_{n,\left\{k-1,\xi -1\right\}{a}_k^{\prime }},\hfill \end{array}$when $a_{k-\xi },\dots ,a_k,a_k^{\prime }\in A;a_k\ne a_k;n=0,\dots ,n_{\mathrm{cut}}$, for

15

$\begin{array}{c}\hfill t_{\left\{k-1,\xi -1\right\},a_k,a_k^{\prime },n}^{\pm }=G_{\pm }\left({\tilde{h}}_{n,\left\{k-1,\xi -1\right\}{a}_k},{\tau }_{a_k,a_k^{\prime },n}^{\xi }\right)\\ \hfill -G_{\pm }^{\prime }\left({\tilde{h}}_{n,\left\{k-1,\xi -1\right\}{a}_k},{\tau }_{a_k,a_k^{\prime },n}^{\xi }\right){\tilde{h}}_{n,\left\{k-1,\xi -1\right\}{a}_k},\\ \hfill s_{\left\{k-1,\xi -1\right\},a_k,a_k^{\prime },n}^{\pm }=G_{\pm }^{\prime }\left({\tilde{h}}_{n,\left\{k-1,\xi -1\right\}{a}_k},{\tau }_{a_k,a_k^{\prime },n}^{\xi }\right),\end{array}$where the reference values $\left\{{\tilde{h}}_{n,\left\{k,\xi \right\}}\right\}$ are given in section ‘Reference values for the linear CS constraints’.

Putting it all together, Eq. (10) and Eq. (12) plus the boundary constraints 0 ≤ y_{n,{k, ξ}} ≤ 1 compose the round-dependent linear programme of the Z-basis detection statistics, whose objective function is given by

16

$\begin{array}{ccc}\hfill Z_{1,\mu }^{\left(k\right)\mathrm{L}}& \hfill =\hfill & q_{\mathrm{Z}}^2{p}_{\mu }\sum _{a_{k-\xi }\in A}\dots \sum _{a_{k-1}\in A}\left(\prod _{i=k-\xi }^{k-1}{p}_{a_i}\right)\hfill \\ \hfill & \hfill \hfill & \times p_1^{\mathrm{L}}{\mid }_{\left\{k-1,\xi -1\right\}\mu }{y}_{1,\left\{k-1,\xi -1\right\}\mu }.\hfill \end{array}$Namely, a lower bound on the total probability of a Z-basis match and a signal-setting single-photon click occurring in round k (averaged over all possible records {k − 1, ξ − 1}).

Importantly, an equivalent LP for the X-basis detection statistics is obtained by replacing X by Z where convenient. The objective function of such LP, defined analogously as $Z_{1,\mu }^{\left(k\right)\mathrm{L}}$ but for the X basis, is denoted as $X_{1,\mu }^{\left(k\right)\mathrm{L}}$.

Analogously, Eq. (11) and Eq. (14) plus the boundary constraints 0 ≤ h_{n,{k, ξ}} ≤ 1 compose the round-dependent linear programme of the X-basis error statistics, whose objective function is given by

17

$\begin{array}{ccc}\hfill E_{1,\mu }^{\left(k\right)\mathrm{L}}& \hfill =\hfill & q_{\mathrm{Z}}^2{p}_{\mu }\sum _{a_{k-\xi }\in A}\dots \sum _{a_{k-1}\in A}\left(\prod _{i=k-\xi }^{k-1}{p}_{a_i}\right)\hfill \\ \hfill & \hfill \hfill & \times p_1^{\mathrm{U}}{\mid }_{\left\{k-1,\xi -1\right\}\mu }{h}_{1,\left\{k-1,\xi -1\right\}\mu }.\hfill \end{array}$Namely, an upper bound on the total probability of an X-basis match and a signal-setting single-photon bit error occurring in round k (averaged over all possible records {k − 1, ξ − 1}).

To finish with, the relevant round-independent linear programmes follow by summing over all N protocol rounds and dividing by N in both the objective functions and the constraints.

Calculation of correlation parameters

On the practical side, there are two challenges in implementing the enhanced decoy-state method. First, the photon number constraint requires to accurately estimate tight upper and lower bounds on the photon number statistics. Second, the CS constraints require to determine lower bounds on the squared inner product between suitable quantum states that only differ on the last intensity setting, which we have defined as ${\tau }_{a_k,a_k^{\prime },n}^{\xi }$ in Eq. (8). As our solution, by including a local monitor as an integral component of the protocol, Alice can obtain a tight estimation of these parameters. Precisely, she measures a fraction of every PRWCP she sends to the quantum channel. By categorizing the measurement results according to the setting sequences {k, ξ}, Alice obtains the gain for each category. Using this information, together with the accurate relative efficiency η_m between the pulse entering the channel and the pulse entering the detection module in the local monitor, she can estimate the different photon-number probabilities for each category, further determining ${\tau }_{a_k,a_k^{\prime },n}^{\xi }$.

In this section, we elucidate how our protocol calculates these parameters by means of using a local monitor situated at Alice’s side^23,30. Remarkably, if compared with the existing analyses^24,25, our estimates result in a significant mitigation of the damaging effect that intensity correlations have on the secret key rate.

The working principle of the local monitor is depicted in Fig. 1, although we remark that alternative configurations may fulfil the same purpose as well. By means of a laser diode (LD) and an Encoder, Alice prepares PRWCPs which are then split into two components by a balanced BS. While the transmitted pulses enter the quantum channel after undergoing attenuation by a variable optical attenuator (VOA), the reflected pulses are directed to the local monitor, which consists of a VOA (VOA-2) and a single-photon detector (SPD-M). The relative attenuation of the local monitor is defined as ${\eta }_{\mathrm{m}}={\eta }_{{\mathrm{d}}_{\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{i}}}{\eta }_{\mathrm{V}\mathrm{O}\mathrm{A}-2}/{\eta }_{\mathrm{V}\mathrm{O}\mathrm{A}-1}$ where ${\eta }_{{\mathrm{d}}_{\mathrm{moni}}}$ is the detection efficiency of SPD-M and η_{V OA−i} is the attenuation of VOA − i. The reason of choosing SPD is mentioned in section ‘Measurement of correlation parameters’.

[See PDF for image]

Fig. 1

Schematic of Alice’s module.

Quan. Mod. quantum module, LD laser diode, Encoder quantum state and decoy-state intensity modulator, BS beam splitter, VOA variable optical attenuator, SPD single-photon detector.

Due to the expression of the correlation parameters, Eq. (8), it is imperative to tightly estimate the conditional photon-number statistics for each possible record of intensity settings. Combining a standard characterization of single-photon detectors with the general assumptions of section ‘Assumptions’, it follows that the detection probability of the local SPD in the kth round conditioned on a record of settings {k, ξ} is given by

18

$D_{\left\{k,\xi \right\}}=1-\left(1-p_d\right)\left(1-p_{ap}\right){\int }_{{\delta }_{\left\{k,\xi \right\}}^{-}}^{{\delta }_{\left\{k,\xi \right\}}^{+}}{g}_{\left\{k,\xi \right\}}\left(\delta \right)e^{-{\eta }_{\mathrm{m}}{\overline{\alpha }}_{\left\{k,\xi \right\}}\left(1+\delta \right)}\mathrm{d}\delta ,$where p_d and p_ap respectively denote the dark count rate and the afterpulse rate of the local SPD. Defining ${\mathcal{D}}_{\left\{k,\xi \right\}}=\left(1-D_{\left\{k,\xi \right\}}\right)/\left[\left(1-p_d\right)\left(1-p_{ap}\right)\right]$, it follows that

19

${\mathcal{D}}_{\left\{k,\xi \right\}}={\int }_{{\delta }_{\left\{k,\xi \right\}}^{-}}^{{\delta }_{\left\{k,\xi \right\}}^{+}}{g}_{\left\{k,\xi \right\}}\left(\delta \right)e^{-{\eta }_{\mathrm{m}}{\overline{\alpha }}_{\left\{k,\xi \right\}}\left(1+\delta \right)}\mathrm{d}\delta .$Then, employing the Taylor expansion of function e^x and combine with the fact that the first-order term of δ vanishes (see Eq. (2)), we can calculate the upper and lower bond of ${\overline{\alpha }}_{\left\{k,\xi \right\}}$ (the detailed derivation is described in section ‘Measurement of correlation parameters’), satisfies

20

$\begin{array}{c}{\overline{\alpha }}_{\left\{k,\xi \right\}}^{\mathrm{U}}=\frac{1-\sqrt{1-2\left(1-{\mathcal{D}}_{\left\{k,\xi \right\}}\right)\left(1+{\zeta }_{\left\{k,\xi \right\}}\right)}}{{\eta }_{\mathrm{m}}\left(1+{\zeta }_{\left\{k,\xi \right\}}\right)}\hfill \\ {\overline{\alpha }}_{\left\{k,\xi \right\}}^{\mathrm{L}}=\frac{1-{\mathcal{D}}_{\left\{k,\xi \right\}}}{{\eta }_{\mathrm{m}}}+\frac{{\left(1-{\mathcal{D}}_{\left\{k,\xi \right\}}\right)}^2}{2{\eta }_{\mathrm{m}}}-\frac{{{\eta }_{\mathrm{m}}}^2{\left({\overline{\alpha }}_{\left\{k,\xi \right\}}^{\mathrm{U}}\right)}^3}{6}\hfill \end{array},$where ${\zeta }_{\left\{k,\xi \right\}}=\max \left\{{\left({\delta }_{\left\{k,\xi \right\}}^{-}\right)}^2,{\left({\delta }_{\left\{k,\xi \right\}}^{+}\right)}^2\right\}$.

Then we can tightly bound the conditional photon-number statistics, p_n∣_{{k, ξ}}. To this end, we re-use the trick of Taylor-expanding around the conditional expectation ${\overline{\alpha }}_{\left\{k,\xi \right\}}$. Similarly, performing a Taylor expansion on the $\exp \left(-{\overline{\alpha }}_{\left\{k,\xi \right\}}\left(1+\delta \right)\right)$ and ${\left[{\overline{\alpha }}_{\left\{k,\xi \right\}}\left(1+\delta \right)\right]}^n$ terms in Eq. (3) and utilizing the fact that the first-order term of δ is zero (the detailed derivation is described in section ‘Measurement of correlation parameters’), we obtain that

21

$\begin{array}{c}\hfill \underset{x,y}{\min }{f}_n\left(x,y\right)\le p_n{\mid }_{\left\{k,\xi \right\}}\le \underset{x,y}{\max }{f}_n\left(x,y\right),\\ \hfill \mathrm{for}x\in \left[{\delta }_{\left\{k,\xi \right\}}^{-},{\delta }_{\left\{k,\xi \right\}}^{+}\right],y\in \left[{\overline{\alpha }}_{\left\{k,\xi \right\}}^{\mathrm{L}},{\overline{\alpha }}_{\left\{k,\xi \right\}}^{\mathrm{U}}\right],\end{array}$where

22

$f_n\left(x,y\right)=\frac{y^n{\mathrm{e}}^{-y}}{n!}\left[{\mathrm{e}}^{-xy}{\left(1+x\right)}^n-\left(n-y\right)x\right].$On the contrary, invoking the monotonicity of the Poissonian photon-number probabilities in both ${\overline{\alpha }}_{\left\{k,\xi \right\}}$ and δ, it readily follows that

23

$\begin{array}{c}{\mathrm{e}}^{-{\overline{\alpha }}_{\left\{k,\xi \right\}}^{\mathrm{L}}\left(1+{\delta }_{\left\{k,\xi \right\}}^{-}\right)}\frac{{\left[{\overline{\alpha }}_{\left\{k,\xi \right\}}^{\mathrm{L}}\left(1+{\delta }_{\left\{k,\xi \right\}}^{-}\right)\right]}^n}{n!}\le p_n{\mid }_{\left\{k,\xi \right\}}\hfill \\ \le {\mathrm{e}}^{-{\overline{\alpha }}_{\left\{k,\xi \right\}}^{\mathrm{U}}\left(1+{\delta }_{\left\{k,\xi \right\}}^{+}\right)}\frac{{\left[{\overline{\alpha }}_{\left\{k,\xi \right\}}^{\mathrm{U}}\left(1+{\delta }_{\left\{k,\xi \right\}}^{+}\right)\right]}^n}{n!},\hfill \end{array}$for n ≥ 1. This latter approach overlooks the fact that the first-order term in δ vanishes upon integration, resulting in looser bounds than those of Eq. (21).

Since the tighter bounds of Eq. (21) further rely on explicit numerical optimization, we use these bounds only up to a certain threshold photon number n_th, and for n > n_th we apply the monotonicity bounds of Eq. (23) directly. In either case, we refer to these lower and upper bounds as $p_n{\mid }_{\left\{k,\xi \right\}}^{\mathrm{L}}$ and $p_n{\mid }_{\left\{k,\xi \right\}}^{\mathrm{U}}$, respectively.

Let us now calculate the lower bounds on the correlation parameters —defined in Eq. (8)— that arise from the $p_n{\mid }_{\left\{k,\xi \right\}}^{\mathrm{L}}$ of the previous subsection. Splitting ${\sum }_{m=0}^{\infty }\sqrt{p_m^{\mathrm{L}}{\mid }_{a_{i-\xi }\dots a_k\dots a_i}{p}_m^{\mathrm{L}}{\mid }_{a_{i-\xi }\dots a_k^{\prime }\dots a_i}}$ into two sums to separately make use of Eqs. (21) and (23) one can readily show that

24

$\begin{array}{c}\sum _{m=0}^{\infty }\sqrt{p_m^{\mathrm{L}}{\mid }_{a_{i-\xi }\dots a_k\dots a_i}{p}_m^{\mathrm{L}}{\mid }_{a_{i-\xi }\dots a_k^{\prime }\dots a_i}}\hfill \\ \ge \sum _{m=0}^{n_{\mathrm{th}}}\sqrt{\underset{x_1,y_1}{\min }{f}_m\left(x_1,y_1\right)\underset{x_2,y_2}{\min }{f}_m\left(x_2,y_2\right)}\hfill \\ +\sum _{m=n_{\mathrm{th}+1}}^{\infty }\left[{\mathrm{e}}^{-a_k\left(1-{\delta }_{\max ,a_k}^{\mathrm{rand}}\right)\left(1-{\delta }_{\max ,a_k}^{\mathrm{corr}}\right)}\right)\hfill \\ \times \left(\frac{{\left[a_k\left(1-{\delta }_{\max ,a_k}^{\mathrm{rand}}\right)\left(1-{\delta }_{\max ,a_k}^{\mathrm{corr}}\right)\right]}^m}{m!}\right]\hfill \\ =1+\sum _{m=0}^{n_{\mathrm{th}}}\left[\sqrt{\underset{x_1,y_1}{\min }{f}_m\left(x_1,y_1\right)\underset{x_2,y_2}{\min }{f}_m\left(x_2,y_2\right)}-{\mathrm{e}}^{-a_k^{\mathrm{L}}}\frac{{\left(a_k^{\mathrm{L}}\right)}^m}{m!}\right],\hfill \end{array}$where $x{\left(y\right)}_i\in \left[x{\left(y\right)}_i^{\mathrm{L}},x{\left(y\right)}_i^{\mathrm{U}}\right]$, $x_1^{\mathrm{L}}={\delta }_{a_{i-\xi }\dots a_k\dots a_i}^{-}$, $x_1^{\mathrm{U}}={\delta }_{a_{i-\xi }\dots a_k\dots a_i}^{+}$, $x_2^{\mathrm{L}}={\delta }_{a_{i-\xi }\dots a_k^{\prime }\dots a_i}^{-}$, $x_2^{\mathrm{U}}={\delta }_{a_{i-\xi }\dots a_k^{\prime }\dots a_i}^{+}$, $y_1^{\mathrm{L}}={\overline{\alpha }}_{a_{i-\xi }\dots a_k\dots a_i}^{\mathrm{L}}$, $y_1^{\mathrm{U}}={\overline{\alpha }}_{a_{i-\xi }\dots a_k\dots a_i}^{\mathrm{U}}$, $y_2^{\mathrm{L}}={\overline{\alpha }}_{a_{i-\xi }\dots a_k^{\prime }\dots a_i}^{\mathrm{L}}$, $y_2^{\mathrm{U}}={\overline{\alpha }}_{a_{i-\xi }\dots a_k^{\prime }\dots a_i}^{\mathrm{U}}$ and $a_k^{\mathrm{L}}=a_k\left(1-{\delta }_{\max ,a_k}^{\mathrm{rand}}\right)\left(1-{\delta }_{\max ,a_k}^{\mathrm{corr}}\right)$. Plugging this into Eq. (8) yields

25

$\begin{array}{c}{\tau }_{a_k,a_k^{\prime },n}^{\xi }\ge \left\{\sum _{a_{k+1}\in A}\dots \sum _{a_{k+\xi }\in A}\left[\left(\prod _{i=k+1}^{k+\xi }{p}_{a_i}\right)\right)\right)\hfill \\ {\left(\left\{1+\left(\sum _{n=0}^{n_{\mathrm{th}}}\left[\sqrt{\underset{x_1,y_1}{\min }{f}_n\left(x_1,y_1\right)\underset{x_2,y_2}{\min }{f}_n\left(x_2,y_2\right)}-{\mathrm{e}}^{-a_k^{\mathrm{L}}}\frac{{\left(a_k^{\mathrm{L}}\right)}^n}{n!}\right]\right\}\right]\right\}}^2.\hfill \end{array}$

Asymptotic secret key rate and simulation results

In the simulations, we follow the asymptotic analysis. For large enough N, one can identify the round-averaged observables with their expectations to approximate the secret key rate via

26

$R\approx Z_{1,\mu }^{\mathrm{L}}\left[1-H\left(\frac{E_{1,\mu }^{\mathrm{U}}}{X_{1,\mu }^{\mathrm{L}}}\right)\right]-f_{\mathrm{EC}}{Z}_{\mu }H\left(E_{\mathrm{tol}}\right),$where $Z_{1,\mu }^{\mathrm{L}}={\sum }_{k=0}^N{Z}_{1,\mu }^{\left(k\right)\mathrm{L}}/N$, $X_{1,\mu }^{\mathrm{L}}={\sum }_{k=0}^N{X}_{1,\mu }^{\left(k\right)\mathrm{L}}/N$, $E_{1,\mu }^{\mathrm{U}}={\sum }_{k=0}^N{E}_{1,\mu }^{\left(k\right)\mathrm{U}}/N$, f_EC is the error-correction efficiency, H(x) is the binary entropy function, E_tol is the expected bit error rate in the Z-basis, and $Z_{\mu }={\sum }_{k=0}^N{Z}_{\mu }^{\left(k\right)}/N$ for

27

$Z_{\mu }^{\left(k\right)}=\sum _{a_{k-\xi }\in A}\dots \sum _{a_{k-1}\in A}{Z}_{\left\{k-1,\xi -1\right\}\mu }.$Note that Z_μ is the expected ratio of rounds where both parties select the Z-basis, Alice selects the signal setting μ, and a detection event is registered. Also, we remark that $Z_{1,\mu }^{\mathrm{L}}$ incorporates the sifting factor $q_{\mathrm{Z}}^2$ by definition (see Eq. (16)). Of course, alternatively, one could consider a symmetric scenario where the parties extract key from both bases.

In the numerical simulations, we use the same channel and detector model deployed in refs. ^24,25, which is consistent with the reference parameters derived in section ‘Reference values for the linear CS constraints’. For ease of comparison with prior work, this model does not incorporate intensity correlations. The parameters of the model are ${\eta }_{\det }$ (detector efficiency of Bob’s detectors), p_d (dark count probability of Bob’s detectors), ${\eta }_{\mathrm{ch}}=10^{-{\alpha }_{\mathrm{att}}L/10}$ (channel transmittance for a lab-to-lab distance of L kilometres and an attenuation coefficient of α_att dB/km), and δ_A (misalignment error in the channel). The model reads

28

$\begin{array}{ccc}\hfill \frac{Z_{\left\{k,\xi \right\}}}{q_Z^2{\prod }_{i=k-\xi }^k{p}_{a_i}}& \hfill =\hfill & \frac{X_{\left\{k,\xi \right\}}}{q_X^2{\prod }_{i=k-\xi }^k{p}_{a_i}}=1-{\left(1-p_d\right)}^2{\mathrm{e}}^{-\eta a_k},\hfill \\ \hfill \frac{E_{\left\{k,\xi \right\}}}{q_X^2{\prod }_{i=k-\xi }^k{p}_{a_i}}& \hfill =\hfill & \frac{p_d^2}{2}+p_d\left(1-p_d\right)\left[1+h_{\eta ,{\delta }_A}\left(a_k\right)\right]\hfill \\ \hfill & \hfill \hfill & +{\left(1-p_d\right)}^2\left[\frac{1}{2}+h_{\eta ,{\delta }_A}\left(a_k\right)-\frac{1}{2}{\mathrm{e}}^{-\eta a_k}\right],\hfill \end{array}$for $\eta ={\eta }_{\mathrm{ch}}{\eta }_{\det }$ and $h_{\eta ,{\delta }_A}\left(a_k\right)=\left({\mathrm{e}}^{-\eta a_k{\cos }^2{\delta }_A}-{\mathrm{e}}^{-\eta a_k{\sin }^2{\delta }_A}\right)/2$. Naturally as well, E_tol = E_{{k, ξ}}/X_{{k, ξ}} because of the symmetry between the Z and the X bases.

In Fig. 2, we present numerical simulations of the asymptotic secret key rate attainable with our enhanced decoy-state method. For ease of comparison with prior work^24,25, we consider an asymmetric decoy-state BB84 protocol with an active receiver. The figure reveals that this method offers a great tolerance to variations in the magnitude and range of the correlations when compared to previous studies. In fact, it consistently achieves higher secret key rates for all considered scenarios.

[See PDF for image]

Fig. 2

Simulation of the asymptotic secret key rate in an asymmetric protocol with an active receiver.

a Secret key rate when Bob uses single-photon detectors (SPDs) with detection efficiency ${\eta }_{\det }=0.2$ and dark count rate p_d = 4.2 × 10⁻⁶. b Secret key rate when Bob uses superconducting nanowire SPDs (SNSPDs) with detection efficiency ${\eta }_{\det }=0.608$ and dark count rate p_d = 9.5 × 10⁻⁸. Solid lines, secret key rate obtained with our method; Dashed lines, secret key rate obtained with the most efficient previous method²⁵. We set ${\delta }_{\max ,a_k}^{\mathrm{corr}}={\delta }_{\max ,a_k}^{\mathrm{rand}}={\delta }_{\max }/2$ for ${\delta }_{\max }\in \left\{10^{-2},10^{-4}\right\}$, and contemplate two correlation ranges, ξ ∈ {1, 3}. We focus on the (asymptotically optimal) limit where p_μ = 1, p_ν = 0, p_ω = 0 and q_Z = 1, where in this case q_Z determines Bob’s probability to actively select the Z-basis as well. Also, we set the fibre attenuation coefficient to α = 0.2 dB/km, the misalignment to δ_A = 0.08 and the error correction efficiency to f_EC = 1.16.

Note that, in this work, we use the model here presented for the simulations of both Figs. 2 and 4. In this respect, we remark that, although the model is specifically suited to an active BB84 receiver, it can also be used to approximately describe the detection statistics of a passive receiver.

Experimental setup and implementation

The QKD setup consists of three modules. On the source side, as usual, there is a quantum module responsible for the BB84 and decoy-state encoding. Subsequently, the monitor module monitors the quantum module, capturing the intensity correlation information. The third module is Bob’s detection setup. Precisely, in the quantum module, for each protocol round k = 1, 2,…, N, Alice selects an intensity setting a_k ∈ A = {μ, ν, ω} with probability $p_{a_k}$, and a bit-and-basis setting r_k ∈ R = {Z0, Z1, X0, X1} with probability $p_{r_k}$. The module encodes these settings on PRWCPs that are attenuated by a beam splitter (BS) and a variable optical attenuator (VOA). Due to the effect of intensity correlations, the average photon-number of the k-th PRWCP is not exactly a_k, but α_k according to the general model described above.

The output of the quantum module consists of two parts, the first one being sent to Bob and the second one being sent to the monitor module. The monitor module contains a VOA and an SPD, and it classifies the detection results according to the record of the previous ξ settings, for the later calculation of (i) the bounds of the conditional n-photon probabilities of the different records, and (ii) the correlation parameters ${\tau }_{a_k,a_k^{\prime },n}^{\xi }$. For this purpose, the method described in the section ‘Calculation of correlation parameters’ is used.

In the detection module, Bob measures the received signals and announces whether a successful result occurs in his measurement unit (MU). He stores the rounds of successful measurements and the corresponding bases and raw key bits. For successful measurements, Alice stores her original key bit and basis, while she stores her decoy settings for all rounds to recognize the records {k, ξ} of every successful measurement. Box 0.1 illustrates the steps of the protocol, where the additional operations introduced to address the intensity correlation problem are highlighted in bold.

To validate our approach, we conducted experiments using the apparatus shown in Fig. 3. We use a polarization coding system where a gain-switched laser diode emits PRWCPs with a pulse width of 50 ps and a repetition rate of 1 GHz. The decoy-state intensity modulator (DS-IM) module includes an isolator, a Sagnac interferometer (SI)²⁰, and a commercial IM made from an integrated LiNbO₃ Mach-Zehnder interferometer (MZI). The SI employs a phase modulator (PM) to modulate clockwise and counterclockwise pulses, causing interference that leads to maximum and minimum attenuation when the relative phase is 0 and π, respectively. Due to the SI’s unique characteristics²⁰, it has two stable operating points during constructive and destructive interference, achieving smaller fluctuations and correlations compared to the commercial IM. SI is used to adjust the signal and decoy state, while the vacuum state is modulated using MZI, which also has stable operating points during constructive and destructive interference^21,27. By combining the SI and MZI IM, each with two stable operating points at maximum and minimum attenuation, we achieve the three required stable operating points for the DS-IM module. The specific interference combinations are shown in Table 1.

[See PDF for image]

Fig. 3

Depiction of the experimental setup.

LD laser diode, ISO isolator, PM phase modulator, IM intensity modulator, Pol.-Enc. polarization encoder, BS beam splitter, EVOA electronic variable optical attenuator, SPD single-photon detector, PC polarization controller, PBS polarization beam splitter.

Table 1. Relation between the intensity settings and the state of the interferometers

	Signal	Decoy	Vacuum
SI	Cons.	Des.	Des.
MZI	Cons.	Cons.	Des.

Cons. constructive interference, Des. destructive interference.

The polarization encoder encodes the DS-IM modulated pulses. A BS then divides the encoded states into two pulses, one of them entering the local monitor —composed of an electronic variable optical attenuator (EVOA-2) and a single-photon detector (SPD-M)— and the other one entering the fibre-based quantum channel after being attenuated with EVOA-1.

At the detection side, we opt for a passive BB84 receiver to measure the arriving polarization states. The detection module comprises one symmetric BS and two MUs. Each MU consists of a polarization controller (PC), a polarization BS (PBS), and two SPDs. The SPDs work on gated mode³¹ at a frequency of 1 GHz, with a dark count rate p_d = 4.2 × 10⁻⁶ and a detection efficiency ${\eta }_{\det }=20\%$. The PC-Z of MU-Z is properly adjusted so that SPD-H (V) measures horizontal (vertical) polarization states. In a similar fashion, MU-X measures the X basis polarization states.

The attenuation of Alice’s local monitor is controlled precisely to achieve the desired value of η_m. In our experiment, we set η_m to 10⁻³. SPD-M records whether or not each round clicks and categorizes all data into 3^ξ+1 types based on the record of settings. To illustrate our method, we have assumed that ξ = 3²⁷ for the data analysis, thus calculating the detection rates of 81 records of settings in total (see the Supplementary Note 2 for the specific calibration method).

Before performing the experiment, we conduct simulations to determine the intensities μ, ν and ω and their selection probabilities. We adjust the splitting ratio of the BS within SI to fix the signal-to-decoy intensities constrained to μ/ν = 5. We optimize both the intensities and their probabilities simultaneously in the simulation, while imposing the constraint $p_{a_k}\ge 0.15$ for all a_k to ensure that all data sets are sufficiently large. For consistency with the experiment, Alice’s Z-basis probability is set to q_Z = 0.5, and a symmetric passive receiver is considered at Bob’s side. In addition, we estimate the magnitude of the intensity correlations in advance^23,27 as ${\delta }_{\max ,\mu }^{\mathrm{corr}}=3.0\times 10^{-3}$, ${\delta }_{\max ,\nu }^{\mathrm{corr}}=2.0\times 10^{-3}$, ${\delta }_{\max ,\omega }^{\mathrm{corr}}=1.6\times 10^{-2}$ and ${\delta }_{\max ,\mu \left(\nu ,\omega \right)}^{\mathrm{rand}}=3.0\times 10^{-2}$, and we set the remaining parameters as in Fig. 2. Notably, prior knowledge of the device’s correlations and fluctuations is not essential for implementing the protocol, but it allows to optimize the secret key rate for each distance. Furthermore, this knowledge can be acquired by Alice’s monitor module before initiating the quantum communication.

The results of the simulation are shown in Fig. 4. We successfully demonstrate the feasibility of our approach at 50 km, 60 km, and 70 km of channel length. As depicted in the figure, the asymptotic key rates are respectively given by 1.82 × 10⁻⁴, 8.13 × 10⁻⁵, and 1.91 × 10⁻⁵, suitably aligned with the corresponding simulation outcomes (see the Supplementary Table 1 to 3 and Supplementary Figs. 1 to 3 in Supplementary Note 1 for detail data).

[See PDF for image]

Fig. 4

Simulation and experimental performance of the intensity-correlation-tolerant QKD protocol in Box 0.1.

Solid lines, secret key rate simulated for the experiment; Dashed lines, secret key rate of the case without correlation and fluctuation. The table inside the figure provides the key parameters of the experiments, where Q_Z, E_Z, y₁ and e_ph correspond, respectively, to the Z basis gain, the Z basis QBER, the Z basis single-photon yield and the phase error rate.

Box 0.1 Protocol description

State preparation: In the k-th round (k = 1, 2,…,N), Alice selects an intensity setting ak ∈ A = {μ, ν, ω} with probability $p_{a_k}$, and a bit-and-basis setting rk ∈ R = {Z0, Z1, X0, X1} with probability $p_{r_k}$, where pZ0 = pZ1 = qZ/2 and pX0 = pX1 = qX/2 for qZ + qX = 1. Then, she prepares a PRWCP accordingly, and a local monitor measures a fraction of the generated PRWCP with a photodetector with relative efficiency ηm. Measurement: Bob measures the received states and announces the successful measurement results. Alice and Bob record these later rounds and their corresponding bases and raw key bits. Also, Alice stores the intensity settings for all rounds. Sifting: Alice and Bob broadcast their basis choices for each round. They keep the raw key bits for the Z(X)-basis coincidences to form the sifted key (parameter estimation data). Parameter estimation: Bob discloses his measurement outcomes for the parameter estimation data. Alice classifies the results recorded by her local monitor in the state preparation step into 3ξ+1 categories —matching the possible records {k, ξ}— and computes the correlation parameters${\tau }_{a_k,a_k^{\prime },n}^{\xi }$. With this information, she uses the enhanced decoy-state method to estimate the Z and X basis single-photon yields, together with the phase error rate in the Z basis. Key distillation: Alice and Bob perform error correction and privacy amplification based on the results of the parameter estimation step, and then use the data in the Z basis to generate the secret keys.

Discussion

In summary, we have proposed an efficient solution to the intensity correlation problem in practical QKD. Our approach achieves high tolerance to this imperfection through the addition of a monitor module at Alice’s side at minimal cost, leading to higher secret key rates and extended maximum achievable distances compared to previous proposals. Leveraging a stable DS-IM module, we have conducted the first experiment addressing intensity correlations in a practical QKD setup without requiring hardware-based suppression methods²⁷, or limiting software-based methods restricted to the nearest-neighbours scenario²¹. Although we have considered a finite correlation range, our results could be promoted to the unbounded range setting using the tools in ref. ³². For future research, it is essential from an experimental perspective to improve the IM module in order to achieve lower levels of fluctuations and correlation. The current design still requires a reduction of fluctuations and correlation magnitude by two to three orders of magnitude to attain a key rate comparable to the ideal case. In addition, developing feedback methods based on low-noise photondiodes to account for time-dependent imperfections is another important research direction. On the theoretical side, analyzing potential security vulnerabilities arising from the source monitoring module would be valuable. Furthermore, it is meaningful to incorporate additional imperfections, such as side-channel effects, into the analysis.

Given that this type of correlation poses a considerable security concern in QKD systems, our work constitutes a significant advance for the practical security of QKD and facilitates the adoption of QKD in large-scale security applications.

Methods

Intensity correlation parameter

In this section, we outline the calculation of these inner products, and for this purpose, we restore the explicit notation of the intensity-setting subscripts for clarity. For further technical explanations, the reader is referred to ref. ²⁴.

In the entanglement-based picture, the global input state of all protocol rounds can be described as

29

$\mid \mathrm{\Psi }\left(\right.⟩=\left[{\sum }_{a_1\dots a_N}{\sum }_{r_1\dots r_N}\left(\prod _{i=1}^N\sqrt{p_{a_i}{q}_{r_i}}\right)\left(\underset{i=1}{\overset{N}{⨂}}{\mid a_i\left(\right.⟩}_{A_i}{\mid r_i\left(\right.⟩}_{{A_i}^{\prime }}{\mid {\psi }_{a_1\dots a_i}^{r_i}\left(\right.⟩}_{B_i{C}_i}\right)\right]\otimes {\mid 0\left(\right.⟩}_E,$where ${\mid a_i\left(\right.⟩}_{A_i}$ is a virtual ancilla storing the intensity setting in round i, ${\mid r_i\left(\right.⟩}_{{A_i}^{\prime }}$ is a virtual ancilla storing the encoded BB84 state in round i (i.e., the bit and basis information), $p_{a_i}$ is the probability of choosing the intensity setting a_i in round i, $q_{r_i}$ is the probability of choosing the BB84 state r_i in round i, and ${\mid 0\left(\right.⟩}_E$ is the vacuum state of Eve’s system. Also, for any given round i, we have defined

30

${\mid {\psi }_{a_1\dots a_i}^{r_i}\left(\right.⟩}_{B_i{C}_i}=\sum _{n_i=0}^{\infty }\sqrt{p_{n_i}{\mid }_{a_1\dots a_i}}{\mid t_{n_i}\left(\right.⟩}_{C_i}{\mid n_i^{r_i}\left(\right.⟩}_{B_i},$where $p_{n_i}{\mid }_{a_1\dots a_i}$ is the conditional n_i-photon probability given the record a₁…a_i, ${\mid t_{n_i}\left(\right.⟩}_{C_i}$ is a virtual ancilla storing the photon number n_i, and ${\mid n_i^{r_i}\left(\right.⟩}_{B_i}$ is a Fock state with n_i photons encoding the BB84 state r_i. Note that, given the finite range ξ of the correlations, for all i > ξ, we can replace $p_{n_i}{\mid }_{a_1\dots a_i}$ by $p_{n_i}{\mid }_{a_{i-\xi }\dots a_i}$ and ${\mid {\psi }_{a_1\dots a_i}^{r_i}\left(\right.⟩}_{B_i{C}_i}$ by ${\mid {\psi }_{a_{i-\xi }\dots a_i}^{r_i}\left(\right.⟩}_{B_i{C}_i}$.

In principle, for any given round k, one could establish CS constraints between the detection statistics of any two arbitrary records of settings. Nevertheless, this exhaustive approach would result in a number of CS constraints that increase exponentially with the correlation range, which is computationally prohibitive. Instead, we follow the simpler analysis presented in ref. ²⁵, where only the tightest CS constraints are considered. Particularly, for round k, one only relates the yields/error yields of records of intensity settings exclusively differing in a_k. Let a_k−ξ,…,a_k−1a_k and $a_{k-\xi }\dots a_{k-1}{a}_k^{\prime }$ be any two such records. As shown in refs. ^24,25, it turns out that the relevant overlap to compute in order to evaluate the CS constraints between these two records is

31

$∣⟨{\psi }_{a_{k-\xi }\dots a_{k-1}{a}_k,n}\mid {\psi }_{a_{k-\xi }\dots a_{k-1}{a}_k^{\prime },n}⟩∣,$where

32

$\mid {\psi }_{a_{k-\xi }\dots a_{k-1}{a}_k,n}\left(\right.⟩=\frac{\mid {\tilde{\psi }}_{a_{k-\xi }\dots a_{k-1}{a}_k,n}\left(\right.⟩}{∥\mid {\tilde{\psi }}_{a_{k-\xi }\dots a_{k-1}{a}_k,n}\left(\right.⟩∥},$for $\mid {\tilde{\psi }}_{a_{k-\xi }\dots a_{k-1}{a}_k,n}\left(\right.⟩={\left.⟨\right)a_{k-\xi }\mid }_{A_{k-\xi }}\otimes \dots \otimes {\left.⟨\right)a_{k-1}\mid }_{A_{k-1}}\otimes {\left.⟨\right)a_k\mid }_{A_k}\otimes {\left.⟨\right)t_n\mid }_{C_k}\mid \mathrm{\Psi }\left(\right.⟩$. Explicit calculation of this overlap follows identically as in ref. ²⁵ and yields

33

$\begin{array}{c}∣\left.⟨\right){\psi }_{a_{k-\xi }\dots a_{k-1}{a}_k,n}\mid \mid {\psi }_{a_{k-\xi }\dots a_{k-1}{a}_k^{\prime },n}\left(\right.⟩∣\hfill \\ =\sum _{a_{k+1}\in A}\dots \sum _{a_{k+\xi }\in A}\left(\prod _{i=k+1}^{k+\xi }{p}_{a_i}\right)\hfill \\ \left(\sum _{m=0}^{\infty }\sqrt{p_m{\mid }_{a_{i-\xi }\dots a_k\dots a_i}{p}_m{\mid }_{a_{i-\xi }\dots a_k^{\prime }\dots a_i}}\right).\hfill \end{array}$As one would expect, the overlap is dependent on the correlation function g_{{k−1, ξ−1}} through the conditional photon-number statistics. Since g_{{k−1, ξ−1}} is unknown, one must derive correlation-function-independent lower bounds on the photon-number statistics —say, $p_n^{\mathrm{L}}{\mid }_{a_{k-\xi }\dots a_k}$ for the specific record a_k−ξ,…,a_k— in order to possibly lower bound the overlap in Eq. (33). In this respect, we recall that lower bounds on the overlaps are necessary to reach loosened but correlation-function-independent CS constraints²⁵. In fact, such loosened CS constraints rely on lower bounds on the squared overlaps, which we refer to as the intensity correlation parameters,

34

$\begin{array}{c}{\tau }_{a_k,a_k^{\prime },n}^{\xi }:=\left\{\sum _{a_{k+1}\in A}\dots \sum _{a_{k+\xi }\in A}\right)\hfill \\ {\left(\left(\prod _{i=k+1}^{k+\xi }{p}_{a_i}\sum _{m=0}^{\infty }\sqrt{p_m^{\mathrm{L}}{\mid }_{a_{i-\xi }\dots a_k\dots a_i}{p}_m^{\mathrm{L}}{\mid }_{a_{i-\xi }\dots a_k^{\prime }\dots a_i}}\right)\right\}}^2\hfill \end{array},$for any given round k, any pair of distinct settings a_k and $a_{k^{\prime }}$, and any photon number n.

Reference values for the linear CS constraints

In this section, we present the reference parameters ${\tilde{y}}_{n,\left\{k-1,\xi -1\right\}{a}_k}$ and ${\tilde{h}}_{n,\left\{k-1,\xi -1\right\}{a}_k}$ that we utilize in the CS constraints. As mentioned earlier, although these parameters can be chosen arbitrarily, they have an impact on the tightness of the CS constraints. In accordance with ref. ²⁴, here we select them following a rather typical channel and detection model free of intensity correlations. The model is depicted in Fig. 5, and the main elements involved are described in the caption.

[See PDF for image]

Fig. 5

Illustration of the channel model that we use to compute the reference parameters for the CS constraints.

PRWCP phase-randomized weak coherent pulse, BS beam splitter, η overall transmittance (including channel transmittance and detection efficiency); σ_a, misalignment; D_err, detector associated to the “error” clicks; D_cor, detector associated to the “correct” (or “no error'') clicks.

The reference n-photon yield in the model is simply given by

35

${\tilde{y}}_{n,\left\{k-1,\xi -1\right\}{a}_k}=1-{\left(1-p_d\right)}^2{\left(1-\eta \right)}^n,$where p_d is the dark count probability of Bob’s detectors and η is the overall transmittance of the system.

In order to write down the reference n-photon error probabilities, we average over all four combinations of genuine detection events (i.e., disregarding dark counts for the moment). These have associated probabilities p₁ = (1−η)ⁿ (no detector clicks), $p_2={\left(1-{\cos }^2{\sigma }_a\eta \right)}^n-p_1$ (only D_cor clicks), $p_3={\left(1-{\sin }^2{\sigma }_a\eta \right)}^n-p_1$ (only D_err clicks) and p₄ = 1 − p₁ − p₂ − p₃ (both detectors click). Assuming that double clicks are randomly assigned to a specific detection outcome, one can readily identify the specific dark count combinations that map each of these four events to a bit error. In doing so, one obtains the following conditional error probabilities for all four events: $p_1^e=p_d\left(1-p_d\right)+\frac{1}{2}{p}_d^2$, $p_2^e={\left(1-p_d\right)}^2+\frac{3}{2}{p}_d\left(1-p_d\right)+\frac{1}{2}{p}_d^2$, $p_3^e=\frac{1}{2}{p}_d\left(1-p_d\right)+\frac{1}{2}{p}_d^2$ and $p_4^e=\frac{1}{2}{\left(1-p_d\right)}^2+p_d\left(1-p_d\right)+\frac{1}{2}{p}_d^2$, and therefore we conclude that

36

${\tilde{h}}_{n,\left\{k-1,\xi -1\right\}{a}_k}=\sum _{i=1}^4{p}_i{p}_i^e.$

Measurement of correlation parameters

The SPD is chosen for the calculation of correlation parameters based on several key reasons. First, Alice’s local monitoring requires security parameters, specifically intensity correlation parameters, which SPDs provide more accurately due to lower noise than conventional PDs. Second, the threshold response model of SPDs can be used for a tighter estimation of the photon-number statistics, which will be shown below. Despite limited detection efficiency, SPDs can still accurately measure intensity correlations by grouping click outcomes based on previous round settings and calculating response rates.

Employing the Taylor expansion e^x < 1 + x + x²/2 (−1 < x < 0) in Eq. (19), we obtain

37

$\begin{array}{c}{\mathcal{D}}_{\left\{k,\xi \right\}}<{\int }_{{\delta }_{\left\{k,\xi \right\}}^{-}}^{{\delta }_{\left\{k,\xi \right\}}^{+}}{g}_{\left\{k,\xi \right\}}\left(\delta \right)\hfill \\ \times \left[1-{\eta }_{\mathrm{m}}{\overline{\alpha }}_{\left\{k,\xi \right\}}\left(1+\delta \right)+\frac{{\eta }_{\mathrm{m}}^2{\overline{\alpha }}_{\left\{k,\xi \right\}}^2{\left(1+\delta \right)}^2}{2}\right]\mathrm{d}\delta ,\hfill \end{array}$and since the first-order term vanishes (see Eq. (2)), we have that

38

${\mathcal{D}}_{\left\{k,\xi \right\}}<1-{\eta }_{\mathrm{m}}{\overline{\alpha }}_{\left\{k,\xi \right\}}+{\eta }_{\mathrm{m}}^2{\overline{\alpha }}_{\left\{k,\xi \right\}}^2\frac{{\zeta }_{\left\{k,\xi \right\}}+1}{2},$for ${\zeta }_{\left\{k,\xi \right\}}=\max \left\{{\left({\delta }_{\left\{k,\xi \right\}}^{-}\right)}^2,{\left({\delta }_{\left\{k,\xi \right\}}^{+}\right)}^2\right\}$. Solving this quadratic equation for ${\overline{\alpha }}_{\left\{k,\xi \right\}}$ we find the upper bound

39

${\overline{\alpha }}_{\left\{k,\xi \right\}}<{\overline{\alpha }}_{\left\{k,\xi \right\}}^{\mathrm{U}}=\frac{1-\sqrt{1-2\left(1-{\mathcal{D}}_{\left\{k,\xi \right\}}\right)\left(1+{\zeta }_{\left\{k,\xi \right\}}\right)}}{{\eta }_{\mathrm{m}}\left(1+{\zeta }_{\left\{k,\xi \right\}}\right)}.$Similarly, keeping the third order we have that e^x > 1 + x + x²/2 + x³/6 (−1 < x < 0), and considering the practical situation of QKD, we have $1-{\eta }_m{\overline{\alpha }}_{\left\{k,\xi \right\}}\left(1-\frac{\delta }{3}\right)>0$ when $0<{\eta }_m,∣\delta ∣\ll 1$. This means that

40

$\begin{array}{c}{\mathcal{D}}_{\left\{k,\xi \right\}}>{\int }_{{\delta }_{\left\{k,\xi \right\}}^{-}}^{{\delta }_{\left\{k,\xi \right\}}^{+}}{g}_{\left\{k,\xi \right\}}\left(\delta \right)\left[1-{\eta }_{\mathrm{m}}{\overline{\alpha }}_{\left\{k,\xi \right\}}\left(1+\delta \right)\right)\hfill \\ +\left(\frac{{\eta }_{\mathrm{m}}^2{\overline{\alpha }}_{\left\{k,\xi \right\}}^2{\left(1+\delta \right)}^2}{2}-\frac{{\eta }_{\mathrm{m}}^3{\overline{\alpha }}_{\left\{k,\xi \right\}}^3{\left(1+\delta \right)}^3}{6}\right]\mathrm{d}\delta \hfill \\ >1-{\eta }_{\mathrm{m}}{\overline{\alpha }}_{\left\{k,\xi \right\}}+\frac{{{\eta }_{\mathrm{m}}}^2{\overline{\alpha }}_{\left\{k,\xi \right\}}^2}{2}-\frac{{{\eta }_{\mathrm{m}}}^3{\overline{\alpha }}_{\left\{k,\xi \right\}}^3}{6}.\hfill \end{array}$Coming next, the second order term can be lower-bounded via ${\overline{\alpha }}_{\left\{k,\xi \right\}}\ge \left(1-{\mathcal{D}}_{\left\{k,\xi \right\}}\right)/{\eta }_{\mathrm{m}}$ —which results from plugging the relation e^x > 1 + x (−1 < x < 0) in Eq. (19) — and the third order term can be lower bounded via ${\overline{\alpha }}_{\left\{k,\xi \right\}}<{\overline{\alpha }}_{\left\{k,\xi \right\}}^{\mathrm{U}}$. Altogether, this yields

41

${\overline{\alpha }}_{\left\{k,\xi \right\}}>{\overline{\alpha }}_{\left\{k,\xi \right\}}^{\mathrm{L}}=\frac{1-{\mathcal{D}}_{\left\{k,\xi \right\}}}{{\eta }_{\mathrm{m}}}+\frac{{\left(1-{\mathcal{D}}_{\left\{k,\xi \right\}}\right)}^2}{2{\eta }_{\mathrm{m}}}-\frac{{{\eta }_{\mathrm{m}}}^2{\left({\overline{\alpha }}_{\left\{k,\xi \right\}}^{\mathrm{U}}\right)}^3}{6}.$The above discussion leads to the conclusion in Eq. (20).

Starting from Eq. (3) and expanding both $\exp \left(-{\overline{\alpha }}_{\left\{k,\xi \right\}}\left(1+\delta \right)\right)$ and ${\left[{\overline{\alpha }}_{\left\{k,\xi \right\}}\left(1+\delta \right)\right]}^n$ in δ yields

42

$\begin{array}{c}{p}_n{\mid }_{\left\{k,\xi \right\}}=\frac{{\overline{\alpha }}_{\left\{k,\xi \right\}}^n{\mathrm{e}}^{-{\overline{\alpha }}_{\left\{k,\xi \right\}}}}{n!}{\int }_{{\delta }_{\left\{k,\xi \right\}}^{-}}^{{\delta }_{\left\{k,\xi \right\}}^{+}}{g}_{\left\{k,\xi \right\}}\left(\delta \right)\left[1-{\overline{\alpha }}_{\left\{k,\xi \right\}}\delta +o\left({\delta }^2\right)\right]\hfill \\ \left[1+n\delta +o\left({\delta }^2\right)\right]\mathrm{d}\delta \hfill \\ =\frac{{\overline{\alpha }}_{\left\{k,\xi \right\}}^n{\mathrm{e}}^{-{\overline{\alpha }}_{\left\{k,\xi \right\}}}}{n!}{\int }_{{\delta }_{\left\{k,\xi \right\}}^{-}}^{{\delta }_{\left\{k,\xi \right\}}^{+}}{g}_{\left\{k,\xi \right\}}\left(\delta \right)\left[1+\left(n-{\overline{\alpha }}_{\left\{k,\xi \right\}}\right)\delta +o\left({\delta }^2\right)\right]\mathrm{d}\delta .\hfill \end{array}$Since the first-order term $\left(n-{\overline{\alpha }}_{\left\{k,\xi \right\}}\right)\delta$ vanishes upon integration, one can explicitly suppress it from the definition of p_n∣_{{k, ξ}} to obtain

43

$\begin{array}{c}{p}_n{\mid }_{\left\{k,\xi \right\}}\hfill \\ =p_n{\mid }_{\left\{k,\xi \right\}}-\frac{{\overline{\alpha }}_{\left\{k,\xi \right\}}^n{\mathrm{e}}^{-{\overline{\alpha }}_{\left\{k,\xi \right\}}}}{n!}{\int }_{{\delta }_{\left\{k,\xi \right\}}^{-}}^{{\delta }_{\left\{k,\xi \right\}}^{+}}{g}_{\left\{k,\xi \right\}}\left(\delta \right)\left(n-{\overline{\alpha }}_{\left\{k,\xi \right\}}\right)\delta \mathrm{d}\delta \hfill \\ =\frac{{\overline{\alpha }}_{\left\{k,\xi \right\}}^n{\mathrm{e}}^{-{\overline{\alpha }}_{\left\{k,\xi \right\}}}}{n!}\hfill \\ {\int }_{{\delta }_{\left\{k,\xi \right\}}^{-}}^{{\delta }_{\left\{k,\xi \right\}}^{+}}{g}_{\left\{k,\xi \right\}}\left(\delta \right)\left[{\mathrm{e}}^{-\delta {\overline{\alpha }}_{\left\{k,\xi \right\}}}{\left(1+\delta \right)}^n-\left(n-{\overline{\alpha }}_{\left\{k,\xi \right\}}\right)\delta \right]\mathrm{d}\delta \hfill \\ ={\int }_{{\delta }_{\left\{k,\xi \right\}}^{-}}^{{\delta }_{\left\{k,\xi \right\}}^{+}}{g}_{\left\{k,\xi \right\}}\left(\delta \right)f_n\left(\delta ,{\overline{\alpha }}_{\left\{k,\xi \right\}}\right)\mathrm{d}\delta ,\hfill \end{array}$for

44

$f_n\left(x,y\right)=\frac{y^n{\mathrm{e}}^{-y}}{n!}\left[{\mathrm{e}}^{-xy}{\left(1+x\right)}^n-\left(n-y\right)x\right].$This methodology, which is also exploited in ref. ²³, allows for a more accurate estimation of the photon-number statistics than the monotonicity arguments provided in refs. ^24,25. Particularly, it follows from Eq. (43) that

45

$\begin{array}{c}\hfill \underset{x,y}{\min }{f}_n\left(x,y\right)\le p_n{\mid }_{\left\{k,\xi \right\}}\le \underset{x,y}{\max }{f}_n\left(x,y\right),\\ \hfill \mathrm{for}x\in \left[{\delta }_{\left\{k,\xi \right\}}^{-},{\delta }_{\left\{k,\xi \right\}}^{+}\right],y\in \left[{\overline{\alpha }}_{\left\{k,\xi \right\}}^{\mathrm{L}},{\overline{\alpha }}_{\left\{k,\xi \right\}}^{\mathrm{U}}\right].\end{array}$The above discussion leads to the conclusion in Eq. (21).

However, although the analysis in section ‘Calculation of correlation parameters’ and the above analysis assumes that η_m can be accurately characterized, if calibration errors of η_m do occur, we can still analyze and estimate the conditional photon-number statistics and correlation parameters. If the range of η_m can be determined, it can be incorporated into Eq. (20) to estimate new upper and lower bounds for ${\overline{\alpha }}_{\left\{k,\xi \right\}}$. By combining this with Eq. (25), we can obtain updated intensity correlation parameters. In subsequent analysis, calculations would be based on the worst-case result within the range of η_m. Additionally, we have numerically verified that an exact characterization of η_m is not essential. A relative deviation in η_m only results in a similarly scaled relative deviation in the secret key rate.

Acknowledgements

This work has been supported by the National Natural Science Foundation of China (Grant nos. 62271463, 62301524, 62105318, 61961136004, 62171424), the Fundamental Research Funds for the Central Universities, the China Postdoctoral Science Foundation (Grant nos. 2022M723064), Natural Science Foundation of Anhui (No. 2308085QF216), the Innovation Programme for Quantum Science and Technology (Grant no. 2021ZD0300700), and the Youth Innovation Fund of the University of Science and Technology of China. V.Z. and M.C. acknowledge support from 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 the grant No. PID2020-118178RB-C21, 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, the European Union’s Horizon Europe Framework Programme under the Marie Sklodowska-Curie Grant No. 101072637 (Project QSI) and the project “Quantum Security Networks Partnership” (QSNP, grant agreement no. 101114043). We would also like to express our sincere gratitude to Matej Pivoluska and his team for their attention to detail in our manuscript, which has contributed to improving its rigor.

Author contributions

J.L. and F.L. start the project. J.L., V.Z., F.L. and Z.Y. provide the theoretical analysis. F.L. Z.W. and J.L. design and perform the experiment and complete the data analysis. J.L., V.Z., M.C. and F.L. mainly prepared the manuscript. Z.Y., M.C., S.W., W.C., D.H., G.G. and Z.H. supervise the project.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request. And the detailed experimental data is provided within the Supplementary Material.

Competing interests

The authors declare no competing interests.

Supplementary information

The online version contains supplementary material available at https://doi.org/10.1038/s41534-025-01059-0.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Bennett, C. H. & Brassard, G. Quantum cryptography: public key distribution and coin tossing. In Conf. on Computers, Systems and Signal Processing, 175 (Bangalore, 1984).

2. Pirandola, S et al. Advances in quantum cryptography. Adv. Opt. Photon.; 2020; 12, pp. 1012-1236. [DOI: https://dx.doi.org/10.1364/AOP.361502]

3. Lo, H-K; Chau, HF. Unconditional security of quantum key distribution over arbitrarily long distances. Science; 1999; 283, pp. 2050-2056. [DOI: https://dx.doi.org/10.1126/science.283.5410.2050]

4. Shor, PW; Preskill, J. Simple proof of security of the BB84 quantum key distribution protocol. Phys. Rev. Lett.; 2000; 85, pp. 441-444. [DOI: https://dx.doi.org/10.1103/PhysRevLett.85.441]

5. Scarani, V et al. The security of practical quantum key distribution. Rev. Mod. Phys.; 2009; 81, pp. 1301-1350. [DOI: https://dx.doi.org/10.1103/RevModPhys.81.1301]

6. Renner, R. Security of quantum key distribution. Int. J. Quantum Inf.; 2008; 06, pp. 1-127. [DOI: https://dx.doi.org/10.1142/S0219749908003256]

7. Hwang, W-Y. Quantum key distribution with high loss: Toward global secure communication. Phys. Rev. Lett.; 2003; 91, 057901. [DOI: https://dx.doi.org/10.1103/PhysRevLett.91.057901]

8. Wang, X-B. Beating the photon-number-splitting attack in practical quantum cryptography. Phys. Rev. Lett.; 2005; 94, 230503. [DOI: https://dx.doi.org/10.1103/PhysRevLett.94.230503]

9. Lo, H-K; Ma, X; Chen, K. Decoy state quantum key distribution. Phys. Rev. Lett.; 2005; 94, 230504. [DOI: https://dx.doi.org/10.1103/PhysRevLett.94.230504]

10. Takesue, H et al. Quantum key distribution over a 40-dB channel loss using superconducting single-photon detectors. Nat. Photonics; 2007; 1, pp. 343-348. [DOI: https://dx.doi.org/10.1038/nphoton.2007.75]

11. Lucamarini, M et al. Efficient decoy-state quantum key distribution with quantified security. Opt. Express; 2013; 21, pp. 24550-24565. [DOI: https://dx.doi.org/10.1364/OE.21.024550]

12. Yuan, Z et al. 10-Mb/s quantum key distribution. J. Lightw. Technol.; 2018; 36, pp. 3427-3433. [DOI: https://dx.doi.org/10.1109/JLT.2018.2843136]

13. Grünenfelder, F et al. Fast single-photon detectors and real-time key distillation enable high secret-key-rate quantum key distribution systems. Nat. Photonics; 2023; 17, pp. 422-426. [DOI: https://dx.doi.org/10.1038/s41566-023-01168-2]

14. Islam, NT; Lim, CCW; Cahall, C; Kim, J; Gauthier, DJ. Provably secure and high-rate quantum key distribution with time-bin qudits. Sci. Adv.; 2017; 3, e1701491. [DOI: https://dx.doi.org/10.1126/sciadv.1701491]

15. Boaron, A et al. Secure quantum key distribution over 421 km of optical fiber. Phys. Rev. Lett.; 2018; 121, 190502. [DOI: https://dx.doi.org/10.1103/PhysRevLett.121.190502]

16. Grünenfelder, F; Boaron, A; Rusca, D; Martin, A; Zbinden, H. Performance and security of 5 GHz repetition rate polarization-based quantum key distribution. Appl. Phys. Lett.; 2020; 117, 144003. [DOI: https://dx.doi.org/10.1063/5.0021468]

17. Li, W et al. High-rate quantum key distribution exceeding 110 Mb s⁻¹. Nat. Photonics; 2023; 17, pp. 416-421. [DOI: https://dx.doi.org/10.1038/s41566-023-01166-4]

18. Xu, F; Ma, X; Zhang, Q; Lo, H-K; Pan, J-W. Secure quantum key distribution with realistic devices. Rev. Mod. Phys.; 2020; 92, 025002.4122934 [DOI: https://dx.doi.org/10.1103/RevModPhys.92.025002]

19. Kobayashi, T; Tomita, A; Okamoto, A. Evaluation of the phase randomness of a light source in quantum-key-distribution systems with an attenuated laser. Phys. Rev. A; 2014; 90, 032320. [DOI: https://dx.doi.org/10.1103/PhysRevA.90.032320]

20. Roberts, GL et al. Patterning-effect mitigating intensity modulator for secure decoy-state quantum key distribution. Opt. Lett.; 2018; 43, pp. 5110-5113. [DOI: https://dx.doi.org/10.1364/OL.43.005110]

21. Yoshino, K-i et al. Quantum key distribution with an efficient countermeasure against correlated intensity fluctuations in optical pulses. npj Quantum Inf.; 2018; 4, [DOI: https://dx.doi.org/10.1038/s41534-017-0057-8] 8.

22. Lu, F-Y et al. Intensity modulator for secure, stable, and high-performance decoy-state quantum key distribution. npj Quantum Inf.; 2021; 7, [DOI: https://dx.doi.org/10.1038/s41534-021-00418-x] 75.

23. Lu, F-Y et al. Intensity tomography method for secure and high-performance quantum key distribution. J. Lightw. Technol.; 2023; 41, pp. 4895-4900. [DOI: https://dx.doi.org/10.1109/JLT.2023.3247766]

24. Zapatero, V; Navarrete, Á; Tamaki, K; Curty, M. Security of quantum key distribution with intensity correlations. Quantum; 2021; 5, 602. [DOI: https://dx.doi.org/10.22331/q-2021-12-07-602]

25. Sixto, X; Zapatero, V; Curty, M. Security of decoy-state quantum key distribution with correlated intensity fluctuations. Phys. Rev. Appl.; 2022; 18, 044069. [DOI: https://dx.doi.org/10.1103/PhysRevApplied.18.044069]

26. Pereira, M; Kato, G; Mizutani, A; Curty, M; Tamaki, K. Quantum key distribution with correlated sources. Sci. Adv.; 2020; 6, eaaz4487. [DOI: https://dx.doi.org/10.1126/sciadv.aaz4487]

27. Kang, X et al. Patterning-effect calibration algorithm for secure decoy-state quantum key distribution. J. Lightw. Technol.; 2023; 41, pp. 75-82. [DOI: https://dx.doi.org/10.1109/JLT.2022.3211442]

28. Huang, A; Mizutani, A; Lo, H-K; Makarov, V; Tamaki, K. Characterization of state-preparation uncertainty in quantum key distribution. Phys. Rev. Appl.; 2023; 19, 014048. [DOI: https://dx.doi.org/10.1103/PhysRevApplied.19.014048]

29. Xie, H-B et al. Optically injected intensity-stable pulse source for secure quantum key distribution. Opt. Express; 2019; 27, pp. 12231-12240. [DOI: https://dx.doi.org/10.1364/OE.27.012231]

30. Wang, X-B. Decoy-state quantum key distribution with large random errors of light intensity. Phys. Rev. A; 2007; 75, 052301. [DOI: https://dx.doi.org/10.1103/PhysRevA.75.052301]

31. He, D-Y et al. Sine-wave gating InGaAs/InP single photon detector with ultralow afterpulse. Appl. Phys. Lett.; 2017; 110, 111104. [DOI: https://dx.doi.org/10.1063/1.4978599]

32. Pereira, M. et al. Quantum key distribution with unbounded pulse correlations. Quantum Sci. Technol.10, 015001 (2024).