# Impossibility of bosonic autonomous entanglement engines in the weak-
coupling limit

Bradley Longstaff  Michael G. Jabbour  Jonatan Bohr Brask  Center for
Macroscopic Quantum States (bigQ), Department of Physics, Technical University
of Denmark, 2800 Kongens Lyngby, Denmark

###### Abstract

Entanglement is a fundamental feature of quantum physics and a key resource
for quantum communication, computing and sensing. Entangled states are fragile
and maintaining coherence is a central challenge in quantum information
processing. Nevertheless, entanglement can be generated and stabilised through
dissipative processes. In fact, entanglement has been shown to exist in the
steady state of certain interacting quantum systems subject solely to
incoherent coupling to thermal baths. This has been demonstrated in a range of
bi- and multipartite settings using systems of finite dimension. Here we focus
on the steady state of infinite-dimensional bosonic systems. Specifically, we
consider any set of bosonic modes undergoing excitation-number-preserving
interactions of arbitrary strength and divided between an arbitrary number of
parties that each couple weakly to thermal baths at different temperatures. We
show that a unique steady state is always separable.

##  I Introduction

Entanglement is a central feature of quantum physics and a key resource for
many quantum-information-processing tasks Horodecki  _et al._ (2009), such as
quantum communication, computing and sensing Gisin and Thew (2007);
Giovannetti  _et al._ (2011); Pironio  _et al._ (2016); Pirandola  _et al._
(2020). The generation and stabilisation of entangled states is a challenging
task. As entanglement is generally fragile with respect to loss and noise
Horodecki  _et al._ (2009), considerable effort is usually required to isolate
the system of interest from the environment. Nevertheless, it has been shown
that coupling to an environment can be used to aid entanglement generation
Plenio  _et al._ (1999); Kim  _et al._ (2002); Jakóbczyk (2002); Braun (2002);
Benatti  _et al._ (2003); Burgarth and Giovannetti (2007); Bellomo  _et al._
(2008); Naseem and Müstecaplıoğlu (2022). Steady-state entanglement can also
be generated by dissipation and external driving Diehl  _et al._ (2008);
Verstraete  _et al._ (2009); Kraus  _et al._ (2008); Ticozzi and Viola (2014).
In fact, dissipative steady-state entanglement generation is even possible
without any source of coherence or external control, a setting referred to as
autonomous. This has been demonstrated in a number of different contexts for
finite-dimensional systems Plenio and Huelga (2002); Hartmann  _et al._ (2006,
2007); Žnidarič (2012); Quiroga  _et al._ (2007); Bellomo and Antezza (2013a,
b); Brask  _et al._ (2015); Tavakoli  _et al._ (2018); Man  _et al._ (2019);
Tavakoli  _et al._ (2020); Brask  _et al._ (2022).

Much less is known for infinite-dimensional systems, although they are widely
employed in practice in experiments with quantum optical, optomechanical, and
superconducting setups. In particular, the class of so-called Gaussian states
and operations in bosonic systems are ubiquitous as they are relatively simple
to realise and their theory is well understood Weedbrook  _et al._ (2012).
Gaussian processes, e.g. squeezing and linear interferometers, are used
extensively in continuous-variable quantum information processing and quantum
optics in general. However, it turns out that non-Gaussian resources are
necessary in order to perform various important quantum-information-processing
tasks including, for example, entanglement distillation Eisert  _et al._
(2002); Giedke and Ignacio Cirac (2002); Fiurášek (2002), quantum error
correction Niset  _et al._ (2009), and universal quantum computation Lloyd and
Braunstein (1999); Knill  _et al._ (2001). It is therefore of interest to
determine the combinations of Gaussian and non-Gaussian resources that can
generate steady-state entanglement. In the present work, we demonstrate that
autonomous, steady-state entanglement in bosonic systems with passive Gaussian
interactions is not possible.

We consider a multi-mode bosonic system undergoing quadratic, excitation-
preserving interactions (i.e., passive Gaussian interactions or linear
interferometers) of arbitrary strength, divided between an arbitrary number of
parties. The modes of different parties couple weakly to thermal baths at
different temperatures, with the energy distribution in each bath following
either Bose-Einstein or Fermi-Dirac statistics (both bosonic and spin baths
are allowed). The system-bath interactions are also quadratic but not
necessarily excitation number preserving. We show that in such a setting, no
entanglement can ever be present in the steady state between the parties of
the system.

The paper is organised as follows. In Section II we introduce the model of our
thermal machines. We then present a separability criterion in Section III that
will be satisfied by the steady states. The master equations that describe the
dynamics of our machines are derived in Section IV. We focus on two
situations: a global approach, which is appropriate when the interactions
between subsystems is strong relative to the system-bath couplings; and a
local approximation, which describes a situation where the interactions
between subsystems is weak relative to the system-bath couplings. We then show
that the steady state is separable in both situations. Finally, we conclude in
Section V.

##  II Bosonic autonomous thermal machine

We consider a collection of bosonic modes undergoing excitation-preserving
interactions of arbitrary strength and distributed between an arbitrary number
of parties that each couple weakly to thermal baths at different temperatures
(Fig. 1).

Figure 1: Sketch of the quantum thermal machine. The system consists of $d$
bosonic modes distributed between $N$ subsystems. The system modes are coupled
via excitation-number-preserving interactions. The subsystems are separately
coupled to either bosonic or spin baths. Bosonic baths exchange both energy
and particles with the system while the spin baths exchange only energy. The
system-bath coupling is taken to be weak, while the inter-system couplings can
be of arbitrary strength.

The total system $S$ consists of $d$ bosonic modes with associated creation
and annihilation operators that satisfy the canonical commutation relations
$[\hat{a}_{j},\hat{a}^{\dagger}_{k}]=\delta_{jk}$,
$[\hat{a}_{j},\hat{a}_{k}]=0$ and
$[\hat{a}^{\dagger}_{j},\hat{a}^{\dagger}_{k}]=0$. The corresponding
quadrature operators for each mode are defined as

$\hat{q}_{j}=\frac{1}{\sqrt{2}}\left(\hat{a}^{\dagger}_{j}+\hat{a}_{j}\right),\quad\hat{p}_{j}=\frac{i}{\sqrt{2}}\left(\hat{a}^{\dagger}_{j}-\hat{a}_{j}\right),$
(1)

and fulfil the commutation relations $[\hat{q}_{j},\hat{p}_{k}]=i\delta_{jk}$,
where we use units with $\hbar=1$ throughout.

The Hamiltonian of the system $\hat{H}_{S}$ is taken to be a general stable
($\hat{H}_{S}>0$), quadratic, and excitation-number-preserving Hamiltonian

$\hat{H}_{S}=\sum_{i,j}\mathcal{H}_{ij}\hat{a}^{\dagger}_{i}\hat{a}_{j},$ (2)

where $\mathcal{H}_{ij}$ is a Hermitian matrix. There are $N$ subsystems and
baths, with the $n$th subsystem denoted by $S_{n}$ and coupling to the $n$th
bath $E_{n}$. The total environment is described by the continuous non-
interacting Hamiltonian

