# High-Purity Entanglement of Hot Propagating Modes Using Nonreciprocity

L. Orr  Dahlem Center for Complex Quantum Systems and Fachbereich Physik,
Freie Universität Berlin, 14195 Berlin, Germany  Saeed A. Khan  Department of
Electrical Engineering, Princeton University, Princeton, NJ 08544, USA  N.
Buchholz  Dahlem Center for Complex Quantum Systems and Fachbereich Physik,
Freie Universität Berlin, 14195 Berlin, Germany  S. Kotler  Racah Institute of
Physics, The Hebrew University of Jerusalem, Jerusalem, 91904, Israel  A.
Metelmann  Dahlem Center for Complex Quantum Systems and Fachbereich Physik,
Freie Universität Berlin, 14195 Berlin, Germany  Institute for Theory of
Condensed Matter, Karlsruhe Institute of Technology, 76131 Karlsruhe, Germany
Institute for Quantum Materials and Technology, Karlsruhe Institute of
Technology, 76344 Eggenstein-Leopoldshafen, Germany

###### Abstract

Distributed quantum information processing and communication protocols demand
the ability to generate entanglement among propagating modes. However, thermal
fluctuations can severely limit the fidelity and purity of propagating
entangled states, especially for low-frequency modes relevant for radio-
frequency (RF) signals. Here we propose nonreciprocity as a resource to render
continuous-variable entanglement of propagating modes robust against thermal
fluctuations. By utilising a cold-engineered reservoir we break the symmetry
of reciprocity in a standard two-mode squeezing interaction between a low- and
a high-frequency mode, and show that the rerouting of thermal fluctuations
allows the generation of flying entangled states with high purity. Our
approach requires only pairwise Gaussian interactions and is thus ideal for
parametric circuit QED implementations.

###### pacs:

84.30.Le 03.65.Ta, 42.50.Pq 42.50.Lc

##  I Introduction

Entanglement of propagating photons [1] is a crucial resource for quantum
information processing and communication protocols [2] and is useful for
distributing entanglement amongst components of a quantum network [3, 4].
However, as with other coherent quantum effects, it is remarkably sensitive to
decoherence channels, such as thermal fluctuations. Operating at cryogenic
temperatures allows for the effects of thermal fluctuations to be overcome by
ensuring that $k_{b}T\ll\hbar\omega$, which is possible for mode frequencies
$\omega$ as low as the microwave domain [5]. For even lower frequency bands,
such as the radio-frequency (RF) domain that is ubiquitous in modern
communication, thermal fluctuations remain appreciable even at the lowest
operating temperatures [6, 7, 8], presenting challenges for RF quantum
communication and sensing. With these cryogenic temperature limitations,
alternative approaches to generate entanglement in systems of hot modes must
be considered.

In this paper, we present an approach utilising engineered nonreciprocity to
generate steady-state entangled output fields from a system of interacting hot
modes coupled to a single cold mode. While several proposals and recent
experiments consider the entanglement of the output fields of cold modes
coupled via an intermediate hot (for example, mechanical) mode [9, 10, 11, 12,
13, 14, 15], we consider situations where the fields to be entangled
themselves are effectively coupled to high temperature baths. Continuous
ambient thermal excitations can severely limit the entanglement fidelity of
steady-state emission from such “hot” modes at a given pump power.
Furthermore, these excitations limit the purity of the generated flying
states, demanding the use of complex state purification protocols [16, 17,
18]. We show that nonreciprocity provides a crucial ingredient to alleviate
these effects: the ability to continuously reroute thermal excitations toward
a cold output. This enables the entanglement of propagating photons with
increased robustness to thermal excitations, and with much higher purity than
is possible using a completely reciprocal two-mode entangling interaction
between the hot modes of interest.

The importance of nonreciprocity has already been firmly established in
quantum information processing, enabling the routing of signals in a quantum
network by realising asymmetric scattering matrices, across diverse
architectures from superconducting circuits [19, 20, 21, 22, 51, 24, 25, 26,
27] to optomechanics [28, 29, 30, 31, 32] and beyond [33, 34]. Our work
analyses an aspect of nonreciprocal interactions which is much less explored:
the role of nonreciprocity in manipulating fluctuations in a quantum system,
to route thermal noise while generating entanglement of the scattered fields.
Building on recent progress in the theory of engineered nonreciprocity [19,
35, 36, 37], we consider a system of three dissipative quantum modes
undergoing configurable coherent interactions, and identify the conditions
required for nonreciprocal scattering and directional transmission.
Interestingly, by analysing the complete output state we find that
entanglement can be enhanced at points of “perfect nonreciprocity”, where
scattering in one direction is forbidden. Perhaps equally as importantly, it
is also possible to engineer nonreciprocal scattering between a pair of modes
without entangling their outputs, highlighting the need for a deeper
understanding of the connections between nonreciprocity and entanglement. To
this end, we develop a heuristic picture drawing connections between steady-
state entanglement in nonreciprocal systems and sequential Gaussian circuit
operations, as well as dissipative entanglement schemes [38, 12], and ideal
two-mode squeezing.

With these foundations, our work finally addresses the impact of
nonreciprocity on entanglement in quantum systems experiencing thermal noise.
The control of thermal noise flow [39] has recently garnered renewed interest
even in classical devices, due to its importance in energy harvesting, heat
management, and information transfer using thermal currents [40, 41]. For
quantum systems, which are our specific focus, the control of thermal noise
flow becomes particularly important [42, 43, 44] in order to protect fragile
quantum properties from thermal decoherence. To this end we show how
nonreciprocal scattering can be engineered to continuously reroute incident
thermal excitations away from hot modes, toward the output of the cold
auxiliary mode introduced to break reciprocity. We then show that
nonreciprocity can increase the entanglement fidelity and state purity of
output fields scattered off the hot modes, above values that are possible
using a reciprocal two-mode squeezing interaction at the same strength. Our
heuristic picture shows that this increased robustness is due to a controlled
swap of the input noise incident on the quantum modes at different
temperatures.

The proposed three-mode system can be efficiently realised in parametric
circuit QED (cQED), where time-dependent pump fields enable tunable
interactions to break reciprocity [45]. Furthermore, we demand only pairwise
squeezing and beam-splitter interactions, capabilities for which have been
suitably demonstrated in recent cQED experiments [20, 22], and which can be
achieved using even a single nonlinear element based on a Josephson junction.
Our model can therefore serve as a practical platform for the detailed study
of quantum entanglement in the presence of nonreciprocal interactions and
thermal noise.

The rest of this paper is organised as follows. In Sec. II, starting with a
two-mode squeezing interaction, we introduce the minimal three-mode system
required to render this interaction nonreciprocal, within the context of
standard approaches to nonreciprocity. In Secs. III and IV, we proceed to
analyse the scattering and entangling properties of the three-mode system,
finding conditions for nonreciprocal scattering, and clarifying the connection
between nonreciprocity and entanglement generation. We find that at the
specific points of “perfect nonreciprocity,” the scattering and entangling
properties of the system can be very efficiently explained as a sequence of
simple pairwise linear operations. Finally, in Sec. V we combine this
understanding to explore the impact of thermal fluctuations on entanglement in
nonreciprocal systems and demonstrate how thermal inputs to a hot mode can be
efficiently routed via nonreciprocity to protect the entanglement of scattered
output fields.

##  II Setup

We begin with the standard description of a nondegenerate two-mode squeezing
(TMS) interaction between two harmonic modes (setting $\hbar=1$),

$\displaystyle\hat{H}_{\rm
TMS}=\,\,g_{12}\left(\hat{a}^{\dagger}_{1}\hat{a}^{\dagger}_{2}+\hat{a}_{1}\hat{a}_{2}\right),$
(1)

where $\hat{a}_{j}$ is the bosonic annihilation operator for mode $j$,
satisfying the standard commutation relations
$[\hat{a}_{j},\hat{a}^{\dagger}_{k}]=\delta_{jk}$. In the cQED architecture
such an interaction is typically realised by appropriately pumping nonlinear
Josephson-junction based superconducting elements [46, 47]. Additionally, this
interaction can be used to generate entangled photon pairs [48] and hence two-
mode squeezed light for quantum information processing applications [49, 50].
Nevertheless, the interaction defined by Eq. (1) is reciprocal.

The conditions required to render interactions of the form of Eq. (1)
nonreciprocal have been clarified in recent years [36, 37]. An arbitrary
bidirectional interaction between two systems, $\hat{H}_{\rm
int}\propto(\hat{A}\hat{B}+h.c.)$ governed by operators, $\hat{A}$ and
$\hat{B}$, must be balanced with a corresponding nonlocal dissipative
interaction $\Gamma\mathcal{D}[\hat{z}]$ as depicted in Fig. 1, where
$\mathcal{D}[\hat{z}]\hat{\rho}=\hat{z}\hat{\rho}\hat{z}^{\dagger}-\frac{1}{2}\\{\hat{z}^{\dagger}\hat{z},\hat{\rho}\\}$
is the standard dissipative superoperator, with collapse operator
$\hat{z}=\hat{A}+\eta e^{i\phi}\hat{B}$. An appropriately chosen interaction
strength $\Gamma$, asymmetry $\eta$, and, most crucially, phase $\phi$ [37],
can then be used to render the desired interaction nonreciprocal.

Figure 1: (a) A diagram of an open two-mode squeezer where both modes are
coupled to the common nonlocal dissipator $\Gamma\mathcal{D}[\hat{z}]$. (b) A
diagram of a minimal loop system consisting of three open modes. Modes $a_{1}$
and $a_{2}$ are coupled via a two-mode squeezing interaction, as are modes
$a_{2}$ and $a_{3}$. Modes $a_{1}$ and $a_{3}$ are coupled via a beam-splitter
interaction. The tunable loop phase is placed on the interaction between modes
$a_{2}$ and $a_{3}$. If mode $a_{3}$ is adiabatically eliminated this loop is
equivalent to the system depicted in (a) with the jump operator given by
$\hat{z}=\hat{a}_{1}+\eta e^{i\phi}\hat{a}^{\dagger}_{2}$.

Applying this approach to the interaction defined by Eq. (1), it is clear that
dissipators with either the collapse operator $\hat{z}=\hat{a}_{1}+\eta
e^{i\phi}\hat{a}_{2}^{\dagger}$ or $\hat{z}=\hat{a}_{1}^{\dagger}+\eta
e^{-i\phi}\hat{a}_{2}$ both satisfy the aforementioned form, and thus can be
employed to attain the desired nonreciprocal scattering matrix. Importantly,
it has been shown that either dissipative interaction alone can generate
steady-state entanglement [38, 12]. Hence, by combining either dissipator with
the coherent two-mode squeezing interaction in Eq. (1), one can also take the
point of view that we are analysing the effects of introducing nonreciprocity
in such dissipative entanglement schemes.

A hint as to what may be expected can be found in the fact that the required
operators $\hat{z}$ are non-Hermitian nonlocal collapse operators.
Nonreciprocal interactions mediated by dissipators with Hermitian collapse
operators have been shown to be equivalent to measurement-based feedforward
schemes [37]: a classical observer makes a measurement on system A, and uses
the result to evolve system B, breaking the reciprocity of interaction between
the systems. As such, this evolution is equivalent to performing local
operations and classical communication, and hence cannot generate any
entanglement between the two quantum systems. In contrast, the non-Hermitian
collapse operators required here have no such mapping, and can in principle
generate entanglement.

To realise either nonlocal dissipator and hence render the TMS interaction
nonreciprocal, we must introduce an auxiliary mode $a_{3}$, as depicted in
Fig. 1(b). The open quantum system comprising this loop and the environment
with which it interacts are described by the quantum optical master equation,

$\displaystyle\dot{\hat{\rho}}=\mathcal{L}\hat{\rho}=-i[\hat{H}_{\rm
NRL},\hat{\rho}]+\sum_{j=1,2,3}\kappa_{j}\mathcal{D}[\hat{a}_{j}]\hat{\rho},$
(2)

and its dynamics are governed by the Hamiltonian

$\displaystyle\hat{H}_{\rm NRL}$
$\displaystyle=\left(g_{12}\hat{a}^{\dagger}_{1}\hat{a}^{\dagger}_{2}+g_{13}\hat{a}^{\dagger}_{3}\left[\hat{a}_{1}+\eta
e^{i\phi}\hat{a}^{\dagger}_{2}\right]\right)+h.c.,$ (3)

where $\eta=g_{23}/g_{13}$ accounts for an asymmetric coupling to the
auxiliary mode, written in the interaction frame with respect to the three
modes. Eq. (3) simply describes the original two-mode squeezing interaction
between modes $a_{1}$ and $a_{2}$, but now with an auxiliary third mode that
couples to mode $a_{1}$ via a beam-splitter interaction, and to mode $a_{2}$
via a two-mode squeezing interaction. When the auxiliary mode $a_{3}$ can be
adiabatically eliminated (namely, when its damping rate $\kappa_{3}$ is the
largest system parameter), this configuration realises the nonlocal dissipator
$\hat{z}=\hat{a}_{1}+\eta e^{i\phi}\hat{a}_{2}^{\dagger}$. Engineering of the
alternative nonlocal dissipator introduced earlier leads to equivalent
results. The resulting three-mode system can thus render the interaction of
Eq. (1) unidirectional and can realise a nonreciprocal loop (NRL) via the
general scheme described above. By explicitly including the dynamics of the
auxiliary mode, we are able to explore the routing of both scattered fields
and - importantly - their correlations around the loop.

We now take a moment to discuss quantum optics platforms that can be used to
realise our proposed system, and its practical implications. Eq. (3) requires
only two-mode squeezing and beam-splitter interactions between pairs of modes;
systems where such interactions can be controllably realised in a three-mode
configuration to enable nonreciprocity have already been demonstrated in
parametric cQED [20, 22]. For concreteness, in Appendix A, we detail how a
simple three-mode circuit with a single nonlinear mixing element can realise
the NRL. The required pairwise interactions have also been realised across
other platforms such as optomechanical circuits [23, 52, 53]. Importantly, the
three modes constituting the NRL can have widely distinct frequencies, and
thus experience thermal effects to varying degrees (see Appendix A). The
coupling of modes with such disparate frequencies has been realised in very
recent experiments [54, 55], and is highly relevant to RF-domain quantum
optics. In precisely these cases, the NRL provides a way of enhancing the
entanglement and purity of hot propagating modes using auxiliary cold modes,
which is a key result of our work.

##  III Scattering properties and nonreciprocity

