# Time-extensive classical and quantum correlations in thermal machines

Milton Aguilar<EMAIL_ADDRESS>Departamento de Física, FCEyN, UBA, Pabellón 1,
Ciudad Universitaria, 1428 Buenos Aires, Argentina  Instituto de Física de
Buenos Aires, UBA CONICET, Pabellón 1, Ciudad Universitaria, 1428 Buenos
Aires, Argentina  Juan Pablo Paz<EMAIL_ADDRESS>Departamento de Física, FCEyN,
UBA, Pabellón 1, Ciudad Universitaria, 1428 Buenos Aires, Argentina  Instituto
de Física de Buenos Aires, UBA CONICET, Pabellón 1, Ciudad Universitaria, 1428
Buenos Aires, Argentina

###### Abstract

We study intraenvironmental classical and quantum correlations in a thermal
machine, which is modeled as a driven quantum system coupled with thermal
reservoirs. We compute the mutual information, the quantum discord, and the
entanglement between two parts of the environment formed by oscillators
centered around two different frequencies. We show that there are only two
processes that generate time-extensive correlations in the long-time limit.
First, there is a resonant process which is responsible for the transport of
excitations between different environmental modes due to the absorption (or
emission) of energy from (or into) the driving field. Second, there is a
nonresonant process that transforms the energy from the external driving into
pairs of excitations in two environmental modes. We show that there is a
regime when the mutual information and the quantum discord between the parts
of the environment correlated by these two processes grow quadratically in
time, while entanglement production is time-extensive.

##  I Introduction

The tools and resources of information theory have brought astounding
advancements to our understanding of the quantum world, and have permeated all
the layers of quantum theory. In recent years, the connection between quantum
information [1] and quantum thermodynamics [2] has been extensively explored,
and in the heart of it all lies the study of correlations between physical
systems. This has provided, for example, much necessary insight into the
emergence of thermodynamic irreversibility from unitary dynamics [3, 4], and a
generalization of the laws of thermodynamics that resolve their apparent
violations in correlated scenarios [5] (such as anomalous heat flows from cold
to hot baths [6], and memory erasure accompanied by work extraction instead of
heat dissipation [7]).

In view of their seemingly fundamental role, it is natural to wonder if
correlations can be exploited as a thermodynamic resource. In the context of
single-shot thermodynamics [8], both classical and quantum correlations have
been studied in their role as a resource in state formation or work extraction
[9, 10, 11, 12]. In addition, considering that separability has often been
used synonymously with classicality, entanglement in particular has been put
in the spotlight in order to understand how nonclassical correlations affect
thermal machines [13, 14, 15, 16].

Understanding the nature and creation of correlations is a necessary task to
further advance our knowledge of the foundations of thermodynamics. In this
work we address these issues by studying how and what type of correlations are
built up between different parts of the environment of a thermal machine while
it is continuously driven beyond the transient regime. We will study a
periodically driven system $\mathcal{S}$ in contact with two bosonic
environments, $\mathcal{E}_{R}$ and $\mathcal{E}_{L}$, at different
temperatures. First, we show that there are only two physical processes that
constantly create intraenvironmental time-extensive correlations in the long-
time limit. On the one hand, the driving field can transport an excitation
from one mode to another one. On the other, the energy of the driving can
split and create a pair of excitations that are dumped in different modes. The
first process is none other than the quantum manifestation of the heat current
from classical thermodynamics, that flows either from the hotter reservoir to
the colder one in a heat engine or the reverse in a heat pump. The second one
is at the core of the third law of thermodynamics: it always heats up the
environment, preventing it from reaching absolute zero temperature. As a
consequence of these two processes, time extensive classical and quantum
correlations are established only between parts of the environment centered
around frequencies such that their sum or difference is a multiple of the
driving frequency $\omega_{d}$. And second, we analyze the nature of those
correlations. We present simple, analytical expressions for the quantification
of classical and quantum correlations and analyze their relation and their
interplay. We also show that the “classical” transport of excitations swaps
the entanglement created by the creation of a pair of excitations to other
modes.

The study of intraenvironmental correlations is not plentiful in the
literature of quantum thermodynamics. Recent findings [3, 4] show that
thermodynamic irreversibility is intrinsically linked to the development of
correlations in the environment. Indeed, in Refs. [3, 4] it was shown that
these correlations are responsible for the entropy production in the long-time
limit, using a fermionic system as an example. Our work can also be
interpreted as providing support for this idea. Thus, we obtain similar
results in a bosonic model, which is the paradigmatic model of a thermal
machine. This paper also contains a generalization of previous results
reported in Ref. [13] where the generation of entanglement in the environment
was studied. Here we study other measures of classical and quantum
correlations, such as the mutual information and the quantum discord. We also
present novel results and interpretations about the nature and behavior of the
entanglement generated in thermal machines.

The paper is organized as follows: In Sec. II we give a brief overview of the
model and its solution. It is a generalization of the standard quantum
Brownian motion model including a time-dependant driving field. In Sec. III we
use the previous results to compute the covariance matrix between different
environmental modes and we discuss its properties in the long-time regime. In
Sec. IV we present the measures we will use later to quantify and describe
classical and quantum correlations, including entanglement. In Sec. V we study
classical and quantum correlations between environmental bands at an arbitrary
environmental temperature. We provide simple analytical expressions to compute
the mutual information, the quantum discord and the entanglement in terms of
physically meaningful parameters. We summarize our results in Sec. VI.

##  II The model

We consider a generalization of the usual quantum Brownian motion model (QBM)
[17, 18] that was presented and solved in Ref. [13]. It consists of a
parametric oscillator, which we will refer to as system $\mathcal{S}$, coupled
with an environment $\mathcal{E}$ formed by $N$ independent oscillators. The
environment is then divided into two pieces $\mathcal{E}_{L}$ and
$\mathcal{E}_{R}$ by preparing them in different thermal states with
temperatures $T_{L}$ and $T_{R}$, respectively. This model represents the
physical situation shown in Fig. (1).

Figure 1: A parametric oscillator $\mathcal{S}$ driven by $V\left(t\right)$ is
coupled with two environments $\mathcal{E}_{L}$ and $\mathcal{E}_{R}$ at
temperatures $T_{L}$ and $T_{R}$, respectively.

Since the solution of this model has been treated in detail in [13], here we
will only present a brief outlook of it. The dynamics are governed by the
total Hamiltonian
$H_{T}=H_{\mathcal{S}}\otimes\mathbb{1}_{\mathcal{E}}+\mathbb{1}_{\mathcal{S}}\otimes
H_{\mathcal{E}}+H_{\mathcal{S},\mathcal{E}}$. The system’s Hamiltonian is
$H_{\mathcal{S}}=p^{2}/2m+mV(t)x^{2}/2$ while the environmental and
interaction terms are, respectively,
$H_{\mathcal{E}}=\sum_{i}\left(p_{i}^{2}/2m_{i}+m_{i}\omega_{i}^{2}q_{i}^{2}/2\right)$
and $H_{\mathcal{S},\mathcal{E}}=x\sum_{i}\lambda_{i}q_{i}$. The solution of
the Heisenberg equations of motion is

$\displaystyle q_{i}=q_{i}^{h}+K_{ij}^{\left(1\right)}\ast
q_{j}^{h}+K^{\left(2\right)}_{i}\ast x_{h},$ (1) $\displaystyle
x=x_{h}+K^{\left(3\right)}_{j}\ast q_{j}^{h},$

where the notation $F\ast
f=\int_{0}^{t}dt^{\prime}F\left(t,t^{\prime}\right)f\left(t^{\prime}\right)$
is used.
$q_{i}^{h}\left(t\right)=q_{i,0}\cos\left(\omega_{i}t\right)+p_{i,0}\sin\left(\omega_{i}t\right)/m_{i}\omega_{i}$
are the free Heisenberg operators of the environmental modes, with $q_{i,0}$
and $p_{i,0}$ Schröedinger operators, and $x_{h}\left(t\right)$ is a dressed
operator for $\mathcal{S}$ that satisfies the linear equation

$\ddot{x}_{h}+V_{R}\left(t\right)x_{h}+\gamma\ast\dot{x}_{h}=0.$ (2)

Above, $\gamma$ is the dissipation kernel $\gamma(t)=\int d\omega
I(\omega)\text{cos}(\omega t)/m\omega$, and
$V_{R}\left(t\right)=V\left(t\right)-\gamma\left(0\right)$ is the renormalized
potential. The explicit expressions for the kernels $K^{\left(1,2,3\right)}$
can be found in [13], but it is sufficient to know that they are functionals
of the Green function $G$ of Eq. (2). Thus, knowing $G$ means solving the full
Heisenberg equations. For a periodic driving with frequency $\omega_{d}$,
$V\left(t\right)=\sum_{k}V_{k}e^{ik\omega_{d}t}$, the Green function $G$ can
always be written as
$G\left(t,t^{\prime}\right)=\textstyle{\sum_{k}}A_{k}\left(t-t^{\prime}\right)e^{ik\omega_{d}t}$,
where $A_{k}$ vanishes for negative arguments. For $G$ to be a Green’s
function of Eq. (2), the Laplace transform of the functions $A_{k}(t)$,
$\tilde{A}_{k}\left(s\right)$, must satisfy a linear set of algebraic
equations:

$\tilde{g}^{-1}\left(s+ik\omega_{d}\right)\tilde{A}_{k}\left(s\right)+\textstyle{\sum_{n\neq
0}}V_{n}\,\tilde{A}_{k-n}\left(s\right)=\delta_{k0},$ (3)

where $\tilde{g}$ is the Laplace transform of the static Green function. The
above system of equations can be simply solved by using a perturbative series
expansion which is valid when the Fourier coefficients of the potential
$\lvert V_{k}\rvert$ are small and the frequency of the driving $\omega_{d}$
is detuned from the parametric resonance (that is, from the renormalized
frequency $\omega\textasciiacute_{r}=V_{0}-\gamma(0)$ [19]). In that case, the
solution of (3) satisfies the following recurrence relation:

$A_{k}^{\left(m\right)}\\!\left(s\right)=\tilde{g}\left(s+ik\omega_{d}\right)[\delta_{k0}-\textstyle{\sum_{n\neq
0}}V_{n}\,\tilde{A}_{k-n}^{\left(m-1\right)}\left(s\right)],$ (4)

for $m\geq 1$, with
$A_{k}^{\left(0\right)}\left(s\right)=\tilde{g}\left(s+ik\omega_{d}\right)\,\delta_{k0}$.

We will focus our attention on Gaussian initial states with first moments
equal to zero. Since the total Hamiltonian is quadratic, states remain
Gaussian for all times. Also, our results are independent of the specific form
of the spectral density $I_{R,L}(\omega)$ of the environments, but we will
require that it is a smooth function of $\omega$ with a frequency-independent
damping constant $\gamma_{0}$. That is, $I(\omega)=\gamma_{0}\,f(\omega)$ with
$f$ smooth.