$\hat{H}_{E}=\sum_{n=1}^{N}\hat{H}_{E_{n}}=\sum_{n=1}^{N}\int\epsilon_{n}(q)\,\hat{c}^{\dagger}_{n}(q)\hat{c}_{n}(q)\,dq,$
(3)

and we assume that the bath spectra are non negative with $\epsilon_{n}(q)\geq
0$. In general, each bath may be bosonic or consist of spins. For the bosonic
baths, the bath creation and annihilation operators satisfy the canonical
commutation relations

$\displaystyle[\hat{c}_{n}(k),\hat{c}^{\dagger}_{m}(q)]$
$\displaystyle=\delta_{nm}\delta(k-q),$ (4)
$\displaystyle[\hat{c}_{n}(k),\hat{c}_{m}(q)]$
$\displaystyle=[\hat{c}^{\dagger}_{n}(k),\hat{c}^{\dagger}_{m}(q)]=0,$ (5)

while spins in different spin baths commute and spins within the same bath
obey the anticommutation relations

$\displaystyle\\{\hat{c}_{n}(k),\hat{c}^{\dagger}_{n}(q)\\}$
$\displaystyle=\delta(k-q),$ (6)
$\displaystyle\\{\hat{c}_{n}(k),\hat{c}_{n}(q)\\}$
$\displaystyle=\\{\hat{c}^{\dagger}_{n}(k),\hat{c}^{\dagger}_{n}(q)\\}=0.$ (7)

Finally, the system-bath interaction is a quadratic Hamiltonian of the form
$\hat{H}_{SE}=\sum_{n=1}^{N}\hat{V}_{n}$, with

$\hat{V}_{n}=\sum_{j\in S_{n}}\int
g_{j}(q)\left(\hat{a}^{\dagger}_{j}+\eta_{j}\hat{a}_{j}\right)\hat{c}_{n}(q)\,dq+h.c..$
(8)

Here the $j$-sum is taken over the modes of subsystem $S_{n}$, the complex
coefficients $g_{j}(q)$ quantify the strength of each interaction, and the
parameters $\eta_{j}\in\\{0,1\\}$ allow for both excitation-number-conserving
($\eta_{j}=0$) and ‘position-like’ ($\eta_{j}=1$) coupling to the environment
.

The bosonic baths exchange energy and particles with the system, and are thus
characterised by an inverse temperature $\beta_{n}=1/k_{B}T_{n}$ and chemical
potential $\mu_{n}$. On the other hand, the spin baths only exchange energy
with the system. We assume that the environment is initially in the thermal
state

$\hat{\rho}_{E}=\bigotimes_{n=1}^{N}\frac{\mathrm{e}^{-(\hat{H}_{E_{n}}-\mu_{n}\hat{N}_{n})\beta_{n}}}{\tr\left[\mathrm{e}^{-(\hat{H}_{E_{n}}-\mu_{n}\hat{N}_{n})\beta_{n}}\right]},$
(9)

where $\hat{N}_{n}=\int\hat{c}^{\dagger}_{n}(q)\hat{c}_{n}(q)\,dq$ is the
number operator of the $n$th bath and the chemical potential should only be
included for bosonic baths. In this case, the two-point correlation functions
satisfy

$\displaystyle\langle\hat{c}^{\dagger}_{n}(k)\hat{c}_{m}(q)\rangle=p_{n}(\epsilon_{n}(k))\,\delta_{nm}\,\delta(k-q),$
(10)
$\displaystyle\langle\hat{c}_{n}(k)\hat{c}_{m}(q)\rangle=0,\quad\langle\hat{c}^{\dagger}_{n}(k)\hat{c}^{\dagger}_{m}(q)\rangle=0,$
(11)

where $\langle\cdot\rangle$ is the expectation value in the state
$\hat{\rho}_{E}$ and

$p_{n}(\epsilon)=\left[\xi_{n}+\mathrm{e}^{(\epsilon-\mu_{n})\beta_{n}}\right]^{-1}$
(12)

is the Fermi-Dirac distribution for spin baths ($\xi_{n}=1$) or the Bose-
Einstein distribution for bosonic baths ($\xi_{n}=-1$).

##  III A criterion for separability

Our goal is to derive the steady state of the system in the model above.
Before proceeding, however, we first consider a criterion for separability. As
the steady state is found to be Gaussian (see Section IV.3), separability can
be determined from the covariance matrix. By definition, a Gaussian state
$\hat{\rho}$ has a Wigner function of the form

$W(x)=\frac{1}{(2\pi)^{d}}\frac{1}{\sqrt{\det\Sigma}}\mathrm{e}^{-\frac{1}{2}(x-\bar{x})^{T}\Sigma^{-1}(x-\bar{x})}.$
(13)

Here $x=(q_{1},p_{1},\ldots,q_{d},p_{d})\in\mathbb{R}^{2d}$ are canonical
phase-space coordinates, $\bar{x}$ is a vector of first moments
$\bar{x}_{j}=\tr[\hat{x}_{j}\hat{\rho}]$ with
$\hat{x}=(\hat{q}_{1},\hat{p}_{1},\ldots,\hat{q}_{d},\hat{p}_{d})$, and the
covariance matrix $\Sigma$ is a real, symmetric and positive-definite matrix
with elements

$\Sigma_{jk}=\frac{1}{2}\langle\hat{x}_{j}\hat{x}_{k}+\hat{x}_{j}\hat{x}_{k}\rangle-\langle\hat{x}_{j}\rangle\langle\hat{x}_{k}\rangle,$
(14)

where $\langle\cdot\rangle=\tr[\cdot\hat{\rho}]$. The covariance matrix of a
valid quantum state satisfies the Robertson-Schrödinger uncertainty relation

$\Sigma+\frac{i}{2}\Omega\geq 0,$ (15)

where $\Omega$ is the symplectic matrix

$\Omega=\bigoplus_{j=1}^{d}\begin{pmatrix}0&1\\\ -1&0\end{pmatrix}.$ (16)

As shown by Werner and Wolf Werner and Wolf (2001), a bipartite Gaussian state
of systems $A$ and $B$ with covariance matrix $\Sigma_{AB}$ is separable if
and only if there exist valid covariance matrices $\Sigma_{A}$ and
$\Sigma_{B}$ for systems $A$ and $B$ such that
$\Sigma_{AB}\geq\Sigma_{A}\oplus\Sigma_{B}$. We would like to apply this
result to our steady states. However, it turns out that the steady-state
covariance matrices derived below in Sec. IV are more conveniently expressed
in coordinates different to $(q_{1},p_{1},\ldots,q_{d},p_{d})$. Specifically,
in the complex coordinates
$(q_{1}+ip_{1},\ldots,q_{d}+ip_{d},q_{1}-ip_{1},\ldots,q_{d}-ip_{d})/\sqrt{2}$
they have the simple block structure

$\Sigma^{\prime}=\begin{pmatrix}\sigma^{\ast}&0\\\ 0&\sigma\end{pmatrix},$
(17)

