Sending-or-not-sending twin-filed quantum key distribution with discrete phase modulation

The implementation process of the SNS protocol with discrete phase modulation WCS sources is similar to that of the original SNS protocol Wang et al. (2018). Here we first introduce the 4-intensity protocol as follows. Obviously, the special case that the intensity of signal state equals to that of one of the decoy state in the 4-intensity makes the 3-intensity protocol.

For each time window, Alice (Bob) randomly decides it is a decoy window or a signal window. If it is a decoy window, Alice (Bob) randomly chooses to prepare a pulse of state $|0\rangle$, $|e^{2m\pi i/N}\sqrt{\mu_{x}}\rangle$ or $|e^{2m^{\prime}\pi i/N}\sqrt{\mu_{y}}\rangle$ ($|0\rangle$, $|e^{2n\pi i/N}\sqrt{\mu_{x}}\rangle$ or $|e^{2n^{\prime}\pi i/N}\sqrt{\mu_{y}}\rangle$), where $m,m^{\prime},n,n^{\prime}$ are randomly chosen from $\{0,1,2,\cdots,N-1\}$ and $N$ is assumed to be an even number. If it is a signal window, Alice (Bob) randomly chooses to prepare a vacuum pulse or a pulse of state $|e^{2l\pi i/N}\sqrt{\mu_{z}}\rangle$ with probabilities $1-\epsilon$ and $\epsilon$ respectively, where $l$ is randomly chosen from $\{0,1,2,\cdots,N-1\}$. If Alice (Bob) decides to send out a vacuum pulse, that is to say sending nothing or not sending, Alice (Bob) takes the corresponding classical bit value as $0(1)$. If Alice (Bob) decides to send out a pulse of state $|e^{2l\pi i/N}\sqrt{\mu_{z}}\rangle$, Alice (Bob) takes the corresponding classical bit value as $1(0)$.

Then Alice and Bob send their pulses to Charlie, Charlie is assumed to perform interferometric measurements on the received pulses. If only one of the two detectors clicks, Charlie would announce this pulse pair causes a click and whether the left detector or right detector clicks. Alice and Bob take it as a one-detector heralded event.

After Alice and Bob repeat the above process for many times, they acquire a series of data which are used to perform the data post-processing.

The first step of data post-processing is sifting. Alice and Bob first announce the types of each time window they have decided. For a window that both Alice and Bob have decided a signal window, it is a $Z$ window. The corresponding bits of the one-detector heralded event of the $Z$ windows, which are also called as the sifted key, are used to extract the final key. Except for the $Z$ windows, Alice and Bob announce the intensities and phases they have chosen in each window. For a window that both Alice and Bob have decided a decoy window, and the intensity of the pulse is $\mu_{x}$ and their phases satisfy

$$1-|\cos(\frac{2m\pi}{N}-\frac{2n\pi}{N})|\leq\delta,$$

(1)

where $\delta$ is a positive number close to 0, it is an $X_{1}$ window. Here in our calculation we shall simply set

$$|\cos(\frac{2m\pi}{N}-\frac{2n\pi}{N})|=1$$

(2)

for the post selection condition of $X_{1}$ windows above. The data of $X_{1}$ windows are used to estimate the phase-flip error rate of untagged bits. The data of other windows are used for decoy-state analysis.

As the phases of WCS pulses in the $Z$ windows are never announced in the public channel, the density matrix of those WCS pulses is

$$\rho_{z}=\frac{1}{N}\sum_{l=0}^{N-1}|e^{2l\pi i/N}\sqrt{\mu_{z}}\rangle\langle e^{-2l\pi i/N}\sqrt{\mu_{z}}|.$$

(3)

For convenience, we define the approximated $j$-photon state $|\lambda_{j}\rangle$ in the following form Cao et al. (2015):

$$|\lambda_{j}\rangle=\frac{1}{\sqrt{P_{j}(\mu_{z})}}\sum_{k=0}^{\infty}\frac{(\sqrt{\mu_{z}})^{kN+j}}{\sqrt{(kN+j)!}}|kN+j\rangle,$$