##  III The covariance matrix

In order to make this paper self-contained, in this section we include some
material from our previous work [13]. We added new and better explanations to
aid the understanding of the broader topics treated here.

Our goal is to study environmental correlations between two parts of the
environment, so we use Eq. (1) to compute all correlation functions between
two environmental bands: one of them consists of oscillators whose frequencies
are distributed around $\omega_{i}\in\mathcal{E}_{R}$ with a bandwidth
$\Delta\omega$, and the other one is centered around
$\omega_{j}\in\mathcal{E}_{L}$. All expressions will depend on the product
$I(\omega)\Delta\omega$ which, for sufficiently small values of
$\Delta\omega$, plays the role of an effective coupling strength between the
reservoir band and $\mathcal{S}$ (since
$I(\omega)\Delta\omega\approx\lambda^{2}/m\omega$). With the correlation
functions we build the two-mode covariance matrix as
$\boldsymbol{\sigma}_{ab}(t)=\langle\\{z_{a}(t),z_{b}(t)\\}\rangle/2$ where
$\vec{z}=(q_{i},p_{i},q_{j},p_{j})$, which determines all properties of
Gaussian states. For the sake of simplicity we avoid writing all correlators
here (in Appendix A we include the general form of the position and momentum
correlators), but we can prove that, for a continuous environment (the
continuous hypothesis is discussed in Ref. [13]),

$\boldsymbol{\sigma}(t)=\boldsymbol{\nu}_{i}\otimes\boldsymbol{\nu}_{j}+\boldsymbol{\sigma}_{0}(t)+\boldsymbol{\sigma}_{osc}(t)+\boldsymbol{\sigma}_{lin}\times
t,$ (5)

where $\boldsymbol{\nu}_{i}\otimes\boldsymbol{\nu}_{j}$ is the initial thermal
state of both bands, $\boldsymbol{\sigma}_{0}(t)$ depends on the initial state
of $\mathcal{S}$ and $\boldsymbol{\sigma}_{0}(t)\to 0$ exponentially fast as
$t\to\infty$, $\boldsymbol{\sigma}_{osc}(t)$ oscillates in time with
frequencies $\omega_{d}$ and its higher harmonics, and
$\boldsymbol{\sigma}_{lin}\times t$ is linear in time. We note that the fact
that $\boldsymbol{\sigma}_{0}(t)\to 0$ is not trivial and rest on the
assumption that a stable stationary regime exists. The existence of such
regime requires the energy pumped into $\mathcal{S}$ to be dissipated by
$\mathcal{E}$. This can be achieved for small driving amplitudes provided that
$\omega_{d}$ is detuned from the parametric resonance.

We are interested in the physical processes that continuously create
correlations beyond the transient regime. Therefore, we will only keep the
terms that do not equal zero after averaging $\boldsymbol{\sigma}$ over a
driving cycle in the long-time limit. Thus, we will work with

$\boldsymbol{\sigma}_{av}(t)=\boldsymbol{\nu}_{i}\otimes\boldsymbol{\nu}_{j}+\boldsymbol{\sigma}_{lin}\times
t.$ (6)

The linear term in time, $\boldsymbol{\sigma}_{lin}$, is a block diagonal
matrix (that is, the cross-correlators are zero) unless the center frequencies
of the bands satisfy a precise relation: either
$\omega_{i}+\omega_{j}=k\omega_{d}$ or $\omega_{j}=\omega_{i}+k\omega_{d}$,
with $k\in\mathbb{Z}$. This is the signature of the two main processes that
creates correlations between the bands. In order to understand their nature,
we compute the energy $E_{i}$ stored in the band $i\in\mathcal{E}_{R}$, which
varies due to the same two processes. $E_{i}$ is proportional to the sum of
the two first diagonal correlators of $\boldsymbol{\sigma}_{av}$
($E_{i}\propto\boldsymbol{\sigma}_{av,11}+\boldsymbol{\sigma}_{av,22}$, see
Appendix A for a brief derivation) and has a simple expression:
$E_{i}(t)=\left[1/2+n_{R}\left(\omega_{i}\right)\right]\omega_{i}+\dot{\mathcal{Q}}_{i}\times
t$, with

$\displaystyle\frac{\dot{\mathcal{Q}}_{i}}{\Delta\omega}=\sum_{k}\sum_{\alpha=R,L}\omega_{i}\left\\{\Theta\left(\omega_{i,k}\right)p_{R,\alpha}^{(k)}\left(\omega_{i}\right)\left[n_{\alpha}\left(\omega_{i,k}\right)-n_{R}\left(\omega_{i}\right)\right]\right.$
(7)
$\displaystyle\left.+\Theta\left(-\omega_{i,k}\right)p_{R,\alpha}^{(k)}\left(\omega_{i}\right)\left[n_{\alpha}\left(\left\lvert\omega_{i,k}\right\rvert\right)+n_{R}\left(\omega_{i}\right)+1\right]\right\\},$

where $\omega_{i,k}=\omega_{i}-k\omega_{d}$, $\Theta$ is the step function,
$n_{\alpha}\left(\omega\right)$ is the Planck distribution with temperature
$T_{\alpha}$, and
$p_{R,\alpha}^{(k)}\left(\omega_{i}\right)=\pi\,I_{R}\left(\omega_{i}\right)I_{\alpha}\left(\left\lvert\omega_{i,k}\right\rvert\right)\lvert\tilde{A}_{k}\left(i\omega_{i,k}\right)\rvert^{2}/2m^{2}$
is the probability that mode $\omega_{i}$ in $\mathcal{E}_{R}$ interacts with
mode $\omega_{i,k}$ in $\mathcal{E}_{\alpha}$ through $\mathcal{S}$
($I_{\alpha}$ is the spectral density of $\mathcal{E}_{\alpha}$). In Eq. (7)
we can see the two main processes that generate correlations between the
bands. In the long-time regime the energy stored in each band varies due to
these two processes and their relative importance depends on the temperatures
$T_{R}$ and $T_{L}$. As the first term on the right-hand side of Eq. (7)
shows, at nonzero environmental temperature, the resonant absorption (or
emission) from (or into) the driving field can transport an excitation in a
mode with frequency $\omega_{i}$ to a mode with frequency
$\omega_{i}-k\omega_{d}$. This process is associated with the classical heat
current that flows from one environment to another and, naturally, it is not
present at zero temperature where there are no excitations to be transported
around. Thus, at very low temperature, $E_{i}$ varies because of a different
process. In the second term of Eq. (7), we can see that the energy of the
driving is dumped into two modes whose frequencies add up to a multiple of
$\omega_{d}$. This splitting of the driving energy between two modes generates
entanglement between them, and is interpreted as the nonresonant creation of a
pair of excitations (one in each mode) [13, 20, 21]. This second term is
positive and therefore it is always associated with heating—more specifically,
with the third law of thermodynamics.

We have seen that the only time-extensive correlations are generated between
environmental bands that satisfy either $\omega_{i}+\omega_{j}=k\omega_{d}$ or
$\omega_{j}=\omega_{i}+k\omega_{d}$. The former are correlated by the
nonresonant creation of a pair of excitations, while the latter are correlated
by the resonant transport of excitations. In following sections we will study
the nature of the correlations generated by these two processes and their
relation with entanglement. Before that, we note that in the time-extensive
regime, $\boldsymbol{\sigma}_{av}$ is symplectically equivalent to the
standard form (see Lemma I in Ref. [22])

$\boldsymbol{\sigma}_{sf}=\begin{pmatrix}a&0&c_{1}&0\\\ 0&a&0&c_{2}\\\
c_{1}&0&b&0\\\ 0&c_{2}&0&b\end{pmatrix}$ (8)

where $c_{1}=c_{2}$ for the bands correlated by the resonant transport of
excitations and $c_{1}=-c_{2}$ for the ones correlated by the nonresonant
creation of a pair of excitations. This is the subclass of squeezed-thermal
states. Therefore the driving might be creating squeezed thermal states
between the bands in the environment.

##  IV The toolbox

In order to study the cause and nature of intraenvironmental correlations in
thermal machines, we need:

  1. (i)

a measure of the total amount of correlations,

  2. (ii)

a measure of quantum correlations, and

  3. (iii)

a measure of entanglement.

For the first item, the immediate choice is the mutual information
$\mathcal{I}$. For a bipartite system divided in $A$ and $B$, it is defined as
$\mathcal{I}(A:B)=S(A)+S(B)-S(AB)$, where $S$ is the von Neumann entropy. It
is worth noticing that, whereas the mutual information for classical
distributions is bounded from above by the entropy of either one of them, the
mutual information for quantum systems can be as large as
$2\,\text{min}\\{S(A),S(B)\\}$ for entangled states (see Araki-Lieb
inequality—Corollary of Theorem 2 in Ref. [23]), which reflects the existence
of quantum correlations beyond the classical ones. In order to quantify these
quantum correlations, we will use a measure called quantum discord (QD). QD
was proposed [24, 25] in an effor to adress a series of developments in Refs.
[26, 27, 28] that challenged the belief that entanglement is the only form of
quantum correlation. QD is defined by quantizing the mismatch between two
classically equivalent measures of mutual information:
$\overleftarrow{D}(A:B)=\mathcal{I}(A:B)-\mathcal{J}(A:B)$, where
$\mathcal{I}(A:B)=H(A)+H(B)-H(A,B)$ and $\mathcal{J}(A:B)=H(A)-H(A\rvert B)$,
with $H$ the Shannon entropy. The discrepancy between $\mathcal{I}$ and
$\mathcal{J}$ in the quantum case occurs because, in quantum theory, the
conditional entropy $H(A\rvert B)$ involves an specific choice of basis to
perform a measurement on $B$ to infer the state of $A$. Thus, the quantization
of $\mathcal{J}$ is not as straightforward as the one of $\mathcal{I}$, which
only requires to replace the Shannon entropy for the von Neumann one. In order
to do this, a minimization is performed over all positive operator valued
measures (POVMs) corresponding to a measurement in $B$ so as to find the one
that disturbs least the overall quantum state and that, at the same time,
allows to extract the most information about $B$. Hence, the quantum analog of
$\mathcal{J}$ is defined as

$\overleftarrow{\mathcal{J}}(A:B)=S(A)-\inf_{\\{\Pi_{n}\\}}\left[\sum_{n}\,p_{n}\,S(A|n)\right],$
(9)