where
$\sigma_{jk}=\frac{1}{2}\langle\hat{a}^{\dagger}_{j}\hat{a}_{k}+\hat{a}_{k}\hat{a}^{\dagger}_{j}\rangle-\langle\hat{a}^{\dagger}_{j}\rangle\langle\hat{a}_{k}\rangle$
and the prime denotes the complex form. The transformation to complex
coordinates is achieved by first applying the permutation $P$ that maps
$(q_{1},p_{1},\ldots,q_{d},p_{d})$ to
$(q_{1},\ldots,q_{d},p_{1},\ldots,p_{d})$, and then the unitary matrix

$U=\frac{1}{\sqrt{2}}\begin{pmatrix}I_{d}&iI_{d}\\\
I_{d}&-iI_{d}\end{pmatrix},$ (18)

where $I_{d}$ is the $d\times d$ identity matrix. Thus, the similarity
transformation $UP$ yields $\Sigma^{\prime}=UP\Sigma P^{\dagger}U^{\dagger}$.

From the steady-state covariance matrix (17), we observe that
$\langle\hat{a}_{j}\hat{a}_{k}\rangle-\langle\hat{a}_{j}\rangle\langle\hat{a}_{k}\rangle=0$,
and so the corresponding Gaussian steady state is not squeezed. Separability
of the steady state then follows from the fact that every entangled Gaussian
state is squeezed, see e.g. Wolf  _et al._ (2003). For completeness, we
provide a proof that any $\Sigma^{\prime}$ of the form (17) corresponds to a
$\Sigma$ that satisfies the separability condition across any bipartition. In
complex coordinates, the physicality condition (15) becomes
$\Sigma^{\prime}+\frac{1}{2}K\geq 0$, where

$K=\begin{pmatrix}I_{d}&0\\\ 0&-I_{d}\end{pmatrix},$ (19)

as may be seen by applying the transformation $UP$. Since $\Sigma^{\prime}$ in
(17) comes from a valid steady state, it follows that
$\sigma-\frac{1}{2}I_{d}\geq 0$. Taking the complex conjugate, we also have
$\sigma^{\ast}-\frac{1}{2}I_{d}\geq 0$. These two results then imply that
$\Sigma^{\prime}\geq\frac{1}{2}I_{2d}$. Transforming back to the real
coordinates we obtain $\Sigma\geq\frac{1}{2}I_{2d}$, as the eigenvalues are
preserved under the transformation. For any bipartition of the total system we
have that $I_{2d}=I_{A}\oplus I_{B}$, with $I_{A}$ and $I_{B}$ identity
matrices corresponding to subsystems A and B, and $\frac{1}{2}I$ is always a
valid covariance matrix (the covariance matrix of the vacuum state). Hence, a
Gaussian state with covariance matrix of the form (17) is separable across
every bipartition.

##  IV Open-system dynamics and steady states

In this section we derive master equations for the open-system dynamics of the
setup in Fig. 1. It is shown that the steady state is a Gaussian state with a
covariance matrix of the form (17), and is therefore separable by the argument
of Sec. III. We consider weak system-bath couplings and large baths, and
assume that the initial composite state of the system and baths factorises. In
this case, one can apply a perturbative approach in the system-bath coupling
to derive a master equation for the time evolution of the system. A time-local
(Markovian) description of the system evolution, which does not depend on the
entire history of the system state, is obtained provided that correlations in
the baths decay much faster than variations of the state of the system
$\tilde{\rho}(t)$ in the interaction picture (denoted by a tilde).

Following the standard procedure, one finds (see, e.g., Refs. Rivas and Huelga
(2012); Potts (2019) for an in-depth derivation)

$\displaystyle\frac{d}{dt}\tilde{\rho}(t)$ (20)
$\displaystyle=-\int_{0}^{\infty}du\,\textnormal{tr}_{E}\bigg{(}\left[\hat{H}_{SE}(t),\left[\hat{H}_{SE}(t-u),\tilde{\rho}(t)\otimes\hat{\rho}_{E}\right]\right]\bigg{)}.$

Substituting in the expressions for $\hat{H}_{SE}$ and $\hat{\rho}_{E}$ leads
to

$\displaystyle\frac{d}{dt}\tilde{\rho}(t)=\sum_{n=1}^{N}\sum_{j,k\in
S_{n}}\int_{0}^{\infty}du\,\bigg{(}\tr_{E}\left[\hat{B}^{\dagger}_{jn}(u)\hat{B}_{kn}\hat{\rho}_{E}\right]\left(\hat{A}_{j}^{\dagger}(t-u)\tilde{\rho}(t)\hat{A}_{k}(t)-\hat{A}_{k}(t)\hat{A}_{j}^{\dagger}(t-u)\tilde{\rho}(t)\right)$
(21)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\quad+\tr_{E}\left[\hat{B}_{jn}(u)\hat{B}^{\dagger}_{kn}\hat{\rho}_{E}\right]\left(\hat{A}_{j}(t-u)\tilde{\rho}(t)\hat{A}_{k}^{\dagger}(t)-\hat{A}_{k}^{\dagger}(t)\hat{A}_{j}(t-u)\tilde{\rho}(t)\right)+h.c.\bigg{)},$

where we have introduced the system operators
$\hat{A}_{j}(t)=\hat{a}_{j}(t)+\eta_{j}\hat{a}^{\dagger}_{j}(t)$ to simplify
the equation, and in the interaction picture the bath operators
$\hat{B}_{jn}=\int g^{\ast}_{j}(q)\hat{c}_{n}(q)\,dq$ have the form

$\hat{B}_{jn}(u)=\int
g^{\ast}_{j}(q)\hat{c}_{n}(q)\mathrm{e}^{-i\epsilon_{n}(q)u}\,dq.$ (22)

It is also convenient to define the one-sided Fourier transforms of the bath
correlation functions,

$\displaystyle C^{(n,1)}_{jk}(\omega)$
$\displaystyle=\int_{0}^{\infty}\tr_{E}\left[\hat{B}^{\dagger}_{jn}(u)\hat{B}_{kn}\hat{\rho}_{E}\right]\mathrm{e}^{-i\omega
u}du,$ (23) $\displaystyle C^{(n,2)}_{jk}(\omega)$
$\displaystyle=\int_{0}^{\infty}\tr_{E}\left[\hat{B}_{jn}(u)\hat{B}^{\dagger}_{kn}\hat{\rho}_{E}\right]\mathrm{e}^{i\omega
u}du,$ (24)

and split them into their Hermitian and anti-Hermitian parts

$C^{(n,l)}_{jk}(\omega)=\frac{1}{2}\gamma^{(n,l)}_{jk}(\omega)+is^{(n,l)}_{jk}(\omega).$
(25)

The following expressions will also be useful when we consider separability of
the steady state,

