# Thermal equilibrium in Gaussian dynamical semigroups

Fabricio Toscano<EMAIL_ADDRESS>Instituto de Física, Universidade
Federal do Rio de Janeiro, 21941-972, Rio de Janeiro, Brazil  Fernando Nicacio
<EMAIL_ADDRESS>Instituto de Física, Universidade Federal do Rio de
Janeiro, 21941-972, Rio de Janeiro, Brazil  Universität Wien, NuHAG, Fakultät
für Mathematik, A-1090 Wien, Austria.

###### Abstract

We characterize all Gaussian dynamical semigroups in continuous variables
quantum systems of $n$-bosonic modes which have a thermal Gibbs state as a
stationary solution. This is performed through an explicit relation between
the diffusion and dissipation matrices, which characterize the semigroup
dynamics, and the covariance matrix of the thermal equilibrium state. We also
show that Alicki’s quantum detailed-balance condition, based on a Gelfand-
Naimark-Segal inner product, allows the determination of the temperature
dependence of the diffusion and dissipation matrices and the identification of
different Gaussian dynamical semigroups which share the same thermal
equilibrium state.

##  I Introduction

In modern Quantum Information Theory for continuous variable systems, i.e,
systems described by $n$-bosonic modes, Gaussian channels are the standard
models in most of the quantum communication protocols Holevo et al. (1999);
Holevo and Werner (2001); Holevo (2002); Cerf et al. (2007); Caruso et al.
(2008); Holevo (2019). These channels are defined as those bosonic channels
which transform Gaussian states into Gaussian states Weedbrook et al. (2012);
Adesso et al. (2014). Further, these states have an exceptional role in
quantum communication as, for example, they are optimal for the transmission
of classical information through Gaussian bosonic quantum channels with
additive capacity Wolf et al. (2006).

The most general form of a one-parameter Gaussian channel of $n$-bosonic modes
is a Gaussian Dynamical Semigroup (GDS) Heinosaari et al. (2010); Toscano et
al. (2021), constituting thus the tool to describe the dynamics of all
memoryless continuous-in-time Gaussian quantum channels. This is the reason
why they are widely used to describe noisy quantum channels in continuous
variable systems Heinosaari et al. (2010); Giovannetti et al. (2010); De Palma
et al. (2016). Further,  GDSs are able to describe all processes which can
formally be written as decomposition and production of noninteracting
particles (or quasi-particles) which can be treated at least approximately as
bosons Alicki and Lendi (2007). In this context,  GDSs are known as quasi-free
completely-positive semigroups Vanheuverzwijn (1978, 1979); Demoen et al.
(1979), which happens, for example, in damped collective modes in deep
inelastic collisions Sǎndulescu and Scutaru (1987).

The dynamics of any Quantum Dynamical Semigroup (QDS), not necessarily
Gaussian, also known as Quantum Markov Semigroups Carlen and Maas (2017), is
described by a master equation in the Lindblad form Lindblad (1976); Gorini et
al. (1976). In an analogous way, a  GDS verifies a Lindblad master equation
where the unitary part of the evolution is set by a quadratic Hamiltonian
which describes the $n$-bosonic modes. By a quadratic Hamiltonian, we mean one
composed by products of any two canonical conjugate operators, positions and
momenta, or the set of self-adjoint operators that constitutes a
representation of the Heisenberg canonical commutation relations of the
particular bosonic system considered Tarasov (2008). Meanwhile, the non-
unitary part is given by Lindblad operators corresponding to complex linear
functions of positions and momenta, which simplifies the non-unitary dynamics
to be described only in terms of two real $2n\times 2n$ matrices, the
diffusion and dissipation matrices Nicacio et al. (2016); Toscano et al.
(2021).

Alternatively, the Weyl-Wigner representation Ozorio de Almeida (1998) for the
master equation of a  GDSs can be employed and corresponds to a linear Fokker-
Planck equation for the evolution of the Wigner function of the evolved states
Nicacio et al. (2010, 2016); Toscano et al. (2021). This is the quantum
counterpart of the classical channels corresponding to Ornstein-Uhlenbeck
processes, thus also known as Bose-Ornstein–Uhlenbeck semigroups Carlen and
Maas (2017), whose evolved probability distributions satisfy exactly the same
Fokker-Planck equation as the  GDSs.

Notably when the  GDS dynamics has a steady state, this will be a dynamically
invariant Gaussian state which attracts, over long times, the evolution of any
initial condition Frigerio (1977, 1978); Carmichael (1999). Thus, for a given
quadratic Hamiltonian, the characterization of stationary situations in  GDSs
corresponds to find all the diffusion and dissipation matrices that allow a
stationary state Toscano et al. (2021). Of particular importance are the
stationary states which also corresponds to thermal equilibrium states,
characterized by a Gibbs state of a quadratic Hamiltonian. This is the main
subject of our study.

Thermal equilibrium in  QDSs in finite-dimensional quantum systems has its own
well-established theory based on the so called Quantum Detailed Balance
Condition (QDBC), as first stated by Alicki Alicki (1976); Alicki and Lendi
(2007), see also Carlen and Maas (2017)). Of note are the results in Carlen
and Maas (2017) where the authors prove, using the  QDBC, that the evolution
governed by a  QDS is a gradient flow, in a particular Riemannian metric on
the set of states, for the relative entropy of a state with respect to a Gibbs
state. On the other hand, the study of thermal equilibrium in  GDSs is scarce,
which is particularly true for the multimode scenario; a prominent exception
is the extension of the results for  QDSs in finite-dimensional quantum
systems to the case of a one mode  GDS performed in Carlen and Maas (2017).

Here we fill this gap and give a complete characterization of all $n$-mode
GDSs with a thermal equilibrium state. To this aim, we employ the Fokker-
Planck equation for the evolution of the Wigner function and show that the
thermal probability current is always null in every phase space point. This
enable us to conclude that  GDSs with thermal equilibrium are characterized by
a set of three commuting matrices. Two of the them are the Hamiltonian
matrices associated to the covariance matrix of the thermal state and the
diffusion matrix of the  GDS. The third one is the skew-Hamiltonian matrix
associated to the dissipation matrix of the  GDS. By another side, this
condition neglects that different  GDSs, characterized by the diffusion and
dissipation matrices, may share the same thermal equilibrium state, as a
consequence of the fact that the relation among that matrices does not set
their temperature dependence. To circumvent this, we show that these
characterizations are possible by extending Alicki’s  QDBC to bosonic-mode
systems.

The extended  QDBC leads to a master equation for a  GDSs with the form of a
Quantum Master Optical Equation (QOME) Gardiner and Zoller (2000) and the
temperature dependence of the diffusion and dissipation matrices is
established for this type of  GDSs. From this temperature characterization, we
establish that all  GDS that leads to thermal equilibrium satisfy a  QDBC if
we allow an arbitrary temperature dependence for coupling constants between
the system and the environment. Finally, we discriminate the Hamiltonians
corresponding to the unitary part of the dynamics of a  GDS which allow the
occurrence of the thermalization. Although, we show that the thermalization
process itself is not affected by these Hamiltonians.

The paper is organized as follows. In Sec.II we introduce the  GDSs, their
action on Gaussian states (Sec.II.1) and the Weyl-Wigner formalism to describe
its dynamics (Sec.II.2). In the introduction of Sec.III, we establish the
general time dependence of the first and second order moments in  GDSs with
stationary solutions. Then, in Sec.III.1 we set up the problem of having
thermal equilibrium as stationary solutions and describe general properties
that must be satisfied by  GDSs with thermal equilibrium. Section III.2
contains one of our main results: a theorem that characterizes all the  GDSs
with a thermal equilibrium state. The extension of the  QDBC to $n$-bosonic
mode systems is placed in Sec. IV, where five theorems are presented. These
theorems completely characterize all  GDSs satisfying the detailed balance. In
Section V, we show that the master equation of  GDSs satisfying the  QDBC
always corresponds to a  QOME; this section finishes with a discussion about
entanglement properties of its thermal equilibrium state solution. We further
explore the characterization of thermal equilibrium states in  GDSs that
satisfy a  QDBC in Sec. VI, where the temperature dependence of the diffusion
and dissipation matrices is developed in Sec. VI.1, the high and low
temperature limits are described in Sec. VI.2, and in Sec. VI.3 we explain the
pure diffusive regime, where the stationary solution is lost. In Sec. VII we
describe the necessary structure of a quadratic Hamiltonian, governing the
unitary part of the  GDS, that has a thermal equilibrium state. In this
section we also clarify the role of this Hamiltonian in the process of
thermalization. Finally, we summarize our findings in Sec.VIII. Some auxiliary
calculations and technical proofs are presented in the Appendixes A, B, C, D,
E, and H.

##  II Gaussian Dynamical Semigroups

In the Schrödinger picture, a  QDS is ruled by the Lindblad master equation
(LME) Lindblad (1976); Gorini et al. (1976)

$\derivative{\hat{\rho}_{t}}{t}=\mathcal{L}[\hat{\rho}_{t}]=\mathcal{L}_{\mathrm{U}}[\hat{\rho}_{t}]+\mathcal{L}_{\mathrm{NU}}[\hat{\rho}_{t}],$
(1)

where

$\displaystyle\mathcal{L}_{\mathrm{U}}[\,\cdot\,]=-\frac{\imath}{\hbar}[\hat{H}_{\mathrm{eff}},\cdot\,]$
(2a)
$\displaystyle\mathcal{L}_{\mathrm{NU}}[\,\cdot\,]=\frac{1}{2\hbar}\sum_{k=1}^{K}\left(2\hat{L}_{k}\cdot\hat{L}_{k}^{\dagger}-\hat{L}_{k}^{\dagger}\hat{L}_{k}\cdot-\,\cdot\hat{L}_{k}^{\dagger}\hat{L}_{k}\right)$
(2b)

are, respectively, the infinitesimal generators of the unitary and non-unitary
parts of the evolution. The operator $\hat{H}_{\mathrm{eff}}$ is the effective
Hamiltonian of the system and $\hat{L}_{k}(k=1,...,K)$ are the Lindblad
operators.

With the help of the Hilbert-Schmidt inner product,
$\langle{\hat{A},\hat{B}}\rangle={\Tr}(\hat{A}^{\dagger}\hat{B})$, the adjoint
$\bar{\mathcal{L}}$ of the superoperator $\mathcal{L}$ is defined by

$\langle{\mathcal{L}[\hat{A}],\hat{B}}\rangle=\langle{\hat{A},\bar{\mathcal{L}}[\hat{B}]}\rangle,$
(3)

which enables us to write the Heisenberg picture of Eq.(1) for an observable
$\hat{O}_{t}$ Tarasov (2008):

$\derivative{\hat{O}_{t}}{t}=\bar{\mathcal{L}}[\hat{O}_{t}]=\bar{\mathcal{L}}_{\mathrm{U}}[\hat{O}_{t}]+\bar{\mathcal{L}}_{\mathrm{NU}}[\hat{O}_{t}],$
(4)

where

$\displaystyle\bar{\mathcal{L}}_{\mathrm{U}}[\,\cdot\,]$ $\displaystyle=$
$\displaystyle\frac{\imath}{\hbar}[\hat{H}_{\mathrm{eff}},\cdot\,],$ (5a)
$\displaystyle\bar{\mathcal{L}}_{\mathrm{NU}}[\cdot]$ $\displaystyle=$
$\displaystyle\frac{1}{2\hbar}\sum_{k=1}^{K}\left(2\hat{L}_{k}^{\dagger}\cdot\hat{L}_{k}-\,\cdot\hat{L}_{k}^{\dagger}\hat{L}_{k}-\hat{L}_{k}^{\dagger}\hat{L}_{k}\cdot\,\right).$
(5b)

Since $\mathcal{L}$ is time-independent, the solution of (1) is formally given
by $\hat{\rho}_{t}=e^{t\mathcal{L}}\hat{\rho}_{0}$, which is the evolution of
an initial condition $\hat{\rho}_{0}$, and the set
$\\{\Lambda_{t}=e^{t\mathcal{L}}\\}_{t\geq 0}$ is properly the  QDS in the
Schrödinger picture Breuer et al. (2002). The solution of (4) is formally
written as $\hat{O}_{t}=e^{t\bar{\mathcal{L}}}\hat{O}_{0}$, which gives the
evolution of an initial condition $\hat{O}_{0}$, where the set
$\\{\bar{\Lambda}_{t}=e^{t\bar{\mathcal{L}}}\\}_{t\geq 0}$ is the Heisenberg
picture version of the  QDS.

The kinematics of a system with $n$-bosonic modes is described by a
$2n$-dimensional column vector of canonical operators,

$\hat{\mathbf{x}}=(\hat{q}_{1},\ldots,\hat{q}_{n},\hat{p}_{1},\ldots,\hat{p}_{n})^{\intercal},$
(6)

satisfying the canonical commutation relation
$[\hat{x}_{j},\hat{x}_{k}]=\imath\hbar\mathsf{J}_{jk}\hat{1}$, where

$\mathsf{J}=\begin{pmatrix}0&\mathbb{1}\\\
-\mathbb{1}&0\end{pmatrix},\,\,\,\mathsf{J}^{-1}=-\mathsf{J}=\mathsf{J}^{\intercal},$
(7)

is a $2n\times 2n$ real antisymmetric symplectic matrix and $\mathbb{1}$ is
the $n\times n$ identity matrix.

In a $n$-bosonic-mode system a Gaussian dynamical semigroup is a  QDS with the
Hamiltonian and Lindblad operators given by