where $\\{\Pi_{n}\\}$ is a POVM, $p_{n}=\text{tr}(\rho_{B}\Pi_{n})$ is the
probability of obtaining the result $n$, and $S(A|n)$ is the von Neumann
entropy of the reduced state of $A$ after obtaining this result. Quantum
discord has been shown to be a property held by almost all quantum states [29]
(the set of states with zero quantum discord has measure zero and is nowhere
dense) and has attracted considerable attention [30, 31, 32, 33, 34].
Initially defined only in finite dimensional systems, the concept of QD was
extended to continuous-variable systems, specifically to the case of two-mode
Gaussian states [35, 36]. When the POVMs are restricted to the set of Gaussian
measurements, it is called Gaussian quantum discord (GQD). Until recently it
was thought to be an upper bound of the QD for continuous-variable systems but
it has been proved that, in fact, QD and GQD are equal for Gaussian states
[37]. Finally, as a measure of entanglement we will use the logarithmic
negativity $E_{\mathcal{N}}=\text{log}\|\boldsymbol{\tilde{\rho}}\|_{1}$,
where $\tilde{\rho}$ is the partially transposed density matrix.
$E_{\mathcal{N}}$ quantifies the violation of the Peres-Horodecki criterion
and it is monotone under local operations and classical communication [38]. In
the case of two-mode Gaussian states, this criterion is a necessary and
sufficient condition for the composite system to be separable.

Gaussian states are completely determined by their covariance matrix
$\boldsymbol{\sigma}$, so it is not surprising that all three quantities
discussed above can be computed using only $\boldsymbol{\sigma}$. For example,
the logarithmic negativity is computed as
$E_{\mathcal{N}}=\text{max}\\{0,-\text{ln}(2\tilde{\lambda}_{-})\\}$, where
$\tilde{\lambda}_{-}$ is the lowest symplectic eigenvalue of the covariance
matrix $\boldsymbol{\tilde{\sigma}}$ corresponding to the partially transposed
density matrix $\boldsymbol{\tilde{\rho}}$ (which differs from
$\boldsymbol{\sigma}$ by just a sign flip in the off-diagonal block matrices).
Expressions for the mutual information and quantum discord in terms of
$\boldsymbol{\sigma}$ can be respectively found in Ref. [39] and in Ref. [36].

##  V Correlations between environmental bands

In this section we present a study of classical and quantum correlations
between two environmental bands which are respectively centered around
frequencies $\omega_{i}\in\mathcal{E}_{R}$ and $\omega_{j}\in\mathcal{E}_{L}$,
in increasing order of quantumness. We will use the toolbox described in the
previous section. First, we will compute the mutual information $\mathcal{I}$
for the total correlations, second, the quantum discord $\overleftarrow{D}$,
which measures the quantumness of correlations, and third, as a measure of
entanglement, the logarithmic negativity $E_{\mathcal{N}}$. We will analyze
the relationship between these quantities and present simple analytic
expressions which are valid in the limit of small driving and weak coupling
($\gamma_{0}/\omega_{i,j}\ll 1$). We will examine and compare two distinct
physical processes which, as discussed above, are the main sources of
generation of time-extensive correlations in the long-time regime: the non-
reasonant creation of a pair of excitations, satisfying the condition
$\omega_{i}+\omega_{j}=k\omega_{d}$, and the resonant transport of
excitations, with $\omega_{j}=\omega_{i}+k\omega_{d}$.

Our formulas are written in terms of only three physically meaningful
symplectic invariants. These are: the individual purities
$\mu_{i}=\text{tr}(\rho_{\alpha}^{2})=1/2\sqrt{\text{det}(\boldsymbol{\alpha})}$
and $\mu_{j}=1/2\sqrt{\text{det}(\boldsymbol{\beta})}$, and the determinant
$\Gamma(t)=4\,\text{det}(\boldsymbol{\gamma})$, which is a measure of the
generation of cross correlations (see below). Here, $\boldsymbol{\alpha}$,
$\boldsymbol{\beta}$ and $\boldsymbol{\gamma}$ are the three $2\times 2$
submatrices in the two-mode covariance block matrix $\boldsymbol{\sigma}$. If
we write $\boldsymbol{\sigma}$ in its standard form as in Eq. (8), then
$\boldsymbol{\alpha}=a\,\mathbb{1}$, $\boldsymbol{\beta}=b\,\mathbb{1}$, and
$\boldsymbol{\gamma}=\text{diag}(c_{1},c_{2})$.

Although the individual purities depend on time (they decrease as time goes
by), in the weak-coupling limit their value is mostly determined by the
temperature of the environments. That is, $\mu_{i,j}\to 1^{-}$ indicates low
temperature whereas $\mu_{i,j}\to 0^{+}$ indicates high temperature. On the
other hand, $\Gamma(t)$ is a quadratic function of time,
$\Gamma_{\pm}(t)=\Gamma_{\pm}\times t^{2}$ (from now on, the plus sign will
indicate the nonresonant case and the minus sign, the resonant one). As we
mentioned above, $\Gamma_{\pm}(t)$ can be interpreted as the generator of
cross-correlations: pair creation or transport of excitations, depending on
the case. This is because the Kullback-Leibler divergence between the Wigner
function of $\rho_{\boldsymbol{\sigma}}$ and the Wigner function of the
product of the marginals
$\rho_{\boldsymbol{\alpha}}\otimes\rho_{\boldsymbol{\beta}}$ can be written as

$D_{KL}(W_{\rho_{\boldsymbol{\sigma}}}||W_{\rho_{\boldsymbol{\alpha}}\otimes\rho_{\boldsymbol{\beta}}})\simeq\mu_{i}\,\mu_{j}\,\lvert\Gamma_{\pm}(t)\rvert$
(10)

(see Ref. [40]). $D_{KL}$ measures the phase-space distinguishability between
$\rho_{\boldsymbol{\sigma}}$ and
$\rho_{\boldsymbol{\alpha}}\otimes\rho_{\boldsymbol{\beta}}$. As we can see
from Eq. (10), while the purities decrease with time,
$\lvert\Gamma_{\pm}(t)\rvert$ increases quadratically making
$\rho_{\boldsymbol{\sigma}}$ more and more distinguishable from the completely
uncorrelated state
$\rho_{\boldsymbol{\alpha}}\otimes\rho_{\boldsymbol{\beta}}$. Thus,
$\Gamma_{\pm}(t)$ is indeed responsible for generating correlations.

In Appendix B we include explicit expressions for the individual purities in
Eq. (37) and for $\Gamma_{\pm}$ in Eqs. (38) and (39), respectively. In
Appendix C we include a brief derivation of the equations presented below.

###  V.1 Mutual information

We begin our study of correlations by computing the mutual information between
the bands. For both the resonant and nonresonant cases, $\mathcal{I}$ has the
same simple form

$\mathcal{I}_{\pm}(\boldsymbol{\sigma}_{av})\simeq
f_{\pm}(\mu_{i},\mu_{j})\,\lvert\Gamma_{\pm}(t)\rvert$ (11)

with

$f_{\pm}(\mu_{i},\mu_{j})=\mu_{i}\mu_{j}[\text{atanh}(\mu_{i})\pm\text{atanh}(\mu_{j})]/(\mu_{i}\pm\mu_{j}).$
(12)

Although Eq. (11) may give the impression that mutual information grows
quadratically in time like $\,\lvert\Gamma_{\pm}(t)\rvert$, this is not the
case. The reason is that, while $\,\lvert\Gamma_{\pm}(t)\rvert$ increases,
$f_{\pm}$ is a monotonically decreasing function of the individual purities.
In fact, for the resonant case it can be shown (see Appendix B) that

$\lvert\Gamma_{-}(t)\rvert\leq(1-\mu_{i})(1-\mu_{j})/\mu_{i}\mu_{j},$ (13)

and, consequently, the mutual information is bounded: $\mathcal{I}_{-}\leq 1$.
Thus, the quadratic growth in time must eventually stop and $\mathcal{I}_{-}$
reaches a saturation value. It seems that the generation of correlations by
means of the transport of excitations cannot increase the mutual information
between the bands to arbitrarily large values. On the contrary, in the
nonresonant case the bound in Eq. (13) does not hold. This is a natural result
taking into account that this process continually creates correlations in
time. However, it exists a regime in which the mutual information
$\mathcal{I}_{\pm}$ grows relatively quadratically in time. For not so long-
times such that the individual purities $\mu_{i,j}$ stay approximately
constant, $f_{\pm}$ stays constant too and the mutual information grows as
$\,\lvert\Gamma_{\pm}(t)\rvert$.

With this in mind, Eq. (11) has a clear physical interpretation: mutual
information increases due to the creation of cross-correlations between the
bands by $\Gamma_{\pm}(t)$, but this increment is countered by the influence
of the thermal agitations in the environment, which are represented by
$f_{\pm}$. Certainly, as time passes or temperature increases the bands
correlate progressively more with the rest of the environment, decreasing its
mutual dependence. Furthermore, since $f_{+}<f_{-}$, the effects of the
thermal agitations in the environment are more detrimental to the bands that
are related by the nonresonant pair creation than to the ones that are related
by the more “classical” resonant transport of excitations. Mutual information
is plotted in Figs. 2 and 3 for different environmental temperatures, in the
nonresonant and resonant cases, respectively.

An interesting feature can be seen from the form of $\Gamma_{-}$:

$\Gamma_{-}\propto\lvert
n_{R}(\omega_{i})\tilde{A}_{k}(i\omega_{i})-n_{L}(\omega_{j})\tilde{A}_{-k}^{\ast}(i\omega_{j})\rvert^{2}$
(14)

(the complete expression can be found in Eq. (39) in Appendix B). In the
resonant case, the transport of excitations can go both ways depending if a
quantum of energy is either absorbed or emitted from the driving field.
Because excitations act as identical quasiparticles, an interference effect
takes place and the correlations generated by the absorption tend to cancel
the ones generated by the emission. Indeed, the case where no excitations have
been transported is indistinguishable from the one where the same number of
excitations have been transferred in both directions. In particular, when the
bands have the same occupation number (i.e.,
$n_{R}(\omega_{i})=n_{L}(\omega_{j})$), absorption and emission are equally
likely and mutual information tends to vanish. This is reflected in the fact
that $\Gamma_{-}$ becomes vanishingly small, as can be seen in Eq. (14). This
is illustrated in Fig. 3(c) (note the scale of the vertical axis).

Seeing that the pair creation mechanism is the dominant process at lower
temperatures, physical intuition tells us that the mutual information between
the bands correlated by it must be bigger than the ones correlated by the
resonant transport of excitations. This is confirmed by an analysis of
$\Gamma_{\pm}$. As opposed to the nonresonant case, in the resonant case,
$\Gamma_{-}\to 0$ as $T\to 0$ (see Eq. (14)). At zero environmental
temperature the bands are not populated and, therefore, the driving cannot
transport excitations from one into the other.