$\displaystyle\tr_{E}\left[\hat{B}^{\dagger}_{jn}(u)\hat{B}_{kn}\hat{\rho}_{E}\right]$
$\displaystyle=\int
g_{j}(q)g^{\ast}_{k}(q)\mathrm{e}^{i\epsilon_{n}(q)u}p_{n}(\epsilon_{n}(q))\,dq,$
(26)
$\displaystyle\tr_{E}\left[\hat{B}_{jn}(u)\hat{B}^{\dagger}_{kn}\hat{\rho}_{E}\right]$
$\displaystyle=\int g^{\ast}_{j}(q)g_{k}(q)\mathrm{e}^{-i\epsilon_{n}(q)u}$
$\displaystyle\qquad\qquad\times[1-\xi_{n}p_{n}(\epsilon_{n}(q))]\,dq.$ (27)

In order to make further progress with the evolution equation (21), we need to
deal with the $u$ dependence appearing in the interaction-picture creation and
annihilation operators. We consider two situations separately and derive a
master equation for each case. First we take the global approach, in which the
baths interact with the eigenmodes of the total system Hamiltonian. This is
appropriate when the interaction between subsystems is strong relative to the
system-bath interactions and results in a so-called global master equation
Hofer  _et al._ (2017). Second, we make a local approximation, where each bath
interacts locally with the eigenmodes of the subsystem to which it couples.
This is appropriate when the interaction between subsystems is weaker than the
system-bath couplings (or comparable, see Hofer  _et al._ (2017)), and results
in a local master equation. Note that the validity regimes of the local and
global master equations have some overlap Hofer  _et al._ (2017), and hence
treating both enables us to apply our separability result in the entire range
from very weak to strong inter-system interactions. In both cases we apply the
full secular approximation to obtain a physical (completely positive trace-
preserving) master equation. The condition for this approximation to be valid
differs between the local and global approaches, and therefore the details are
provided separately below.

###  IV.1 Global master equation

Let us first derive a global master equation. To this end, we follow
D’Abbruzzo and Rossini (2021) and make use of the Bogoliubov transformation.
The system Hamiltonian (2) can be written in the form
$\hat{H}_{S}=\mathbf{\hat{a}^{\dagger}}\mathcal{H}\mathbf{\hat{a}}$, where
$\mathbf{\hat{a}^{\dagger}}=(\hat{a}^{\dagger}_{1},\ldots,\hat{a}^{\dagger}_{d})$
and $\mathcal{H}$ is a Hermitian matrix that can be diagonalised with a
unitary matrix $U$. The transformation $\mathbf{\hat{a}}=U\mathbf{\hat{b}}$
then brings the system Hamiltonian into the form

$\hat{H}_{S}=\sum_{j}\omega_{j}\hat{b}^{\dagger}_{j}\hat{b}_{j},$ (28)

where $\omega_{j}>0$ are the eigenvalues of $\mathcal{H}$, and the new
creation and annihilation operators $\hat{b}_{j},\hat{b}^{\dagger}_{j}$
satisfy the usual bosonic canonical commutation relations.

Starting from (21) and applying a full-secular approximation yields a
Markovian master equation for the system state $\hat{\rho}(t)$ (in the
Schrödinger picture) that is completely positive and trace preserving,

$\displaystyle\frac{d}{dt}\hat{\rho}(t)=-i[\hat{H}_{S}+\hat{H}_{LS},\hat{\rho}(t)]+$
$\displaystyle\sum_{\lambda=1}^{D}\sum_{u,v\in\Lambda_{\lambda}}\left[\Gamma_{uv}^{(1)}(\omega_{\lambda})\mathcal{D}(\hat{b}_{u}^{\dagger},\hat{b}_{v}^{\dagger})+\Gamma_{uv}^{(2)}(\omega_{\lambda})\mathcal{D}(\hat{b}_{u},\hat{b}_{v})\right][\hat{\rho}(t)].$
(29)

Here we have introduced the superoperator
$\mathcal{D}(\hat{X}_{i},\hat{X}_{j})$ that, given the two operators
$\hat{X}_{i},\hat{X}_{j}$, acts on a quantum state $\hat{\rho}$ as

$\mathcal{D}(\hat{X}_{i},\hat{X}_{j})[\hat{\rho}]=\hat{X}_{i}\hat{\rho}\hat{X}_{j}^{\dagger}-\frac{1}{2}\\{\hat{X}_{j}^{\dagger}\hat{X}_{i},\hat{\rho}\\}.$
(30)

The system has $D$ different energy eigenspaces labelled by $\lambda$, where
$\Lambda_{\lambda}$ is the set of normal modes associated with the $\lambda$th
eigenspace corresponding to eigenvalue $\omega_{\lambda}$, and

$\displaystyle\Gamma^{(1)}_{uv}(\omega)$
$\displaystyle=\sum_{n=1}^{N}\sum_{i,j\in
S_{n}}\gamma^{(n,1)}_{ij}(\omega)U^{\ast}_{iu}U_{jv},$ (31)
$\displaystyle\Gamma^{(2)}_{uv}(\omega)$
$\displaystyle=\sum_{n=1}^{N}\sum_{i,j\in
S_{n}}\gamma^{(n,2)}_{ij}(\omega)U_{iu}U^{\ast}_{jv},$ (32)

where the $\gamma^{(n,k)}_{ij}(\omega)$ are defined in (25). The unitary part
of the dynamics is generated by the system Hamiltonian

$\hat{H}_{S}=\sum_{\lambda=1}^{D}\sum_{u\in\Lambda_{\lambda}}\omega_{\lambda}\hat{b}^{\dagger}_{u}\hat{b}_{u}$
(33)

and the so-called Lamb-shift Hamiltonian

$\hat{H}_{LS}=\sum_{\lambda=1}^{D}\sum_{u,v\in\Lambda_{\lambda}}\varphi_{uv}(\omega_{\lambda})\hat{b}^{\dagger}_{u}\hat{b}_{v},$
(34)

where

$\displaystyle\varphi_{uv}(\omega)$ $\displaystyle=\sum_{n=1}^{N}\sum_{i,j\in
S_{n}}\left[\left(s_{ij}^{(n,1)}(\omega)+s_{ji}^{(n,2)}(\omega)\right)\right.$
(35)
$\displaystyle+\left.\left(s^{(n,1)}_{ji}(-\omega)+s^{(n,2)}_{ij}(-\omega)\right)\eta_{i}\eta_{j}\right]U^{\ast}_{iu}U_{jv},$

and the $s^{(n,k)}_{ij}(\omega)$ are also defined in (25). The Lamb-shift
Hamiltonian commutes with $\hat{H}_{S}$ and so produces a shift in the energy
levels of $\hat{H}_{S}$.

The full-secular approximation, which was required to arrive at the master
equation (IV.1), is valid provided that the differences $\nu$ of the normal-
mode eigenvalues satisfy
$\min_{\nu\neq\nu^{\prime}}|\nu-\nu^{\prime}|>1/\tau_{B}$, where $\tau_{B}$ is
the largest correlation time of any of the baths. Furthermore, this
approximation has removed all terms from the master equation that could
generate squeezing, i.e., only passive terms remain.