When analysing the three-mode loop, we restrict the initial state
$\hat{\rho}_{\mathrm{in}}$ of the system to Gaussian states. Each mode is
equipped with a pair of quadrature operators,
$\hat{X}_{j}=(\hat{a}^{\dagger}_{j}+\hat{a}_{j})/\sqrt{2}$ and
$\hat{P}_{j}=i(\hat{a}^{\dagger}_{j}-\hat{a}_{j})/\sqrt{2}$ where
$j\in\\{1,2,3\\}$. We can also define quadrature modes for the input or output
states by replacing $j$ with $(j,\\{\mathrm{in},\mathrm{out}\\})$,
respectively. By declaring the initial state $\hat{\rho}_{\mathrm{in}}$ to be
Gaussian, we mean that it is completely characterised by the first and second
moments of these quadrature operators. Since the dynamics induced by the
Hamiltonian in Eq. (3) are entirely linear, the state of the system at any
time will remain Gaussian, including the output state
$\hat{\rho}_{\mathrm{out}}$ [56].

Our approach to analysing the steady-state scattering properties of the system
is standard: we solve the linear Heisenberg-Langevin equations in frequency
space, and use quantum input-output theory [57] (see Appendix B). From this,
we obtain the scattering matrix $\mathbf{S}[\omega]$ relating output-field
operators to the input fields,

$\vec{R}_{\mathrm{out}}[\omega]=\mathbf{S}[\omega]\vec{R}_{\mathrm{in}}[\omega],$
(4)

where $\vec{R}_{\\{\mathrm{in},\mathrm{out}\\}}[\omega]$ is a vector of
quadrature operators for the input and output modes, respectively.

We are interested in the scattering behaviour on resonance, that is, when
$\omega=0$ in this frame, which describes the response of fields resonant with
the individual modes comprising the system. The scattering matrix in this
simpler case, $\mathbf{S}[0]\equiv\mathbf{S}$, can be expressed purely in
terms of the cooperativities
$\mathcal{C}_{jk}=4g_{jk}^{2}/\kappa_{j}\kappa_{k}$ parameterising
interactions between modes $a_{j}$ and $a_{k}$ in relation to their individual
decay rates. By balancing the cooperativities and adjusting the loop phase
$\phi$, scattering between any pair of modes in the system can be rendered
nonreciprocal. The full form of the scattering matrix is still unwieldy (see
Appendix C). We therefore introduce measures that allow us to quantify the
scattering properties of the system more compactly.

For nonreciprocal systems, we are primarily interested in the asymmetry of
scattering between modes $j$ and $k$. To more precisely quantify this
asymmetry, we introduce the normalised degree of nonreciprocity
$\mathcal{N}^{(j,k)}$ as

$\displaystyle\mathcal{N}^{(j,k)}=\frac{||{\rm abs}~{}\mathbf{S}_{jk}-{\rm
abs}~{}\mathbf{S}_{kj}||}{||\mathbf{S}_{jk}||+||\mathbf{S}_{kj}||}$ (5)

where $\mathbf{S}_{jk}$ is the two-mode block of the scattering matrix
corresponding to modes $j$ and $k$, ${\rm abs}~{}\mathbf{O}$ is an element-
wise absolute value operation, and $||\cdot||$ is the Frobenius norm.
$\mathcal{N}^{(j,k)}$ is a quantity that remains bounded within $[0,1]$, and
measures the difference in amplitude (and not phase) of scattering between a
pair of modes. As a result, it vanishes for scattering that is reciprocal in
amplitude but that may differ in phase. As discussed in Appendix C, it can be
shown that $\mathcal{N}^{(j,k)}=0~{}\forall~{}j\neq k$ only when $\phi=0$.

For all other values of $\phi$, the three-mode system exhibits nonreciprocal
scattering properties to varying degrees. To understand how this
nonreciprocity influences the quantum properties of scattered input fields, we
can analyse the covariance matrix of the output fields. However, before
analysing the general case, we find that several key ideas can be understood
via a simple heuristic picture that is valid when $\mathcal{N}^{(j,k)}=1$.
From Eq. (5), this corresponds to perfectly asymmetric scattering, where
either $||\mathbf{S}_{jk}||$ or $||\mathbf{S}_{kj}||$ vanishes. We refer to
these as points of perfect nonreciprocity and the required conditions are
summarised below. The arrows denote the direction in which signal transmission
is allowed; scattering matrix elements in the reverse direction vanish
exactly:

$\begin{array}[]{|c||c|c|}\hline\cr&\phi=-\pi/2&\phi=+\pi/2\\\
\hline\cr\hline\cr\mathbf{\mathcal{C}_{12}=\mathcal{C}_{13}\mathcal{C}_{23}}&\bm{a_{1}\rightarrow
a_{2}}&\bm{a_{1}\leftarrow a_{2}}\\\
\hline\cr\mathcal{C}_{23}=\mathcal{C}_{12}\mathcal{C}_{13}&a_{2}\rightarrow
a_{3}&a_{2}\leftarrow a_{3}\\\
\hline\cr\mathcal{C}_{13}=\mathcal{C}_{12}\mathcal{C}_{23}&a_{1}\rightarrow
a_{3}&a_{1}\leftarrow a_{3}\\\ \hline\cr\end{array}.$ (6)

In this work, we pay particular attention to scattering between modes $a_{1}$
and $a_{2}$, and hence set
$\mathcal{C}_{12}=\mathcal{C}_{13}\mathcal{C}_{23}$. Plotting
$\mathcal{N}^{(1,2)}$ in Fig. 2(a), we see that scattering between modes
$a_{1}$ and $a_{2}$ can be rendered perfectly nonreciprocal,
$\mathcal{N}^{(1,2)}=1$, when $\phi=\pm\pi/2$.

###  III.1 Circuit decomposition at points of perfect nonreciprocal scattering

The scattering matrix describes a potentially complicated set of linear
operations on Gaussian states. In order to more easily understand the
scattering behaviour of the system, a variety of decomposition schemes can be
used to represent this set of operations more efficiently. One such prominent
example is the Bloch-Messiah decomposition [58], a special case of the
singular value decomposition for the group of real symplectic matrices, of
which the steady-state scattering matrix $\mathbf{S}$ is an element,
$\mathbf{S}\in Sp(2n,\mathbb{R})$.

We find that for the perfectly nonreciprocal system under consideration, the
less commonly used polar decomposition [59] proves to be simpler. This
decomposition allows us to write a symplectic matrix in the form
$\mathbf{RU}$, referred to as its left polar decomposition, where
$\mathbf{U}\in Sp(2n,\mathbb{R})\cap O(2n,\mathbb{R})$ is a real symplectic
orthogonal matrix, and $\mathbf{R}\in
Sp(2n,\mathbb{R})\cap\mathrm{Sym}^{+}(2n)$ is a real symplectic symmetric
positive definite matrix. Physically, the matrix $\mathbf{U}$ represents
passive optical transformations, namely beam-splitters and phase shifters. The
matrix $\mathbf{R}$ then includes any single and two-mode squeezing
interactions. An equivalent form $\mathbf{UR}$ of the scattering matrix is
provided by the right polar decomposition, where $\mathbf{R}$ and $\mathbf{U}$
are in general distinct from the left polar decomposition.

In general, the polar decomposition leads to dense and complicated matrices
$\mathbf{R}$ and $\mathbf{U}$. Remarkably, we find that the scattering matrix
describing our three-mode system has extremely simple forms for the left
and/or right polar decomposition at points of perfect nonreciprocal
scattering, $\mathcal{N}^{(1,2)}=1$, up to a global change in phase on the
scattering matrix, $-\mathbf{S}$. Due to the structure of the covariance
matrix [see Eq. (10)] this change of phase will not affect the resulting
covariances. In these simple cases, the $\mathbf{R}$ and $\mathbf{U}$ matrices
involve only a single interaction between one pair of modes. The polar
decomposition when $\mathcal{N}^{(1,2)}=1$ is then comprised of symplectic
matrices corresponding to the following two unitary operations:

$\displaystyle\mathbf{U}^{(1,3)}\leftrightarrow$
$\displaystyle\exp\left[2i\arctan\Big{(}\\!\sqrt{\mathcal{C}_{13}}\Big{)}\left(\hat{a}^{\dagger}_{1}\hat{a}_{3}+\hat{a}_{1}\hat{a}^{\dagger}_{3}\right)\right]$
$\displaystyle\mathbf{R}^{(2,3)}\leftrightarrow$
$\displaystyle\exp\left[-2\,\mathrm{artanh}\Big{(}\\!\sqrt{\mathcal{C}_{23}}\Big{)}\left(\hat{a}^{\dagger}_{2}\hat{a}^{\dagger}_{3}-\hat{a}_{2}\hat{a}_{3}\right)\right].$
(7)

Crucially, we find that both left and right polar decompositions can provide
useful complementary insights into the action of the three-mode system at
points of perfect nonreciprocal scattering.

Figure 2: (a) The scattering properties of the nonreciprocal loop as a
function of the loop phase $\phi$. $\mathcal{N}^{(1,2)}$ is plotted in blue
for the impedance-matched case only; the general case is qualitatively similar
and is omitted for clarity. $||\mathbf{S}_{11}||$ is plotted in orange for the
impedance-matched (solid) and general (dashed) cases. (b) Circuit descriptions
of the nonreciprocal system at points of perfect nonreciprocal scattering,
where $\mathcal{N}^{(1,2)}=1$. The form of these circuits holds regardless of
the value of $||\mathbf{S}_{11}||$. The left diagram corresponds to
$\phi=-\pi/2$ which yields perfect nonreciprocal scattering in the
$a_{1}\rightarrow a_{2}$ direction, while the right diagram corresponds to
$\phi=\pi/2$ where scattering is only allowed in the $a_{1}\leftarrow a_{2}$
direction. The input is initially uncorrelated and moves from bottom to top.

For the loop phase $\phi=-\pi/2$, it is the left polar decomposition that
takes on a simple form:

$-\mathbf{S}=\mathbf{R}^{(2,3)}\,\mathbf{U}^{(1,3)}\text{ when
}\phi=-\frac{\pi}{2}.$ (8)

As illustrated in Fig. 2, modes $a_{1}$ and $a_{3}$ interact first via a beam-
splitter $\mathbf{U}^{(1,3)}$ which acts to exchange input from $a_{1}$ to
$a_{3}$, and vice versa. This operation is followed by a two-mode squeezer
$\mathbf{R}^{(2,3)}$ between modes $a_{2}$ and $a_{3}$, where the output of
mode $a_{2}$ then becomes dependent on the input of mode $a_{1}$, thus
realising directional transmission from $a_{1}\rightarrow a_{2}$. The output
of mode $a_{1}$ cannot have any dependence on the input of $a_{2}$ because it
can only arrive at the output of mode $a_{3}$ via the same two-mode squeezer.

For the opposite sign of the loop phase, $\phi=+\pi/2$, the right polar
decomposition yields

$-\mathbf{S}=\mathbf{U}^{(1,3)}\,\mathbf{R}^{(2,3)}\text{ when
}\phi=\frac{\pi}{2},$ (9)

which describes the same component operations as Eq. (8) but applied in
reverse order. As a result, the scattering behaviour is reversed and still
nonreciprocal, allowing transmission from $a_{1}\leftarrow a_{2}$.

The relatively simple form of these circuits provides a heuristic picture of
nonreciprocal scattering in this system, where changing the direction of the
nonreciprocal scattering is equivalent to changing the order of operations in
the circuit. However, note that the nonreciprocal behaviour is not explained
just by the sequential beam-splitters, but also involves two-mode squeezing
interactions. This already hints at the possibility of generating nontrivial
quantum correlations in scattered output fields and thus connects to
entangling properties of the three-mode system, as we will see in Sec. IV.

##  IV Output entanglement and purity

The output state of the three-mode system, and all of its entangling
capabilities, is completely characterised by the generally frequency-dependent
covariance matrix of the output-field quadrature operators,
$\mathbf{V}[\omega]$, which can be calculated using Eq. (4) and the known
correlation relations of the input-field operators. The output covariance
matrix can then be written in terms of the scattering matrix

$\mathbf{V}[\omega]=\frac{1}{2}\left(\mathbf{S}[\omega]\mathbf{V}_{\mathrm{in}}\mathbf{S}^{T}[-\omega]+\mathbf{S}[-\omega]\mathbf{V}_{\mathrm{in}}\mathbf{S}^{T}[\omega]\right),$
(10)

where $\mathbf{V}_{\mathrm{in}}$ is the matrix of correlations of the input
fields. Assuming that the input noise for different modes is uncorrelated,
$\mathbf{V}_{\mathrm{in}}$ contains variances determined by the thermal
occupation number $n_{j}^{\mathrm{th}}$ and vacuum fluctuations (see Appendix
B),

$\mathbf{V}_{\mathrm{in}}=\bigoplus_{j=1}^{n}\left(n_{j}^{\mathrm{th}}+\frac{1}{2}\right)\mathbf{I}$
(11)

where $n=3$ is the number of modes and $\mathbf{I}$ is the $2\times 2$
identity matrix. Once again considering the response on resonance, we set
$\omega=0$ in Eq. (10) so that the covariance matrix of interest,
$\mathbf{V}[0]\equiv\mathbf{V}$, takes the simple form

$\mathbf{V}=\mathbf{S}\mathbf{V}_{\mathrm{in}}\mathbf{S}^{T}.$ (12)

From this covariance matrix, we aim to calculate useful entanglement metrics
for different bipartitions of the output fields in order to investigate the
effects of nonreciprocal scattering. These include the Simon-Peres-Horodecki
criterion for the separability of Gaussian states [60], as well as the
logarithmic negativity $E_{N}^{(j,k)}$, an entanglement monotone for the
shared output state of modes $a_{j}$ and $a_{k}$ [61, 62, 56] for the shared
output state of modes $a_{j}$ and $a_{k}$. The latter may be calculated from
the minimum symplectic eigenvalue $\nu_{-}^{(j,k)}$ of the partial transpose
of the corresponding two-mode block $\mathbf{V}^{(j,k)}$ from the total
covariance matrix $\mathbf{V}$ via

$E^{(j,k)}_{N}=\begin{cases}0&\quad\text{for }2\nu_{-}^{(j,k)}\geq 1\\\
-\log\left(2\nu_{-}^{(j,k)}\right)&\quad\text{for
}2\nu_{-}^{(j,k)}<1.\end{cases}$ (13)

Similarly, we can define the marginal purity $\mu^{(j,k)}$ of a given
bipartition of the output field of modes using

$\mu^{(j,k)}=\frac{1}{4\sqrt{\det\mathbf{V}^{(j,k)}}}.$ (14)

This measure has a maximum value of $\mu^{(j,k)}=1$ only for pure states; for
mixed states the purity will be $\mu^{(j,k)}<1$.

###  IV.1 Entanglement and purity in a nonreciprocal system

