Proposal for distribution of multi-photon entanglement with optimal rate-distance scaling

TL;DR

This paper introduces a protocol for distributing multi-photon entanglement over long distances using current integrated quantum photonic technology, enabling secure quantum communication and metrology with optimal rate-distance scaling.

Contribution

It presents a novel scheme utilizing squeezed vacuum states and photon-number-resolving detectors for long-range entanglement distribution with high loss tolerance and optimal rate scaling.

Findings

Achieves near-maximal entanglement distribution over long distances.

Supports loophole-free Bell tests for certification.

Provides measurement-device-independent security with potential for full device independence.

Abstract

We propose a protocol to perform long-range distribution of near-maximally entangled multiphoton states, allowing versatile applications such as quantum key distribution (QKD) and quantum metrology which can provide alternatives to state-of-the-art protocols. Our scheme uses resources available within the current integrated quantum photonic technology: squeezed vacuum states and photon-number-resolving detectors. The distributed entanglement can be certified by Bell tests which have the potential to be loophole free, and may be directly used in well established QKD protocols. Generally, this provides measurement-device-independent (MDI) levels of security, which may be upgraded to fully device-independent (DI) security if the Bell test is loophole-free. In both cases, the protocol is robust to extremely high transmission losses, matching the optimal $O(\sqrt{\eta})$ scaling of key rate with channel transmittance.

Tables

Table 1. Table 1: Example states | Ψ out ( k , σ ) ⟩ ket superscript subscript Ψ out 𝑘 𝜎 \ket{\Psi_{\text{out}}^{(k,\sigma)}} for σ ≤ 3 𝜎 3 \sigma\leq 3 .

$k$	$\sigma$	$|{\Psi }_{\text{out}}^{(k,\sigma )}⟩$
0	1	$\frac{1}{\sqrt{2}}\left(|0,1⟩-i|1,0⟩\right)$
1	1	$\frac{1}{\sqrt{2}}\left(|0,1⟩+i|1,0⟩\right)$
0	2	$\frac{1}{2}|0,2⟩-\frac{i}{\sqrt{2}}|1,1⟩-\frac{1}{2}|2,0⟩$
1	2	$\frac{1}{\sqrt{2}}\left(|0,2⟩+|2,0⟩\right)$
2	2	$\frac{1}{2}|0,2⟩+\frac{i}{\sqrt{2}}|1,1⟩-\frac{1}{2}|2,0⟩$
0	3	$\frac{1}{\sqrt{8}}\left(|0,3⟩-i\sqrt{3}|1,2⟩-\sqrt{3}|2,1⟩+i|3,0⟩\right)$
1	3	$\frac{1}{\sqrt{8}}\left(\sqrt{3}|0,3⟩-i|1,2⟩+|2,1⟩-i\sqrt{3}|3,0⟩\right)$
2	3	$\frac{1}{\sqrt{8}}\left(\sqrt{3}|0,3⟩+i|1,2⟩+|2,1⟩+i\sqrt{3}|3,0⟩\right)$
3	3	$\frac{1}{\sqrt{8}}\left(|0,3⟩+i\sqrt{3}|1,2⟩-\sqrt{3}|2,1⟩-i|3,0⟩\right)$

Table 2. Table 2: Optimal settings for the Bell test for each state | Ψ out ( k , σ ) ⟩ ket subscript superscript Ψ 𝑘 𝜎 out \ket{\Psi^{(k,\sigma)}_{\textrm{out}}} . | B | 𝐵 |B| is the value of the Bell expression for zero idler losses. | B | 𝐵 |B| (40dB) is the value for 40dB idler losses (80dB losses in the Alice-Bob channel). A value of | B | > 2 𝐵 2 |B|>2 shows non-locality, with a blank entry signifying that there is no observed non-locality ( | B | ≤ 2 𝐵 2 |B|\leq 2 ). The detector efficiency shown is that required when there are 40dB idler losses.

$k$	$\sigma$	variables	$[{\delta }_{\alpha 1},{\delta }_{\alpha 2}],[{\delta }_{\beta 1},{\delta }_{\beta 2}]$	$\left|B\right|$	$\left|B\right|$ (40dB)	Alice/Bob detector efficiency (40dB)
0	1	zero/non-zero	[0.17, -0.56], [0.17$i$, -0.56$i$]	2.69	2.63	$\ge 85$%
1	1	zero/non-zero	[0.17, -0.56], [-0.17$i$, 0.56$i$]	2.69	2.63	$\ge 85$%
0	2	even/odd	[0.06, -0.28], [0.06$i$, -0.28$i$]	2.35	2.17	$\ge 98$%
1	2	zero/non-zero	[0.00, 0.63], [0.00, 0.63]	2.07	2.03	$\ge 96$%
2	2	even/odd	[0.06, -0.28], [-0.06$i$, 0.28$i$]	2.35	2.17	$\ge 98$%
0	3	even/odd	[0.06, -0.24], [0.06$i$, -0.24$i$]	2.40	2.16	$\ge 99$%
1	3	even/odd	[0.02, -0.22], [0.02$i$, -0.22$i$]	2.10	—	—
2	3	even/odd	[0.02, -0.22], [-0.02$i$, 0.22$i$]	2.10	—	—
3	3	even/odd	[0.06, -0.24], [-0.06$i$, 0.24$i$]	2.40	2.16	$\ge 99$%

Table 3. Table 3: Optimal settings for the QKD protocol in the lossless case.

$k$	$\sigma$	$[{\delta }_{\alpha 0},{\delta }_{\alpha 1},{\delta }_{\alpha 2}],[{\delta }_{\beta 1},{\delta }_{\beta 2}]$	$\left|B\right|$	Bit error rate, $q$	Raw key rate, $K$	Key rate per pulse, $R$
0	1	[-0.17, 0.17, -0.56], [0.17$i$, -0.56$i$]	2.688	$4\times {10}^{-4}$	$\ge 0.556$	$\ge 1.39\times {10}^{-2}$
1	1	[-0.17, 0.17, -0.56], [-0.17$i$, 0.56$i$]	2.688	$4\times {10}^{-4}$	$\ge 0.556$	$\ge 1.39\times {10}^{-2}$

Table 4. Table 4: Optimal settings for the QKD protocol in the case of 40dB idler losses (80dB losses in the Alice-Bob channel).

$k$	$\sigma$	$[{\delta }_{\alpha 0},{\delta }_{\alpha 1},{\delta }_{\alpha 2}],[{\delta }_{\beta 1},{\delta }_{\beta 2}]$	$\left|B\right|$	Bit error rate, $q$	Raw key rate, $K$	Key rate per pulse, $R$
0	1	[-0.17, 0.17, -0.56], [0.17$i$, -0.56$i$]	2.635	$0.01$	$\ge 0.705$	$\ge 1.12\times {10}^{-6}$
1	1	[-0.17, 0.17, -0.56], [-0.17$i$, 0.56$i$]	2.635	$0.01$	$\ge 0.705$	$\ge 1.12\times {10}^{-6}$

Equations

$$|\Psi ⟩=\sum _{n=0}^{\infty }\sqrt{{\lambda }_n}{|n⟩}_1{|n⟩}_2,$$

$${|n⟩}_{a_2}=\frac{(a_2^{†}{)}^n}{\sqrt{n!}}|0⟩,{|\sigma -n⟩}_{b_2}=\frac{(b_2^{†}{)}^{\sigma -n}}{\sqrt{(\sigma -n)!}}|0⟩,$$

$${\mathcal{A}}_{\sigma }(k,n)=$$

$$=$$

$$|{\Psi }_{\text{out}}^{(k,\sigma )}⟩=$$

$$=$$

$$\times ⟨k,\sigma -k|U_{\text{BS}}{|n,m⟩}_{a_2,b_2},$$

$$|{\Psi }_{\text{out}}^{(k,\sigma )}⟩=\sum _{n=0}^{\sigma }{\mathcal{A}}_{\sigma }(k,n){|n,\sigma -n⟩}_{a_1,b_1}.$$

$$\displaystyle E_{\mathcal{N}}\big{(}\ket{\Psi_{\text{out}}^{(k,\sigma)}}\bra{\Psi_{\text{out}}^{(k,\sigma)}}\big{)}={}$$

$$\displaystyle\qquad{}=2\log_{2}\Bigg{\{}\sum_{n=0}^{\sigma}\Big{|}\phi_{k}(n-\tfrac{\sigma}{2},\sigma)\Big{|}\Bigg{\}},$$

$${\rho }_{\text{out}}^{(k,\sigma )}=\sum _{S=\sigma }^{\infty }{p}_{S|\sigma }{\rho }_{\text{int}}^{(\sigma ,k,S)},$$

$$p_{(k,\sigma )}=\frac{{\lambda }_{\sigma }}{{\cosh }^{2}g}\cdot \frac{(1-r{)}^{\sigma }}{(1-r{\tanh }^{2}g{)}^{\sigma +2}},$$

$$|{\Psi }_{\text{out}}^{(k,\sigma )}⟩=\sum _{n=0}^{\sigma }{\mathcal{A}}_{\sigma }(k,n){|n,\sigma -n⟩}_{a_1,b_1},$$

$${\mathcal{A}}_{\sigma }(k,n)=⟨k,\sigma -k|U_{\text{BS}}|n,\sigma -n⟩=i^{k-n}{\varphi }_k(n-\sigma /2,\sigma ).$$

$${\varphi }_k(n-\sigma /2,\sigma )=(-1{)}^k\sqrt{\frac{1}{2^{\sigma }}\left(\genfrac{}{}{0ex}{}{\sigma }{n}\right)\left(\genfrac{}{}{0ex}{}{\sigma }{k}\right)}{}_2{F}_1(-k,-n;-\sigma ;2),$$

$$\sum _{n=0}^{\sigma }{\varphi }_k(n-\sigma /2,\sigma ){\varphi }_{k^{\prime }}(n-\sigma /2,\sigma )={\delta }_{k,k^{\prime }}.$$

$$⟨{\Psi }_{\text{out}}^{(k,\sigma )}⟩{\Psi }_{\text{out}}^{(k^{\prime },\sigma )}={\delta }_{k,k^{\prime }}.$$

$$E_{\mathcal{N}}(\rho_{\text{out}}^{(k,\sigma)})=2\log_{2}\left(\sum_{n=0}^{\sigma}\bigl{\lvert}\mathcal{A}_{\sigma}(k,n)\bigr{\rvert}\right).$$

$${\rho }_{\text{out}}^{(k,\sigma )}$$

$$=\sum _{n,m=0}^{\sigma }\begin{array}{rl}& \underset{{\mathcal{A}}_{\sigma }(k,n){\mathcal{A}}_{\sigma }^{\ast }(k,m)}{\underbrace{[\rho {]}_{n,\sigma -n;m,\sigma -m}}}\times \\ & \times |n,\sigma -n⟩⟨m,\sigma -m|,\end{array}$$

$${\left({\rho }_{\text{out}}^{(k,\sigma )}\right)}^{{\Gamma }_2}$$

$${\left({\left({\rho }_{\text{out}}^{(k,\sigma )}\right)}^{{\Gamma }_2}\right)}^{†}$$

$${\left({\left({\rho }_{\text{out}}^{(k,\sigma )}\right)}^{{\Gamma }_2}\right)}^{†}\cdot {\left({\rho }_{\text{out}}^{(k,\sigma )}\right)}^{{\Gamma }_2}$$

$$=\sum _{n,m,p,q=0}^{\sigma }[\rho {]}_{n,\sigma -n;m,\sigma -m}^{\ast }[\rho {]}_{p,\sigma -p;q,\sigma -q}{\delta }_{n,p}{\delta }_{m,q}|m,\sigma -n⟩⟨q,\sigma -p|$$

$$=\sum _{n,m=0}^{\sigma }[\rho {]}_{n,\sigma -n;m,\sigma -m}^{\ast }[\rho {]}_{n,\sigma -n;m,\sigma -m}|m,\sigma -n⟩⟨m,\sigma -n|$$

$$\displaystyle{}=\sum_{n,m=0}^{\sigma}\bigl{\lvert}[\rho]_{n,\sigma-n;m,\sigma-m}\bigr{\rvert}^{2}\ket{m}\bra{m}\otimes\ket{\sigma-n}\bra{\sigma-n},$$

$${\left({\left({\rho }_{\text{out}}^{(k,\sigma )}\right)}^{{\Gamma }_2}\right)}^{†}\cdot {\left({\rho }_{\text{out}}^{(k,\sigma )}\right)}^{{\Gamma }_2}$$