Figure 2: Dependence of the mutual information, quantum discord and
logarithmic negativity on the frequency $\omega_{i}$ for the nonresonant case,
with $\omega_{i}+\omega_{j}=\omega_{d}$. (a) When both environments are at
zero temperature, entanglement is maximized and
$\overleftarrow{D}_{+}/\mathcal{I}_{+}\to 1/2$. (b,c) Temperature is raised in
one of the environments to $T=7.5\gamma_{0}$, while the other one is kept at
$T=0$. This illustrates the asymmetry of the quantum discord. (d) Both
temperatures are raised to $T=7.5\gamma_{0}$ and the only entanglement left is
present at the center frequencies. Here,
$V(t)=\omega_{r}^{2}+V\text{cos}(\omega_{d}t)$ and
$I(\omega)=2m\gamma_{0}\omega\Lambda^{2}/\pi(\omega^{2}+\Lambda^{2})$. The
plots are normalized using $E_{0}=\gamma_{0}\Delta\omega Vt/\omega_{r}^{3}$
(we plot $\mathcal{I}_{+}/E_{0}^{2}$, $\overleftarrow{D}_{+}/E_{0}^{2}$ and
$E_{\mathcal{N}}/E_{0}$). The parameters used are
$\omega_{d}=\omega_{r}/\sqrt{11}$, $\omega_{r}=800\gamma_{0}$,
$V=\omega_{r}^{2}/32$, $m=10m_{i}$, $t=20\gamma_{0}$ and $\gamma_{0}=0.005$.
Figure 3: Dependence of the mutual information and quantum discord on the
frequency $\omega_{i}$ for the resonant case, with
$\omega_{j}=\omega_{i}+\omega_{d}$. (a,b) Temperature is raised in one of the
environments to $T=1500\gamma_{0}$, while the other one is kept at $T=0$. This
illustrates the limits shown in Eqs. (19) and (20), respectively. (c) Both
temperatures are raised to $T=7500\gamma_{0}$, showing that $\mathcal{I}_{-}$
becomes vanishingly small when occupation numbers are similar (note that the
scale of the vertical axis is 100 times smaller than the previous plots).
Here, $V(t)=\omega_{r}^{2}+V\text{cos}(\omega_{d}t)$ and
$I(\omega)=2m\gamma_{0}\omega\Lambda^{2}/\pi(\omega^{2}+\Lambda^{2})$. The
plots are normalized using $E_{0}=\gamma_{0}\Delta\omega Vt/\omega_{r}^{3}$
(we plot $\mathcal{I}_{-}/E_{0}^{2}$ and $\overleftarrow{D}_{-}/E_{0}^{2}$).
The parameters used are $\omega_{d}=\omega_{r}/\sqrt{11}$,
$\omega_{r}=800\gamma_{0}$, $V=\omega_{r}^{2}/32$, $m=10m_{i}$,
$t=20\gamma_{0}$ and $\gamma_{0}=0.005$.

###  V.2 Quantum discord

In this section we will compute the ratio $\overleftarrow{D}/\mathcal{I}$
which represents the amount of quantum correlations among the totality of
correlations. The result is

$\overleftarrow{D}_{\pm}(\boldsymbol{\sigma}_{av})/\mathcal{I}_{\pm}(\boldsymbol{\sigma}_{av})\simeq
1-g_{\pm}(\mu_{i},\mu_{j}),$ (15)

with

$g_{\pm}(\mu_{i},\mu_{j})=\frac{1}{1+\mu_{j}}\left(1\pm\frac{\mu_{j}}{\mu_{i}}\right)\left[1\pm\frac{\text{atanh}(\mu_{j})}{\text{atanh}(\mu_{i})}\right]^{-1}.$
(16)

We note that $g_{\pm}$ is such that $0\leq g_{\pm}\leq 1$ and it is a
monotonically increasing function of the purities. Since $g_{+}>g_{-}$, it
follows that
$\overleftarrow{D}_{+}/\mathcal{I}_{+}<\overleftarrow{D}_{-}/\mathcal{I}_{-}$.
That is, the fraction of quantum correlations in the mutual information
between the bands correlated by the nonresonant pair creation is less than the
fraction between the bands correlated by the resonant transport of
excitations. Nevertheless, it is possible for the quantum discord in the
nonresonant case to be greater than the one in the resonant one (i.e.,
$\overleftarrow{D}_{+}>\overleftarrow{D}_{-}$). For example, that is the case
in the low temperature limit where $\mathcal{I}_{-}\to 0$ but
$\mathcal{I}_{+}\nrightarrow 0$. As a consequence of the mutual information in
the resonant case being bounded, the quantum discord is bounded too:
$\overleftarrow{D}_{-}\leq 1$. This result was previously obtained in Refs.
[36, 35] when studying quantum discord in squeezed thermal states. Thus, it is
consistent with the idea that that is the type of state generated by the
resonant process. Now we will provide some examples that illustrate the
previous results in two different regimes.

We will analyze the case where environmental temperatures are similar (a more
precise statement would be $\mu_{i}\simeq\mu_{j}$, but in the limits of high
and low temperatures both conditions are equivalent). We will express all
quantities in terms of $\bar{\mu}=(\mu_{i}+\mu_{j})/2$ and
$\Delta\mu=(\mu_{i}-\mu_{j})/2$, and we will work up to first order in
$\lvert\Delta\mu\rvert/\bar{\mu}\ll 1$. In this case, the ratio
$\overleftarrow{D}/\mathcal{I}$ is

$\overleftarrow{D}_{+}(\boldsymbol{\sigma}_{av})/\mathcal{I}_{+}(\boldsymbol{\sigma}_{av})\simeq\frac{\bar{\mu}}{1+\bar{\mu}}$
(17)

and

$\overleftarrow{D}_{-}(\boldsymbol{\sigma}_{av})/\mathcal{I}_{-}(\boldsymbol{\sigma}_{av})\simeq
1-\frac{1}{\bar{\mu}}(1-\bar{\mu})\,\text{atanh}(\bar{\mu}).$ (18)

As we can see from Eqs. (17) and (18) above, at high temperatures, when
$\bar{\mu}\to 0^{+}$, both expressions go to zero (i.e.,
$\overleftarrow{D}/\mathcal{I}\to 0$). In the other limit, when $\bar{\mu}\to
1^{-}$, we have $\overleftarrow{D}_{+}/\mathcal{I}_{+}\to 1/2$ and
$\overleftarrow{D}_{-}/\mathcal{I}_{-}\to 1$. For the nonresonant case, it
means that half of the total correlations are quantum in origin. This can be
interpreted as follows. For the quantum system, mutual information is bounded
by $\mathcal{I}_{+}\leq 2\,S_{\bar{\mu}}$ (Araki-Lieb inequality with
$S_{\bar{\mu}}$ the entropy of an “average band” represented by $\bar{\mu}$),
while for a classical system the mutual information satisfies
$\mathcal{I}_{+}\leq S_{\bar{\mu}}$. Therefore, the above limit
($\overleftarrow{D}_{+}/\mathcal{I}_{+}\to 1/2$) seems to indicate that the
pair creation mechanism at zero temperature saturates the Araki-Lieb
inequality. Thus, $\mathcal{I}_{+}\to 2\,S_{\bar{\mu}}$, and, qualitatively,
$0<\mathcal{I}_{+}\leq S_{\bar{\mu}}$ corresponds to classical correlations
and $S_{\bar{\mu}}<\mathcal{I}_{+}\leq 2\,S_{\bar{\mu}}$ to quantum ones. This
is shown in Fig. 2 (a). On the other hand, in the resonant case the situation
is different. When $\bar{\mu}\to 1^{-}$, the mechanism creating correlations
tends to vanish (remember that in this case,
$\mathcal{I}_{-}\ll\mathcal{I}_{+}$) but the small amount of correlations that
are present will have nonzero quantum discord and satisfy
$\overleftarrow{D}_{-}\simeq\mathcal{I}_{-}$. As time passes or temperature
increases, classical correlations start to develop and the quantum discord
begins to decrease.

Since the definition of the quantum discord is based on measurements over one
of the subsystems, it is naturally asymmetric with respect to an
$\\{i,R\\}\leftrightarrow\\{j,L\\}$ interchange. It is interesting to explore
how this asymmetry manifests itself in terms of the temperatures of the
environments. For example, let us use the resonant case (the nonresonant one
is analog), and suppose one of the environments is at low temperature and the
other one is at an arbitrary (but higher) one. If $\mathcal{E}_{L}$ is the one
at low temperature, then Eq. (15) can be approximated as

$\overleftarrow{D}_{-}(\boldsymbol{\sigma}_{av})/\mathcal{I}_{-}(\boldsymbol{\sigma}_{av})\simeq
1-\frac{\mu_{j}}{2\,\text{atanh}(\mu_{j})}.$ (19)

Equation (19) shows that most correlations are quantum ones in that limit. In
contrast, if $\mathcal{E}_{R}$ is the one at low temperature, we have

$\overleftarrow{D}_{-}(\boldsymbol{\sigma}_{av})/\mathcal{I}_{-}(\boldsymbol{\sigma}_{av})\simeq
2\,\frac{\mu_{j}}{\mu_{i}}.$ (20)

Now if we increase the temperature of $\mathcal{E}_{L}$, quantum discord
reaches its minimum value: all correlations are classical. As we can see from
the results above, the quantumness of the correlations measured by the quantum
discord highly depends on which system is being observed. This difference is
shown in Figs. 2 (b) and (c) for the nonresonant case, and Figs. 3 (a) and (b)
for the resonant one (in this case we explicitly show the limits of Eqs. (19)
and (20)).

###  V.3 Entanglement

The generation of entanglement in the nonresonant case was addressed in our
previous work [13]. Here we rederive its main equations but this time written
in terms of symplectic invariants. We expand on its interpretation and the
connection between entanglement, the phase-space entropy and the generation of
correlations.

In the resonant case we have $\text{det}(\boldsymbol{\gamma})>0$, and
therefore the logarithmic negativity is zero. That is, bands with center
frequencies such that $\omega_{j}=\omega_{i}+k\omega_{d}$ are not entangled.
Therefore, we focus our attention on the bands correlated by the pair creation
mechanism (i.e, such that $\omega_{i}+\omega_{j}=k\omega_{d}$). We can show
that these bands are entangled and the logarithmic negativity
$E_{\mathcal{N}}$ is the maximum between zero and

$\text{E}_{\mathcal{N}}(t)\simeq-S_{ij}+\Gamma_{\mathcal{N}}\times t$ (21)

