# Long-distance distribution of qubit-qubit entanglement using Gaussian-
correlated photonic beams

J. Agustí1 Y. Minoguchi1 J. M. Fink2 P. Rabl1 1 Vienna Center for Quantum
Science and Technology, Atominstitut, TU Wien, 1020 Vienna, Austria
2 Institute of Science and Technology Austria, am Campus 1, 3400
Klosterneuburg, Austria

###### Abstract

We investigate the deterministic generation and distribution of entanglement
in large quantum networks by driving distant qubits with the output fields of
a non-degenerate parametric amplifier. In this setting, the amplifier produces
a continuous Gaussian two-mode squeezed state, which acts as a quantum-
correlated reservoir for the qubits and relaxes them into a highly entangled
steady state. Here we are interested in the maximal amount of entanglement and
the optimal entanglement generation rates that can be achieved with this
scheme under realistic conditions taking, in particular, the finite amplifier
bandwidth, waveguide losses and propagation delays into account. By combining
exact numerical simulations of the full network with approximate analytic
results, we predict the optimal working point for the amplifier and the
corresponding qubit-qubit entanglement under various conditions. Our findings
show that this passive conversion of Gaussian- into discrete-variable
entanglement offers a robust and experimentally very attractive approach for
operating large optical, microwave or hybrid quantum networks, for which
efficient parametric amplifiers are currently developed.

##  I Introduction

The distribution of entanglement between separated nodes of small- and large-
scale quantum networks is a fundamental task for many quantum communication
and quantum information processing applications Nielsen2000 ; Kimble2008 ;
Wehner2018 ; Kim2014 ; Lukin2007 ; Bianchi2020 . Once established,
entanglement can be locally purified Bennett1996 ; Duer1999 ; Macchiavello1999
; Massar1999 ; Cleland2022  and used for quantum teleportation and remote gate
operation protocols that then require classical communication only Chuang1999
. The availability of a large number of highly entangled qubit pairs, shared
among different nodes, is thus a universal and in practice one of the most
essential resources for quantum network applications.

Existing protocols for distributing entanglement in realistic systems are
mostly based on one of the following two strategies Northup2014 : A first
approach is to generate a pair of entangled qubits locally and then transfer
one of the states to the distant node, for example, through a controlled
emission and reabsorption of optical or microwave photons Cirac1997 ;
Kurpiers2018 ; Schoelkopf2018 ; Wallraff2020 ; Cleland2021 . This protocol is
fully deterministic but requires a sufficiently high level of (synchronous)
control on both nodes of the network. A second strategy is to generate only
weakly correlated qubit-photon pairs from which highly entangled qubit states
can be distilled through measurements and postselection Duan2001 ; Barrett2005
; Matsukevich2006 ; Moehring2007 ; Hensen2015 ; Riedinger2018 . This approach
has the advantage that it requires only little local control and that it is
intrinsically robust with respect to losses. It is, however, only
probabilistic and for many implementations Matsukevich2006 ; Moehring2007 ;
Hensen2015 ; Riedinger2018  the scalability of this approach is limited by the
low success probabilities, which also decrease exponentially with the number
of qubit pairs that must be entangled simultaneously in this way.

Figure 1: Setup. A non-degenerate parametric amplifier is used to produce two
Gaussian-correlated beams of photons, which drive two qubits located in
separate notes of the network into an entangled state. See text for more
details.

A complimentary third strategy for distributing entanglement is to drive
distant qubits with continuous beams of quantum-correlated photons, as shown
in Fig. 1. In the simplest case, entangled beams of optical and/or microwave
fields can be generated in a parametric down-conversion process
Textbooks_CarmichaelOptics , which produces as an output a propagating two-
mode squeezed state. This process only requires a weak intrinsic nonlinearity
as it occurs, for example, in nonlinear optical crystals Zeilinger1995 ;
Wong2006  or in driven Josephson junctions in the microwave regime Eichler2011
; Menzel2012 ; Huard2012 ; Devoret2016 ; Winkel2020 ; Peugeot2021 . Therefore,
parametric down-conversion is currently the most common method to generate
entangled pairs of optical photons Walmsley2011 ; Couteau2018 , however,
usually in a probabilistic and postselected manner. In contrast, here we are
interested in the regime where the parametric amplifier is strongly pumped,
such that the output fields contain many photons on average. By using these
correlated photons to drive two spatially separated qubits, the Gaussian
entanglement can be mapped onto an entangled qubit state. This scheme has the
obvious benefit that it only relies on an externally pumped
$\chi^{(2)}$-nonlinearity for generating the entanglement, which is typically
much easier to realize than strong few-photon interactions or high-fidelity
qubit-qubit gates. At the same time, this approach does not rely on
postselection and can be used to distribute entanglement deterministically.

The basic idea of using two-mode squeezed states of light to drive qubits into
highly entangled states has been originally discussed in a paper by Kraus and
Cirac Kraus2004 , and since then several related schemes for engineering
correlated reservoirs for qubits have been described Paternostro2004 ;
Serafini2006 ; Muschik2011 ; Felicetti2014 ; Blais2017 ; Mirrahimi20018 ;
You2018 ; Govia2021 . Importantly, these works already show that under ideal
conditions a maximally entangled qubit state can be prepared in this way. In
view of the vastly different nature of entangled continuous-variable states
and entangled discrete-variable states, this result is a priori not apparent.
However, all these predictions rely on the crucial assumption that the
correlated photon source can be treated as an effective Markovian reservoir,
which ignores the finite bandwidth of the parametric amplifier, propagation
delays and other effects that are relevant for practical networking
applications. Therefore, a systematic analysis of the performance of this
approach under realistic conditions, for example, in terms of speed,
robustness, or the fidelity of the produced entangled states, is still
wanting.

In this paper, we study the entanglement properties of two distant qubits that
are driven by the output of a realistic two-mode parametric amplifier. In this
analysis we go beyond the usual Markov approximation and use both exact
numerical simulations as well as approximate analytic methods to evaluate the
dynamics and steady states of the qubits, taking the finite bandwidth of the
amplifier and other relevant sources of imperfections fully into account. This
analysis allows us to predict optimal operation points for the parametric
amplifier and assess the maximal fidelities and entanglement generation rates
that can be reached with this scheme under given experimental constraints.
These findings do not rely on any specific implementation and can thus be
applied for the optimization of remote entanglement distribution schemes in
various optical, microwave or hybrid Xiang2013 ; Kurizki2015  quantum
networks.

The remainder of the paper is organized as follows: In Sec. II we introduce
the setup and the individual network components. In Sec. III we first briefly
review the generation of qubit-qubit entanglement in the infinite-bandwidth
limit and discuss the relationship between the achievable amount of
entanglement and the basic characteristics of the squeezed reservoir. In Sec.
IV we go beyond the Markov approximation and use numerical and approximate
analytic results to evaluate the steady-state entanglement in realistic
networks. Finally, in Sec. V and Sec. VI we discuss the performance of the
scheme in asymmetric networks as well as the optimal entanglement distribution
rates in a pulsed operation mode. Finally, in Sec. VII we summarize our main
findings and place them into context with respect to the current experimental
state of the art.

##  II Setup

We consider the generic setting depicted in Fig. 1, where two distant qubits
with frequencies $\omega_{q1}$ and $\omega_{q2}$ are driven by the output
fields of a non-degenerate parametric amplifier. The amplifier consists of two
distinct bosonic modes with frequencies $\omega_{1}$ and $\omega_{2}$ and
annihilation operators $a_{1}$ and $a_{2}$, respectively. These modes are
driven into a correlated two-mode squeezed state via an externally pumped
$\chi^{(2)}$-process, and each of them is connected to one of the qubits
through a unidirectional channel. In practice, such a scenario can be realized
with the help of circulators, directional couplers or chiral waveguides (see,
for example, Refs. Cirac1997 ; Stannigel2014 ; Lodahl2017 ). The following
analysis is kept deliberately general and applies to implementations in the
optical and microwave domain, as well as mixed scenarios, where, for example,
correlated pairs of optical and microwave photons are generated via the
electro-optical effect Rueda2019 ; Krastanov2021 . However, throughout our
analysis we will assume a matching of the corresponding qubit and photon
frequencies, i.e., $\omega_{qi}\approx\omega_{i}$, and that all parts of the
network are sufficiently cold such that thermal excitations can be neglected
Huard2015 .

###  II.1 Qubits

We model the qubits as simple two-level systems with a ground state
$|0\rangle$ and an excited state $|1\rangle$. The qubits decay into the
waveguide with rates $\gamma_{i}$ and undergo dephasing with rate
$\gamma_{\phi}$. In the absence of the amplifier and by changing into a
rotating frame with respect to the qubit frequencies $\omega_{q1}$ and
$\omega_{q2}$, the resulting dynamics for the qubit density operator
$\rho_{q}$ is given by the master equation

$\dot{\rho}_{q}=\sum_{i=1,2}\gamma_{i}\mathcal{D}[\sigma^{-}_{i},\sigma_{i}^{+}]\rho_{q}+\frac{\gamma_{\phi}}{2}\mathcal{D}[\sigma^{z}_{i},\sigma^{z}_{i}]\rho_{q}\equiv\mathcal{L}_{q}\rho_{q}.$
(1)

Here, $\sigma^{-}=(\sigma^{+})^{\dagger}=|0\rangle\langle 1|$ and
$\sigma^{z}=|1\rangle\langle 1|-|0\rangle\langle 0|$ and we have introduced
the short notation

$\mathcal{D}[A,B]\rho=A\rho B-\frac{1}{2}\left(BA\rho+\rho BA\right)$ (2)

for the Liouville superoperators describing incoherent decay and dephasing.

###  II.2 Parametric amplifier

Parametric amplification relies on a $\chi^{(2)}$-type process between three
bosonic modes $(\hbar=1)$ Textbooks_CarmichaelOptics ,

$H_{\rm\chi}=i\chi(a_{1}^{\dagger}a_{2}^{\dagger}a_{0}-a_{1}a_{2}a_{0}^{\dagger}),$
(3)

where the nonlinearity $\chi$ is small, but one of the modes, $a_{0}$, is
strongly pumped. Under this assumption, the pumped mode can be treated as a
classical field, $a_{0}\rightarrow\alpha_{0}\in\mathbbm{C}$, and $H_{\rm\chi}$
reduces to a two-mode squeezing interaction. This process becomes most
effective when the resonance condition
$\omega_{1}+\omega_{2}\approx\omega_{0}$ between the three modes is fulfilled
(see Fig. 1) and is frequently employed in nonlinear optical crystals to
produce entangled photon pairs. Similar interactions also occur in
superconducting circuits, where driven Josephson junctions are used to
generate two-mode squeezed microwave beams, and in many other devices. In the
following, we fix the pump frequency to a value of
$\omega_{0}=\omega_{q1}+\omega_{q2}$, in order to maximize the resulting
qubit-qubit correlations discussed below. The residual photon-qubit detunings
are denoted by $\Delta_{i}=\omega_{i}-\omega_{qi}$.

The parametrically generated pairs of photons in modes $a_{1}$ and $a_{2}$
decay into the waveguides with rates $\kappa_{1}$ and $\kappa_{2}$,
respectively. By changing into a rotating frame with respect to $\omega_{q1}$
and $\omega_{q2}$, the full dynamics of the two-mode photonic state $\rho_{p}$
is described by the master equation

$\dot{\rho}_{p}=-i\left[H_{p},\rho_{p}\right]+\sum_{i=1,2}\kappa_{i}\mathcal{D}[a_{i},a_{i}^{\dagger}]\rho_{p}\equiv\mathcal{L}_{p}\rho_{p},$
(4)

where

$H_{p}=\sum_{i=1,2}\Delta_{i}a_{i}^{\dagger}a_{i}+i\frac{\sqrt{\kappa_{1}\kappa_{2}}\epsilon}{2}\left(a_{1}^{\dagger}a_{2}^{\dagger}-a_{1}a_{2}\right)$
(5)

and $\epsilon\sim\chi|\alpha_{0}|$ is the dimensionless pump parameter, which
we can assume to be real. For $\Delta_{i}=0$, the value of $\epsilon=1$ marks
the onset of the parametric instability, beyond which our linearized
description of the amplifier is no longer valid. For the remainder of this
work, we focus mainly on the resonant case and restrict the pumping parameter
to $\epsilon\in[0,1)$.

###  II.3 Cascaded photon-qubit interactions

In the considered cascaded setting, photons emitted by the parametric
amplifier will drive the qubits located further along the waveguides, while
there is no backaction of the qubits on the photonic modes. By assuming a
waveguide with a broad and linear dispersion relation, this scenario can be
modelled in terms of the cascaded master equation Carmichael1992 ;
Gardiner1992