$$\displaystyle{}=\left(\sum_{m=0}^{\sigma}\bigl{\lvert}\mathcal{A}_{\sigma}(k,m)\bigr{\rvert}^{2}\,\ket{m}\bra{m}\right)\otimes\left(\sum_{n=0}^{\sigma}\bigl{\lvert}\mathcal{A}_{\sigma}(k,n)\bigr{\rvert}^{2}\,\ket{\sigma-n}\bra{\sigma-n}\right).$$

$$\sqrt{{\left({\left({\rho }_{\text{out}}^{(k,\sigma )}\right)}^{{\Gamma }_2}\right)}^{†}\cdot {\left({\rho }_{\text{out}}^{(k,\sigma )}\right)}^{{\Gamma }_2}}$$

$${\Vert {\left({\rho }_{\text{out}}^{(k,\sigma )}\right)}^{{\Gamma }_2}\Vert }_1=\mathrm{Tr}\sqrt{{\left({\left({\rho }_{\text{out}}^{(k,\sigma )}\right)}^{{\Gamma }_2}\right)}^{†}\cdot {\left({\rho }_{\text{out}}^{(k,\sigma )}\right)}^{{\Gamma }_2}}$$

$$\displaystyle{}=\left(\sum_{n=0}^{\sigma}\bigl{\lvert}\mathcal{A}_{\sigma}(k,n)\bigr{\rvert}\right)^{2}.$$

$$|\Psi ⟩=\sum _{n=0}^{\infty }\sqrt{{\lambda }_n}{|n⟩}_1{|n⟩}_2,$$

$${\rho }_{\text{in}}=$$

$$\rho^{(k,\sigma)}_{\text{out}}=\mathcal{N}^{2}\!\!\!\sum_{n,n^{\prime},m,m^{\prime}=0}^{\infty}\!\!\!\begin{aligned} &\delta_{\sigma,n+m}\,\delta_{\sigma,n^{\prime}+m^{\prime}}\times{}\\ &\qquad\times\sqrt{\lambda_{n}\,\lambda_{n^{\prime}}\,\lambda_{m}\,\lambda_{m^{\prime}}}\times{}\\ &\qquad\times\mathcal{A}_{\sigma}(k,n)\,\bigl{(}\mathcal{A}_{\sigma}(k,n^{\prime})\bigr{)}^{*}\times{}\\ &\qquad\times\ket{n,m}\bra{n^{\prime},m^{\prime}}_{a_{1},b_{1}},\end{aligned}$$

$${\rho }_{\text{out}}^{(k,\sigma )}$$

$$=|{\Psi }_{\text{out}}^{(k,\sigma )}⟩⟨{\Psi }_{\text{out}}^{(k,\sigma )}|.$$

$${\mathrm{Tr}}_r\left\{U_{\text{BS}}^{(r)}|n,0⟩⟨n,0|{\left(U_{\text{BS}}^{(r)}\right)}^{†}\right\}$$

$$=\sum _{p=0}^n\left(\genfrac{}{}{0ex}{}{n}{p}\right)(1-r{)}^{n-p}{r}^p|n-p⟩{⟨n-p|}_t,$$

Full text

I Introduction

Distribution of photonic entanglement is key to building quantum networks which facilitate secure long-distance quantum communications [1], distributed quantum computation [2] and sensing [3, 4], leading finally to a full realization of the quantum internet [5]. A fundamental problem is that entanglement becomes easily corrupted by losses in the transmitting channels, which results in low transmission rates of entangled states. Since amplification of quantum signals is impossible [6], alternative remedies are in high demand. A popular option is to use quantum repeaters [7, 8]. However, this approach necessitates numerous intermediate stations, quantum memories and multiple two-photon Bell pairs, a resource that is often created non-deterministically. Despite rapid advances in the field [9, 10, 11, 12, 13, 14, 15, 16, 17, 18], experimental realizations of quantum repeaters currently show deterioration with distance which is insufficient for most quantum communications applications [19]. Another possibility is to use satellites to distribute two-photon Bell pairs, which is advantageous thanks to photons traveling the majority of the distance in free space and thus experiencing less loss than in atmospheric links or optical fibers, while connecting arbitrarily distant locations. This approach has seen much progress recently [31, 32, 33, 24, 25, 26, 34, 27, 28, 29, 30, 19, 20, 21, 22, 23], although the losses are still detrimental, with typically only approximately 1 in $10^{6}$ of produced Bell pairs surviving the trip to the detection stations on the ground [34]. Thus the distribution of entanglement in a way that is resource-efficient, verifiable and well suited to the existing quantum-photonic technologies continues to be an open problem.

Quantum key distribution (QKD) uses quantum correlations to distribute cryptographic keys between remote parties which enables secure communication. In the long distance regime, the current state-of-the-art is twin-field (TF)-QKD [35], which displays a ground breaking improvement over previous methods, exemplified by the key-rate scaling $O(\sqrt{\eta})$ with channel transmittance $\eta$, allowing efficient demonstrations of QKD over 600km [36, 37, 38]. There, weak coherent pulses are phase-encoded and sent to an untrusted node for measurement. This is an example of measurement-device-independent (MDI)-QKD [39], which closes all detector side channels but remains susceptible to eavesdropping attacks at the source. For example, photon-number-splitting attacks [40] which must be carefully avoided via decoy-state methods [41]. Entanglement-based versions, which benefit from the same key-rate scaling, were soon theorized [42, 43]. Generally, entanglement-based approaches make security easier to verify and have the potential for loophole-free Bell tests which could close all possible side channels. This is called device-independent (DI)-QKD [44, 23]. It is clear that for entanglement distribution protocols to remain competitive for QKD at long distances, they must achieve the same optimal $O(\sqrt{\eta})$ key-rate scaling.

Here, we propose a protocol that uses multiphoton bipartite entanglement and photon-number-resolving (PNR) detectors to establish long-distance near-maximal entanglement distribution, allowing entanglement-based QKD protocols that match the $O(\sqrt{\eta})$ scaling. It is robust to high transmission losses which, remarkably, deteriorate only its efficiency, but not the amount of generated entanglement, a feature until now seen only in protocols based on Bell pairs. As a result, the protocol has the potential to outperform other methods of entanglement distribution which have transmission rates scaling as $O(\eta)$. This is particularly true in high loss scenarios, for example satellite-based Earth-space channels and long fiber networks. Once distributed, the multiphoton entanglement of the states may be verified by a CHSH Bell inequality [45], forming a basis for establishing a cryptographic key. Depending on the detection efficiencies available, this can be achieved either by DI-QKD [42], or entanglement-based phase-matching QKD [43]. The higher photon number states generated are additionally useful for quantum metrology which has been extensively covered in our previous work [46].

Bipartite entanglement between two physical systems such as, e.g. two photons, can be revealed by their Bell nonlocality, i.e. the fact that outcomes of local measurements on these two subsystems exhibit correlations which cannot be described by local hidden variable models [47]. In this paper we focus on multiphoton bipartite entanglement, engineered from two-mode squeezed-vacuum (SV) states, which are used as our primary resource. This Gaussian quantum state is naturally and deterministically produced from a pulsed coherent light via the spontaneous parametric down-conversion process (SPDC), the most ubiquitous source of quantum light [48]. It carries entanglement between the quadratures of the electric field (continuous variables, CVs) but also between photon numbers in the two modes (discrete variables, DVs), which are unbounded [49]. The latest advances in experimental integrated quantum photonics facilitate precise generation and manipulation of near-single-mode SV states in optical chips [50], as well as detection of their photon number statistics using state-of-the-art PNR detectors [51]. The most successful PNR detectors are transition edge sensors (TES) which have negligible dark count rates and quantum efficiencies which can exceed $95\%$ [52, 53]. Integrated photonic devices have also contributed to the development of practical quantum communications scenarios [54, 55].

In our protocol, a verifiable DV entanglement of high local dimension is created from two sources of an SV state using an entangling measurement at a remote station. The measurement is realized by multiphoton quantum interference on a beam splitter (BS) followed by PNR detectors. The protocol can also be carried out in a delayed-choice scheme which frees the parties from using quantum memories and allows them to share near-maximally entangled states in realistic implementations. These states have been proven to be useful not only for QKD, as specified above, but also find an immediate application in quantum metrology. These are multiphoton generalized Holland–Burnett (HB) states, which have been experimentally proven to allow near-optimal quantum-enhanced optical phase estimation in a lossy environment [46].

II Results

II.1 Resources

We employ two copies of an SV state $\ket{\Psi}$ coming from a pulsed source as our input, $\ket{\Psi_{\text{in}}}=\ket{\Psi}_{a}\otimes\ket{\Psi}_{b}$. In its Schmidt basis, $\ket{\Psi}$ takes the form of a superposition of $n$-photon pairs with real-valued probability amplitudes $\sqrt{\lambda_{n}}=\frac{\tanh^{n}g}{\cosh g}$

[TABLE]

where $g$ is the parametric gain, which is the key parameter characterizing an SPDC source, setting the mean photon number in $\ket{\Psi}$ to $2\sinh^{2}g$ for each pulse. Quantum correlations in $\ket{\Psi}$ are manifested by equal photon numbers in modes 1 and 2, called the signal and idler, which can be spatially resolved. The subsequent photon number contributions to the SV state: the vacuum, single-photon and higher-order ($n>1$) emissions, occur with a probability which follows a geometric progression with common ratio $\tfrac{\lambda_{n+1}}{\lambda_{n}}=\tanh^{2}g$. Thus, the SPDC generates a considerable amount of multiphoton events which are utilized to our benefit in our protocol.

A schematic explanation of the proposed entanglement distribution protocol is shown in Fig. 1. The setup consists of two pulsed SPDC sources, one held by Alice and one by Bob, each generating $\ket{\Psi}$. The idler beams emitted into modes $a_{2}$ and $b_{2}$ are sent to Charlie at a remote station, who interferes them on a balanced ($50:50$) BS and then detects them with PNR detectors. This is an entangling measurement, whereby the signal modes $a_{1}$ and $b_{1}$ become photon-number correlated.

II.2 Protocol in ideal circumstances

Let us consider the lossless case first. The detection of $\sigma$ photons in total in the entangling measurement means that $\sigma$ photons were distributed between the two idler beams entering the BS and thus that the total number of photons, $S$, in the output state shared by Alice and Bob is $S=\sigma$. The BS performs a linear operation on the input idler creation operators, $U_{\text{BS}}^{\dagger}\,a_{2}^{\dagger}\,U_{\text{BS}}=\tfrac{1}{\sqrt{2}}\,(a_{2}^{\dagger}-ib_{2}^{\dagger})$, $U_{\text{BS}}^{\dagger}\,b_{2}^{\dagger}\,U_{\text{BS}}=\tfrac{1}{\sqrt{2}}\,(-ia_{2}^{\dagger}+b_{2}^{\dagger})$, while the signal modes $a_{1}$ and $b_{1}$ are intact. Applying this operation to Fock states

[TABLE]

requires taking powers of the transformed operators. This results in a transformation on $a_{2}^{\dagger}$ and $b_{2}^{\dagger}$ governed by a binomial distribution and an output state described by an arcsine probability distribution. Through such a multiphoton Hong–Ou–Mandel (HOM) effect, entanglement between the BS output modes is generated [56]. The probability amplitudes of detecting $k$ and $\sigma\!-\!k$ photons behind the BS are equal to [57]

[TABLE]

where $k=0,\dots,\sigma$ and $\phi_{k}$ are symmetric Kravchuk functions – orthonormal discrete polynomials which converge to Hermite–Gauss polynomials for large $\sigma$ [58]. The output state in our protocol is therefore

[TABLE]

where $\mathcal{N}$ is the normalization and $\ket{\Psi}$ is defined in Eq. (1). Since $n+m=\sigma$ must hold true,

[TABLE]

We now see that $\ket{\Psi_{\text{out}}^{(k,\sigma)}}$ lives in a Hilbert space of dimension $(\sigma\!+\!1)$ and is parameterized by the entangling measurement readouts $k$ and $\sigma\!-\!k$ which define its photon number statistics. The probability of detecting $n$ photons in mode $a_{1}$ (and $\sigma-n$ in mode $b_{1}$) is $\lvert\mathcal{A}_{\sigma}(k,n)\rvert^{2}=\lvert\phi_{k}(n-\frac{\sigma}{2},\sigma)\rvert^{2}$. Example photon number distributions, and details of the derivation, are given in Figs. 6-7 and in Appendix A.

This protocol can also be realized in a delayed-choice scheme. Then, the measurements taken by Alice and Bob on the signal modes $a_{1}$ and $b_{1}$ precede the ones performed on the idler modes $a_{2}$ and $b_{2}$ at the remote station. By the no-signaling principle, the photon-number statistics observed by Alice and Bob is again determined by $\lvert\mathcal{A}_{\sigma}(k,n)\rvert^{2}$.