We begin by examining the entanglement properties and purity of output fields
under vacuum input, $n^{\mathrm{th}}_{j}=0$, for all modes. In this instance
the initial matrix of correlations is comprised of vacuum noise and is
therefore the identity matrix, $\mathbf{V}_{\mathrm{in}}=\mathbf{I}_{6}/2$, so
the output covariance matrix is
$\mathbf{V}=\frac{1}{2}\mathbf{S}\mathbf{S}^{T}$. Balancing the
cooperativities $\mathcal{C}_{12}=\mathcal{C}_{13}\mathcal{C}_{23}$, we plot
the logarithmic negativity between the output of modes $a_{1}$ and $a_{2}$,
$E_{N}^{(1,2)}$, as well as $E_{N}^{(2,3)}$, in Fig. 3(a), where we note the
strong dependence on the value of the loop phase. The simple form for the
circuit decomposition means that the behaviour of the entanglement at the
points of perfect nonreciprocal scattering between modes $a_{1}$ and $a_{2}$,
corresponding to the points where $\phi=\pm\pi/2$ in Fig. 3, can be explained
simply as follows.

The covariance matrix for the $a_{1}\rightarrow a_{2}$ direction of perfect
nonreciprocal scattering (where $\phi=-\pi/2$) may be written using the form
of the scattering matrix from Eq. (8):

$\mathbf{V}=\frac{1}{2}\mathbf{R}^{(2,3)}\big{(}\mathbf{R}^{(2,3)}\big{)}^{T}$
(15)

where the beam-splitter component does not appear because it is an orthogonal
transformation. The entangling behaviour for this direction of the
nonreciprocal scattering is therefore equivalent to a two-mode squeezer
between modes $a_{2}$ and $a_{3}$; hence we must have $E^{(2,3)}_{N}>0$. Since
modes $a_{1}$ and $a_{2}$ do not share any squeezing in this representation,
there will be no entanglement generated between these two modes, which is
evident in Fig. 3.

We can then ask whether there are other operating points where the
entanglement between the output of modes $a_{1}$ and $a_{2}$ vanishes. It is
possible to determine this for all system parameters, and not just at the
points of perfect nonreciprocity, by examining the Simon-Peres-Horodecki
criterion for the output of modes $a_{1}$ and $a_{2}$:

$\sqrt{\frac{\mathcal{C}_{12}}{\mathcal{C}_{13}\mathcal{C}_{23}}}+\sqrt{\frac{\mathcal{C}_{13}\mathcal{C}_{23}}{\mathcal{C}_{12}}}\leq-2\sin\phi.$
(16)

The output fields of modes $a_{1}$ and $a_{2}$ are separable so long as the
above inequality is satisfied, which only occurs for one set of parameters:
$\mathcal{C}_{12}=\mathcal{C}_{13}\mathcal{C}_{23}$ and $\phi=-\pi/2$, which
is the point of perfect nonreciprocal scattering where $a_{1}\rightarrow
a_{2}$. This is, then, the only point of operation where the entanglement of
the output of modes $a_{1}$ and $a_{2}$ vanishes: $E^{(1,2)}_{N}=0$. Away from
this point, the output for these two modes will always be entangled.

For the reverse direction, $a_{1}\leftarrow a_{2}$ (where $\phi=+\pi/2$), we
can use Eq. (9) to write the covariance matrix as

$\mathbf{V}=\frac{1}{2}\mathbf{U}^{(1,3)}\mathbf{R}^{(2,3)}\big{(}\mathbf{R}^{(2,3)}\big{)}^{T}\\!\big{(}\mathbf{U}^{(1,3)}\big{)}^{T}.$
(17)

Importantly, the beam-splitter between modes $a_{1}$ and $a_{3}$ appears and
therefore plays an important role in the entanglement generation here. The
two-mode squeezer acts first to entangle modes $a_{2}$ and $a_{3}$, while the
later action of the beam-splitter swaps some of these squeezed correlations
from mode $a_{3}$ to $a_{1}$, generating entanglement between the output of
modes $a_{1}$ and $a_{2}$. This resembles the dissipative entanglement
protocol [38, 12], where the entanglement of two modes is mediated by a
strongly damped auxiliary mode.

Moreover, at $\phi=\pi/2$, there exists only one specific mode of operation
where the additional entanglement between mode $a_{2}$ and $a_{3}$ vanishes:
when $\mathcal{C}_{13}=1$, the beam-splitter in Eq. (III.1) acts to perfectly
swap all squeezing from mode $a_{3}$ to $a_{1}$ and swap all the uncorrelated
vacuum noise from mode $a_{1}$ to $a_{3}$. The result of this perfect swap is
that $a_{1}$ and $a_{2}$ now form a two-mode squeezed vacuum state, so the
value of $E_{N}^{(1,2)}$ will be equivalent to the value achieved by a TMS, as
seen in Fig. 3. Since modes $a_{2}$ and $a_{3}$ no longer share any squeezed
correlations their measure of entanglement must vanish: $E_{N}^{(2,3)}=0$.

Figure 3: (a) The logarithmic negativity $E_{N}^{(j,k)}$ and (b) the purity
$\mu^{(j,k)}$ of the stationary output states when
$\mathcal{C}_{12}=\mathcal{C}_{13}\mathcal{C}_{23}$ as a function of the loop
phase. The results for the joint states of the outputs of modes $a_{1}$ and
$a_{2}$ (blue) and $a_{2}$ and $a_{3}$ (orange) are shown, for two parameter
choices: “symmetric,” $\mathcal{C}_{23}=\mathcal{C}_{12},\mathcal{C}_{13}=1$
(solid) and “asymmetric,”
$\mathcal{C}_{23}=\mathcal{C}_{12}/2,\mathcal{C}_{13}=2$ (dashed), with
$\mathcal{C}_{12}=0.5$ in both cases. All baths are in the vacuum state. The
dashed black lines at the top of each plot give results for the output state
of a reciprocal two-mode squeezed system with vacuum input (TMSV),
$\mathcal{C}_{13}=\mathcal{C}_{23}=0$.

In order to discuss the behaviour of $E^{(2,3)}_{N}$ in more detail we again
use the Simon-Peres-Horodecki criterion; the output of $a_{2}$ and $a_{3}$
will be separable so long as the following inequality is satisfied:

$\sqrt{\frac{\mathcal{C}_{23}}{\mathcal{C}_{12}\mathcal{C}_{13}}}+\sqrt{\frac{\mathcal{C}_{12}\mathcal{C}_{13}}{\mathcal{C}_{23}}}\leq
2\sin\phi.$ (18)

Referring to Eq. (6), it is evident that this is only satisfied when there is
nonreciprocal scattering with direction $a_{2}\leftarrow a_{3}$. The required
phase here is $\phi=+\pi/2$ which is the opposite phase requirement from Eq.
(16).

It is evident that the degree and direction of nonreciprocal scattering, and
therefore the value of the phase $\phi$, plays a crucial role in the behaviour
of the output-field entanglement, as depicted in Fig. 3. In particular, when
$\mathcal{C}_{12}=\mathcal{C}_{23}$ and $\mathcal{C}_{13}=1$,
$\mathcal{N}^{(1,2)}=1$ and $\mathcal{N}^{(2,3)}=1$, so the scattering
processes for both pairs of modes are perfectly nonreciprocal. $E_{N}^{(1,2)}$
and $E_{N}^{(2,3)}$ reach maximum values for this parameter regime, where they
can both realise the same entangling power of a reciprocal two-mode squeezer;
however, the maxima are achieved at different values of the phase. In addition
this is also the only operating regime where both $E_{N}^{(1,2)}$ and
$E_{N}^{(2,3)}$ reach the absolute minimum value of 0.

Furthermore, in this parameter regime, where perfect swapping is also
observed, the scattering of modes $a_{1}$ and $a_{3}$ is impedance matched in
both cases (see Appendix C), so the input noise is not reflected in the output
fields. For later convenience, we refer to the regime where
$\mathcal{C}_{13}=1$ as the “symmetric” case, since
$\mathcal{C}_{12}=\mathcal{C}_{23}$. We therefore label the regime where
$\mathcal{C}_{13}\neq 1$ the “asymmetric” case. The degree of impedance on
mode $a_{1}$ is presented in Fig. 2(a), where we see that the reflection of
mode $a_{1}$ vanishes, i.e., $||\mathbf{S}_{11}||=0$, only in the symmetric
case.

The circuit decomposition also allows for a heuristic explanation of the
behaviour of the marginal purities, seen in Fig. 3(b). Since the initial
covariance matrix is $\mathbf{V}_{\mathrm{in}}=\mathbf{I}_{6}/2$ for vacuum
inputs, the marginal purities for the input states will be $\mu^{(1,2)}=1$ and
$\mu^{(2,3)}=1$. These purities will remain unchanged in the output state
provided that the corresponding two-mode block of the output covariance matrix
can be reached by a symplectic transformation. This follows since the
determinant of a symplectic transformation is $\mathrm{det}(\mathbf{S})=1$ and
therefore
$\mathrm{det}(\mathbf{S}\mathbf{O}\mathbf{S}^{T})=\mathrm{det}(\mathbf{O})$
for any matrix $\mathbf{O}$.

Since the covariance matrix for $\phi=-\pi/2$ [see Eq. (15)] simply describes
a two-mode squeezing interaction between modes $a_{2}$ and $a_{3}$, the
marginal purity for their outputs will always remain the same,
$\mu^{(2,3)}=1$.. On the other hand, the marginal purity of the output state
between modes $a_{1}$ and $a_{2}$ is below one, $\mu^{(1,2)}<1$, in this case,
as the two-mode block $\mathbf{V}^{(1,2)}$ cannot be reached by any
combination of symplectic transformations. For the opposite phase
$\phi=+\pi/2$, this same reasoning applies to the two-mode block
$\mathbf{V}^{(2,3)}$, so that now $\mu^{(2,3)}<1$. However, $\mu^{(1,2)}$ is
still generally less than one, with a notable exception: the symmetric case,
where $\mathbf{V}^{(1,2)}$ becomes equivalent to the covariance matrix for a
two-mode squeezed vacuum state. In fact, the condition that saturates the
marginal purity $\mu^{(1,2)}=1$ is the same condition for which only the
output of modes $a_{1}$ and $a_{2}$ are entangled, as $E_{N}^{(2,3)}=0$ [see
Eq. (18)].

We have therefore observed and explained how entanglement arises in a
nonreciprocal system. Crucially, we find that the points of perfect
nonreciprocity play a special role in maximising the achievable output
entanglement. We also note how the special symmetric case allows for the
purity of the output states to be maximised. We are now in a position to
analyse the role of nonreciprocity in entanglement generation in the presence
of thermal fluctuations.

##  V Output entanglement in the presence of thermal noise

Figure 4: The logarithmic negativity (solid) and purity (dashed) of NRL output
states as a function of the loop phase in the presence of thermal noise. Here,
$n^{\rm th}_{1}=10$ and $n^{\rm th}_{2}=n^{\rm th}_{3}=0$. The results shown
are for the symmetric case where $\mathcal{C}_{12}=\mathcal{C}_{23}=0.5$ and
$\mathcal{C}_{13}=1$, describing the joint output states of modes $a_{1}$ and
$a_{2}$ (blue), and modes $a_{2}$ and $a_{3}$ (orange). We also show the
results for a two-mode squeezed state where $n_{1}^{\mathrm{th}}=10$ and
$n_{2}^{\mathrm{th}}=0$ (green). The dashed black line indicates both the
logarithmic negativity and purity for a two-mode squeezed vacuum state (TMSV).

Thermal noise is an unwanted feature when attempting to generate entanglement.
In the case of a reciprocal two-mode squeezer, thermal noise incident on one
or both modes will only serve to degrade the logarithmic negativity. While
this can be overcome by increasing the strength of the two-mode squeezing
interaction (e.g., in parametric cQED, using stronger pump strengths), the
same is not true for the purity of the generated output state. More precisely,
the purity for a two-mode squeezed system where the thermal noise at the
inputs for both modes is $n_{1}^{\mathrm{th}}$ and $n_{2}^{\mathrm{th}}$ is
given by $\mu^{(1,2)}=1/(2n_{1}^{\mathrm{th}}+1)(2n_{2}^{\mathrm{th}}+1)$,
which is independent of the degree of squeezing.

One might expect that reciprocally coupling an auxiliary cold mode to the hot
modes of interest would help in mitigating this impact. However, while such a
coupling can reduce the internal occupation of the hot modes, it is unable to
continuously route thermal inputs in a specified direction: away from the
propagating output fields. Combining this cold auxiliary mode with
nonreciprocity enables unidirectional scattering of coherent input signals
which extends to the routing of thermal fluctuations, while also allowing for
the output of the target modes to be entangled. While it is also possible for
a three-mode reciprocal system to route thermal fluctuations in a similar
manner, no entanglement can be generated between the outputs of the target
modes (see Appendix F).

Figure 5: The entanglement between the stationary output of the modes $a_{1}$
and $a_{2}$ as measured by the logarithmic negativity, as a function of the
strength of the thermal noise input $n_{1}^{\rm th}$ on mode $a_{1}$. We
compare the results from an open two-mode squeezed system (TMS, grey) with our
nonreciprocal three-mode loop (NRL, blue). In both cases, the squeezing
interaction between modes $a_{1}$ and $a_{2}$ has cooperativity
$\mathcal{C}_{12}=0.5$. The symmetric case is shown for the nonreciprocal
loop, so $\mathcal{C}_{23}=0.5$, $\mathcal{C}_{13}=1$, and $\phi=+\pi/2$. The
thermal occupation of mode $a_{2}$ is taken to be $n^{\rm th}_{2}=0$ (solid)
and $n^{\rm th}_{2}=1$ (dashed). Mode $a_{3}$ is always taken to have vacuum
input. The dashed black lines correspond to the logarithmic negativity for a
TMS where $n_{1}^{\rm th}\rightarrow\infty,n_{2}^{\rm th}=0$ and $n_{1}^{\rm
th}\rightarrow\infty,n_{2}^{\rm th}=1$.

###  V.1 Rerouting thermal fluctuations using nonreciprocity

The nonreciprocal loop provides a way for us to avoid these effects on the
state shared between the output of modes $a_{1}$ and $a_{2}$ by setting the
parameters to the symmetric case, where the scattering of modes $a_{1}$ and
$a_{3}$ are both impedance matched. Provided that the input of mode $a_{3}$ is
vacuum noise, if we set $\phi=+\pi/2$, thermal noise in the input of mode
$a_{1}$ can be rerouted to the output of mode $a_{3}$. Due to the presence of
thermal noise, the initial covariance matrix is no longer proportional to the
identity matrix and so the output covariance matrix has the following form:

$\mathbf{V}=\mathbf{U}^{(1,3)}_{\mathrm{swap}}\mathbf{R}^{(2,3)}\mathbf{V}_{\mathrm{in}}\big{(}\mathbf{R}^{(2,3)}\big{)}^{T}\\!\big{(}\mathbf{U}^{(1,3)}_{\mathrm{swap}}\big{)}^{T},$
(19)