(4)

and $j=0,1,2\cdots,N-1$. It is easy to see that $\rho_{z}$ is a classical mixture of different $|\lambda_{j}\rangle$. Explicitly, we have

$$\rho_{z}=\sum_{j=0}^{N-1}P_{j}(\mu_{z})|\lambda_{j}\rangle\langle\lambda_{j}|,$$

(5)

where

$$P_{j}(\mu_{z})=\sum_{k=0}^{\infty}\frac{\mu_{z}^{kN+j}e^{-\mu_{z}}}{(kN+j)!}.$$

(6)

For the pulse pairs in $X_{1}$ windows, if the phases of the pulse pair satisfy $n\equiv(m+q){\rm mod}N$ where $q$ is a constant integer, the density matrix of those pulse pairs is

$$\rho_{X1}(q)=\frac{1}{N}\sum_{m=0}^{N-1}|e^{2m\pi i/N}\sqrt{\mu_{x}}\rangle_{A}|e^{2n\pi i/N}\sqrt{\mu_{x}}\rangle_{B}\langle e^{-2m\pi i/N}\sqrt{\mu_{x}}|_{A}\langle e^{-2n\pi i/N}\sqrt{\mu_{x}}|_{B},$$

(7)

where the subscript $A$ and $B$ indicate Alice and Bob respectively. After a simple calculation, one can find that $\rho_{X1}(q)$ is actually the classical mixture of the state $|\varphi_{j}^{q}\rangle$,

$$\rho_{X1}(q)=\sum_{j=0}^{N-1}PX_{j}(\mu_{x})|\varphi_{j}^{q}\rangle\langle\varphi_{j}^{q}|$$

(8)

where

$$|\varphi_{j}^{q}\rangle=\frac{e^{-\mu_{x}}}{\sqrt{PX_{j}(\mu_{x})}}\sum_{k=0}^{\infty}\sum_{k_{1}=0}^{kN+j}\frac{(\sqrt{\mu_{x}})^{kN+j}}{\sqrt{k_{1}!(kN+j-k_{1})!}}e^{\frac{2\pi i}{N}q(kN+j-k_{1})}|k_{1};kN+j-k_{1}\rangle,$$

(9)

and the probability to obtain the state $|\varphi_{j}^{q}\rangle$ is

$$PX_{j}(\mu_{x})=\sum_{k=0}^{\infty}\sum_{k_{1}=0}^{kN+j}\frac{{\mu_{x}}^{kN+j}e^{-2\mu_{x}}}{{k_{1}!(kN+j-k_{1})!}}=\sum_{k=0}^{\infty}\frac{{(2\mu_{x})}^{kN+j}e^{-2\mu_{x}}}{{(kN+j)!}}=P_{j}(2\mu_{x}).$$

(10)

Follow the security proof in Ref. Wang et al. (2018), let’s first consider the equivalent entanglement protocol of the SNS protocol.

In the entanglement protocol, for each time window, Alice and Bob pre-share the entanglement state

$$\begin{split}|\Psi_{1}\rangle=&\frac{1}{\sqrt{2}}(|0\lambda_{1}\rangle\otimes|\tilde{0}\tilde{0}\rangle+|\lambda_{1}0\rangle\otimes|\tilde{1}\tilde{1}\rangle)\\ =&\frac{1}{\sqrt{2}}[\frac{1}{\sqrt{2}}(|0\lambda_{1}\rangle+|\lambda_{1}0\rangle)\otimes\frac{1}{\sqrt{2}}(|\tilde{0}\tilde{0}\rangle+|\tilde{1}\tilde{1}\rangle)+\frac{1}{\sqrt{2}}(|0\lambda_{1}\rangle-|\lambda_{1}0\rangle)\otimes\frac{1}{\sqrt{2}}(|\tilde{0}\tilde{0}\rangle-|\tilde{1}\tilde{1}\rangle)],\end{split}$$