To quantify entanglement in $\ket{\Psi_{\text{out}}^{(k,\sigma)}}$ we employ the logarithmic negativity $E_{\mathcal{N}}=\log_{2}\lvert\lvert\rho^{\Gamma}\rvert\rvert_{1}$, where $\rho$ denotes a density operator, $\Gamma$ is the partial transpose operation and $\lvert\lvert\cdot\rvert\rvert_{1}$ is the trace norm. This entanglement measure is easily computable and gives an upper bound to the distillable entanglement, reflecting the maximal number of Bell states that can be extracted [59, 60]. Inserting $\rho=\ket{\Psi_{\text{out}}^{(k,\sigma)}}\bra{\Psi_{\text{out}}^{(k,\sigma)}}$ we obtain

[TABLE]

see Appendix B for the full derivation. Since the readouts $k$ and $\sigma\!-\!k$ uniquely define the state $\ket{\Psi_{\text{out}}^{(k,\sigma)}}$, they also determine the amount of entanglement in it. The maximum amount of entanglement which can be created in Hilbert space of dimension $(\sigma+1)$ is $E_{\mathcal{N}_{\text{max}}}=\log_{2}(\sigma+1)$. Our protocol allows one to achieve values close to this maximum, as shown in the lossless case for $\sigma=4$ in Fig. 2, with the maximal value $E_{\mathcal{N}_{\text{max}}}\approx 2.3$. The decrease in entanglement for $k=\frac{\sigma}{2}$ results from the multiphoton HOM effect, whereupon all odd photon-number components disappear. This leads to an effective reduction in the Hilbert space dimension to $(\frac{\sigma}{2}+1)$ where the maximal entanglement is $\log_{2}(\frac{\sigma}{2}+1)$, or $\approx 1.6$ at $\sigma=4$. For comparison, Bell states have logarithmic negativity, $E_{\mathcal{N}_{\text{Bell}}}=1$. This is the same as we obtain for $\sigma=1$.

II.3 Operation in lossy conditions

Depending on the application of our protocol, the amount of losses in the photon transmission may vary, ranging from 0.2dB/km for optical fibers to approximately 40dB for an uplink ground-to-satellite channel [61]. In addition, the efficiency of a PNR detection system may drop to 50–60% (by 2–3dB) if one includes uncorrelated background counts and inefficient coupling [46].

Losses are modeled by inserting additional beam splitters into the pathways of the photons, whose reflectivity $r$ quantifies the amount of loss. In the case of symmetric idler mode losses, the state shared by Alice and Bob is no longer the pure state $\ket{\Psi_{\text{out}}^{(k,\sigma)}}$, but a mixed state described by the density operator

[TABLE]

where $\rho^{(\sigma,k,S)}_{\text{int}}$ is a density operator of a state with $S$ photons in total in the signal modes which are assumed lossless. The full representation can be found in Appendix C. Unlike in the lossless case, $S$ is not always equal to $\sigma$, and instead the probability for Alice and Bob to detect $S$ photons between them is given by $p_{S|\sigma}=\tilde{\mathcal{N}}^{2}\chi_{\sigma,S}$ where $\chi_{\sigma,S}=r^{S-\sigma}\lambda_{S}\,\binom{S+1}{\sigma+1}$ and $\tilde{\mathcal{N}}$ is a normalization constant independent of $S$. As shown in Fig. 3, in the limit $r\tanh^{2}(g)\ll 1$, the most likely event is $S=\sigma$, with additional components $S>\sigma$ becoming increasingly less relevant as $\chi_{\sigma,S}$ drops rapidly towards zero. The primary component of the mixed state is then $\rho^{(\sigma,k,\sigma)}_{\text{int}}$, which is identical to the lossless state $\ket{\Psi_{\text{out}}^{(k,\sigma)}}\bra{\Psi_{\text{out}}^{(k,\sigma)}}$. Provided that the gain $g$ is at typical low levels, $g^{2}\ll 1$, the limit $r\tanh^{2}(g)\ll 1$ holds even as $r\to 1$, so that the entanglement of the final states remains high even in the presence of arbitrarily high losses in the idler modes.

To understand how realistic conditions affect the distributed entanglement we consider different scenarios of losses. Fig. 2 shows numerically calculated logarithmic negativity at $g=0.1$ in the ideal case, the case of large symmetric losses in the idler beams ($r_{a_{2},b_{2}}\!=\!99.99\%$, i.e., 40dB), the case of unbalanced losses in the idler beams with relative transmittance $\epsilon=\frac{1-r_{b_{2}}}{1-r_{a_{2}}}=0.5$ ($r_{a_{2}}\!=\!99.98\%$, $r_{b_{2}}\!=\!99.99\%$, i.e., 37, 40dB, respectively) and the case of non-ideal signal modes with 20% losses in those arms. Losses in dB are given by $-10\log_{10}(1-r)$. We do not assume perfect entangling measurements, instead setting realistic losses at the detectors ($r_{d_{1},d_{2}}\!=\!50\%$), however, we found that their influence is negligible since they can be viewed as a minor addition to the already large idler losses. For more examples of scenarios involving losses see Figs. 8-10 and 12-15, and Appendices C-E. As can be seen in Fig. 2, the amount of entanglement obtained in the case of large symmetric losses is almost exactly the same as in the ideal case. In fact, the entangling measurement can be viewed as a $g^{(2)}$ correlation function measurement, which is known to be immune to losses [62].

As for the asymmetrical transmission losses, the logarithmic negativity stays within 90% of the lossless values provided that the relative transmittance of the channels $\epsilon$ stays above 0.4, see Fig. 13 and Appendix D for an explanation of this effect. One should note that relative transmittance $\epsilon<0.4$ is unlikely to occur in reality [29, 61] and artificial losses may always be introduced to balance the channels. Another consideration is that losses for each channel will vary with time [61], especially with atmospheric losses when applied to a satellite-based scheme. However, losses fluctuate on time-scales an order of magnitude higher than the time of flight of photons between the ground and satellite, therefore it suffices to take the average of Alice and Bob’s output logarithmic negativity. The effect of fluctuations on the average is not detrimental, as shown in Fig. 16.

In contrast to idler mode losses, detection inefficiencies or losses in signal modes $a_{1}$ and $b_{1}$ are critical since they spoil generated entanglement, see Fig. 10 for more examples of how signal losses affect the entanglement. However, these losses may be significantly reduced by using a delayed-choice scheme. In such a scheme, Alice and Bob would perform their measurements immediately before waiting for the detection by Charlie at the remote station. By the no-signaling principle, all measurement outcomes are independent of this modification. If Alice, Bob and Charlie’s measurements have space-like separation from each other, this can allow Bell tests which simultaneously close the locality and detection loopholes, as in “event-ready” Bell tests [63]. In our case the event-ready signal is Charlie’s successful measurement. The delayed-choice scheme also mitigates the effect of phase-damping. After Alice and Bob perform their near-immediate measurement of the signal beams, from their reference frame, the idlers are projected onto Fock states which are unaffected by phase-damping. Remaining synchronization errors and phase-damping may be further reduced by using two reference pulses, one to measure the delay in the arrival times of the two idler beams, and one to measure the rate of phase-shifting. These errors can then be corrected in postprocessing [36, 37, 38].

II.4 Protocol performance

The probability to generate the state $\rho^{(k,\sigma)}_{\text{out}}$ is equal to the probability for Charlie to obtain the measurement outcomes $k$ and $\sigma-k$ photons after the BS

[TABLE]

which is independent of $k$, see Appendix F for the calculation. This probability is useful to calculate the efficiency of our protocol, which we define as the probability of Charlie receiving at least one photon. When that happens, we are certain that Alice and Bob share entangled states, which may be multiphoton. The protocol’s efficiency, $1-p_{(0,0)}\approx 2(1-r)\lambda_{1}/\cosh^{2}(g)$, increases with the parametric gain $g$ and decreases with idler losses $r$. Investigating the dependence on transmission losses more closely, we see that it is proportional to the transmittance $t=1-r$ between Alice/Bob and Charlie. It is customary to instead use the transmittance $\eta=t^{2}$ of the summed distance from Alice to Bob [35], assuming that Charlie is approximately at the mid-point between them. In this context, it is clear that our protocol’s efficiency scales as $O(\sqrt{\eta})$. This should be contrasted to schemes where Charlie directly distributes Bell pairs to Alice and Bob [34, 29], which have efficiencies $O(\eta)$. The reason for the quadratic improvement is that only a single photon has to survive the trip between Alice/Bob and Charlie, since an entangled state is established even with $\sigma=1$. In contrast, distributing a Bell pair requires two photons to survive the trip, therefore losses are more significant.

Fig. 4 depicts the protocol efficiency and generated entanglement as a function of the idler losses for three values of the parametric gain: $g=\{0.1,0.2,0.4\}$. Provided that $g^{2}\ll 1$, as can be seen for the case $g=0.1$, the idler losses influence only the efficiency of the protocol without changing its principle of operation or the generated entanglement. This feature is characteristic to protocols based on two-photon Bell pair entanglement, and now we have shown it for all multiphoton schemes based on our setup. If $g$ is increased, the efficiency of the protocol is improved but at the cost of reducing the generated entanglement, as can be seen for $g=0.2$ and $g=0.4$.

II.5 Applications

II.5.1 Quantum key distribution

Depending on the detection efficiencies available, two different Bell tests are possible for Alice and Bob to verify the amount of entanglement in their shared state, with each test being associated with a well-established QKD protocol: single-photon entanglement-based DI-QKD [64, 42], and single-photon entanglement-based phase-matching (SEPM)-QKD [43]. Both of these protocols display the optimal $O(\sqrt{\eta})$ scaling of key rate with channel transmittance $\eta$. They rely on Alice and Bob interfering their signal modes with coherent states $\ket{\alpha}$ and $\ket{\beta}$ on local beam splitters with reflectivities $r_{a}$ and $r_{b}$ respectively, and then measuring the output photon numbers.

The first test is specifically suited to cases where Alice and Bob’s detector efficiencies and signal mode transmittances are high enough for the test to be performed without any postselection. Alice and Bob interfere their signal modes with coherent states to perform the displacement operations $D(-\delta_{\alpha})=\exp(\delta_{\alpha}^{*}a_{1}-\delta_{\alpha}a_{1}^{\dagger})$ and $D(-\delta_{\beta})=\exp(\delta_{\beta}^{*}b_{1}-\delta_{\beta}b_{1}^{\dagger})$, where $\delta_{\alpha}=i\alpha\sqrt{r_{a}}$ and $\delta_{\beta}=i\beta\sqrt{r_{b}}$, and finally measure their transmitted photon numbers with PNR detectors. They evaluate a CHSH inequality by assigning ${+}1/{-}1$ either to vacuum/non-vacuum events, or to even/odd photon numbers. By choosing randomly between two local measurement settings each, $\delta_{\alpha 1},\delta_{\alpha 2}$ and $\delta_{\beta 1},\delta_{\beta 2}$, they can then observe violation of the CHSH inequality $|B|\leq 2$ where $B=E(\delta_{\alpha 1},\delta_{\beta 1})+E(\delta_{\alpha 1},\delta_{\beta 2})+E(\delta_{\alpha 2},\delta_{\beta 1})-E(\delta_{\alpha 2},\delta_{\beta 2})$, and $E(\delta_{\alpha},\delta_{\beta})$ is the correlation between their dichotomized variables. The greatest violation $|B|=2.63$ is found for the $\sigma=1$ states, even with large transmission losses in excess of 40dB, where they comprise the vast majority (99.9999%) of the generated states. For these states we found the optimal choice of measurement settings to be $\delta_{\alpha 1}=+0.17,\delta_{\alpha 2}=-0.56,\delta_{\beta 1}=\pm 0.17i,\delta_{\beta 2}=\mp 0.56i$, with the upper/lower signs corresponding to the $k=0/1$ states. This test remains effective $|B|>2$ as long as Alice and Bob’s local detector efficiencies, including signal mode transmittances, are above 85%. It does not employ any postselection to filter out any vacuum or multiphoton events, closing the detection loophole and opening up the possibility of performing a loophole-free Bell test [65, 66, 67]. Bell inequality violation is also observed for the states $\sigma\geq 2$, e.g., for $\sigma=2$ one obtains $|B|=2.35$ in the ideal case, but it is not as tolerant to losses.