###  IV.2 Local master equation

Now suppose that the $N$ subsystems are weakly interacting with each other. In
this case the local approach may be justified, where the coupling between the
subsystems is neglected when calculating the effects of the environment. We
thus begin by writing the system Hamiltonian in the form
$\hat{H}_{S}=\sum_{n}\hat{H}_{n}+\hat{H}_{C}$, where $\hat{H}_{n}$ is the
local Hamiltonian for subsystem $S_{n}$ and $\hat{H}_{C}$ accounts for the
coupling between the different subsystems. Each $\hat{H}_{n}$ is excitation-
number-conserving and, as before, can be put into diagonal form

$\hat{H}_{n}=\sum_{j=1}^{M_{n}}\omega_{nj}\hat{b}_{nj}^{\dagger}\hat{b}_{nj}$
(36)

with a Bogoliubov transformation

$\hat{a}_{i}=\sum_{j=1}^{M_{n}}U^{n}_{ij}\hat{b}_{nj}.$ (37)

Here $M_{n}$ is the number of modes in subsystem $S_{n}$, $\omega_{nj}>0$ are
the eigenvalues of $\mathcal{H}_{n}$, $\hat{b}_{nj}$ are the Bogoliubov
operators associated with subsystem $S_{n}$, and $i\in S_{n}$. The components
$U_{ij}^{n}$ can be viewed as the elements of a unitary matrix $U^{n}$ after a
relabelling of the $i$ index.

These transformations do not diagonalise $\hat{H}_{S}$ and will typically
result in a complicated coupling term $\hat{H}_{C}$. However, this does not
complicate the derivation, as the coupling between the subsystems is ignored
when calculating the influence of the environment. With this in mind, starting
from (21) and applying a full-secular approximation yields a Markovian master
equation that is again completely positive and trace preserving

$\displaystyle\frac{d}{dt}\hat{\rho}(t)=-i[\hat{H}_{S}+\hat{H}_{LS},\hat{\rho}(t)]+$
$\displaystyle\sum_{n=1}^{N}\sum_{\lambda=1}^{D_{n}}\sum_{u,v\in
n_{\lambda}}\bigg{[}\Gamma_{uv}^{(n,1)}(\omega_{n\lambda})\mathcal{D}(\hat{b}_{nu}^{\dagger},\hat{b}_{nv}^{\dagger})+$
$\displaystyle\hskip
80.00012pt\Gamma_{uv}^{(n,2)}(\omega_{n\lambda})\mathcal{D}(\hat{b}_{nu},\hat{b}_{nv})\bigg{]}[\hat{\rho}(t)].$
(38)

Here $D_{n}$ is the number of energy eigenspaces associated with the subsystem
$S_{n}$, and $n_{\lambda}$ is the set of normal modes associated with the
$\lambda$th eigenspace corresponding to the eigenvalue $\omega_{n\lambda}$.
Furthermore,

$\displaystyle\Gamma^{(n,1)}_{uv}(\omega)$ $\displaystyle=\sum_{i,j\in
S_{n}}\gamma^{(n,1)}_{ij}(\omega){U^{n}_{iu}}^{\ast}U^{n}_{jv},$ (39)
$\displaystyle\Gamma^{(n,2)}_{uv}(\omega)$ $\displaystyle=\sum_{i,j\in
S_{n}}\gamma^{(n,2)}_{ij}(\omega)U^{n}_{iu}{U^{n}_{jv}}^{\ast},$ (40)

where the $\gamma^{(n,k)}_{ij}(\omega)$ are defined in (25).

While the mode transformations simplify each local subsystem Hamiltonian
$\hat{H}_{n}$, their effect on the coupling Hamiltonian $\hat{H}_{C}$ leads to
a system Hamiltonian of the form

$\hat{H}_{S}=\sum_{n=1}^{N}\sum_{\lambda=1}^{D_{n}}\sum_{u\in
n_{\lambda}}\omega_{n\lambda}\hat{b}^{\dagger}_{nu}\hat{b}_{nu}+\hat{H}_{C}^{\prime}.$
(41)

The contribution $\hat{H}_{C}^{\prime}$ contains terms of the form
$\hat{b}_{nu}^{\dagger}\hat{b}_{mv}$, which involve Bogoliubov operators for
different subsystems and energy eigenspaces. However, the full details of the
coefficients appearing in this part are not needed in what follows. The
coupling of the system to the baths produces an additional contribution to the
unitary part of the dynamics

$\hat{H}_{LS}=\sum_{n=1}^{N}\sum_{\lambda=1}^{D_{n}}\sum_{u,v\in
n_{\lambda}}\varphi^{(n)}_{uv}(\omega_{n\lambda})\hat{b}^{\dagger}_{nu}\hat{b}_{nv},$
(42)

where

$\displaystyle\varphi^{(n)}_{uv}(\omega)$ $\displaystyle=\sum_{i,j\in
S_{n}}\left[\left(s_{ij}^{(n,1)}(\omega)+s_{ji}^{(n,2)}(\omega)\right)\right.$
(43)
$\displaystyle+\left.\left(s^{(n,1)}_{ji}(-\omega)+s^{(n,2)}_{ij}(-\omega)\right)\eta_{i}\eta_{j}\right]{U^{n}_{iu}}^{\ast}U^{n}_{jv},$

and the $s^{(n,k)}_{ij}(\omega)$ are also defined in (25). In contrast to the
global approach, this term does not commute with the system Hamiltonian and
therefore does not simply shift the energies of the system Hamiltonian.

In the local approach, the validity of the full-secular approximation is
determined by the local Hamiltonians of each subsystem. For each subsystem,
the differences $\nu_{n}$ of the normal-mode eigenvalues associated with the
local Hamiltonian $\hat{H}_{n}$ must satisfy
$\min_{\nu_{n}\neq\nu_{n}^{\prime}}|\nu_{n}-\nu_{n}^{\prime}|>1/\tau_{B}$.
Again, this approximation results in a master equation containing only passive
terms that cannot generate squeezing.

###  IV.3 The steady-state solution

We now show that whenever the steady state of the global or local master
equation is unique, it is a Gaussian state with a covariance matrix of the
form (17). Thus, for any initial system state the steady state will be a
separable Gaussian state. We first provide a short proof that assumes the
steady state is unique, and then derive a condition for the steady state to be
unique. Both derivations make use of the fact that the global and local master
equations take Gaussian states to Gaussian states. This is because the
Hamiltonian part is quadratic in the Bogoliubov operators, and the Lindblad
operators in the dissipative part are linear in the Bogoliubov operators.