$\dot{\rho}=\left(\mathcal{L}_{q}+\mathcal{L}_{p}+\mathcal{L}_{\rm
cas}\right)\rho,$ (6)

where $\rho$ is the density operator of the whole system, including both the
qubits and the photonic degrees of freedom. The last term in Eq. (6),

$\mathcal{L}_{\rm
cas}\rho=\sum_{i=1,2}\sqrt{\eta\gamma_{i}\kappa_{i}}\left([a_{i}\rho,\sigma_{i}^{+}]+[\sigma^{-}_{i},\rho
a_{i}^{\dagger}]\right),$ (7)

accounts for the cascaded, i.e. unidirectional coupling between the photons
and the qubits. Here we have introduced the additional phenomenological
parameter $\eta\in[0,1]$. A value of $\eta<1$ can either be used to model
losses along the waveguide, in which case $(1-\eta)$ is the probability that a
propagating photon is lost between the amplifier and the qubit, or additional
decay channels for the qubits other than the emission into the waveguide.

Note that in writing Eq. (6) we have already reabsorbed all propagation phases
into a redefinition of the qubit states. For now, we have also neglected the
finite propagation times $\tau_{i}=d_{i}/c_{i}$, which it takes for photons
with group velocity $c_{i}$ to travel the distance $d_{i}$ between the
amplifier and the $i$-th qubit. This is usually a valid assumption for small
on-chip networks, but propagation delays can become very relevant when
discussing entanglement distribution in larger networks, in particular when
$d_{1}\neq d_{2}$. This issue will be addressed in more detail in Sec. V
below.

###  II.4 Reduced qubit state

Starting from the cascaded master equation of the full network, Eq. (6), our
main goal in the following is to study the dynamics and stationary states of
the reduced qubit state, $\rho_{q}(t)={\rm Tr}_{p}\\{\rho(t)\\}$, as a
function of the pumping strength $\epsilon$ and other system parameters. From
the structure of the superoperator $\mathcal{L}_{\rm cas}$ it is
straightforward to show that the cascaded coupling doesn’t affect the dynamics
of the reduced state of the amplifier, i.e., $\dot{\rho}_{p}={\rm
Tr}_{q}\\{\mathcal{L}_{\rm cas}\rho\\}=0$. However, this doesn’t imply that
the system can be factorized into a photonic and a qubit part. In contrast, it
can be explicitly shown that in such cascaded systems a large amount of
entanglement between the individual subsystems can emerge Stannigel2014 .
Therefore, in general, it is necessary to solve the full master equation in
order to obtain accurate predictions for $\rho_{q}$.

##  III Qubits coupled to a two-mode squeezed reservoir

As a starting point, it is instructive to revisit the idealized limit of a
broadband amplifier, $\kappa_{i}\rightarrow\infty$, where the dynamics of the
amplifier modes is much faster than the qubit evolution and can be
adiabatically eliminated. In this case, it is indeed possible to obtain a
reduced master equation for the qubit state $\rho_{q}$ only, which reads

$\begin{split}\dot{\rho}_{q}&=\mathcal{L}_{q}\rho_{q}+\sum_{i=1,2}\sum_{s=\pm}\gamma_{i}N_{i}\mathcal{D}[\sigma_{i}^{s},\sigma_{i}^{-s}]\rho_{q}\\\
&-\sqrt{\gamma_{1}\gamma_{2}}\left\\{M\left(\mathcal{D}[\sigma_{1}^{+},\sigma_{2}^{+}]\rho_{q}+\mathcal{D}[\sigma_{2}^{+},\sigma_{1}^{+}]\rho_{q}\right)+\rm
H.c.\right\\}.\end{split}$ (8)

Here the occupation numbers
$N_{i}=2\eta\mathrm{Re}\\{I_{a_{i}^{\dagger}a_{i}}(0)\\}$ and the squeezing
term $M=\eta[I_{a_{1}a_{2}}(0)+I_{a_{2}a_{1}}(0)]$ are determined by the
steady-state correlation spectra of the amplifier modes,
$I_{A_{i}B_{j}}(\omega)=\sqrt{\kappa_{i}\kappa_{j}}\int_{0}^{\infty}\mathrm{d}t\,\braket{A_{i}(t)B_{j}(0)}e^{-i\omega
t}$, evaluated at frequency $\omega=0$. A detailed derivation of this master
equation and explicit expressions for the two-time correlation functions are
given in Appendix A and Appendix B, respectively.

In the symmetric case, $\kappa_{i}=\kappa$, and for $\Delta_{i}=0$ we obtain
the simple results

$\displaystyle
2\mathrm{Re}\\{I_{a_{i}^{\dagger}a_{i}}(\omega)\\}=\epsilon\left[\Gamma_{-}(\omega)-\Gamma_{+}(\omega)\right],$
(9)

and

$\displaystyle
I_{a_{1}a_{2}}(\omega)+I_{a_{2}a_{1}}(-\omega)=\epsilon\left[\Gamma_{-}(\omega)+\Gamma_{+}(\omega)\right],$
(10)

where

$\Gamma_{\pm}(\omega)=\frac{\kappa^{2}}{\kappa^{2}(1\pm\epsilon)^{2}+4\omega^{2}}$
(11)

are Lorentzian functions. The parameters $N=N_{i}$ and $M$ in Eq. (8) are then
given by

$\displaystyle N$ $\displaystyle=$
$\displaystyle\epsilon\eta\left(\frac{1}{\left(1-\epsilon\right)^{2}}-\frac{1}{\left(1+\epsilon\right)^{2}}\right),$
(12) $\displaystyle M$ $\displaystyle=$
$\displaystyle\epsilon\eta\left(\frac{1}{\left(1-\epsilon\right)^{2}}+\frac{1}{\left(1+\epsilon\right)^{2}}\right).$
(13)

###  III.1 Squeezing and purity

To obtain additional intuition, we focus on the symmetric configuration
$\kappa_{i}=\kappa$, where we can reinterpret $N=\langle
a^{\dagger}_{i}a_{i}\rangle_{\rho_{\rm eff}}$ and $M=\langle
a_{1}a_{2}\rangle_{\rho_{\rm eff}}$ as expectation values of an effective two-
mode squeezed state

$\rho_{\rm eff}=S(r_{\rm eff})\rho_{\rm th}(\overline{n}_{\rm
eff})S^{\dagger}(r_{\rm eff}).$ (14)

Here, $\rho_{\rm th}(\overline{n}_{\rm eff})$ is a thermal two-mode state with
occupation number $\overline{n}_{\rm eff}$ and

$S(r)=e^{r(a_{1}^{\dagger}a_{2}^{\dagger}-a_{1}a_{2})}$ (15)

is the two-mode squeezing operator. Instead of the thermal occupation number,
we can also work with the purity of the effective state, $\mu_{\rm
eff}=1/(2\overline{n}_{\rm eff}+1)$, in which case we obtain the relations
Serafini04_TwoModeGaussianNoise

$\displaystyle r_{\rm eff}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\tanh^{-1}{\left[\frac{2|M|}{2N+1}\right]},$ (16)
$\displaystyle\mu_{\rm eff}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{(2N+1)^{2}-4|M|^{2}}}.$ (17)

From the expressions given in Eq. (12) and Eq. (13) we find that
$|M|^{2}=N(N+\eta)$ and therefore $\mu_{\rm eff}=1$ for a lossless channel.
This implies that in the infinite-bandwidth limit, Eq. (8) describes two
qubits that are coupled to a two-mode squeezed zero-temperature reservoir with
a squeezing parameter

$r_{\rm eff}(\epsilon)=2\tanh^{-1}(\epsilon),$ (18)

which becomes arbitrarily large when $\epsilon\rightarrow 1$. However, this is
a consequence of the Markov approximation, where the environment is only
probed at a single frequency. Below we will show that this is no longer true
when a nonvanishing ratio $\kappa/\gamma$ is taken into account, in which case
averaging over a finite frequency window $\Delta\omega\neq 0$ translates into
an impure effective state with reduced squeezing. The same is already true in
the Markovian limit in the presence of transmission losses, $\eta<1$.

###  III.2 Steady state entanglement

From the reduced master equation in Eq. (8) we can derive the steady state of
the two qubits, which we express in general as

$\rho_{q}^{0}=\sum_{\boldsymbol{s},\boldsymbol{s}^{\prime}}\rho^{0}_{\boldsymbol{s},\boldsymbol{s}^{\prime}}|\boldsymbol{s}\rangle\langle\boldsymbol{s}^{\prime}|,$
(19)

where $\boldsymbol{s}=(s_{1},s_{2})$ with $s_{i}=0,1$ labels the two-qubit
states. For a symmetric parametric amplifier, $\kappa_{1}=\kappa_{2}$, and for
symmetric qubits, $\gamma_{1}=\gamma_{2}$, there are only six nonvanishing
matrix elements,

$\displaystyle\Lambda\rho^{0}_{00,00}$ $\displaystyle=$
$\displaystyle(1+N)^{2}(2\Gamma_{\phi}+1+2N)-|M|^{2}(3+2N),$
$\displaystyle\Lambda\rho^{0}_{10,10}$ $\displaystyle=$ $\displaystyle
2\Gamma_{\phi}N(N+1)+(2N+1)(N(N+1)-|M|^{2}),$
$\displaystyle\Lambda\rho^{0}_{11,11}$ $\displaystyle=$ $\displaystyle
N^{2}(1+2\Gamma_{\phi}+2N)+|M|^{2}(1-2N),$
$\displaystyle\Lambda\rho^{0}_{11,00}$ $\displaystyle=$ $\displaystyle M,$

and $\rho^{0}_{01,01}=\rho^{0}_{10,10}$ and
$\rho^{0}_{00,11}=(\rho^{0}_{11,00})^{*}$. Here we introduced the
normalization constant
$\Lambda=(1+2N)[1+2\Gamma_{\phi}+4(N(N+1+\Gamma_{\phi})-|M|^{2}))]$ and the
normalized dephasing rate $\Gamma_{\phi}=\gamma_{\phi}/\gamma$. In the Markov
limit and assuming also otherwise ideal conditions, $\gamma_{\phi}=0$ and
$\eta=1$, we have $|M|^{2}=N(N+1)$, $\rho^{0}_{10,10}=0$ and
$|\rho^{0}_{11,00}|^{2}=\rho^{0}_{11,11}\rho^{0}_{00,00}$. This means that the
steady state is the pure state $\rho_{q}^{0}=\ket{\Psi_{q}}\bra{\Psi_{q}}$,
where Cirac1997

$\ket{\Psi_{q}}=\sqrt{\frac{N+1}{1+2N}}\ket{00}+e^{i\theta}\sqrt{\frac{N}{1+2N}}\ket{11},$
(20)

and we defined $M=|M|e^{i\theta}$. On resonance, $\Delta_{i}=0$, we obtain
$\theta=0$ and the steady state approaches the maximally entangled triplet
state $\ket{\Phi^{+}}=\frac{1}{\sqrt{2}}(\ket{00}+\ket{11})$ for $N\gg 1$.

Before we proceed, let us provide some additional insights about the emergence
of such a pure entangled steady state, by considering the coupling of the
qubits to two isolated modes $a_{1}$ and $a_{2}$ via a Jaynes-Cummings
interaction of the form

$H_{\rm int}\sim
i\left(\sigma_{1}^{-}a_{1}^{\dagger}-\sigma_{1}^{+}a_{1}+\sigma_{2}^{-}a_{2}^{\dagger}-\sigma_{2}^{+}a_{2}\right).$
(21)

An ideal two-mode squeezed state of those two modes can be written as

$\ket{\Psi_{\rm
TMS}}=\frac{1}{\sqrt{1-x^{2}}}\sum_{n=0}^{\infty}x^{n}|n\rangle_{1}|n\rangle_{2},$
(22)

where we can set $x=\sqrt{N/(N+1)}$ to match the convention from above. This
expression shows that the number of photons in the two modes are perfectly
correlated, suggesting that the qubits are only excited and de-excited
pairwise. However, this argument is too naive, since, for example, the action
of $H_{\rm int}$ on the state $|00\rangle\ket{\Psi_{\rm TMS}}$ would also
generate singly-excited states $\sim|10\rangle,|01\rangle$. To explain the
existence of the steady state given in Eq. (20) it is thus important to take
into account the coherence between the $|n\rangle_{1}|n\rangle_{2}$
components, which leads to the following relations

$\displaystyle a_{1}\ket{\Psi_{\rm TMS}}$ $\displaystyle=$ $\displaystyle
xa_{2}^{\dagger}\ket{\Psi_{\rm TMS}},$ (23) $\displaystyle a_{2}\ket{\Psi_{\rm
TMS}}$ $\displaystyle=$ $\displaystyle xa_{1}^{\dagger}\ket{\Psi_{\rm TMS}}.$
(24)