If a loophole-free Bell test can be established, DI-QKD schemes become possible which are immune to all possible side channels and eavesdropping strategies [44]. The $\sigma=1$ states are particularly well suited for this due to their tolerance to losses. These states are generally multiphoton (see Fig. 3a) but have a dominant single-photon contribution $\frac{1}{\sqrt{2}}(\ket{0}_{A}\ket{1}_{B}\mp i\ket{1}_{A}\ket{0}_{B})$. DI-QKD using these states can then be viewed as a form of single-photon DI-QKD [42]. Once a state has been successfully distributed to Alice and Bob, they perform the same procedure as described in the Bell test above but now Alice has an additional choice for her measurement setting, $\delta_{\alpha 0}=-0.17$. The pair $\{\delta_{\alpha 0}$, $\delta_{\beta 1}\}$ is chosen such that Alice and Bob’s transmitted photon numbers are highly correlated for these measurements, and these results are kept secret and used to generate the raw key. Results for the combination of settings $\{\delta_{\alpha 1},\delta_{\alpha 2},\delta_{\beta 1},\delta_{\beta 2}\}$ are communicated publicly and used to violate the CHSH inequality as described above, which bounds the information that is leaked to Eve [44]. The key rate is proportional to the efficiency of our entanglement distribution protocol, which we have seen scales as $O(\sqrt{\eta})$. Thus single-photon DI-QKD based on our distributed entanglement can achieve the optimal scaling with channel transmittance. In fact, as shown in Fig. 5, the key rates are almost identical to those obtained in TF-QKD if both protocols are run in ideal conditions, i.e., 100% detector efficiencies, 100% error correction efficiency, no dark counts and no background or stray photons. The key rate is calculated assuming collective attacks, following the general procedure outlined in [44], although the security may be improved to protect against completely general/coherent attacks [68, 69]. See Appendix G for details of the calculation and a sketch of the protocol in Fig. 18.

The higher photon-number states $\sigma\geq 2$ also violate Bell inequalities and thus can be used for DI-QKD. However, while the entanglement distribution protocol as a whole has an efficiency scaling as $O(\sqrt{\eta})$, from Eq. (8) the efficiency of extracting a specific state $\rho^{(k,\sigma)}_{\text{out}}$ scales as $O(\eta^{\sigma/2})$. Thus these higher states don’t benefit from such scaling but could still be useful in shorter range networks with lower losses.

If high detection efficiencies and signal mode transmittances ($>85\%$) are not available, a different Bell test may be carried out using postselection and the fair-sampling assumption [43]. Alice and Bob now set $r_{a}=r_{b}=\frac{1}{2}$, and assign a value ${+}1/{-}1$ when photons are detected in the transmitted/reflected channels, ignoring outcomes where neither or both detectors click. Near maximum violation of the CHSH inequality may be obtained if Alice chooses the phase of her coherent state, $\theta_{a}$, randomly from the set $\{0,\frac{\pi}{2}\}$, while Bob chooses his phase, $\theta_{b}$, randomly from $\{\frac{\pi}{4}-\theta,-\frac{\pi}{4}-\theta\}$, where $\theta=-\frac{\pi}{2}$ when $k=0$ and $\theta=\frac{\pi}{2}$ when $k=1$. The amplitudes of the coherent states are chosen to be equal $|\alpha|=|\beta|$. This Bell inequality is no longer loophole-free so it cannot provide DI levels of security, but it can nevertheless be used to perform SEPM-QKD, first described in [43]. If the phase-matching condition $\theta_{a}-\theta_{b}-\theta=\pm\pi$ is satisfied, their results are highly correlated and can be used to extract a shared key, while results for other phases are used to evaluate the CHSH inequality, placing a limit on the information leaked to an eavesdropper Eve. Since the fair-sampling assumption is used, Alice and Bob must now trust their local sources and measurement apparatus, but still need not trust the measurement at the remote station performed by Charlie/Eve, placing the security within the framework of measurement-device-independent (MDI)-QKD [39]. The protocol differs slightly in our case since the entanglement is now heralded by Charlie’s measurement rather than being sent to Alice and Bob directly, but is otherwise identical in operation. The key rate in this case is once again proportional to the efficiency of the entanglement distribution, and thus has the same scaling $O(\sqrt{\eta})$. To compare SEPM-QKD based upon our entanglement distribution to TF-QKD we have considered both in the same realistic scenario discussed in [35], i.e., 30% detector efficiencies and an error correction coefficient of 1.15. Although TES detectors have a negligible dark count rate [70, 52], for a fairer comparison we have also added a small probability of false detection events, $10^{-8}$ per gate interval, which could be due to either background (stray) photons or dark counts if the TES were replaced with single-photon detectors. The results are presented in Fig. 5, and we see that the two key rates remain within an order of magnitude, with SEPM-QKD achieving a slightly larger maximum distance. In the case of SEPM-QKD, this maximum distance results when the false detection probability becomes comparable with the protocol efficiency, and can in principle be increased by using extremely low dark count TES detectors and sufficiently shielding from background photons [71]. Our version of SEPM-QKD may also be viewed as a realistic implementation of the ideal entanglement-based protocol presented in [72]. See Appendix H for more information.

II.5.2 Earth-space communications with $O(\sqrt{\eta})$ scaling

The robustness to extreme losses of our protocol makes it ideal for deployment in a satellite-based quantum communications scheme where entanglement is distributed between two remote locations on Earth and employed for quantum key distribution. Our protocol has the potential to outperform state-of-the-art demonstrations based on the distribution of two-photon Bell pairs, such as those reported in [34] and [29]. There, an SPDC source pumped by a continuous wave laser generated space-borne Bell pairs at random times, which were next sent by a downlink channel to two stations on Earth separated by nearly 1200km. The entanglement between the stations was verified by a CHSH Bell test under the fair sampling hypothesis. Moreover, only approximately two Bell pairs out of $5.9\times 10^{6}$ reached the stations per second. The receiving stations had to be synchronized by the beacon laser and classical two-way communication between Alice and Bob was necessary to discard the cases when only one of the stations received a photon.

In contrast, in our protocol both Alice and Bob create locally an SV state. Their sources can be synchronized offline similarly to entanglement swapping schemes. A third party, Charlie, is located on a satellite and performs conditional state preparation. He acts as a central authority who broadcasts via a classical downlink if Alice and Bob share an entangled state. Thus, direct communication between Alice and Bob is not required to generate the entanglement, although for QKD some one-way communication is still necessary to perform information reconciliation. Most importantly, in a satellite-based approach, the improved scaling $O(\sqrt{\eta})$ is more significant. For example, in [29], two Bell pairs reach the ground stations per second, whereas for similar atmospheric losses (40 dB between ground and satellite, 80 dB summed losses) and $g=0.1$ in our scenario, we estimate 160 successful entangled states to be generated per second.

Our protocol operating in an Earth-space environment is also original in that it uses an uplink rather than a downlink to send the quantum optical signals. Although downlinks have been shown to perform better than uplinks [61], from an experimental point of view uplinks provide better control and flexibility of the quantum source when it is located on the ground rather than in space, and thus are the better choice to study global-scale quantum communications implementations [61]. Turbulence and background radiation affect uplinks more significantly, however, the pointing accuracy is a more consequential factor in our scheme, due to the necessity of time synchronization of the ground sources. In this aspect uplinks have in fact been shown to be more precise [61]. Further analysis on the feasibility of uplink channels has been done in [73] and [74].

II.5.3 Quantum metrology

The entanglement distribution protocol described above can be used to generate states with applications to quantum metrology. $\ket{\Psi^{(k,\sigma)}_{\text{out}}}$ is a generalized HB state, which becomes an exact HB state for $k=\sigma/2$ [75]. It offers quantum-enhanced optical phase estimation even in the presence of significant optical loss and approximates the performance of an optimal probe. Details of this analysis together with an experimental demonstration can be found in [46].

III Discussion

We have proposed an entanglement distribution protocol based on an optimal use of the existing integrated quantum optical components, namely SPDC sources and PNR detectors. An experiment in which our family of states was successfully produced has been reported in [57] and [46], manifesting the feasibility of our protocol. Other distinctive features of our protocol include robustness to arbitrarily high transmission losses, ability of choosing the local dimension of the generated multiphoton entangled state, the possibility of performing a loophole-free Bell test by adopting a preselection instead of a postselection procedure and the flexibility in the choice of the security framework depending on the available detector efficiencies. These allow our protocol to be applicable to quantum metrology, as already shown in [46], Earth-space quantum communications and QKD. This would require space-ready TES detection, which is already being developed [76]. We demonstrated that our protocol is capable of achieving a quadratic improvement in transmission rates compared to the distribution of polarization-entangled photon pairs, such as in the recently deployed Earth-space setup [29]. This result enabled us to employ our distributed entangled states in well-established QKD protocols that match the key rates of state-of-the-art protocols. At smaller distances, such as in metropolitan networks, the generation of robust and near-maximally entangled multiphoton states also opens the possibility of performing high-dimensional QKD, which has recently produced record breaking results in other platforms [77].

IV Acknowledgments

We would like to thank Stefanie Barz, Christine Silberhorn and Tim Bartley for discussions. M.E.M., T.M., A.B. and M.S. were supported by the Foundation for Polish Science “First Team” project No. POIR.04.04.00-00-220E/16-00 (originally FIRST TEAM/2016-2/17) and the National Science Centre ‘Sonata Bis’ project No. 2019/34/E/ST2/00273. Numerical computations were performed in the ACK “Cyfronet” AGH computer centre (Zeus cluster).

Appendix A Photon number statistics

of the output state $\ket{\Psi^{(k,\sigma)}_{\text{out}}}$

In the lossless case, the family of output states generated by our protocol is given by

[TABLE]

where $\sigma$ is the total number of photons distributed between the signal modes and the probability amplitudes $\mathcal{A}_{\sigma}(k,n)$ are expressed using symmetric Kravchuk functions $\phi_{k}$ [57, 58, 78]

[TABLE]

Here $U_{\text{BS}}=U_{\text{BS}}^{(1/2)}$ is the unitary for a beam splitter (BS) of reflectivity $r=\frac{1}{2}$ where $U_{\text{BS}}^{(r)}=\exp(-iH_{\text{BS}}\theta)$ with $H_{\text{BS}}=\frac{1}{2}(a^{\dagger}b+ab^{\dagger})$ and $\theta=2\arcsin(\sqrt{r})$. The Kravchuk functions $\phi_{k}$ are defined in terms of the Gaussian hypergeometric function ${}_{2}F_{1}$ as

[TABLE]

and form an orthonormal set

[TABLE]

The states $\ket{\Psi_{\text{out}}^{(k,\sigma)}}$ then form a complete orthonormal set in the $(\sigma+1)$-dimensional Hilbert space of fixed photon number $\sigma$

[TABLE]

These states are listed in Table I below for $\sigma=1,2,3$, while Figs. 6 and 7 depict photon number distributions $\lvert\mathcal{A}_{\sigma}(k,n)\rvert^{2}=\lvert\phi_{k}(n-\sigma/2,\sigma)\rvert^{2}$ for some higher states, $\sigma=4$ and $\sigma=10$ respectively. Notice that for $k=\sigma/2$, the odd photon number components vanish due to multi-photon Hong-Ou-Mandel interference, effectively reducing the Hilbert space dimension to $(\frac{\sigma}{2}+1)$. Additionally, the probability distribution $\lvert\mathcal{A}_{\sigma}(\sigma/2,n)\rvert^{2}$ possesses an envelope $f(n)=\frac{4}{\pi\sigma\sqrt{1-(2n/\sigma-1)^{2}}}$ which is a probability density function of an arcsine distribution [79].

Appendix B Logarithmic negativity of $\ket{\Psi^{(k,\sigma)}_{\text{out}}}$

The logarithmic negativity of a physical state is a measure of its amount of entanglement. Given a state described by a density operator $\rho$, it is defined by $E_{\mathcal{N}}(\rho)=\log_{2}\left\lVert\rho^{\Gamma}\right\lVert_{1}$, where $\Gamma$ denotes the partial transposition operation and $\left\lVert\rho^{\Gamma}\right\lVert_{1}$ is the trace norm of $\rho^{\Gamma}$. In the lossless case our output state is $\rho_{\text{out}}^{(k,\sigma)}=\ket{\Psi_{\text{out}}^{(k,\sigma)}}\bra{\Psi_{\text{out}}^{(k,\sigma)}}$ and the logarithmic negativity is

[TABLE]

This formula can be derived in the following way

[TABLE]

[TABLE]

The above operator is diagonal and therefore, it is easy to compute the square root of it by taking the square root of its eigenvalues

[TABLE]

In these formulas, $\Gamma_{2}$ denotes the partial transpose with respect to the second subsystem, and $\delta_{n,m}$ denotes the Kronecker delta, equal to $1$ for $n=m$ and [math] otherwise.

Appendix C Protocol performance in case of symmetric losses in idler modes $a_{2}$ and $b_{2}$

C.1 Ideal case

A two-mode squeezed vacuum state is given by

[TABLE]

so that the density operator for the input state $\ket{\Psi_{\text{in}}}=\ket{\Psi}^{\otimes 2}$ is as follows

[TABLE]