For the simple proof, note that the right hand side of both the global (IV.1)
and local (IV.2) master equations can be written in terms of a Liouvillian
$\mathcal{L}$ as $\mathcal{L}\hat{\rho}(t)$. In each case $\mathcal{L}$
commutes with the superoperator
$\mathcal{N}\boldsymbol{\cdot}=\left[\hat{N},\boldsymbol{\cdot}\right]$, where
the total excitation-number operator is
$\hat{N}=\sum_{j}\hat{b}_{j}^{\dagger}\hat{b}_{j}$ in the global approach and
$\hat{N}=\sum_{n,j}\hat{b}_{nj}^{\dagger}\hat{b}_{nj}$ in the local approach.
Therefore, any state that is initially block diagonal in the $\hat{N}$
eigenbasis, with no coherences between different number sectors, will remain
block diagonal for all time. Such a state has no squeezing and thus a
covariance matrix of the separable form (17). This proves our result provided
that the steady state is unique because an initial Gaussian state with
covariance matrix of the form (17) must then remain a Gaussian state with a
covariance matrix of this form for all time.

We now derive a condition for the steady state to be unique. Recall that any
initial system state $\hat{\rho}(0)$ can be written as

$\hat{\rho}(0)=\frac{1}{\pi^{2d}}\int\langle{\alpha}|\hat{\rho}(0)|{\beta}\rangle|{\alpha}\rangle\langle{\beta}|d^{2}{\alpha}d^{2}{\beta},$
(44)

where $|{\alpha}\rangle$ is a $d$-mode Glauber coherent state with amplitudes
${\alpha}=(\alpha_{1},\ldots,\alpha_{d})$. By the linearity of the master
equations, if we obtain the long-time evolution of each term
$|{\alpha}\rangle\langle{\beta}|$ in the decomposition (44) we can then deduce
the form of the steady state.

Let us first deal with the coherent states $|{\alpha}\rangle\langle{\alpha}|$
appearing in (44). These are Gaussian states, each with the same covariance
matrix $I_{2d}/2$ but different centres. We now show that in the long-time
limit every coherent state $|{\alpha}\rangle\langle{\alpha}|$ tends to the
same Gaussian state $\hat{\rho}_{G}$ with zero first moments and a covariance
matrix of the form (17). Because the master equations map Gaussian states to
Gaussian states, the initial coherent states $|{\alpha}\rangle$ will remain
Gaussian for all time. It is then enough to work out the steady-state
expectation values of the first and second moments.

We focus on the global master equation (IV.1), as the derivation for the local
master equation proceeds along the same lines. From (IV.1) it follows that the
time evolution of an operator $\hat{O}$ is determined by

$\frac{d}{dt}\hat{O}=i[\hat{H}_{S}+\hat{H}_{LS},\hat{O}]+\sum_{\lambda=1}^{D}\sum_{u,v\in\Lambda_{\lambda}}\Gamma_{uv}^{(1)}(\omega_{\lambda})\left(\hat{b}_{v}\hat{O}\hat{b}^{\dagger}_{u}-\frac{1}{2}\\{\hat{b}_{v}\hat{b}^{\dagger}_{u},\hat{O}\\}\right)+\Gamma_{uv}^{(2)}(\omega_{\lambda})\left(\hat{b}^{\dagger}_{v}\hat{O}\hat{b}_{u}-\frac{1}{2}\\{\hat{b}^{\dagger}_{v}\hat{b}_{u},\hat{O}\\}\right).$
(45)

The equations of motion for the one and two-point correlators
$\langle\hat{b}_{i}\rangle$ and $\langle\hat{b}_{i}\hat{b}_{j}\rangle$ are
thus obtained by setting $\hat{O}$ to $\hat{b}_{i}$ and
$\hat{b}_{i}\hat{b}_{j}$ respectively. Now recall that every coherent state
$|{\alpha}\rangle$ has the same covariance matrix, $I_{2d}/2$, and that the
Bogoliubov operators are related to the original operators by the unitary
transformation $\mathbf{\hat{a}}=U\mathbf{\hat{b}}$. From this we see that
$\langle\hat{b}_{i}\hat{b}_{j}\rangle-\langle\hat{b}_{i}\rangle\langle\hat{b}_{j}\rangle=0$
for all $i,j$ at time $t=0$. By plugging these initial conditions into the
equations of motion for
$\beta_{ij}=\langle\hat{b}_{i}\hat{b}_{j}\rangle-\langle\hat{b}_{i}\rangle\langle\hat{b}_{j}\rangle$
it is straightforward to check that $\beta_{ij}=0$ for all time (see Appendix
A). Therefore,
$\langle\hat{a}_{i}\hat{a}_{j}\rangle-\langle\hat{a}_{i}\rangle\langle\hat{a}_{j}\rangle=0$
at every time for all $i,j$. Because the equations of motion for the
(co)variances do not couple to the first moments, the covariance matrix of
every initial coherent state $|{\alpha}\rangle$ tends to the same steady-state
covariance matrix of separable form (17).

All that remains is to show that the centre of every initial coherent state
goes to the same value in the long-time limit. For this we need the equations
of motion for the first moments

$\displaystyle\frac{d}{dt}\langle\hat{b}_{j}\rangle$
$\displaystyle=-i\omega_{j}\langle\hat{b}_{j}\rangle-i\sum_{u\in\Lambda_{j}}\varphi_{ju}(\omega_{j})\langle\hat{b}_{u}\rangle$
$\displaystyle\qquad+\frac{1}{2}\sum_{u\in\Lambda_{j}}\left(\Gamma^{(1)}_{ju}(\omega_{j})-\Gamma^{(2)}_{uj}(\omega_{j})\right)\langle\hat{b}_{u}\rangle.$
(46)

For simplicity we have redefined $\omega_{j}$ as the eigenvalue associated
with the mode $\hat{b}_{j}$ and $\Lambda_{j}$ as the eigenspace to which the
eigenvalue $\omega_{j}$ belongs. This set of linear equations can be written
in the form $d{b}/dt=R{b}$, where $b_{i}=\langle\hat{b}_{i}\rangle$ and $R$ is
a non-Hermitian matrix. The Hermitian part of $R$ is negative semidefinite
(see Appendix B), and therefore, if $R$ is invertible, then ${b}={0}$ is the
only steady-state solution. If this condition is satisfied, it follows that
every initial state $|{\alpha}\rangle$ tends to the same Gaussian state
$\hat{\rho}_{G}$ with zero first moments and a covariance matrix of the form
(17).

The long-time dynamics of the cross terms $|{\alpha}\rangle\langle{\beta}|$
(${\alpha}\neq{\beta}$) appearing in (44) require a little more work. The
Wigner-Weyl transformation of $|{\alpha}\rangle\langle{\beta}|$ is a complex
Gaussian in phase space that is highly oscillatory. A familiar example where
terms like this appear is in the Wigner function of a cat state. For Lindblad
master equations with linear Lindblad operators and quadratic Hamiltonians,
the coupling to the environment exponentially damps the oscillatory parts of
these phase-space functions. This is a manifestation of decoherence.
Therefore, the complex oscillatory parts cannot persist in the steady state.
The dynamics of complex Gaussian wave packets generated by equations of
Lindblad type were analysed in Graefe  _et al._ (2018). Using these results,
we find that each $|{\alpha}\rangle\langle{\beta}|$ tends to
$N_{\alpha\beta}\hat{\rho}_{G}$. Here $\hat{\rho}_{G}$ is the same Gaussian
steady state we found before, and the $N_{\alpha\beta}$ are normalisation
factors that guarantee the steady state $\hat{\rho}_{\infty}$ is normalised.
We can then conclude that the steady state $\hat{\rho}_{\infty}$ is the
separable Gaussian state $\hat{\rho}_{G}$.