Equivalently, there exists a unique dark state of the interaction,

$H_{\rm int}(|00\rangle+x|11\rangle)\ket{\Psi_{\rm TMS}}=0.$ (25)

Therefore, once the system reaches the state $\ket{\Psi_{q}}$ in Eq. (20), the
emission of a photon by one qubit interferes destructively with the absorption
of a photon in the other mode. When applied to the original setting, this
argument shows that the absence of any components $\sim|01\rangle$ or
$\sim|10\rangle$ in $|\Psi_{q}\rangle$ is a consequence of a nonlocal
interference effect between two well-separated, but correlated parts of the
network.

###  III.3 Non-ideal squeezed reservoirs

To investigate the performance of this entanglement distribution scheme also
under non-ideal conditions, for example, when $|M|^{2}<N(N+1)$ or when a
reduced master equation for $\rho_{q}$ is not available, we consider the
fidelity

$\mathcal{F}=|\braket{\Phi^{+}}{\rho^{0}_{q}}{\Phi^{+}}|^{2}\leq\frac{(1+\epsilon)^{4}}{2(1+6\epsilon^{2}+\epsilon^{4})}$
(26)

as a measure of how close the qubit state approaches the maximally entangled
triplet state. In addition, we can use the concurrence of the reduced state,
$\mathcal{C}(\rho_{q}^{0})$ Horodeki2009 , to quantify directly the amount of
qubit-qubit entanglement. Note that the upper bound in Eq. (26) derived for
the ideal Markovian limit shows that already for moderate driving strengths,
$\epsilon\gtrsim 0.5$, fidelities of about $\mathcal{F}\approx 0.99$ can be
reached.

In Fig. 2 we show a plot of the steady state fidelity $\mathcal{F}$ for the
steady state of master equation (8), under the assumption that
$\Gamma_{\phi}=0$, but allowing for arbitrary values of $N_{i}=N$ and
$|M|^{2}\leq N(N+1)$. For later convenience, these parameters are in turn
expressed in terms of the effective squeezing parameter $r_{\rm eff}$ and the
effective purity $\mu_{\rm eff}$, as defined in Eq. (16) and Eq. (17). We see
that for $\mu_{\rm eff}\simeq 1$, the fidelity approaches unity as the
squeezing parameter is increased, consistent with the bound stated in Eq.
(26). For values of $r_{\rm eff}\gtrsim 1$, the fidelity becomes almost
independent of the squeezing parameter, but decreases as
$\mathcal{F}\simeq\frac{1}{4}(1+3\mu^{2}_{\rm eff})$ for impure reservoirs.
Entanglement is present as long as $\mathcal{F}>0.5$, which corresponds to a
minimal purity of about $\mu_{\rm eff}\simeq 1/\sqrt{3}\simeq 0.57$. This
result shows that improving the purity of the effective photonic bath will be
most relevant for this entanglement distribution scheme.

Figure 2: Contour plot of the steady state fidelity $\mathcal{F}$ defined in
Eq. (26) as a function of the effective squeezing parameter $r_{\rm eff}$ and
the effective purity $\mu_{\rm eff}$. The solid lines indicate the path in
this parameter space that one obtains by increasing $\epsilon$ from 0 to 1 for
different values of $\beta=\kappa/\gamma$ and $\Gamma_{\phi}=0$ and $\eta=1$
(see discussion in Sec. IV). The dashed line marks the boundary of vanishing
concurrence $\mathcal{C}=0$, above which the reduced qubit state is entangled.

##  IV Entanglement distribution in realistic networks

Let us now address the more realistic scenario, where finite waveguide losses,
decoherence of the qubits and, in particular, the finite bandwidth of the
parametric amplifier are taken into account. In this case, it is in general no
longer possible to eliminate the photon modes and the full cascaded master
equation in Eq. (6) must be solved numerically. In view of the large Hilbert
space required to represent the two-mode squeezed state, these simulations
become very demanding when approaching the parametric instability,
$\epsilon=1$. Thus, in the following, we restrict all numerical simulations to
values of $\epsilon\leq 0.8$, where a convergence of the results can still be
ensured by truncating the Hilbert space of each photon mode to $n_{\rm
trunc}\lesssim 25$ states and by performing finite size scaling analysis (see
Appendix C). In addition, we introduce a filtered-mode approximation (FMA) in
order to obtain also an analytic dependence of the fidelity on all the
relevant system parameters.

###  IV.1 Filtered mode approximation

In the Markov limit studied in the previous section, the state of the qubits
obeys a master equation, where the characteristic bath parameters $N_{i}$ and
$M$ are determined by the output fields of the amplifier at a single
frequency, $\omega=0$. To go beyond this approximation, we must take into
account that the qubits will be affected by photons within a finite region of
the spectra $S_{a_{i}^{\dagger}a_{j}}(\omega)$ and $S_{a_{i}a_{j}}(\omega)$
that cannot be associated with a pure squeezed state. The relevant bandwidth
of frequencies will be determined by the dynamics of the qubits themselves and
will be roughly given by the decay rates $\gamma_{i}$. Based on this
intuition, we introduce the two filtered modes

$a_{f,i}(t)=\sqrt{\gamma_{i}\kappa_{i}\eta}\int_{-\infty}^{t}{\rm
d}s\,e^{-\gamma_{i}(t-s)/2}a_{i}(s-\tau_{i}),$ (27)

where, for a later generalization, we have already included the propagation
delays $\tau_{i}$. These modes represent the output of the two-mode amplifier,
but delayed by $\tau_{i}$ and filtered by the response of the qubits.

We can now use these filtered modes to define an adjusted set of parameters
for the qubit master equation in Eq. (8),

$N_{i}=\langle a^{\dagger}_{f,i}a_{f,i}\rangle,\qquad
M=\langle\mathcal{T}a_{f,1}a_{f,2}\rangle,$ (28)

where $\mathcal{T}$ denotes the time-ordering operator applied to the
amplifier modes $a_{1}$ and $a_{2}$. These parameters include the
characteristic timescales of the qubits and of the photons on an equal
footing. Specifically, we obtain

$N_{i}=2\eta\gamma_{i}\int_{-\infty}^{\infty}\frac{\mathrm{d}\omega}{2\pi}\frac{I_{a_{i}^{\dagger}a_{i}}(\omega)}{\gamma_{i}^{2}/4+\omega^{2}}$
(29)

for the occupation numbers and

$\begin{split}M=\sqrt{\gamma_{1}\gamma_{2}}\eta\int_{-\infty}^{\infty}\frac{\mathrm{d}\omega}{2\pi}\frac{\left[I_{a_{1}a_{2}}(\omega)+I_{a_{2}a_{1}}(-\omega)\right]e^{i\omega(\tau_{2}-\tau_{1})}}{(\gamma_{1}/2+i\omega)(\gamma_{2}/2-i\omega)}\end{split}$
(30)

for the correlation parameter. These expressions explicitly show how the
effective photonic reservoir seen by the qubits depends on the amplifier
correlations within finite frequency windows set by the decay rates
$\gamma_{i}$. The Markovian limit discussed in Sec. III is recovered in the
limit $\gamma_{i}\rightarrow 0$ and $\tau_{i}\rightarrow 0$.

###  IV.2 Effective squeezing parameters for non-ideal amplifiers

While the derivation of a master equation in terms of the filtered modes is
only an approximation, it becomes exact in the regime of low qubit
excitations, i.e. for small values of $\epsilon$ (see Appendix D). Even for
moderate and large driving strengths it still results in a considerable
improvement over the conventional master equation discussed in Sec. III. In
particular, the FMA allows us to account for the effects of a finite amplifier
bandwidth. For example, by setting $\kappa_{i}=\kappa$, $\gamma_{i}=\gamma$
and assuming $\tau_{i}=0$ for now, we obtain

$\displaystyle N$ $\displaystyle=$
$\displaystyle\frac{2\epsilon^{2}\beta(1+2\beta)\eta}{[(\beta+1)^{2}-\beta^{2}\epsilon^{2}](1-\epsilon^{2})},$
(31) $\displaystyle M$ $\displaystyle=$
$\displaystyle\frac{2\epsilon\beta(\epsilon^{2}\beta+\beta+1)\eta}{[(\beta+1)^{2}-\beta^{2}\epsilon^{2}](1-\epsilon^{2})},$
(32)

where $\beta=\kappa/\gamma$ is the ratio between the amplifier bandwidth and
the qubit decay.

Figure 3: Dependence of the effective squeezing parameter and the effective
purity on the driving strength within the FMA. The different curves are
evaluated for different values of $\beta$ and assuming $\eta=1$ and symmetric
conditions, $\kappa_{i}=\kappa$ and $\gamma_{i}=\gamma$. On the right axis, we
plot the squeezing factor $S_{\rm eff}=e^{2r_{\rm eff}}$ in $\mathrm{dB}$.
Note that $r_{\rm eff}\rightarrow\infty$ and $\mu_{\rm eff}\rightarrow 0$ when
$\epsilon\rightarrow 1$ for all values of $\beta$.

As discussed in Sec. III.1 and shown explicitly in Fig. 3, these parameters
can be reexpressed in terms of an effective squeezing parameter $r_{\rm eff}$
and an effective purity $\mu_{\rm eff}$. In this way, the fidelity of the
resulting steady state can be read off directly from the general plot in Fig.
2.

We may also use the analytic expressions for $N$ and $M$ from above to
evaluate the dependence of these quantities up to the first order corrections
in the inverse bandwidth ratio $1/\beta$,

$\displaystyle r_{\rm eff}$ $\displaystyle\simeq$ $\displaystyle
2\tanh^{-1}(\epsilon)-\frac{1}{\beta}\frac{2\epsilon}{(1-\epsilon^{2})^{2}},$
(33) $\displaystyle\mu_{\rm eff}$ $\displaystyle\simeq$ $\displaystyle
1-\frac{1}{\beta}\frac{4\epsilon^{2}}{(1-\epsilon^{2})^{2}}.$ (34)

which is valid for $\beta(1-\epsilon)^{2}\gg 1$. We see that finite-bandwidth
corrections get strongly amplified as one approaches the parametric
instability. In particular, the purity of the effective squeezed reservoir,
which is most relevant for the entanglement of the reduced qubit state,
decreases significantly. Thus, even for $\beta\gg 1$, it is not possible to
assess the achievable amount of entanglement using a purely Markovian
description.

###  IV.3 Optimal fidelities

Figure 4: (a) Plot of the Bell-state fidelity $\mathcal{F}$ as a function of
the driving strength and for different amplifier bandwidths,
$\beta=\kappa/\gamma$. The solid lines represent the results obtained from the
numerical solution of the full cascaded master equation, Eq. (6), which are
compared with the predictions under the FMA (dashed lines), assuming
$\Gamma_{\phi}=0$ and $\eta=1$. (b) Plot of the optimal fidelity
$\mathcal{F}_{\rm op}$ and the corresponding optimal driving strength
$\epsilon_{\rm op}$ (inset) under the same conditions. The dotted line shows
the analytic approximation given in Eq. (37). (c) Plot of the optimal fidelity
as a function of the qubit dephasing $\Gamma_{\phi}$ for $\eta=1$ and (d) as a
function of the channel transmissivity $\eta$ for $\Gamma_{\phi}=0$. In both
plots, we assume a value of $\beta=10^{2}$ for the upper curves and $\beta=1$
for the lower curves, where the dashed lines represent the respective FMA
results. For all plots in this figure we have set $\Delta_{i}=0$,
$\kappa_{1,2}=\kappa$ and $\gamma_{1,2}=\gamma$.

Figure 4 summarizes the performance of the entanglement distribution scheme in
realistic settings. First of all, Fig. 4 (a) shows the dependence of the
fidelity $\mathcal{F}$ on the driving strength $\epsilon$ for different ratios
$\beta=\kappa/\gamma$. Here we compare the results from a simulation of the
full master equation with the predictions obtained from the FMA. In both
cases, we find the expected maximum for intermediate values of $\epsilon$,
which results from an increase in squeezing on the one hand and from the loss
of purity on the other hand. While for moderate and large driving strengths we
see a deviation of the approximate results from exact numerics, the
qualitative trends are still accurately captured. Importantly, in all the
investigated parameter regimes we find that the FMA agrees with or
underestimates the exact fidelity and can thus be reliably used to predict
lower bounds for the achievable amount of entanglement.

In the limit of low pump values, the FMA becomes exact and we obtain a simple
analytic expression for the fidelity,