Modes $a_{2}$ and $b_{2}$ interfere on a beam splitter and the resulting output modes are projected onto $\ket{k,\sigma-k}_{d_{2},d_{1}}$. In the ideal case, the total number of photons registered behind the BS, $\sigma$, is equal to the total number of photons, $S$, in the remaining modes $a_{1}$ and $b_{1}$. The state created in modes $a_{1}$ and $b_{1}$ takes the form

[TABLE]

where $\mathcal{N}^{2}=\frac{\cosh^{4}g}{\tanh^{2\sigma}g}$. This simplifies to

[TABLE]

C.2 Symmetric idler losses

Here we assume no losses in signal modes $a_{1}$ and $b_{1}$, as well as ideal detection performed by photon-number-resolving (PNR) detectors. More general cases including losses in idler and signal modes are studied in subsequent sections.

Losses can be modelled by a beam splitter with reflectivity $r\in[0,1]$ which quantifies the amount of loss, and then tracing over the reflected mode. For a Fock state $\ket{n}$ it reduces to

[TABLE]

where $t$ and $r$ denote the transmitted and reflected mode, respectively, $r+t=1$, and

[TABLE]

For $r=0$ it reduces to

[TABLE]

Applying this procedure to the idler modes of the input state $\ket{\Psi_{\text{in}}}$ leads to

[TABLE]

where $r_{a}$ and $r_{b}$ denote losses at Alice’s and Bob’s idler modes. As we assume that there are no losses in modes $a_{1}$ and $b_{1}$, the total number of photons in these modes equals $S=n+m=n^{\prime}+m^{\prime}$. The lossy modes $a_{2}$ and $b_{2}$ interfere on a balanced beam splitter and are projected onto $\ket{k,\sigma-k}_{d_{2},d_{1}}$, where $\sigma\leq S$. This results in the creation of the following state

[TABLE]

where

[TABLE]

We now assume that the losses are equal, i.e. $r_{a}=r_{b}=r$. Then,

[TABLE]

This allows us to write

[TABLE]

where we have absorbed the factor $(1-r)^{\sigma}$ into $\tilde{\mathcal{N}}^{2}$. Since $S=n+m=n^{\prime}+m^{\prime}$, we notice that $n,n^{\prime}\leq S$ and $\sigma=S-p-q$. Since $q=S-\sigma-p$, we are able to remove the sum over $q$, which creates constraints on $p$: $p\leq S-\sigma$ and $p\geq n-\sigma$ as well as $p\geq n^{\prime}-\sigma$. The state takes the form

[TABLE]

where $\chi_{\sigma,S}=\frac{r^{S-\sigma}\lambda_{S}}{\tilde{\mathcal{N}}^{2}_{\text{int}}}$. Here we have used $\lambda_{n}\lambda_{S-n}=\frac{\lambda_{S}}{\cosh^{2}g}$ and absorbed the factor $\frac{1}{\cosh^{2}g}$ into $\tilde{\mathcal{N}}^{2}$. Eq. (36) then simplifies to

[TABLE]

The internal matrix is as follows

[TABLE]

where

[TABLE]

By shifting the summation index $p\to n-p$ and then using the Chu-Vandermonde identity $\sum_{n=0}^{S}\binom{n}{p}\binom{S-n}{\sigma-p}=\binom{S+1}{\sigma+1}$, this may be evaluated as $\frac{1}{\tilde{\mathcal{N}}^{2}_{\text{int}}}=\binom{S+1}{\sigma+1}$. This allows one to simplify the form of $\chi_{\sigma,S}$ to

[TABLE]

Furthermore, Eq. (42) now simplifies to

[TABLE]

where we have used the identity $\sum_{S=\sigma}^{\infty}x^{S}\binom{S+1}{\sigma+1}=(1-x)^{-\sigma-2}x^{\sigma}$. The probability to obtain $S$ photons, given a measurement $\ket{k,\sigma-k}$ at the remote station, is then equal to

[TABLE]

and Eq. (41) may be rewritten as

[TABLE]

For $S=\sigma$, as expected, $\rho^{(\sigma,k,S)}_{\text{int}}$ reproduces the density operator for the lossless case

[TABLE]

The probability $p_{S|\sigma}$ quantifies the contribution of $\rho^{(\sigma,k,S)}_{\text{int}}$ to the final density operator. Provided we are in the limit $r\tanh^{2}g\ll 1$, the term $(r\tanh^{2}g)^{S-\sigma}$ decreases with $S$ more rapidly than $\binom{S+1}{\sigma+1}$ increases, so that $p_{S|\sigma}$ is largest for $S=\sigma$ and rapidly decreases. In the realistic limit $g^{2}\ll 1$, this is satisfied for arbitrarily high idler losses $r\to 1$, leading to the output state being largely unchanged by the losses $\rho^{(k,\sigma)}_{\text{out}}\approx\ket{\Psi_{\text{out}}^{(k,\sigma)}}\bra{\Psi_{\text{out}}^{(k,\sigma)}}$.

C.3 Numerical computations

We have computed the logarithmic negativity for the density operator in Eq. (40) assuming losses $r$ in idler modes $a_{2}$ and $b_{2}$ to be symmetric. The signal mode and detector losses, $r_{s}$ and $r_{d}$, have been set to zero. As is customary, the reflectivities in the results below and in the following subsections are displayed in percentages (%).

We have changed $r$ between [math] and $90\%$. The case of $r=99.99\%$ ($40\kern 2.5pt\mathrm{dB}$), which is typical of Earth-to-space scenarios, has been also examined. The results are shown in Fig. 8. The logarithmic negativity is not deteriorated by increasing values of $r$ and even for $r=99.99\%$ it is close to the ideal case $r=0$.

C.4 Numerical computations including lossy detection

We have repeated the above computations for the density operator in Eq. (40) with non-zero losses $r_{d}$ at the PNRs located behind the beam splitter included in the numerical programme. Fig. 9 depicts the values of $E_{\mathcal{N}}$ for $r_{d}$ equal to $5\%$ and $21\%$, which correspond to detector efficiency of $95\%$ and $79\%$ in case of on-chip detection scheme [80]. We have also set the range of $r\in[0,99.99\%]$. Similarly to the previous case, the results are not deteriorated by the losses.

C.5 Comparison with the case of symmetric losses in signal modes $a_{1}$ and $b_{1}$

In order to check the influence of losses $r_{s}$ in signal modes $a_{1}$ and $b_{1}$ we have set them to symmetric values in the range between [math] and $50\%$. The losses at the detectors have been set to $5\%$ and $21\%$ while we assumed no losses at idler modes. The results are depicted in Fig. 10. Non-zero values of $r_{s}$ significantly lower the obtained logarithmic negativity, which drops below $1.0$ for $r_{s}=30\%$, in the case of $g=0.1$, $\sigma=4$. This is not a major problem since the signal beams only travel a short distance to Alice / Bob’s local detectors, compared to the large distance travelled by the idler beams. The signal beams may be measured before the idler beams (the delayed-choice scheme) without affecting the measurement statistics, by the no-signaling principle.

Appendix D Protocol performance for asymmetric losses in idler modes $a_{2}$ and $b_{2}$

D.1 Analytic derivations

In this subsection we assume no losses in signal modes as well as ideal detection performed by all PNR detectors. Moreover, since analytical derivations become quite complex for a general case of losses occurring in our setup, we will first analyse the case when Bob’s idler mode, $b_{2}$, is lossless ($t_{b}=1$ or $r_{b}=0$) and Alice’s idler mode’s transmittance, $t_{a}$, is a fraction $\epsilon$ of Bob’s mode’s transmittance ($t_{a}=\epsilon$). This approach will allow us to understand the impact of asymmetric losses on the protocol’s performance.

We repeat the steps as in Appendix C, but with $r_{a}=1-\epsilon$ to arrive at the form equivalent to Eq. (35)

[TABLE]

[TABLE]

where we could remove the sum over $q$, as only the $q=0$ term would contribute to the sum for $t_{b}=1$.

The Kronecker delta functions tell us that $n+m=n^{\prime}+m^{\prime}=S$, we can then set $m=S-n,m^{\prime}=S-n^{\prime}$ and replace the sum over $m,m^{\prime}$ with a single sum over $S$. They also tell us that $p=S-\sigma$, so that the sum over $p$ can also be removed. Thus, the density matrix takes the following form

[TABLE]

Finally, we notice that $\lambda_{n}\lambda_{S-n}=\frac{\lambda_{S}}{\cosh^{2}g}$, and again absorb $\frac{1}{\cosh^{2}g}$ into $\tilde{\mathcal{N}}^{2}$

[TABLE]

where the normalization factor equals

[TABLE]

We notice that, unlike in Eq. (41), it is impossible to decompose the matrix into $\epsilon$-independent components. Instead we obtain

[TABLE]

where \chi_{\sigma,S}(\epsilon)=\tfrac{\lambda_{S}}{\tilde{\mathcal{N}}(\epsilon)^{2}}\big{(}\tfrac{1-\epsilon}{\epsilon}\big{)}^{S-\sigma} and the density matrix components are

[TABLE]

with

[TABLE]

We find that, similar to the case of symmetric losses, the probability $\tilde{\mathcal{N}}^{2}\chi_{\sigma,S}(\epsilon)$ is strongly peaked at $S=\sigma$ when $g^{2}\ll 1$, so that $\rho^{(k,\sigma)}_{\text{out}}\approx\rho^{(\sigma,k,\sigma)}_{\text{int}}$. However $\rho^{(\sigma,k,\sigma)}_{\text{int}}$ is no longer equal to the lossless state $\ket{\Psi_{\text{out}}^{(k,\sigma)}}\bra{\Psi_{\text{out}}^{(k,\sigma)}}$ but is instead given by $\rho^{(\sigma,k,\sigma)}_{\text{int}}=\ket{\Psi_{\epsilon}^{(k,\sigma)}}\bra{\Psi_{\epsilon}^{(k,\sigma)}}$ where

[TABLE]

The asymmetry lowers the number of photons, $n$, measured by Alice, breaking the $n\to\sigma-n$ symmetry of the state. This is shown for the state $k=0,\sigma=4$ and $g=0.1$ in Fig. 11, where the probability for Alice, Bob to measure $n,\sigma-n$ photons is plotted for $\epsilon\in\{0,0.33,0.66,1\}$. The suppression of high $n$ states leads to an effective reduction of the Hilbert space dimension, lowering the entanglement. Following the same steps as in Note 2, the logarithmic negativity of $\ket{\Psi_{\epsilon}^{(k,\sigma)}}$ is

[TABLE]

This correctly reduces to Eq. (14) in the symmetric limit $\epsilon\to 1$, while in the limit of extreme asymmetry $\epsilon\to 0$ it becomes $E_{\mathcal{N}}=0$, since the state reduces to the unentangled $\ket{0,\sigma}$.

D.2 Numerical computations

To back up our analytic derivations we have computed the logarithmic negativity for the full state given in Eq. (55) for various values of transmittances $t_{a}=\epsilon\cdot 100\%$, while keeping Bob’s mode ideal ($t_{b}=100\%$). As before, we take $g=0.1$ and display transmittances in percentages (%) in the following results. A representative figure showing the behavior of logarithmic negativity as the transmittance of idler modes is made progressively more unequal is shown in Fig. 12. We find that a small asymmetry $\epsilon\gtrsim 70\%$ has almost no effect but thereafter there is a significant reduction in the amount of entanglement.

Further analysis reveals that the entanglement is best preserved for the case where all photons emerge from the same output of the beam splitter, i.e. for $k=0$ and $k=\sigma$. The entanglement deteriorates quickest when the outputs of the beam splitter are equally populated, $k=\sigma/2$. These results are shown in Fig. 13 where the dependence of the logarithmic negativity on $\epsilon=\tfrac{t_{a}}{t_{b}}$ is plotted for $\sigma=4$. The numerically calculated values (points) match almost exactly the approximate analytical result (lines) calculated from Eq. (61). From these calculations it can be inferred that the protocol maintains its quality of output (up to 90% of the maximal value of entanglement) down to $\epsilon=0.7$ for $k=2$ and $\epsilon=0.4$ for $k=0$.

D.3 Numerical computations with both idler modes lossy

We include computations for the more general cases, when Bob’s mode is no longer assumed lossless. The results for $\sigma=4$ are shown in Fig. 14. $t_{b_{2}}$ takes values from 90% to 10% and for each $t_{b_{2}}$, $t_{a_{2}}$ goes from 5% to $t_{b_{2}}$ in steps of 5%. When the transmittances of the two idler modes are equal, the amount of entanglement is close to maximal, given by $\log_{2}(5)\approx 2.32$. As the symmetry is broken, entanglement decreases. In fact, Eq. (61) gives a very good analytical approximation even in this more general case, showing that the entanglement depends only on $\epsilon=\frac{t_{a}}{t_{b}}$.