where
$\Gamma_{\mathcal{N}}=(\mu_{i}+\mu_{j})e^{-2S_{ij}}\sqrt{\lvert\Gamma_{+}\rvert}/2\mu_{i}\mu_{j}$
and
$S_{ij}=\text{ln}\left[\left(\mu_{i}^{2}+\mu_{j}^{2}\right)/2\mu_{i}^{2}\mu_{j}^{2}\right]/2\geq
0$. We note that, just like it happens with the mutual information and the
quantum discord, $\text{E}_{\mathcal{N}}$ is not linear in time as Eq. (21)
may suggest. Entanglement production is time-extensive as long as the
individual purities stay approximately constant. From Eq. (21) we can see that
the creation of entanglement is a competition between the generation of
correlations by the pair creation mechanism and the thermal agitations of the
environment that try to destroy them. These thermal agitations have two
distinct effects on the logarithmic negativity. First, they reduce the rate of
generation $\Gamma_{\mathcal{N}}$
($\Gamma_{\mathcal{N}}\leq\sqrt{\lvert\Gamma_{+}\rvert}$).
$\Gamma_{\mathcal{N}}$ is a function of the individual purities and it reaches
its maximum value when $\mu_{i,j}\to 1^{-}$
($\Gamma_{\mathcal{\mathcal{N}}}\to\sqrt{\lvert\Gamma_{+}\rvert}$). And
second, they provide a threshold $S_{ij}$ that the pair creation mechanism has
to overcome to entangle the bands. Logarithmic negativity is plotted alongside
mutual information and quantum discord in Figs. 2 (a)-(d) for different
environmental temperatures. Note that, as temperature is increased,
entanglement disappears from the side that corresponds to the lower
frequencies in the higher temperature environment.

It is worth noticing that $S_{ij}$ is obtained from the Shannon entropy of the
phase-space Wigner distributions corresponding to $\rho_{\boldsymbol{\alpha}}$
and $\rho_{\boldsymbol{\beta}}$. Thus, it is a measure of the disorder of the
bands. In fact, it can be written as

$S_{ij}=\bar{S}_{2}+\text{ln}\left[\text{cosh}\left(\Delta
S_{2}\right)\right]/2,$ (22)

where
$\bar{S}_{2}=\left[S_{2}(\boldsymbol{\alpha})+S_{2}(\boldsymbol{\beta})\right]/2$
and $\Delta S_{2}=S_{2}(\boldsymbol{\alpha})-S_{2}(\boldsymbol{\beta})$, with
$S_{2}(\rho)=-\text{ln}[\text{tr}(\rho^{2})]$ the Renyi-2 entropy. For
Gaussian states, $S_{2}$ coincides with the Shannon entropy of the
corresponding phase-space Wigner distribution (up to an additive constant)
[40]. As Eq. (22) shows, not only the entropy of the bands affects the
generation of entanglement, but their difference too. For example, since for
high temperatures $S_{2}\sim\text{ln}(T)$, it is harder to entangle two bands
at different temperatures than two at the same one. $S_{ij}$ in turn fixes a
latency time $t_{ent}=S_{ij}/\Gamma_{\mathcal{N}}$, which is the time it takes
for the bands to be entangled.

As we previously mentioned, the generator of intra environmental correlations
is $\Gamma_{\pm}(t)$ (see Eq. (10)). In view of the above results, in the
nonresonant case, we can rewrite this generator as

$\lvert\Gamma_{+}(t)\rvert=\frac{\lvert\Gamma_{+}\rvert}{\Gamma_{\mathcal{N}}^{2}}[E_{\mathcal{N}}(t)+S_{ij}]^{2}\qquad
t>t_{ent}.$ (23)

Equation (23) shows that, from the moment the bands are entangled onwards, the
pair creation mechanism can be interpreted as having the dual effect of
increasing that entanglement and generating entropy in the environment. It
must be emphasized that, since $S_{ij}$ is mostly constant through time, the
predominant effect is the production of entanglement. As a consequence of Eq.
(23), the mutual information in Eq. (11) and the quantum discord in Eq. (15)
can be cast in terms of the logarithmic negativity and, for fixed
$E_{\mathcal{N}}$, both are nonmonotonic functions of the purities, in
accordance with previous results [35]. We also note that the second term in
Eq. (21), which represents the generation of entanglement, is proportional to
the number of entangled pairs $\Delta\omega$, as it should be since the
logarithmic negativity is additive. This is in opposition to the mutual
information in Eq. (11) (and the quantum discord), which is proportional to
$\Delta\omega^{2}$ and not additive.

As we mentioned in the beginning of this section, bands such that
$\omega_{j}=\omega_{i}+k\omega_{d}$ are not entangled. That is, the resonant
process cannot create entanglement. Nevertheless, it may play a role in the
entanglement present between bands correlated by the nonresonant process. Let
us suppose the excitation in band $i$ that the driving absorbs to be
transported to band $j$ was previously put there by the pair creation
mechanism. In that case, there must exist a $j^{\prime}$ band such that
$\omega_{i}+\omega_{j^{\prime}}=k^{\prime}\omega_{d}$. Thus, the $j$ and
$j^{\prime}$ bands are now related by
$\omega_{j}+\omega_{j^{\prime}}=(k+k^{\prime})\omega_{d}$, meaning that they
are, in fact, entangled. From this we conclude that it is possible for the
resonant process to swap the entanglement created by the nonresonant process
to other bands. This should be a higher order effect in the damping constant
$\gamma_{0}$ and thus not visible in our expressions above.

##  VI Conclusions

In this paper we presented a complete analysis of the origin and nature of the
time-extensive correlations that are present in the reservoirs of a thermal
machine (which we divided in the so-called bands). By studying a
generalization of the usual QBM model that includes a time-dependent driving
enforced on the system, we were able to show that there are only two processes
that generate intraenvironmental correlations beyond the transient regime: the
resonant transport of excitations, responsible for the classical heat flow,
and the nonresonant pair creation, linked to the third law of thermodynamics.
We would like to stress the most important results:

  1. (i)

There is a regime where mutual information and quantum discord between the
bands grow approximately quadratically in time. This happens as long as the
individual purities stay approximately constant. Furthermore, in the resonant
case (with center frequencies such that $\omega_{j}=\omega_{i}+k\omega_{d}$)
both quantities are bounded from above (showing that this regime must
eventually end).

  2. (ii)

Quantum discord is always present between the bands, independently of the
process that correlates them. As shown in Eq. (15), the fraction of quantum
correlations is negatively affected by the thermal agitations of the
environment. This effect is worse on the bands correlated by the nonresonant
process than the ones correlated by the resonant one.

  3. (iii)

Entanglement is only present on the bands that are connected by the
nonresonant process (with center frequencies such that
$\omega_{i}+\omega_{j}=k\omega_{d}$). We showed that the logarithmic
negativity grows approximately linear in time and that the pair creation
mechanism has to overcome a threshold $S_{ij}$ to entangle the bands that is a
measure of the disorder of them (see Eq. (22)).

  4. (iv)

The resonant process cannot create entanglement — at most it can swap it.
Since it can only transport excitations between bands, if the band on which
the excitation is absorbed was already entangled with another one, it can swap
that entanglement to the band on which the excitation is dumped. This is a
higher order effect which was not included in our formulas above.

Our results are consistent with previous studies of classical and quantum
correlations in Gaussian states [35, 36]. Using the inequality presented in
Eq. (13) we showed that the quantum discord between the bands correlated by
the resonant process is bounded from above. Additionally, using Eq. (23), we
showed that the mutual information (and the quantum discord) between the ones
correlated by the nonresonant process written in terms of the entanglement is
a nonmonotonic function of the purities.

##  Appendix A Position and momentum correlation functions

Here we present the general expression for the position and momentum
correlation functions for two modes $i\in\mathcal{E}_{R}$ and
$j\in\mathcal{E}_{L}$. In order to obtain the correlation function for just
one mode (e.g.
$\left\langle\left\\{q_{i}\left(t\right),q_{i}\left(t\right)\right\\}\right\rangle$
instead of
$\left\langle\left\\{q_{i}\left(t\right),q_{j}\left(t\right)\right\\}\right\rangle$)
we just need to make the replacement
$\left\\{i,R\right\\}\rightarrow\left\\{j,L\right\\}$. The position correlator
at an arbitrary environmental temperature is

$\displaystyle\left\langle\left\\{\hat{q}_{i}\left(t\right),\hat{q}_{j}\left(t\right)\right\\}\right\rangle$
$\displaystyle=\frac{1}{m_{i}\omega_{i}}\,\left[2\,n_{R}\left(\omega_{i}\right)+1\right]\,\delta_{ij}$
(24)
$\displaystyle+\frac{1}{2m}\frac{1}{\sqrt{m_{i}\omega_{i}}}\frac{1}{\sqrt{m_{j}\omega_{j}}}\Delta\omega\sqrt{I_{R}\left(\omega_{i}\right)I_{L}\left(\omega_{j}\right)}\,\left[2\,n_{R}\left(\omega_{i}\right)+1\right]$
$\displaystyle\times\text{Im}\left[\mathcal{J}\left(\omega_{i},\omega_{j},t\right)e^{-i\left(\omega_{i}-\omega_{j}\right)t}-\mathcal{J}\left(\omega_{i},-\omega_{j},t\right)e^{-i\left(\omega_{i}+\omega_{j}\right)t}\right]$
$\displaystyle+\frac{1}{2m}\frac{1}{\sqrt{m_{i}\omega_{i}}}\frac{1}{\sqrt{m_{j}\omega_{j}}}\Delta\omega\sqrt{I_{R}\left(\omega_{i}\right)I_{L}\left(\omega_{j}\right)}\,\left[2\,n_{L}\left(\omega_{j}\right)+1\right]$
$\displaystyle\times\text{Im}\left[\mathcal{J}\left(\omega_{j},\omega_{i},t\right)e^{i\left(\omega_{i}-\omega_{j}\right)t}-\mathcal{J}\left(\omega_{j},-\omega_{i},t\right)e^{-i\left(\omega_{i}+\omega_{j}\right)t}\right]$
$\displaystyle+\frac{1}{\sqrt{m_{i}\omega_{i}}}\frac{1}{\sqrt{m_{j}\omega_{j}}}\,\Delta\omega\sqrt{I_{R}\left(\omega_{i}\right)I_{L}\left(\omega_{j}\right)}\int_{0}^{t}dt_{1}\int_{0}^{t}dt_{2}\,\text{sin}\left[\omega_{i}\left(t-t_{1}\right)\right]\,\text{sin}\left[\omega_{j}\left(t-t_{2}\right)\right]\left\langle\left\\{x^{h}\left(t_{1}\right),x^{h}\left(t_{2}\right)\right\\}\right\rangle$
$\displaystyle+\frac{1}{4m^{2}}\frac{1}{\sqrt{m_{i}\omega_{i}}}\frac{1}{\sqrt{m_{j}\omega_{j}}}\Delta\omega\sqrt{I_{R}\left(\omega_{i}\right)I_{L}\left(\omega_{j}\right)}\sum_{\alpha}\int_{0}^{\infty}d\omega\,I_{\alpha}\left(\omega\right)\,\left[2\,n_{\alpha}\left(\omega\right)+1\right]$
$\displaystyle\times\text{Re}\left[\mathcal{J}\left(\omega,\omega_{i},t\right)\mathcal{J}^{\ast}\left(\omega,\omega_{j},t\right)e^{i\left(\omega_{i}-\omega_{j}\right)t}+\mathcal{J}\left(\omega,-\omega_{i},t\right)\mathcal{J}^{\ast}\left(\omega,-\omega_{j},t\right)e^{-i\left(\omega_{i}-\omega_{j}\right)t}\right.$
$\displaystyle\left.-\mathcal{J}\left(\omega,\omega_{i},t\right)\mathcal{J}^{\ast}\left(\omega,-\omega_{j},t\right)e^{i\left(\omega_{i}+\omega_{j}\right)t}-\mathcal{J}\left(\omega,-\omega_{i},t\right)\mathcal{J}^{\ast}\left(\omega,\omega_{j},t\right)e^{-i\left(\omega_{i}+\omega_{j}\right)t}\right].$