$\mathcal{F}(\epsilon\ll
1)\simeq\frac{1}{2}+\frac{2\beta\eta}{(1+\beta)(1+2\Gamma_{\phi})}\epsilon.$
(35)

Therefore, the fidelity increases linearly with the pump strength [see Fig. 4
(a)], with a slope that depends on all the different sources of imperfections.
To estimate the maximally achievable fidelities, we assume that this maximum
is reached for a pumping strength $\epsilon\approx 1$ and expand $\mathcal{F}$
to lowest order in $(1-\epsilon)$, $1/\beta$, $\Gamma_{\phi}$ and $(1-\eta)$,

$\mathcal{F}(\epsilon\approx 1)\simeq
1-\frac{(1-\epsilon)^{4}}{16}-\frac{3}{(1-\epsilon)^{2}}\left[\frac{1}{2\beta}+\Gamma_{\phi}+(1-\eta)\right].$
(36)

By optimizing this result with respect to the driving strength we obtain

$\mathcal{F}^{\rm app}_{\rm op}\simeq
1-\frac{3\sqrt[3]{9}}{4}\left[\frac{1}{2\beta}+\Gamma_{\phi}+(1-\eta)\right]^{\frac{2}{3}}.$
(37)

Although the result in Eq. (37) is based on various crude approximations, it
still gives a good estimate for the overall scaling of the maximal fidelity
that is achievable with this scheme in the presence of imperfections. In
particular, as shown in Fig. 4 (b), for the parameter regimes of interest,
$\mathcal{F}_{\rm op}>\mathcal{F}^{\rm app}_{\rm op}$, where $\mathcal{F}_{\rm
op}$ is the exact optimized fidelity evaluated numerically. In Figs. 4 (c) and
(d) $\mathcal{F}_{\rm op}$ is also shown for different non-ideal settings.

It is worth emphasizing that a fidelity of $\mathcal{F}>0.5$, and therefore a
finite amount of steady-state entanglement, can be reached with this scheme
even for rather high waveguide losses. This is in contrast to a deterministic
state transfer scheme, where losses of $(1-\eta)\geq 0.5$ would not permit a
finite amount of entanglement.

###  IV.4 The role of nonlinearities

All the results presented in this work are based on a linearized treatment of
the parametric amplifier. For driving strengths up to $\epsilon=0.8$, the
average photon number in each mode, $\langle
a_{i}^{\dagger}a_{i}\rangle\lesssim 1$, is still rather low and a linearized
description should be justified in most systems of interest. However, to
estimate the impact of a finite nonlinearity in the system, we can study the
steady state fidelity $\mathcal{F}$ as a function of the Hilbertspace
truncation number $n_{\rm trunc}$ in our numerical simulation. Such an
investigation is presented in Appendix C and shows that for small and moderate
$\epsilon\lesssim 0.5$, the amplifier must be linear only up to photon numbers
of $n_{\rm trunc}\approx 4-6$ in order to reach the same optimal fidelities as
in the fully linear case. As we approach the instability, the requirements on
the linearity of the amplifier become more stringent, but only when the
bandwidth is sufficiently high. For example, $n_{\rm trunc}\approx 20$ is
sufficient to reach $\mathcal{F}>0.98$ for $\epsilon=0.8$ and $\beta=10^{3}$.

##  V Steady state entanglement in asymmetric networks

In our discussion so far we have neglected the finite time that it takes for
the photons to propagate from the down-conversion source to the qubits. This
is not a crucial assumption in situations where both qubits are located
approximately equally far away from the amplifier. The photons then simply
take a time $\tau_{1}\approx\tau_{2}$ to propagate from the source to the
qubits, which is irrelevant for the steady state correlations. However, in
situations where $\tau_{1}\neq\tau_{2}$, the qubits are driven by photons that
have been emitted at two different times. When the time lag
$\tau=\tau_{2}-\tau_{1}$ is too long, correlations between these photons are
lost Fedorov2018 . Such situations can occur, for example, when the amplifier
and the first qubit are located in the same laboratory, while the second beam
is sent via an optical fibre to a qubit at a remote location.

###  V.1 Quantum networks with propagation delays

For larger networks with non-negligible propagation delays, different choices
for the definition of entanglement and correlations can be considered. For
quantum key distribution schemes or similar applications, where the quantum
states are only used once, one is typically interested in correlations between
measurements that are delayed by the respective propagation times of the
transmitting photons. However, for other applications, where quantum states
are redistributed within the network multiple times, the more relevant
question is how much entanglement there exists between different nodes at a
given point in time. In the following, we are interested in this second type
of scenario, where signal delays become relevant.

In the presence of finite propagation delays, the dynamics of the full density
operator $\rho(t)$ can no longer be described by the time-local master
equation in Eq. (6), which was the basis for our exact numerical simulations
so far. However, as discussed in more detail in Appendix D, this master
equation can still be used to evaluate non-equal time correlation functions,
which can be related to the steady state expectation values of the actual
network with time delays. More precisely, given an arbitrary product of two
qubit operators $O_{1}\equiv O_{1}\otimes\mathbbm{1}_{2}$ and
$O_{2}\equiv\mathbbm{1}_{1}\otimes O_{2}$, its steady state expectation value
can be computed as

$\langle O_{1}O_{2}\rangle_{0}=\langle
O_{1}(\tau_{2}-\tau_{1})O_{2}\rangle_{0}\Big{|}_{\rm loc},$ (38)

assuming that $\tau_{2}>\tau_{1}$. Here, $\left.\langle
O_{1}(\tau)O_{2}\rangle_{0}\right|_{\rm loc}$ denotes a non-equal time
correlation function, which is evaluated with the help of the quantum
regression theorem using the time-local cascaded master equation in Eq. (6).

Based on this relation, we express the reduced steady state of the two qubits
as

$\rho_{q}^{0}=\sum_{\mu,\nu=0,x,y,z}^{3}\braket{\sigma_{1}^{\mu}\sigma_{2}^{\nu}}_{0}\sigma_{1}^{\mu}\sigma_{2}^{\nu},$
(39)

where $\sigma^{0}\equiv\mathbbm{1}$. In this way, we can employ Eq. (38) to
evaluate the full two-qubit density matrix of a time-delay network through
numerical simulations of the time-local master equation given in Eq. (6). In
addition, we can again use the FMA to derive a time-local master equation for
$\rho_{q}$ only. In this approach all the propagation delays are already
included in the parameters $N_{i}$ and $M$ in Eq. (29) and Eq. (30), which are
derived from the retarded field operators in Eq. (27).

###  V.2 Fidelities and entanglement in the presence of propagation delays

Figure 5: (a) Time delayed effective parameters $r_{\rm eff}(\tau)$ and (b)
$\mu_{\rm eff}(\tau)$ for $\epsilon=0.4$ (solid) and $\epsilon=0.9$ (dotted)
and for different amplifier bandwidths $\beta=1$ and $\beta=10^{3}$. (c)
Steady state concurrence of a network with time delays $\mathcal{C}(\tau)$ for
$\epsilon=0.5$ and different values of $\beta$. The solid lines represent the
results obtained from a full numerical simulation based on Eq. (39), while the
dashed lines indicate the FMA predictions. (d) Plot of the entanglement time
$\tau_{\rm ent}$, i.e., the minimal delay time beyond which the entanglement
vanishes, $\mathcal{C}(\tau_{\rm ent})=0$. In all plots, $\eta=1$,
$\Gamma_{\phi}=0$ and a symmetric setup with $\Delta_{i}=0$ have been assumed.

In Fig. 5 we consider a network with different propagation times
$\tau_{1}\neq\tau_{2}$ for the two beams and study the (equal-time) steady
state entanglement as a function of $\tau=\tau_{2}-\tau_{1}$. In Fig. 5 (a)
and (b) we plot, first of all, the effective squeezing $r_{\rm eff}(\tau)$ and
the purity $\mu_{\rm eff}(\tau)$, as obtained from the FMA for different
values of $\beta$ and $\epsilon$. Within this approximation and in the limit
of $\beta\gg 1$ we find that

$M(\tau)\simeq\frac{2\epsilon(1+\epsilon^{2})}{(1-\epsilon^{2})^{2}}e^{-\gamma\tau/2}.$
(40)

This shows that in this limit the correlations decay on the timescale set by
$\gamma$, and not by the relaxation rates of the amplifier,
$\kappa_{\pm}=\kappa(1\pm\epsilon)$. Therefore, a finite amount of squeezing
prevails up to delay times of about $\tau\approx\gamma^{-1}$ for all driving
strengths. In contrast, the behaviour of the effective purity depends more
strongly on the driving strength. This is because a larger value of
$M(\tau=0)$ implies a larger absolute change of $M(\tau)$, as relevant for the
purity $\mu_{\rm eff}$ defined in Eq. (17). The timescale that determines the
decay of the purity, and therefore the entanglement, can become considerably
shorter than $\gamma^{-1}$ for large driving strength. This conclusion is
consistent with the entanglement decay in delayed two-mode squeezed states
reported in Ref. Fedorov2018 . In the opposite regime $\beta\lesssim 1$, the
squeezing and purity parameters are smaller to begin with, but they are more
robust and decay only after a delay $\tau>\gamma^{-1}$, roughly independently
of $\epsilon$.

In Fig. 5 (c) we evaluate the actual steady state concurrence of this system
using Eq. (39) for a moderate driving strength of $\epsilon=0.5$. We see that
the dependence of $\mathcal{C}(\tau)$ captures the overall trend inferred from
$\mu_{\rm eff}(\tau)$. However, the exact simulations not only predict
consistently higher values for $\mathcal{C}(\tau=0)$, they also show that the
entanglement of the qubits is considerably more robust with respect to time
delays than the entanglement of the filtered modes. In Fig. 5 (d) we define
the delay time $\tau_{\rm ent}$ as the maximal delay time for which a finite
amount of steady-state entanglement can still be distributed. This timescale
is roughly given by $\tau_{\rm ent}\sim\gamma^{-1}$, but can be significantly
reduced for very large driving strengths.

##  VI Entanglement distribution rates

In the previous sections we have focused on the amount of entanglement that
can be reached under stationary driving conditions. However, for practical
applications it is equally important to know how fast this entangled state can
be reached. In particular, given the possibility of distilling a highly
entangled state from many copies of a state with a low amount of entanglement
Bennett1996 ; Duer1999 ; Macchiavello1999 ; Massar1999 ; Cleland2022 , it can
be more favourable to optimize the generation rate rather than the fidelity.
These considerations are specifically relevant for the current entanglement
distribution scheme since close to the parametric instability, where the
correlations are maximized, the relaxation time of the parametric amplifier
diverges. Therefore, even for an ideal, broadband amplifier, operating close
to the threshold might not be the optimal choice Govia2021 .

For a given application, the optimal compromise between the entanglement
generation speed and the achievable fidelities depends on many details, in
particular also on the local resources that are available to carry out
entanglement purification protocols. To avoid such an application-specific
discussion, we consider here only the following rudimentary scenario: The two-
mode squeezing source is running continuously, while the qubits in each node
are initialized in state $|0\rangle$. At time $t=0$ the coupling between the
qubits and the photonic channels is switched on for a duration $T$, after
which the qubits are decoupled again and stored in a local register. This
process is then repeated with a fresh pair of qubits and so on, such that an
entangled two-qubit state $\rho_{q}(T)$ is distributed between the two nodes
every time interval $T$.

To analyze this scenario, we consider equal qubit decay rates
$\gamma_{i}=\gamma$ and introduce the normalized entanglement distribution
rate

$\mathcal{R}=\frac{E_{F}(T)}{\gamma T}.$ (41)

Here $E_{F}(T)\equiv E_{F}(\rho_{q}(T))$ is the entanglement of formation,
which is related to the concurrence by Horodeki2009

$E_{F}=h\left(\frac{1+\sqrt{1-\mathcal{C}^{2}}}{2}\right),$ (42)

where $h(p)=-p\log_{2}(p)-(1-p)\log_{2}(1-p)$ is the Shannon entropy function.
The entanglement of formation has the following meaning
Wootters97_EntanglementFormation . It quantifies the number of pure singlet
states that are needed on average to generate the state $\rho(T)$ through LOCC
operations only. This more intuitive interpretation, while still being easily
computable, makes $E_{F}$ well suited for a comparative study of entanglement
rates. Note, however, that $E_{F}$ is only an upper bound Horodeki2009  on the
number of singlet states that can be extracted from multiple copies of
$\rho_{q}(T)$, which depends on the available purification protocols and many
other details that go beyond the scope of this analysis.