where $\mathbf{U}_{\rm swap}^{(1,3)}$ is the beam-splitter operation from Eq.
(III.1) when $\mathcal{C}_{13}=1$, which describes a perfect swap. However,
the circuit description shown in Fig. 2 still holds: modes $a_{2}$ and $a_{3}$
are entangled and the subsequent beam-splitter acts as a perfect swap between
modes $a_{1}$ and $a_{3}$. Since $n_{1}^{\mathrm{th}}\neq 0$, the output of
mode $a_{3}$ will receive the unwanted thermal noise while the output of mode
$a_{1}$ forms a two-mode squeezed state with mode $a_{2}$. Provided that the
input for mode $a_{2}$ is also vacuum noise, then the shared state for the
output of modes $a_{1}$ and $a_{2}$ will be a two-mode squeezed vacuum state
with maximum purity and entanglement, unaffected by the value of
$n_{1}^{\mathrm{th}}$, as seen in Fig. 4 (for details of the scatttering
behaviour, see Appendix D).

At this point of operation, a complementary circuit description can be
obtained using the right polar decomposition of the scattering matrix instead,
which also takes on a simple form,

$-\mathbf{S}=\mathbf{R}^{(1,2)}\,\mathbf{U}^{(1,3)}_{\mathrm{swap}}\text{ when
}\phi=+\frac{\pi}{2}$ (20)

where $\mathbf{R}^{(1,2)}$ is a two-mode squeezing operation between modes
$a_{1}$ and $a_{2}$:

$\mathbf{R}^{(1,2)}\leftrightarrow\exp\left[2i\,\mathrm{artanh}\Big{(}\\!\sqrt{\mathcal{C}_{23}}\Big{)}\left(\hat{a}^{\dagger}_{1}\hat{a}^{\dagger}_{2}+\hat{a}_{1}\hat{a}_{2}\right)\right].$
(21)

Then, we can write the covariance matrix described by Eq. (19) in an
equivalent form:

$\mathbf{V}=\mathbf{R}^{(1,2)}\mathbf{U}^{(1,3)}_{\mathrm{swap}}\mathbf{V}_{\mathrm{in}}\big{(}\mathbf{U}^{(1,3)}_{\mathrm{swap}}\big{)}^{T}\\!\big{(}\mathbf{R}^{(1,2)}\big{)}^{T}.$
(22)

Here, the action of the beam-splitter on $\mathbf{V}_{\rm in}$ can be seen
explicitly: it swaps the thermal noise from mode $a_{1}$ with the input from
mode $a_{3}$, which is vacuum noise. This is followed by a two-mode squeezer
acting directly to entangle modes $a_{1}$ and $a_{2}$, creating a state with
maximum purity and entanglement. It is important to note that mode $a_{1}$ is
not cooled using this scheme, and that the nonreciprocal loop only allows for
the thermal noise to be rerouted so as to not appear in the output field.

If we tune $\phi$ away from operating points of perfect nonreciprocity and
impedance matching, the entanglement and purity of modes $a_{1}$ and $a_{2}$
degrade as before. However, comparing Fig. 4 and the previous results when
only considering vacuum inputs (see Fig. 3) we note that the degradation is
more pronounced when the input to mode $a_{1}$ is thermal. This observation
further highlights the importance of nonreciprocity in implementing perfect
swaps of thermal inputs.

Finally, while thermal noise in one mode is detrimental to both entanglement
and purity, the effects are compounded when both modes contain some thermal
noise input. Fig. 5 demonstrates the effects of incident thermal noise on both
modes of an entangled pair. The nonreciprocal loop has the benefit that
regardless of the amount of thermal noise incident on mode $a_{2}$, the
thermal noise incident on mode $a_{1}$ is always swapped away in the output.
The usual two-mode squeezed state, on the other hand, will experience extra
degradation of the entanglement and purity as $n_{1}^{\mathrm{th}}$ increases
for even relatively small values of $n_{2}^{\mathrm{th}}$. In addition,
provided that the noise incident on mode $a_{1}$ is at a higher temperature
than the noise incident on mode $a_{3}$, $n_{1}^{\rm{th}}>n_{3}^{\rm{th}}$,
then the NRL will always improve the fidelity of the entangled output state
for modes $a_{1}$ and $a_{2}$ when compared to the usual TMS at the same
interaction strength. This is true even when other internal (unmonitored) loss
channels are present (see Appendix G).

###  V.2 Entangling the output fields of two hot modes

Figure 6: The schematic for a potential entanglement-swapping scheme to
entangle the output fields of two hot modes. The scheme begins with two
nonreciprocal loops, each comprising a hot primary mode (red), a secondary
mode (green) that we wish to optimally entangle with the hot mode, and a cold
auxiliary mode (blue). As with the setup covered in this paper, the hot and
cold auxiliary modes are coupled via a beam-splitter interaction. The other
interactions must then be two-mode squeezers. To entangle the output of both
hot modes, the output state of the secondary modes is passed through a 50/50
beam-splitter and then measured. The displacement of the output for the hot
modes is then conditioned on the measurements; their outputs will then be
entangled.

The nonreciprocal loop allows for the output fields of one hot mode and one
cold mode to be entangled with maximum purity. As seen in Fig. 5, when using
this system to realise such an entangled state between the output fields of
two hot modes, the logarithmic negativity is reduced, since it is only
possible to reroute the thermal noise from one of the input modes.

However, it is possible to entangle the outputs of two hot modes if we use two
nonreciprocal loops and employ the entanglement-swapping protocol [63, 64]. A
possible setup is shown in Fig. 6, with two realizations of the NRL ($A$ and
$B$, respectively). Two hot modes, $a_{1}$ and $b_{1}$, are coupled to cold
auxiliary modes, $a_{3}$ and $b_{3}$, respectively, with secondary modes,
$a_{2}$ and $b_{2}$, completing the corresponding loops. Both loops are
operated under symmetric configurations and are thus parameterised by a single
two-mode squeezing cooperativity each ($\mathcal{C}_{A}$ and $\mathcal{C}_{B}$
respectively). In light of our previous results, we also operate at points of
perfect nonreciprocal scattering and impedance matching. The loop phase is
assumed to be tuned such that the thermal noise incident on the hot mode is
routed to the output of the cold auxiliary mode. Provided that the input for
each secondary mode is also in the vacuum state, the output for the hot and
secondary modes from each loop will be optimally entangled as a two-mode
squeezed vacuum state with maximum purity.

In the swapping protocol, the outputs from the secondary modes are mixed via a
beam-splitter and then measured. The outputs for the hot modes are then
displaced, conditioned on the result of these measurements. Perfect
application of this entanglement-swapping protocol would allow for the output
fields of the two hot modes to be combined to produce another two-mode
squeezed vacuum state with purity $\mu=1$. The entanglement for the resulting
state is then [64]

$E_{N}=2r\text{ where
}\tanh(r)=\frac{4\sqrt{\mathcal{C}_{A}\mathcal{C}_{B}}}{(1+\mathcal{C}_{A})(1+\mathcal{C}_{B})}$
(23)

which is independent of the thermal noise incident on both hot modes. We
therefore see that the entanglement of the propagating fields from hot modes
$a_{1}$ and $b_{1}$, and the purity of the generated flying states, can be
rendered robust against their thermal inputs using nonreciprocity.

##  VI Conclusions and outlook

In this paper, we analyse how engineered nonreciprocity in quantum systems
influences their ability to entangle propagating fields incident on their
constituent quantum modes. This requires analysing the role of nonreciprocity
beyond asymmetric scattering for signal routing - a notion that can be defined
completely classically - and, in particular, exploring its influence on
quantum correlations and steady-state entanglement of fields. Using a minimal
system consisting of a two-mode squeezer where each mode is coupled to a third
auxiliary mode in a closed-loop configuration, we show that the generated
entanglement of output fields depends strongly on the direction of
nonreciprocal scattering. It is not a priori obvious that it should be
possible to entangle outputs from two modes for which signal flow is only
unidirectional. However, we show that this is indeed possible, given the right
configuration of the system.

To explain this somewhat surprising entanglement behaviour, we develop a
heuristic picture that is based on a polar decomposition of the scattering
matrix. This description maps nonreciprocal scattering to sequential Gaussian
circuit operations, including pairwise beam splitters and, more importantly,
two-mode squeezers which are necessary for generating entanglement. This
picture helps us to explain a second key result: that engineered
nonreciprocity can be used to reroute thermal fluctuations from a hot
propagating mode toward the output of the cold auxiliary mode (via the beam-
splitter component), while simultaneously allowing the entanglement of
propagating output fields (via the two-mode squeezer component). This renders
output-field entanglement much more robust to thermal fluctuations when
compared to a reciprocal two-mode squeezing interaction.

Our work is relevant to the generation of stationary entanglement of itinerant
low-frequency modes, where thermal occupations can be appreciable even at
cryogenic temperatures. Our analysis also brings to light the possible uses of
nonreciprocity in entanglement generation. With recent interest in
multipartite entanglement in quantum systems of increasing scale, our work
invites the exploration of whether engineered nonreciprocity can be a useful
resource in improving robustness of multipartite entanglement in low-frequency
modes.

##  VII Acknowledgements

We thank C. M. Wilson and C. W. Sandbo Chang for the useful discussions. A.M.
and L.O. acknowledge funding by the Deutsche Forschungsgemeinschaft through
the project CRC 910 and the Emmy Noether program (Grant No. ME 4863/1-1). S.K.
acknowledges funding by the Israeli Science Foundation (Grant No. 1794/21) and
the Israeli Innovation Authority (Grant No. 73754). We acknowledge support by
the KIT-Publication Fund of the Karlsruhe Institute of Technology.

##  Appendix A Experimental implementation in parametric cQED

In this appendix section we detail how the model we consider in this work, Eq.
(3), can be realised with standard techniques in parametric cQED.

###  A.1 Circuit Lagrangian

The minimal system we require consists of three modes coupled via a single
nonlinear mixing element. The circuit Lagrangian for such a system is given
by:

$\displaystyle\mathbb{L}=\sum_{m=1,2,3}\\!\left(\frac{1}{2}C_{m}\dot{\tilde{\Phi}}_{m}^{2}-\frac{\tilde{\Phi}_{m}^{2}}{2L_{m}}\right)-U_{J}(\\{\tilde{\Phi}_{m}\\})$
(24)

where $U_{J}(\\{\tilde{\Phi}_{m}\\})$ describes the energy due to the
inductance of superconducting circuit elements incorporating Josephson
junctions. Expressing $U_{J}(\\{\tilde{\Phi}_{m}\\})$ using a Taylor expansion
around an equilbrium point, $\\{\tilde{\Phi}_{m}\\}=0$, gives the following
form:

$\displaystyle
U_{J}=\sum_{mn}\tilde{c}^{(2)}_{mn}\tilde{\Phi}_{m}\tilde{\Phi}_{n}+\sum_{mnr}\tilde{c}^{(3)}_{mnr}\tilde{\Phi}_{m}\tilde{\Phi}_{n}\tilde{\Phi}_{r}+\ldots$
(25)

where

$\displaystyle\tilde{c}^{(2)}_{mn}$
$\displaystyle=\frac{1}{2!}\frac{\partial^{2}U_{\rm
NL}}{\partial\tilde{\Phi}_{m}\partial\tilde{\Phi}_{n}}\Bigg{|}_{\\{\tilde{\Phi}_{m}\\}\to
0},$ (26a) $\displaystyle\tilde{c}^{(3)}_{mnr}$
$\displaystyle=\frac{1}{3!}\frac{\partial^{3}U_{\rm
NL}}{\partial\tilde{\Phi}_{m}\partial\tilde{\Phi}_{n}\partial\tilde{\Phi}_{r}}\Bigg{|}_{\\{\tilde{\Phi}_{m}\\}\to
0}.$ (26b)

The Josephson term provides quadratic contributions that serve to renormalise
the bare linear modes, as we will soon show. Earlier work [22] has considered
the time-modulation of these quadratic terms as a means of generating tunable
Gaussian interactions.

While this is a viable approach, in this appendix section we will instead
explore means of realising the NRL by coherent pumping of nonlinear mixing
terms, demonstrated in several recent works [65, 66]. The lowest-order
nonlinear contribution is defined as

$\displaystyle U_{\rm
NL}=\sum_{mnr}\tilde{c}^{(3)}_{mnr}\tilde{\Phi}_{m}\tilde{\Phi}_{n}\tilde{\Phi}_{r}.$
(27)

Note that higher-order nonlinear contributions cannot always be neglected, as
they can contribute non-rotating terms, such as Kerr terms. However, recent
work [46, 47] has shown methods to engineer Kerr-free nonlinear potentials
which can substantially suppress such contributions.

Now, by introducing the vector of mode fluxes
$\bm{\tilde{\Phi}}=(\tilde{\Phi}_{1},\tilde{\Phi}_{2},\tilde{\Phi}_{3})$, the
circuit Lagrangian of Eq. (24) can be written in the compact matrix form

$\displaystyle\mathbb{L}=\frac{1}{2}\dot{\bm{\tilde{\Phi}}}{\\!\\!~{}}^{T}\mathbf{C}\dot{\bm{\tilde{\Phi}}}-\frac{1}{2}\bm{\tilde{\Phi}}^{T}\mathbf{L}^{-1}\bm{\tilde{\Phi}}-U_{\rm
NL}(\\{\tilde{\Phi}_{m}\\})$ (28)

where we have introduced the capacitance matrix $\mathbf{C}$ and inductance
matrix $\mathbf{L}$, whose matrix elements are defined as:

$\displaystyle\mathbf{C}_{nm}=C_{n}\delta_{nm},~{}\mathbf{L}^{-1}_{nm}=L_{n}^{-1}\delta_{nm}+\tilde{c}^{(2)}_{nm}$
(29)

where $\delta_{nm}$ is the Kronecker $\delta$ function. It now proves useful
to diagonalise the quadratic part of the circuit Lagrangian. To do this, we
first write the Euler-Lagrange equations for the quadratic part of the
Lagrangian,

$\displaystyle\mathbf{C}\ddot{\bm{\tilde{\Phi}}}=-\mathbf{L}^{-1}\tilde{\bm{\Phi}}.$
(30)

The Euler-Lagrange equations allow us to introduce dimensionless circuit
eigenmodes $\\{\bm{\varphi}^{(j)}\\}$ with eigenfrequencies
$\\{\omega_{j}\\}$, which satisfy the generalised eigenproblem

$\displaystyle\omega_{j}^{2}\mathbf{C}\bm{\varphi}^{(j)}=\mathbf{L}^{-1}\bm{\varphi}^{(j)}$
(31)

and obey the orthogonality relations