The same procedure can be followed with the local master equation (IV.2) to
obtain the same result. The only difference, barring additional indices, is
the coupling term $\hat{H}_{C}^{\prime}$ appearing in the system Hamiltonian
(41). However, this does not pose any complications when following the steps
above, as $\hat{H}_{C}^{\prime}$ only includes terms of the form
$\hat{b}_{nu}^{\dagger}\hat{b}_{mv}$ and cannot introduce any squeezing.

##  V Conclusion

We have shown that autonomous, steady-state entanglement generation is
impossible in bosonic systems of arbitrary size undergoing quadratic
excitation-number-preserving interactions (i.e., passive Gaussian interactions
or linear interferometers) and weakly coupled to thermal bosonic and/or spin
baths at different temperatures. Our result applies to both the local and
global regimes, which cover the entire range from weak to strong inter-system
coupling. Our result holds provided that: the system Hamiltonian is stable
($\hat{H}_{S}>0$); a secular approximation is justified; and the steady state
is unique. This result contrasts with previous findings for finite-dimensional
autonomous thermal machines with excitation-preserving interactions, which can
generate steady-state entanglement strong enough to exhibit both steering and
nonlocality. Passive Gaussian interactions, even combined with non-bosonic
baths, are sufficiently restrictive to preclude any steady-state entanglement.
We note that the condition for the steady state to be unique ($R$ invertible)
will most likely be satisfied in all practical situations, where the model
parameter values will never be completely symmetric.

Having ruled out steady-state entanglement in this case, it is interesting to
ask whether bosonic steady-state entanglement can be generated in more
complicated settings. One natural next step would be to include higher-order
inter-system interactions, such as three- or four-wave mixing. The challenge
here may be to determine whether the (generally non-Gaussian) steady-state is
separable or not. One could also attempt to go beyond the weak-coupling limit
by applying, for example, the reaction coordinate method Nazir and Schaller
(2018).

###### Acknowledgements.

We thank Patrick Potts for enlightening discussions and feedback. This work
was supported by a research grant (40864) from VILLUM FONDEN, as well as by
the Carlsberg Foundation CF19-0313 and the Independent Research Fund Denmark
7027-00044B.

## References

  * Horodecki  _et al._ (2009) R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Rev. Mod. Phys. 81, 865 (2009).
  * Gisin and Thew (2007) N. Gisin and R. Thew, Nature Photonics 1, 165 (2007).
  * Giovannetti  _et al._ (2011) V. Giovannetti, S. Lloyd, and L. Maccone, Nat Photon 5, 222 (2011).
  * Pironio  _et al._ (2016) S. Pironio, V. Scarani, and T. Vidick, New Journal of Physics 18, 100202 (2016).
  * Pirandola  _et al._ (2020) S. Pirandola, U. L. Andersen, L. Banchi, M. Berta, D. Bunandar, R. Colbeck, D. Englund, T. Gehring, C. Lupo, C. Ottaviani, J. L. Pereira, M. Razavi, J. S. Shaari, M. Tomamichel, V. C. Usenko, G. Vallone, P. Villoresi, and P. Wallden, Adv. Opt. Photon. 12, 1012 (2020).
  * Plenio  _et al._ (1999) M. B. Plenio, S. F. Huelga, A. Beige, and P. L. Knight, Phys. Rev. A 59, 2468 (1999).
  * Kim  _et al._ (2002) M. S. Kim, J. Lee, D. Ahn, and P. L. Knight, Phys. Rev. A 65, 040101(R) (2002).
  * Jakóbczyk (2002) L. Jakóbczyk, Journal of Physics A: Mathematical and General 35, 6383 (2002).
  * Braun (2002) D. Braun, Phys. Rev. Lett. 89, 277901 (2002).
  * Benatti  _et al._ (2003) F. Benatti, R. Floreanini, and M. Piani, Phys. Rev. Lett. 91, 070402 (2003).
  * Burgarth and Giovannetti (2007) D. Burgarth and V. Giovannetti, Phys. Rev. A 76, 062307 (2007).
  * Bellomo  _et al._ (2008) B. Bellomo, R. Lo Franco, S. Maniscalco, and G. Compagno, Phys. Rev. A 78, 060302(R) (2008).
  * Naseem and Müstecaplıoğlu (2022) M.-T. Naseem and Ö. Müstecaplıoğlu, Quantum Sci. Technol. 7, 045012 (2022).
  * Diehl  _et al._ (2008) S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Büchler, and P. Zoller, Nature Physics 4, 878 (2008).
  * Verstraete  _et al._ (2009) F. Verstraete, M. M. Wolf, and J. Ignacio Cirac, Nature Physics 5, 633 (2009).
  * Kraus  _et al._ (2008) B. Kraus, H. P. Büchler, S. Diehl, A. Kantian, A. Micheli, and P. Zoller, Phys. Rev. A 78, 042307 (2008).
  * Ticozzi and Viola (2014) F. Ticozzi and L. Viola, Quantum Information and Computation 14, 0265 (2014).
  * Plenio and Huelga (2002) M. B. Plenio and S. F. Huelga, Phys. Rev. Lett. 88, 197901 (2002).
  * Hartmann  _et al._ (2006) L. Hartmann, W. Dür, and H.-J. Briegel, Phys. Rev. A 74, 052304 (2006).
  * Hartmann  _et al._ (2007) L. Hartmann, W. Dür, and H.-J. Briegel, New Journal of Physics 9, 230 (2007).
  * Žnidarič (2012) M. Žnidarič, Phys. Rev. A 85, 012324 (2012).
  * Quiroga  _et al._ (2007) L. Quiroga, F. J. Rodríguez, M. E. Ramírez, and R. París, Phys. Rev. A 75, 032308 (2007).
  * Bellomo and Antezza (2013a) B. Bellomo and M. Antezza, EPL (Europhysics Letters) 104, 10006 (2013a).
  * Bellomo and Antezza (2013b) B. Bellomo and M. Antezza, New Journal of Physics 15, 113052 (2013b).
  * Brask  _et al._ (2015) J. B. Brask, G. Haack, N. Brunner, and M. Huber, New Journal of Physics 17, 113029 (2015).
  * Tavakoli  _et al._ (2018) A. Tavakoli, G. Haack, M. Huber, N. Brunner, and J. B. Brask, Quantum 2, 73 (2018).
  * Man  _et al._ (2019) Z.-X. Man, A. Tavakoli, J. B. Brask, and Y.-J. Xia, Physica Scripta 94, 075101 (2019).
  * Tavakoli  _et al._ (2020) A. Tavakoli, G. Haack, N. Brunner, and J. B. Brask, Phys. Rev. A 101, 012315 (2020).
  * Brask  _et al._ (2022) J. B. Brask, F. Clivaz, G. Haack, and A. Tavakoli, Quantum 6, 672 (2022).
  * 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, Rev. Mod. Phys. 84, 621 (2012).
  * Eisert  _et al._ (2002) J. Eisert, S. Scheel, and M. B. Plenio, Phys. Rev. Lett. 89, 137903 (2002).
  * Giedke and Ignacio Cirac (2002) G. Giedke and J. Ignacio Cirac, Phys. Rev. A 66, 032316 (2002).
  * Fiurášek (2002) J. Fiurášek, Phys. Rev. Lett. 89, 137904 (2002).
  * Niset  _et al._ (2009) J. Niset, J. Fiurášek, and N. J. Cerf, Phys. Rev. Lett. 102, 120501 (2009).
  * Lloyd and Braunstein (1999) S. Lloyd and S. L. Braunstein, Phys. Rev. Lett. 82, 1784 (1999).
  * Knill  _et al._ (2001) E. Knill, R. Laflamme, and G. J. Milburn, Nature 409, 46 (2001).
  * Werner and Wolf (2001) R. F. Werner and M. M. Wolf, Phys. Rev. Lett. 86, 3658 (2001).
  * Wolf  _et al._ (2003) M. M. Wolf, J. Eisert, and M. B. Plenio, Phys. Rev. Lett. 90, 047904 (2003).
  * Rivas and Huelga (2012) Á. Rivas and S. F. Huelga,  _Open Quantum Systems_ (Springer Berlin, Heidelberg, 2012).
  * Potts (2019) P. Potts, ArXiv e-prints (2019), 10.48550/arXiv.1906.07439, arXiv:1906.07439 [quant-ph] .
  * Hofer  _et al._ (2017) P. P. Hofer, M. Perarnau-Llobet, L. D. M. Miranda, G. Haack, R. Silva, J. B. Brask, and N. Brunner, New Journal of Physics 19, 123037 (2017).
  * D’Abbruzzo and Rossini (2021) A. D’Abbruzzo and D. Rossini, Phys. Rev. A 103, 052209 (2021).
  * Graefe  _et al._ (2018) E. M. Graefe, B. Longstaff, T. Plastow, and R. Schubert, J. Phys. A: Math. Theor. 51, 365203 (2018).
  * Nazir and Schaller (2018) A. Nazir and G. Schaller, “The reaction coordinate mapping in quantum thermodynamics,” in  _Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions_, edited by F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso (Springer International Publishing, Cham, 2018) pp. 551–577.