and the momentum correlator is

$\displaystyle\left\langle\left\\{p_{i}\left(t\right),p_{j}\left(t\right)\right\\}\right\rangle$
$\displaystyle=m_{i}\,\omega_{i}\,\left[2\,n_{R}\left(\omega_{i}\right)+1\right]\,\delta_{ij}$
(25)
$\displaystyle+\frac{1}{2m}\sqrt{m_{i}\omega_{i}}\sqrt{m_{j}\omega_{j}}\Delta\omega\sqrt{I_{R}\left(\omega_{i}\right)I_{L}\left(\omega_{j}\right)}\,\left[2\,n_{R}\left(\omega_{i}\right)+1\right]$
$\displaystyle\times\text{Im}\left[\mathcal{J}\left(\omega_{i},\omega_{j},t\right)e^{-i\left(\omega_{i}-\omega_{j}\right)t}+\mathcal{J}\left(\omega_{i},-\omega_{j},t\right)e^{-i\left(\omega_{i}+\omega_{j}\right)t}\right]$
$\displaystyle+\frac{1}{2m}\sqrt{m_{i}\omega_{i}}\sqrt{m_{j}\omega_{j}}\Delta\omega\sqrt{I_{R}\left(\omega_{i}\right)I_{L}\left(\omega_{j}\right)}\,\left[2\,n_{L}\left(\omega_{j}\right)+1\right]$
$\displaystyle\times\text{Im}\left[\mathcal{J}\left(\omega_{j},\omega_{i},t\right)e^{i\left(\omega_{i}-\omega_{j}\right)t}+\mathcal{J}\left(\omega_{j},-\omega_{i},t\right)e^{-i\left(\omega_{i}+\omega_{j}\right)t}\right]$
$\displaystyle+\sqrt{m_{i}\omega_{i}}\sqrt{m_{j}\omega_{j}}\,\Delta\omega\sqrt{I_{R}\left(\omega_{i}\right)I_{L}\left(\omega_{j}\right)}\int_{0}^{t}dt_{1}\int_{0}^{t}dt_{2}\,\text{cos}\left[\omega_{i}\left(t-t_{1}\right)\right]\,\text{cos}\left[\omega_{j}\left(t-t_{2}\right)\right]\left\langle\left\\{x^{h}\left(t_{1}\right),x^{h}\left(t_{2}\right)\right\\}\right\rangle$
$\displaystyle+\frac{1}{4m^{2}}\sqrt{m_{i}\omega_{i}}\sqrt{m_{j}\omega_{j}}\Delta\omega\sqrt{I_{R}\left(\omega_{i}\right)I_{L}\left(\omega_{j}\right)}\sum_{\alpha}\int_{0}^{\infty}d\omega\,I_{\alpha}\left(\omega\right)\,\left[2\,n_{\alpha}\left(\omega\right)+1\right]$
$\displaystyle\times\text{Re}\left[\mathcal{J}\left(\omega,\omega_{i},t\right)\mathcal{J}^{\ast}\left(\omega,\omega_{j},t\right)e^{i\left(\omega_{i}-\omega_{j}\right)t}+\mathcal{J}\left(\omega,-\omega_{i},t\right)\mathcal{J}^{\ast}\left(\omega,-\omega_{j},t\right)e^{-i\left(\omega_{i}-\omega_{j}\right)t}\right.$
$\displaystyle\left.+\mathcal{J}\left(\omega,\omega_{i},t\right)\mathcal{J}^{\ast}\left(\omega,-\omega_{j},t\right)e^{i\left(\omega_{i}+\omega_{j}\right)t}+\mathcal{J}\left(\omega,-\omega_{i},t\right)\mathcal{J}^{\ast}\left(\omega,\omega_{j},t\right)e^{-i\left(\omega_{i}+\omega_{j}\right)t}\right],$

where the $\mathcal{J}$ function is defined as

$\mathcal{J}\left(\omega,\omega_{i},t\right)=\sum_{k}\int_{0}^{t}dt^{\prime}\,e^{i\left(\omega-\omega_{i}+k\omega_{d}\right)t^{\prime}}\int_{0}^{t^{\prime}}dt^{\prime\prime}\,A_{k}\left(t^{\prime\prime}\right)\,e^{-i\omega
t^{\prime\prime}},$ (26)

which, by integrating the exponentials with a change of the order of
integration, can be formally solved as

$\mathcal{J}\left(\omega,\omega_{i},t\right)=\sum_{k}\left[t\,\text{sinc}\left[\left(\omega-\omega_{i}+k\omega_{d}\right)t/2\right]\,a_{k}\left(i\omega\right)\,e^{i\left(\omega-\omega_{i}+k\omega_{d}\right)t/2}+F_{k}\left(\omega,\omega_{i}\right)\right],$
(27)

where we used the notation

$\displaystyle a_{k}\left(i\omega\right)$
$\displaystyle=\int_{0}^{t}dt^{\prime}\,A_{k}\left(t^{\prime}\right)\,e^{-i\omega
t^{\prime}}$ (28) $\displaystyle F_{k}\left(\omega,\omega_{i}\right)$
$\displaystyle=\frac{a_{k}\left(i\omega\right)-a_{k}\left[i\left(\omega_{i}-k\omega_{d}\right)\right]}{i\left(\omega-\omega_{i}+k\omega_{d}\right)}$

Both $a_{k}$ and $F_{k}$ are functions of the variable $t$, but we do not
write its explicit dependence in an effort to keep the notation simple. When
computing correlators and quantities related to them, it is sometimes useful
to use a relation which is a direct consequence of the unitary evolution of
the operators $q_{i}$ and $p_{i}$:

$\displaystyle\text{Im}\left[\mathcal{J}\left(\omega_{i},\omega_{i},t\right)\right]$
$\displaystyle=-\frac{1}{4m}\sum_{\alpha}\int_{0}^{\infty}d\omega\,I_{\alpha}\left(\omega\right)\left\lvert\mathcal{J}\left(\omega,\omega_{i},t\right)\right\rvert^{2}+\frac{1}{4m}\sum_{\alpha}\int_{0}^{\infty}d\omega\,I_{\alpha}\left(\omega\right)\left\lvert\mathcal{J}\left(\omega,-\omega_{i},t\right)\right\rvert^{2}$
(29)
$\displaystyle-\frac{1}{2}m\int_{0}^{t}dt_{1}\int_{0}^{t}dt_{2}\,\text{sin}\left[\omega_{i}\left(t_{1}-t_{2}\right)\right]\left\langle
i\left[x^{h}\left(t_{1}\right),x^{h}\left(t_{2}\right)\right]\right\rangle.$

The previous equality is obtained by imposing that
$[q_{i}(t),p_{i}(t)]=i\mathbb{1}$ for all times to the solutions of the
equations of motion shown in the main text. We can identify the different
parts of the covariance matrix as written in Eq. (5) by looking at, for
example, the position correlator in Eq. (24). The initial thermal state is

$\boldsymbol{\nu}_{i}\otimes\boldsymbol{\nu}_{j}\rightarrow\frac{1}{2}\,\left[2\,n_{R}\left(\omega_{i}\right)+1\right]\,\delta_{ij}.$
(30)

The term which involves the initial state of $\mathcal{S}$ is

$\boldsymbol{\sigma}_{0}(t)\rightarrow\frac{1}{2}\,\Delta\omega\sqrt{I_{R}\left(\omega_{i}\right)I_{L}\left(\omega_{j}\right)}\int_{0}^{t}dt_{1}\int_{0}^{t}dt_{2}\,\text{sin}\left[\omega_{i}\left(t-t_{1}\right)\right]\,\text{sin}\left[\omega_{j}\left(t-t_{2}\right)\right]\left\langle\left\\{x^{h}\left(t_{1}\right),x^{h}\left(t_{2}\right)\right\\}\right\rangle.$
(31)

The other two terms, $\boldsymbol{\sigma}_{osc}(t)$ and
$\boldsymbol{\sigma}_{lin}\times t$, are given by the rest of Eq. (24).

Using the correlators in Eq. (24) and Eq. (25) we can obtain $E_{i}$, which is
the energy stored in the band $i\in\mathcal{E}_{R}$ presented in Eq. (7), as

$E_{i}=\frac{1}{4m_{i}}\left\langle\left\\{p_{i}\left(t\right),p_{i}\left(t\right)\right\\}\right\rangle+\frac{1}{4}m_{i}\omega_{i}^{2}\left\langle\left\\{q_{i}\left(t\right),q_{i}\left(t\right)\right\\}\right\rangle.$
(32)

A direct computation using relation (29) shows that

$\displaystyle E_{i}\left(t\right)$
$\displaystyle=\frac{1}{2}\omega_{i}\,\left[2\,n_{R}\left(\omega_{i}\right)+1\right]$
(33)
$\displaystyle+\frac{1}{2}\Delta\omega\,\omega_{i}\,I_{R}\left(\omega_{i}\right)\left[2\,n_{R}\left(\omega_{i}\right)+1\right]\int_{0}^{t}dt_{1}\int_{0}^{t}dt_{2}\,\left\langle\hat{x}^{h}\left(t_{1}\right)\hat{x}^{h}\left(t_{2}\right)\right\rangle\,e^{-i\omega_{i}\left(t_{1}-t_{2}\right)}$
$\displaystyle-\frac{1}{2}\Delta\omega\,\omega_{i}\,I_{R}\left(\omega_{i}\right)n_{R}\left(\omega_{i}\right)\int_{0}^{t}dt_{1}\int_{0}^{t}dt_{2}\,\text{cos}\left[\omega_{i}\left(t_{1}-t_{2}\right)\right]\left\langle\left\\{x^{h}\left(t_{1}\right),x^{h}\left(t_{2}\right)\right\\}\right\rangle$
$\displaystyle+\frac{1}{4m^{2}}\Delta\omega\,\omega_{i}\,I_{R}\left(\omega_{i}\right)\sum_{\alpha=R,L}\int_{0}^{\infty}d\omega\,I_{\alpha}\left(\omega\right)\left[n_{\alpha}\left(\omega\right)-n_{R}\left(\omega_{i}\right)\right]\left\lvert\mathcal{J}\left(\omega,\omega_{i},t\right)\right\rvert^{2}$
$\displaystyle+\frac{1}{4m^{2}}\Delta\omega\,\omega_{i}\,I_{R}\left(\omega_{i}\right)\sum_{\alpha=R,L}\int_{0}^{\infty}d\omega\,I_{\alpha}\left(\omega\right)\left[n_{\alpha}\left(\omega\right)+n_{R}\left(\omega_{i}\right)+1\right]\left\lvert\mathcal{J}\left(\omega,-\omega_{i},t\right)\right\rvert^{2}.$