Figure 6: (a) Plot of the fidelity $\mathcal{F}$ as function of the pulse
length $T$ and (b) the corresponding entanglement distribution rate
$\mathcal{R}$ for different $\epsilon$. For both plots, $\beta=10$,
$\Gamma_{\phi}=0$ and $\eta=1$. (c) Dependence of the maximally achievable
rate, $\mathcal{R}_{\rm max}$, on the driving strength $\epsilon$ for
different values of $\beta$. (d) Plot of the entanglement rate
$\mathcal{R}_{\rm max}$, optimized with respect to both $T$ and $\epsilon$,
${\rm max}_{T,\epsilon}\\{\mathcal{R}\\}$, for different amplifier bandwidths
and waveguide losses.

The behaviour of the entanglement distribution rate as a function of different
parameters is summarized in Fig. 6. First of all, in Fig. 6 (a) we show the
evolution of the fidelity $\mathcal{F}$ as a function of the interaction time
$T$. While for small driving strength we find a featureless increase to the
steady state value discussed in the previous sections, we observe additional
small modulations for large $\epsilon$. These can be traced back to the
appearance of Rabi oscillations between the states $|00\rangle$ and
$|11\rangle$ in this regime. However, these oscillations are not very
pronounced and the fidelity stays well below unity in this transient regime.
Therefore, while the corresponding entanglement rates shown in Fig. 6 (b)
reach their optimum for smaller $T$ when the driving strength is large, higher
overall values for $\mathcal{R}$ can be reached for moderate $\epsilon\lesssim
0.5$. This is shown in more detail in Fig. 6 (c), where the maximal rate
$\mathcal{R}_{\rm max}={\rm max}_{T}\\{\mathcal{R}(T)\\}$ is plotted as
function of $\epsilon$ for different amplifier bandwidths $\beta$. We see that
the rate is maximized at values of $\epsilon\approx 0.3-0.4$.

Finally, in Fig. 6 (d) we plot to maximal entanglement rate as a function of
$\beta$ and the waveguide losses $\eta$, which are the most relevant network
parameters, considering otherwise ideal qubits. Again we find appreciable
rates of $\mathcal{R}\gtrsim 0.1$ for a large range of parameters, including
regimes with significant losses and a comparably low amplifier bandwidth.
While under these conditions Bell-state fidelities close to one are not
achievable, the scheme still performs very well for applications, where the
distribution of a large number of weakly entangled qubits is sufficient.

##  VII Discussion and Conclusions

In summary, we have presented a detailed analysis of a robust entanglement
distribution scheme between distant qubits, which only requires moderately
entangled Gaussian photonic beams. Specifically, we have investigated the
performance of this scheme under realistic experimental conditions, taking a
finite amplifier bandwidth, waveguide losses and also propagation delays into
account. This analysis shows that while reaching extremely high fidelities of
the entangled state still requires close-to-ideal conditions, the scheme is
very efficient in distributing non-ideal entangled states, where it remains
rather robust with respect to common experimental imperfections. For example,
compared to deterministic state-transfer schemes, a finite amount of
entanglement can be distributed with this scheme over very lossy channels. At
the same time, the entanglement distribution rates remain high,
$\mathcal{R}\sim\gamma^{-1}$, and don’t suffer from the weak-driving
requirement and low success probabilities of probabilistic entanglement
distribution schemes.

The analysis presented in this work was explicitly carried out for parametric
amplifiers as one of the most common sources for entangled photon pairs, but
it can be readily generalized for other methods Huard2021  for producing two-
mode squeezed beams. Therefore, all the results can be applied to various
quantum technology platforms, where such sources are currently developed. In
particular, in the field of circuit QED, highly entangled and very pure
sources of entangled beams of microwave photons have been demonstrated in
recent years Eichler2011 ; Menzel2012 ; Huard2012 ; Devoret2016 ; Winkel2020 ;
Peugeot2021 , and used for quantum illumination Fink2020  and quantum
teleportation Fedorov2021  applications. With demonstrated purities as large
as $\mu_{\rm eff}\approx 0.9$ Huard2021 , such sources can be directly applied
to entangle superconducting qubits over microwave channels. Similarly, in the
optical domain two-mode entangled states with more than 10 dB squeezing and
above-threshold purities $\mu_{\rm eff}\geq 1/\sqrt{3}$ Schnabel2007 ;
Schnabel2013 ; Polzik2021  have already been demonstrated and can be further
developed into efficient sources for entangling distant atoms, quantum dots or
defect centers.

An important area for applications for the presented entanglement distribution
scheme are hybrid quantum networks connecting qubits in the microwave and the
optical domain. To entangle qubits with such vastly different frequencies,
optomechanical, electro-optical and other types of quantum transducers Han2021
are currently developed. Most of these transducers can be readily operated as
a parametric amplifier by simply changing the driving frequency and be used to
produce highly entangled beams of microwave and optical photons. The
generation and transduction of entanglement can then be combined and optimized
within a single device, which can significantly boost efficiencies and reduce
the control complexity of such hybrid network architectures.

###### Acknowledgements.

We thank Themis Mavrogordatos and Daoquan Zhu for initial contribution on the
presented topic and Kirill Fedorov for stimulating discussions on entangled
microwave beams. This work was supported by the Austrian Science Fund (FWF)
through Grant No. P32299 (PHONED) and the European Union’s Horizon 2020
research and innovation programme under grant agreement No. 899354
(SuperQuLAN). Most of the computational results presented were obtained using
the CLIP cluster WebCLIP .

##  Appendix A Derivation of the master equation for the qubit state

We start by considering the cascaded master equation for qubits and the
parametric amplifier, which is spelt out in Eq. (6). We first change to an
interaction picture with respect to the free evolution,
$\mathcal{L}_{0}=\mathcal{L}_{p}+\mathcal{L}_{q}$, by writing
$\rho(t)=e^{\mathcal{L}_{0}t}\rho_{I}(t)$. In this representation, the master
equation reads

$\begin{split}\dot{\rho}_{I}(t)=\,&\mathcal{L}_{\rm casc}(t)\rho_{I}(t),\\\
=&\sum_{i=1,2}\sum_{s=\pm}\sqrt{\eta\gamma_{i}\kappa_{i}}
e^{-\mathcal{L}_{0}t}\mathcal{K}_{i}^{s}e^{\mathcal{L}_{0}t}\rho_{I}(t),\end{split}$
(43)

where we defined the cascaded coupling in the interaction picture,
$\mathcal{L}_{\rm casc}(t)=e^{-\mathcal{L}_{0}t}\mathcal{L}_{\rm
casc}e^{\mathcal{L}_{0}t}$, and the superoperators
$\mathcal{K}_{i}^{+}\rho=[a_{i}\rho,\sigma_{i}^{+}]$ and
$\mathcal{K}_{i}^{-}\rho=[\sigma_{i}^{-},\rho a^{\dagger}_{i}]$.

By assuming that the relaxation time of the amplifier is fast compared to the
timescale of the qubit dynamics, we can approximate the full density operator
by the product $\rho_{I}(t)\simeq\rho_{p}^{0}\otimes\rho_{q,I}(t)$, where
$\rho_{p}^{0}$ is the steady state of the photonic system,
$\mathcal{L}_{p}\rho_{p}^{0}=0$. Up to second order in $\mathcal{L}_{\rm
casc}(t)$ and under the validity of the usual Born-Markov approximation, we
can then derive a master equation for the reduced qubit state
$\rho_{q,I}(t)=\mathrm{Tr}_{p}\\{\rho_{I}(t)\\}$, which is given by

$\dot{\rho}_{q,I}(t)=\int_{-\infty}^{t}\mathrm{d}\tau\,\mathrm{Tr}_{p}\\{\mathcal{L}_{\rm
casc}(t)\mathcal{L}_{\rm casc}(\tau)\rho_{p}^{0}\otimes \rho_{q,I}(t)\\}.$
(44)

After inserting the explicit expression for $\mathcal{L}_{\rm casc}(t)$ and
undoing the transformation to the interaction picture, the result can be
written in a compact form as

$\begin{split}\dot{\rho}_{q}(t)=&\mathcal{L}_{q}\rho_{q}(t)+\sum_{i,j=1}^{2}\eta\sqrt{\gamma_{i}\gamma_{j}}\sum_{s,s^{\prime}=\pm}\int_{-\infty}^{t}\mathrm{d}\tau\,C_{ij}^{ss^{\prime}}(t-\tau)\\\
&\times
ss^{\prime}[\sigma^{-s}_{i},e^{\mathcal{L}_{q}(t-\tau)}[\sigma_{j}^{-s^{\prime}},e^{-\mathcal{L}_{q}(t-\tau)}\rho_{q}(t)]].\end{split}$
(45)

Here we have introduced the bosonic correlation functions

$C_{ij}^{ss^{\prime}}(\tau)=\sqrt{\kappa_{i}\kappa_{j}}\langle:a_{i}^{s}(\tau)a_{j}^{s^{\prime}}(0):\rangle,$
(46)

where we identified $a_{i}^{+}\equiv a_{i}^{\dagger}$ and $a_{i}^{-}\equiv
a_{i}$ and assumed the normal ordering prescription
$\langle:a^{s}(\tau)a(0):\rangle=\mathrm{Tr}_{p}\\{a^{s}e^{\mathcal{L}_{p}\tau}(a\rho_{p}^{0})\\}=\langle
a^{s}(\tau)a(0)\rangle$, while
$\langle:a^{s}(\tau)a^{\dagger}(0):\rangle=\mathrm{Tr}_{p}\\{a^{s}e^{\mathcal{L}_{p}\tau}(\rho_{p}^{0}a^{\dagger})\\}=\langle
a^{\dagger}(0)a^{s}(\tau)\rangle$.

Consistent with the Markov approximation, we also neglect the slow dynamics of
the qubits in Eq. (45), which leaves us with

$\begin{split}\dot{\rho}_{q}(t)=\mathcal{L}_{q}\rho_{q}(t)+&\sum_{i,j=1}^{2}\eta\sqrt{\gamma_{i}\gamma_{j}}\\\
\times&\sum_{s,s^{\prime}=\pm}ss^{\prime}I_{ij}^{ss^{\prime}}(0)[\sigma^{-s}_{i},[\sigma_{j}^{-s^{\prime}},\rho_{q}(t)]],\end{split}$
(47)

where
$I_{ij}^{ss^{\prime}}(\omega)=\int_{0}^{\infty}\mathrm{d}\tau\,C_{ij}^{ss^{\prime}}(\tau)e^{-i\omega\tau}$
is evaluated at resonance $\omega=0$. To proceed, we make use of the general
relation $I_{ij}^{s,s}=(I_{ij}^{-s,-s})^{*}$ and additional simplifications
$I_{ii}^{s,s}=I_{12}^{s,-s}=I_{21}^{s,-s}=0$, which we derive in Appendix B
for the setup at hand. With these relations and identifying
$N_{i}=2\eta\mathrm{Re}\\{I_{ii}^{+-}(0)\\}$ and
$M=\eta(I_{12}^{--}(0)+I_{21}^{--}(0))$, we end up with the master equation in
Eq. (8) together with a small Hamiltonian correction

$H^{\prime}=-\eta\sum_{i}\gamma_{i}{\rm Im}\\{I_{ii}^{+-}(0)\\}
\sigma_{i}^{z}.$ (48)

These frequency shifts can be compensated by local qubit detunings and we also
find that ${\rm Im}\\{I_{ii}^{+-}(0)\\}=0$ for the symmetric parameters
considered in most of our examples. Therefore, this corrections is neglected
in the main part of the paper.

##  Appendix B Correlation Functions and Spectra of the Parametric Amplifier

In this appendix, we summarize the general results for the spectra of the two
amplifier modes for arbitrary parameters. Since in the operation regime of
interest the amplifier is linear, this can be done most conveniently by
converting the master equation in Eq. (4) into an equivalent set of quantum
Langevin equations Textbooks_ZollerNoise ; Collett1984  for the Heisenberg
operators $a_{1}(t)$ and $a_{2}(t)$ (see also Appendix D). These equations
read

$\displaystyle\dot{a}_{1}=-\left(i\Delta_{1}+\frac{\kappa_{1}}{2}\right)a_{1}+\frac{\sqrt{\kappa_{1}\kappa_{2}}\epsilon}{2}a_{2}^{\dagger}-\sqrt{\kappa_{1}}
f_{{\rm in},1},$ (49)
$\displaystyle\dot{a}_{2}=-\left(i\Delta_{2}+\frac{\kappa_{2}}{2}\right)a_{2}+\frac{\sqrt{\kappa_{1}\kappa_{2}}\epsilon}{2}a_{1}^{\dagger}-\sqrt{\kappa_{2}}
f_{{\rm in},2},$ (50)