(11)

where $|\tilde{0}\tilde{0}\rangle$ and $|\tilde{1}\tilde{1}\rangle$ are local states that only exist in Alice’s and Bob’s labs and $|0\lambda_{1}\rangle$ and $|\lambda_{1}0\rangle$ are real states that are sent to Charlie. According to the measurement results announced by Charlie, Alice and Bob can get a series of almost perfect entanglement state by applying entanglement purification to the local states. Two parameters are needed in the entanglement purification: the first is the bit-flip error rate in the $Z$ basis, $e_{z}$, and the second is the bit-flip error rate in the $X$ basis, $e^{ph}$, which is also the phase-flip error rate in the $Z$ basis. Here the $Z$ basis means $\{|\tilde{0}\rangle,|\tilde{1}\rangle\}$, and the $X$ basis means $\{\frac{1}{\sqrt{2}}(|\tilde{0}\rangle+|\tilde{1}\rangle),\frac{1}{\sqrt{2}}(|\tilde{0}\rangle-|\tilde{1}\rangle)\}$. Finally, by measuring the local states in the $Z$ basis, Alice and Bob can get secure final keys.

As Alice and Bob only concern about the secure final keys, they needn’t have to measure their local states after Charlie announces his measurement results, but they can just measure their local states before they send the real states to Charlie. If Alice and Bob each measures their local qubits in the $Z$ basis, it is equivalent to that Alice and Bob randomly send the pulse of state $|0\lambda_{1}\rangle$ or $|\lambda_{1}0\rangle$ to Charlie. If Alice and Bob each measures their local qubits in the $X$ basis, it is equivalent to that Alice and Bob randomly send the pulse of state $|\chi_{0}\rangle=\frac{1}{\sqrt{2}}(|0\lambda_{1}\rangle+|\lambda_{1}0\rangle)$ or $|\chi_{1}\rangle=\frac{1}{\sqrt{2}}(|0\lambda_{1}\rangle-|\lambda_{1}0\rangle)$ to Charlie. As shown in Ref. Wang et al. (2018), a phase error happens if Alice and Bob send $|\chi_{0}\rangle$ to Charlie and Charlie announces the right detector clicks or Alice and Bob send $|\chi_{1}\rangle$ to Charlie and Charlie announces the left detector clicks.

Denote $T_{0}^{R}$ as the probability that Charlie announces the right detectors clicks while Alice and Bob send $|\chi_{0}\rangle$ to Charlie. Denote $T_{1}^{L}$ as the probability that Charlie announces the left detectors clicks while Alice and Bob send $|\chi_{1}\rangle$ to Charlie. Denote $s_{1}$ as the yield of $\frac{1}{2}(|0\lambda_{1}\rangle\langle 0\lambda_{1}|+|\lambda_{1}0\rangle\langle\lambda_{1}0|)$. We can calculate the phase-flip error rate by

$$e^{ph}=\frac{T_{0}^{R}+T_{1}^{L}}{2s_{1}}.$$

(12)

According to the above discussion and the tagged-model, we can define the untagged bits in the real protocol as the bits in the $Z$ windows that Alice decides not to send and Bob actually sends a pulse of state $|\lambda_{1}\rangle\langle\lambda_{1}|$ while Bob decides to send a WCS pulse with intensity $\mu_{z}$, or Bob decides not to send and Alice actually sends a pulse of state $|\lambda_{1}\rangle\langle\lambda_{1}|$ while Alice decides to send a WCS pulse with intensity $\mu_{z}$. Finally, we can get the secure final key rate by

$$R=2\epsilon(1-\epsilon)P_{1}(\mu_{z})s_{1}[1-H(e^{ph})]-S_{z}fH(E),$$

(13)

where $H(x)=-x\log_{2}(x)-(1-x)\log_{2}(1-x)$ is the Shannon entropy, $S_{z}$ is the yield of the events in the $Z$ windows, $f$ is error correction efficiency factor, and $E$ is the bit error rate of the sifted keys.