In order to arrive at the final expression, we need the long-time behavior of
the $\mathcal{J}$ functions. For this, we will use the fact that

$\lim_{t\to\infty}t\,\text{sinc}\left(\omega\,t/2\right)=\lim_{t\to\infty}t\,\text{sinc}^{2}\left(\omega\,t/2\right)=2\pi\,\delta\left(\omega\right),$
(34)

and $\lim_{t\to\infty}a_{k}(i\omega)=\tilde{A}_{k}(i\omega)$. Thus,

$\left\lvert\mathcal{J}\left(\omega,\pm\omega_{i},t\right)\right\rvert^{2}\to
2\pi
t\,\sum_{k}\delta(\omega\mp\omega_{i}+k\omega_{d})\lvert\tilde{A}_{k}(i\omega)\rvert^{2},$
(35)

where we only kept the linear terms in time. Using this, and that
$\tilde{A}_{k}^{\ast}(i\omega)=\tilde{A}_{-k}(-i\omega)$ (which is a
consequence of $G$ being a real valued function), we arrive at

$E_{i}(t)\to\left[1/2+n_{R}\left(\omega_{i}\right)\right]\omega_{i}+\dot{\mathcal{Q}}_{i}\times
t$ (36)

with $\dot{\mathcal{Q}}_{i}$ as shown in Eq. (7) (remember that the terms
involving the initial state of $\mathcal{S}$ decay exponentially fast for
long-times).

##  Appendix B Individual purities and $\Gamma_{\pm}$

Here we present explicit expressions for the individual purities $\mu_{i,j}$
and for $\Gamma_{\pm}$, in the long-time regime (that is, given the covariance
matrix $\boldsymbol{\sigma}_{av}$ in Eq. (6)). As discussed above, $\mu_{i,j}$
are related to $\boldsymbol{\sigma}_{av}$ by
$\mu_{i}=1/2\sqrt{\text{det}(\boldsymbol{\alpha})}$ and
$\mu_{j}=1/2\sqrt{\text{det}(\boldsymbol{\beta})}$. Since the purities do not
depend on the relation between the frequencies $\omega_{i}$ and $\omega_{j}$,
we just show

$\displaystyle 2\sqrt{\text{det}(\boldsymbol{\alpha})}$
$\displaystyle=1+2n_{R}\left(\omega_{i}\right)$ (37)
$\displaystyle+\frac{\pi}{m^{2}}\,\Delta\omega\,I_{R}\left(\omega_{i}\right)\sum_{\alpha=R,L}\sum_{q^{\prime}}I_{\alpha}\left(\omega_{i}-q^{\prime}\omega_{d}\right)\left[n_{\alpha}\left(\omega_{i}-q^{\prime}\omega_{d}\right)-n_{R}\left(\omega_{i}\right)\right]\left\lvert\tilde{A}_{q^{\prime}}\left[i\left(\omega_{i}-q^{\prime}\omega_{d}\right)\right]\right\rvert^{2}\times
t$
$\displaystyle+\frac{\pi}{m^{2}}\,\Delta\omega\,I_{R}\left(\omega_{i}\right)\sum_{\alpha=R,L}\sum_{q^{\prime\prime}}I_{\alpha}\left(q^{\prime\prime}\omega_{d}-\omega_{i}\right)\left[n_{\alpha}\left(q^{\prime\prime}\omega_{d}-\omega_{i}\right)+n_{R}\left(\omega_{i}\right)+1\right]\left\lvert\tilde{A}_{-q^{\prime\prime}}\left[i\left(q^{\prime\prime}\omega_{d}-\omega_{i}\right)\right]\right\rvert^{2}\times
t$

where $q^{\prime}$ and $q^{\prime\prime}$ are integers such that
$\omega_{i}-q^{\prime}\omega_{d}\geq 0$ and
$q^{\prime\prime}\omega_{d}-\omega_{i}\geq 0$, respectively. In order to
obtain $2\sqrt{\text{det}(\boldsymbol{\beta})}$ we just need to make the
replacement $\left\\{i,R\right\\}\rightarrow\left\\{j,L\right\\}$. On the
contrary, $\Gamma_{\pm}$ depends on whether
$\omega_{i}+\omega_{j}=k\omega_{d}$ or $\omega_{j}=\omega_{i}+k\omega_{d}$.
Thus, for the nonresonant case ($\omega_{i}+\omega_{j}=k\omega_{d}$) we have

$\Gamma_{+}=-\frac{1}{4m^{2}}\Delta\omega^{2}I_{R}\left(\omega_{i}\right)I_{L}\left(\omega_{j}\right)\left\lvert\left[2n_{R}\left(\omega_{i}\right)+1\right]\tilde{A}_{-k}^{\ast}\left(i\omega_{i}\right)+\left[2n_{L}\left(\omega_{j}\right)+1\right]\tilde{A}_{-k}^{\ast}\left(i\omega_{j}\right)\right\rvert^{2},$
(38)

and for the resonant one ($\omega_{j}=\omega_{i}+k\omega_{d}$),

$\Gamma_{-}=\frac{1}{m^{2}}\Delta\omega^{2}I_{R}\left(\omega_{i}\right)I_{L}\left(\omega_{j}\right)\left\lvert
n_{R}\left(\omega_{i}\right)\tilde{A}_{k}\left(i\omega_{i}\right)-n_{L}\left(\omega_{j}\right)\tilde{A}_{-k}^{\ast}\left(i\omega_{j}\right)\right\rvert^{2}.$
(39)

It is clear from Eqs. (37) and (39) that the inequality in Eq. (13) holds. To
see this, we can write the inequality as

$\lvert\Gamma_{-}(t)\rvert\leq(1/\mu_{i}-1)(1/\mu_{j}-1).$ (40)

From Eq. (39), we know that when temperature goes to zero, the left hand side
of the inequality goes to zero too. On the contrary, from the third line of
Eq. (37), the right hand side does not go to zero. In the general case, the
left hand side of the inequality is of order
$(\Delta\omega^{2}\gamma_{0}^{2})n(\omega)$, while the right hand side is of
order $n(\omega)^{2}+(\Delta\omega\gamma_{0}^{2})n(\omega)$ which is much
bigger in the weak coupling limit. So, as temperature increases, the right
hand side is always greater than the left hand side.

##  Appendix C Symplectic eigenvalues and measures of correlations

Here we will briefly show how to obtain the expressions of the mutual
information, quantum discord and logarithmic negativity presented above. In
order to do this, first we will compute the symplectic eigenvalues of
$\boldsymbol{\sigma}_{av}$. For the sake of clarity we will use the standard
notation: $A=\text{det}(\boldsymbol{\alpha})=1/4\mu_{i}^{2}$,
$B=\text{det}(\boldsymbol{\beta})=1/4\mu_{j}^{2}$,
$C=\text{det}(\boldsymbol{\gamma})=\Gamma(t)/4$ and
$D=\text{det}(\boldsymbol{\sigma}_{av})$. We remind the reader that all
computations are valid in the weak coupling regime and for long times.
Therefore, in many steps we will use Taylor expansions to first nontrivial
order in $\lvert C\rvert$ when comparing $\lvert C\rvert\sim
O(\gamma_{0}^{2})$ to $\lvert A\rvert,\lvert B\rvert\sim O(1)$ (see Eqs. (37),
(38) and (39)). The symplectic eigenvalues are

$\lambda_{1,2}^{2}=\frac{1}{2}(\Delta\pm\sqrt{\Delta^{2}-4\,D}),$ (41)

where $\Delta=A+B+2C$ [38, 39]. In the nonresonant case we have
$\Delta=A+B-2\lvert C\rvert$ and $D=(\sqrt{AB}-\lvert C\rvert)^{2}$ (see Eq.
(8)). Thus, after a Taylor expansion we obtain

$\sqrt{\Delta^{2}-4\,D}\simeq(A-B)-2\frac{\sqrt{A}-\sqrt{B}}{\sqrt{A}+\sqrt{B}}\lvert
C\rvert.$ (42)

Therefore, after another expansion, the result is

$\displaystyle\lambda_{1}\simeq\sqrt{A}-\frac{\lvert
C\rvert}{\sqrt{A}+\sqrt{B}}$ (43)
$\displaystyle\lambda_{2}\simeq\sqrt{B}-\frac{\lvert
C\rvert}{\sqrt{A}+\sqrt{B}}$

On the other hand, in the resonant case we have $\Delta=A+B+2C$ and
$D=(\sqrt{AB}-C)^{2}$. Then, after a Taylor expansion we get

$\sqrt{\Delta^{2}-4\,D}\simeq(A-B)+2\frac{\sqrt{A}+\sqrt{B}}{\sqrt{A}-\sqrt{B}}C.$
(44)

In order to arrive at the previous expression, it was used that
$(\sqrt{A}-\sqrt{B})^{2}>4C$. This in an immediate consequence of the
inequality in Eq. (13). The symplectic eigenvalues in this case are

$\displaystyle\lambda_{1}\simeq\sqrt{A}+\frac{C}{\sqrt{A}-\sqrt{B}}$ (45)
$\displaystyle\lambda_{2}\simeq\sqrt{B}-\frac{C}{\sqrt{A}-\sqrt{B}}$

Now we can proceed to compute the mutual information and quantum discord. The
mutual information is defined as [39]

$\mathcal{I}(\boldsymbol{\sigma})=f(\sqrt{A})+f(\sqrt{B})-f(\lambda_{1})-f(\lambda_{2}),$
(46)

where $f(x)=(x+1/2)\,\text{ln}(x+1/2)-(x-1/2)\,\text{ln}(x-1/2)$. Using a
first order Taylor expansion in $\lvert C\rvert/(\sqrt{A}\pm\sqrt{B})$, we get