D.4 General case which assumes losses in all modes and in detection

Finally, we include computations for losses present in all modes and detectors, and asymmetry in idler losses. We have set the losses in the idler modes $a_{2}$ and $b_{2}$ to 37 and 40dB, the detector losses to 50%, while changing losses in signal modes $a_{1}$ and $b_{1}$ between [math] and $50\%$. The results are depicted in Fig. 15, and we see that there is still a high amount of entanglement, even in these unfavourable conditions.

Appendix E Losses varying with time

As has been discussed in ref. [61], in Earth-space scenarios, losses in each mode will vary over time on the timescale of 10-100ms for losses coming from atmospheric turbulence and 0.1-1s for losses deriving from jitter in the telescopes. A single pulse travelling towards a satellite on a low-Earth orbit will take approximately 3ms to reach the satellite. Therefore any fluctuations will only affect the time-averaged logarithmic negativity Alice and Bob compute for all successful events.

Let us assume that Bob’s modes’s transmittance is fixed and the Alice’s has a normal distribution with variance 1dB. We then generate 500 random values from that distribution and calculate the resulting logarithmic negativity. We repeat the process for different values of transmittance $t_{b_{2}}$. An example is shown in Fig. 16 for $\sigma=4$ and typical attenuation of 40dB.

Appendix F Efficiencies and success rates

F.1 Ideal case

Interference of modes $a_{2}$ and $b_{2}$ on the BS alters the input state $\ket{\Psi_{\text{in}}}=\ket{\Psi}^{\otimes 2}$ in the following way

[TABLE]

After tracing out modes $a_{1}$ and $b_{1}$ we obtain

[TABLE]

Thus, the probability of detecting $\ket{k,\sigma-k}$ behind the BS, generating the state $\rho_{\textrm{out}}^{(k,\sigma)}$, equals

[TABLE]

Since the beam splitter interaction is particle number conserving, $\sigma=n+m$ must hold true. In addition, we note that $\lambda_{a}\lambda_{b}=\frac{\lambda_{a+b}}{\cosh^{2}g}$ and therefore,

[TABLE]

This probability depends solely on the parametric gain $g$ and $\sigma$. Since there are $\sigma+1$ values of $k$, the probability to generate any $\sigma$ photon state is $(\sigma+1)p_{(k,\sigma)}$. Fig. 17 shows $p_{(k,\sigma)}$ for $g=0.1$.

The success rate of obtaining a particular total photon number $\sigma$ is then calculated from multiplying $(\sigma+1)p_{(k,\sigma)}$ and the repetition rate of the pump, $f_{\mathrm{rep}}$. For example, assuming $g=0.1$ and $f_{\mathrm{rep}}=80$MHz, we would observe events having $\sigma=4$ photons in total with frequency 3.82Hz and events having $\sigma=2$ with frequency 23.2kHz.

F.2 Symmetric idler losses

Including idler losses, the probability $p_{(k,\sigma)}$ may be calculated similarly by

[TABLE]

where $\rho_{\text{in}}$ is given by Eq. (34) and $U_{\text{BS}}$ is understood to act on the modes $a_{2},b_{2}$. The state we are tracing over has already been calculated in Appendix C, for the case of symmetric losses $r_{a}=r_{b}=r$. The trace should then be given by the normalization constant $\frac{1}{\tilde{\mathcal{N}}^{2}}$, but we subsequently absorbed a factor $\frac{(1-r)^{\sigma}}{\cosh^{2}(g)}$ into $\tilde{\mathcal{N}}^{2}$. The probability is then

[TABLE]

which correctly reduces to Eq. (69) in the lossless case $r=0$.

The protocol efficiency is proportional to the probability of receiving at least one photon at the remote station, as then we are certain that Alice and Bob share entangled states. The complement of that event is that no photons reach the satellite, equivalent to $\sigma=0$. Setting $\sigma=0$ in Eq. (71) we find

[TABLE]

The equivalent expression in the case of asymmetric losses is calculated as follows

[TABLE]

The overall success rate is then $f_{\mathrm{rep}}(1-p_{(0,0)})$. For a typical parametric gain ($g=0.1$) and Earth-space losses ($r_{a}=r_{b}=40$dB) we obtain $1-p_{(0,0)}=2\times 10^{-6}$. Assuming a pulse repetition rate of $80$MHz, our success rate will be $160$Hz, i.e. 160 entangled states generated per second.

Appendix G Loophole-free Bell test and DI-QKD

G.1 Loophole-free Bell test

The entanglement of Alice and Bob’s state $\rho^{(k,\sigma)}_{\textrm{out}}$ can be verified by a Bell test which was first introduced for generic two-mode states of light in ref. [81]. Alice and Bob interfere their individual modes with coherent states $\ket{\alpha}$ and $\ket{\beta}$, on beam splitters of reflectivities $r_{a}$ and $r_{b}$ (Alice and Bob’s idler losses are now assumed equal and given by $r$). These operations can be approximated by displacement operations $D(-\delta_{\alpha})$ and $D(-\delta_{\beta})$, where $\delta_{\alpha}=i\alpha\sqrt{r_{a}}$, $\delta_{\beta}=i\beta\sqrt{r_{b}}$, which become exact in the limit $r_{a},r_{b}\to 0$ and $|\alpha|^{2},|\beta|^{2}\to\infty$ such that $\delta_{\alpha}$ and $\delta_{\beta}$ are finite. They then measure their transmitted photon numbers $i$ and $j$, separating their results into two categories (dichotomized variables). This can be done by assigning $+1/-1$ either to vacuum/non-vacuum events, or to even/odd photon numbers. By choosing randomly between two local measurement settings each, $\delta_{\alpha 1},\delta_{\alpha 2}$ and $\delta_{\beta 1},\delta_{\beta 2}$, they can observe violation of the CHSH inequality $|B|\leq 2$ where $B=E(\delta_{\alpha 1},\delta_{\beta 1})+E(\delta_{\alpha 1},\delta_{\beta 2})+E(\delta_{\alpha 2},\delta_{\beta 1})-E(\delta_{\alpha 2},\delta_{\beta 2})$, and $E(\delta_{\alpha},\delta_{\beta})$ is the correlation between their dichotomized variables.

First let us consider the zero/non-zero case. It is convenient to write $B$ in the Clauser–Horne form

[TABLE]

Here $p(i=0,j=0|\delta_{\alpha},\delta_{\beta})$ are the probabilities that they both detect zero photons, given measurement settings $\delta_{\alpha}$, $\delta_{\beta}$. The use of TES detectors (transition edge sensors) as PNRs make the detection of vacuum events much easier by significantly reducing the dark count rate. Alternatively, swapping the labelling of the outputs has no effect on $B$, so an equivalent form may be obtained by replacing $i=0,j=0$ with $i\neq 0,j\neq 0$ throughout this equation, which may be a more practical form for experiments. In the case of perfect detectors, the probability $p(i=0,j=0|\delta_{\alpha},\delta_{\beta})$ is equal to the Husimi Q function of the state

[TABLE]

where following [81] we have dropped unnecessary factors of $\pi$ from the usual definition. Similarly, the marginal probabilities $p(i=0|\delta_{\alpha})$ and $p(j=0|\delta_{\beta})$ are given by the marginal Q functions

[TABLE]

and similarly for Bob, where $\textrm{Tr}_{b}$ indicates a partial trace over Bob’s modes. With zero idler losses, the Q functions of the state $\ket{\Psi^{(k,\sigma)}_{\textrm{out}}}$ are

[TABLE]

These forms are extremely useful for finding the optimal measurement settings that maximize the Bell inequality violation.

Similarly, in the even/odd test we can express the Bell expression in terms of the Wigner function. The correlation may be written as an expectation value E(\delta_{\alpha},\delta_{\beta})=\textrm{Tr}\Big{(}\rho^{(k,\sigma)}_{\textrm{out}}\hat{W}(\delta_{\alpha},\delta_{\beta})\Big{)} of the operator

[TABLE]

where $\hat{\Pi}^{(+)}_{a}$ and $\hat{\Pi}^{(-)}_{a}$ are measurement operators for Alice to measure even and odd photon numbers respectively

[TABLE]

and similarly for Bob. Using the parity operator

[TABLE]

$\hat{W}(\delta_{\alpha},\delta_{\beta})$ may be simplified to

[TABLE]

Thus $E(\delta_{\alpha},\delta_{\beta})$ is the expectation value of this displaced parity operator, which is one of the equivalent definitions of the Wigner function [81]. To make this clear we write $E(\delta_{\alpha},\delta_{\beta})=W(\delta_{\alpha},\delta_{\beta})$, which again is true only for perfect photon detection. With zero idler losses, the Wigner function of the state $\ket{\Psi^{(k,\sigma)}_{\textrm{out}}}$ is

[TABLE]

where $L_{k}$ are Laguerre polynomials. This may be calculated by noting that $\ket{\Psi_{\text{out}}^{(k,\sigma)}}=U_{\textrm{BS}}\ket{k}\ket{\sigma-k}$ (this follows from the symmetry $A_{\sigma}(n,k)=A_{\sigma}(k,n)$), and evolving the Wigner function for the two-mode Fock state.

For each state we have numerically maximized the Bell inequality with respect to the four measurement settings $[\delta_{\alpha 1},\delta_{\alpha 2}],[\delta_{\beta 1},\delta_{\beta 2}]$. Table II lists the optimal choice of dichotomized variables and settings for each state, and the resulting value of $|B|$. We have also seen how this value is affected by idler losses. Since the entanglement decreases and plateaus after some finite amount of losses (see Fig. 4), the amount of Bell violation displays the same behaviour. For definiteness we have investigated typical Earth-space idler losses of 40dB (80dB summed losses in the Alice-Bob channel), with $g=0.1$, although due to the plateau effect these losses can be considered the limiting behaviour as $r\to\infty$. We have also determined the minimum efficiency of Alice and Bob’s local photon detectors required to still observe some violation.

The $\sigma=1$ states have the most robust violation in the presence of losses, $|B|=2.63$ and require detector efficiencies $>85\%$. States equivalent to our $\sigma=1$ states in the lossless limit have been explored in ref. [81], although they used a less optimal choice of settings, obtaining $|B|\approx 2.4$. The higher states $\sigma\geq 2$ also violate these Bell inequalities, but the violation is fragile with respect to signal mode losses and detector efficiencies. In ideal conditions, the states $k=\frac{\sigma}{2}$ have the lowest amount of Bell violation due to the lower amount of entanglement in these states. Some states e.g. $(k=1,\sigma=3)$ and $(k=2,\sigma=3)$ violate inequalities for zero idler losses but fail to do so with arbitrarily high losses. For these specific examples, the maximum amount of losses is around 50%, or 3dB. It is expected that the states $\sigma\geq 2$ violate some higher dimensional Bell inequalities beyond the simple CHSH scenario, which will be more loss tolerant, but this is still an active area of research.

G.2 Device-independent (DI)-QKD

If Alice and Bob share entangled states and can violate some Bell inequalities, it is possible for them to extract a shared cryptographic key. If the Bell test is loophole-free, they can place a limit on the amount of information obtained by an eavesdropper Eve and ensure that their key is secure, immune to all possible side channels and eavesdropping strategies. This is called device-independent quantum key distribution (DI-QKD). We follow closely the formulation presented in ref. [44], where we can think of Charlie sending Alice and Bob the state $\rho_{\textrm{out}}^{(k,\sigma)}$ depending on their measurement, $\ket{k,\sigma-k}$ at the remote station. The protocol is secure even if Charlie is replaced by the adversarial Eve.

As discussed in the main text, we focus on long-distance QKD and thus $\sigma=1$ states, since they display success rates scaling with transmittance as $O(\sqrt{\eta})$, and are more tolerant to signal mode losses and detector inefficiencies. In the long distance limit, these states form the vast majority of generated states, e.g. for our typical Earth-space parameters with 40dB losses and $g=0.1$, $p_{(k,1)}/(1-p_{(0,0)})=99.9999\%$. In the absence of losses these states are exactly the maximally entangled single-photon states $\frac{1}{\sqrt{2}}(\ket{0,1}_{a_{1},b_{1}}\mp i\ket{1,0}_{a_{1},b_{1}})$. These have been considered in a single-photon entanglement DI-QKD scheme described in ref. [42], the protocol developed here can be considered equivalent, but using the realistic state $\rho_{\textrm{out}}^{(k,\sigma)}$ rather than the idealized states, and with improved parameters and key rates.