To clearly show how to apply the decoy-state method to this protocol, we denote Alice’s sources $|0\rangle$, $|e^{2m\pi i/N}\sqrt{\mu_{x}}\rangle$, and $|e^{2m^{\prime}\pi i/N}\sqrt{\mu_{y}}\rangle$ by $o_{A}$, $x_{A}$, and $y_{A}$ respectively, and we also denote Bob’s sources $|0\rangle$, $|e^{2n\pi i/N}\sqrt{\mu_{x}}\rangle$, and $|e^{2n^{\prime}\pi i/N}\sqrt{\mu_{y}}\rangle$ by $o_{B}$, $x_{B}$, and $y_{B}$. We denote the source of pulse pairs by $\kappa_{A}\varsigma_{B}$, where $\kappa,\varsigma=o,x,y$. And for simplicity, we omit the subscripts. For example, source $ox$ represents that Alice uses the source $o_{A}$ and Bob uses the source $x_{B}$. Without phase post-selection, the density matrices of sources $x_{A},y_{A},x_{B},y_{B}$ have similar form with Eq. (4). Specifically, we have

$$\rho_{w}=\sum_{j=0}^{N-1}P_{j}(\mu_{w})|\lambda_{j}^{w}\rangle\langle\lambda_{j}^{w}|,\quad(w=x,y),$$

(14)

where

$$|\lambda_{j}^{w}\rangle=\frac{1}{\sqrt{P_{j}(\mu_{w})}}\sum_{k=0}^{\infty}\frac{(\sqrt{\mu_{w}})^{kN+j}}{\sqrt{(kN+j)!}}|kN+j\rangle,$$

(15)

and $P_{j}(\mu_{w})$ is defined in Eq. (6). Here $\rho_{x}$ is the density matrix of sources $x_{A}$ and $x_{B}$, and $\rho_{y}$ is the density matrix of sources $y_{A}$ and $y_{B}$.

Denote the yield of sources $\kappa\varsigma$ by $S_{\kappa\varsigma}$. Denote the yield of states $|0\lambda_{j}^{w}\rangle\langle 0\lambda_{j}^{w}|$ and $|\lambda_{j}^{w}0\rangle\langle\lambda_{j}^{w}0|$ by $Y_{vj}^{w}$ and $Y_{jv}^{w}$. In the original SNS protocol Wang et al. (2018), as the continuously modulated phase-randomized WCS sources are used, we have

$$Y_{vj}^{x}=Y_{vj}^{y},\quad Y_{jv}^{x}=Y_{jv}^{y},$$

(16)

but this equality no longer holds in this protocol. Consider the properties of trace distance Cao et al. (2015), we have

$$|Y_{vj}^{x}-Y_{vj}^{y}|\leq\sqrt{1-\left(F_{xy}^{j}\right)^{2}},\quad|Y_{jv}^{x}-Y_{jv}^{y}|\leq\sqrt{1-\left(F_{xy}^{j}\right)^{2}},$$

(17)

where

$$F_{xy}^{j}=\frac{|\langle\lambda_{j}^{x}|\lambda_{j}^{y}\rangle|}{\sqrt{\langle\lambda_{j}^{x}|\lambda_{j}^{x}\rangle\langle\lambda_{j}^{y}|\lambda_{j}^{y}\rangle}}=\frac{\sum_{k=0}^{\infty}\frac{(\mu_{x}\mu_{y})^{(kN+j)/2}}{(kN+j)!}}{\sqrt{\sum_{k=0}^{\infty}\frac{\mu_{x}^{kN+j}}{(kN+j)!}\sum_{k=0}^{\infty}\frac{\mu_{y}^{kN+j}}{(kN+j)!}}},$$

(18)

is the fidelity of states $|\lambda_{j}^{x}\rangle$ and $|\lambda_{j}^{y}\rangle$.