##  Appendix A Time evolution of the covariance matrix elements $\beta_{jk}$

The equations of motion for the covariance matrix elements
$\beta_{jk}=\langle\hat{b}_{j}\hat{b}_{k}\rangle-\langle\hat{b}_{j}\rangle\langle\hat{b}_{k}\rangle$
can be obtained from the operator evolution equation (45). We first compute
the equations for $\langle\hat{b}_{j}\rangle$ and
$\langle\hat{b}_{j}\hat{b}_{k}\rangle$, and then combine the results to find

$\displaystyle\frac{d}{dt}\beta_{jk}=$
$\displaystyle-i\left(\omega_{j}+\omega_{k}\right)\beta_{jk}+$
$\displaystyle\sum_{u\in\Lambda_{j}}f_{ju}(\omega_{j})\beta_{ku}+\sum_{u\in\Lambda_{k}}f_{ku}(\omega_{k})\beta_{ju}.$
(47)

For simplicity we have redefined $\omega_{j}$ as the eigenvalue associated
with mode $\hat{b}_{j}$ and $\Lambda_{j}$ as the eigenspace to which
$\omega_{j}$ belongs. The functions $f_{iu}(\omega)$ are defined as

$f_{iu}(\omega)=\frac{1}{2}\left[\Gamma_{iu}^{(1)}(\omega)-\Gamma_{ui}^{(2)}(\omega)\right]-i\varphi_{iu}(\omega),$
(48)

in terms of functions appearing in the main text. The key result of this
Appendix is that when each element $\beta_{jk}=0$ at $t=0$, then each
$\beta_{jk}$ is zero for all time. This follows immediately from (A).

##  Appendix B The Hermitian part of $R$ is negative semidefinite

Let the Hermitian part of $R$ be denoted by $R_{H}=(R+R^{\dagger})/2$. From
the equations of motion for the first moments in the main text (IV.3) we
obtain the matrix elements

$[R_{H}]_{ij}=\frac{1}{2}\left(\Gamma^{(1)}_{ij}(\omega_{i})-\Gamma^{(2)}_{ji}(\omega_{i})\right)\delta_{\omega_{i},\omega_{j}},$
(49)

where $\omega_{i}$ is the eigenvalue associated with the Bogoliubov mode
$\hat{b}_{i}$. We permute the rows and columns of $R_{H}$ to bring it into
block diagonal form, where each block corresponds to the same eigenvalue.
After a relabelling of indices, the block $B$ corresponding to the eigenvalue
$\omega$ has the elements

$B_{ij}=\frac{1}{2}\left(\Gamma^{(1)}_{ij}(\omega)-\Gamma^{(2)}_{ji}(\omega)\right).$
(50)

We now show that this block is negative semidefinite.

Making use of the expressions (23)-(27) and (31)-(32) we find that the matrix
elements in (50) can be written as

$\displaystyle B_{ij}$
$\displaystyle=\pi\sum_{n=1}^{N}\int\tilde{g}_{in}(q)\tilde{g}^{\ast}_{jn}(q)\delta(\epsilon_{n}(q)-\omega)$
(51)
$\displaystyle\qquad\qquad\quad\times\left((1+\xi_{n})p_{n}(\epsilon_{n}(q))-1\right)\,dq,$

where we have defined

$\tilde{g}_{in}(q)=\sum_{j\in S_{n}}U^{\ast}_{ji}g_{j}(q)$ (52)

in terms of the unitary matrix $U$ that diagonalises the system Hamiltonian
and the coupling parameters between the system and environment $g_{j}(q)$. For
every complex vector $x$

$\displaystyle x^{\dagger}Bx$
$\displaystyle=\pi\sum_{n=1}^{N}\int|x_{n}(q)|^{2}\delta(\epsilon_{n}(q)-\omega)$
(53)
$\displaystyle\qquad\qquad\quad\times\left((1+\xi_{n})p_{n}(\epsilon_{n}(q))-1\right)\,dq,$

with $x_{n}(q)=\sum_{j}x_{j}\tilde{g}^{\ast}_{jn}(q)$. For the bosonic baths,
$\xi_{n}=-1$ and the contribution to (53) is less than or equal to zero. For
the spin baths, $\xi_{n}=+1$ but $p_{n}(\omega)<1/2$ since $\omega>0$. The
contribution from the spin baths is therefore also less than or equal to zero,
and it follows that $x^{\dagger}Bx\leq 0$. This holds for every block of
$R_{H}$, and so the Hermitian part of $R$ is negative semidefinite.