$\displaystyle\hat{H}_{\mathrm{eff}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\hat{\mathbf{x}}^{\intercal}\mathbb{B}^{\prime}\hat{\mathbf{x}}+\hat{\mathbf{x}}^{\intercal}\mathsf{J}\boldsymbol{\xi}^{\prime},$
(8a) $\displaystyle\hat{L}_{k}$ $\displaystyle=$
$\displaystyle\mathbf{l}_{k}^{\intercal}\mathsf{J}\hat{\mathbf{x}},\,\,\,(k=1,...,K),$
(8b)

where $\mathbb{B}^{\prime}=(\mathbb{B}^{\prime})^{\intercal}$ is the Hessian
matrix of the Hamiltonian, $\boldsymbol{\xi}^{\prime}$ is an
$(2n)$-dimensional real column vector, and the $\mathbf{l}_{k}$’s are
$(2n)$-dimensional complex column vectors. In this case, according to Toscano
et al. (2021), the superoperator in (2b) becomes111From now on, we will use
$\tr$ to denote the trace of a matrix and $\Tr$ to denote the trace of an
operator.

$\displaystyle\mathcal{L}_{\mathrm{NU}}^{\mathrm{G}}[\hat{\rho}_{t}]$
$\displaystyle=$
$\displaystyle\tr\\!\left(\mathbb{C}\mathsf{J}\partialderivative{\hat{\mathbf{x}}}\,[\hat{\mathbf{x}}^{\intercal}\hat{\rho}_{t}]+\frac{1}{2}\mathbb{D}\partialderivative{\hat{\mathbf{x}}}\partialderivative{\hat{\mathbf{x}}^{\intercal}}\hat{\rho}_{t}\right),$
(9)

where
$\partialderivative{\hat{\mathbf{x}}}\cdot=\frac{\imath}{\hbar}[(\mathsf{J}\hat{\mathbf{x}}),\cdot]$
is a column vector operator whose components are
$\partialderivative{x_{k}}\cdot=\frac{\imath}{\hbar}[(\mathsf{J}\hat{\mathbf{x}})_{k},\cdot]$
with $k=1,\ldots,2n$. The matrices

$\displaystyle\mathbb{D}=\hbar\real(\mathbb{\Gamma})\quad\text{and}\quad\mathbb{C}=\imaginary(\mathbb{\Gamma}),$
(10)

are, respectively, the diffusion and dissipation matrices, both defined
through the decoherence matrix

$\mathbb{\Gamma}=\sum_{k=1}^{K}\mathbf{l}_{k}\mathbf{l}_{k}^{\dagger},$ (11)

which is composed by the vectors in the Lindblad operators (8b). The adjoint
generator $\bar{\mathcal{L}}_{\mathrm{NU}}$ in (4), using Eq.(5b) for the
present case, is

$\displaystyle\bar{\mathcal{L}}_{\mathrm{NU}}^{\mathrm{G}}[\hat{O}_{t}]$
$\displaystyle=$
$\displaystyle\tr\\!\left(\mathbb{C}\mathsf{J}\hat{\mathbf{x}}\partialderivative{\hat{\mathbf{x}}^{\intercal}}\hat{O}_{t}+\frac{1}{2}\mathbb{D}\partialderivative{\hat{\mathbf{x}}}\partialderivative{\hat{\mathbf{x}}^{\intercal}}\hat{O}_{t}\right).$
(12)

According to the definitions in (10), $\mathbb{D}=\mathbb{D}^{\intercal}\geq
0$ and $\mathbb{C}=-\mathbb{C}^{\intercal}$. Note also that, according to
(11), $\mathbb{\Gamma}\geq 0$, and thus

$\hbar\mathbb{\Gamma}=\mathbb{D}+\imath\hbar\mathbb{C}\geq 0,$ (13)

which can be interpreted as a generalized fluctuation–dissipation relation
Wiseman and Milburn (2009). For the next sections, a useful result concerning
a relation between the matrices in (13) is the following Lemma, which is
proved in the Appendix A.

###### Lemma 1

If $\det\mathbb{C}\neq 0$ in (13), i.e., $\mathbb{C}$ is invertible, then both
the diffusion and the decoherence matrices are invertible and strictly
positive-definite, that is, $\mathbb{D}>0$ and $\mathbb{\Gamma}>0$.

Two quantities of main importance for establishing the results of this work
are the mean-value vector

$\expectationvalue{\hat{\mathbf{x}}}_{t}=\Tr(\hat{\rho}_{t}\hat{\mathbf{x}})$
(14)

and the (dimensionless) covariance matrix

$\mathbb{V}_{\\!t}=\frac{1}{2\hbar}\Tr(\hat{\rho}_{t}(\hat{\mathbf{x}}-\expectationvalue{\hat{\mathbf{x}}}_{t})(\hat{\mathbf{x}}-\expectationvalue{\hat{\mathbf{x}}}_{t})^{\intercal}).$
(15)

Despite the evolution of the system state through a  GDS can be analytically
determined Carmichael (1999), the description for the system behavior is
improved when analyzing the evolution of these two moments. Taking the
temporal derivative of above equations and using the  LME in (1) for the  GDS,
i.e., with the operators in (8), the cyclicity of the trace together with the
canonical commutation relation yield Nicacio et al. (2016)

$\derivative{\expectationvalue{\hat{\mathbf{x}}}_{t}}{t}=\mathbb{A}\expectationvalue{\hat{\mathbf{x}}}_{t}-\boldsymbol{\xi},$
(16)

and

$\derivative{\mathbb{V}_{\\!t}}{t}=\left(\mathbb{A}\mathbb{V}_{\\!t}+\mathbb{V}_{\\!t}\mathbb{A}^{\\!\intercal}\right)+\frac{\mathbb{D}}{\hbar},$
(17)

where we defined the drift matrix

$\mathbb{A}=\mathsf{J}\mathbb{B}^{\prime}-\mathbb{C}\mathsf{J},$ (18)

for $\mathbb{B}^{\prime}$ from (8a) and $\mathbb{C}$ from (10).

By direct integration, the solutions of Eqs.(16) and (17) are, respectively,

$\displaystyle\\!\\!\expectationvalue{\hat{\mathbf{x}}}_{t}$ $\displaystyle=$
$\displaystyle
e^{\mathbb{A}t}\expectationvalue{\hat{\mathbf{x}}}_{0}-\int_{0}^{t}\\!dt^{\prime}e^{\mathbb{A}t^{\prime}}\boldsymbol{\xi},$
(19a) $\displaystyle\\!\\!\mathbb{V}_{\\!t}$ $\displaystyle=$ $\displaystyle
e^{\mathbb{A}t}\mathbb{V}_{0}e^{\mathbb{A}{\\!{}^{\intercal}}t}+\frac{1}{\hbar}\int_{0}^{t}\\!dt^{\prime}\;e^{\mathbb{A}(t-t^{\prime})}\;\mathbb{D}\;e^{\mathbb{A}^{\\!\intercal}(t-t^{\prime})}.$
(19b)

If the matrix $\mathbb{A}$ is invertible, the integral in (19a) can be
explicitly performed and this solution becomes

$\expectationvalue{\hat{\mathbf{x}}}_{t}=e^{\mathbb{A}t}(\expectationvalue{\hat{\mathbf{x}}}_{0}-\mathbb{A}^{-1}\boldsymbol{\xi})+\mathbb{A}^{-1}\boldsymbol{\xi}.$
(20)

###  II.1 Gaussian States

The formalism presented so far describes the action of a  GDS on a generic
quantum state. However, the term “Gaussian” in the acronym “ GDS” refers to
the fact that this kind of dynamics is a quantum channel that preserves the
Gaussian character of an initial Gaussian state throughout the whole
evolution.

The density operator $\hat{\sigma}_{t}$ of a Gaussian state can be expressed
as Holevo (2019); Banchi et al. (2015)

$\hat{\sigma}_{t}=\frac{{e^{-\frac{1}{2\hbar}(\hat{\mathbf{x}}-\expectationvalue{\hat{\mathbf{x}}}_{t})^{\intercal}\mathbb{U}_{t}(\hat{\mathbf{x}}-\expectationvalue{\hat{\mathbf{x}}}_{t})}}}{\sqrt{{\det}(\mathbb{V}_{\\!t}+\tfrac{\imath}{2}\mathsf{J})}},$
(21)

which is completely determined only by the moments in Eqs.(14) and (15), where
the mean-value is
$\expectationvalue{\hat{\mathbf{x}}}_{t}=\Tr(\hat{\sigma}_{t}\hat{\mathbf{x}})$
and the matrix $\mathbb{U}_{t}$ is given by

$\mathbb{U}_{t}=2\imath\mathsf{J}\coth^{-1}\left(2\imath\mathbb{V}_{\\!t}\mathsf{J}\right).$
(22)

Note that $\mathbb{V}_{\\!t}+\tfrac{\imath}{2}\mathsf{J}\geq 0$ is the bona
fide condition of a covariance matrix of a quantum state Simon et al. (1994),
thus the determinant in the denominator (21) is never negative. When subjected
to a  GDS, the evolved state is like (21) with
$\expectationvalue{\hat{\mathbf{x}}}_{t}$ and $\mathbb{V}_{\\!t}$ given in
(19).

The relation between the matrices in (22) can be strengthened, which will be
necessary for our future results. It is immediate from (22) that

$\mathsf{J}\mathbb{U}_{t}\mathbb{V}_{\\!t}=\mathbb{V}_{\\!t}\mathbb{U}_{t}\mathsf{J}\Leftrightarrow[\mathsf{J}\mathbb{U}_{t},\mathbb{V}_{\\!t}\mathsf{J}]=0;$
(23)

however, we will prove this relation using well-known results in order to
establish methods and notations for several future occasions. First, we use
the Williamson theorem Simon et al. (1994); Nicacio (2021a) which establishes
that for every $2n\times 2n$ real symmetric and positive-definite matrix
$\mathbb{V}_{\\!t}$, i.e.,
$\mathbb{V}_{\\!t}^{\intercal}=\mathbb{V}_{\\!t}>0$, there exists a symplectic
matrix $\mathsf{S}_{t}\in{\rm Sp}(2n,\mathds{R})$ such that

$\mathsf{S}_{t}\mathbb{V}_{\\!t}\mathsf{S}_{t}^{\intercal}=\mathbb{k}_{t}\oplus\mathbb{k}_{t},$
(24)

where $\mathbb{k}_{t}=\mbox{diag}(\kappa_{1}(t),\ldots,\kappa_{n}(t))$ is the
symplectic spectra of $\mathbb{V}_{\\!t}$ and $\kappa_{j}(t)\geq
1/2\,(j=1,...,n)$ are the symplectic eigenvalues. Next, we use the following
Lemma, also a consequence of the Williamson theorem.

###### Lemma 2

A Hamiltonian matrix222A real $2n\times 2n$ matrix $\mathbb{M}$ is said
Hamiltonian matrix iff $\mathsf{J}\mathbb{M}$ (or equivalently
$\mathbb{M}\mathsf{J}$) is symmetric, where $\mathsf{J}$ is in (7).
$\mathbb{O}\mathsf{J}$, where $\mathbb{O}$ is symmetric and positive-definite
and $\mathsf{J}$ is in (7), is diagonalized by the similarity transformation

$(\mathbb{Q}\mathsf{S})\,\mathbb{O}\mathsf{J}\,(\mathbb{Q}\mathsf{S})^{-1}=(\imath\mathbb{o})\oplus(-\imath\mathbb{o}),$
(25)

where $\mathbb{o}={\rm diag}(o_{1},\ldots,o_{n})$, $o_{j}>0$ ($j=1,\ldots,n$)
are the symplectic eigenvalues of $\mathbb{O}$ through $\mathsf{S}$, i.e.,
$\mathsf{S}\mathbb{O}\mathsf{S}^{\intercal}=\mathbb{o}\oplus\mathbb{o}$, and
$\mathbb{Q}$ is the complex matrix 333The matrix $\mathbb{Q}$ is a member of
the compact symplectig group ${\rm Sp}(n):={\rm Sp}(2n,\mathds{C})\;\cap\;{\rm
SU}(2n)$.

$\mathbb{Q}=\mathbb{Q}^{\intercal}=\frac{1}{\sqrt{2}}\begin{pmatrix}\mathbb{1}&-\imath\mathbb{1}\\\
-\imath\mathbb{1}&\mathbb{1}\\\ \end{pmatrix}$ (26)

such that $\mathbb{Q}^{\intercal}\mathsf{J}\mathbb{Q}=\mathsf{J}$ and
$\mathbb{Q}^{-1}=\mathbb{Q}^{\dagger}$. Equivalently, we have

$(\mathbb{Q}\mathsf{S}^{-\intercal})\,\mathsf{J}\mathbb{O}\,(\mathbb{Q}\mathsf{S}^{-\intercal})^{-1}=(\imath\mathbb{o})\oplus(-\imath\mathbb{o}).$
(27)

Returning to the proof of Eq.(23), the above Lemma can be used to diagonalize
the matrix $\mathbb{V}_{\\!t}\mathsf{J}$, i.e.,

$\mathbb{Q}\mathsf{S}_{t}\,\mathbb{V}_{\\!t}\mathsf{J}(\mathbb{Q}\mathsf{S}_{t})^{-1}=(\imath\mathbb{k}_{t})\oplus(-\imath\mathbb{k}_{t})=:(\mathbb{V}_{\\!t}\mathsf{J})_{\rm
d}.$ (28)

As a useful notation, the diagonal matrix $(\mathbb{V}_{\\!t}\mathsf{J})_{\rm
d}$ defined above will be called the canonical form of
$\mathbb{V}_{\\!t}\mathsf{J}$. From (22), we write
$\mathsf{J}\mathbb{U}_{t}=-2\imath\,g(2\imath\mathbb{V}_{\\!t}\mathsf{J})$,
where

$g(x)=2\coth^{-1}(2x)>0$ (29)

is a continuous function for $x>1/2$. Consequently, employing Eq.(28) and
noting that $g(x)=-g(-x)$, we attain

$\displaystyle\mathbb{Q}\mathsf{S}_{t}\,\mathsf{J}\mathbb{U}_{t}\,(\mathbb{Q}\mathsf{S}_{t})^{-1}=[\imath
g(\mathbb{k}_{t})]\oplus[-\imath
g(\mathbb{k}_{t})]=:(\mathsf{J}\mathbb{U}_{t})_{\rm d},$ (30)

which is the canonical form of the Hamiltonian matrix
$\mathsf{J}\mathbb{U}_{t}$. Therefore, the matrix
$(\mathbb{Q}\mathsf{S}_{t})^{-1}$ simultaneously diagonalizes the matrices
$\mathbb{V}_{\\!t}\mathsf{J}$ and $\mathsf{J}\mathbb{U}_{t}$ and they must
commute, as we wanted to prove.

Noteworthy, the matrix $\mathbb{U}_{t}$ is positive-definite so the Williamson
theorem can be applied. The symplectic diagonalization can be easily inferred
using (27) with $\mathbb{O}=\mathbb{U}_{t}$ and
$\mathsf{S}^{-\intercal}=\mathsf{S}_{t}$, so we get the canonical form in
(30). Therefore according to the hypothesis of the Lemma the matrix that
diagonalizes symplectically $\mathbb{U}_{t}$ is $\mathsf{S}_{t}^{-\intercal}$,
i.e.,

$\mathsf{S}_{t}^{-\intercal}\mathbb{U}_{t}\mathsf{S}_{t}^{-1}=g(\mathbb{k}_{t})\oplus
g(\mathbb{k}_{t}),$ (31)

where
$g(\mathbb{k}_{t})=\mbox{diag}(g(\kappa_{1}(t)),\ldots,g(\kappa_{n}(t)))$ and
$g(\kappa_{j})>0\,(j=1,...,n)$ are the symplectic eigenvalues of
$\mathbb{U}_{t}$.

From (31), we see that the matrix $\mathbb{U}_{t}$ is finite whenever
$\kappa_{j}(t)>1/2$ for all $j=1,\ldots,n$ and for each fixed value of $t$. In
this case the density operator $\hat{\sigma}_{t}$ in (21) corresponds to a
full-rank mixed-state. Still, the representation of Gaussian states as in (21)
is also valid in the limit $\kappa_{j}(t)\to 1/2\,\forall j$, where both the
matrix $\mathbb{U}_{t}$ and
${\det}(\mathbb{V}_{\\!t}+\tfrac{\imath}{2}\mathsf{J})$ diverge. In this
limit, we have $\hat{\sigma}_{t}\rightarrow\outerproduct{\Psi_{t}}{\Psi_{t}}$,
where $\ket{\Psi_{t}}$ is an $n$-mode pure-Gaussian-state. When some but not
all symplectic eigenvalues are such that $\kappa_{j}(t)=1/2$, the same
divergences happen, and $\hat{\sigma}_{t}$ in (21) represents a rank-deficient
mixed quantum-state in the limit $\kappa_{j}(t)\to 1/2$. In conclusion,
$\hat{\sigma}_{t}$ in (21) is a valid representation of the density operator
of any Gaussian state.

###  II.2 Wigner Function and Fokker-Planck Equation

In continuous variable systems, the following sets
$\\{\hat{T}_{\boldsymbol{\xi}}=e^{\frac{\imath}{\hbar}\hat{\mathbf{x}}^{\intercal}\mathsf{J}\boldsymbol{\xi}}\,\,|\,\,\boldsymbol{\xi}\\!\in\\!\mathds{R}^{2n}\\}$
and $\\{\hat{R}_{\mathbf{x}}=(4\pi\hbar)^{-n}\int
d\boldsymbol{\xi}\,e^{\frac{\imath}{\hbar}\boldsymbol{\xi}^{\intercal}\mathsf{J}\mathbf{x}}\,\hat{T}_{\boldsymbol{\xi}}\,\,|\,\,\mathbf{x}\\!\in\\!\mathds{R}^{2n}\\}$
are basis of a vector space constituted by operators acting on the separable
infinite-dimensional Hilbert space ${\cal H}=\otimes_{j=1}^{n}{\cal H}_{j}$ of
the $n$-mode bosonic system Nicacio (2021b). The elements of these sets are
called translations and reflections, respectively, and refer to their action
on the vector operator $\hat{\mathbf{x}}$ in (6), corresponding to the
Heisenberg picture, namely
$\hat{T}_{\boldsymbol{\xi}}^{\dagger}\,\hat{\mathbf{x}}\,\hat{T}_{\boldsymbol{\xi}}=\hat{\mathbf{x}}+\boldsymbol{\xi}\hat{1}$
and
$\hat{R}_{\mathbf{x}}^{\dagger}\,\hat{\mathbf{x}}\,\hat{R}_{\mathbf{x}}=-\hat{\mathbf{x}}+2\mathbf{x}\hat{1}$
Ozorio de Almeida (1998). Translation operators are unitary
$\hat{T}^{\dagger}=\hat{T}_{\boldsymbol{\xi}}^{-1}=\hat{T}_{-\boldsymbol{\xi}}$
and reflection operators are unitary and Hermitian, i.e., involutory
operators, $\hat{R}_{\mathbf{x}}^{2}=\hat{1}$. The operator
$\hat{T}_{\boldsymbol{\xi}}$ is also known as Weyl operator and
$\hat{R}_{\mathbf{x}}$ as Wigner operator Tarasov (2008).

When dealing with continuous variable systems, the existence of unbounded
operators and operators with continuous spectra may cause some mathematical
difficulties. In particular, it is often difficult to find the algebra of
operators that defines the domain of applicability of a given formalism. In
this work, we circumvent this difficulty by applying our formalism to the
algebra of all operators with a Weyl and Wigner representation. This will be
particularly important for the demonstration developed in Appendixes B and C.

The Weyl and Wigner representations of an operator $\hat{A}$ are,
respectively, the Hilbert-Schimidt inner products
$A(\boldsymbol{\xi})=\langle\hat{A},\hat{T}_{\boldsymbol{\xi}}\rangle$ and
$A(\mathbf{x})=\langle\hat{A},\hat{R}_{\mathbf{x}}\rangle$, which are the
coefficients of the expansion (also called symbols of $\hat{A}$) in one of the
mentioned bases through the Bochner integrals de Gosson (2006):

$\hat{A}=\int\frac{d\boldsymbol{\xi}}{(2\pi\hbar)^{n}}\,A(\boldsymbol{\xi})\,\hat{T}_{\boldsymbol{\xi}}=\int\frac{d\mathbf{x}}{(\pi\hbar)^{n}}\,A(\mathbf{x})\,\hat{R}_{\mathbf{x}}.$
(32)

These facts are consequences of the orthogonality relations
$\langle\hat{T}_{\boldsymbol{\xi}},\hat{T}_{\boldsymbol{\xi}^{\prime}}\rangle=2^{n}\langle\hat{R}_{\boldsymbol{\xi}},\hat{R}_{\boldsymbol{\xi}^{\prime}}\rangle=(2\pi\hbar)^{n}\delta(\boldsymbol{\xi}^{\prime}-\boldsymbol{\xi})$.
In particular, the Wigner representation of the operator vector
$\hat{\mathbf{x}}$ in (6),

$\mathbf{x}=\expectationvalue*{\hat{\mathbf{x}},\hat{R}_{\mathbf{x}}}=(q_{1},\ldots,q_{n},p_{1},\ldots,p_{n})^{\intercal},$
(33)

is a real vector in phase-space. In the following, we will use the Wigner
representation of the master equation for a  GDS, which is nothing more than
an alternative description of the system evolution. Through this
representation we will establish our first result in Theorem 1 in Section
III.2. Also, this will be important to prove the results in Appendixes B and
C.

The Wigner representation of the  LME in (1) for the operators in (8), see
Carmichael (1999); Nicacio et al. (2010); Wiseman and Milburn (2009), is the
Fokker-Planck equation Risken (1996)

$\derivative{W(\mathbf{x},t)}{t}=-\partialderivative{\mathbf{x}^{\intercal}}\,[{\bf
v}_{\mathrm{U}}(\mathbf{x},t)+{\bf v}_{\mathrm{NU}}(\mathbf{x},t)]$ (34)

for the Wigner function444 The Wigner function is proportional to the Wigner
symbol of the density operator, this proportionality guarantees the
normalization of that function as a quasi-probability density, see Carmichael
(1999); Nicacio (2021b); Ozorio de Almeida (1998); de Gosson (2006), for
instance.
$W(\mathbf{x})=\frac{1}{(\pi\hbar)^{n}}\expectationvalue*{\hat{\rho},\hat{R}_{\mathbf{x}}}$.
In the above equation, we identify the Fokker-Planck current vectors:

$\displaystyle{\bf v}_{\mathrm{U}}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle(\mathsf{J}\mathbb{B}^{\prime}\;\mathbf{x}-\boldsymbol{\xi}^{\prime})\,W({\bf
x},t),\,\,\,$ (35a) $\displaystyle{\bf v}_{\mathrm{NU}}(\mathbf{x},t)$
$\displaystyle=$
$\displaystyle-\tfrac{1}{2}\mathbb{D}\frac{\partial}{\partial{\bf x}}W({\bf
x},t)-\mathbb{C}\mathsf{J}\mathbf{x}\;,$ (35b)

corresponding, respectively, to the unitary (reversible) and non-unitary
(irreversible) contributions to the evolution of $W(\bf x)$. Note that the
first term in Eq.(34) is the Poisson bracket

$-\partialderivative{\mathbf{x}^{\intercal}}\,[{\bf
v}_{\mathrm{U}}(\mathbf{x},t)]=[H_{\mathrm{eff}},W(\mathbf{x},t)]_{\rm cl}$
(36)

between the Wigner function and the Hamiltonian
$H_{\mathrm{eff}}=\frac{1}{2}\mathbf{x}^{\intercal}\mathbb{B}^{\prime}\mathbf{x}+\mathbf{x}^{\intercal}\mathsf{J}\boldsymbol{\xi}^{\prime}$,
which is the Wigner symbol of (8a).

When the initial state is a Gaussian state, the evolved Wigner function of
$\hat{\sigma}_{t}$ in (21) is

$W_{\rm
G}(\mathbf{x},t)=\frac{\expectationvalue*{\hat{\sigma}_{t},\hat{R}_{\mathbf{x}}}}{(\pi\hbar)^{n}}=\frac{e^{-\frac{1}{2\hbar}(\mathbf{x}-\expectationvalue{\hat{\mathbf{x}}}_{t})^{\intercal}\mathbb{V}_{\\!t}^{-1}(\mathbf{x}-\expectationvalue{\hat{\mathbf{x}}}_{t})}}{(2\pi\hbar)^{n}\sqrt{{\det}\mathbb{V}_{\\!t}}},$
(37)

for $\expectationvalue{\hat{\mathbf{x}}}_{t}$ and $\mathbb{V}_{\\!t}$ in (19).
This function is a multivariate Gaussian distribution, so the Fokker-Planck
currents in (35) are true probability currents given by

$\displaystyle\\!\\!\\!\\!\\!\\!{\bf v}^{\rm G}_{\mathrm{U}}(\mathbf{x},t)$
$\displaystyle=$
$\displaystyle(\mathsf{J}\mathbb{B}^{\prime}\;\mathbf{x}-\boldsymbol{\xi}^{\prime})\,W_{\rm
G}({\bf x},t),\,\,\,$ (38a) $\displaystyle\\!\\!\\!\\!\\!\\!{\bf v}^{\rm
G}_{\mathrm{NU}}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle[\tfrac{1}{2\hbar}\mathbb{D}\mathbb{V}_{\\!t}^{-1}(\mathbf{x}-\expectationvalue{\hat{\mathbf{x}}}_{t})-\mathbb{C}\mathsf{J}\mathbf{x}]W_{\rm
G}(\mathbf{x},t).$ (38b)

##  III Thermal equilibrium in GDS: general considerations

For  GDSs, if there exists a stationary state $\hat{\sigma}^{\rm S}$, it will
be unique for any initial state $\hat{\rho}_{0}$, i.e.,
$\lim_{t\to\infty}\hat{\rho}_{t}=\hat{\sigma}^{\rm S}$ Frigerio (1977, 1978);
Carmichael (1999). In particular, starting with an initial Gaussian state, the
evolved state remains Gaussian throughout the whole evolution, therefore
$\hat{\sigma}^{\rm S}$ is necessarily Gaussian.

The first moments and covariance matrix of the stationary state can be
determined through the asymptotic behavior of Eqs.(19) and for that we resort
to the Lyapunov theory of stability Horn and Johnson (1991). Note that the
only way to erase all the information about any initial condition in Eqs.(19)
is to admit a matrix $\mathbb{A}$ with all of its eigenvalues with negative
real part, which is the same as saying that $\mathbb{A}$ is a Hurwitz matrix.
Consequently, the covariance matrix of $\hat{\sigma}^{\rm S}$ is a solution of
the Lyapunov equation

$\displaystyle\derivative{\mathbb{V}^{\rm
S}}{t}=\left(\mathbb{A}\mathbb{V}^{\rm S}+\mathbb{V}^{\rm
S}\mathbb{A}^{\\!\intercal}\right)+\frac{\mathbb{D}}{\hbar}=0.$ (39)

Recalling that $\mathbb{V}^{\rm S}$ is strictly positive-definite
$\mathbb{V}^{\rm S}>0$ and together with the Hurwitz condition over
$\mathbb{A}$, the Lyapunov theorem Horn and Johnson (1991) sets that
$\mathbb{D}>0$. In this case, the linear dynamical system in (16) is said
globally asymptotically stable (AS)Horn and Johnson (1991) and the solutions
in (19) attains the asymptotic values

$\displaystyle\expectationvalue{\hat{\mathbf{x}}}^{\rm S}$ $\displaystyle=$
$\displaystyle\mathbb{A}^{-1}\boldsymbol{\xi}^{\prime},$ (40a)
$\displaystyle\mathbb{V}^{\rm S}$ $\displaystyle=$
$\displaystyle\frac{1}{\hbar}\int_{0}^{\infty}dt\;e^{\mathbb{A}t}\;\mathbb{D}\;e^{\mathbb{A}^{\intercal}t},$
(40b)

where it is clear that any trace of the initial state disappears.

###  III.1 Gibbs States as Stationary States

A stationary Gaussian state in a  GDS corresponds to a thermal equilibrium
state when $\hat{\sigma}^{\rm S}=\hat{\sigma}^{\rm{th}}$ is the Gibbs state

$\hat{\sigma}^{\mathrm{th}}=\frac{\hat{\bar{\sigma}}^{\mathrm{th}}}{{\cal
Z}^{\mathrm{th}}},\quad\hat{\bar{\sigma}}^{\mathrm{th}}=e^{-\beta\hat{H}},\quad{\cal
Z}^{\mathrm{th}}={\Tr}(\hat{\bar{\sigma}}^{\mathrm{th}}),$ (41)

where $\beta$ is the “inverse temperature” and
$\hat{H}=\hat{H}_{\mathrm{eff}}$ is the quadratic Hamiltonian in (8a) or, more
generically, another quadratic Hamiltonian such that
$[\hat{H},\hat{H}_{\mathrm{eff}}]=0$. In the following we will establish
necessary conditions over the quadratic Hamiltonians, $\hat{H}_{\mathrm{eff}}$
and $\hat{H}$, which allows a  GDS to have an equilibrium thermal state.

Without loss of generality, we can set the origin of the phase-space
coordinates $\mathbf{x}$ in (33) such
$\expectationvalue{\hat{\mathbf{x}}}^{\mathrm{th}}=\Tr(\hat{\sigma}^{\mathrm{th}}\hat{\mathbf{x}})=-\mathbb{B}^{-1}\mathsf{J}\boldsymbol{\xi}=0$,
where the vector $\boldsymbol{\xi}$ is associated with a possible linear term
of $\hat{H}$. Therefore, the quadratic Hamiltonian $\hat{H}$ can be chosen as

$\hat{H}=\tfrac{1}{2}\hat{\mathbf{x}}^{\intercal}\mathbb{B}\hat{\mathbf{x}},$
(42)

i.e., with $\boldsymbol{\xi}=0$. The Hessian matrix of the Hamiltonian has to
be positive definite, $\mathbb{B}>0$, in order to fulfill the normalization
condition ${\Tr}(\hat{\sigma}^{\mathrm{th}})=1$ Nicacio (2021b). From
$\hat{H}_{\mathrm{eff}}$ in (8a) and $\hat{H}$ in (42), we have

$\displaystyle[\hat{H},\hat{H}_{\mathrm{eff}}]$ $\displaystyle=$
$\displaystyle-\tfrac{\imath\hbar}{2}\hat{\mathbf{x}}^{\intercal}\mathsf{J}[\mathsf{J}\mathbb{B},\mathsf{J}\mathbb{B}^{\prime}]\hat{\mathbf{x}}-\imath\hbar(\boldsymbol{\xi}^{\prime})^{\intercal}\mathbb{B}\hat{\mathbf{x}}=0$
(43) $\displaystyle\Leftrightarrow$
$\displaystyle[\mathsf{J}\mathbb{B},\mathsf{J}\mathbb{B}^{\prime}]=0\;\;\text{and}\;\;\boldsymbol{\xi}^{\prime}=0.$

Therefore, the Hamiltonian of the free evolution of a  GDS with an equilibrium
thermal state must also be of the form

$\hat{H}_{\mathrm{eff}}=\tfrac{1}{2}\hat{\mathbf{x}}^{\intercal}\mathbb{B}^{\prime}\hat{\mathbf{x}}.$
(44)

Comparing the general form of the density operator of a Guassian state,
Eq.(21), with the thermal state $\hat{\sigma}^{\mathrm{th}}$ in (41) with
$\hat{H}$ in (42), we arrive to

$\mathbb{U}^{\mathrm{th}}=\hbar\beta\mathbb{B}.$ (45)

Applying condition (23), we get

$\mathsf{J}\mathbb{B}\mathbb{V}^{\mathrm{th}}=\mathbb{V}^{\mathrm{th}}\mathbb{B}\mathsf{J}\Leftrightarrow[\mathsf{J}\mathbb{B},\mathbb{V}^{\mathrm{th}}\mathsf{J}]=0,$
(46)

where $\mathbb{V}^{\mathrm{th}}$ is the covariance matrix of the Gibbs
Gaussian state in (41). Using (22), it is clear that
$\mathbb{V}^{\mathrm{th}}\mathsf{J}$ is a function of $\mathsf{J}\mathbb{B}$,
viz,

$\mathbb{V}^{\mathrm{th}}\mathsf{J}=-\frac{\imath}{2}\coth\left(\frac{\imath\hbar\beta}{2}\mathsf{J}\mathbb{B}\right).$
(47)

Taking into account the relation
$[\mathsf{J}\mathbb{B},\mathsf{J}\mathbb{B}^{\prime}]=0$ from (43) and that
$\mathbb{V}^{\mathrm{th}}\mathsf{J}$ is a function of $\mathsf{J}\mathbb{B}$,
it is also true that

$[\mathsf{J}\mathbb{B}^{\prime},\mathbb{V}^{\mathrm{th}}\mathsf{J}]=0\Leftrightarrow\mathsf{J}\mathbb{B}^{\prime}\mathbb{V}^{\mathrm{th}}=\mathbb{V}^{\mathrm{th}}\mathbb{B}^{\prime}\mathsf{J}.$
(48)

According to Williamson theorem it is possible to find a symplectic matrix
${\mathsf{S}}^{\mathrm{th}}$ such that

$\displaystyle({\mathsf{S}}^{\mathrm{th}})\mathbb{V}^{\mathrm{th}}({\mathsf{S}}^{\mathrm{th}})^{\intercal}=\mathbb{k}^{\mathrm{th}}\oplus\mathbb{k}^{\mathrm{th}},$
(49)

and using Eq.(25) of Lemma 2, the matrix
$(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})$ diagonalizes
$\mathbb{V}^{\mathrm{th}}\mathsf{J}$:

$\mathbb{Q}{\mathsf{S}}^{\mathrm{th}}\,\mathbb{V}^{\mathrm{th}}\mathsf{J}(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{-1}=(\imath\mathbb{k}^{\mathrm{th}})\oplus(-\imath\mathbb{k}^{\mathrm{th}}).$
(50)

However, due to (47), the same matrix $(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})$
also diagonalizes $\mathsf{J}\mathbb{B}$. Using Eq.(31) one realizes that
${({{\mathsf{S}}^{\mathrm{th}}})}^{-\intercal}$ is the symplectic matrix that
diagonalizes $\mathbb{B}$, i.e.,

${({{\mathsf{S}}^{\mathrm{th}}})}^{-\intercal}\mathbb{B}({\mathsf{S}}^{\mathrm{th}})^{-1}=\mathbb{w}\oplus\mathbb{w},\,\,\,\mathbb{w}=\text{diag}(\omega_{1},\ldots,\omega_{n}),$
(51)

where $\omega_{j}>0$ ($j=1,\ldots,n$) are the symplectic eigenvalues of the
Hessian matrix $\mathbb{B}>0$ in (42). It is also possible to define, again
according to Lemma 2, the canonical form of $\mathsf{J}\mathbb{B}$:

$(\mathsf{J}\mathbb{B})_{\rm
d}=(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})\,\mathsf{J}\mathbb{B}\,(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{-1}=(\imath\mathbb{w})\oplus(-\imath\mathbb{w}).$
(52)

Finally, from (47), the relations between the symplectic spectra of
$\mathbb{V}^{\mathrm{th}}$ and $\mathbb{B}$ is

$\displaystyle\mathbb{k}^{\mathrm{th}}=\frac{1}{2}\coth\left(\frac{\hbar\beta\mathbb{w}}{2}\right),$
(53)

or equivalently $\mathbb{w}=\frac{1}{\hbar\beta}g(\mathbb{k}^{\mathrm{th}})$,
see Eq.(29).

For a phase space described by $\mathbf{x}$ in (33), the classical counterpart
of the Hamiltonian $\hat{H}$ in (42) coincides with its Wigner symbol, i.e.,

$H=\tfrac{1}{2}\mathbf{x}^{\intercal}\mathbb{B}\mathbf{x}=\langle\hat{H},\hat{R}_{\mathbf{x}}\rangle$
(54)

and the solution of the Hamilton equation
$\dot{\mathbf{x}}=\mathsf{J}\partialderivative{H}{\mathbf{x}}$ is given by
$\mathbf{x}(t)=\tilde{\mathsf{S}}_{t}\mathbf{x}(0)$ with

$\tilde{\mathsf{S}}_{t}=e^{\mathsf{J}\mathbb{B}t}=(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{-1}\,e^{(\mathsf{J}\mathbb{B})_{\rm
d}t}\,(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})$ (55)

and $(\mathsf{J}\mathbb{B})_{\rm d}$ in (52). The matrix
$\tilde{\mathsf{S}}_{t}$ defined above will be important in Section III.2, but
here it is worth to note that it generates a Hamiltonian flow around the
elliptical fixed point $\mathbf{x}=0$ and that, from (51), $\omega_{j}$ are
the eigenfrequencies of the Hamiltonian (54). This is a direct consequence of
the positive-definiteness of $\mathbb{B}$ and, by this reason, we call
positive-elliptic all the Hamiltonians in (42) with $\mathbb{B}>0$.
Consequently, a positive-elliptic Hamiltonian in (42) is a necessary condition
for a  GDS to have an $n$-mode equilibrium thermal-state, since it is
necessary for the convergence of the partition function ${\cal
Z}^{\mathrm{th}}$ in (41).

Up to this point we were describing some properties of thermal states
associated to quadratic Hamiltonians. In the next section, we will show the
conditions over the diffusion and dissipation matrices $\mathbb{D}$ and
$\mathbb{C}$, respectively, which define a  GDS with a thermal equilibrium
state.

###  III.2 Diffusion and dissipation matrices for thermal equilibrium

The Lyapunov equation (39) for the covariance matrix of a  GDS thermal-
equilibrium-state $\mathbb{V}^{\mathrm{th}}$, where $\mathbb{A}$ is given in
(18), attains a simpler form through condition (48):

$\mathbb{C}\mathsf{J}\mathbb{V}^{\mathrm{th}}+\mathbb{V}^{\mathrm{th}}\mathsf{J}\mathbb{C}=\frac{\mathbb{D}}{\hbar},$
(56)

whose unique formal solution, see Eq.(40b), is

$\mathbb{V}^{\mathrm{th}}=\frac{1}{\hbar}\int_{0}^{\infty}dt\;e^{-\mathbb{C}\mathsf{J}t}\;\mathbb{D}\;e^{-\mathsf{J}\mathbb{C}t},$
(57)

where according to Lyapunov theorem Horn and Johnson (1991), the matrix
$-\mathbb{C}\mathsf{J}$ must be Hurwitz, since $\mathbb{V}^{\mathrm{th}}>0$
and $\mathbb{D}>0$. Therefore, the dissipation matrix $\mathbb{C}$ must be
invertible and according to Lemma 1, $\mathbb{D}>0$ and $\mathbb{\Gamma}>0$,
i.e., the diffusion and decoherence matrices, in (10) and (11) respectively,
of a  QDS with a thermal equilibrium state must be positive definite. In the
following we show that an explicit solution of the integral (57) can be
obtained from the stationary condition over the Fokker-Planck equation (34)
corresponding to a  QDS with a thermal equilibrium state.

When a  GDS has a stationary state, this is a Gaussian state
$\hat{\sigma}^{s}$ and from the Fokker-Planck equation (34), we obtain the
condition

$-\partialderivative{\mathbf{x}^{\intercal}}\,[{\bf
v}^{G}_{\mathrm{U}}(\mathbf{x})+{\bf v}^{G}_{\mathrm{NU}}(\mathbf{x})]=0,$
(58)

with ${\bf v}^{G}_{\mathrm{U}}(\mathbf{x})$ and ${\bf
v}^{G}_{\mathrm{NU}}(\mathbf{x})$ in Eqs.(38). If this stationary state is a
thermal equilibrium state, $\hat{\sigma}^{s}=\hat{\sigma}^{\mathrm{th}}$,
condition (58) simplifies to

$-\partialderivative{\mathbf{x}^{\intercal}}\,[{\bf
v}^{\mathrm{th}}_{\mathrm{NU}}(\mathbf{x})]=0,$ (59)

since, due to Eq.(36), one has

$\displaystyle[H_{\mathrm{eff}}(\mathbf{x}),W^{\mathrm{th}}(\mathbf{x})]_{\mathrm{cl}}$
$\displaystyle=$
$\displaystyle\tr\left(\mathsf{J}\mathbb{B}^{\prime}\mathbf{x}\mathbf{x}^{\intercal}(\mathbb{V}^{\mathrm{th}})^{-1}\right)=0,$
(60)

where we employed Eq.(48) and the fact that
$\tr(\mathbb{A})=\tr(\mathbb{A}^{\intercal})$ for any matrix $\mathbb{A}$. Now
we can establish the following theorem that characterizes a  QDS with a
thermal equilibrium state:

###### Theorem 1

A  QDS has a thermal equilibrium state iff

${\bf v}^{\mathrm{th}}_{\mathrm{NU}}(\mathbf{x})=0.$ (61)

The covariance matrix of such state is given by

$\mathsf{J}\mathbb{V}^{\mathrm{th}}=\frac{1}{2\hbar}\mathsf{J}\mathbb{D}(\mathsf{J}\mathbb{C})^{-1},$
(62)

where

$[\mathsf{J}\mathbb{D},\mathsf{J}\mathbb{C}]=0.$ (63)

In order to prove our theorem, we use the Divergence Theorem555See, for
example, Appendix A of Valero-Toranzo et al. (2018). and (59), both enable us
to relate the divergence of the vector field ${\bf
v}^{\mathrm{th}}_{\mathrm{NU}}(\mathbf{x})$ with the flux through the boundary
$\partial\Omega$ of the region $\Omega\in\mathds{R}^{2n}$,

$\int_{\partial\Omega}\,{\bf n}^{\intercal}{\bf
v}^{\mathrm{th}}_{\mathrm{NU}}(\mathbf{x})\,ds=\int_{\Omega}\partialderivative{\mathbf{x}^{\intercal}}\,[{\bf
v}^{\mathrm{th}}_{\mathrm{NU}}(\mathbf{x})]\,d\mathbf{x}^{2n}=0,$ (64)

where ${\bf n}$ is the $2n-$dimensional real vector normal to the surface
$\partial\Omega$. Since $\Omega$ has arbitrary volume, the necessary and
sufficient condition in (61) is proved. From (38b), ${\bf
v}^{\mathrm{th}}_{\mathrm{NU}}(\mathbf{x})=\left(\frac{1}{2\hbar}\mathbb{D}(\mathbb{V}^{\mathrm{th}})^{-1}-\mathbb{C}\mathsf{J}\right)\mathbf{x}\,W^{\mathrm{th}}(\mathbf{x})=0$
for any $\mathbf{x}$, ending up with Eq.(62). The relation in Eq.(63) follows
from the requirement
$\mathbb{V}^{\mathrm{th}}=(\mathbb{V}^{\mathrm{th}})^{\intercal}$. Note that,
if $[\mathsf{J}\mathbb{D},\mathsf{J}\mathbb{C}]=0$ then
$\mathbb{D}\;e^{-\mathsf{J}\mathbb{C}t}=e^{-\mathbb{C}\mathsf{J}t}\;\mathbb{D}$
and the integration in (57) can be explicitly performed to obtain exactly the
expression in (62).

It is worth to note that the covariance matrix $\mathbb{V}^{\mathrm{th}}$ of
the thermal equilibrium state $\hat{\sigma}^{\mathrm{th}}$ in (41) is
completely determined by the Hessian $\mathbb{B}$ of the Hamiltonian $\hat{H}$
in (42). This is shown in Eqs.(49) and (53), where
${({{\mathsf{S}}^{\mathrm{th}}})}^{-\intercal}$ symplectically diagonalizes
$\mathbb{B}$, whose symplectic spectrum is contained in the diagonal matrix
$\mathbb{w}$. So, the expression in Eq.(62) simply establishes the connection
between the fixed matrix $\mathbb{V}^{\mathrm{th}}$ and the dynamics of the
GDS, that is, the one determined by the matrices $\mathbb{C}$ and
$\mathbb{D}$, which at the end determines $\hat{\sigma}^{\mathrm{th}}$ as an
equilibrium state. However, the Lyapunov equation (56) has common solutions
Bialas and Gora (2015), i.e., there are different matrices $\mathbb{D}$ and
$\mathbb{C}$ which are able to give the same covariance matrix
$\mathbb{V}^{\mathrm{th}}$ in (62). Each pair $(\mathbb{D},\mathbb{C})$
corresponds to a different  GDS which has as steady state the same thermal
state $\hat{\sigma}^{\mathrm{th}}$. In the next section we show that the  QDBC
determines common solutions.

Also, according to Theorem 1, a  QDS has a thermal equilibrium state iff the
set
$\\{\mathsf{J}\mathbb{V}^{\mathrm{th}},\mathsf{J}\mathbb{D},\mathsf{J}\mathbb{C}\\}$
is a commuting set of matrices. Therefore, there is a matrix that
simultaneously diagonalizes the three matrices
$\mathsf{J}\mathbb{V}^{\mathrm{th}}$, $\mathsf{J}\mathbb{D}$, and
$\mathsf{J}\mathbb{C}$ (see Appendix H). Using Lemma 2, the Hamiltonian matrix
$\mathsf{J}\mathbb{V}^{\mathrm{th}}$ is diagonalized by
$(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{\dagger}$ so, for a  GDS with a
thermal equilibrium state, we can always write

$\mathsf{J}\mathbb{V}^{\mathrm{th}}=(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{\dagger}\;(\mathsf{J}\mathbb{V}^{\mathrm{th}})_{\rm
d}\;((\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{\dagger})^{-1},$ (65)

with $(\mathsf{J}\mathbb{V}^{\mathrm{th}})_{\rm
d}=\imath\mathbb{k}^{\mathrm{th}}\oplus\left(-\imath\mathbb{k}^{\mathrm{th}}\right)$
and $\mathbb{k}^{\mathrm{th}}$ in (53). Equivalently, applying the same Lemma
to the diagonalization of $\mathsf{J}\mathbb{D}$ we arrive at

$\mathsf{J}\mathbb{D}=(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{\dagger}\;(\mathsf{J}\mathbb{D})_{\rm
d}\;((\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{\dagger})^{-1},$ (66)

where $(\mathsf{J}\mathbb{D})_{\rm
d}=(\imath\hbar\mathbb{d})\oplus(-\imath\hbar\mathbb{d})$ with
$\hbar\mathbb{d}$ the $\mathds{R}^{n\times n}$ diagonal matrix with the
symplectic spectrum of $\mathbb{D}$ in its diagonal666It is worth to note that
${\mathsf{S}}^{\mathrm{th}}$ simultaneously diagonalizes simplectically
$\mathbb{V}^{\mathrm{th}}$ and $\mathbb{D}$.. Because $\mathsf{J}\mathbb{C}$
is a skew-Hamiltonian matrix 777A real $2n\times 2n$ matrix $\mathbb{M}$ is
said a skew-Hamiltonian matrix iff $\mathsf{J}\mathbb{M}$ (or equivalently
$\mathbb{M}\mathsf{J}$) is a skew-symmetric matrix., its eigenvalues are real
and with at least multiplicity equal to two Fassbender et al. (1999). So we
can write

$\mathsf{J}\mathbb{C}=(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{\dagger}\;(\mathsf{J}\mathbb{C})_{\rm
d}\;((\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{\dagger})^{-1},$ (67)

where $(\mathsf{J}\mathbb{C})_{\rm d}=\mathbb{jc}\oplus\mathbb{jc}$ with
$\mathbb{jc}$ a $\mathds{R}^{n\times n}$ diagonal matrix. Therefore, using
(62) we arrive to $\frac{1}{2\hbar}(\mathsf{J}\mathbb{D})_{\rm
d}(\mathsf{J}\mathbb{C})_{\rm
d}^{-1}=(\mathsf{J}\mathbb{V}^{\mathrm{th}})_{\rm
d}=(\imath\mathbb{k}^{\mathrm{th}})\oplus(-\imath\mathbb{k}^{\mathrm{th}})$,
or equivalently to

$\mathbb{d}(\mathbb{jc})^{-1}=2\mathbb{k}^{\mathrm{th}}=\coth\left(\frac{\hbar\beta\mathbb{w}}{2}\right)=\frac{\mathbb{1}+e^{-\hbar\beta\mathbb{w}}}{\mathbb{1}-e^{-\hbar\beta\mathbb{w}}}.$
(68)

However, this relation says nothing about the dependence on $\beta$, the
inverse temperature, of the matrices $\mathbb{d}$ and $\mathbb{jc}$ composed
by the eigenvalues of the matrices $\mathsf{J}\mathbb{D}$ and
$\mathsf{J}\mathbb{C}$, respectively. In the next section we will show that a
QDBC allows the determination of this dependence.

##  IV  GDSs satisfying a detailed balance condition

The notion of detailed balance is the principle governing the way thermal
equilibrium is attained by classical Markov processes Risken (1996). It has
several different quantum versions (see Carlen and Maas (2017) and references
therein) and in the context of  QDS for finite-dimensional systems, the one
due to Alicki Alicki (1976) stands out because it allows the extension of
time-reversal invariance of classical equilibrium to the quantum realm Carlen
and Maas (2017).

Inspired by the classical case in Markov processes, where the time-reversal
invariance of transition probabilities is related to a particular definition
of an inner product, Alicki’s definition for quantum detailed balance is based
on the $\hat{\bar{\rho}}$-Gelfand-Naimark-Segal ($\hat{\bar{\rho}}$-GNS) inner
product in finite dimension Hilbert spaces:

$\langle{\hat{A},\hat{B}}\rangle_{\text{GNS}}=\Tr(\hat{\bar{\rho}}\,\hat{A}^{\dagger}\hat{B}),$
(69)

where $\hat{\bar{\rho}}$ is a positive operator888The operator
$\hat{\bar{\rho}}$ in Alicki’s work Alicki (1976) is a full rank density
operator in finite dimensional systems. However, it is more convenient to
extend the definition of a $\hat{\bar{\rho}}$- GNS inner product for
unnormalized density operators $\hat{\bar{\rho}}$ and, in particular, to
unnormalized Gibbs states like $\hat{\bar{\sigma}}^{\mathrm{th}}$ in (41), see
Carlen and Maas (2017)., $\hat{\bar{\rho}}=\hat{\bar{\rho}}^{\dagger}>0$, and
the operators $\hat{A}$ and $\hat{B}$ belong to a finite-dimensional
C∗-algebra. Thus, Alicki’s  QDBC relies on the notion of self-adjointness with
respect to the $\hat{\bar{\rho}}^{\rm th}$- GNS inner product, where
$\hat{\bar{\rho}}^{\rm th}=e^{-\beta\hat{H}}$ is an unnormalized Gibbs state
with Hamiltonian $\hat{H}$. A superoperator $\Lambda$ is said to be self-
adjoint with respect to the $\hat{\bar{\rho}}$- GNS inner product if

$\displaystyle\expectationvalue{\Lambda[\hat{A}],\hat{B}}_{\text{GNS}}=\expectationvalue{\hat{A},\Lambda[\hat{B}]}_{\text{GNS}},$
(70)

for any operators $\hat{A}$ and $\hat{B}$ in the C∗-algebra.

The extension of Alicki’s approach for continuous variable systems relies on a
definition for the set of operators where the  GNS-inner-product is well
defined. In this regard, we consider operators acting on the separable Hilbert
space of $n$-bosonic modes, ${\cal H}=\otimes_{j=1}^{n}{\cal H}_{j}$, of
infinite dimension. In our case, the $\hat{\bar{\sigma}}^{\mathrm{th}}$- GNS
inner products are computed with $\hat{\bar{\sigma}}^{\mathrm{th}}$ in (41)
and $\hat{H}$ being the quadratic Hamiltonian in (42). Also, regardless of
whether the operators $\hat{A}$ and $\hat{B}$ are bounded or unbounded, having
continuous spectra or not, the domain of applicability of the
$\hat{\bar{\sigma}}^{\mathrm{th}}$- GNS inner product in (69) with
$\hat{\bar{\rho}}=\hat{\bar{\sigma}}^{\mathrm{th}}$ is over all operators such
that the trace on this formula is finite999The existence of the trace in (69)
can be checked, for example, using
$\langle{\hat{A},\hat{B}}\rangle_{\text{GNS}}=\int_{\mathds{R}^{2n}}\,d\boldsymbol{\xi}\,\bar{\sigma}^{\rm
th}(\boldsymbol{\xi})(A^{\dagger}B)(\boldsymbol{\xi})=\int_{\mathds{R}^{2n}}\,d\mathbf{x}\,\bar{\sigma}^{\rm
th}(\mathbf{x})(A^{\dagger}B)(\mathbf{x})$, where in the integrands we have
the Weyl and Wigner symbols of $\hat{\bar{\sigma}}^{\mathrm{th}}$ and
$\hat{A}^{\dagger}\hat{B}$, respectively. Notwithstanding, any other
representation that could be more convenient can be used..

For the definition of the  QDBC in the context of  GDSs, we recall the
notation of Sec.II, where $\bar{\Lambda}_{t}=e^{t\bar{\mathcal{L}}}$ with
$\bar{\mathcal{L}}$ in (4) represents the superoperator that generates a  GDS
in the Heisenberg picture and ${\Lambda}_{t}=e^{t\mathcal{L}}$ with
$\mathcal{L}$ in (1), the one in the Schrodinger picture.

###### Definition 1

Consider the  GDS $\\{\bar{\Lambda}_{t}=e^{t\bar{\mathcal{L}}}\\}_{t\geq 0}$
with the infinitesimal generator
$\bar{\mathcal{L}}=\bar{\mathcal{L}}_{\mathrm{U}}+\bar{\mathcal{L}}_{\mathrm{NU}}^{\mathrm{G}}$,
where $\bar{\mathcal{L}}_{\mathrm{U}}$ is defined in (5a) for the quadratic
Hamiltonian in (44) and $\bar{\mathcal{L}}_{\mathrm{NU}}^{\mathrm{G}}$ is
defined in (12). This  GDS satisfies the  QDBC with respect to
$\hat{\bar{\sigma}}^{\mathrm{th}}$ in (41) if $\bar{\Lambda}_{t}$ is self-
adjoint with respect to the $\hat{\bar{\sigma}}^{\mathrm{th}}$- GNS inner
product for all $t$. In such case, we say that $\bar{\Lambda}_{t}$ satisfies a
$\hat{\bar{\sigma}}^{\mathrm{th}}$-DBC.

The connection between the  QDBC and a steady state of the  GDS is in the
following theorem, which is proved in Appendix B.

###### Theorem 2

If a  GDS $\bar{\Lambda}_{t}$ (Heisenberg picture) satisfies a
$\hat{\bar{\sigma}}^{\mathrm{th}}$- DBC, then
$\hat{\bar{\sigma}}^{\mathrm{th}}$ is invariant under the  GDS $\Lambda_{t}$
(Schrödinger picture), i.e.,
$\Lambda_{t}[\hat{\bar{\sigma}}^{\mathrm{th}}]=\hat{\bar{\sigma}}^{\mathrm{th}}$
or $\mathcal{L}[\hat{\bar{\sigma}}^{\mathrm{th}}]=0$, equivalently.

In the theorem above, the statement
$\mathcal{L}[\hat{\bar{\sigma}}^{\mathrm{th}}]=0$ follows from
$d\Lambda_{t}[\hat{\bar{\sigma}}^{\mathrm{th}}]/dt=\mathcal{L}[\Lambda_{t}[\hat{\bar{\sigma}}^{\mathrm{th}}]]$.
So, if a  GDS satisfies a $\hat{\bar{\sigma}}^{\mathrm{th}}$- DBC, the quantum
state
$\hat{\sigma}^{\mathrm{th}}=\hat{\bar{\sigma}}^{\mathrm{th}}/\Tr(\hat{\bar{\sigma}}^{\mathrm{th}})$
is a stationary state of the evolution. Therefore, in order to attain thermal
equilibrium, it is enough that the superoperator $\bar{\Lambda}_{t}$ of a  GDS
satisfies a $\hat{\bar{\sigma}}^{\mathrm{th}}$- DBC for (the unnormalized
Gibbs state) $\hat{\bar{\sigma}}^{\mathrm{th}}$ defined in (41) with $\hat{H}$
in (42).

The equilibrium properties of a  GDS can be extracted from its relation with
the so called modular automorphism group Carlen and Maas (2017):

$\Xi_{t}[\hat{A}]=e^{\frac{\imath}{\hbar}\hat{H}t}\,\hat{A}\,e^{-\frac{\imath}{\hbar}\hat{H}t},$
(71)

with $\hat{H}$ in (42) and $t\in\mathds{C}$. Of particular relevance will be
the elements of the group given by the superoperator
$\Xi_{-\imath\hbar\beta}[\cdot]=(\hat{\bar{\sigma}}^{\mathrm{th}})^{-1}\,\cdot\,\hat{\bar{\sigma}}^{\mathrm{th}}$,
which existence is guaranteed for any finite value of $\beta$. The relation
between the $\hat{\bar{\sigma}}^{\mathrm{th}}$- DBC for an  GDS and the above
defined modular group is established in the following theorem:

###### Theorem 3

If a  GDS satisfies a $\hat{\bar{\sigma}}^{\mathrm{th}}$- DBC, then
$\bar{\Lambda}_{t}=e^{t\bar{\mathcal{L}}}$ and $\bar{\mathcal{L}}$ both
commute with $\Xi_{t}$ for all values of $t\in\mathds{C}$.

This theorem was proved in Carlen and Maas (2017) for  QDS in finite-
dimensional unital C∗-algebras and our demonstration for  GDSs follows almost
the same lines, see Appendix C.

Now, due to the commutation relation between the Hamiltonians (42) and (44),
$[\hat{H},\hat{H}_{\mathrm{eff}}]=0$, we have that
$\bar{\mathcal{L}}_{\mathrm{U}}$ commutes with $\Xi_{t}$. Therefore, the
commutation of
$\bar{\mathcal{L}}=\bar{\mathcal{L}}_{\mathrm{U}}+\bar{\mathcal{L}}_{\mathrm{NU}}^{\mathrm{G}}$
with the automorphism $\Xi_{t}$ in (71) is equivalent to the following
statement: $\bar{\mathcal{L}}_{\mathrm{NU}}^{\mathrm{G}}$ commutes with
$\Xi_{t}$. The generator $\bar{\mathcal{L}}_{\mathrm{NU}}^{\rm G}$ in (12) for
a  GDS is an explicit function of both diffusion and dissipation matrices,
$\mathbb{D}$ and $\mathbb{C}$ respectively. The properties of theses matrices
that stems from the fact that $\bar{\mathcal{L}}_{\mathrm{NU}}^{\mathrm{G}}$
commutes with $\Xi_{t}$ is settled by the following theorem:

###### Theorem 4

A  GDS satisfies a $\hat{\bar{\sigma}}^{\mathrm{th}}$- DBC if and only if the
diffusion and dissipation matrices, defined in Eqs.(10), are such that

$\displaystyle\mathbb{D}=\tilde{\mathsf{S}}_{t}\mathbb{D}\tilde{\mathsf{S}}_{t}^{\intercal}\,\,\,{\text{and}}\,\,\,\mathbb{C}=\tilde{\mathsf{S}}_{t}\mathbb{C}\tilde{\mathsf{S}}_{t}^{\intercal},$
(72)

for $\tilde{\mathsf{S}}_{t}$ ($t\in\mathds{R}$) in Eq.(55).

This means that both matrices are invariant under a congruence relation
through the symplectic matrix $\tilde{\mathsf{S}}_{t}$.

We begin the proof first noting that, for real values of $t$, the operator
$e^{\frac{\imath}{\hbar}\hat{H}t}$ in (71), with $\hat{H}$ in (42), belongs to
the metaplectic group Mp($2n$,$\mathds{R}$) of unitary operators and,
consequently, is associated with the symplectic matrix
$\tilde{\mathsf{S}}_{t}^{-1}=\tilde{\mathsf{S}}_{-t}$ defined by (55) Ozorio
de Almeida (1998); de Gosson (2006). So, in the Heisenberg picture, the action
of these operators on the vector (6) is described by

$\displaystyle\Xi_{-t}[\hat{\mathbf{x}}]=\tilde{\mathsf{S}}_{t}^{-1}\hat{\mathbf{x}}.$
(73)

Note that the above equation is equivalent to
$\Xi_{t}[\hat{\mathbf{x}}]=\tilde{\mathsf{S}}_{t}\hat{\mathbf{x}}$. Using
these actions, in Appendix D, we prove that

$\Xi_{-t}\left[\partialderivative{\hat{\mathbf{x}}^{\intercal}}[\cdot]\right]=\partialderivative{\hat{\mathbf{x}}^{\intercal}}[\Xi_{-t}[\cdot]]\,\tilde{\mathsf{S}}_{t},$
(74)

where $\Xi_{-t}=\Xi_{t}^{-1}$. This relation can be equivalently rewritten as
$\Xi_{-t}[\partial/\partial\hat{\mathbf{x}}[\cdot]]=\tilde{\mathsf{S}}_{t}^{\intercal}\,\partial/\partial\hat{\mathbf{x}}[\Xi_{-t}[\cdot]]$.
Now, using Eq.(74), we get

$\displaystyle\Xi_{-t}\left[\mathsf{J}\hat{\mathbf{x}}\partialderivative{\hat{\mathbf{x}}}\,[\Xi_{t}[\cdot]]\right]$
$\displaystyle=$
$\displaystyle\tilde{\mathsf{S}}_{t}^{\intercal}\,\mathsf{J}\hat{\mathbf{x}}\partialderivative{\hat{\mathbf{x}}}[\cdot]\,\tilde{\mathsf{S}}_{t},$
(75a)
$\displaystyle\Xi_{-t}\left[\partialderivative{\hat{\mathbf{x}}}\partialderivative{\hat{\mathbf{x}}^{\intercal}}[\Xi_{t}[\cdot]]\right]$
$\displaystyle=$
$\displaystyle\tilde{\mathsf{S}}_{t}^{\intercal}\partialderivative{\hat{\mathbf{x}}}\partialderivative{\hat{\mathbf{x}}^{\intercal}}[\cdot]\tilde{\mathsf{S}}_{t},$
(75b)

where in (75a) we used the symplectic condition
$\tilde{\mathsf{S}}_{t}^{-1}=-\mathsf{J}\tilde{\mathsf{S}}_{t}^{\intercal}\mathsf{J}$.
Finally, inserting Eqs.(75) in (12) we attain

$\displaystyle\Xi_{-t}[\bar{\mathcal{L}}_{\mathrm{NU}}^{\rm
G}[\Xi_{t}[\cdot]]]=$
$\displaystyle=\tr\\!\left(\tilde{\mathsf{S}}_{t}\mathbb{C}\tilde{\mathsf{S}}_{t}^{\intercal}\mathsf{J}\hat{\mathbf{x}}\partialderivative{\hat{\mathbf{x}}^{\intercal}}\,[\hat{O}_{t}]+\frac{1}{2}\tilde{\mathsf{S}}_{t}\mathbb{D}\tilde{\mathsf{S}}_{t}^{\intercal}\partialderivative{\hat{\mathbf{x}}}\partialderivative{\hat{\mathbf{x}}^{\intercal}}\hat{O}_{t}\right),$
(76)

which is equal to $\bar{\mathcal{L}}_{\mathrm{NU}}^{\rm G}[\cdot]$ in (12) iff
$\mathbb{D}$ and $\mathbb{C}$ satisfy Eqs.(72). In summary, all these prove
that the Eqs.(72) are equivalent to the statement that
$\bar{\mathcal{L}}_{\mathrm{NU}}^{\rm G}$ commutes with $\Xi_{t}$ for any real
value of $t$. However, due to
$\Xi_{t^{*}}=\left(\Xi_{t}\right)^{{t^{*}}/{t}}$, $\bar{\mathcal{L}}^{\rm
G}_{\mathrm{NU}}$ must also commute with $\Xi_{t^{*}}$ for any complex value
$t^{*}$. We finish the demonstration of Theorem 4 noting that the commutation
between $\bar{\mathcal{L}}_{\mathrm{NU}}^{\rm G}$ and $\Xi_{t}$ with
$t\in\mathds{C}$ is tantamount to say that a  GDS verifies a
$\hat{\bar{\sigma}}^{\mathrm{th}}$- DBC, according to Theorem 3.

As a consequence of the results in Theorem 4, the Lindblad operators in (8b)
are restricted to a particular structure, since the congruence relations in
(72) are extended to the decoherence matrix $\mathbb{\Gamma}$ in (11), due to
the definitions in (10), that is,
$\mathbb{\Gamma}=\tilde{\mathsf{S}}_{t}\mathbb{\Gamma}\tilde{\mathsf{S}}_{t}^{\intercal}$
for $\tilde{\mathsf{S}}_{t}$ ($t\in\mathds{R}$) in Eq.(55). Explicitly, the
continuous-variable version of Theorem 3 from Alicki’s work Alicki (1976),
which deals with  QDS in discrete Hilbert spaces, is a mere reformulation of
our Theorem 4:

###### Theorem 5

A  GDS satisfies a $\hat{\bar{\sigma}}^{\mathrm{th}}$- DBC iff the Lindblad
operators describing the  GDS are eigenoperators of the automorphism group
$\Xi_{-t}$, i.e.,

$\displaystyle\Xi_{-t}[\hat{L}_{j}]$ $\displaystyle=$ $\displaystyle
e^{\imath\omega_{j}t}\hat{L}_{j},$ (77a)
$\displaystyle\Xi_{-t}[\hat{L}_{n+j}]$ $\displaystyle=$ $\displaystyle
e^{-\imath\omega_{j}t}\hat{L}_{n+j}=e^{-\imath\omega_{j}t}e^{-\frac{1}{2}\hbar\beta\omega_{j}}\hat{L}_{j}^{\dagger},$
(77b)

with $j=1,\ldots,n$ and $\omega_{j}>0$ are the eigenfrequencies (symplectic
eigenvalues) of the Hessian matrix $\mathbb{B}$ which defines, through the
Hamiltonian Eq.(42), the thermal equilibrium state
$\hat{\sigma}^{\mathrm{th}}$ of the  GDS.

The proof for this theorem stands on Theorem 4 and on Eq.(73), and some
technical details are placed in Appendix E. In this appendix we prove that
Eqs.(72) are equivalent to write the decoherence matrix $\mathbb{\Gamma}$ in
the following characteristic form:

$\displaystyle\mathbb{\Gamma}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{n}\left(|s_{j}|^{2}\bar{\mathbf{l}}_{j}\bar{\mathbf{l}}_{j}^{\dagger}+|r_{j}|^{2}\bar{\mathbf{l}}^{*}_{j}(\bar{\mathbf{l}}^{*}_{j})^{\dagger}\right)$
(78) $\displaystyle=$
$\displaystyle(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{-1}(\mathbb{s}\oplus\mathbb{r})(\mathbb{s}^{*}\oplus\mathbb{r}^{*})(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{-\dagger},$

where $\\{\bar{\mathbf{l}}_{j}\\}_{j=1,\ldots,n}$ are the eigenvectors of
$\tilde{\mathsf{S}}_{t}$, i.e.,

$\tilde{\mathsf{S}}_{t}\,\bar{\mathbf{l}}_{j}=e^{\imath\omega_{j}t}\bar{\mathbf{l}}_{j}\quad(j=1,\ldots,n),\\\
$ (79)

the matrix $\mathbb{r}:=\text{diag}(r_{1},\ldots,r_{n})$ is the diagonal
matrix satisfying

$|\mathbb{r}|^{2}=\text{diag}(|r_{1}|^{2},\ldots,|r_{n}|^{2})=e^{-\frac{1}{2}\hbar\beta\mathbb{w}}\,|\mathbb{s}|^{2},$
(80)

see Eq.(124) in Appendix E; the diagonal matrix $\mathbb{w}$ is defined in
(51) and contains the symplectic spectrum of the Hessian matrix $\mathbb{B}$
of the Hamiltonian in (42); the matrix
$|\mathbb{s}|^{2}=\mbox{diag}(|s_{1}|^{2},\ldots,|s_{n}|^{2})$ is a real
diagonal matrix with a particular temperature dependence. Although Theorem 5
implies this dependence for any  GDS satisfying its conditions, we will keep
our track on the proof postponing the analysis of $|\mathbb{s}|^{2}$ to
Sec.VI, see Eq.(94).

Comparing the canonical form (78) with (11), the Lindblad operators
$\hat{L}_{k}=\mathbf{l}_{k}^{\intercal}\mathsf{J}\hat{\mathbf{x}}$ in (8b),
with $k=1,\ldots,K=2n$, correspond to the vectors

$\displaystyle\mathbf{l}_{j}=s_{j}\bar{\mathbf{l}}_{j}\,,\quad\mathbf{l}_{n+j}$
$\displaystyle=$ $\displaystyle
r_{j}\bar{\mathbf{l}}^{\ast}_{j}\quad(j=1,\ldots,n).$ (81)

From (79),
$\tilde{\mathsf{S}}_{t}\mathbf{l}_{n+j}=\tilde{\mathsf{S}}_{t}\bar{\mathbf{l}}_{j}^{\ast}=e^{-\imath\omega_{j}t}r_{j}\bar{\mathbf{l}}_{j}^{\ast}$,
consequently
$\Xi_{-t}[\hat{L}_{j}]=\mathbf{l}_{j}^{\intercal}\tilde{\mathsf{S}}_{t}^{\intercal}\mathsf{J}\hat{\mathbf{x}}=s_{j}\bar{\mathbf{l}}_{j}^{\intercal}\tilde{\mathsf{S}}_{t}^{\intercal}\mathsf{J}\hat{\mathbf{x}}=e^{\imath\omega_{j}t}s_{j}\bar{\mathbf{l}}_{j}^{\intercal}\mathsf{J}\hat{\mathbf{x}}=e^{\imath\omega_{j}t}\hat{L}_{j}$
and
$\Xi_{-t}[\hat{L}_{n+j}]=r_{j}\bar{\mathbf{l}}^{\dagger}_{j}\tilde{\mathsf{S}}_{t}^{\intercal}\mathsf{J}\hat{\mathbf{x}}=e^{\imath\omega_{j}t}\,r_{j}\,\bar{\mathbf{l}}^{\dagger}_{j}\mathsf{J}\hat{\mathbf{x}}=e^{-\imath\omega_{j}t}e^{-\frac{1}{2}\hbar\beta\omega_{j}}\hat{L}^{\dagger}_{j}$,
where we used (73) and (80) and the symplectic condition
$\mathsf{J}\tilde{\mathsf{S}}_{t}^{-1}=\tilde{\mathsf{S}}_{t}^{\intercal}\mathsf{J}$.
With all these, we finish the proof of Theorem 5.

Here, two important observations are in order. First, the decoherence matrix
(78) could be written as
$\mathbb{\Gamma}=\mathbb{\Upsilon}\mathbb{\Upsilon}^{\dagger}$, where

$\mathbb{\Upsilon}=(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{-1}(\mathbb{s}\oplus\mathbb{s}e^{-\frac{1}{2}\hbar\beta\mathbb{w}})\mathbb{W}$
(82)

and $\mathbb{W}\in\mathds{C}^{2n\times 2n}$ is an arbitrary unitary matrix.
Therefore, alternatively we can use the columns vectors
$\mathbf{l}^{\prime}_{j}$ of the matrix $\mathbb{\Upsilon}$ to define new
Lindblad operators
$\hat{L}_{k}^{\prime}=(\mathbf{l}^{\prime}_{k})^{\intercal}\mathsf{J}\hat{\mathbf{x}}$
with $k=1,\ldots,2n$. Notwithstanding, it is straightforward to check that
transformation (82) corresponds to

$\hat{L}_{k}\longrightarrow\hat{L}^{\prime}_{k}=\sum^{2n}_{j=1}\mathbb{W}_{kj}\,\hat{L}_{j}$
(83)

and the arbitrariness introduced by $\mathbb{W}$ in (82) is equivalent to a
well known symmetry of any  QDS (see e.g., (Breuer et al., 2002, Sec.3.2.2)):
the  LME (1) is invariant under the unitary transformation in (83) of the
Lindblad operators. Therefore, this symmetry also holds for a  GDS. Thus, the
semigroup dynamic associated to the new set of Lindblad operators,
$\hat{L}_{k}^{\prime}$, is exactly the same as the one generated by the old
set, i.e., $\hat{L}_{k}$ in (8b) with $\bar{\mathbf{l}}_{k}$ in (81).
Secondly, Theorem 5 shows that it is enough to consider only $n$ Lindblad
operators to describe an  GDS which satisfies a
$\hat{\bar{\sigma}}^{\mathrm{th}}$- DBC, although we are dealing with infinite
dimensional quantum systems. In this sense,  GDS are like  QDS in finite-
dimensional systems.

Properties of the environment, characterized by the matrices $\mathbb{D}$ and
$\mathbb{C}$, can be extracted from the canonical form of the decoherence
matrix $\mathbb{\Gamma}$ in (78). To this end, we conveniently rewrite
$\mathbb{\Gamma}$ in (78) as

$\mathbb{\Gamma}=\frac{\mathbb{D}}{\hbar}+\imath\mathbb{C}=({\mathsf{S}}^{\mathrm{th}})^{-1}\left(\mathbb{\Lambda}_{r}-\imath\mathsf{J}\mathbb{\Lambda}_{i}\right)({\mathsf{S}}^{\mathrm{th}})^{-\intercal},$
(84)

where we define the diagonal positive defined matrices

$\displaystyle\mathbb{\Lambda}_{r}$ $\displaystyle=$
$\displaystyle\frac{|\mathbb{s}|^{2}+|\mathbb{r}|^{2}}{2}\oplus\frac{|\mathbb{s}|^{2}+|\mathbb{r}|^{2}}{2}$
(85a) $\displaystyle=$
$\displaystyle\frac{|\mathbb{s}|^{2}}{2}(\mathbb{1}+e^{-\hbar\beta\mathbb{w}})\oplus\frac{|\mathbb{s}|^{2}}{2}(\mathbb{1}+e^{-\hbar\beta\mathbb{w}}),$
$\displaystyle\text{and}\;\,\mathbb{\Lambda}_{i}$ $\displaystyle=$
$\displaystyle\frac{|\mathbb{s}|^{2}-|\mathbb{r}|^{2}}{2}\oplus\frac{|\mathbb{s}|^{2}-|\mathbb{r}|^{2}}{2}$
(85b) $\displaystyle=$
$\displaystyle\frac{|\mathbb{s}|^{2}}{2}(\mathbb{1}-e^{-\hbar\beta\mathbb{w}})\oplus\frac{|\mathbb{s}|^{2}}{2}(\mathbb{1}-e^{-\hbar\beta\mathbb{w}}),$

using the notation
$|\mathbb{s}|^{2}=\text{diag}(|s_{1}|^{2},\ldots,|s_{n}|^{2})$,
$|\mathbb{r}|^{2}=\text{diag}(|r_{1}|^{2},\ldots,|r_{n}|^{2})$ and the matrix
relation in (80). Therefore, the diffusion and dissipation matrices are

$\displaystyle\mathbb{D}$ $\displaystyle=$
$\displaystyle({\mathsf{S}}^{\mathrm{th}})^{-1}\hbar\mathbb{\Lambda}_{r}({\mathsf{S}}^{\mathrm{th}})^{-\intercal},$
(86a) $\displaystyle\mathbb{C}$ $\displaystyle=$
$\displaystyle({\mathsf{S}}^{\mathrm{th}})^{-1}\mathsf{J}^{\intercal}\mathbb{\Lambda}_{i}({\mathsf{S}}^{\mathrm{th}})^{-\intercal}.$
(86b)

Note that each pair of matrices $\mathbb{D}$ and $\mathbb{C}$ uniquely
determine a  GDS. Since these only depend on the real matrix
$|\mathbb{s}|^{2}$, it is enough to choose a real matrix
$|\mathbb{s}|=\text{diag}(|s_{1}|,\ldots,|s_{n}|)$ instead of a complex matrix
$\mathbb{s}$ in (78) with $\mathbb{r}$ in (80). In this way the relations in
(81) change to

$\displaystyle\mathbf{l}_{j}=|s_{j}|\bar{\mathbf{l}}_{j},\,\,\mathbf{l}_{n+j}=|s_{j}|\,e^{-\frac{1}{2}\hbar\beta\omega_{j}}\bar{\mathbf{l}}^{\ast}_{j}\,\,\,\,(j=1,\ldots,n),$
(87)

which gives the expressions for the coefficients of the Lindblad operators in
(8b).

The structure of $\mathbb{D}$ and $\mathbb{C}$ in (86) is determined by the
matrices ${\mathsf{S}}^{\mathrm{th}}$, $\mathbb{w}$, and $|\mathbb{s}|$. It is
worth to note that, according to Eq.(52), the matrices
${\mathsf{S}}^{\mathrm{th}}$ and $\mathbb{w}$ can be extracted from the
ordinary diagonalization of the Hamiltonian matrix $\mathsf{J}\mathbb{B}$,
where $\mathbb{B}$ is the Hessian matrix of the Hamiltonian $\hat{H}$ in (42),
that defines the thermal equilibrium state $\hat{\sigma}^{\mathrm{th}}$ in
(41). However, we will see in Sec.VI that $|\mathbb{s}|$ depends on the
coupling constants of the system and the environment. Thus, each matrix
$|\mathbb{s}|$, corresponding to different coupling constants, defines one
different pair of diffusion and dissipation matrices through Eqs.(86). The
dynamics of the  GDSs associated to these matrices is different because each
one corresponds to different Lindblad operators, which are not associated with
the symmetry in (83). Nonetheless, all these  GDSs have the same thermal
equilibrium state $\hat{\sigma}^{\mathrm{th}}$. This is checked in Appendix E
through the symplectic diagonalization of $\mathbb{V}^{\mathrm{th}}$ in (121),
where the matrix
$\frac{1}{2}\mathbb{\Lambda}_{r}\mathbb{\Lambda}_{i}^{-1}=\mathbb{k}^{\mathrm{th}}\oplus\mathbb{k}^{\mathrm{th}}$
is the direct sum of the symplectic spectrum of $\mathbb{V}^{\mathrm{th}}$
that does not depend on $|\mathbb{s}|$.

##  V The master equation of a GDS satisfying a
$\hat{\bar{\sigma}}^{\mathrm{th}}$-DBC.

Here, we prove that the master equation of a  GDS that satisfies a
$\hat{\bar{\sigma}}^{\mathrm{th}}$- DBC, i.e., satisfying all theorems in last
section, has the form of the Quantum Optical Master Equation ( QOME), see
Eq.(5.4.14) in Gardiner and Zoller (2000):

$\displaystyle\derivative{\hat{\rho}_{t}}{t}$ $\displaystyle=$
$\displaystyle-\frac{\imath}{\hbar}[\hat{H}_{\mathrm{OME}},\hat{\rho}_{t}]$
(88)
$\displaystyle+\sum_{k}\frac{\bar{\gamma}_{k}}{2}(\bar{n}_{k}+1)\left(2\hat{X}^{-}_{k}\hat{\rho}_{t}\hat{X}^{+}_{k}-\\{\hat{X}^{+}_{k}\hat{X}^{-}_{k},\hat{\rho}_{t}\\}\right)$
$\displaystyle+\sum_{k}\frac{\bar{\gamma}_{k}}{2}\bar{n}_{k}\left(2\hat{X}^{+}_{k}\hat{\rho}_{t}\hat{X}^{-}_{k}-\\{\hat{X}^{-}_{k}\hat{X}^{+}_{k},\hat{\rho}_{t}\\}\right),$

where $\bar{\gamma}_{k}$ are the coupling constants between the bath and the
system, $\bar{n}_{k}=(e^{\hbar\beta\tilde{\omega}_{k}}-1)^{-1}$ is the Planck
factor with $\tilde{\omega}_{k}>0$, and the operators $\hat{X}^{\pm}_{k}$ are
eigenoperators of $\hat{H}_{\mathrm{OME}}$, i.e.,

$[\hat{H}_{\mathrm{OME}},\hat{X}^{\pm}_{k}]=\pm\hbar\tilde{\omega}_{k}\hat{X}^{\pm}_{k},$
(89)

so $(\hat{X}^{-}_{k})^{\dagger}=\hat{X}^{+}_{k}$.

In order to rewrite the master equation of a  GDS satisfying a
$\hat{\bar{\sigma}}^{\mathrm{th}}$- DBC, first we rewrite the Lindblad
operators in (8b) as

$\displaystyle\hat{L}_{j}$ $\displaystyle=$ $\displaystyle|s_{j}|\hat{X}^{-},$
(90a) $\displaystyle\hat{L}_{n+j}$ $\displaystyle=$ $\displaystyle
e^{-\frac{1}{2}\hbar\beta\omega_{j}}\hat{L}^{\dagger}_{j}=|s_{j}|e^{-\frac{1}{2}\hbar\beta\omega_{j}}\hat{X}^{+}_{j}$
(90b)

with

$\hat{X}^{-}_{j}=\bar{\mathbf{l}}_{j}^{\intercal}\mathsf{J}\hat{\mathbf{x}}\quad\text{and}\quad\hat{X}^{+}_{j}=(\hat{X}^{-})^{\dagger}=(\bar{\mathbf{l}}^{\ast}_{j})^{\intercal}\mathsf{J}\hat{\mathbf{x}},$
(91)

when employing Eq.(87) for $j=1,...,n$. Therefore, using these expressions for
the Lindblad operators, the master equation in (1) becomes

$\displaystyle\\!\\!\\!\\!\derivative{\hat{\rho}_{t}}{t}$ $\displaystyle=$
$\displaystyle-\frac{\imath}{\hbar}[\hat{H}_{\mathrm{eff}},\hat{\rho}_{t}]$
$\displaystyle+$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\sum_{j=1}^{n}\frac{|s_{j}|^{2}}{2\hbar}(2\hat{X}^{-}_{j}\hat{\rho}_{t}\hat{X}^{+}_{j}-\\{\hat{X}^{+}_{j}\hat{X}^{-}_{j},\hat{\rho}_{t}\\})$
(92a) $\displaystyle+$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\sum_{j=1}^{n}\frac{|s_{j}|^{2}}{2\hbar}e^{-\hbar\beta\omega_{j}}\\!(2\hat{X}^{+}_{j}\hat{\rho}_{t}\hat{X}^{-}_{j}-\\{\hat{X}^{-}_{j}\hat{X}^{+}_{j},\hat{\rho}_{t}\\}).$
(92b)

In Appendix F, we prove that $\hat{X}_{j}^{\pm}$ are eigenoperators of
$\hat{H}$ in (42), i.e.,

$[\hat{H},\hat{X}_{j}^{\pm}]=\pm\hbar\omega_{j}\hat{X}_{j}^{\pm}\,\,\,\,\,\,(j=1,\ldots,n),$
(93)

while in Appendix G there is a proof for

$|s_{j}|^{2}=\hbar\bar{\gamma}_{j}\,\bar{n}_{j}\,e^{\hbar\beta\omega_{j}}=\hbar\bar{\gamma}_{j}\,(\bar{n}_{j}+1),$
(94)

where now the Planck factor is

$\bar{n}_{k}=(e^{\hbar\beta\omega_{k}}-1)^{-1}.$ (95)

With all these, we show that the master equation (92) has the structure of a
QOME like the one in (88).

Using Eq.(51) and writing $\hat{U}_{{\mathsf{S}}^{\mathrm{th}}}$ for the
metaplectic operator associated with the symplectic matrix
${\mathsf{S}}^{\mathrm{th}}$, i.e.,
$\hat{U}_{{\mathsf{S}}^{\mathrm{th}}}^{\dagger}\hat{\mathbf{x}}\hat{U}_{{\mathsf{S}}^{\mathrm{th}}}={\mathsf{S}}^{\mathrm{th}}\hat{\mathbf{x}}$,
one can find

$\hat{H}=\frac{1}{2}\hat{\mathbf{x}}^{\intercal}\mathbb{B}\hat{\mathbf{x}}=\hat{U}_{{\mathsf{S}}^{\mathrm{th}}}^{\dagger}\hat{H}_{{\rm
ho}}\hat{U}_{({\mathsf{S}}^{\mathrm{th}})},$ (96)

where

$\hat{H}_{{\rm
ho}}:=\frac{1}{2}\hat{\mathbf{x}}^{\intercal}(\mathbb{w}\oplus\mathbb{w})\hat{\mathbf{x}}=\sum_{j=1}^{n}\frac{\omega_{j}}{2}(\hat{p}_{j}^{2}+\hat{q}_{j}^{2})$

is the Hamiltonian of a multimode harmonic oscillator. The eigenequation for
$\hat{H}_{\rm ho}$ is

$\hat{H}_{\rm ho}\ket{\bf n}=E_{\bf n}\ket{\bf n},\,\,\,E_{\bf
n}=\sum_{l=1}^{n}\hbar\omega_{l}\left(n_{l}+\frac{1}{2}\right),$ (97)

where $\ket{\bf n}=\ket{n_{1}}\otimes\ldots\otimes\ket{n_{n}}$, ${\bf
n}=(n_{1},\ldots,n_{n})$ with $n_{l}=0,\ldots,+\infty$. Therefore, due to the
similarity relation between $\hat{H}_{{\rm ho}}$ and $\hat{H}$, the
eigenequation of $\hat{H}$ is $\hat{H}\ket{\phi_{\bf n}}=E_{\bf
n}\ket{\phi_{\bf n}}$, with

$\ket{\phi_{\bf n}}=\hat{U}_{{\mathsf{S}}^{\mathrm{th}}}^{\dagger}\ket{{\bf
n}}.$ (98)

Now, from (93), it is straightforward to verify that $\hat{X}^{+}_{j}$ and
$\hat{X}^{-}_{j}$ are the creation and annihilation operators of a quantum
$\hbar\omega_{j}$ associated to the states $\ket{\phi_{\bf n}}$; that is,
$\hat{X}^{\pm}_{j}\ket{\phi_{\bf n}}=\alpha_{{\bf n}_{j}^{\pm}}|\phi_{{\bf
n}_{j}^{\pm}}\rangle$ is an eigenvector of $\hat{H}$ with energy $E_{{\bf
n}_{j}^{\pm}}$ such that ${\bf n}_{j}^{\pm}=(n_{1},\ldots,n_{j}\pm
1,\ldots,n_{n})$. Since $\hat{X}_{j}^{\pm}$ satisfies the commutation relation
(93), the operator $\hat{X}^{+}_{j}\hat{X}^{-}_{j}$ is the number operator in
the eigenbasis $\\{\ket{\phi_{\bf n}}\\}$, i.e.,
$\hat{X}^{+}_{j}\hat{X}^{-}_{j}\ket{\phi_{\bf n}}=n_{j}\ket{\phi_{\bf n}}$,
therefore $\bar{\alpha}_{{\bf n}_{j}^{+}}=\sqrt{n_{j}+1}$ and
$\bar{\alpha}_{{\bf n}_{j}^{-}}=\sqrt{n_{j}}$ 101010This can be checked
through a lengthy, but not difficult, calculation using the definition of
$\hat{X}_{j}^{\pm}$ in (91) and $\ket{\phi_{\bf n}}$ in (98)..

It is worth to note that we consider the same mode structure for $\hat{H}_{\rm
ho}$ and $\hat{H}$, what changes is the nature of the stationary states, i.e.,
while $\ket{\bf n}$ are separable states, $\ket{\phi_{\bf n}}$ could be
entangled with respect to the considered mode structure. In an analogous way,
we can rewrite the thermal equilibrium state in (41) as

$\displaystyle\hat{\sigma}^{\mathrm{th}}$ $\displaystyle=$
$\displaystyle\frac{1}{{\cal
Z}^{\mathrm{th}}}e^{-\frac{\beta}{2}\hat{\mathbf{x}}^{\intercal}({\mathsf{S}}^{\mathrm{th}})^{\intercal}(\mathbb{w}\oplus\mathbb{w}){\mathsf{S}}^{\mathrm{th}}\hat{\mathbf{x}}}$
(99) $\displaystyle=$
$\displaystyle\hat{U}_{{\mathsf{S}}^{\mathrm{th}}}^{\dagger}\hat{\sigma}^{\mathrm{th}}_{\rm
ho}\hat{U}_{{\mathsf{S}}^{\mathrm{th}}},$

where we recognize in

$\hat{\sigma}^{\mathrm{th}}_{\rm ho}=\frac{1}{{\cal
Z}^{\mathrm{th}}}e^{-\frac{\beta}{2}\hat{\mathbf{x}}^{\intercal}(\mathbb{w}\oplus\mathbb{w})\hat{\mathbf{x}}}$
(100)

the Gibbs’s state of the multimode harmonic oscillator with Hamiltonian
$\hat{H}_{\rm ho}$ in (96). The state $\hat{\sigma}^{\mathrm{th}}_{\rm ho}$ is
manifestly separable, so the possible entanglement of the state
$\hat{\sigma}^{\mathrm{th}}$ is due to the action of the unitary operation
$\hat{U}_{({\mathsf{S}}^{\mathrm{th}})}$, where ${\mathsf{S}}^{\mathrm{th}}$
is the matrix of the symplectic diagonalization of the Hessian matrix
$\mathbb{B}$, as shown in Eq.(51).

Noteworthy that usually in the derivation of the  QOME, it is assumed that the
commutation relation (89) is valid Gardiner and Zoller (2000). Here, the
quantum detailed balance condition for a  GDS shows that commutation relation
(93) must be valid instead. When ${\mathsf{S}}^{\mathrm{th}}$ is the identity
and $\hat{H}_{\mathrm{eff}}=\hat{H}=\hat{H}_{\rm ho}$, we have that
$\hat{X}^{-}_{j}=\hat{a}_{j}$ and $\hat{X}^{+}_{j}=\hat{a}_{j}^{\dagger}$ are
the usual creation and annihilation operators of the multimode quantum
harmonic oscillator. Then, in this case the  QOME (92) is the well known
master equation of a multimode damped harmonic oscillator Gardiner and Zoller
(2004).

##  VI Additional properties of GDSs satisfying a
$\hat{\bar{\sigma}}^{\mathrm{th}}$-DBC

###  VI.1 The temperature dependence of $\mathbb{D}$ and $\mathbb{C}$

The final form for the diffusion and dissipation matrices of GDSs that satisfy
a $\hat{\bar{\sigma}}^{\mathrm{th}}$- DBC is obtained replacing (94) into
(86), so

$\displaystyle\mathbb{D}$ $\displaystyle=$
$\displaystyle\hbar^{2}({\mathsf{S}}^{\mathrm{th}})^{-1}(\bar{\mathbb{g}}\mathbb{k}^{\mathrm{th}}\oplus\bar{\mathbb{g}}\mathbb{k}^{\mathrm{th}})\;({\mathsf{S}}^{\mathrm{th}})^{-\intercal},$
(101a) $\displaystyle\mathbb{C}$ $\displaystyle=$
$\displaystyle\frac{\hbar}{2}({\mathsf{S}}^{\mathrm{th}})^{-1}\;\mathsf{J}^{\intercal}(\bar{\mathbb{g}}\oplus\bar{\mathbb{g}})({\mathsf{S}}^{\mathrm{th}})^{-\intercal},$
(101b)

where $\mathbb{k}^{\mathrm{th}}$ is given in (53) and $\bar{\mathbb{g}}$ is
the diagonal matrix containing the coupling constants, i.e.,

$\bar{\mathbb{g}}=\text{diag}(\bar{\gamma}_{1},\ldots,\bar{\gamma}_{n}).$
(102)

It is worth to remember here that ${\mathsf{S}}^{\mathrm{th}}$ and
$\mathbb{k}^{\mathrm{th}}$ are ultimately determined by $\mathbb{B}$, the
Hessian matrix of $\hat{H}$ that defines the thermal equilibrium state.
However, the coupling constants can take different values, thus defining
different pair of matrices, $\mathbb{D}$ and $\mathbb{C}$, all having the same
equilibrium state. These are the multiples solutions of Eq.(62).

Note that here we consider that the coupling constants do not depend on
temperature. This is consistent with the fact that the Fokker-Planck equation
(34) for Gaussian states corresponds to an Ornstein–Uhlenbeck process Risken
(1996), where the drift term
$\partialderivative{\mathsf{J}\mathbb{C}\mathbf{x}}{\mathbf{x}^{\intercal}}=\tr(\mathsf{J}\mathbb{C})$
corresponds to a drift force $-\mathsf{J}\mathbb{C}\mathbf{x}$ which does not
depend on the temperature Toscano et al. (2021). From Eq.(66), we can write
$\mathbb{D}=\hbar({\mathsf{S}}^{\mathrm{th}})^{-1}(\mathbb{d}\oplus\mathbb{d})({\mathsf{S}}^{\mathrm{th}})^{-\intercal}$,
where we use that ${\mathsf{S}}^{\mathrm{th}}$ is symplectic and that
$\mathbb{Q}^{\dagger}\,(\imath\mathbb{d}\oplus(-\imath\mathbb{d}))\,\mathbb{Q}=\hbar\mathsf{J}(\mathbb{d}\oplus\mathbb{d})$.
Analogously, from Eq.(67), we can write
$\mathbb{C}=({\mathsf{S}}^{\mathrm{th}})^{-1}\mathsf{J}^{\intercal}(\mathbb{jc}\oplus\mathbb{jc})({\mathsf{S}}^{\mathrm{th}})^{-\intercal}$,
where
$\mathbb{Q}^{\dagger}(\mathbb{jc}\oplus\mathbb{jc})\mathbb{Q}=\mathbb{jc}\oplus\mathbb{jc}$.
Comparing with Eqs.(101), we have
$\mathbb{d}=\hbar\bar{\mathbb{g}}\,\mathbb{k}^{\mathrm{th}}$ and
$\mathbb{jc}=\frac{\hbar}{2}\bar{\mathbb{g}}$, so
$\mathbb{d}=2\mathbb{jc}\,\mathbb{k}^{\mathrm{th}}$. Therefore, the diffusion
and dissipation matrices of a  GDS satisfying a  QDBC verify the condition
(68); this condition guarantees that the  GDS has a thermal equilibrium state.

For the validity of the reciprocal implication (every  GDS that attains
thermal equilibrium must satisfy a  QDBC), we need to allow an arbitrary
temperature dependence for the coupling constants in the matrix
$\bar{\mathbb{g}}$. This arbitrariness is clear taking into account (68),
where $\mathbb{d}=2\mathbb{jc}\,\mathbb{k}^{\mathrm{th}}$; now $\mathbb{jc}$
can depend on temperature still matching (101b) with
$\mathbb{jc}=\frac{\hbar}{2}\bar{\mathbb{g}}$ for an arbitrary dependence of
$\bar{\mathbb{g}}$ on temperature.

###  VI.2 The high and low temperature limits

The high temperature limit is obtained when considering
$\hbar\beta\|\mathbb{w}\|\ll 1$ in (101a), where
$\|\mathbb{w}\|={\max}\\{\omega_{1},...,\omega_{n}\\}$, so one can write
$\mathbb{k}^{\mathrm{th}}\approx(\hbar\beta\mathbb{w})^{-1}+\frac{1}{12}\hbar\beta\mathbb{w}+{\cal
O}(\hbar\beta\mathbb{w})^{3}$ and neglecting higher order terms, we write

$\mathbb{D}\approx\hbar\beta^{-1}({\mathsf{S}}^{\mathrm{th}})^{-1}(\bar{\mathbb{g}}\mathbb{w}^{-1}\oplus\bar{\mathbb{g}}\mathbb{w}^{-1})({\mathsf{S}}^{\mathrm{th}})^{-\intercal}.$
(103)

Inserting this into Eq.(62) with $\mathbb{C}$ in (101b), using the symplectic
condition for ${\mathsf{S}}^{\mathrm{th}}$, and the symplectic diagonalization
in (51), one obtains

$\mathbb{V}\approx({\hbar\beta}\mathbb{B})^{-1}.$ (104)

This limit is also called classical limit and the Wigner function in (37) for
this covariance matrix is the classical Boltzmann factor $\exp[-\beta H]$ for
the classical version (Wigner symbol) of Hamiltonian (42) Nicacio (2021b).

The low temperature limit corresponds to $\hbar\beta\|\mathbb{w}\|\gg 1$ and
in this case $\mathbb{k}^{\mathrm{th}}=\frac{1}{2}\mathbb{1}+{\cal
O}(e^{-\hbar\beta\mathbb{w}})$. From (101),

$\mathbb{D}\approx\hbar\frac{1}{2}({\mathsf{S}}^{\mathrm{th}})^{-1}({\mathsf{S}}^{\mathrm{th}})^{-\intercal}\mathsf{J}\mathbb{C}.$
(105)

Consequently,

$\mathbb{V}^{\mathrm{th}}\approx\frac{1}{2}({\mathsf{S}}^{\mathrm{th}})^{-1}({\mathsf{S}}^{\mathrm{th}})^{-\intercal},$
(106)

and the thermal equilibrium state corresponds to a pure Gaussian state,
$\hat{\sigma}^{\mathrm{th}}\approx\outerproduct{\phi_{{\bf n}={\bf
0}}}{\phi_{{\bf n}={\bf 0}}}$, where $\ket{\phi_{{\bf n}={\bf 0}}}$ is the
ground state of the Hamiltonian $\hat{H}$ in (42).

###  VI.3 The pure diffusive regime

This regime corresponds to the limit $\bar{n}_{j}\rightarrow+\infty$ (i.e.,
$\beta\rightarrow 0$) together with $\bar{\gamma}_{j}\rightarrow 0$, such that
$\bar{n}_{j}\bar{\gamma}_{j}=\bar{c}_{j}$ are constant values for
$j=1,\ldots,n$. In such limit, the coeficients in (92a) become equal to the
ones in (92b) for the master equation of the  GDS, see also Eq.(88). For the
difussion and dissipation matrices, respectively, we have

$\displaystyle\mathbb{D}$ $\displaystyle=$
$\displaystyle\hbar^{2}({\mathsf{S}}^{\mathrm{th}})^{-1}\,\bar{\mathbb{c}}\oplus\bar{\mathbb{c}}\,({\mathsf{S}}^{\mathrm{th}})^{-\intercal},$
(107a) $\displaystyle\mathbb{C}$ $\displaystyle=$ $\displaystyle 0,$ (107b)

where $\bar{\mathbb{c}}=\text{diag}(\bar{c}_{1},\ldots,\bar{c}_{n})$. In this
case the  GDS has no thermal equilibrium solution because $\mathbb{D}$ and
$\mathbb{C}$ do not satisfy Eq.(62).

It is worth to note that the master equation in (92) has no pure dissipative
regime, i.e., $\mathbb{D}=0$, since Eq. (13) is violated. However, if $\hbar$
represents an effective Planck constant in the master equation of the  GDS in
(92), the semiclassical limit $\hbar\ll 1$ and the low temperature condition
$\hbar\beta\|\mathbb{w}\|\gg 1$ guarantee the validity of (105), thus the
contributions from $\mathbb{D}$ can be effectively neglected when compared to
those coming from $\mathbb{C}$; in this case the evolution is thus dominated
by dissipation Toscano et al. (2005).

##  VII The role of $\hat{H}_{\rm eff}$ in thermalization

So far, we completely answer the question of which kind of environments,
characterized by the diffusion and dissipation matrices $\mathbb{D}$ and
$\mathbb{C}$, allows the existence of a given thermal equilibrium state in a
GDS. This state is given by a Gaussian Gibbs state,
$\hat{\sigma}^{\mathrm{th}}=e^{-\beta\hat{H}}/{\cal Z}^{\mathrm{th}}$,
corresponding to a quadratic Hamiltonian $\hat{H}$ of positive-elliptic type,
which is completely characterized by its covariance matrix
$\mathbb{V}^{\mathrm{th}}$ such that $\mathsf{J}\mathbb{V}^{\mathrm{th}}$
belongs to the commuting set
$\\{\mathsf{J}\mathbb{V}^{\mathrm{th}},\mathsf{J}\mathbb{D},\mathsf{J}\mathbb{C}\\}$.

However, it is interesting to analyze this result from the perspective of a
fixed  GDS that governs the system to thermal equilibrium. In this case, the
commuting matrices $\mathsf{J}\mathbb{D}$ and $\mathsf{J}\mathbb{C}$ determine
$\mathsf{J}\mathbb{B}$ through the relation (62) and (47). Then, it is clear
that the effective action of the environment, through the non-unitary part of
the evolution in (9), is to confine the system in phase space (with
coordinates $\mathbf{x}$) around the the equilibrium point
$\expectationvalue{\hat{\mathbf{x}}}^{\mathrm{th}}=0$. The energetic balance
between the environment and the system is described by the quantum positive-
elliptic Hamiltonian $\hat{H}$, which is dynamically stable (with spectrum
necessarily discrete with bounded eigenstates). If we admit that every  GDS
that leads to thermal equilibrium satisfies a  QDBC, the form of the master
equation in (92) together with (93) and (94), helps to understand in detail
the energetic balance of the thermal equilibrium process. Indeed, this process
is determined by the action of the eigenoperators $\hat{X}_{j}^{\pm}$ of
$\hat{H}$, which promote transitions between the energy levels of this
Hamiltonian as a consequence of the interaction with the environment. This
contrasts with the usual derivation Gardiner and Zoller (2000) of the  QOME in
(88), where $\\{\hat{X}_{j}^{\pm}\\}$ is considered as a basis of
eigenoperators of $\hat{H}_{\mathrm{OME}}$, which expands the algebra of
operators acting on the system Gardiner and Zoller (2000).

The induced environmental confinement process of the system in thermal
equilibrium is not influenced by the unitary dynamics generated by
$\hat{H}_{\mathrm{eff}}$, since the commuting set of matrices
$\\{\mathsf{J}\mathbb{V}^{\mathrm{th}},\mathsf{J}\mathbb{D},\mathsf{J}\mathbb{C}\\}$,
which ultimately determines the thermal equilibrium through (62), do not
depend on the Hessian matrix $\mathbb{B}^{\prime}$ of
$\hat{H}_{\mathrm{eff}}$. However, there are some of these Hamiltonians that
allows the process of thermalization to occur, these are determined by the
condition $[\hat{H},\hat{H}_{\mathrm{eff}}]=0$. Using Eq.(43) and the results
in Appendix H, all these Hamiltonians are written as

$\hat{H}_{\mathrm{eff}}=\frac{1}{2}\hat{\mathbf{x}}^{\intercal}\mathbb{B}^{\prime}\hat{\mathbf{x}}=\hat{U}_{{\mathsf{S}}^{\mathrm{th}}}^{\dagger}\hat{H}^{\prime}_{{\rm
ho}}\hat{U}_{({\mathsf{S}}^{\mathrm{th}})},$ (108)

where we employed the same reasoning performed in Eq.(96) with
$\mathbb{B}^{\prime}=({\mathsf{S}}^{\mathrm{th}})^{\intercal}(\mathbb{x}\oplus\mathbb{x}){\mathsf{S}}^{\mathrm{th}}$.
The Hamiltonian $\hat{H}^{\prime}_{{\rm ho}}$ is the one of a multimode
harmonic oscillator, i.e.,

$\hat{H}^{\prime}_{{\rm
ho}}:=\frac{1}{2}\hat{\mathbf{x}}^{\intercal}(\mathbb{x}\oplus\mathbb{x})\hat{\mathbf{x}}=\sum_{j=1}^{n}\frac{\lambda_{j}}{2}(\hat{p}_{j}^{2}+\hat{q}_{j}^{2})$

with $\mathbb{x}=\text{diag}(\lambda_{1},\ldots,\lambda_{n})$ and
$\lambda_{j}\in\mathds{R}$, $j=1,...,n$.

Likewise $\hat{H}$, the Hamiltonian $\hat{H}_{\mathrm{eff}}$ in (108) is also
elliptic, or dynamically stable. However, contrary to the multimode harmonic
oscillator $\hat{H}_{{\rm ho}}$ in the expression (96), the arbitrary
frequencies $\lambda_{j}$ of $\hat{H}^{\prime}_{{\rm ho}}$ can be positive,
negative, or even null. Consequently, $\hat{H}_{\mathrm{eff}}$ can be the null
Hamiltonian and the system still thermalizes to the same state
$\hat{\sigma}^{\mathrm{th}}$. Note, however, that taking a null Hamiltonian is
quite different of considering the interaction representation of the  QOME in
(88), where the density operator turns to
$\hat{\rho}_{t}^{I}=e^{\frac{\imath}{\hbar}\hat{H}_{\mathrm{OME}}t}\hat{\rho}_{t}e^{-\frac{\imath}{\hbar}\hat{H}_{\mathrm{OME}}t}$,
since the frequencies $\tilde{\omega}_{k}$ still correspond to
$\hat{H}_{\mathrm{OME}}$ through the Planck factor $\bar{n}_{k}$.

It is worth to note that the operators $\hat{X}^{\pm}_{k}$ are eigenoperators
of $\hat{H}$ when the master equation (92) of a  GDS thermalizes, see Eq.
(93), but they are also eigenoperators of $\hat{H}_{\mathrm{eff}}$, i.e.,
$[\hat{H},\hat{X}_{j}^{\pm}]=\pm\hbar\lambda_{j}\hat{X}_{j}^{\pm}$ with
$j=1,\ldots,n$. Nevertheless, the frequencies $\omega_{k}$ in the master
equation (92) are associated with $\hat{H}$. Thus, the energetic balance of
the thermalization process is associated with transitions in the Hamiltonian
$\hat{H}$ associated with the confinement of the system.

##  VIII Conclusions

When a  GDS has a thermal equilibrium state $\hat{\sigma}^{\mathrm{th}}$,
every initial condition of the evolution ends up on this state. Without lost
of generality, we can choose
$\expectationvalue{\hat{\mathbf{x}}}^{\mathrm{th}}=0$, so the thermal state
$\hat{\sigma}^{\mathrm{th}}$ must have the form in (41) with $\hat{H}$ in
(42), and the Hamiltonian of the unitary part of the evolution for the GDS
given by $\hat{H}_{\mathrm{eff}}$, of the form in (44) such that
$\hat{H}=\hat{H}_{\mathrm{eff}}$ or $[\hat{H},\hat{H}_{\mathrm{eff}}]=0$.
Another prerequisite for $\hat{\sigma}^{\mathrm{th}}$ came from the
normalization of a density operator: in order to have unity trace, the Hessian
matrix $\mathbb{B}$ of $\hat{H}$ must be positive definite Nicacio (2021b).
This means that $\hat{H}$ is elliptic, i.e., its classical counterpart in (54)
generates a Hamiltonian flow given by the symplectic matrix
$\tilde{\mathsf{S}}_{t}$ in (55) around the elliptic fixed point
$\mathbf{x}=0$. Therefore, $\hat{H}$ is dynamically stable having a discrete
energy spectrum corresponding to bound states.

The non-unitary part of the evolution of a  GDS is determined by the diffusion
and dissipation matrices, $\mathbb{D}$ and $\mathbb{C}$, respectively. These
are the real and imaginary parts, respectively, of the decoherence matrix
$\mathbb{\Gamma}$, see Eq.(11). In (62), we show the relation between the
covariance matrix $\mathbb{V}^{\mathrm{th}}$ of $\hat{\sigma}^{\mathrm{th}}$
and the matrices $\mathbb{D}$ and $\mathbb{C}$ for every  GDS with a thermal
equilibrium state, which has appeared before in Toscano et al. (2021) in a
narrower context. Further, we show that $\mathbb{V}^{\mathrm{th}}$ is
completely determined by the symplectic matrix ${\mathsf{S}}^{\mathrm{th}}$ in
(49), which diagonalizes symplectically the Hessian matrix $\mathbb{B}$ as
shown in (51), and by the symplectic spectrum of $\mathbb{B}$, which contains
the eigenfrequencies $\omega_{j}>0$ of the system Hamiltonian.

However, Eq.(62) does not establish neither the temperature dependence of
$\mathbb{D}$ and $\mathbb{C}$, nor characterize all the different GDSs that
have the same thermal equilibrium state $\hat{\sigma}^{\mathrm{th}}$. This
deficiency can be fixed extending Alicki’s  QDBC (developed for systems
described by discrete Hilbert spaces) to $n$-mode bosonic systems. We have
shown how to implement this extension and the result is that all GDSs
verifying a  QDBC have a thermal equilibrium state whose dynamics is
characterized by symplectically invariant diffusion and dissipation matrices,
see Eq.(72). This condition allows the characterization of the Lindblad
operators of a  GDS satifying a  QDBC in Theorem 5.

The corresponding master equation has the structure of a Quantum Optical
Master Equation (QOME), see Eq. (92) together with Eqs. (93) and (94). As a
consequence, the matrices $\mathbb{D}$ and $\mathbb{C}$ have the specific
structure in (101), which sets explicitly their dependence on temperature.
These expressions also show that different GDSs with the same thermal
equilibrium state differ from each other only on the value of the coupling
constants, namely $\bar{\gamma}_{j}$ in (94), because
${\mathsf{S}}^{\mathrm{th}}$ is determined by the symplectic diagonalization
of $\mathbb{B}$, i.e., the Hessian matrix of $\hat{H}$. Usually the coupling
constants $\bar{\gamma}_{j}$ do not depend on temperature, however we show
that if we allow these to have an arbitrary temperature dependence, then it is
possible to say that every  GDS with thermal equilibrium satisfies a  QDBC.

We also have shown that it is possible to define a pure diffusive regime for
any  GDS satisfying the  QDBC, where the  GDS lost its stationary solution.
This is a specific example of the necessity to balance diffusion and
dissipation to achieve a stationary solution in any quantum dynamical
semigroup Toscano et al. (2021).

Finally, we also show that the contribution of $\hat{H}_{\mathrm{eff}}$ to the
dynamics of a  GDS has no influence in the thermalization process. However, by
determining all $\hat{H}_{\mathrm{eff}}$ such that
$[\hat{H},\hat{H}_{\mathrm{eff}}]=0$, we have determined all Hamiltonians that
allows thermalization. These Hamiltonians are also dynamically stable as
$\hat{H}$, however with arbitrary eigenenergy frequencies. Therefore, the
energetic balance of the thermalization process, determined by $\hat{H}$, is
completely independent of $\hat{H}_{\mathrm{eff}}$. This marks a fundamental
difference between the  QOME for a  GDS satisfying a  QDBC and the usual  QOME
Gardiner and Zoller (2000).

###### Acknowledgements.

The authors are members of the Brazilian National Institute of Science and
Technology for Quantum Information [CNPq INCT-IQ (465469/2014-0)]. FN
acknowledges partial financial support from the Brazilian agency CAPES
[PrInt2019 (88887.468382/2019-00)]. The authors are grateful for the referees’
comments which allow a significant improvement of the manuscript.

##  Appendix A Proof of Lemma 1

Let us prove first that $\mathbb{D}>0$ under the hypothesis of the Lemma.
Since $\mathbb{C}$ is antisymmetric ($\mathbb{C}=-\mathbb{C}^{\intercal}$) and
invertible, the spectral theorem Horn and Johnson (1990) guarantees that there
is an unitary diagonalizing matrix $\mathbb{U}$, such that
$\mathbb{U}\mathbb{C}\mathbb{U}^{\dagger}=(-\imath\mathbb{c})\oplus(\imath\mathbb{c})=\mathbb{C}^{\prime}$,
where $\mathbb{c}=\mbox{diag}(c_{1},\ldots,c_{n})$, with $c_{i}>0$, for
$i=1,\ldots,n$. Thus, using (13) we have
$\mathbb{U}\hbar\mathbb{\Gamma}\mathbb{U}^{\dagger}=\mathbb{D}^{\prime}+\mathbb{C}^{\prime}\geq
0$, where $\mathbb{D}^{\prime}=\mathbb{U}\mathbb{D}\mathbb{U}^{\dagger}$ and
$\mathbb{C}^{\prime}=(\hbar\mathbb{c})\oplus(-\hbar\mathbb{c})$. Now consider
the matrix $\mathbb{\Gamma}^{\intercal}$ and apply the same unitary
transformation, reaching
$\mathbb{U}\mathbb{\Gamma}^{\intercal}\mathbb{U}^{\dagger}=\mathbb{D}^{\prime}-\mathbb{C}^{\prime}\geq
0$, which is positive-semidefinite, since $\mathbb{\Gamma}^{\intercal}$ is
Hermitian and has the same eigenvalues of $\mathbb{\Gamma}$. For a complex
$2n$-dimensional vector $\bf z$, the above matrix inequalities are equivalent
to ${\bf z}^{\dagger}\mathbb{D}^{\prime}{\bf z}\pm{\bf
z}^{\dagger}\mathbb{C}^{\prime}{\bf z}\geq 0$, for all vectors ${\bf z}$,
which implies ${\bf z}^{\dagger}\mathbb{D}^{\prime}{\bf z}>0$, thus
$\mathbb{D}>0$.

To prove that $\mathbb{\Gamma}>0$, it is enough to use that ${\bf
z}^{\dagger}\mathbb{C}{\bf z}=0$ for all vectors $\bf z$, since $\mathbb{C}$
is antisymmetric. Consequently, $\hbar{\bf z}^{\dagger}\mathbb{\Gamma}{\bf
z}={\bf z}^{\dagger}\mathbb{D}{\bf z}$, since $\mathbb{D}>0$ under the
hypothesis of the in Sec.VIII, thus $\mathbb{\Gamma}>0$. This finishes the
proof.

Observing the particular case where $\mathbb{C}=0$, one has
$\mathbb{\Gamma}>0$ iff $\mathbb{D}>0$, which shows that the converse of the
Lemma is not true.

##  Appendix B Proof of Theorem 2

Since $\bar{\Lambda}_{t}$ is a  GDS and satisfies a
$\hat{\bar{\sigma}}^{\mathrm{th}}$-DBC, we have

$\displaystyle\expectationvalue{\bar{\Lambda}_{t}[\hat{A}],\hat{B}}_{\text{GNS}}=\expectationvalue{\hat{A},\bar{\Lambda}_{t}[\hat{B}]}_{\text{GNS}}.$
(109)

Evaluating for $\hat{A}=\hat{1}$, we have $\bar{\Lambda}_{t}[\hat{1}]=\hat{1}$
and thus

$\expectationvalue{\hat{1},\hat{B}}_{\text{GNS}}=\expectationvalue{\hat{1},\bar{\Lambda}_{t}[\hat{B}]}_{\text{GNS}}.$
(110)

However,
$\expectationvalue{\hat{1},\hat{B}}_{\text{GNS}}=\expectationvalue{\hat{\bar{\sigma}},\hat{B}}$
and
$\expectationvalue{\hat{1},\bar{\Lambda}_{t}[\hat{B}]}_{\text{GNS}}=\expectationvalue{\hat{\bar{\sigma}},\bar{\Lambda}_{t}[\hat{B}]}=\expectationvalue{\Lambda_{t}[\hat{\bar{\sigma}}],\hat{B}}$,
where we use that $\Lambda_{t}$ is the adjoint of $\bar{\Lambda}_{t}$ with
respect to the Hilbert-Schmidt scalar product $\langle{\cdot}\rangle$.
Therefore, we arrive to

$\displaystyle\expectationvalue{\hat{\bar{\sigma}}^{\mathrm{th}},\hat{B}}=\expectationvalue{\Lambda_{t}[\hat{\bar{\sigma}}^{\mathrm{th}}],\hat{B}},$
(111)

valid for all operators $\hat{B}$ such the Hilbert-Schmidt scalar product on
both sides is finite. In particular, we can use
$\hat{B}=\hat{R}_{\mathbf{x}}$, a reflection operator defined in Section II.2,
therefore, (111) it is just the equality between the Wigner symbols of
$\hat{\bar{\sigma}}^{\mathrm{th}}$ and
$\Lambda_{t}[\hat{\bar{\sigma}}^{\mathrm{th}}]$, which implies that
$\Lambda_{t}[\hat{\bar{\sigma}}^{\mathrm{th}}]=\hat{\bar{\sigma}}^{\mathrm{th}}$.

##  Appendix C Proof of Theorem 3

Let us start writing

$\displaystyle\expectationvalue{\bar{\Lambda}_{t}[\Xi_{\tau}[\hat{A}]],\hat{B}}_{\text{GNS}}=\Tr(\hat{\bar{\sigma}}^{\mathrm{th}}(\bar{\Lambda}_{t}[\Xi_{\tau}[\hat{A}]])^{\dagger}\hat{B})$
$\displaystyle=\Tr(\hat{\bar{\sigma}}^{\mathrm{th}}(\Xi_{\tau}[\hat{A}])^{\dagger}\bar{\Lambda}_{t}[\hat{B}])=\Tr(\hat{\bar{\sigma}}^{\mathrm{th}}\Xi_{{\tau^{\ast}}}[\hat{A}^{\dagger}]\bar{\Lambda}_{t}[\hat{B}])$
$\displaystyle=\Tr(\hat{\bar{\sigma}}^{\mathrm{th}}\,e^{\frac{\imath}{\hbar}\hat{H}{\tau^{\ast}}}\,\hat{A}^{\dagger}e^{-\frac{\imath}{\hbar}\hat{H}{\tau^{\ast}}}\,\bar{\Lambda}_{t}[\hat{B}])$
$\displaystyle=\Tr(e^{-\frac{\imath}{\hbar}\hat{H}{\tau^{\ast}}}\,\bar{\Lambda}_{t}[\hat{B}^{\dagger}]^{\dagger}\,e^{\frac{\imath}{\hbar}\hat{H}{\tau^{\ast}}}\hat{\bar{\sigma}}^{\mathrm{th}}\hat{A}^{\dagger}),$
(112)

where in the first line we used that $\bar{\Lambda}_{t}$ satisfies a
$\hat{\bar{\sigma}}^{\mathrm{th}}$- DBC and the notation $\tau^{*}$ means
complex conjugation of $\tau$; in the third line, we used that
$[e^{\pm\frac{\imath}{\hbar}\hat{H}\tau},\hat{\bar{\sigma}}^{\mathrm{th}}]=0$
and that $\bar{\Lambda}_{t}$ is a real superoperator, i.e.,
$\bar{\Lambda}_{t}[\hat{B}^{\dagger}]^{\dagger}=\bar{\Lambda}_{t}[\hat{B}]$.

Following Carlen and Maas (2017), we specialize (112) for
${\tau^{\ast}}=-\imath\hbar\beta$, therefore
$e^{\frac{\imath}{\hbar}\hat{H}(-\imath\hbar\beta)}=(\hat{\bar{\sigma}}^{\mathrm{th}})^{-1}$,
then

$\displaystyle\expectationvalue{\bar{\Lambda}_{t}[\Xi_{\imath\hbar\beta}[\hat{A}]],\hat{B}}_{\text{GNS}}=\Tr(\hat{\bar{\sigma}}^{\mathrm{th}}\,\bar{\Lambda}_{t}[\hat{B}^{\dagger}]^{\dagger}\,\hat{A}^{\dagger})$
$\displaystyle=\Tr(\hat{\bar{\sigma}}^{\mathrm{th}}\,(\hat{B}^{\dagger})^{\dagger}\,\bar{\Lambda}_{t}[\hat{A}^{\dagger}])=\Tr(\hat{\bar{\sigma}}^{\mathrm{th}}\,\hat{B}\,\bar{\Lambda}_{t}[\hat{A}]^{\dagger})$
$\displaystyle=\Tr(\bar{\Lambda}_{t}[\hat{A}]^{\dagger}\,\hat{\bar{\sigma}}^{\mathrm{th}}\,\hat{B})=\Tr(\hat{\bar{\sigma}}^{\mathrm{th}}\,(\hat{\bar{\sigma}}^{\mathrm{th}})^{-1}\bar{\Lambda}_{t}[\hat{A}]^{\dagger}\,\hat{\bar{\sigma}}^{\mathrm{th}}\,\hat{B})$
$\displaystyle=\Tr(\hat{\bar{\sigma}}^{\mathrm{th}}\,{(\hat{\bar{\sigma}}^{\mathrm{th}}\bar{\Lambda}_{t}[\hat{A}]\,(\hat{\bar{\sigma}}^{\mathrm{th}})^{-1})}^{\dagger}\,\hat{B})$
$\displaystyle=\Tr(\hat{\bar{\sigma}}^{\mathrm{th}}\,(\Xi_{\imath\hbar\beta}[\bar{\Lambda}_{t}[\hat{A}]])^{\dagger}\,\hat{B})$
$\displaystyle=\expectationvalue{\Xi_{\imath\hbar\beta}[\bar{\Lambda}_{t}[\hat{A}]],\hat{B}}_{\text{GNS}}.$
(113)

Since this identity is valid for any operator $\hat{B}$, we can set it equal
to the reflection operator, i.e., $\hat{B}=\hat{R}_{\mathbf{x}}$. Then, in
this case, Eq.(113) just restates that the Wigner symbols of the operator
$\hat{\bar{\sigma}}^{\mathrm{th}}\bar{\Lambda}_{t}[\Xi_{\imath\hbar\beta}[\hat{A}]]$
and of
$\hat{\bar{\sigma}}^{\mathrm{th}}\Xi_{\imath\hbar\beta}[\bar{\Lambda}_{t}[\hat{A}]]$
are identical, therefore also the operators themselves. For any finite value
$\beta$, the operator $\hat{\bar{\sigma}}^{\mathrm{th}}$ has an inverse, then
it must be true that

$\bar{\Lambda}_{t}[\Xi_{\imath\hbar\beta}[\hat{A}]]=\Xi_{\imath\hbar\beta}[\bar{\Lambda}_{t}[\hat{A}]],$
(114)

for all $\hat{A}$, which means that $\bar{\Lambda}_{t}$ commutes with
$\Xi_{\imath\hbar\beta}$ for all $t$. Noting that
$(\hat{\bar{\sigma}}^{\mathrm{th}})^{\pm\frac{\imath
t}{\hbar\beta}}=(e^{-\beta\hat{H}})^{\pm\frac{\imath
t}{\hbar\beta}}=e^{{\mp}\frac{\imath}{\hbar}t\hat{H}}$, therefore we can
define $(\Xi_{\imath\hbar\beta})^{{-\frac{\imath
t}{\hbar\beta}}}[\cdot]={(\hat{\bar{\sigma}}^{\mathrm{th}})^{-\frac{\imath
t}{\hbar\beta}}\,\cdot\,((\hat{\bar{\sigma}}^{\mathrm{th}})^{-1})^{-\frac{\imath
t}{\hbar\beta}}}=e^{\frac{\imath}{\hbar}\hat{H}t}\,\cdot\,e^{-\frac{\imath}{\hbar}\hat{H}t}=\Xi_{t}[\cdot]$.
But, according to (114), $\bar{\Lambda}_{t}$ commutes with
$\Xi_{\imath\hbar\beta}$ and thus it must commute with any function of
$\Xi_{\imath\hbar\beta}$, then with
$(\Xi_{\imath\hbar\beta})^{\frac{\imath}{\hbar}t}=\Xi_{t}$ for all
$t\in\mathds{C}$. Taking the time derivative on both sides of Eq.(114) and
using
$d\bar{\Lambda}_{t}[\cdot]/dt=\bar{\mathcal{L}}[\bar{\Lambda}_{t}[\cdot]]$, we
get
$\bar{\mathcal{L}}[\bar{\Lambda}_{t}[\Xi_{\imath\hbar\beta}[\hat{A}]]]=\Xi_{\imath\hbar\beta}[\bar{\mathcal{L}}\bar{\Lambda}_{t}[[\hat{A}]]]$,
which shows that $\bar{\mathcal{L}}$ commute with $\Xi_{\imath\hbar\beta}$
when $t=0$. Following the same reasoning as before, $\bar{\mathcal{L}}$
commutes with $\Xi_{t}$.

##  Appendix D Demonstration of Eq.(74)

Let us start from the left-hand side of (74):

$\displaystyle\Xi_{-t}\left[\partialderivative{\hat{\mathbf{x}}^{\intercal}}[\cdot]\right]$
$\displaystyle=$
$\displaystyle\Xi_{-t}\left[\frac{\imath}{\hbar}[(\mathsf{J}\hat{\mathbf{x}})^{\intercal},\cdot]\right]=\frac{\imath}{\hbar}\left(\Xi_{-t}\left[(\mathsf{J}\hat{\mathbf{x}})^{\intercal}\,\cdot\right]-\frac{\imath}{\hbar}\Xi_{-t}\left[\cdot\,(\mathsf{J}\hat{\mathbf{x}})^{\intercal}\right]\right)=\frac{\imath}{\hbar}\left(\Xi_{-t}\left[(\mathsf{J}\hat{\mathbf{x}})^{\intercal}\right]\Xi_{-t}\left[\cdot\right]-\frac{\imath}{\hbar}\Xi_{-t}\left[\cdot\right]\Xi_{-t}\left[(\mathsf{J}\hat{\mathbf{x}})^{\intercal}\right]\right)$
(115) $\displaystyle=$
$\displaystyle\frac{\imath}{\hbar}\left(-\Xi_{-t}\left[\hat{\mathbf{x}}^{\intercal}\right]\,\mathsf{J}\,\Xi_{-t}\left[\cdot\right]+\frac{\imath}{\hbar}\Xi_{-t}\left[\cdot\right]\,\Xi_{-t}\left[\hat{\mathbf{x}}^{\intercal}\right]\,\mathsf{J}\right)=\frac{\imath}{\hbar}\left(-\hat{\mathbf{x}}^{\intercal}\tilde{\mathsf{S}}_{t}^{-\intercal}\,\mathsf{J}\,\Xi_{-t}\left[\cdot\right]+\frac{\imath}{\hbar}\Xi_{-t}\left[\cdot\right]\,\hat{\mathbf{x}}^{\intercal}\tilde{\mathsf{S}}_{t}^{-\intercal}\,\mathsf{J}\right)$
$\displaystyle=$
$\displaystyle\frac{\imath}{\hbar}\left(-\hat{\mathbf{x}}^{\intercal}\,\mathsf{J}\,\tilde{\mathsf{S}}_{t}\,\Xi_{-t}\left[\cdot\right]+\frac{\imath}{\hbar}\Xi_{-t}\left[\cdot\right]\,\hat{\mathbf{x}}^{\intercal}\,\mathsf{J}\tilde{\mathsf{S}}_{t}\right)=\frac{\imath}{\hbar}\left((\mathsf{J}\hat{\mathbf{x}})^{\intercal}\Xi_{-t}\left[\cdot\right]-\frac{\imath}{\hbar}\Xi_{-t}\left[\cdot\right]\,(\mathsf{J}\,\hat{\mathbf{x}})^{\intercal}\right)\tilde{\mathsf{S}}_{t}=\partialderivative{\hat{\mathbf{x}}^{\intercal}}\left[\Xi_{-t}[\cdot]\right]\tilde{\mathsf{S}}_{t},$

where we used that
$\Xi_{-t}[\hat{\mathbf{x}}^{\intercal}]=\hat{\mathbf{x}}\,\tilde{\mathsf{S}}_{t}^{-\intercal}$
and
$\tilde{\mathsf{S}}_{t}^{-\intercal}\mathsf{J}=\mathsf{J}\tilde{\mathsf{S}}_{t}$
because $\tilde{\mathsf{S}}_{t}$ is a symplectic matrix.

##  Appendix E Demonstration of Eq.(78)

In order to prove Eq.(78), note that the decoherence matrix $\mathbb{\Gamma}$,
that stems from Theorem 4, must be invariant under a congruence with
$\tilde{\mathsf{S}}_{t}$ in (55). This follows using (72) into the definition
of the matrices $\mathbb{D}$ and $\mathbb{C}$ in (10), and the fact that
$\tilde{\mathsf{S}}_{t}$ is a real matrix.

Now, using Eq.(55) we can write
$\tilde{\mathsf{S}}_{t}(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{-1}=(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{-1}e^{(\mathsf{J}\mathbb{B})_{\rm
d}\,t}$, so the column vectors
$(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{-1}|_{k}=\bar{\mathbf{l}}_{k}$ for
$k=1,\ldots,2n$ are the complex eigenvectors of $\tilde{\mathsf{S}}_{t}$.
Respecting the block order in the matrix in (52), we can write

$\displaystyle\tilde{\mathsf{S}}_{t}\,\bar{\mathbf{l}}_{j}$ $\displaystyle=$
$\displaystyle e^{\imath\omega_{j}t}\bar{\mathbf{l}}_{j},$ (116a)
$\displaystyle\tilde{\mathsf{S}}_{t}\bar{\mathbf{l}}_{n+j}$ $\displaystyle=$
$\displaystyle\tilde{\mathsf{S}}_{t}\bar{\mathbf{l}}_{j}^{*}=e^{-\imath\omega_{j}t}\bar{\mathbf{l}}_{j}^{*},$
(116b)

where $j=1,\ldots,n$ and $\omega_{j}>0$ are the frequencies corresponding to
the symplectic spectra of $\mathbb{B}$ in (51). According to (116), it is
clear that the canonical form of $\mathbb{\Gamma}$ in (78) is manifestly
invariant under a congruence through $\tilde{\mathsf{S}}_{t}$ and is positive-
definite, as it must be.

Let us now prove that the matrices $\mathbb{s}$ and $\mathbb{r}$ in (78)
satisfy the relation in Eq.(80). We first note that the expression of
$\mathbb{\Gamma}$ in (78) can be rewritten as in (84), which allows us to show
that the matrix $({\mathsf{S}}^{\mathrm{th}})^{\intercal}\mathbb{Q}^{\dagger}$
diagonalizes simultaneously the Hamiltonian matrix $\mathsf{J}\mathbb{D}$ and
the skew Hamiltonian matrix $\mathsf{J}\mathbb{C}$, therefore, we must have
$[\mathsf{J}\mathbb{D},\mathsf{J}\mathbb{C}]=0$. Indeed, from (84), we have

$\displaystyle\mathsf{J}\mathbb{D}$ $\displaystyle=$
$\displaystyle\hbar\mathsf{J}({\mathsf{S}}^{\mathrm{th}})^{-1}\mathbb{\Lambda}_{r}({\mathsf{S}}^{\mathrm{th}})^{-\intercal}=\hbar({\mathsf{S}}^{\mathrm{th}})^{\intercal}\mathsf{J}\mathbb{\Lambda}_{r}({\mathsf{S}}^{\mathrm{th}})^{-\intercal}$
(117) $\displaystyle=$
$\displaystyle\hbar({\mathsf{S}}^{\mathrm{th}})^{\intercal}\mathbb{Q}^{\dagger}\mathbb{Q}\mathsf{J}\mathbb{\Lambda}_{r}\mathbb{Q}^{\dagger}\mathbb{Q}({\mathsf{S}}^{\mathrm{th}})^{-\intercal}$
$\displaystyle=$
$\displaystyle({\mathsf{S}}^{\mathrm{th}})^{\intercal}\mathbb{Q}^{\dagger}(\mathsf{J}\mathbb{D})_{\rm
d}(({\mathsf{S}}^{\mathrm{th}})^{\intercal}\mathbb{Q}^{\dagger})^{-1},$

where $\mathbb{Q}$ is the unitary matrix defined in (26). In the above steps,
the symplectic condition
$\mathsf{J}({\mathsf{S}}^{\mathrm{th}})^{-1}=({\mathsf{S}}^{\mathrm{th}})^{\intercal}\mathsf{J}$
and the fact that
$\hbar\mathbb{Q}\mathsf{J}\mathbb{\Lambda}_{r}\mathbb{Q}^{\dagger}$ is a
diagonal matrix were employed. Since Eq.(117) is a similarity transformation
the matrix $\hbar\mathbb{Q}\mathsf{J}\mathbb{\Lambda}_{r}\mathbb{Q}^{\dagger}$
has the eigenvalues of $\mathsf{J}\mathbb{D}$ in its diagonal. Therefore, we
can write

$(\mathsf{J}\mathbb{D})_{\rm
d}=\frac{\imath\hbar}{2}\left((|\mathbb{s}|^{2}+|\mathbb{r}|^{2})\oplus\left(-|\mathbb{s}|^{2}-|\mathbb{r}|^{2}\right)\right).$
(118)

It is worth to note that, according to Lemma 2, the diagonal matrix
$\mathbb{d}=\tfrac{1}{2}(|\mathbb{s}|^{2}+|\mathbb{r}|^{2})$ contains the
symplectic spectrum of $\mathbb{D}$. In an analogous way, we have

$\displaystyle\mathsf{J}\mathbb{C}$ $\displaystyle=$
$\displaystyle\mathsf{J}({\mathsf{S}}^{\mathrm{th}})^{-1}(-\mathsf{J}\mathbb{\Lambda}_{i})({\mathsf{S}}^{\mathrm{th}})^{-\intercal}=({\mathsf{S}}^{\mathrm{th}})^{\intercal}\mathbb{\Lambda}_{i}({\mathsf{S}}^{\mathrm{th}})^{-\intercal}$
(119) $\displaystyle=$
$\displaystyle({\mathsf{S}}^{\mathrm{th}})^{\intercal}\mathbb{Q}^{\dagger}\mathbb{\Lambda}_{i}\mathbb{Q}({\mathsf{S}}^{\mathrm{th}})^{-\intercal}$
$\displaystyle=$
$\displaystyle({\mathsf{S}}^{\mathrm{th}})^{\intercal}\mathbb{Q}^{\dagger}(\mathsf{J}\mathbb{C})_{\rm
d}(({\mathsf{S}}^{\mathrm{th}})^{\intercal}\mathbb{Q}^{\dagger})^{-1},$

where the eigenvalue matrix corresponding to $\mathsf{J}\mathbb{C}$ is

$(\mathsf{J}\mathbb{C})_{\rm
d}=\mathbb{Q}^{\dagger}\mathbb{\Lambda}_{i}\mathbb{Q}=\mathbb{\Lambda}_{i}=\mathbb{jc}.$
(120)

From (117) and (119), we immediately realize that $\mathsf{J}\mathbb{D}$
commutes with $\mathsf{J}\mathbb{C}$, as we wanted to prove.

The condition $[\mathsf{J}\mathbb{D},\mathsf{J}\mathbb{C}]=0$ in Theorem 1
(proved in the main text) determines univocally the covariance matrix
$\mathbb{V}^{\mathrm{th}}$ of a thermal equilibrium state
$\hat{\sigma}^{\mathrm{th}}$. Therefore, using the expressions (62) and (86),
we can write

$\mathbb{V}^{\mathrm{th}}=({\mathsf{S}}^{\mathrm{th}})^{-1}\frac{1}{2}\mathbb{\Lambda}_{r}\mathbb{\Lambda}_{i}^{-1}({\mathsf{S}}^{\mathrm{th}})^{-\intercal},$
(121)

where we used the symplectic condition
$\mathsf{J}({\mathsf{S}}^{\mathrm{th}})=({\mathsf{S}}^{\mathrm{th}})^{-\intercal}\mathsf{J}$.
Comparing with Eq.(49), we immediately recognize in Eq.(121) the symplectic
spectrum of $\mathbb{V}^{\mathrm{th}}$, i.e.,

$\frac{1}{2}\mathbb{\Lambda}_{r}\mathbb{\Lambda}_{i}^{-1}=\mathbb{k}^{\mathrm{th}}\oplus\mathbb{k}^{\mathrm{th}}.$
(122)

Equating the matrix elements on both sides of (122), we obtain
$\frac{1}{2}\frac{1+a_{j}}{1-a_{j}}=\kappa_{j}^{\rm th}$ with
$a_{j}=\frac{|r_{j}|^{2}}{|s_{j}|^{2}}$ and $j=1,\ldots,n$. Inverting these
equations and rewriting in a matrix structure, we get

$\frac{|\mathbb{r}|^{2}}{|\mathbb{s}|^{2}}=\frac{2\mathbb{k}^{\mathrm{th}}-\mathbb{1}}{2\mathbb{k}^{\mathrm{th}}+\mathbb{1}}.$
(123)

Now using (31), (45), and (51), we arrive to the equality
$g(\mathbb{k}^{\mathrm{th}})=\hbar\beta\mathbb{w}$. Then, employing the
identity $\exp[-2\coth^{-1}(2x)]=\frac{2x-1}{2x+1}$ for $x\geq\frac{1}{2}$, we
can rewrite Eq.(123) as

$\frac{|\mathbb{r}|^{2}}{|\mathbb{s}|^{2}}=e^{-\hbar\beta\mathbb{w}},$ (124)

which finally implies the relation in (80).

##  Appendix F Proof of Eq.(93)

Let us prove an equivalent statement of Eq.(93), i.e., in a  GDS satisfying a
$\hat{\bar{\sigma}}^{\mathrm{th}}$-DBC we have

$[\hat{H},\hat{L}_{j}]=-\hbar\omega_{j}\hat{L}_{j}\Rightarrow[\hat{H},\hat{L}^{\dagger}_{j}]=\hbar\omega_{j}\hat{L}^{\dagger}_{j},$
(125)

for $j=1,\ldots,n$, $\hat{H}$ given in (42), and the Lindblad operators
$\hat{L}_{j}$ in (90).

For a generic symmetric matrix $\mathbb{X}$ and a generic complex vector
$\mathbf{l}$, it is straightforward to show that

$[\tfrac{1}{2}\hat{\mathbf{x}}^{\intercal}\mathbb{X}\hat{\mathbf{x}},\mathbf{l}^{\intercal}\mathsf{J}\hat{\mathbf{x}}]=\imath\hbar\mathbf{l}^{\intercal}\mathbb{X}\hat{\mathbf{x}}=\imath\hbar\mathbf{l}^{\intercal}\mathbb{X}\mathsf{J}^{\intercal}\mathsf{J}\hat{\mathbf{x}},$
(126)

which is the commutation relation for a generic quadratic Hamiltonian with a
linear Lindblad operator. But, if the  GDS satisfies a
$\hat{\bar{\sigma}}^{\mathrm{th}}$- DBC according to Eq.(55),
$\bar{\mathbf{l}}_{j}$ appearing in (91) is an eigenvector of the matrix
$\mathsf{J}\mathbb{B}$, i.e.,
$\mathsf{J}\mathbb{B}\bar{\mathbf{l}}_{j}=\imath\omega_{j}\bar{\mathbf{l}}_{j}$,
which is the same as say that
$\bar{\mathbf{l}}_{j}=(\mathbb{Q}{\mathsf{S}}^{\mathrm{th}})^{-1}|_{j}$, see
Appendix E. Therefore, using (126) with
$\mathbf{l}_{j}^{\intercal}\mathbb{B}\mathsf{J}=\imath\omega_{j}\mathbf{l}_{j}^{\intercal}$
we arrive to (125).

##  Appendix G Proof of Eq.(94)

Let us start noting that for the thermal state $\hat{\sigma}^{\mathrm{th}}$ in
(41), or equivalently in (99), the eigenenergy distribution is the usual
(multimode) Planck distribution:

$\displaystyle\\!\\!P({\bf
n})=\expectationvalue{\hat{\sigma}^{\mathrm{th}}}{\phi_{\bf
n}}=\expectationvalue{\hat{\sigma}_{\rm ho}^{\mathrm{th}}}{{\bf
n}}=\prod_{j=1}^{n}P_{n_{j}}(\bar{n}_{j}),$ (127)

which is nothing but a consequence of Eqs.(96) and (98). In the above
equation,

$\displaystyle
P_{n_{j}}(\bar{n}_{j})=\frac{1}{\bar{n}_{j}+1}\left(\frac{\bar{n}_{j}}{\bar{n}_{j}+1}\right)^{n_{j}},$
(128)

where $\bar{n}_{j}$ is in (95). Note that
$\bar{n}_{j}(\bar{n}_{j}+1)^{-1}=e^{-\hbar\beta\omega_{j}}$.

This distribution can be recovered considering the stationary regime
$\derivative{\hat{\rho}_{t}}{t}=0$ in (92), subsequently taking the diagonal
matrix elements in the eigenbasis {$\ket{\phi_{\bf n}}\\}$ and using the
commutation relation $[\hat{X}^{-}_{j},\hat{X}^{+}_{j}]=1$, indeed

$\displaystyle\sum_{j=1}^{n}\frac{|s_{j}|^{2}}{\hbar}\left((n_{j}+1)\,{P}_{n{{}_{j}}+1}-n_{j}{P}_{n{{}_{j}}}\right)+$
$\displaystyle\sum_{j=1}^{n}\frac{|s_{j}|^{2}}{\hbar}e^{-\hbar\beta\omega_{j}}\left(n_{j}\;{P}_{n{{}_{j}}-1}-(n_{j}+1){P}_{n{{}_{j}}}\right)=0,$

where ${P}_{n{{}_{j}}}=\langle\phi_{\bf
n}|\hat{\sigma}^{\mathrm{th}}|\phi_{\bf n}\rangle$ and ${P}_{n{{}_{j}}\pm
1}=\langle\phi_{{\bf n}_{j}^{\pm}}|\hat{\sigma}^{\mathrm{th}}|\phi_{{\bf
n}_{j}^{\pm}}\rangle$. As can be directly checked, the solution is

${P}_{n{{}_{j}}}=\frac{\hbar\bar{\gamma}_{j}}{|s_{j}|^{2}}e^{-\hbar\beta\omega_{j}n_{j}}.$
(129)

The constant $\bar{\gamma}_{j}$ is included in order to
${\hbar\bar{\gamma}_{j}}/{|s_{j}|^{2}}$ be dimensionless. Comparing Eqs. (129)
and (128) we arrive to (94).

##  Appendix H The algebra of commuting elliptic Hamiltonian matrices

In this appendix we will prove the following theorem:

###### Theorem 6

Consider a positive definite symmetric matrix $\mathbb{B}$ and a symmetric
matrix $\mathbb{B}^{\prime}$ such that
$[\mathsf{J}\mathbb{B},\mathsf{J}\mathbb{B}^{\prime}]=0$, then

$\mathsf{J}\mathbb{B}^{\prime}=(\mathbb{Q}\mathsf{S})^{-1}\,\imath\mathbb{x}\oplus(-\imath\mathbb{x})\,(\mathbb{Q}\mathsf{S}),$
(130)

where $\mathbb{Q}$ is defined in (26),
$\mathbb{x}:={\rm{diag}}(\lambda_{1},\ldots,\lambda_{n})$ with
$\lambda_{i}\in\mathds{R}$ for $i=1,\ldots,n$, and $\mathsf{S}$ is a
symplectic matrix that diagonalizes symplectically $\mathbb{B}$, i.e.,

$\mathsf{S}^{-\intercal}\mathbb{B}\mathsf{S}=\mathbb{w}\oplus\mathbb{w},\,\,\,\mathbb{w}={\rm{diag}}(\omega_{1},\ldots,\omega_{n})$
(131)

for $\omega_{i}>0$ and $i=1,\ldots,n$. Further, Eq.(130) is equivalent to a
“symplectic diagonalization” of $\mathbb{B}^{\prime}$:

$\mathsf{S}^{-\intercal}\mathbb{B}^{\prime}\mathsf{S}^{-1}=\mathbb{x}\oplus\mathbb{x}.$
(132)

Let $\tilde{\mathsf{S}}$ be a symplectic matrix such that
$\tilde{\mathsf{S}}^{-\intercal}\mathbb{B}\tilde{\mathsf{S}}=\mathbb{w}\oplus\mathbb{w}$,
then

$\displaystyle[\mathsf{J}\mathbb{B},\mathsf{J}\mathbb{B}^{\prime}]$
$\displaystyle=$
$\displaystyle[\mathsf{J}{\tilde{\mathsf{S}}}^{\intercal}(\mathbb{w}\oplus\mathbb{w})\tilde{\mathsf{S}},\mathsf{J}\mathbb{B}^{\prime}]$
(133) $\displaystyle=$
$\displaystyle\mathsf{J}\tilde{\mathsf{S}}^{\intercal}\,[\mathsf{J}(\mathbb{w}\oplus\mathbb{w}),{\tilde{\mathsf{S}}}^{-\intercal}\mathbb{B}^{\prime}{\tilde{\mathsf{S}}}^{-1}]\,{\tilde{\mathsf{S}}},$

where we used the symplectic condition
${\tilde{\mathsf{S}}}\mathsf{J}=\mathsf{J}{\tilde{\mathsf{S}}}^{-\intercal}$
and
$\mathsf{J}(\mathbb{w}\oplus\mathbb{w})=(\mathbb{w}\oplus\mathbb{w})\mathsf{J}$.
Writing

$\tilde{\mathbb{B}}:={\tilde{\mathsf{S}}}^{-\intercal}\mathbb{B}^{\prime}{\tilde{\mathsf{S}}}^{-1}=\begin{pmatrix}\tilde{\mathbb{x}}&\tilde{\mathbb{y}}\\\
\tilde{\mathbb{y}}^{\intercal}&\tilde{\mathbb{z}}\\\ \end{pmatrix},$

where $\tilde{\mathbb{x}},\tilde{\mathbb{y}}$ and $\tilde{\mathbb{z}}$ are
$n\times n$ real matrices such that
$\tilde{\mathbb{x}}=\tilde{\mathbb{x}}^{\intercal}$ and
$\mathbb{z}=\mathbb{z}^{\intercal}$, and using the commutation relation
$[\mathsf{J}\mathbb{B},\mathsf{J}\mathbb{B}^{\prime}]=0$ in (133) one attains

$\tilde{\mathbb{y}}^{\top}=-\mathbb{w}^{-1}\tilde{\mathbb{y}}\mathbb{w}=-\mathbb{w}\tilde{\mathbb{y}}\mathbb{w}^{-1}\,\,\,{\text{and}}\,\,\,\tilde{\mathbb{z}}=\mathbb{w}^{-1}\tilde{\mathbb{x}}\mathbb{w}=\mathbb{w}\tilde{\mathbb{x}}\mathbb{w}^{-1};$

consequently,
$[\tilde{\mathbb{y}},\mathbb{w}]=[\tilde{\mathbb{x}},\mathbb{w}]=0$, thus
$\tilde{\mathbb{x}}=\tilde{\mathbb{z}}$ and $\tilde{\mathbb{y}}$ is skew-
symmetric $\tilde{\mathbb{y}}^{\intercal}=-\tilde{\mathbb{y}}$. All of these
relations enable us to write

$\tilde{\mathbb{B}}:={\tilde{\mathsf{S}}}^{-\intercal}\mathbb{B}^{\prime}{\tilde{\mathsf{S}}}^{-1}=\begin{pmatrix}\tilde{\mathbb{x}}&\tilde{\mathbb{y}}\\\
-\tilde{\mathbb{y}}&\tilde{\mathbb{x}}\\\ \end{pmatrix}.$ (134)

Multiplying by $\mathsf{J}$ from left, considering the symplecticity of
$\tilde{\mathsf{S}}$, and using $\mathbb{Q}$, last equation is equivalently
rewritten as

$(\mathbb{Q}\tilde{\mathsf{S}})\,\mathsf{J}\mathbb{B}^{\prime}\,(\mathbb{Q}\tilde{\mathsf{S}})^{-1}=\mathbb{Q}\mathsf{J}\tilde{\mathbb{B}}\mathbb{Q}^{\dagger}=(\imath\tilde{\mathbb{x}}-\tilde{\mathbb{y}})\oplus(-\imath\tilde{\mathbb{x}}-\tilde{\mathbb{y}}).$
(135)

The above particular block structure is a consequence of the degenerated
structure of the diagonal matrix $\mathbb{w}\oplus\mathbb{w}$, where each
diagonal element is at least doubly-degenerated.

The two blocks in the matrix of rightmost equality in (135) are skew-Hermitian
and moreover they are complex conjugate of each other. Recalling that a skew-
Hermitian matrix has pure imaginary (possibly null) eigenvalues and is
unitarily diagonalizable Horn and Johnson (2013), then there exists an unitary
matrix $\mathbb{u}$ such that
$\mathbb{u}(\imath\tilde{\mathbb{x}}-\tilde{\mathbb{y}})\mathbb{u}^{\dagger}=\imath\mathbb{x}$,
where $\mathbb{x}=\text{diag}(\lambda_{1},\ldots,\lambda_{n})$ is the diagonal
matrix containing the eigenvalues $\lambda_{j}\in\mathds{R}$ of the matrix
$\imath\tilde{\mathbb{x}}-\tilde{\mathbb{y}}$.

Defining $\mathsf{S}=\mathsf{R}\,\tilde{\mathsf{S}}$ with

$\mathsf{R}=\mathbb{Q}\begin{pmatrix}\mathbb{u}&0\\\ 0&\mathbb{u}^{*}\\\
\end{pmatrix}\mathbb{Q}^{\dagger}\,\,\in{\rm Sp}(2n,\mathds{R})\;\cap\;{\rm
O}(2n)$ (136)

a real symplectic and orthogonal matrix, from (135) we can write

$\mathbb{Q}\mathsf{S}\;\mathsf{J}\mathbb{B}^{\prime}\;(\mathbb{Q}\mathsf{S})^{-1}=\,\imath\mathbb{x}\oplus(-\imath\mathbb{x}),$
(137)

which is Eq.(130). Since $\mathsf{S}$ is symplectic and using that
$\mathbb{Q}^{\dagger}\,(\imath\mathbb{x}\oplus(-\imath\mathbb{x}))\mathbb{Q}=\mathsf{J}(\mathbb{x}\oplus\mathbb{x})$,
we immediately recover Eq.(132) from (137).

However it remains to prove that $\mathsf{S}$ satisfies (131). If the
symplectic spectrum in $\mathbb{w}$ is non-degenerate, conditions
$[\tilde{\mathbb{y}},\mathbb{w}]=[\tilde{\mathbb{x}},\mathbb{w}]=0$,
$\tilde{\mathbb{x}}^{\intercal}=\tilde{\mathbb{x}}$, and
$\tilde{\mathbb{y}}^{\intercal}=-\tilde{\mathbb{y}}$ imply
$\tilde{\mathbb{x}}=\mathbb{x}=\text{diag}(\lambda_{1},\ldots,\lambda_{n})$
and $\mathbb{y}=0$, so one can choose $\mathbb{u}=\mathbb{1}$ such that
$\mathsf{R}=\mathbb{1}$ and $\mathsf{S}=\tilde{\mathsf{S}}$; consequently
$\mathsf{S}$ satisfies (131), as claimed. When the symplectic spectrum in
$\mathbb{w}$ is degenerate, $\tilde{\mathbb{x}}$ is diagonal outside the
degenerate subspaces of $\mathbb{w}$, while $\mathbb{y}$ is null outside the
same subspaces. Therefore, the unitary matrix $\mathbb{u}$ have to diagonalize
$\imath\tilde{\mathbb{x}}-\tilde{\mathbb{y}}$ only inside the degenerate
subspaces of $\mathbb{w}$. This is possible choosing $\mathbb{u}$ block
diagonal such
$\mathbb{u}\imath\tilde{\mathbb{x}}-\tilde{\mathbb{y}}\mathbb{u}^{\dagger}=\imath\mathbb{x}$.
In this case we also have that
$\mathsf{R}^{\intercal}\mathbb{w}\oplus\mathbb{w}\,\mathsf{R}=\mathbb{w}\oplus\mathbb{w}$,
and therefore $\mathsf{S}=\mathsf{R}\,\tilde{\mathsf{S}}$ also diagonalize
symplectically $\mathbb{B}$, i.e.,

$\mathbb{B}=\mathsf{S}^{\intercal}\mathbb{w}\oplus\mathbb{w}\mathsf{S}=\tilde{\mathsf{S}}^{\intercal}\mathsf{R}^{\intercal}\mathbb{w}\oplus\mathbb{w}\,\mathsf{R}\tilde{\mathsf{S}}=\tilde{\mathsf{S}}^{\intercal}\mathbb{w}\oplus\mathbb{w}\,\tilde{\mathsf{S}},$
(138)

which shows that $\mathsf{S}$ satisfies (131) for the degenerate case, with
this we finish the proof of the theorem.

## References

  * Holevo et al. (1999) A. S. Holevo, M. Sohma, and O. Hirota, Phys. Rev. A 59, 1820 (1999).
  * Holevo and Werner (2001) A. S. Holevo and R. F. Werner, Phys. Rev. A 63, 032312 (2001).
  * Holevo (2002) A. S. Holevo, Sending quantum information with gaussian states, in  Quantum Communication, Computing, and Measurement 2, edited by P. Kumar, G. M. D’Ariano, and O. Hirota (Springer US, Boston, MA, 2002) p. 75–82.
  * Cerf et al. (2007) N. J. Cerf, G. Leuchs, and E. S. Polzik,  Quantum Information with Continuous Variables of Atoms and Light (Imperial College Press, London, 2007) .
  * Caruso et al. (2008) F. Caruso, J. Eisert, V. Giovannetti, and A. S. Holevo, New Journal of Physics 10, 083030 (2008).
  * Holevo (2019) A. S. Holevo,  Quantum Systems, Channels, Information (De Gruyter, 2019).
  * 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, Reviews of Modern Physics 84, 621 (2012).
  * Adesso et al. (2014) G. Adesso, S. Ragy, and A. R. Lee, Open Systems & Information Dynamics 21, 1440001 (2014).
  * Wolf et al. (2006) M. M. Wolf, G. Giedke, and J. I. Cirac, Physical Review Letters 96, 080502 (2006).
  * Heinosaari et al. (2010) T. Heinosaari, A. Holevo, and M. Wolf, Quantum Information and Computation 10, 619–635 (2010).
  * Toscano et al. (2021) F. Toscano, G. M. Bosyk, S. Zozor, and M. Portesi, Phys. Rev. A 104, 062207 (2021).
  * Giovannetti et al. (2010) V. Giovannetti, A. S. Holevo, S. Lloyd, and L. Maccone, Journal of Physics A 43, 415305 (2010).
  * De Palma et al. (2016) G. De Palma, D. Trevisan, and V. Giovannetti, IEEE Transactions on Information Theory 62, 2895 (2016).
  * Alicki and Lendi (2007) R. Alicki and K. Lendi, Quantum dynamical semigroups and applications, Lecture notes in physics (Springer, Berlin, 2007).
  * Vanheuverzwijn (1978) P. Vanheuverzwijn, Annales de l’I.H.P. Physique théorique 29, 123 (1978).
  * Vanheuverzwijn (1979) P. Vanheuverzwijn, Annales de l’I.H.P. Physique théorique 30, 83 (1979).
  * Demoen et al. (1979) B. Demoen, P. Vanheuverzwijn, and A. Verbeure, Reports on Mathematical Physics 15, 27–39 (1979).
  * Sǎndulescu and Scutaru (1987) A. Sǎndulescu and H. Scutaru, Annals of Physics 173, 277 (1987).
  * Carlen and Maas (2017) E. A. Carlen and J. Maas, Journal of Functional Analysis 273, 1 (2017).
  * Lindblad (1976) G. Lindblad, Communications in Mathematical Physics 48, 119 (1976).
  * Gorini et al. (1976) V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Journal of Mathematical Physics 17, 821 (1976).
  * Tarasov (2008) V. Tarasov, Quantum mechanics of non-hamiltonian and dissipative systems (Elsevier, 2008).
  * Nicacio et al. (2016) F. Nicacio, M. Paternostro, and A. Ferraro, Physical Review A 94, 052129 (2016).
  * Ozorio de Almeida (1998) A. M. Ozorio de Almeida, Physics Reports 295, 265 (1998).
  * Nicacio et al. (2010) F. Nicacio, R. N. P. Maia, F.Toscano, and R. O. Vallejos, Physics Letters A 374, 4385 (2010).
  * Frigerio (1977) A. Frigerio, Letters in Mathematical Physics 2, 79–87 (1977).
  * Frigerio (1978) A. Frigerio, Communications in Mathematical Physics 63, 269–276 (1978).
  * Carmichael (1999) J. H. Carmichael, Statistical Methods in Quantum Optics 1, Master Equations and Fokker-Plank Equations (Springer-Verlag, 1999).
  * Alicki (1976) R. Alicki, Reports on Mathematical Physics 10, 249–258 (1976).
  * Gardiner and Zoller (2000) C. W. Gardiner and P. Zoller, Quantum noise, 2nd ed., Springer series in synergetics (Springer, Berlin ; Heidelberg [u.a.], 2000).
  * Breuer et al. (2002) H. Breuer, P. Breuer, F. Petruccione, and S. Petruccione,  The Theory of Open Quantum Systems (Oxford University Press, 2002).
  * Wiseman and Milburn (2009) H. M. Wiseman and G. J. Milburn,  Quantum Measurement and Control (Cambridge University Press, 2009).
  * Banchi et al. (2015) L. Banchi, S. L. Braunstein, and S. Pirandola, Physical Review Letters 115, 260501 (2015).
  * Simon et al. (1994) R. Simon, N. Mukunda, and B. Dutta, Phys. Rev. A 49, 1567 (1994).
  * Nicacio (2021a) F. Nicacio, American Journal of Physics 89, 1139 (2021a), https://doi.org/10.1119/10.0005944 .
  * Nicacio (2021b) F. Nicacio, Journal of Physics A: Mathematical and Theoretical 54, 055004 (2021b).
  * de Gosson (2006) M. de Gosson,  Symplectic Geometry and Quantum Mechanics, Operator Theory: Advances and Applications (Birkhäuser Basel, 2006).
  * Risken (1996) H. Risken, The Fokker-Planck Equation: Methods of Solution and Applications (Springer, 1996).
  * Horn and Johnson (1991) R. A. Horn and C. R. Johnson,  Topics in Matrix Analysis (Cambridge University Press, 1991).
  * Valero-Toranzo et al. (2018) I. Valero-Toranzo, S. Zozor, and J. Brossier, IEEE Transactions on Information Theory 64, 6743 (2018).
  * Bialas and Gora (2015) S. Bialas and M. Gora, Bulletin of the Polish Academy of Sciences Technical Sciences 63, 163–168 (2015).
  * Fassbender et al. (1999) H. Fassbender, D. S. Mackey, N. Mackey, and H. Xu, Linear Algebra and its Applications 287, 125–159 (1999).
  * Gardiner and Zoller (2004) C. W. Gardiner and P. Zoller, Quantum noise : a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics, 3rd ed. (Springer, 2004).
  * Toscano et al. (2005) F. Toscano, R. L. de Matos Filho, and L. Davidovich, Physical Review A 71, 010101(R) (2005).
  * Horn and Johnson (1990) R. A. Horn and C. R. Johnson, Topics in Matrix Analysis (Cambridge University Press, New York, 1994).
  * Horn and Johnson (2013) R. A. Horn and C. R. Johnson, Matrix Analysis, 2nd ed. (Cambridge University Press, Cambridge; New York, 2013).