Denote the yield of states $|0\lambda_{1}\rangle$ and $|\lambda_{1}0\rangle$ by $s_{01}$ and $s_{10}$ respectively. And we have $s_{1}=\frac{1}{2}(s_{01}+s_{10})$. We can get the lower bound of $s_{1}$ by either analytical formula or the linear programming.

Combining the equations

$$S_{ox}=\sum_{j=0}^{N-1}P_{j}(\mu_{x})Y_{vj}^{x}=\sum_{j=0}^{N-1}P_{j}(\mu_{x})Y_{vj}^{y}+\Delta,\quad S_{oy}=\sum_{j=0}^{N-1}P_{j}(\mu_{y})Y_{vj}^{y},$$

(19)

where

$$\Delta=\sum_{j=0}^{N-1}P_{j}(\mu_{x})(Y_{vj}^{x}-Y_{vj}^{y}).$$

(20)

we have

$$Y_{v1}^{y}=\frac{P_{2}(\mu_{y})S_{ox}-P_{2}(\mu_{x})S_{oy}-[P_{0}(\mu_{x})P_{2}(\mu_{y})-P_{0}(\mu_{y})P_{2}(\mu_{x})]Y_{v0}^{y}-P_{2}(\mu_{y})\Delta-\xi}{P_{1}(\mu_{x})P_{2}(\mu_{y})-P_{1}(\mu_{y})P_{2}(\mu_{x})},$$

(21)

where

$$\xi=\sum_{j=3}^{N-1}[P_{j}(\mu_{x})P_{2}(\mu_{y})-P_{j}(\mu_{y})P_{2}(\mu_{x})]Y_{vj}^{y}.$$

(22)

It is easy to check that $\xi\leq 0$ if the following condition holds

$$\frac{P_{1}(\mu_{x})}{P_{1}(\mu_{y})}\geq\frac{P_{2}(\mu_{x})}{P_{2}(\mu_{y})}\geq\frac{P_{j}(\mu_{x})}{P_{j}(\mu_{y})},\quad j=3,4,\cdots,N-1,$$

(23)

which can be easily examined given values of $\mu_{x}$ and $\mu_{y}$. And in our numerical simulation, we found Eq. (23) always holds. With

$$Y_{v0}^{y}\leq S_{oo}+\sqrt{1-F_{0}^{2}},\quad\Delta\leq\Delta^{U}=\sum_{j=0}^{N-1}P_{j}(\mu_{x})\sqrt{1-\left(F_{xy}^{j}\right)^{2}},\quad s_{01}\geq Y_{v1}^{y}-\sqrt{1-F_{1}^{2}},$$

(24)

where

$$F_{0}=\frac{1}{\sqrt{\sum_{k=0}^{\infty}\frac{\mu_{y}^{kN}}{(kN)!}}},\quad F_{1}=\frac{\sum_{k=0}^{\infty}\frac{(\mu_{y}\mu_{z})^{(kN+1)/2}}{(kN+1)!}}{\sqrt{\sum_{k=0}^{\infty}\frac{\mu_{y}^{kN+1}}{(kN+1)!}\sum_{k=0}^{\infty}\frac{\mu_{y}^{kN+1}}{(kN+1)!}}}.$$

(25)

We have

$$s_{01}\geq s_{01}^{L}=\frac{P_{2}(\mu_{y})S_{ox}-P_{2}(\mu_{x})S_{oy}-[P_{0}(\mu_{x})P_{2}(\mu_{y})-P_{0}(\mu_{y})P_{2}(\mu_{x})]\left(S_{oo}+\sqrt{1-F_{0}^{2}}\right)-P_{2}(\mu_{y})\Delta^{U}}{P_{1}(\mu_{x})P_{2}(\mu_{y})-P_{1}(\mu_{y})P_{2}(\mu_{x})}-\sqrt{1-F_{1}^{2}}.$$

(26)

By similar method, we can prove that