where $f_{{\rm in},i}$ are independent white noise operators satisfying
$[f_{{\rm in},i}(t),f^{\dagger}_{{\rm
in},j}(t^{\prime})]=\delta_{ij}\delta(t-t^{\prime})$. By introducing the
vectors $\boldsymbol{v}=(a_{1},a_{2}^{\dagger},a_{2},a_{1}^{\dagger})^{\top}$
and $\boldsymbol{f}=(\sqrt{\kappa_{1}}f_{{\rm in},1},\sqrt{\kappa_{2}}f_{{\rm
in},2}^{\dagger},\sqrt{\kappa_{2}}f_{{\rm in},2},\sqrt{\kappa_{1}}f_{{\rm
in},1}^{\dagger})^{\top}$, these equations can be written in a compact form as

$\dot{\boldsymbol{v}}=\mathcal{M}\boldsymbol{v}-\boldsymbol{f},$ (51)

where the matrix $\mathcal{M}$ is given by

$\mathcal{M}=\begin{pmatrix}-i\Delta_{1}-\frac{\kappa_{1}}{2}&\frac{\epsilon\sqrt{\kappa_{1}\kappa_{2}}}{2}&0&0\\\
\frac{\epsilon\sqrt{\kappa_{1}\kappa_{2}}}{2}&i\Delta_{2}-\frac{\kappa_{2}}{2}&0&0\\\
0&0&-i\Delta_{2}-\frac{\kappa_{2}}{2}&\frac{\epsilon\sqrt{\kappa_{1}\kappa_{2}}}{2}\\\
0&0&\frac{\epsilon\sqrt{\kappa_{1}\kappa_{2}}}{2}&i\Delta_{1}-\frac{\kappa_{1}}{2}\\\
\end{pmatrix}.$ (52)

For long times, $t\rightarrow\infty$, the formal solution of Eq. (51) is

$\boldsymbol{v}(t)=-\int_{-\infty}^{t}d\tau\,e^{\mathcal{M}(t-\tau)}\boldsymbol{f}(\tau).$
(53)

From this result we obtain the full covariance matrix in steady state,
$\mathcal{V}_{0}=\langle\boldsymbol{v}\boldsymbol{v}^{\dagger}\rangle(t\rightarrow\infty)$,
as

$\mathcal{V}_{0}=\int_{0}^{\infty}ds\,e^{\mathcal{M}s}\mathcal{R}e^{\mathcal{M}^{\dagger}s},$
(54)

where $\mathcal{R}={\rm diag}(\kappa_{1},0,\kappa_{2},0)$ is a diagonal
matrix. Since $\mathcal{M}$ is block-diagonal, the matrix exponential and
therefore the individual entries of the covariance matrix can be solved
analytically. We obtain

$\displaystyle\langle
a^{\dagger}_{i}a_{i}\rangle=\frac{(\bar{\kappa}-\kappa_{i})\bar{\kappa}\epsilon^{2}}{4\bar{\Delta}^{2}+(1-\epsilon^{2})\bar{\kappa}^{2}},$
(55) $\displaystyle\langle
a_{1}a_{2}\rangle=\frac{\sqrt{\kappa_{1}\kappa_{2}}(-2i\bar{\Delta}+\bar{\kappa})\epsilon}{4\bar{\Delta}^{2}+(1-\epsilon^{2})\bar{\kappa}^{2}},$
(56)

where $\bar{\kappa}=\kappa_{1}+\kappa_{2}$ and
$\bar{\Delta}=\Delta_{1}+\Delta_{2}$. All other expectation values vanish,
i.e., $\langle a_{1}^{\dagger}a_{2}\rangle=\langle a_{i}^{2}\rangle=0$.

To evaluate the spectra we make use of the quantum regression theorem

$\partial_{\tau}\langle\boldsymbol{v}(\tau)\boldsymbol{v}^{\dagger}\rangle=\mathcal{M}\langle\boldsymbol{v}(\tau)\boldsymbol{v}^{\dagger}\rangle,$
(57)

and define

$\begin{split}\mathcal{I}(\omega)=&\int_{0}^{\infty}d\tau\,\langle\boldsymbol{v}(\tau)\boldsymbol{v}^{\dagger}\rangle
e^{-i\omega\tau}=\frac{1}{i\omega\mathbbm{1}-\mathcal{M}}\mathcal{V}_{0}.\end{split}$
(58)

The spectra for the photon occupation numbers and correlations can be obtained
elementwise as
$I_{a_{1}^{\dagger}a_{1}}(\omega)=\kappa_{1}[\mathcal{I}(\omega)]_{44}$,
$I_{a_{2}^{\dagger}a_{2}}(\omega)=\kappa_{2}[\mathcal{I}(\omega)]_{22}$,
$I_{a_{1}a_{2}}(\omega)=\sqrt{\kappa_{1}\kappa_{2}}[\mathcal{I}(\omega)]_{12}$
and
$I_{a_{2}a_{1}}(\omega)=\sqrt{\kappa_{1}\kappa_{2}}[\mathcal{I}(\omega)]_{34}$.
The general expressions for the spectra are rather lengthy, but for the
symmetric case, $\kappa_{i}=\kappa$ and $\Delta_{i}=\Delta$, we obtain

$I_{a^{\dagger}_{i}a_{i}}(\omega)=\frac{2\epsilon^{2}\kappa^{3}(\kappa+i\omega)}{(\kappa^{2}-\delta^{2})(\kappa^{2}-\delta^{2}+4i\kappa\omega-4\omega^{2})},$
(59)

and

$I_{a_{1}a_{2}}(\omega)=\frac{\epsilon\kappa^{2}(\delta^{2}+\kappa^{2}-4i\Delta(\kappa+i\omega)+2i\omega\kappa)}{(\kappa^{2}-\delta^{2})(\kappa^{2}-\delta^{2}+4i\kappa\omega-4\omega^{2})},$
(60)

where $\delta^{2}=\epsilon^{2}\kappa^{2}-4\Delta^{2}$. These results reduce to
the expressions in Eq. (12) and Eq. (13) for $\Delta=0$. Note that for a
symmetric amplifier ${\rm Im}\\{I_{a^{\dagger}_{i}a_{i}}(0)\\}=0$ and there is
no Lamb-shift contribution. In general this is not the case.

##  Appendix C Hilbertspace truncation

In our numerical simulations of the full cascaded master equation we truncate
the Hilbertspace of each photon mode and only include a finite set of number
states $|n\rangle$ with $n<n_{\rm trunc}$. In Fig. 7 (a) and (b) we plot the
entanglement fidelity for different truncation numbers $n_{\rm trunc}$ and two
very different values of $\beta$. Fig. 7 (c) shows the corresponding average
photon number in each of the amplifier modes. We see that even though for
driving strengths of $\epsilon\lesssim 0.8$ the average number of photons
remains small, $\langle a^{\dagger}_{i}a_{i}\rangle<1$, in the case of large
$\beta$ a substantially larger $n_{\rm trunc}$ must be used in order to
capture the qubit-qubit entanglement accurately. For even larger $\epsilon$
the average photon number grows considerably and exact simulations of the full
master equation become very demanding. For this reasons all our simulations
are restricted to the parameter regime $\epsilon\leq 0.8$ and $\beta\leq
10^{3}$, where a convergence to the required level of accuracy can be achieved
with $n_{\rm trunc}=25$ [see Fig. 7 (d)].

Figure 7: Dependence of the numerically evaluated fidelity $\mathcal{F}$ on
the photon truncation number $n_{\rm trunc}$ for (a) $\beta=1$ and (b)
$\beta=10^{3}$. (c) Plot of the average number of photons in each of the
parametric amplifier modes as a function of the driving strength and for
different $n_{\rm trunc}$. (d) Convergence of the fidelity $\mathcal{F}$ with
increasing photon truncation number for the case of $\beta=10^{3}$ at
$\epsilon=0.8$. For all plots, $\Delta_{i}=0$, $\eta=1$, $\Gamma_{\phi}=0$ and
otherwise symmetric conditions have been assumed.

##  Appendix D Quantum Langevin equations for cascaded networks with
propagation delays

Instead of working with a master equation for the system density operator
$\rho(t)$, an equivalent description of open quantum systems can be obtained
from a set of quantum Langevin equations for Heisenberg operators
Textbooks_ZollerNoise . For the considered setup, these quantum Langevin
equations can be derived from a unidirectional system-waveguide interaction of
the form

$H_{\rm
int}=i\sum_{i=1,2}\left[\sqrt{\kappa_{i}}F^{\dagger}_{i}(0)a_{i}+\sqrt{\gamma_{i}}F^{\dagger}_{i}(d_{i})\sigma_{i}^{-}-{\rm
H.c.}\right],$ (61)

where

$F_{i}(z)=\frac{1}{\sqrt{2\pi}}\int_{0}^{\infty}d\omega\,e^{i\omega
z/c_{i}}b_{i}(\omega)$ (62)

is the field operator for the $i$-th channel at position $z$ and
$[b_{i}(\omega),b^{\dagger}_{j}(\omega^{\prime})]=\delta_{ij}\delta(\omega-\omega^{\prime})$.
The equation of motion for an arbitrary system operator $O(t)$ then obeys

$\begin{split}\dot{O}=&\,i[H_{\rm
sys},O]-\sum_{i=1,2}\left(\sqrt{\kappa_{i}}[O,a_{i}^{\dagger}]F_{i}(0,t)+{\rm
H.c.}\right)\\\ &-\sum_{i=1,2}\left(\sqrt{\gamma_{i}}[O,\sigma_{i}^{+}]
F_{i}(d_{i},t)+{\rm H.c.}\right),\end{split}$ (63)

where $H_{\rm sys}$ is the bare system Hamiltonian. In turn, under the
validity of the usual Born-Markov approximation, the field operator is given
by Textbooks_ZollerNoise

$\begin{split}F_{i}(z,t)=&\,F_{{\rm
in},i}(z,t)+\sqrt{\kappa_{i}}\Theta(z)a_{i}(t-z/c_{i})\\\
&+\sqrt{\gamma_{i}}\Theta(z-d_{i})\sigma_{i}^{-}(t-z/c_{i}+\tau_{i}),\end{split}$
(64)

where we have omitted nonessential propagation phases and $\Theta(x)$ denotes
the unit step function. The field $F_{{\rm in},i}(z,t)$ is the free field
operator in the $i$-th channel and obeys

$F_{{\rm in},i}(z,t+\tau)=F_{{\rm in},i}(z-c_{i}\tau,t).$ (65)

Since we neglect thermal excitations, we find that $F_{{\rm
in},i}(z,t)\rho_{\rm full}^{0}=\rho_{\rm full}^{0}F^{\dagger}_{{\rm
in},i}(z,t)=0$, where $\rho_{\rm full}^{0}$ is the initial state of the full
network. We identify $f_{{\rm in},i}(t)=F_{{\rm in},i}(z=0,t)$ with the noise
operators used in Appendix B.

Note that for every position $z$, the fields $F_{i}(z,t)$ represent
independent degrees of freedom and thus

$\displaystyle[F_{i}(z,t),O(t)]$ $\displaystyle=$ $\displaystyle 0,$ (66)
$\displaystyle\left[F_{i}(z,t),F_{j}(z^{\prime},t)\right]$ $\displaystyle=$
$\displaystyle c_{i}\delta_{ij}\delta(z-z^{\prime}).$ (67)

Together with Eq. (64) and Eq. (65), these relations can be use to derive all
the non-equal time correlation functions used below.

###  D.1 Weak driving limit

From Eq. (63) we can readily derive expectation values of qubit observables in
the weak excitation limit, $\epsilon\ll 1$. In this limit we can approximate
$\sigma^{z}_{i}(t)\simeq-1$ and obtain

$\begin{split}\dot{\sigma}_{i}^{-}(t)\simeq-\frac{\gamma_{i}}{2}
\sigma_{i}^{-}(t)-\sqrt{\gamma_{i}\eta}f_{{\rm out},i}(t-\tau_{i})\\\
-\sqrt{\gamma_{i}(1-\eta)} f_{{\rm in},i}^{\prime}(t).\end{split}$ (68)

Here $f_{{\rm out},i}(t)\equiv F_{{\rm out},i}(z=0,t)$, where $F_{{\rm
out},i}(z,t)=F_{{\rm in},i}(z,t)+\sqrt{\kappa_{i}}a_{i}(t-z/c_{i})$, and
$f_{{\rm in},i}^{\prime}(t)$ is an independent noise operator, which we have
included to account for waveguide losses. In the limit $t\rightarrow\infty$ we
obtain

$\sigma_{i}^{-}(t)=-\sqrt{\gamma_{i}\eta}\int_{-\infty}^{t}{\rm
d}s\,e^{-\gamma_{i}(t-s)/2}f_{{\rm out},i}(s-\tau_{i}),$ (69)