The setup of the QKD protocol and a sequence diagram are displayed in Fig. 18. Upon receipt of the entangled state $\rho^{(0,1)}_{\textrm{out}}$ or $\rho^{(1,1)}_{\textrm{out}}$, Alice and Bob each perform a measurement as described in the previous section, by performing displacement operations with parameters $\delta_{\alpha}$ and $\delta_{\beta}$ respectively, and measuring the output photon numbers $i$ and $j$. Alice randomly chooses her setting $\delta_{\alpha}$ from three possible choices $[\delta_{\alpha 0},\delta_{\alpha 1},\delta_{\alpha 2}]$ while Bob chooses randomly from two choices $[\delta_{\beta 1},\delta_{\beta 2}]$. They repeat this process a number of times and then communicate publicly which settings they chose for each trial.

The pair $\{\delta_{\alpha 0}$, $\delta_{\beta 1}\}$ are chosen such that $i$ and $j$ are highly correlated for these measurements, and these results are kept secret and used to generate the raw key. The combination of settings $\{\delta_{\alpha 1},\delta_{\alpha 2},\delta_{\beta 1},\delta_{\beta 2}\}$ are chosen such that the measurements violate the CHSH inequality. Alice and Bob communicate through an authenticated channel and reveal the results $i$ and $j$ for all combinations of these settings and ensure that the inequality Eq. (74) is violated. This ensures that these measurement results are random i.e. that a potential eavesdropper Eve cannot predict the results with certainty. The amount of randomness increases with increasing Bell violation, with Eve’s probability of guessing the result bounded by $p_{\textrm{guess}}\leq\frac{1}{2}(1+\sqrt{2-B^{2}/4})$ [82]. Since $p(j|\delta_{\beta 1})$ doesn’t depend on Alice’s setting by the non-signalling principle, Eve also cannot have complete information about the results for $\{\delta_{\alpha 0}$, $\delta_{\beta 1}\}$ or the raw key generated by them.

Alice and Bob now both have raw keys of $n$ bits, but they may not be perfectly correlated, and Eve may have a small amount of information about them. They thus carry out a process of classical information reconciliation (error correction) and privacy amplification, discarding a number of bits to obtain the final shared key of $\ell$ bits. The raw key rate is defined as the limiting ratio

[TABLE]

and depends on the type of attacks we consider from Eve. Eve is assumed to control the source (for us, she is in control of the remote station measurement), and also to have fabricated Alice and Bob’s measurement devices. Following [44], we consider collective attacks whereby Eve interacts a fresh system with Alice and Bob’s idler modes each time they reach the remote station. She applies the same unitary operation each time, and can perform a coherent measurement on all her systems at any time. The security may be improved to protect against completely general/coherent attacks by using the results from refs. [68, 69]. However this modification will not significantly reduce the key rate, which is dominated by the efficiency of the entanglement distribution, which remains unchanged.

The security proof presented in ref. [44] applies provided that the measurements are described by Hermitian operators with eigenvalues $\pm 1$. This is true in our case as we now prove. The event where Alice measures zero transmitted photons after interfering with the coherent state $\ket{\alpha}$ is assigned the value $+1$ and has probability

[TABLE]

where $\rho_{a}=\operatorname{Tr}_{b}(\rho_{\textrm{out}}^{(k,\sigma)})$ is Alice’s reduced density matrix. Thus $\ket{\delta_{\alpha}}\bra{\delta_{\alpha}}$ is a projection operator (and POVM element) for the outcome $+1$, and the corresponding operator for outcome $-1$ is $1-\ket{\delta_{\alpha}}\bra{\delta_{\alpha}}$. The measurement can then be described by the Hermitian operator

[TABLE]

and the security proof in ref. [44] then directly follows.

The raw key rate for collective attacks may be calculated by the Devetak–Winter formula [83]

[TABLE]

where $I(A_{0}:B_{1})$ is the mutual information between Alice and Bob when they use settings $\{\delta_{\alpha 0},\delta_{\beta 1}\}$, and $\chi(B_{1}:E)$ is the Holevo quantity which upper-bounds the mutual information between Bob and Eve. This equation assumes the classical communication is performed one way from Bob to Alice. This is because $\chi(A_{0}:E)>\chi(B_{1}:E)$ holds for this protocol, i.e. Eve has more information about Alice’s results than Bob’s, the asymmetry arising from the fact Alice uses three settings while Bob uses only two. The mutual information $\chi(B_{1}:E)$ is upper-bounded [44] by

[TABLE]

where $h(p)$ is the Binary entropy function. Thus we have a lower bound for the raw key rate

[TABLE]

and all that is left to do is calculate this quantity for our protocol and maximize over the five measurement settings.

The mutual information is

[TABLE]

where $H(A_{0})$ and $H(B_{1})$ are the entropies of Alice and Bob’s marginal probabilities

[TABLE]

and $H(A_{0},B_{1})$ is the entropy of the joint probabilities

[TABLE]

In the case of perfect photon detectors, these can be conveniently described in terms of the $Q$ functions using Eqs. (75) and (76), and writing

[TABLE]

The Bell expression $B$ is given by Eq. (74) as before, which can similarly be expressed in terms of the $Q$ functions.

One can now simply perform a numerical maximization of $K$ as a function of our five measurement settings, but we will take a more instructive approach. As before, we will calculate the optimal parameters for the lossless case and then consider the effect of losses after. The parameters $\{\delta_{\alpha 1},\delta_{\alpha 2},\delta_{\beta 1},\delta_{\beta 2}\}$ must be chosen to maximize the Bell parameter $|B|$ so as to decrease Eve’s knowledge about the key $\chi(B_{1}:E)$. At the same time $\{\delta_{\alpha 0}$, $\delta_{\beta 1}\}$ must be chosen to maximize the correlations between Alice and Bob’s key bits $I(A_{0}:B_{1})$. These cannot be done independently since $\delta_{\beta 1}$ appears in both, but we can get a very good approximation by first maximizing $B$ to find $\delta_{\beta 1}$ and then choosing $\delta_{\alpha 0}$. The parameters $\{\delta_{\alpha 1},\delta_{\alpha 2},\delta_{\beta 1},\delta_{\beta 2}\}$ have thus already been calculated and given in Table 2.

There are two ways for the results for $\{\delta_{\alpha 0}$, $\delta_{\beta 1}\}$ to be perfectly correlated. Either $p(i=0,j=0)+p(i\neq 0,j\neq 0)=1$, so that Alice and Bob always have equal results, or $p(i=0,j=0)+p(i\neq 0,j\neq 0)=0$, so that they always have opposite results. Let’s first look at the state $(k=0,\sigma=1)$, the probabilities are given by our Q function expressions

[TABLE]

Thanks to the interference term, the first of these can be set to zero by setting $\delta_{\alpha 0}=i\delta_{\beta 1}$, which also reduces the second probability, simplifying to

[TABLE]

where we have expanded around small $|\delta_{\beta 1}|$. Reading from Table 2 we have $\delta_{\beta 1}=0.17i$, so we should choose $\delta_{\alpha 0}=-0.17$ and thus we find $p(i=0,j=0)+p(i\neq 0,j\neq 0)=4\times 10^{-4}$. Thus Alice and Bob’s results are almost perfectly correlated, if Alice measures $i\neq 0$, Bob always measures $j=0$ and vice-versa, one party just needs to swap their results to obtain the same key bit. We thus define the quantum bit error rate $q$ as

[TABLE]

which will become non-zero once we include losses. Interestingly, Alice and Bob’s marginal probabilities are also nearly perfectly balanced between the two outcomes

[TABLE]

so that they have near maximum entropy $H(A_{0})=H(B_{1})=1-10^{-7}$. The joint entropy and mutual information are $H(A_{0},B_{1})=1.005$ and $I(A_{0}:B_{1})=0.995$. Inserting our value for the Bell expression $B=2.688$ into Eq. (89) we find the Holevo bound $\chi(B_{1}:E)\leq 0.290$ and thus raw key rate for the lossless case

[TABLE]

The key rate per pulse $R$ (number of secure key bits per second divided by the repetition rate of the pump $f_{\textrm{rep}}$) is found by multiplying the dimensionless key rate $K$ by the number of raw key bits generated per pulse. The probabilities of using the settings $\{\delta_{\alpha 1},\delta_{\alpha 2},\delta_{\beta 2}\}$ can be made very small such that the majority of Alice and Bob’s measurements create raw key bits using the settings $\{\delta_{\alpha 0}$, $\delta_{\beta 1}\}$. This leaves a small fraction of other measurements, but still enough to reliably calculate the Bell inequality. Thus the number of raw key bits generated per pulse is approximately equal to the protocol efficiency $(1-p_{(0,0)})$. Thus $R$ is given by

[TABLE]

with $K$ given by Eq. (90), and we have inserted $p_{(0,0)}$ from Eq. (72). In the lossless case we obtain

[TABLE]

The calculation for the state $(k=1,\sigma=1)$ proceeds in the same way except that we take $\delta_{\alpha 0}=-i\delta_{\beta 1}$. The results are summarized in Table 3.

We will see that the raw key rate $K$ is not significantly affected by idler losses, but the true rate $R$ decreases with the same scaling as the protocol efficiency, $O(\sqrt{\eta})=O(1-r)$. Although we have calculated the optimal parameters $[\delta_{\alpha 0},\delta_{\alpha 1},\delta_{\alpha 2}],[\delta_{\beta 1},\delta_{\beta 2}]$ assuming the lossless scenario, we find that they remain very close to optimal when losses are included.

Inserting idler losses of 40 dB, Alice and Bob’s marginals are still nearly perfectly balanced $H(A_{0})=H(B_{1})=1-10^{-4}$, while the joint entropy and mutual information are only slightly modified $H(A_{0},B_{1})=1.073$, $I(A_{0}:B_{1})=0.927$. This corresponds to a quantum bit error rate of $q=0.01$. We already calculated $|B|$ for 40 dB losses and found $|B|=2.635$, which yields a Holevo bound $\chi(B_{1}:E)\leq 0.371$. Thus for 40 dB idler losses and perfect photon detectors we have a raw key rate

[TABLE]

and multiplying by the protocol efficiency, $1-p_{(0,0)}=2\times 10^{-6}$, the key rate per pulse is

[TABLE]

The results are summarized in Table 4.

So far we have assumed perfect efficiency of Alice and Bob’s local detectors. Non-zero detector losses reduce $K$ by the reduction of the Bell violation (increasing $\chi(B_{1}:E)$), and by an increase of the quantum bit error rate $q$ (decreasing $I(A_{0}:B_{1})$). Assuming 95% efficiency, we obtain marginal entropies $H(A_{0})=H(B_{1})=0.999$, joint entropy $H(A_{0},B_{1})=1.290$ and mutual information $I(A_{0}:B_{1})=0.708$, corresponding to a quantum bit error rate $q=0.055$. The Bell violation remains large $|B|=2.430$, yielding $\chi(B_{1}:E)\leq 0.622$. Thus we find, in the presence of 40 dB idler losses (80 dB Alice-Bob losses) and 95% detector efficiency

[TABLE]

Fig. 19 shows how the key rate decreases as a function of transmission losses and detector efficiencies.

Appendix H Bell test with post-selection and SEPM-QKD

H.1 Bell test with post-selection

If high detection efficiencies are not available for Alice and Bob’s local detectors, a different Bell test may be carried out using post-selection [43]. The setup is the same as before (Fig. 19a), with Alice and Bob interfering with coherent states but they now measure both the transmitted and reflected photon numbers with detectors $A_{1},A_{2}$ and $B_{1},B_{2}$. We again focus on the $\sigma=1$ states due to the increased efficiency of their generation.

After measuring $\sigma=1$ at the remote station, Charlie has effectively distributed the following state to Alice and Bob

[TABLE]

where $\theta=\mp\pi/2$ for $k=0/1$. In our analytical results we will use this state generated with zero idler losses, but in our numerical results we will use the full state $\rho_{\textrm{out}}^{(k,1)}$, with the understanding that in the limit $g^{2}\ll 1$ they become identical. Alice interferes her mode with a coherent state $\ket{\alpha}$, where $\alpha=|\alpha|e^{i\theta_{a}}$, on a 50/50 beam splitter $U_{A}$. Bob similarly interferes his mode with a coherent state $\ket{\beta}$, where $\beta=|\beta|e^{i\theta_{b}}$ on a 50/50 beam splitter $U_{B}$. After the interactions they then share the state

[TABLE]

They each assign a result $+1$ to events where photons are detected exclusively in the transmitted channel (detector $A_{1}$/$B_{1}$ clicks for Alice/Bob respectively), and a result $-1$ to events where photons are detected exclusively in the reflected channel (detector $A_{2}$/$B_{2}$ clicks). They ignore outcomes where neither or both detectors click i.e., they apply the fair-sampling assumption in the usual way.