$$s_{10}\geq s_{10}^{L}=\frac{P_{2}(\mu_{y})S_{xo}-P_{2}(\mu_{x})S_{yo}-[P_{0}(\mu_{x})P_{2}(\mu_{y})-P_{0}(\mu_{y})P_{2}(\mu_{x})]\left(S_{oo}+\sqrt{1-F_{0}^{2}}\right)-P_{2}(\mu_{y})\Delta^{U}}{P_{1}(\mu_{x})P_{2}(\mu_{y})-P_{1}(\mu_{y})P_{2}(\mu_{x})}-\sqrt{1-F_{1}^{2}},$$

(27)

if Eq. (23) holds. Finally, we have $s_{1}\geq\frac{1}{2}(s_{01}^{L}+s_{10}^{L})$.

The remaining task is to estimate the upper bound of phase-flip error rate, $e^{ph}$, which is equivalent to estimate the upper bounds of $T_{0}^{R}$ and $T_{1}^{L}$.

Denote $T_{+}^{R}$ as the probability that Charlie announces the right detector clicks while Alice and Bob send out the pulse pairs in the $X_{1}$ window and their phases satisfy $m=n$. Denote $T_{-}^{L}$ as the probability that Charlie announces the left detector clicks while Alice and Bob send out the pulse pairs in the $X_{1}$ window and their phases satisfy $m=(n+N/2){\rm mod}N$. Given the discussion above, we know that the pulse pairs in the $X_{1}$ window whose phases satisfy $m=n$ are the classical mixture of the state $|\varphi_{j}^{0}\rangle$. Thus by denoting $t_{j}^{R}$ as the probability that Charlie announces the right detector clicks while Alice and Bob send out the pulse pairs of state $|\varphi_{j}^{0}\rangle$, we have

$$T_{+}^{R}=\sum_{j=0}^{N-1}PX_{j}(\mu_{x})t_{j}^{R}.$$

(28)

Denote $T_{00}$ and $T_{00}^{\prime}$ as the probability that Charlie announces the right detector clicks and left detector clicks while Alice and Bob send out a vacuum pulse pair. Consider the properties of trace distance, we have

$$|t_{0}^{R}-T_{00}|\leq\sqrt{1-F_{00}^{2}},\quad|t_{1}^{R}-T_{0}^{R}|\leq\sqrt{1-\left[F_{11}^{+}(0)\right]^{2}},$$

(29)

where

$$F_{00}=\frac{|\langle\varphi_{0}^{0}|0\rangle|}{\sqrt{\langle\varphi_{0}^{0}|\varphi_{0}^{0}\rangle\langle 0|0\rangle}}=\frac{1}{\sqrt{\sum_{k=0}^{\infty}\frac{{(2\mu_{x})}^{kN}}{{(kN)!}}}},$$

(30)

and

$$\begin{split}&F_{11}^{+}(q)=\frac{|\langle\varphi_{1}^{q}|\chi_{0}\rangle|}{\sqrt{\langle\varphi_{1}^{q}|\varphi_{1}^{q}\rangle\langle\chi_{0}|\chi_{0}\rangle}}=\frac{\sqrt{Re_{+}^{2}+Im^{2}}}{\sqrt{\sum_{k=0}^{\infty}\frac{{(2\mu_{x})}^{kN+1}}{{(kN+1)!}}\sum_{k=0}^{\infty}\frac{\mu_{z}^{kN+1}}{(kN+1)!}}},\\ &Re_{+}=\frac{1}{\sqrt{2}}\sum_{k=0}^{\infty}\frac{(\mu_{x}\mu_{z})^{(kN+1)/2}}{(kN+1)!}[1+\cos\frac{2\pi}{N}q(kN+1)],\\ &Im=\frac{1}{\sqrt{2}}\sum_{k=0}^{\infty}\frac{(\mu_{x}\mu_{z})^{(kN+1)/2}}{(kN+1)!}\sin\frac{2\pi}{N}q(kN+1).\end{split}$$

(31)

Combine Eqs. (28) and (29), we have