$\displaystyle\bm{\varphi}^{(j)T}\mathbf{C}\bm{\varphi}^{(k)}=C_{S}\,\delta_{jk},~{}\bm{\varphi}^{(j)T}\mathbf{L}^{-1}\bm{\varphi}^{(k)}=C_{S}\,\omega_{k}^{2}\delta_{jk}$
(32)

where $C_{S}$ is a scaling capacitance introduced to ensure that the circuit
eigenmodes are dimensionless. Note that $C_{S}\,\omega_{j}^{2}$ has units of
inverse inductance, as required by the above expression.

The eigenmodes form a complete basis, which allows us to expand the circuit
flux variables in terms of flux variables corresponding to the eigenmodes,

$\displaystyle\bm{\tilde{\Phi}}=\sum_{j}\Phi_{j}\bm{\varphi}^{(j)},~{}\dot{\bm{\tilde{\Phi}}}=\sum_{j}\dot{\Phi}_{j}\bm{\varphi}^{(j)},$
(33)

Substituting the above into Eq. (28) and making use of the orthogonality of
the circuit eigenmodes, we immediately arrive at the Lagrangian,

$\displaystyle\mathbb{L}=\frac{1}{2}\sum_{j}C_{S}\dot{\Phi}_{j}^{2}-\frac{1}{2}\sum_{j}C_{S}\omega_{j}^{2}\Phi_{j}^{2}-U_{\rm
NL}(\\{\Phi_{j}\\})$ (34)

where, using Eq. (27), the nonlinear contribution in terms of eigenmode fluxes
takes the form:

$\displaystyle U_{\rm
NL}(\\{\Phi_{j}\\})=\sum_{mnr}\tilde{c}^{(3)}_{mnr}\sum_{jkl}\varphi_{m}^{(j)}\varphi_{n}^{(k)}\varphi_{r}^{(l)}\Phi_{j}\Phi_{k}\Phi_{l}$
$\displaystyle=\sum_{jkl}\\!\left[\sum_{mnr}\tilde{c}^{(3)}_{mnr}\varphi_{m}^{(j)}\varphi_{n}^{(k)}\varphi_{r}^{(l)}\right]\\!\\!\Phi_{j}\Phi_{k}\Phi_{l}\equiv\sum_{jkl}c^{\prime}_{jkl}\Phi_{j}\Phi_{k}\Phi_{l}.$
(35)

In this diagonalised form, the conjugate momenta are simply given by
$\frac{\partial\mathbb{L}}{\partial\dot{\Phi}_{j}}=Q_{j}=C_{S}\dot{\Phi}_{j}$,
so that the Legendre transformation defining the Hamiltonian,
$\mathbb{H}=\sum_{j}Q_{j}\dot{\Phi}_{j}-\mathbb{L}$, may be carried out
straightforwardly. We finally obtain the circuit Hamiltonian,

$\displaystyle\mathbb{H}=\frac{1}{2}\sum_{j}\frac{Q_{j}^{2}}{C_{S}}+\frac{1}{2}\sum_{j}C_{S}\omega_{j}^{2}\Phi_{j}^{2}+U_{\rm
NL}(\\{\Phi_{j}\\}).$ (36)

We can now obtain the quantum Hamiltonian by promoting the canonical position
and momenta $\Phi_{j},Q_{j}$ to operators, then writing them in the basis of
creation and annihilation operators satisfying the usual commutation relations
$[\hat{d}_{j},\hat{d}^{\dagger}_{k}]=\delta_{jk}$,

$\displaystyle\hat{\Phi}_{j}$
$\displaystyle=\sqrt{\frac{1}{2C_{S}\omega_{j}}}\left(\hat{d}_{j}+\hat{d}_{j}^{\dagger}\right)\equiv\sqrt{\frac{\hbar
Z_{j}}{2}}\left(\hat{d}_{j}+\hat{d}_{j}^{\dagger}\right)$ (37a)
$\displaystyle\hat{Q}_{j}$
$\displaystyle=-i\sqrt{\frac{C_{S}\omega_{j}}{2}}\left(\hat{d}_{j}-\hat{d}_{j}^{\dagger}\right)\equiv-i\sqrt{\frac{\hbar}{2Z_{j}}}\left(\hat{d}_{j}-\hat{d}_{j}^{\dagger}\right)$
(37b)

where $Z_{j}$ defines the effective impedance of the $j$th circuit eigenmode.

In terms of mode creation and annihilation operators, the system Hamiltonian
then takes the form

$\displaystyle\hat{\mathbb{H}}/\hbar=\sum_{j}\omega_{j}\hat{d}_{j}^{\dagger}\hat{d}_{j}+\sum_{jkl}c_{jkl}(\hat{d}_{j}+\hat{d}_{j}^{\dagger})(\hat{d}_{k}+\hat{d}_{k}^{\dagger})(\hat{d}_{l}+\hat{d}_{l}^{\dagger})$
(38)

where $c_{jkl}=\sqrt{\hbar Z_{j}Z_{k}Z_{l}/8}~{}c^{\prime}_{jkl}$.

Note that throughout this derivation we have made no assumptions regarding the
frequencies $\\{\omega_{j}\\}$ of the modes. In particular, Eq. (38) can
describe a system with a single “hot” low-frequency (RF) mode coupled to two
cooler higher frequency (microwave) modes, as analysed in the main text for
rerouting of thermal fluctuations.

###  A.2 NRL Hamiltonian using parametric drives

We now introduce parametric drives to Eq. (38) that allow us to realise the
NRL proposed in the main text. In what follows, we set $\hbar=1$ and consider
three coherent pump tones at frequencies $\\{\nu_{j}\\}$ with generally
complex pump amplitudes $|\alpha_{j}|e^{-i\phi_{pj}}$, applied to the $j$th
mode. In the steady-state, the pump tones lead to a coherent displacement of
the system modes which can be accounted for via a standard (and exact)
displacement transformation

$\displaystyle\hat{d}_{j}=|\alpha_{j}|e^{-i\phi_{pj}}e^{-i\nu_{j}t}+\hat{a}_{j}e^{-i\phi_{j}},$
(39)

where $\phi_{j}$ defines an arbitrary phase for the $j$th mode operator that
leaves commutation relations unchanged. We simultaneously move to an
interaction picture via the following unitary transformation:

$\displaystyle\hat{\mathcal{U}}=\prod_{j}{\rm
exp}\left[-i\omega_{j}\hat{a}_{j}^{\dagger}\hat{a}_{j}t\right],$ (40)

to remove trivial evolution due to the bare mode Hamiltonians. The transformed
Hamiltonian is then given by:

$\displaystyle\hat{H}_{\rm
NRL}=\hat{\mathcal{U}}\hat{\mathbb{H}}\,\hat{\mathcal{U}}^{\dagger}-i\,\dot{\hat{\mathcal{U}}}\,\hat{\mathcal{U}}^{\dagger}$
(41)

We are interested in the quadratic terms in the transformed Hamiltonian
(linear terms only lead to displacements while higher-order nonlinear terms
are suppressed for strong enough pump amplitudes). These take the form:

$\displaystyle\hat{H}_{\rm
NRL}=\sum_{jkl}3c_{jkl}(|\alpha_{j}|e^{-i\phi_{pj}}e^{-i\nu_{j}t}+h.c.)\times$
$\displaystyle(\hat{a}_{k}e^{-i\phi_{k}}e^{-i\omega_{k}t}+\hat{a}_{k}^{\dagger}e^{i\phi_{k}}e^{i\omega_{k}t})(\hat{a}_{l}e^{-i\phi_{l}}e^{-i\omega_{l}t}+\hat{a}_{l}^{\dagger}e^{i\phi_{l}}e^{i\omega_{l}t})$
(42)

Through an appropriate choices of pump frequencies, we can now make specific
interactions in Eq. (42) resonant in the interaction picture. To realise the
NRL, we choose the following pump frequencies:

$\displaystyle\nu_{1}=\omega_{2}+\omega_{3},~{}\nu_{2}=\omega_{3}-\omega_{1},~{}\nu_{3}=\omega_{1}+\omega_{2}$
(43)

It is straightforward to see from Eq. (42) that $\nu_{1}$ now resonantly pumps
a two-mode squeezing interaction between modes $a_{2}$ and $a_{3}$, while
$\nu_{3}$ pumps a two-mode squeezing interaction between the modes $a_{2}$ and
$a_{3}$. In contrast, $\nu_{2}$ pumps a beam-splitter interaction between
modes $a_{1}$ and $a_{3}$. All other interaction terms will be rapidly
oscillating in this frame and can be neglected within the rotating wave
approximation (RWA), which will be discussed later in this section.

Under this choice of pump frequencies, the system Hamiltonian takes the form:

$\displaystyle\hat{H}_{\rm
NRL}=3c_{123}\Big{(}|\alpha_{3}|e^{-i(\phi_{p3}-\phi_{1}-\phi_{2})}\hat{a}_{1}^{\dagger}\hat{a}^{\dagger}_{2}+$
$\displaystyle|\alpha_{2}|e^{-i(\phi_{p2}-\phi_{3}+\phi_{1})}\hat{a}_{3}^{\dagger}\hat{a}_{1}+|\alpha_{1}|e^{-i(\phi_{p1}-\phi_{2}-\phi_{3})}\hat{a}_{2}^{\dagger}\hat{a}_{3}^{\dagger}+h.c.\Big{)}$
(44)

where we have retained only the resonant terms.

We can finally address the question of the phases. The gauge phases
$\\{\phi_{j}\\}$ can be freely chosen to absorb dependencies on pump phases.
In particular, requiring $\phi_{p3}=\phi_{1}+\phi_{2}$,
$\phi_{p2}=\phi_{3}-\phi_{1}$ removes the phase dependence of the first two
terms above, and yields the final system Hamiltonian

$\displaystyle\hat{H}_{\rm
NRL}=\left(g_{12}\hat{a}_{1}^{\dagger}\hat{a}^{\dagger}_{2}+g_{13}\hat{a}_{3}^{\dagger}\hat{a}_{1}+g_{23}e^{i\phi}\hat{a}_{2}^{\dagger}\hat{a}_{3}^{\dagger}\right)+h.c.$
(45)

where:

$\displaystyle
g_{12}=3c_{123}|\alpha_{3}|,~{}g_{13}=3c_{123}|\alpha_{2}|,~{}g_{23}=3c_{123}|\alpha_{1}|$
(46)

and the remaining loop phase is given by $\phi=\phi_{p3}+\phi_{p2}-\phi_{p1}$.
In particular, the loop phase is independent of $\\{\phi_{j}\\}$ and hence
cannot be “gauged away.” Furthermore, $\phi$ can be tuned by rotating the
phase of any of the incident pump fields. We have thus obtained Eq. (3) of the
main text, with interaction strengths tunable via pump amplitudes, having
started from a very general but minimal three-mode circuit Lagrangian with a
single three-wave nonlinear mixing element.

From the above derivation, it is clear that the validity of Eq. (3) hinges on
the RWA. To verify that the RWA is valid, we consider a concrete three-mode
system with frequencies $\omega_{j}/(2\pi)\in\\{0.5,7.5,10.5\\}~{}{\rm GHz}$;
these are representative of modes realised in very recent RF quantum optics
experiments, such as Ref. [54]. Then, the desired pump frequencies are given
by $\nu_{j}/(2\pi)\in\\{18,10,8\\}~{}{\rm GHz}$. Terms we have dropped in
arriving at Eq. (3) include undesired beam splitter interactions pumped by
difference frequencies
$f_{-}/(2\pi)\in\\{\omega_{2}-\omega_{1},\omega_{3}-\omega_{2}\\}/(2\pi)=\\{7,3\\}~{}{\rm
GHz}$, and undesired amplifying interactions pumped by sum frequencies
$f_{+}/(2\pi)\in\\{\omega_{1}+\omega_{3},2\omega_{1},2\omega_{2},2\omega_{3}\\}/(2\pi)=\\{11,1,15,21\\}~{}{\rm
GHz}$. For this particular choice of $\\{\omega_{j}\\}$, the desired pump
frequencies are hence at least $500~{}{\rm Mhz}$ (and typically further) away
from all undesired pumping frequencies and individual mode frequencies. As a
result, all non-resonant interaction terms will be oscillating with a
frequency of at least $500~{}{\rm MHz}$. The RWA is valid provided this
oscillation frequency is much larger than the interaction strengths $g_{jk}$.
Cooperativities $\mathcal{C}_{jk}\sim O(1)$ considered in the main text imply
$g_{jk}\simeq\kappa_{j},\kappa_{k}$, where $\kappa_{j}$ is the $j$th mode
decay rate. For typical decay rates $\kappa_{j}\sim O(1)~{}{\rm MHz}$ in cQED,
this means the fast rotating frequency is at least $O(100)$ times larger than
the required interaction strengths. Hence the RWA can be expected to hold.

##  Appendix B Heisenberg-Langevin equations of motion

We work in the quadrature basis, which consists of the position and momentum
quadrature operators,
$\hat{X}_{j}=(\hat{a}^{\dagger}_{j}+\hat{a}_{j})/\sqrt{2}$ and
$\hat{P}_{j}=i(\hat{a}^{\dagger}_{j}-\hat{a}_{j})/\sqrt{2}$, respectively. The
canonical commutation relations have the usual form
$[{\hat{X}_{j}},{\hat{P}_{k}}]=i\delta_{jk}$. Defining the vector of
quadrature operators for the three-mode system
$\vec{R}=(\hat{X}_{1},\hat{P}_{1},\hat{X}_{2},\hat{P}_{2},\hat{X}_{3},\hat{P}_{3})$,
we can write the Heisenberg-Langevin equations as the follows

$\frac{d}{dt}\vec{R}(t)=\mathbf{M}\vec{R}(t)-\sqrt{\bm{\kappa}}\vec{R}_{\mathrm{in}}(t)$
(47)

where $\mathbf{M}$ is a time independent dynamical matrix,
$\bm{\kappa}=\mathrm{diag}(\kappa_{1},\kappa_{1},\kappa_{2},\kappa_{2},\kappa_{3},\kappa_{3})$
is a diagonal matrix of the mode damping rates, and $\vec{R}_{\mathrm{in}}$ is
the vector of input noise operators in the quadrature basis. The correlators
of the elements of the matrix $\vec{R}_{\mathrm{in}}$ have the following form:

$\displaystyle\left\langle\hat{X}_{j,\mathrm{in}}(t)\hat{X}_{k,\mathrm{in}}(t^{\prime})\right\rangle=\delta_{jk}\left(n^{\rm
th}_{j}+\frac{1}{2}\right)\delta(t-t^{\prime})$
$\displaystyle\left\langle\hat{P}_{j,\mathrm{in}}(t)\hat{P}_{k,\mathrm{in}}(t^{\prime})\right\rangle=\delta_{jk}\left(n^{\rm
th}_{j}+\frac{1}{2}\right)\delta(t-t^{\prime})$
$\displaystyle\left\langle\hat{X}_{j,\mathrm{in}}(t)\hat{P}_{k,\mathrm{in}}(t^{\prime})\right\rangle=\delta_{jk}\frac{i}{2}\delta(t-t^{\prime}).$
(48)