$\mathcal{I}_{\pm}(\boldsymbol{\sigma}_{av})\simeq\pm[f^{\prime}(\sqrt{A})\pm
f^{\prime}(\sqrt{B})]\,\frac{\lvert C\rvert}{\sqrt{A}\pm\sqrt{B}}$ (47)

This is the expression shown in (11). The quantum discord is defined as [36]:

$\overleftarrow{D}(\boldsymbol{\sigma})=f(\sqrt{B})-f(\lambda_{1})-f(\lambda_{2})+f(\sqrt{E_{min}}).$
(48)

In our case, $E_{\min}$ is

$E_{\min}=\frac{2C^{2}+(1/4-B)(A-4D)+2\lvert
C\rvert\sqrt{C^{2}+(1/4-B)(A-4D)}}{4(1/4-B)^{2}}.$ (49)

Replacing $D=(\sqrt{AB}-\lvert C\rvert)^{2}$ in the equation above, the
expression for $E_{\min}$ can be greatly simplified. Indeed, noticing that

$\frac{2\lvert
C\rvert^{2}+(1/4-B)(A-4D)}{4(1/4-B)^{2}}=\left(\sqrt{A}-\frac{\sqrt{B}\lvert
C\rvert}{\lvert 1/4-B\rvert}\right)^{2}+\left(\frac{\lvert C\rvert}{2\lvert
1/4-B\rvert}\right)^{2}$ (50)

and

$\frac{2\lvert C\rvert\sqrt{C^{2}+(1/4-B)(A-4D)}}{4(1/4-B)^{2}}=\frac{\lvert
C\rvert}{\lvert 1/4-B\rvert}\left(\sqrt{A}-\frac{\sqrt{B}\lvert
C\rvert}{\lvert 1/4-B\rvert}\right),$ (51)

we arrive at

$E_{min}=\left(\sqrt{A}-\frac{\lvert C\rvert}{\lvert
1/2+\sqrt{B}\rvert}\right)^{2},$ (52)

where we used the fact that $1/4-B<0$ and $\sqrt{A}>\sqrt{B}\lvert
C\rvert/\lvert 1/4-B\rvert$. Thus, using a first order Taylor expansion in
$\lvert C\rvert/(\lvert 1/2+\sqrt{B}\rvert)$, we get

$f(\sqrt{E_{min}})\simeq f(\sqrt{A})-f^{\prime}(\sqrt{A})\frac{\lvert
C\rvert}{\lvert 1/2+\sqrt{B}\rvert}.$ (53)

Finally, the quantum discord can be written as shown in Eq. (15):

$\overleftarrow{D}_{\pm}(\boldsymbol{\sigma}_{av})\simeq\mathcal{I}_{\pm}(\boldsymbol{\sigma}_{av})-f^{\prime}(\sqrt{A})\frac{\lvert
C\rvert}{\lvert 1/2+\sqrt{B}\rvert}.$ (54)

For the logarithmic negativity, we need the lowest symplectic eigenvalue of
the transposed covariance matrix $\tilde{\boldsymbol{\sigma}}_{av}$. This is,

$\tilde{\lambda}_{2}^{2}=\frac{1}{2}(\tilde{\Delta}-\sqrt{\tilde{\Delta}^{2}-4\,D}),$
(55)

where $\tilde{\Delta}=A+B-2C$ [38, 39]. In this case, we have

$\sqrt{\tilde{\Delta}^{2}-4\,D}\simeq 2\sqrt{\lvert
C\rvert}(\sqrt{A}+\sqrt{B})\left[1+\frac{(\sqrt{A}-\sqrt{B})^{2}}{8\lvert
C\rvert}\right]$ (56)

and therefore,

$\tilde{\lambda}_{2}^{2}\simeq\frac{1}{2}(A+B)\left[1-\frac{\sqrt{A}+\sqrt{B}}{A+B}\sqrt{\lvert
C\rvert}\right].$ (57)

The logarithmic negativity as shown in Eq. (21) is obtained after computing
$\text{max}\\{0,-\text{ln}(2\tilde{\lambda}_{-})\\}$ and using the fact that
$\text{ln}(1+x)\simeq x$ for small $x$.

## References

  * [1] Nielsen M. and Chuang I. Quantum Computation and Quantum Information: 10th Anniversary Edition. Cambridge University Press, 2010
  * [2] Vinjanampathy S. and Anders J. Quantum thermodynamics. Contemporary Physics, 57(4), 2016.
  * [3] Esposito M., Lindenberg K., Van den Broeck C. Entropy production as correlation between system and reservoir. New Journal of Physics, 12(1):013013, 2010
  * [4] Ptaszyński K. and Esposito M. Entropy Production in Open Systems: The Predominant Role of Intraenvironment Correlations. Physical Review Letters, 123(20):200603, 2019
  * [5] Bera M.N., Riera A., Lewenstein M. and Winter, A. Generalized laws of thermodynamics in the presence of correlations. Nature Communications, 8, 2180, 2017.
  * [6] Jennings D. and Rudolph T. Entanglement and the thermodynamic arrow of time. Physical Review E, 81(6):061130, 2010.
  * [7] del Rio L., Åberg J., Renner R., Dahlsten O., and Vedral V. The thermodynamic meaning of negative entropy. Nature, 474, 61-63, 2011.
  * [8] Åberg J. Truly work-like work extraction via a single-shot analysis. Nature Communications, 4, 1925, 2013.
  * [9] Sapienza F., Cerisola F., and Roncaglia A. Correlations as a resource in quantum thermodynamics. Nature Communications, 10, 2492, 2019.
  * [10] Müller M. Correlating Thermal Machines and the Second Law at the Nanoscale. Physical Review X, 8(4):041051, 2018.
  * [11] Perarnau-Llobet M., Hovhannisyan K., Huber M., Skrzypczyk P., Brunner N., and Acín A. Extractable Work from Correlations. Physical Review X, 5(4):041011, 2015.
  * [12] Francica G., Goold J., Plastina F., and Paternostro M. Daemonic ergotropy: enhanced work extraction from quantum correlations. npj Quantum Information, 3, 12, 2017.
  * [13] Aguilar M., Freitas N., and Paz J. P. Entanglement generation in quantum thermal machines. Physical Review A, 102(6):062422, 2020.
  * [14] Brunner N., Huber M., Linden N., Popescu S., Silva R. and Skrzypczyk P. Entanglement enhances cooling in microscopic quantum refrigerators. Physical Review E, 89(3):032115, 2014.
  * [15] Bohr Brask J., Haack G., Brunner N. and Huber M. Autonomous quantum thermal machine for generating steady-state entanglement. New Journal of Physics, 17(11):113029, 2015.
  * [16] Khandelwal S., Palazzo N., Brunner N. and Haack G. Critical heat current for operating an entanglement engine. New Journal of Physics, 22:073039, 2020.
  * [17] Hu B. L., Paz J. P., and Zhang Y. Quantum Brownian motion in a general environment: Exact master equation with nonlocal dissipation and colored noise. Physical Review D, 25(8):2843, 1992.
  * [18] Caldeira A. and Leggett A. Path integral approach to quantum Brownian motion. Physica A, 121(3), 1983.
  * [19] Landau L. and Lifshitz E. Course of Theoretical Physics. Mechanics. Vol. 1 Elsevier, 1982
  * [20] Freitas N. and Paz J. P. Fundamental limits for cooling of linear quantum refrigerators. Physical Review E, 95(1):012146, 2017.
  * [21] Freitas N. and Paz J. P. Cooling a quantum oscillator: A useful analogy to understand laser cooling as a thermodynamical process. Physical Review A, 97(3):032104, 2018.
  * [22] Lu-Ming D., Giedke G., Cirac J. I., and Zoller P. textitInseparability Criterion for Continuous Variable Systems. Physical Review Letters, 84(12):2722, 2000.
  * [23] Araki H. and Lieb, E Entropy inequalities. Communications in Mathematical Physics, 18, 160-170, 1970.
  * [24] Ollivier H., and Zurek W. Quantum Discord: A Measure of the Quantumness of Correlations. Physical Review Letters, 88(1):017901, 2001.
  * [25] Henderson L., and Vedral V. Classical, quantum and total correlations. Journal of Physics A, 34:6899, 2001.
  * [26] Knill E., and Laflamme R. Power of One Bit of Quantum Information. Physical Review Letters, 81(25):5672, 1998.
  * [27] Bennett C., DiVincenzo D., Fuchs C., Mor T., Rains E., Shor P., Smolin J., and Wootters W. Quantum nonlocality without entanglement. Physical Review A., 59(2):1070, 1999.
  * [28] Horodecki M., Sen(De) A., Sen U., and Horodecki K. Local Indistinguishability: More Nonlocality with Less Entanglement. Physical Review Letters, 90(4):047902, 2003.
  * [29] Ferraro A., Aolita L., Cavalcanti D., Cucchietti F., and Acín A. Almost all quantum states have nonclassical correlations. Physical Review A., 81(5):052378, 2010.
  * [30] Datta A., Shaji A., and Caves C. Quantum Discord and the Power of One Qubit. Physical Review Letters, 100(5):050502, 2008.
  * [31] Rodríguez-Rosario C., Modi K., Kuah A., Shaji A., and Sudarshan E. Completely positive maps and classical correlations. Journal of Physics A, 41:205301, 2008.
  * [32] Piani M., Horodecki P., and Horodecki R. No-Local-Broadcasting Theorem for Multipartite Quantum Correlations. Physical Review Letters, 100(9):090502, 2008.
  * [33] Lanyon B., Barbieri M., Almeida M., and White A. Experimental Quantum Computing without Entanglement. Physical Review Letters, 101(9):200501, 2008.
  * [34] Freitas N. and Paz J. P. Dynamics of Gaussian discord between two oscillators interacting with a common environment. Physical Review A, 85(3):032118, 2012.
  * [35] Giorda P., and Paris M. Gaussian Quantum Discord. Physical Review Letters, 105(2):020503, 2010.
  * [36] Adesso G., and Datta A. Quantum versus classical correlations in Gaussian states. Physical Review Letters, 105(3):030501, 2010.
  * [37] Pirandola S., Spedalieri G., Braunstein S., Cerf N., and Lloyd S. Optimality of Gaussian Discord. Physical Review Letters, 113(14):140405, 2014.
  * [38] Adesso G. and Illuminati F. Gaussian measures of entanglement versus negativities: Ordering of two-mode Gaussian states. Physical Review A, 72(3):032334, 2005.
  * [39] Serafini A., Illuminati F., and De Siena S. Symplectic invariants, entropic measures and correlations of Gaussian states. Journal of Physics B, 37(2), 2003.
  * [40] Adesso G., Girolami D., and Serafini A. Measuring Gaussian Quantum Information and Correlations Using the Rényi Entropy of Order 2. Physical Review Letters, 109(19):190502, 2012.