$$T_{0}^{R}\leq T_{0}^{R,U}=\frac{T_{+}^{R}-PX_{0}(\mu_{x})\left(T_{00}-\sqrt{1-F_{00}^{2}}\right)}{PX_{1}(\mu_{x})}+\sqrt{1-\left[F_{11}^{+}(0)\right]^{2}}.$$

(32)

With similar method, we have

$$T_{1}^{L}\leq T_{1}^{L,U}=\frac{T_{-}^{L}-PX_{0}(\mu_{x})\left(T_{00}^{\prime}-\sqrt{1-F_{00}^{2}}\right)}{PX_{1}(\mu_{x})}+\sqrt{1-\left[F_{11}^{-}(\frac{N}{2})\right]^{2}},$$

(33)

where

$$\begin{split}&F_{11}^{-}(q)=\frac{|\langle\varphi_{1}^{q}|\chi_{1}\rangle|}{\sqrt{\langle\varphi_{1}^{q}|\varphi_{1}^{q}\rangle\langle\chi_{1}|\chi_{1}\rangle}}=\frac{\sqrt{Re_{-}^{2}+Im^{2}}}{\sqrt{\sum_{k=0}^{\infty}\frac{{(2\mu_{x})}^{kN+1}}{{(kN+1)!}}\sum_{k=0}^{\infty}\frac{\mu_{z}^{kN+1}}{(kN+1)!}}},\\ &Re_{-}=\frac{1}{\sqrt{2}}\sum_{k=0}^{\infty}\frac{(\mu_{x}\mu_{z})^{(kN+1)/2}}{(kN+1)!}[1-\cos\frac{2\pi}{N}q(kN+1)].\end{split}$$

(34)

and $Im$ has already been shown in Eq. (31). Finally, we have

$$e^{ph}\leq\frac{T_{0}^{R,U}+T_{1}^{L,U}}{s_{01}^{L}+s_{10}^{L}}.$$

(35)

With all those formulas, we can now calculate the secure final key rate according to the observed values in the experiment.

We have obtained explicit formulas above for key rate calculation. In our calculation below, we shall use our analytical formulas above. Definitely, the key rate here can be also calculated through linear programming. It is

$$\begin{split}\min\quad&s_{1}=\frac{1}{2}(s_{01}+s_{10}),\\ s.t.\;\;\;\,\,&{\rm constraints\;with\;observed\;values\;of\;number\;of\;post\;selected\;counts,\;e.g.}\\ \quad&S_{ox}=\sum_{j=0}^{N-1}P_{j}(\mu_{x})Y_{vj}^{x},\quad S_{oy}=\sum_{j=0}^{N-1}P_{j}(\mu_{y})Y_{vj}^{y},\\ &S_{xo}=\sum_{j=0}^{N-1}P_{j}(\mu_{x})Y_{jv}^{x},\quad S_{yo}=\sum_{j=0}^{N-1}P_{j}(\mu_{y})Y_{jv}^{y},\\ &|S_{oo}-Y_{0v}^{y}|\leq\sqrt{1-F_{0}^{2}},\quad|S_{oo}-Y_{v0}^{y}|\leq\sqrt{1-F_{0}^{2}},\\ &|Y_{vj}^{x}-Y_{vj}^{y}|\leq\sqrt{1-(F_{xy}^{j})^{2}},\quad|Y_{jv}^{x}-Y_{jv}^{y}|\leq\sqrt{1-(F_{xy}^{j})^{2}},\\ &|s_{01}-Y_{v1}^{y}|\leq\sqrt{1-F_{1}^{2}},\quad|s_{10}-Y_{1v}^{y}|\leq\sqrt{1-F_{1}^{2}},\end{split}$$

(36)

and so on. After we get the lower bound of $s_{1}$ through linear programming, the remaining process is the same to the above method.

The 3-intensity protocol is simply a special case of our 4-intensity above by setting $\mu_{z}=\mu_{y}$.