These are Gaussian white noise processes, and so have a mean of zero. The
dynamical matrix for the system can be written as follows:

$\mathbf{M}=\begin{pmatrix}-\displaystyle{\frac{\kappa_{1}}{2}}\mathbf{I}&-g_{12}\mathbf{X}&g_{13}\mathbf{J}\\\\[6.99997pt]
-g_{12}\mathbf{X}&-\displaystyle{\frac{\kappa_{2}}{2}}\mathbf{I}&\begin{subarray}{c}\displaystyle{g_{23}(\sin\phi\mathbf{Z}}\\\\[1.99997pt]
\displaystyle{-\cos\phi\mathbf{X})}\end{subarray}&\\\\[10.00002pt]
g_{13}\mathbf{J}&\begin{subarray}{c}\displaystyle{g_{23}(\sin\phi\mathbf{Z}}\\\\[1.99997pt]
\displaystyle{-\cos\phi\mathbf{X})}\end{subarray}&-\displaystyle{\frac{\kappa_{3}}{2}}\mathbf{I}\end{pmatrix}$
(49)

where we have written the matrix in block form using the following $2\times 2$
matrices

$\mathbf{X}=\begin{pmatrix}0&1\\\ 1&0\end{pmatrix}\hskip
11.38109pt\mathbf{Z}=\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix}\hskip
11.38109pt\mathbf{J}=\begin{pmatrix}0&1\\\ -1&0\end{pmatrix}$ (50)

while $\mathbf{I}$ is the $2\times 2$ identity matrix and $\mathbf{0}$ the
$2\times 2$ zero matrix. The linear Heisenberg-Langevin equations of motion,
Eq. (47), can be transformed to frequency space and written in the compact
form

$-i\omega\vec{R}[\omega]=\mathbf{M}\vec{R}[\omega]-\sqrt{\bm{\kappa}}\vec{R}_{\mathrm{in}}[\omega].$
(51)

This linear algebraic system can be straightforwardly solved,

$\vec{R}[\omega]=-(i\omega+\mathbf{M})^{-1}\sqrt{\bm{\kappa}}\vec{R}_{\mathrm{in}}[\omega].$
(52)

Quantum input-output theory [57] relates the output field quadratures to the
input and system fields:

$\vec{R}_{\mathrm{out}}[\omega]=\vec{R}_{\mathrm{in}}[\omega]+\sqrt{\bm{\kappa}}\vec{R}[\omega].$
(53)

Using Eq. (52) then allows us to express the output fields entirely in terms
of the input fields, which defines the scattering matrix $\mathbf{S}[\omega]$,
as introduced in Eq. (4) of the main text. The covariance matrix for the
output modes in frequency space may be written as

$\displaystyle\mathbf{V}[\omega]=$
$\displaystyle\frac{1}{2}\int^{\infty}_{\infty}\\!\\!\\!\left\langle\vec{R}_{\mathrm{out}}[\omega]\vec{R}_{\mathrm{out}}[\omega^{\prime}]^{T}+\vec{R}_{\mathrm{out}}[\omega^{\prime}]\vec{R}_{\mathrm{out}}[\omega]^{T}\right\rangle
d\omega^{\prime}$ (54)

where we take the outer product of the $\vec{R}_{\mathrm{out}}$ vectors. The
above can be rewritten in terms of the input fields and the scattering matrix
using Eq. (4), and the correlators of the input fields may be calculated using
the frequency space version of the correlators from Eq. (B), where the Fourier
transform replaces $\delta(t-t^{\prime})$ with
$\delta(\omega+\omega^{\prime})$. Evaluating the integral in Eq. (54) will
then yield Eq. (10) from the main text.

##  Appendix C Block form of full scattering matrix

Using Eq. (52) in conjunction with Eq. (53) allows us to obtain the full,
frequency-dependent scattering matrix of the three-mode system. However, as
discussed in the main text, we are typically interested in scattering
properties at $\omega=0$. Furthermore, we have also discussed how the simplest
intuitive scattering behaviour can be analysed in the symmetric, impedance-
matched case defined by $\mathcal{C}_{12},\mathcal{C}_{23}\equiv\mathcal{C}$,
and $\mathcal{C}_{13}=1$. In this case, $\mathbf{S}[0]\equiv\mathbf{S}$ is
determined entirely by $\mathcal{C}$ and the loop phase $\phi$, and takes the
form:

$\displaystyle\mathbf{S}=\begin{pmatrix}\mathbf{S}_{11}&\mathbf{S}_{12}&\mathbf{S}_{13}\\\
\mathbf{S}_{21}&\mathbf{S}_{22}&\mathbf{S}_{23}\\\
\mathbf{S}_{31}&\mathbf{S}_{32}&\mathbf{S}_{33}\\\ \end{pmatrix}\equiv
D(\phi)\begin{pmatrix}\mathbf{F}(\phi)\cos\phi&\mathbf{B}(\phi)(1+\sin\phi)&\mathbf{A}(\phi)\\\
\mathbf{B}(\phi-\pi)(1-\sin\phi)&-(1-\mathcal{C}^{2})\mathbf{I}+\mathbf{\bar{F}}(\phi)\cos\phi&\mathbf{C}(\phi)\\\
\mathbf{A}(-\phi)&\mathbf{\bar{C}}(\phi)&\mathbf{F}(\phi)\cos\phi\\\
\end{pmatrix}.$ (55)

Here $D(\phi)=(1-\mathcal{C})^{2}+\mathcal{C}^{2}\cos^{2}\phi$ is an overall
multiplicative factor that does not influence the nonreciprocity of
scattering. The diagonal terms describing reflections take the form:

$\displaystyle\mathbf{F}(\phi)=\mathcal{C}^{2}\cos(\phi)\,\mathbf{I}-(1-\mathcal{C})\mathcal{C}\,\mathbf{J}$
$\displaystyle\mathbf{\bar{F}}(\phi)=\mathcal{C}^{2}\cos(\phi)\,\mathbf{I}+2\mathcal{C}\,\mathbf{J}.$
(56)

The off-diagonal terms that describe transmission between modes $a_{1}$ and
$a_{2}$ take the form:

$\displaystyle\mathbf{B}(\phi)=\,$
$\displaystyle\sqrt{\mathcal{C}}\frac{\cos\phi(1+\mathcal{C}\sin\phi)}{1+\sin\phi}\mathbf{Z}$
$\displaystyle-\sqrt{\mathcal{C}}(2\mathcal{C}-\mathcal{C}\sin\phi-1)\,\mathbf{X}.$
(57)

The interaction between modes $a_{1}$ and $a_{3}$ is a beam-splitter and is
compactly described by a single $\phi$-dependent matrix $\mathbf{A}(\phi)$:

$\displaystyle\mathbf{A}(\phi)=\,\,$
$\displaystyle\mathcal{C}^{2}\cos\phi(1+\sin\phi)\mathbf{I}$
$\displaystyle+\left[-\mathcal{C}^{2}\cos^{2}\phi-(1-\mathcal{C})(1+\mathcal{C}\sin\phi)\right]\mathbf{J}.$
(58)

Finally, the interaction between modes $a_{2}$ and $a_{3}$ is described by:

$\displaystyle\mathbf{C}(\phi)$
$\displaystyle=-\sqrt{\mathcal{C}}(1+\sin\phi)(1-\mathcal{C}\sin\phi)\mathbf{Z}$
$\displaystyle~{}~{}~{}+\sqrt{\mathcal{C}}\cos\phi(1-2\mathcal{C}-\mathcal{C}\sin\phi)\mathbf{X}$
$\displaystyle\mathbf{\bar{C}}(\phi)$
$\displaystyle=+\sqrt{\mathcal{C}}(1-\sin\phi)(1+\mathcal{C}\sin\phi)\mathbf{Z}$
$\displaystyle~{}~{}~{}+\sqrt{\mathcal{C}}\cos\phi(1-2\mathcal{C}+\mathcal{C}\sin\phi)\mathbf{X}.$
(59)

From Eq. (55), it is now straightforward to read off conditions for specific
desired scattering properties. For example, impedance matching of mode $a_{1}$
demands $||\mathbf{S}_{11}||=0$, which clearly requires that $\cos\phi=0$, and
hence $\phi=\pm\pi/2$. Similarly, perfect nonreciprocal scattering between
modes $a_{1}$ and $a_{2}$ as defined in Eq. (5) of the main text,
$\mathcal{N}^{(1,2)}=1$, clearly requires that $(1\pm\sin\phi)=0$, which again
implies that $\phi=\pm\pi/2$. These are the conditions shown in Fig. 2.

##  Appendix D Comparison of the scattering properties of the NRL and TMS

When the NRL is optimised to allow swapping of thermal noise from mode $a_{1}$
to mode $a_{3}$, we use the impedance-matched case from the previous section,
additionally setting the phase to $\phi=\pi/2$. In the quadrature basis the
steady-state scattering matrix will have the following form:

$\mathbf{S}_{\mathrm{NRL}}=\begin{pmatrix}\mathbf{0}&\displaystyle{\frac{2\sqrt{\mathcal{C}}}{1-\mathcal{C}}}\mathbf{X}&-\displaystyle{\frac{1+\mathcal{C}}{1-\mathcal{C}}}\mathbf{J}\\\\[6.00006pt]
\mathbf{0}&-\displaystyle{\frac{1+\mathcal{C}}{1-\mathcal{C}}}\mathbf{I}&-\displaystyle{\frac{2\sqrt{\mathcal{C}}}{1-\mathcal{C}}}\mathbf{Z}&\\\\[6.00006pt]
-\mathbf{J}&\mathbf{0}&\mathbf{0}\end{pmatrix}.$ (60)

Replacing the cooperativity with the following squeezing parameter
$r=\mathrm{artanh}[2\sqrt{\mathcal{C}}/(1+\mathcal{C})]$ we can rewrite the
above scattering matrix as follows:

$\mathbf{S}_{\mathrm{NRL}}\equiv\begin{pmatrix}\mathbf{0}&\sinh(r)\mathbf{X}&-\cosh(r)\mathbf{J}\\\
\mathbf{0}&-\cosh(r)\mathbf{I}&-\sinh(r)\mathbf{Z}&\\\
-\mathbf{J}&\mathbf{0}&\mathbf{0}\end{pmatrix}.$ (61)

This makes clear the behaviour of the system at this point of nonreciprocity.
Modes $a_{1}$ and $a_{2}$ are independent of the input noise on mode $a_{1}$,
which only appears in mode $a_{3}$ showing how the noise is rerouted there.
Meanwhile, modes $a_{2}$ and $a_{3}$ share some squeezed correlations.

We can also calculate the steady-state scattering matrix for an open system
with a TMS Hamiltonian given by
$ig(\hat{a}^{\dagger}_{1}\hat{a}^{\dagger}_{2}-\hat{a}_{1}\hat{a}_{2})$
(choosing the TMS phase to be zero as in Eq. (1) is not appropriate here; this
comes from the polar decomposition):

$\mathbf{S}_{\mathrm{TMS}}=\begin{pmatrix}-\displaystyle{\frac{1+\mathcal{C}}{1-\mathcal{C}}}\mathbf{I}&-\displaystyle{\frac{2\sqrt{\mathcal{C}}}{1-\mathcal{C}}}\mathbf{Z}\\\\[6.00006pt]
-\displaystyle{\frac{2\sqrt{\mathcal{C}}}{1-\mathcal{C}}}\mathbf{Z}&-\displaystyle{\frac{1+\mathcal{C}}{1-\mathcal{C}}}\mathbf{I}\end{pmatrix}.$
(62)

Defining the squeezing in the same way we arrive at the following scattering
matrix:

$\mathbf{S}_{\mathrm{TMS}}\equiv\begin{pmatrix}-\cosh(r)\mathbf{I}&-\sinh(r)\mathbf{Z}\\\
-\sinh(r)\mathbf{Z}&-\cosh(r)\mathbf{I}\end{pmatrix}.$ (63)

While the covariance matrix of the TMS and modes $a_{1}$ and $a_{2}$ in the
NRL will be identical, the scattering behaviour is markedly different. This is
expected given the circuit decomposition for the NRL; since the squeezing is
swapped from mode $a_{3}$ to mode $a_{1}$ by a beam-splitter, the quadratures
are rotated during this swap in a manner that cannot be replicated in a TMS
alone.

##  Appendix E Stability conditions

We provide below the Routh-Hurwitz stability criterion for the three-mode loop
at the phases $\phi=\pm\pi/2$:

$\displaystyle 0<\kappa_{1}+\kappa_{2}+\kappa_{3}$ $\displaystyle
0<1-\mathcal{C}_{12}+\mathcal{C}_{13}-\mathcal{C}_{23}$ $\displaystyle
0<1-\frac{\mathcal{C}_{12}}{(1+\kappa_{3}/\kappa_{1})(1+\kappa_{3}/\kappa_{2})}+\frac{\mathcal{C}_{13}}{(1+\kappa_{2}/\kappa_{1})(1+\kappa_{2}/\kappa_{3})}$
$\displaystyle\hskip
22.76219pt-\frac{\mathcal{C}_{23}}{(1+\kappa_{1}/\kappa_{2})(1+\kappa_{1}/\kappa_{3})}.$
(64)

Away from these phases there are more conditions which must be met, and the
conditions in general take on a much more complicated form. Provided the
cooperativities and dissipation rates are chosen appropriately the system can
be stable for all values of the loop phase.

In case we also apply the condition
$\mathcal{C}_{12}=\mathcal{C}_{13}\mathcal{C}_{23}$ to make the system
nonreciprocal, the second listed condition takes on a much simpler form

$0<(1-\mathcal{C}_{23})(1+\mathcal{C}_{13}).$ (65)

The system is naturally stable for all choices of the beam-splitter
cooperativity, and is therefore limited by the free two-mode squeezing
cooperativity $\mathcal{C}_{23}<1$ which is identical to the stability
criterion for an open two-mode squeezer. The two-mode squeezing cooperativity
fixed by the nonreciprocity condition, $\mathcal{C}_{12}$, can therefore grow
quite large with the system remaining stable. Again, provided the dissipation
rates are chosen correctly, the other stability conditions can be satisfied as
well.

##  Appendix F Thermal rerouting prevents entanglement in a reciprocal system