where contributions from $f_{{\rm in},i}^{\prime}(t)$, which always act on the
vacuum state, have already been omitted. For the evaluation of the expectation
value $\langle\sigma_{i}^{+}\sigma_{i}^{-}\rangle(t\rightarrow\infty)=N_{i}$
we use

$\langle f^{\dagger}_{{\rm out},i}(t)f_{{\rm
out},i}(t^{\prime})\rangle=\kappa_{i}\langle
a^{\dagger}_{i}(t)a_{i}(t^{\prime})\rangle$ (70)

and after some manipulations we obtain the result for $N_{i}$ given in Eq.
(29). For the evaluation of the correlations,
$\langle\sigma_{1}^{-}\sigma_{2}^{-}\rangle(t\rightarrow\infty)=M$, we must
take into account that $f_{{\rm in},1}(t)$ and $a_{2}(t^{\prime})$ do not
commute in general and Collett1985

$\langle f_{{\rm out},1}(t)f_{{\rm
out},2}(t^{\prime})\rangle=\sqrt{\kappa_{1}\kappa_{2}}\langle\mathcal{T}a_{1}(t)a_{2}(t^{\prime})\rangle,$
(71)

where $\mathcal{T}$ denotes the time-ordering operator. Again we find that the
resulting expression for $M$ matches Eq. (30). Therefore, this comparison
shows that the results obtained within the FMA for the steady state of the
qubits become exact in the weak driving limit, $N_{i},|M|\ll 1$.

###  D.2 Propagation delays

To handle finite propagation delays in a more general manner, we introduce a
set of shifted Heisenberg operators

$\tilde{\sigma}_{i}^{k}(t)=\sigma_{i}^{k}(t+\tau_{i}),\qquad\tilde{a}_{i}(t)=a_{i}(t+\tau_{\varepsilon}).$
(72)

Here, $\tau_{\varepsilon}$ is a small time delay, which is negligible on the
timescale of the system dynamics, but must be kept finite when evaluating
commutation relations with the waveguide fields.

The shifted Heisenberg operators obey (assuming $\eta=1$ for simplicity)

$\displaystyle\dot{\tilde{a}}_{i}(t)$ $\displaystyle=$ $\displaystyle
i[\tilde{H}_{p}(t),\tilde{a}_{i}(t)]-\frac{\kappa_{i}}{2}\tilde{a}_{i}(t)$
(73) $\displaystyle-\sqrt{\kappa_{i}}f_{{\rm in},i}(t+\tau_{\varepsilon}),$
$\displaystyle\dot{\tilde{\sigma}}_{i}^{-}(t)$ $\displaystyle=$
$\displaystyle-\frac{\gamma_{i}}{2}\tilde{\sigma}_{i}^{-}(t)+\sqrt{\gamma_{i}\kappa_{i}}\tilde{\sigma}_{i}^{z}(t)a_{i}(t)$
(74) $\displaystyle+\sqrt{\gamma_{i}}\tilde{\sigma}_{i}^{z}(t)f_{{\rm
in},i}(t),$ $\displaystyle\dot{\tilde{\sigma}}_{i}^{z}(t)$ $\displaystyle=$
$\displaystyle-2\gamma_{i}\tilde{\sigma}_{i}^{+}(t)\tilde{\sigma}_{i}^{-}(t)$
(75)
$\displaystyle-2\sqrt{\gamma_{i}\kappa_{i}}\left[\tilde{\sigma}_{i}^{+}(t)a_{i}(t)+a^{\dagger}_{i}(t)\tilde{\sigma}_{i}^{-}(t)\right]$
$\displaystyle-2\sqrt{\gamma_{i}}\left[\tilde{\sigma}_{i}^{+}(t)f_{{\rm
in},i}(t)+f^{\dagger}_{{\rm in},i}(t)\tilde{\sigma}_{i}^{-}(t)\right].$

Since in these equations all $f_{{\rm in},i}(t)$ operators appear to the right
and all $f^{\dagger}_{{\rm in},i}(t)$ operators to the left, we can perform
the expectation value with respect to the initial vacuum state $\rho_{\rm
full}^{0}$ and take the limit $\tau_{\epsilon}\rightarrow 0$ afterwards. As a
result, the expectation values for
$\langle\dot{\tilde{\sigma}}_{i}^{k}(t)\rangle$ and
$\langle\dot{\tilde{a}}_{i}(t)\rangle$ do not explicitly depend on the delay
times $\tau_{i}$ anymore and their expressions are identical to the ones
obtained for $\langle\dot{\sigma}_{i}^{k}(t)\rangle$ and
$\langle\dot{a}_{i}(t)\rangle$ from the cascaded master equation given in Eq.
(6).

As a next step we show that the same is true for arbitrary operator products
$\tilde{S}_{1}(t)\tilde{S}_{2}(t)\tilde{A}(t)$, where $\tilde{S}_{i}(t)$ are
Pauli operators and $\tilde{A}(t)$ is an arbitrary product of operators
$\tilde{a}_{1,2}(t)$ and $\tilde{a}^{\dagger}_{1,2}(t)$. To evaluate the time
derivative of this product, we apply the product rule and use the time
derivatives for the individual operators given in Eqs. (73)-(75). This results
in terms of the form

$\begin{split}&\frac{d}{dt}\tilde{S}_{1}(t)\tilde{S}_{2}(t)\tilde{A}(t)=\dots\\\
&-\sqrt{\gamma_{1}}[\tilde{S}_{1}(t),\tilde{\sigma}_{1}^{+}(t)]F_{1}(d_{1},t+\tau_{1})\tilde{S}_{2}(t)\tilde{A}(t)+\dots\\\
&-\sqrt{\gamma_{2}}\tilde{S}_{1}(t)F_{2}^{\dagger}(d_{2},t+\tau_{2})[\tilde{\sigma}_{2}^{-}(t),\tilde{S}_{2}(t)]\tilde{A}(t)+\dots\\\
&-\sqrt{\kappa_{1}}\tilde{S}_{1}(t)\tilde{S}_{2}(t)F_{1}^{\dagger}(0,t+\tau_{\epsilon})[\tilde{a}_{1}(t),\tilde{A}(t)]+\dots\end{split}$
(76)

Before we can take the expectation value with respect to $\rho_{\rm
full}^{0}$, all $f_{{\rm in},i}(t)$ operators must be commuted to the right
and all $f^{\dagger}_{{\rm in},i}(t)$ operators to the left. To do so we
write, for example,

$F_{1}(d_{1},t+\tau_{1})=\frac{\sqrt{\gamma_{1}}}{2}\tilde{\sigma}_{1}^{-}(t)+F_{{\rm
out},1}(0,t),$ (77)

and use $F_{{\rm out},1}(0,t)=F_{1}(c\tau_{\epsilon},t+\tau_{\epsilon})$ to
show that $[F_{{\rm out},1}(0,t),\tilde{A}(t)]=0$. Further, we write
$\tilde{S}_{2}(t)=S_{2}(t)+\Delta S_{2}(t)$, where $\Delta S_{2}(t)$ depends
on the field $F_{{\rm out},2}(d_{2},t^{\prime})$ for times
$t^{\prime}\in[t,t+\tau_{2}]$ only. Equivalently, it depends on the field
$F_{2}(z,t)$ located in the region $z\in(0,d_{2}]$ and the operator

$F_{{\rm out},2}(0,t)=F_{2}(0,t)+\frac{\sqrt{\kappa_{2}}}{2}a_{2}(t).$ (78)

This implies that also $[\tilde{S}_{2}(t),F_{{\rm out},1}(0,t)]=0$ and
$F_{{\rm out},1}(0,t)=f_{{\rm in},1}(t)+\sqrt{\kappa_{1}}a_{1}(t)$ in the
second line of Eq. (76) can be commuted all the way to the right. Similar
arguments can also be made for all the other terms to achieve the desired
operator ordering.

In summary, from this derivation we obtain a set of coupled equations of
motion for the expectation values of arbitrary operator products
$\langle\tilde{S}_{1}(t)\tilde{S}_{2}(t)\tilde{A}(t)\rangle$, which are
independent of the delay times $\tau_{i}$ and have the same structure as the
corresponding equations of motion derived from the time-local cascaded master
equation. After taking again the limit $\tau_{\epsilon}\rightarrow 0$, this
result implies that

$\langle S_{1}(t+\tau_{1})S_{2}(t+\tau_{2})A(t)\rangle=\langle
S_{1}S_{2}A\rangle(t)\Big{|}_{\rm loc},$ (79)

assuming appropriately matched initial conditions.

In the main text we are interested in equal-time expectation values
$\lim_{t\rightarrow\infty}\langle S_{1}(t)S_{2}(t)\rangle$, for which Eq. (79)
can not be directly applied. Instead we repeat the whole derivation for the
operators $\tilde{S}_{1}(t+t_{0})$ and take the average with respect to the
state $S_{2}(t_{0}+\tau_{2})\rho_{\rm full}(t_{0})$. Since
$S_{2}(t_{0}+\tau_{2})$ depends on $f_{{\rm in},i}(t)$ for $t\leq t_{0}$ only,
it commutes with the relevant noise terms $\sim f_{{\rm in},i}(t+t_{0})$ and
we can still make use of $f_{{\rm
in},i}(t+t_{0})S_{2}(t_{0}+\tau_{2})\rho_{\rm full}(t_{0})=0$. Therefore, this
approach provides us with the relation

$\displaystyle\langle S_{1}(t+\tau_{1}+t_{0})S_{2}(t_{0}+\tau_{2})\rangle$
$\displaystyle=$ (80) $\displaystyle\langle
S_{1}(t+t_{0})S_{2}(t_{0})\rangle\Big{|}_{\rm loc},$

which extends the result from above to more general correlations. By assuming
that $t_{0}$ is long enough such that the system has reached a steady state,
we can redefine $t_{0}\rightarrow t_{0}-\tau_{2}$ and set
$t=\tau_{2}-\tau_{1}$. This leaves us with the result stated in Eq. (38).