Let’s consider the event where $A_{1}$ and $B_{1}$ click. The primary contribution to the probability $p(A_{1},B_{1})$ is from single photon events, and assuming $|\alpha|=|\beta|$ this probability takes the simple form

[TABLE]

If Alice and Bob used PNR detectors they could post-select on these single photon outcomes and this would be the total probability $p(A_{1},B_{1})$. If however they use non-photon-number-resolving detectors, such as superconducting nanowires, there will be extra contributions from multi-photon events. Retaining terms up to $|\alpha|^{4}$ the extra contributions come from

[TABLE]

so that the total probability is

[TABLE]

Similar calculations follow for the other three events $(A_{1},B_{2})$, $(A_{2},B_{1})$ and $(A_{2},B_{2})$. In general, the probability that detectors $A_{i}$ and $B_{j}$ click is

[TABLE]

The term proportional to $[1+(-1)^{i+j}\cos(\theta_{a}-\theta_{b}-\theta)]$ allows for perfect correlations while the remaining term introduces a small error due to multi-photon events.

We now consider the effect of losses. In contrast to ref [43], our states $\rho_{\textrm{out}}^{(k,\sigma)}$ are heralded at Alice and Bob via a pre-selection measurement rather than being sent directly through free space. As we’ve seen, this means the idler losses have a negligible effect on the state so that $\rho_{\textrm{out}}^{(k,\sigma)}\approx\ket{\Psi_{\text{out}}^{(k,\sigma)}}\bra{\Psi_{\text{out}}^{(k,\sigma)}}$, which comes at the cost of reduced efficiency of the protocol. In contrast, signal mode losses $r_{s}$, including detector inefficiencies of Alice and Bob’s detectors, lower the amount of entanglement. With non-perfect signal mode transmittance $\eta_{s}=1-r_{s}$, the state reaching the local beam splitters is then

[TABLE]

Adding the extra contribution from the vacuum interfering with the coherent states

[TABLE]

Eq. (117) is modified to

[TABLE]

A slightly different result is obtained in ref. [43] due to them truncating the coherent states at the single photon level, but this doesn’t affect their conclusions.

The expectation value of the product of their outcomes can then be calculated

[TABLE]

and using the measurement settings $\theta_{a}=\{0,\frac{\pi}{2}\}$ and $\theta_{b}=\{\frac{\pi}{4}-\theta,-\frac{\pi}{4}-\theta\}$, one can evaluate the CHSH expression

[TABLE]

In the limit $|\alpha|^{2}\to 0$ this becomes a maximal violation of $B=2\sqrt{2}$, however the combination of multi-photon events and non-perfect transmittance lowers the violation. The transmittance required to observe $|B|>2$ is

[TABLE]

which for e.g. $|\alpha|=0.35$, is $\eta_{s}>22\%$, much better than the $85\%$ efficiency required in the loophole-free case. This of course comes at the cost of having to post-select on a small amount of events. The fraction of events that are kept is given by

[TABLE]

equal to $4\%$ for $|\alpha|=0.35$ and $\eta_{s}=30\%$. Note that this probability does not take into account the probability of generating the starting state $\ket{\Psi_{\text{out}}^{(k,1)}}$, which is approximately equal to the entanglement distribution protocol efficiency $1-p_{(0,0)}$. This is considered a pre-selection rather than a post-selection, since if Charlie measures $\sigma=0$ and communicates this, Alice and Bob don’t share an entangled state and don’t perform any measurements.

H.2 Single-photon entanglement-based phase-matching (SEPM)-QKD

The above Bell inequality is no longer loophole-free so we do not have DI levels of security, but it can nevertheless be used to limit the information leaked to Eve in a MDI-QKD protocol described in ref. [43]. This protocol is called single-photon entanglement-based phase-matching (SEPM)-QKD. Applying this protocol to our case requires slight modifications due to the states $\rho_{\textrm{out}}^{(k,\sigma)}$ being heralded at Alice and Bob via a pre-selection measurement rather than being sent directly through free space.

The coherent state phases are expressed as $\theta_{a}=\phi_{a}+k_{a}\pi$ and $\theta_{b}=-\theta+\phi_{b}+k_{b}\pi$. The phases $\phi_{a(b)}$ are both randomly chosen from the set $\{-\frac{\pi}{4},0,\frac{\pi}{4},\frac{\pi}{2}\}$ while $k_{a(b)}$ are random bits $\{0,1\}$. The settings used above for the Bell test are contained within this set of measurements, but there are also additional events where the phases $\phi_{a}$ and $\phi_{b}$ match. The shared key bits are extracted from these events. Setting $\phi_{a}=\phi_{b}$ one has $\cos(\theta_{a}-\theta_{b}-\theta)=(-1)^{k_{a}+k_{b}}$, so that if Alice and Bob’s key bits are equal $k_{a}=k_{b}$, their measurement probabilities $p(A_{i},B_{j})$ are

[TABLE]

Similarly, if their key bits are different $k_{a}\neq k_{b}$ we have the opposite scenario

[TABLE]

After the measurements are complete they both publish the phases $\phi_{a(b)}$ and measurement results $i,j$ in each trial but keep the key bits $k_{a(b)}$ secret. From the above equations we see that for events where they obtain the same outcome $i=j$, their key bits are the same with high probability, and for events where they obtain different outcomes $i\neq j$, their key bits are different with high probability, so Bob simply flips his key. The error rate is small and given by (e.g. for $k_{a}=k_{b}$)

[TABLE]

From the published information, Eve can determine $|k_{a}-k_{b}|$ for each trial but cannot determine the exact values $k_{a},k_{b}$. The raw key rate considering attacks by Eve is

[TABLE]

The mutual information between Alice and Bob is $I(A:B)=1-h(q)$ while the Holevo quantity $\chi(AB:E)$ is the maximum amount of information obtainable by Eve. In ref. [43] they determine the optimal collective attack by Eve and find that Eve’s information gain is limited to twice the quantum bit error rate $\chi(AB:E)\leq 2q$. They also consider a beam splitting attack which must be treated separately to collective attacks, however we find that the information gain is negligible, being proportional to $|\alpha|^{4}$, whereas $q$ is proportional to $|\alpha|^{2}$.

The measurement results for the other combinations of settings may be used to evaluate the CHSH inequality and, since $B\approx(1-2q)2\sqrt{2}$, the amount of violation can be used to estimate Eve’s information gain and the quantum bit error rate. After performing error correction, they then have a shared secure key. Since the fair sampling approximation was used to evaluate the Bell inequality, the security is no longer device-independent. Alice and Bob must now trust their local sources and measurement apparatus, but still need not trust the measurement at the remote station performed by Charlie/Eve, placing the security within the framework of measurement-device-independent (MDI)-QKD.

The key rate per pulse is found by multiplying $K$ by the number of raw key bits generated per pulse. There is the probability of post-selection $\sum_{ij}p(A_{i},B_{j})$ but in our case there is also the probability of distributing the entangled state, i.e. pre-selection, $1-p_{(0,0)}$. We then obtain

[TABLE]

where we have allowed for an inefficiency $f$, for the error correction. Due to the pre-selection, the term $\sum_{ij}p(A_{i},B_{j})$ depends only on the signal mode transmittance, with the idler mode transmittance affecting only the term $1-p_{(0,0)}$. This is a key benefit, since in ref. [43] they were restricted to the regime $|\alpha|^{2}<1-r$ where the key rates were much lower than in state-of-the-art protocols. Our pre-selection means we instead have the less strict limit $|\alpha|^{2}<\eta_{s}$, allowing higher values of $|\alpha|$ and higher key rates.

We have numerically simulated this SEPM-QKD protocol and compared it to the case of realistic TF-QKD, with the results being presented in Fig. 5. The simulation parameters were chosen to match those in ref. [35]: 30% detector efficiencies at Alice and Bob’s local detectors and an error correction coefficient $f=1.15$. The coherent state amplitude was chosen as $|\alpha|=0.35$, which we found to be close to optimal. For a more fair comparison we have also included a small probability of false detection events $p_{\textrm{dark}}=10^{-8}$, which could be due to dark counts or background (stray) photons. The false detection rate per unit time is then $p_{\textrm{dark}}f_{\textrm{rep}}$ where $f_{\textrm{rep}}$ is the repetition rate of the pump. The main effect of these false detections is that sometimes Charlie incorrectly declares that the state $\rho_{\textrm{out}}^{(k,1)}$ has been distributed due to a false detection of $\sigma=1$. If the protocol efficiency becomes comparable or lower than the false detection efficiency, $(1-p_{(0,0)})<p_{\textrm{dark}}$, the state can then no longer be distributed and the key rate drops to zero. This is responsible for the maximum distance shown in Fig. 5 The other key rates plotted on this figure are obtained with equations presented in [35]. Lines I-III are three QKD protocols with key rates scaling as $O(\eta)$. I: Ideal BB84, $R=\eta$, II: Decoy state QKD, $R=\frac{\eta}{e}$, III: Decoy-state measurement-device-independent (MDI)-QKD, $R=\frac{\eta}{2e^{2}}$

Finally we note the similarity between this protocol and the ideal entanglement-based one presented in ref. [72] (protocol one). There Alice and Bob each prepare a state $\ket{\Phi}=\sqrt{q}\ket{0,0}_{a1,a2}+\sqrt{1-q}\ket{1,1}_{a1,a2}$ where $0\leq q\leq 1$ and send the idler beams to Charlie who interferes them on a beam splitter and measures the output photon numbers. Alice and Bob then share an entangled state, and perform measurements either in the $z$-basis, $\ket{0},\ket{1}$, or the $x$-basis, $\ket{\pm}=(\ket{0}\pm\ket{1})/\sqrt{2}$. The QKD protocol presented here can be considered equivalent, with realistic interference measurements replacing the ideal $x$-basis qubit measurement, and the two-mode squeezed vacuum state approximating $\ket{\Phi}$.

Appendix I Software

The software that we have developed and used was written in Python 3 using the NumPy package, which was employed for matrix algebra. The program performs operations described by the following algorithm

for both SPDC sources belonging to Alice and Bob, compute $\rho_{\Psi}=\ket{\Psi}\bra{\Psi}$ where, $\ket{\Psi}=\sum_{n=0}^{n_{\text{max}}}\sqrt{\lambda_{n}}\,\lvert n,n\rangle$, and $n_{\text{max}}$ is a sum cutoff found by solving $\lambda_{n_{\text{max}}}/\lambda_{0}<10^{-15}$ for a given $g$ and $10^{-15}$ results from the precision of the numeric data type; usually $n_{\text{max}}\geq 6$. 2. 2.

for all modes leaving the sources, add losses and then trace out unwanted modes, i.e. perform the following operation

[TABLE]

where $r_{s}$ and $r_{i}$ are losses in signal and idler modes, respectively, 3. 3.

construct the density operator matrix for $\rho_{\text{in}}=\rho_{\Psi^{\prime}}\otimes\rho_{\Psi^{\prime}}$, 4. 4.

apply the $50:50$ BS operation, $\rho_{\text{BS}}=U_{\text{BS}}\,\rho_{\text{in}}\,\left(U_{\text{BS}}\right)^{\dagger}$, 5. 5.

add losses $r_{d_{1}}$ and $r_{d_{2}}$ to the modes leaving the BS

[TABLE] 6. 6.

for all $k=0,\dots,\sigma$ perform the following operations:

1. (a)

project $\rho_{\text{out}}$ onto $\ket{k,\sigma-k}$: $\rho_{\text{out}}^{(k,\sigma)}=\bra{k,\sigma-k}\rho_{\text{out}}\ket{k,\sigma-k}$ and renormalize the result, 2. (b)

compute partial transposition of $\rho_{\text{out}}^{(k,\sigma)}$: $\left(\rho_{\text{out}}^{(k,\sigma)}\right)^{\Gamma}$, 3. (c)

compute eigenvalues $\{\alpha_{k}\}$ of $\left(\rho_{\text{out}}^{(k,\sigma)}\right)^{\Gamma}$ using routines built into the NumPy package, 4. (d)

compute the logarithmic negativity

[TABLE]

The advantage of the above algorithm is relatively fast operation at a cost of huge memory requirements because of the large size of the matrices. In the worst case, $\rho_{\text{out}}$ contains $n_{\text{max}}^{8}$ double-precision values (ca. 13MB for $n_{\text{max}}=6$, 800MB for $n_{\text{max}}=10$) and therefore, computations for $n_{\text{max}}$ as large as 20 or 50 would be impossible. This has been solved by noticing that the matrix contains mostly zeros and that the sum of the number of photons in all modes cannot be simultaneously higher than $4n_{\text{max}}$, which allowed us to apply denser packing. Once the final density operator matrix is computed, subsequent computations (step 6) of the algorithm can be performed for all interesting values of $S$ and $k$.