In order to reroute thermal noise away from the output of mode $a_{1}$ in a
three-mode system while still realising entanglement between the outputs of
modes $a_{1}$ and $a_{2}$, it is required that the scattering of mode $a_{1}$
be impedance matched and that the scattering between modes $a_{1}$ and $a_{2}$
be nonreciprocal. It is simple to demonstrate that it is not possible to do
both in a three-mode system where the scattering between modes $a_{1}$ and
$a_{2}$ is reciprocal. We begin with a scattering matrix which can reroute the
thermal excitations away from the output of mode $a_{1}$:

$\mathbf{S}=\begin{pmatrix}\mathbf{0}&\mathbf{0}&\mathbf{S}_{13}\\\
\mathbf{0}&\mathbf{S}_{22}&\mathbf{S}_{23}\\\
\mathbf{S}_{31}&\mathbf{S}_{32}&\mathbf{S}_{33}\end{pmatrix}.$ (66)

Writing the initial covariance matrix as
$\mathbf{V}_{\mathrm{in}}=\mathrm{diag}(\mathbf{V}_{1,\mathrm{in}},\mathbf{V}_{2,\mathrm{in}},\mathbf{V}_{3,\mathrm{in}})$,
the covariance matrix for the output of modes $a_{1}$ and $a_{2}$ is then

$\displaystyle\mathbf{V}^{(1,2)}=$ $\displaystyle\hskip
5.69054pt\begin{pmatrix}\mathbf{S}_{13}\mathbf{V}_{3,\mathrm{in}}\mathbf{S}_{13}^{T}&\mathbf{S}_{13}\mathbf{V}_{3,\mathrm{in}}\mathbf{S}_{23}^{T}\\\
\mathbf{S}_{23}\mathbf{V}_{3,\mathrm{in}}\mathbf{S}_{13}^{T}&\mathbf{S}_{22}\mathbf{V}_{2,\mathrm{in}}\mathbf{S}_{22}^{T}+\mathbf{S}_{23}\mathbf{V}_{3,\mathrm{in}}\mathbf{S}_{23}^{T}\end{pmatrix}$
(67)

which is independent of $\mathbf{V}_{1,\mathrm{in}}$, as desired. Assuming
that Eq. (66) is a valid scattering matrix, it must be symplectic and hence
satisfy the condition
$\mathbf{S}\mathbf{\Omega}\mathbf{S}^{T}=\mathbf{\Omega}$ where
$\mathbf{\Omega}=\mathrm{diag}(\mathbf{J},\mathbf{J},\mathbf{J})$ is the
symplectic form. Using this it is possible to come up with conditions for the
block elements of Eq. (66).

First, it may be determined that
$\mathbf{S}_{13}\mathbf{J}\mathbf{S}_{13}^{T}=\mathbf{J}$ which indicates that
$\mathbf{S}_{13}\in Sp(2,\mathbb{R})$. We also have to satisfy
$\mathbf{S}_{13}\mathbf{J}\mathbf{S}_{23}^{T}=\mathbf{0}$; since
$\mathbf{S}_{13}\mathbf{J}\in Sp(2,\mathbb{R})$ must be invertible it follows
that $\mathbf{S}_{23}^{T}=(\mathbf{S}_{13}\mathbf{J})^{-1}\mathbf{0}$ and so
$\mathbf{S}_{23}=\mathbf{0}$. As a consequence, the off-diagional blocks in
Eq. (67) vanish, indicating that the output of modes $a_{1}$ and $a_{2}$ are
never entangled. It is therefore not possible to realise thermal noise
rerouting and entanglement in a reciprocal three-mode loop.

##  Appendix G Internal losses

In this appendix section we analyse the role of “internal” loss channels that
are not used to direct inputs to, or measure outputs from, the NRL, but still
contribute added noise. These channels can describe unmonitored ports of the
system (for example the undercoupled port of a two-sided cavity) as well as
material losses. To account for these effects, we rewrite the total loss rates
for the system modes as

$\displaystyle\kappa_{j}=\kappa_{j}^{e}+\kappa_{j}^{\rm int}$ (68)

where $\kappa_{j}^{e}$ defines the loss rate via the monitored or external
loss channel, while internal losses are described by $\kappa^{\rm int}_{j}$.
In this way, the total loss rates are the same as those used in our analysis
in the main text. Thus far, our analysis has considered $\kappa_{j}^{\rm
int}=0$. Accounting for couplings to additional loss channels in the case of
non-zero $\kappa_{j}^{\rm int}$, the Heisenberg-Langevin equations can be
written as:

$\displaystyle\frac{d}{dt}\vec{R}=\mathbf{M}\vec{R}(t)-\sqrt{\bm{\kappa}^{e}}\vec{R}_{\rm
in}(t)-\sqrt{\bm{\kappa}^{\rm int}}\vec{R}_{\rm int}(t)$ (69)

where $\bm{\kappa}^{e}$, $\bm{\kappa}^{\rm int}$ are matrices of external and
internal losses respectively, analogous to $\bm{\kappa}$ introduced earlier.
The term $\vec{R}_{\rm int}(t)$ defines additional noise introduced due to
internal or generally unmonitored loss channels, which we again take to be
Gaussian white noise processes with temperature $n_{j}^{\rm th}$ for mode
$a_{j}$.

Figure 7: (a) The entanglement between the stationary output of modes $a_{1}$
and $a_{2}$ and (b) the purity of their joint quantum state as a function of
the ratio of internal losses experienced by each NRL mode. We consider the
symmetric configuration, with $\mathcal{C}=0.5$, and $n^{\rm th}_{1}=10$,
$n^{\rm th}_{2}=n^{\rm th}_{3}=0$.

Transforming to frequency space as before and rearranging, we find:

$\displaystyle\vec{R}[\omega]=(i\omega+\mathbf{M})^{-1}\left(\sqrt{\bm{\kappa}^{e}}\vec{R}_{\rm
in}[\omega]+\sqrt{\bm{\kappa}^{\rm int}}\vec{R}_{\rm int}[\omega]\right).$
(70)

Then, using the modified input-output relations for the monitored modes,

$\vec{R}_{\mathrm{out}}[\omega]=\vec{R}_{\mathrm{in}}[\omega]+\sqrt{\bm{\kappa}^{e}}\vec{R}[\omega],$
(71)

we arrive at

$\displaystyle\vec{R}_{\rm out}[\omega]$
$\displaystyle=\left[\mathbf{I}+\sqrt{\bm{\kappa}^{e}}(i\omega+\mathbf{M})^{-1}\sqrt{\bm{\kappa}^{e}}\right]\vec{R}_{\rm
in}[\omega]$
$\displaystyle~{}~{}~{}~{}+\sqrt{\bm{\kappa}^{e}}(i\omega+\mathbf{M})^{-1}\sqrt{\bm{\kappa}^{\rm
int}}\vec{R}_{\rm int}[\omega]$
$\displaystyle\equiv\mathbf{S}_{e}[\omega]\vec{R}_{\rm
in}[\omega]+\mathbf{T}_{\rm int}[\omega]\vec{R}_{\rm int}[\omega].$ (72)

The scattering matrix for monitored channels $\mathbf{S}_{e}[\omega]$ is in
general distinct from $\mathbf{S}[\omega]$ in Eq. (4) due to some signal being
lost to unmonitored channels. Furthermore, additional noise contributions
appear at monitored output ports via $\mathbf{T}_{\rm int}$ due to noise
incident from these unmonitored channels. The output covariance matrix can
then be calculated (once again on resonance) using Eq. (54),

$\displaystyle\mathbf{V}=\mathbf{S}_{e}\mathbf{V}_{\rm
in}\mathbf{S}_{e}^{T}+\mathbf{T}_{\rm int}\mathbf{V}_{\rm in}\mathbf{T}_{\rm
int}^{T}$ (73)

which follows since $\vec{R}_{\rm in}(t)$ and $\vec{R}_{\rm int}(t)$ have no
cross-correlations.

The effect of internal losses is therefore encapsulated in the structure of
the matrix $\mathbf{T}_{\rm int}$. With all pairwise interactions turned off,
$\mathcal{C}_{jk}=0~{}\forall~{}j,k$, $\mathbf{T}_{\rm int}$ is a block
diagonal matrix with $\mathbf{T}_{\rm int}^{(j,j)}=-2\sqrt{\kappa_{j}^{\rm
int}\kappa_{j}^{e}/\kappa^{2}_{j}}~{}\mathbf{I}$, where $\mathbf{I}$ is the
$2\times 2$ identity matrix. Once interactions between modes are turned on to
realise the NRL, the precise form of $\mathbf{T}_{\rm int}$ may be modified,
but these diagonal blocks will remain non-zero. As a result, noise from
unmonitored channels on mode $a_{j}$ will appear at the monitored port of the
same mode, and can therefore effect the quantum properties of the monitored
modes.

We now analyse the effect of internal losses on the NRL numerically. We will
focus on modes $a_{1}$ and $a_{2}$, and therefore operate at the point of
perfect nonreciprocal scattering ($\mathcal{N}^{(1,2)}=1$) that maximises
entanglement of the output fields of these modes, by choosing $\phi=+\pi/2$.
We also consider the symmetric NRL configuration with $\mathcal{C}=0.5$. The
resulting parameters are therefore the same as analysed in Fig. 5.

Figure 8: (a) The entanglement between the stationary output of the modes
$a_{1}$ and $a_{2}$ and (b) the purity of their joint quantum state, as a
function of the strength of the thermal input noise $n_{1}^{\mathrm{th}}$ on
mode $a_{1}$. We compare outputs of the NRL to the TMS, both for the case of
no internal losses ($\kappa_{j}^{\rm int}/\kappa_{j}=0~{}\forall~{}j$), and
for the case of a nonzero fixed internal loss ratio ($\kappa_{1}^{\rm
int}/\kappa_{1}=10^{-2},\kappa_{2}^{\rm int}/\kappa_{2}=0$). The value of
$\kappa_{3}^{\rm int}/\kappa_{3}$ does not influence performance.

Fig. 7 shows the output field entanglement $E_{N}^{(1,2)}$ and purity
$\mu^{(1,2)}$ for modes $a_{1}$ and $a_{2}$ as a function of the ratio of
internal losses $\kappa_{j}^{\rm int}/\kappa_{j}$ for each mode. Here,
$n_{1}^{\rm th}=10,n_{2}^{\rm th}=n_{3}^{\rm th}=0$. Note that internal losses
in modes $a_{1}$ and $a_{2}$, whose entanglement is being considered, have a
detrimental impact on performance. This is to be expected: noise incident on
unmonitored channels for modes $a_{1}$ and $a_{2}$ appears directly, and
without any quantum correlations, at the monitored ports via $\mathbf{T}_{\rm
int}^{(1,1)}$ and $\mathbf{T}_{\rm int}^{(2,2)}$ respectively, which are
nonzero. This uncorrelated noise reduces the fidelity and purity of the
entangled state for the output of modes $a_{1}$ and $a_{2}$. Also, noise
incident at a higher temperature (for the internal loss channel of mode
$a_{1}$) is more detrimental.

However, we note that $E_{N}^{(1,2)}$ and $\mu^{(1,2)}$ do not depend on the
internal losses of the auxiliary mode $a_{3}$, introduced in the NRL to enable
nonreciprocal routing. There are two reasons for this effect. First, any noise
from the internal loss channels of mode $a_{3}$ appears uncorrelated only at
the monitored port of $a_{3}$ via $\mathbf{T}_{\rm int}^{(3,3)}$; this does
not influence the output of modes $a_{1}$ and $a_{2}$. Therefore, for this
noise to appear at the output of modes $a_{1}$ and $a_{2}$, it must undergo
the entangling and swapping interactions of the NRL. In this process, this
noise in fact seeds output field entanglement of modes $a_{1}$ and $a_{2}$,
and is evidently not detrimental to the performance of the NRL.

This observation has an important implication: since only the internal losses
of modes $a_{1}$ and $a_{2}$ are important, the NRL can outperform a TMS given
the same internal losses and at the same interaction strength. To demonstrate
this, in Fig. 8 we plot $E_{N}^{(1,2)}$ and $\mu^{(1,2)}$ for the NRL and TMS
when $\mathcal{C}=0.5$, as a function of thermal input noise $n_{1}^{\rm th}$
with $n_{2}^{\rm th}=1,n_{3}^{\rm th}=0$. For simplicity, we choose only the
dominant internal loss rate to be nonzero, $\kappa_{1}^{\rm
int}/\kappa_{1}=0.01,\kappa_{2}^{\rm int}/\kappa_{2}=0$. The plotted curves
have no dependence on the value of $\kappa_{3}^{\rm int}/\kappa_{3}$. For
completeness, we also show the case of zero internal losses that was analysed
in the main text. Clearly, while internal losses on the system modes are
always detrimental, and have greater impact at higher thermal noise inputs,
the NRL always provides higher entanglement fidelity and increased purity than
the TMS. The only constraint on the NRL, then, is one we have already
identified as necessary in the main text: the thermal bath fluctuations of the
auxiliary mode must be cooler than those of the hot mode in the NRL,
$n_{3}^{\rm th}<n_{1}^{\rm th}$.