## References

  * (1) M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information, Cambridge University Press, Cambridge (2000).
  * (2) L. Jiang, J. M. Taylor, A. S. Sørensen, and M. D. Lukin, Distributed quantum computation based on small quantum registers, Phys. Rev. A 76, 062323 (2007).
  * (3) H. J. Kimble, The quantum internet, Nature (London) 453, 1023 (2008).
  * (4) C. Monroe, R. Raussendorf, A. Ruthven, K. R. Brown, P. Maunz, L.-M. Duan, and J. Kim, Large-scale modular quantum-computer architecture with atomic memory and photonic interconnects, Phys. Rev. A 89, 022317 (2014).
  * (5) S. Wehner, D. Elkouss, and R. Hanson, Quantum internet: A vision for the road ahead, Science 362, 303 (2018).
  * (6) A. S. Cacciapuoti, M. Caleffi, F. Tafuri, F. S. Cataliotti, S. Gheradini, and G. Bianchi, Quantum Internet: Networking Challenges in Distributed Quantum Computing, IEEE Network 34, 137 (2020).
  * (7) C. H. Bennett, G. Brassard, S. Popescu, B. Schumacher, J. A. Smolin, and W. K. Wootters, Purification of noisy entanglement and faithful teleportation via noisy channels, Phys. Rev. Lett. 76, 722 (1996).
  * (8) W. Dür, H.-J. Briegel, J. Cirac, and P. Zoller, Quantum repeaters based on entanglement purification, Phys. Rev. A 59, 169 (1999).
  * (9) J. I. Cirac, A. K. Ekert, S. F. Huelga, and C. Macchiavello, Distributed quantum computation over noisy channels, Phys. Rev. A. 59, 4249 (1999).
  * (10) A. Kent, N. Linden, and S. Massar, Optimal Entanglement Enhancement for Mixed States, Phys. Rev. Lett. 83, 2656 (1999).
  * (11) H. Yan, Y. Zhong, H. -S. Chang, A. Bienfait, M.-H. Chou, C. R. Conner, E. Dumur, J. Grebel, R. G. Povey, and A. N. Cleland, Entanglement purification and protection in a superconducting quantum network, Phys. Rev. Lett. 128, 080504 (2022).
  * (12) D. Gottesman and I. L. Chuang, Demonstrating the viability of universal quantum computation using teleportation and single-qubit operations, Nature 402, 390 (1999).
  * (13) T. E. Northup and R. Blatt, Quantum information transfer using photons, Nat. Photonics 8, 356 (2014).
  * (14) J. I. Cirac, P. Zoller, H. J. Kimble, and H. Mabuchi, Quantum State Transfer and Entanglement Distribution among Distant Nodes in a Quantum Network, Phys. Rev. Lett. 78, 3221 (1997).
  * (15) P. Kurpiers, P. Magnard, T. Walter, B. Royer, M. Pechal, J. Heinsoo, Y. Salathe, A. Akin, S. Storz, J.-C. Besse, S. Gasparinetti, A. Blais, and A. Wallraff, Deterministic quantum state transfer and remote entanglement using microwave photons, Nature (London) 558, 264 (2018).
  * (16) K. S. Chou, J. Z. Blumoff, C. S. Wang, P. C. Reinhold, C. J. Axline, Y. Y. Gao, L. Frunzio, M. H. Devoret, L. Jiang, and R. J. Schoelkopf, Deterministic teleportation of a quantum gate between two logical qubits, Nature (London) 561, 368 (2018).
  * (17) P. Magnard, S. Storz, P. Kurpiers, J. Schär, F. Marxer, J. Lütolf, T. Walter, J.-C. Besse, M. Gabureac, K. Reuer, A. Akin, B. Royer, A. Blais, and A. Wallraff, Microwave Quantum Link between Superconducting Circuits Housed in Spatially Separated Cryogenic Systems, Phys. Rev. Lett. 125, 260502 (2020).
  * (18) Y. Zhong, H.-S. Chang, A. Bienfait, E. Dumur, M. -H. Chou, C. R. Conner, J. Grebel, R. G. Povey, H. Yan, D. I. Schuster, and A. N. Cleland, Deterministic multi-qubit entanglement in a quantum network, Nature (London) 590, 571 (2021).
  * (19) L.-M. Duan, M. D. Lukin, J. I. Cirac, and P. Zoller, Long-distance quantum communication with atomic ensembles and linear optics, Nature (London) 414, 413 (2001).
  * (20) S. D. Barrett and P. Kok, Efficient high-fidelity quantum computation using matter qubits and linear optics, Phys. Rev. A 71, 060310(R) (2005).
  * (21) D. N. Matsukevich, T. Chaneliere, S. D. Jenkins, S.-Y. Lan, T. A. B. Kennedy, and A. Kuzmich, Entanglement of remote atomic qubits, Phys. Rev. Lett. 96, 030405 (2006).
  * (22) D. L. Moehring, P. Maunz, S. Olmschenk, K. C. Younge, D. N. Matsukevich, L.-M. Duan, and C. Monroe, Entanglement of single-atom quantum bits at a distance, Nature (London) 449, 68 (2007).
  * (23) B. Hensen, H. Bernien, A. E. Dreau, A. Reiserer, N. Kalb, M. S. Blok, J. Ruitenbern, R. F. L. Vermeulen, R. N. Schouten, C. Abellan, W. Amaya, V. Pruneri, M. W. Mitchell, M. Markham, D. J. Twitchen, D. Elkouss, S. Wehner, T. H. Taniniau, and R. Hanson, Loophole-free Bell inequality violation using electronic spins separated by 1.3 kilometers, Nature (London) 526, 682 (2015).
  * (24) R. Riedinger, A. Wallucks, I. Marinkovic, C. Löschnauer, M. Aspelmeyer, S. Hong, and S. Gröblacher, Remote quantum entanglement between two micromechanical oscillators, Nature (London) 556, 473 (2018).
  * (25) H. J. Carmichael, Statistical Methods in Quantum Optics 2: Non-Classical Fields, Springer, Berlin (1999).
  * (26) P. G. Kwiat, K. Mattle, H. Weinfurter, and A. Zeilinger, New High-Intensity Source of Polarization-Entangled Photon Pairs, Phys. Rev. Lett. 75, 4337 (1995).
  * (27) T. Kim, M. Fiorentino, and F. N. C. Wong, Phase-Stable source of polarization-entangled photons using a polarization Segnac interferometer, Phys. Rev. A 73, 012316 (2006).
  * (28) C. Eichler, D. Bozyigit, C. Lang, M. Baur, L. Steffen, J. M. Fink, S. Filipp, and A. Wallraff, Observation of two-mode squeezing in the microwave frequency domain, Phys. Rev. Lett. 107, 113601 (2011).
  * (29) E. P. Menzel, R. Di Candia, F. Deppe, P. Eder, L. Zhong, M. Ihmig, M. Haeberlein, A. Baust, E. Hoffmann, D. Ballester, K. Inomata, T. Yamamoto, Y. Nakamura, E. Solano, A. Marx, and R. Gross, Path Entanglement of Continuous-Variable Quantum Microwaves, Phys. Rev. Lett. 109, 250502 (2012).
  * (30) E. Flurin, N. Roch, F. Mallet, M. H. Devoret, and B. Huard, Generating Entangled Microwave Radiation Over Two Transmission Lines, Phys. Rev. Lett. 109, 183901 (2012).
  * (31) A. Roy and M. Devoret, Introduction to parametric amplification of quantum signals with Josephson circuits, C. R. Phys. 17, 740 (2016).
  * (32) P. Winkel, I. Takmakov, D. Rieger, L. Planat, W. Hasch-Guichard, L. Grünhaupt, N. Maleeva, F. Foroughi, F. Henriques, K. Borisov, J. Ferrero, A. V. Ustinov, W. Wernsdorfer, N. Roch, and I. M. Pop, Nondegenerate Parametric Amplifiers Based on Dispersion-Engineered Josephson-Junction Arrays, Phys. Rev. Applied 13, 024015 (2020).
  * (33) A. Peugeot, G. Menard, S. Dambach, M. Westig, B. Kubala, Y. Mukharsky, C. Altimiras, P. Joyez, D. Vion, P. Roche, D. Esteve, P. Milman, J. Leppäkangas, G. Johansson, M. Hofheinz, J. Ankerhold, and F. Portier, Generating Two Continuous Entangled Microwave Beams Using a dc-Biased Josephson Junction, Phys. Rev. X. 11, 031008 (2021).
  * (34) J. P. Torres, K. Banaszek, and I. A. Walmsley, Engineering Nonlinear Optic Sources of Photonic Entanglement, Progress in Optics 56, 227 (2011).
  * (35) C. Couteau, Spontaneous parametric down-conversion, Cont. Phys. 59, 291 (2018).
  * (36) B. Kraus and J. I. Cirac, Discrete entanglement distribution with squeezed light, Phys. Rev. Lett 92, 013602 (2004).
  * (37) M. Paternostro, W. Son, and M. S. Kim, Complete Conditions for Entanglement Transfer, Phys. Rev. Lett. 92, 197901 (2004).
  * (38) A. Serafini, M. Paternostro, M. S. Kim, and S. Bose, Enhanced dynamical entanglement transfer with multiple qubits, Phys. Rev. A 73, 022312 (2006).
  * (39) C. A. Muschik, E. S. Polzik, and J. I. Cirac, Dissipatively driven entanglement of two macroscopic atomic ensembles, Phys. Rev. A. 83, 052312 (2011).
  * (40) S. Felicetti, M. Sanz, L. Lamata, G. Romero, G. Johansson, P. Delsing, and E. Solano, Dynamical Casimir Effect Entangles Artificial Atoms, Phys. Rev. Lett. 113, 093602 (2014).
  * (41) A. L. Grimsmo and A. Blais, Squeezing and quantum state engineering with Josephson travelling wave amplifiers, npj Quantum Information 3, 20 (2017).
  * (42) N. Didier, J. Guillaud, S. Shankar, and M. Mirrahimi, Remote entanglement stabilization and concentration by quantum reservoir engineering, Phys. Rev. A 98, 012329 (2018).
  * (43) J. You, Z. Liao, S.-W. Li, and M. S. Zubairy, Waveguide quantum electrodynamics in squeezed vacuum, Phys. Rev. A 97, 023810 (2018).
  * (44) L. C. G. Govia, A. Lingenfelter, and A. A. Clerk, Stabilizing two-qubit entanglement by mimicking a squeezed environment, Phys. Rev. Research 4, 023010 (2022).
  * (45) Z.-L. Xiang, S. Ashhab, J. Q. You, and F. Nori, Hybrid quantum circuits: Superconducting circuits interacting with other quantum systems, Rev. Mod. Phys. 85, 623 (2013).
  * (46) G. Kurizki, P. Bertet, Y. Kubo, K. Mølmer, D. Petrosyan, P. Rabl, and J. Schmiedmayer, Quantum technologies with hybrid systems, PNAS 112, 3866 (2015).
  * (47) K. Stannigel, P. Rabl, and P. Zoller, Driven-dissipative preparation of entangled states in cascaded quantum-optical networks, New J. Phys. 14, 063014 (2014).
  * (48) P. Lodahl, S. Mahmoodian, S. Stobbe, P. Schneeweiss, J. Volz, A. Rauschenbeutel, H. Pichler, and P. Zoller, Chiral quantum optics, Nature (London) 541, 473 (2017).
  * (49) A. Rueda, W. Hease, S. Barzanjeh, and J. M. Fink, Electro-optic entanglement source for microwave to telecom quantum state transfer, npj Quantum Information 5, 108 (2019).
  * (50) S. Krastanov, H. Raniwala, J. Holzgrafe, K. Jacobs, M. Loncar, M. J. Reagor, and D. R. Englund, Optically Heralded Entanglement of Superconducting Systems in Quantum Networks, Phys. Rev. Lett. 127, 040503 (2021).
  * (51) E. Flurin, N. Roch, J. D. Pillet, F. Mallet, and B.Huard, Superconducting Quantum Node for Entanglement and Storage of Microwave Radiation, Phys. Rev. Lett. 114, 090503 (2015).
  * (52) H. J. Carmichael, Quantum trajectory theory for cascaded open systems, Phys. Rev. Lett. 70, 2273 (1993).
  * (53) C. W. Gardiner, Driving a Quantum System with the Output Field From Another Driven Quantum System, Phys. Rev. Lett. 70, 2269 (1993).
  * (54) A. Serafini, F. Illuminati, M. G. A. Paris, and S. De Siena, Entanglement and purity of two-mode Gaussian states in noisy channels, Phys. Rev. A. 69, 022318 (2004).
  * (55) R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. 81, 865 (2009).
  * (56) K. G. Fedorov, S. Pogorzalek, U. Las Heras, M. Sanz, P. Yard, P. Eder, M.Fischer, J.Goetz, E. Xie, K. Inomata, Y. Nakamura, R. Di Candia, E. Solano, A. Marx, F. Deppe, and R. Gross, Finite-time quantum entanglement in propagating squeezed microwaves, Sci. Rep. 8, 6416 (2018).
  * (57) S. Hill and W. K. Wootters, Entanglement of a Pair of Quantum Bits, Phys. Rev. Lett. 78, 5022 (1997).
  * (58) R. Dassonneville, R. Arrouly, T. Peronnin, A. A. Clerk, A. Bienfait, and B. Huard, Dissipative Stabilization of Squeezing Beyond 3 dB in a Microwave Mode, PRX Quantum 2, 020323 (2021).
  * (59) J. DiGuglielmo, B. Hage, A. Franzen, J. Fiurášek, and R. Schnabel, Experimental characterization of Gaussian quantum-communication channels, Phys. Rev. A 76, 012323 (2007).
  * (60) T. Eberle, V. Händchen, and R. Schnabel, Stable control of 10 dB two-mode squeezed vacuum states of light, Opt. Express 21, 11546 (2013).
  * (61) T. B. Brasil, V. Novikov, H. Kerdoncuff, M. Lassen, and E. Polzik, Two-color high-purity Einstein-Podolsky-Rosen photonic state, arXiv: 2110.00066 (2021).
  * (62) S. Barzankeh, S. Pirandola, D. Vitali, and J. M. Fink, Microwave quantum illumination using a digital receiver, Sci. Adv. 6, eabb0451 (2020).
  * (63) K. G. Fedorov, M. Renger, S. Pogorzalek, R. Di Candia, Q. Chen, Y. Nojiri, K. Inomata, Y. Nakamura, M. Partanen, A. Marx, R. Gross, and F. Deppe, Experimental quantum teleportation of propagating microwaves, Sci. Adv. 7, abk0891 (2021).
  * (64) X. Han, W. Fu, C.-L. Zou, L. Jiang, and H. X. Tang, Microwave-optical quantum frequency conversion, Optica 8, 1050 (2021).
  * (65) https://clip.science/
  * (66) C. W. Gardiner and P. Zoller, Quantum noise, Springer, Berlin (2004).
  * (67) M. J. Collett and C. W. Gardiner, Squeezing of intracavity and traveling-wave fields produced in parametric amplification, Phys. Rev. A 30, 1386 (1984).
  * (68) C. W. Gardiner and M. J. Collett, Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation, Phys. Rev. A 31, 3761 (1985).