## References

  * Horodecki  _et al._ [2009] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Reviews of Modern Physics 81, 865 (2009), publisher: American Physical Society.
  * Kimble [2008] H. J. Kimble, The quantum internet, Nature 453, 1023 (2008), number: 7198 Publisher: Nature Publishing Group.
  * Cirac  _et al._ [1997] J. I. Cirac, P. Zoller, H. J. Kimble, and H. Mabuchi, Quantum State Transfer and Entanglement Distribution among Distant Nodes in a Quantum Network, Physical Review Letters 78, 3221 (1997), publisher: American Physical Society.
  * Acín  _et al._ [2007] A. Acín, J. I. Cirac, and M. Lewenstein, Entanglement percolation in quantum networks, Nature Physics 3, 256 (2007), number: 4 Publisher: Nature Publishing Group.
  * Eichler  _et al._ [2012] C. Eichler, C. Lang, J. M. Fink, J. Govenius, S. Filipp, and A. Wallraff, Observation of Entanglement between Itinerant Microwave Photons and a Superconducting Qubit, Physical Review Letters 109, 240501 (2012), publisher: American Physical Society.
  * Gely  _et al._ [2019] M. F. Gely, M. Kounalakis, C. Dickel, J. Dalle, R. Vatré, B. Baker, M. D. Jenkins, and G. A. Steele, Observation and stabilization of photonic Fock states in a hot radio-frequency resonator, Science 363, 1072 (2019), publisher: American Association for the Advancement of Science.
  * Li and Gröblacher [2020] J. Li and S. Gröblacher, Stationary quantum entanglement between a massive mechanical membrane and a low frequency LC circuit, New Journal of Physics 22, 063041 (2020), publisher: IOP Publishing.
  * Steele  _et al._ [2021] G. Steele, M. Gely, I. Rodrigues, D. Bothner, M. Kounalakis, C. Dickel, J. Dalle, R. Vatré, B. Baker, and M. Jenkins, Photonics at radio frequency: Quantum sensing of radio-frequency photons and vibrations, in  _Optical and Quantum Sensing and Precision Metrology_, Vol. 11700 (SPIE, 2021) p. 117001B.
  * Abdi  _et al._ [2015] M. Abdi, P. Tombesi, and D. Vitali, Entangling two distant non-interacting microwave modes, Annalen der Physik 527, 139 (2015), publisher: John Wiley & Sons, Ltd.
  * Vitali  _et al._ [2007] D. Vitali, S. Gigan, A. Ferreira, H.R. Böhm, P. Tombesi, A. Guerreiro, V. Vedral, A. Zeilinger, and M. Aspelmeyer, Optomechanical entanglement between a movable mirror and a cavity field, Physical Review Letters 98, 030405 (2007).
  * Genes  _et al._ [2008] C. Genes, A. Mari, P. Tombesi, and D. Vitali, Robust entanglement of a micromechanical resonator with output optical fields, Physical Review A 78, 032316 (2008), publisher: American Physical Society.
  * Wang  _et al._ [2015] Y.-D. Wang, S. Chesi, and A. A. Clerk, Bipartite and tripartite output entanglement in three-mode optomechanical systems, Physical Review A 91, 013807 (2015), publisher: American Physical Society.
  * Bienfait  _et al._ [2019] A. Bienfait, K. J. Satzinger, Y. Zhong, H.-S. Chang, M.-H. Chou, C. R. Conner, É. Dumur, J. Grebel, G. A. Peairs, R. G. Povey, _et al._ , Phonon-mediated quantum state transfer and remote qubit entanglement, Science 364, 368 (2019).
  * Barzanjeh  _et al._ [2019] S. Barzanjeh, E. S. Redchenko, M. Peruzzo, M. Wulf, D. P. Lewis, G. Arnold, and J. M. Fink, Stationary entangled radiation from micromechanical motion, Nature 570, 480 (2019), number: 7762 Publisher: Nature Publishing Group.
  * Chen  _et al._ [2020] J. Chen, M. Rossi, D. Mason, and A. Schliesser, Entanglement of propagating optical modes via a mechanical interface, Nature Communications 11, 943 (2020).
  * Pan  _et al._ [2001] J.-W. Pan, C. Simon, Č. Brukner, and A. Zeilinger, Entanglement purification for quantum communication, Nature 410, 1067 (2001).
  * Duan  _et al._ [2000] L.-M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, Entanglement purification of gaussian continuous variable quantum states, Phys. Rev. Lett. 84, 4002 (2000).
  * Franzen  _et al._ [2006] A. Franzen, B. Hage, J. DiGuglielmo, J. Fiurás̆ek, and R. Schnabel, Experimental demonstration of continuous variable purification of squeezed states, Phys. Rev. Lett. 97, 150505 (2006).
  * Kamal  _et al._ [2011] A. Kamal, J. Clarke, and M. Devoret, Noiseless non-reciprocity in a parametric active device, Nature Physics 7, 311 (2011).
  * Sliwa  _et al._ [2015] K.M. Sliwa, M. Hatridge, A. Narla, S. Shankar, L. Frunzio, R.J. Schoelkopf, and M.H. Devoret, Reconfigurable josephson circulator/directional amplifier, Physical Review X 5, 041020 (2015).
  * Ruesink  _et al._ [2016] F. Ruesink, M.-A. Miri, A. Alu, and E. Verhagen, Nonreciprocity and magnetic-free isolation based on optomechanical interactions, Nature communications 7, 1 (2016).
  * Lecocq  _et al._ [2017] F. Lecocq, L. Ranzani, G.A. Peterson, K. Cicak, R.W. Simmonds, J.D. Teufel, and J. Aumentado, Nonreciprocal microwave signal processing with a field-programmable josephson amplifier, Physical Review Applied 7, 024028 (2017).
  * Peterson  _et al._ [2017a] G. A. Peterson, F. Lecocq, K. Cicak, R. W. Simmonds, J. Aumentado, and J. D. Teufel, Demonstration of efficient nonreciprocity in a microwave optomechanical circuit, Physical Review X 7, 031001 (2017a).
  * Chapman  _et al._ [2017] B. J. Chapman, E. I. Rosenthal, J. Kerckhoff, B. A. Moores, L. R. Vale, J.A.B. Mates, G. C. Hilton, K. Lalumiere, A. Blais, and K.W. Lehnert, Widely tunable on-chip microwave circulator for superconducting quantum circuits, Physical Review X 7, 041043 (2017).
  * Miri  _et al._ [2017] M.-A. Miri, F. Ruesink, E. Verhagen, and A. Alù, Optical nonreciprocity based on optomechanical coupling, Physical Review Applied 7, 064014 (2017).
  * Fang  _et al._ [2017a] K. Fang, J. Luo, A. Metelmann, M. H. Matheny, F. Marquardt, A. A. Clerk, and O. Painter, Generalized non-reciprocity in an optomechanical circuit via synthetic magnetism and reservoir engineering, Nature Physics 13, 465 (2017a).
  * Ando  _et al._ [2020] F. Ando, Y. Miyasaka, T. Li, J. Ishizuka, T. Arakawa, Y. Shiota, T. Moriyama, Y. Yanase, and T. Ono, Observation of superconducting diode effect, Nature 584, 373 (2020).
  * Shen  _et al._ [2016] Z. Shen, Y.-L. Zhang, Y. Chen, C.-L. Zou, Y.-F. Xiao, X.-B. Zou, F.-W. Sun, G.-C. Guo, and C.-H. Dong, Experimental realization of optomechanically induced non-reciprocity, Nature Photonics 10, 657 (2016).
  * Bernier  _et al._ [2017a] N. R. Bernier, L. D. Toth, A. Koottandavida, M. A. Ioannou, D. Malz, A. Nunnenkamp, A. Feofanov, and T. Kippenberg, Nonreciprocal reconfigurable microwave optomechanical circuit, Nature communications 8, 1 (2017a).
  * Barzanjeh  _et al._ [2017] S. Barzanjeh, M. Wulf, M. Peruzzo, M. Kalaee, P. Dieterle, O. Painter, and J. M. Fink, Mechanical on-chip microwave circulator, Nature communications 8, 1 (2017).
  * MercierdeLépinay  _et al._ [2019] L. Mercier de Lépinay, E. Damskägg, C. F. Ockeloen-Korppi, and M. A. Sillanpää, Realization of directional amplification in a microwave optomechanical device, Physical Review Applied 11, 034027 (2019).
  * Xu  _et al._ [2019] H. Xu, L. Jiang, A. Clerk, and J. Harris, Nonreciprocal control and cooling of phonon modes in an optomechanical system, Nature 568, 65 (2019).
  * Dong  _et al._ [2015] C.-H. Dong, Z. Shen, C.-L. Zou, Y.-L. Zhang, W. Fu, and G.-C. Guo, Brillouin-scattering-induced transparency and non-reciprocal light storage, Nature communications 6, 1 (2015).
  * Kim  _et al._ [2015] J. Kim, M. C. Kuzyk, K. Han, H. Wang, and G. Bahl, Non-reciprocal brillouin scattering induced transparency, Nature Physics 11, 275 (2015).
  * Metelmann and Clerk [2015] A. Metelmann and A. A. Clerk, Nonreciprocal photon transmission and amplification via reservoir engineering, Phys. Rev. X 5, 021025 (2015).
  * Ranzani and Aumentado [2015] L. Ranzani and J. Aumentado, Graph-based analysis of nonreciprocity in coupled-mode systems, New Journal of Physics 17, 023024 (2015).
  * Metelmann and Clerk [2017] A. Metelmann and A.A. Clerk, Nonreciprocal quantum interactions and devices via autonomous feedforward, Physical Review A 95, 013837 (2017).
  * Wang and Clerk [2013] Y.-D. Wang and A. A. Clerk, Reservoir-Engineered Entanglement in Optomechanical Systems, Physical Review Letters 110, 253601 (2013), publisher: American Physical Society.
  * Li  _et al._ [2021] Y. Li, W. Li, Han Tiancheng, Z. Xu, J. Li, B. Li, Fan Shanhui, and Q. Cheng-Wei, Transforming heat transfer with thermal metamaterials and devices, Nature Reviews. Materials 6, 488 (2021).
  * Terraneo  _et al._ [2002] M. Terraneo, M. Peyrard, and G. Casati, Controlling the Energy Flow in Nonlinear Lattices: A Model for a Thermal Rectifier, Physical Review Letters 88, 094302 (2002), publisher: American Physical Society.
  * Li  _et al._ [2004] B. Li, L. Wang, and G. Casati, Thermal Diode: Rectification of Heat Flux, Physical Review Letters 93, 184301 (2004), publisher: American Physical Society.
  * Joulain  _et al._ [2016] K. Joulain, J. Drevillon, Y. Ezzahri, and J. Ordonez-Miranda, Quantum thermal transistor, Physical review letters 116, 200601 (2016).
  * Barzanjeh  _et al._ [2018] S. Barzanjeh, M. Aquilina, and A. Xuereb, Manipulating the flow of thermal noise in quantum devices, Physical Review Letters 120, 060601 (2018).
  * Ordonez-Miranda  _et al._ [2019] J. Ordonez-Miranda, Y. Ezzahri, J. A. Tiburcio-Moreno, K. Joulain, and J. Drevillon, Radiative thermal memristor, Physical review letters 123, 025901 (2019).
  * Jalas  _et al._ [2013] D. Jalas, A. Petrov, M. Eich, W. Freude, S. Fan, Z. Yu, R. Baets, M. Popović, A. Melloni, J. D. Joannopoulos, M. Vanwolleghem, C. R. Doerr, and H. Renner, What is — and what is not — an optical isolator, Nature Photonics 7, 579 (2013).
  * Frattini  _et al._ [2017] N. E. Frattini, U. Vool, S. Shankar, A. Narla, K. M. Sliwa, and M. H. Devoret, 3-wave mixing Josephson dipole element, Applied Physics Letters 110, 222603 (2017), publisher: American Institute of Physics.
  * Sivak  _et al._ [2019] V.V. Sivak, N.E. Frattini, V.R. Joshi, A. Lingenfelter, S. Shankar, and M.H. Devoret, Kerr-Free Three-Wave Mixing in Superconducting Quantum Circuits, Physical Review Applied 11, 054060 (2019), publisher: American Physical Society.
  * Walls and Milburn [2008] D. Walls and G. J. Milburn, eds.,  _Quantum Optics_ (Springer, Berlin, Heidelberg, 2008).
  * Liu  _et al._ [2020] G. Liu, X. Cao, T.-C. Chien, C. Zhou, P. Lu, and M. Hatridge, Noise reduction in qubit readout with a two-mode squeezed interferometer, arXiv preprint arXiv:2007.15460 (2020).
  * Qiu  _et al._ [2022] J. Y. Qiu, A. Grimsmo, K. Peng, B. Kannan, B. Lienhard, Y. Sung, P. Krantz, V. Bolkhovsky, G. Calusine, D. Kim, _et al._ , Broadband squeezed microwaves and amplification with a josephson traveling-wave parametric amplifier, arXiv preprint arXiv:2201.11261 (2022).
  * Peterson  _et al._ [2017] G.A. Peterson, F. Lecocq, K. Cicak, R.W. Simmonds, J. Aumentado, and J.D. Teufel, Demonstration of efficient nonreciprocity in a microwave optomechanical circuit, Physical Review X 7, 031001 (2017).
  * Bernier  _et al._ [2017b] N. R. Bernier, L. D. Tóth, A. Koottandavida, M. A. Ioannou, D. Malz, A. Nunnenkamp, A. K. Feofanov, and T. J. Kippenberg, Nonreciprocal reconfigurable microwave optomechanical circuit, Nature Communications 8, 604 (2017b), number: 1 Publisher: Nature Publishing Group.
  * Fang  _et al._ [2017b] K. Fang, J. Luo, A. Metelmann, M. H. Matheny, F. Marquardt, A. A. Clerk, and O. Painter, Generalized non-reciprocity in an optomechanical circuit via synthetic magnetism and reservoir engineering, Nature Physics 13, 465 (2017b), number: 5 Publisher: Nature Publishing Group.
  * Rodrigues  _et al._ [2021] I. C. Rodrigues, D. Bothner, and G. A. Steele, Cooling photon-pressure circuits into the quantum regime, Science Advances 7, eabg6653 (2021).
  * Rodrigues  _et al._ [2022] I. C. Rodrigues, G. A. Steele, and D. Bothner, Parametrically enhanced interactions and nonreciprocal bath dynamics in a photon-pressure Kerr amplifier, Science Advances 8, eabq1690 (2022), publisher: American Association for the Advancement of Science.
  * Weedbrook  _et al._ [2012] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Gaussian quantum information, Reviews of Modern Physics 84, 621 (2012).
  * Gardiner and Collett [1985] C. W. Gardiner and M. J. Collett, Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation, Physical Review A 31, 3761 (1985).
  * Braunstein [2005] S. L. Braunstein, Squeezing as an irreducible resource, Physical Review A 71, 055801 (2005).
  * Gosson [2011] M. A. d. Gosson,  _Symplectic Methods in Harmonic Analysis and in Mathematical Physics_ , 1st ed. (Birkhäuser, Basel, 2011).
  * Simon [2000] R. Simon, Peres-horodecki separability criterion for continuous variable systems, Physical Review Letters 84, 2726 (2000).
  * Vidal and Werner [2002] G. Vidal and R. F. Werner, Computable measure of entanglement, Physical Review A 65, 032314 (2002).
  * Plenio [2005] M. B. Plenio, Logarithmic negativity: A full entanglement monotone that is not convex, Phys. Rev. Lett. 95, 090503 (2005).
  * Hoelscher-Obermaier and van Loock [2011] J. Hoelscher-Obermaier and P. van Loock, Optimal gaussian entanglement swapping, Physical Review A 83, 012319 (2011).
  * van Loock [2002] P. van Loock, Quantum communication with continuous variables, Fortschritte der Physik: Progress of Physics 50, 1177 (2002).
  * Chien  _et al._ [2020] T.-C. Chien, O. Lanes, C. Liu, X. Cao, P. Lu, S. Motz, G. Liu, D. Pekker, and M. Hatridge, Multiparametric Amplification and Qubit Measurement with a Kerr-free Josephson Ring Modulator, Physical Review A 101, 042336 (2020), arXiv:1903.02102 [quant-ph].
  * Zhou  _et al._ [2022] C. Zhou, P. Lu, M. Praquin, T.-C. Chien, R. Kaufman, X. Cao, M. Xia, R. Mong, W. Pfaff, D. Pekker, and M. Hatridge, A modular quantum computer based on a quantum state router, arXiv:2109.06848 [quant-ph] 10.48550/arXiv.2109.06848 (2022).