#  Driven quantum harmonic oscillators: A working medium for thermal machines

Heather Leitch  Centre for Theoretical Atomic, Molecular and Optical Physics,
Queen’s University Belfast, Belfast BT7 1NN, United Kingdom  Nicolò Piccione
Bruno Bellomo  Institut UTINAM, CNRS UMR 6213, Université Bourgogne Franche-
Comté, Observatoire des Sciences de l’Univers THETA, 41 bis avenue de
l’Observatoire, F-25010 Besançon, France  Gabriele De Chiara  Centre for
Theoretical Atomic, Molecular and Optical Physics, Queen’s University Belfast,
Belfast BT7 1NN, United Kingdom

###### Abstract

The study of quantum thermodynamics is key to the development of quantum
thermal machines. In contrast to most of the previous proposals based on
discrete strokes, here we consider a working substance that is permanently
coupled to two or more baths at different temperatures and continuously
driven. To this end, we investigate parametrically driven quantum harmonic
oscillators coupled to heat baths via a collision model. Using a
thermodynamically consistent local master equation, we derive the heat flows
and power of the working device which can operate as an engine, refrigerator
or accelerator and analyze the instantaneous and average efficiencies and
coefficients of performance. Studying the regimes of both slow and fast
driving of the system, we find that an increased driving frequency can lead to
a change of functioning to a dissipator. Finally, we investigate the effect of
squeezing one of the thermal baths: it leads to an apparent higher efficiency
compared to the corresponding Carnot value of an equilibrium bath with the
same temperature and to sustained entanglement between the working substance
oscillators in the limit cycle.

##  I Introduction

The study of quantum thermal machines and their potential advantages over
their classical counterparts has gained much interest in recent years Kosloff
and Levy (2014); Deffner and Campbell (2019); Mitchison (2019), after the
first experimental realizations had been demonstrated Roßnagel  _et al._
(2016); Maslennikov  _et al._ (2019); Peterson  _et al._ (2019); Gluza  _et
al._ (2021). The theory of open quantum systems provides us with the tools to
calculate the dynamics of the working substance interacting with one or more
much larger environments. Studying the energy flows which result from a system
interacting with heat baths of different temperatures is key to understanding
the conditions required for thermal machines such as engines and refrigerators
to operate.

The working substance itself may consist of quantum systems of different
nature, including quantum harmonic oscillators (QHO’s) Kosloff and Rezek
(2017); Singh and Müstecaplıoğlu (2020); Reid  _et al._ (2017); Deffner
(2018); Serafini  _et al._ (2020), qubits Quan  _et al._ (2007); Hewgill  _et
al._ (2018); Solfanelli  _et al._ (2020) or multilevel systems Quan  _et al._
(2005); Campisi and Fazio (2016); Niedenzu and Kurizki (2018); Piccione  _et
al._ (2021). Thermal machines can be normally split into three categories
depending on how the working substance is operated upon: autonomous (where the
Hamiltonian of the working substance is time-independent), discrete-stroke
(e.g. the Otto or Carnot cycle, in which the evolution is split into finite-
time intervals in which the baths are connected and disconnected) and
continuously driven (in which the system is permanently coupled to the baths
and its Hamiltonian is time-dependent). Though the equivalence of these
categories has been proven Uzdin  _et al._ (2015), for practical applications
it is still significant to maintain the separation. Besides, a no-go theorem
exists for autonomous linearly-coupled QHO to operate as an absorption
refrigerator Martinez and Paz (2013).

In this paper, we consider thermal machines of the third category: a network
of driven QHO’s permanently coupled to heat baths. The evolution of the
system, in the limit of weak coupling, is described by an adiabatic Markovian
Lindblad master equation within the local approximation. The thermodynamic
consistency of local master equations is guaranteed by their microscopic model
based on the so-called repeated interactions or collision models, Refs. De
Chiara  _et al._ (2018) and Hewgill  _et al._ (2021) (see Refs.Campbell and
Vacchini (2021) and Ciccarello  _et al._ (2021) for recent reviews).

We focus on working substances made of one and two driven QHO’s interacting
with two baths and find the conditions required for the system to operate as
an engine, a refrigerator or an accelerator. The effect of both slow and fast
driving is explored and we see that the increased driving speed leads to a
larger injection of power into the system causing heat to flow from the system
into both baths. This operation is referred to as a dissipator or heater
Solfanelli  _et al._ (2020). We study the efficiency and the coefficient of
performance (COP) of the device when operating as an engine and as a
refrigerator, respectively. We show that driving the system oscillators causes
the power to oscillate periodically in time. In some cases, this leads to the
instantaneous coefficient of performance of a refrigerator being higher than
the Otto value. We also show how the presence of squeezing in one of the baths
effectively raises the temperature of the bath resulting in an effective
higher Carnot efficiency. Squeezing is also required for the generation of
entanglement in a working medium made of two QHO’s.

This paper is organized as follows. In Sec. II we explain our microscopic
model and derive the relevant thermodynamic quantities. In Sec. III we
illustrate the results for the simple case of a single QHO coupled to two
baths at different temperatures while in Sec. IV we extend our study to a
working substance of two QHO’s. In Sec. V we introduce squeezing in the hotter
bath in the two-oscillator model showing an (apparent) enhanced efficiency and
the sustenance of entanglement between the two oscillators in the long-time
limit. Finally, in Sec. VI we summarize our findings and conclude. Some
details of our analysis can be found in the appendices.

##  II The Model

In this section we describe the general setup we are going to consider in this
paper. The system is composed of $N$ driven quantum harmonic oscillators
(QHO’s) coupled to $N_{B}$ baths via a collision model, as depicted in Fig. 1.

Figure 1: A network of $N$ coupled quantum harmonic oscillators interacting
with $N_{B}$ heat baths via the collision model.

Let us define the position $x_{i}$ and momentum $p_{i}$ operators of the ith
oscillator in the system subject to the time-dependent frequency
$\omega_{i}(t)$ and with unit mass. The time-dependent Hamiltonian for the
system is (setting $\hbar=1$),

$\displaystyle H_{S}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{N}\left[p^{2}_{i}+\omega^{2}_{i}(t)x_{i}^{2}\right]+$
(1) $\displaystyle+$
$\displaystyle\sum_{\begin{subarray}{c}{i<j}\end{subarray}}\lambda_{ij}\left[a_{i}^{\dagger}(t)a_{j}(t)+a_{i}(t)a_{j}^{\dagger}(t)\right],$

where the time-dependent annihilation operators $a_{i}(t)$ are defined as

$a_{i}(t)=\sqrt{\frac{\omega_{i}(t)}{2}}x_{i}+\mathbbm{i}\sqrt{\frac{1}{2\omega_{i}(t)}}p_{i}.$
(2)

In Eq. (1), $\lambda_{ij}$ (assumed to be real) is the strength of the beam-
splitter-type coupling between the ith and jth system oscillator and is equal
to zero if there is no coupling between them. The number-conserving form of
the coupling term ensures that excitations can hop between coupled oscillators
but are not destroyed or created in the process. The presence of counter-
rotating terms in the coupling, arising for instance from a position-position
coupling, tends to hinder the transport efficiency as discussed in Ref. De
Chiara  _et al._ (2018). The beam-splitter-type coupling of Eq. (1) may be
implemented in photonic waveguides as discussed in Ref. Reid  _et al._ (2017).

The oscillators’ frequencies are modulated in time according to

$\omega_{i}(t)=\omega_{0,i}+\delta\omega_{i}\sin(\theta_{i}t),$ (3)

with $\delta\omega_{i}<\omega_{0,i}$ to ensure positivity of the frequencies
and where $\theta_{i}$ is the modulation frequency.

The baths are modeled as an infinite ensemble of QHO’s such that the total
bath Hamiltonian is $H_{B}=\sum_{\alpha=1}^{N_{B}}H_{B,\alpha}$ where each
individual bath Hamiltonian reads

$H_{B,\alpha}=\frac{1}{2}\sum_{j=1}^{\infty}\left(P^{2}_{\alpha,j}+\Omega_{\alpha}^{2}X^{2}_{\alpha,j}\right),$
(4)

where we have defined the corresponding position $X_{\alpha,j}$ and momentum
$P_{\alpha,j}$ operators of the jth oscillator in the $\alpha$th bath. We also
define the corresponding time-independent annihilation operators

$b_{\alpha,j}=\sqrt{\frac{\Omega_{\alpha}}{2}}X_{\alpha,j}+\mathbbm{i}\sqrt{\frac{1}{2\Omega_{\alpha}}}P_{\alpha,j}.$
(5)

Note that we are considering the case where all the oscillators in each bath
have the same frequency. We assume that the oscillators in the baths are
initially prepared in a thermal state at their respective temperature
$1/\beta_{\alpha}$ (setting $k_{B}=1$),

$\rho_{B,\alpha}=\frac{e^{-\beta_{\alpha}H_{B,\alpha}}}{Z_{\alpha}},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;Z_{\alpha}=\text{Tr}\left[e^{-\beta_{\alpha}H_{B,\alpha}}\right],$
(6)

and interact with the system by means of repeated interactions or collisions.
At each time step, one oscillator from each bath interacts with one QHO in the
system for a short time period $\tau$. Each bath interacts always with the
same QHO. This interaction is described by the Hamiltonian

$H_{I}(t)=\sum_{i=1}^{N}\sum_{\alpha=1}^{N_{B}}\sum_{j=1}^{\infty}\frac{g_{i,\alpha}\Theta_{j}(t)}{\sqrt{\tau}}\left[a_{i}^{\dagger}(t)b_{\alpha,j}+a_{i}(t)b_{\alpha,j}^{\dagger}\right],$
(7)

where $g_{i,\alpha}$ is the coupling between the ith system oscillator and the
bath $\alpha$ and where

$\Theta_{j}(t)=\begin{cases}1&\text{for}\;t_{j}\leq t\leq t_{j}+\tau,\\\
0&\text{otherwise},\end{cases}$ (8)

with $t_{j}$ being the time at which the jth bath oscillator starts
interacting with the system. We assume $t_{1}=0$ and that there is no time
between subsequent collisions. In summary, the Hamiltonian of the complete
system and baths is

$H(t)=H_{S}(t)+H_{B}+H_{I}(t).$ (9)

A technical note: the $1/\sqrt{\tau}$ scaling in Eq. 7 ensures the consistency
in the continuous limit ($\tau\to 0$) we are now going to describe and
naturally arises in quantum optical systems Ciccarello  _et al._ (2021).

We now introduce the adiabatic master equation for the evolution of the
system’s density matrix in which non-adiabatic terms arising in the
dissipative part of the evolution are absent. To this end we assume that the
collision time $\tau$ is the shortest timescale in the problem and that the
modulation frequency $\theta_{i}$ corresponds to a much slower timescale
Albash  _et al._ (2012). This assumption results in a master equation that
does not contain derivatives of the system modulated frequencies explicitly.
The derivation then follows closely the one presented in Ref. De Chiara  _et
al._ (2018), where an expansion up to order $\tau^{2}$ for the evolved density
matrix is considered. Making the limit $\tau\to 0$ leads to the adiabatic
Markovian master equation,

$\displaystyle\dot{\rho}_{S}$ $\displaystyle=$
$\displaystyle-\mathbbm{i}\commutator{H_{S}(t)}{\rho_{S}}$ $\displaystyle+$
$\displaystyle\sum_{i=1}^{N}\sum_{\alpha=1}^{N_{B}}g_{i,\alpha}^{2}n_{\alpha}\left[a_{i}^{\dagger}\rho_{S}a_{i}-\frac{1}{2}\left\\{a_{i}a_{i}^{\dagger},\rho_{S}\right\\}\right]$
$\displaystyle+$
$\displaystyle\sum_{i=1}^{N}\sum_{\alpha=1}^{N_{B}}g_{i,\alpha}^{2}(n_{\alpha}+1)\left[a_{i}\rho_{S}a_{i}^{\dagger}-\frac{1}{2}\left\\{a_{i}^{\dagger}a_{i},\rho_{S}\right\\}\right],$

where for brevity we have omitted the time-dependence of $\rho_{S}(t)$ and
$a_{i}(t)$. We have also defined the thermal occupation of each bath
oscillator as

$n_{\alpha}=\frac{1}{e^{\beta_{\alpha}\Omega_{\alpha}}-1}.$ (11)

Note that Eq. (II) is a local master equation (whose validity is well
justified irrespective of the magnitude of the interaction between system
oscillators, see Ref. De Chiara  _et al._ (2018)) which, for a suitable choice
of the parameters, exhibits the presence of exceptional points Scali  _et al._
(2021).

Following Ref. De Chiara  _et al._ (2018), we obtain the heat current to the
bath $\alpha$,

$\displaystyle\dot{Q}_{\alpha}(t)$ $\displaystyle=$
$\displaystyle\lim_{\tau\rightarrow
0}\frac{1}{\tau}\text{Tr}\left[H_{B,\alpha}\left(\rho_{B,\alpha}(t+\tau)-\rho_{B,\alpha}(t)\right)\right]$
(12) $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}g^{2}_{i,\alpha}\Omega_{\alpha}\left[n_{\alpha}-\langle
a_{i}^{\dagger}a_{i}\rangle\right],$

and the power done on the system,

$\displaystyle\dot{W}(t)$ $\displaystyle=$ $\displaystyle\lim_{\tau\rightarrow
0}\frac{1}{\tau}\text{Tr}\left[H_{S}(t+\tau)\rho_{S}(t+\tau)-H_{S}(t)\rho_{S}(t)\right]$
(13) $\displaystyle-$ $\displaystyle\sum_{\alpha}\dot{Q}_{\alpha}(t),$

giving

$\displaystyle\dot{W}(t)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\sum_{\alpha=1}^{N_{B}}g^{2}_{i,\alpha}[\omega_{i}(t)-\Omega_{\alpha}]\left[n_{\alpha}-\langle
a_{i}^{\dagger}a_{i}\rangle\right]+$ (14) $\displaystyle+$
$\displaystyle\sum_{i=1}^{N}\sum_{\alpha=1}^{N_{B}}\dot{\omega}_{i}(t)\frac{2\expectationvalue{a_{i}^{\dagger}a_{i}}+1}{2},$

where $\langle\cdot\rangle$ denotes expectation value with respect to the
state $\rho_{S}(t)$ of the system at time $t$. The heat currents and the power
fulfill the first law of thermodynamics

$\dot{U}_{S}(t)=\dot{W}(t)+\sum_{\alpha=1}^{N_{B}}\dot{Q}_{\alpha}(t),$ (15)

where we have defined the internal energy as

$U_{S}(t)={\rm Tr}\left[H_{S}(t)\rho_{S}(t)\right].$ (16)

Notice that we follow the convention that work done on the system and heat
flowing into the system are positive.

Since we are dealing with linearly coupled harmonic oscillators, we find it
convenient to describe the evolution of the system using the formalism of
Gaussian continuous variable systems Braunstein and van Loock (2005). To this
end, we define a vector $R$ containing the system’s position and momentum
operators,

$R=\\{x_{1},p_{1},x_{2},p_{2},\dots x_{N},p_{N}\\}^{T}.$ (17)

In terms of $R$, we can express the system’s Hamiltonian as a quadratic form,

$H_{S}(t)=\frac{1}{2}RA_{S}(t)R^{T},$ (18)

with

$A_{S}(t)=B(t)+C(t),$ (19)

where

$\displaystyle B(t)=\bigoplus_{i=1}^{N}\begin{pmatrix}\omega_{i}^{2}(t)&0\\\
0&1\end{pmatrix},$ $\displaystyle
C(t)=\begin{pmatrix}0&C_{12}&\ldots&C_{1N}\\\ C_{21}&0&&\vdots\\\
\vdots&&&\vdots\\\ C_{N1}&C_{N2}&\ldots&0\end{pmatrix},$ $\displaystyle
C_{ij}=\begin{pmatrix}\lambda_{i,j}\sqrt{\omega_{i}(t)\omega_{j}(t)}&0\\\
0&\lambda_{i,j}\sqrt{\frac{1}{\omega_{i}(t)\omega_{j}(t)}}\end{pmatrix}.$ (20)

The evolution of Gaussian states can be described completely by the first and
second moments of $R$. The first moments fulfill the linear differential
equation

$\dot{\langle R\rangle}=D(t)\langle R\rangle,$ (21)

where we have defined $D(t)=S_{N}A_{S}(t)-K$, being

$S_{N}=\bigoplus_{i=1}^{N}\left(\begin{array}[]{cc}0&1\\\
-1&0\end{array}\right)$ (22)

the symplectic matrix and

$K=\frac{1}{2}\bigoplus_{i=1}^{N}\sum_{\alpha=1}^{N_{B}}g^{2}_{i,\alpha}\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right)$ (23)

the dissipation term. The second moments can be described by the covariance
matrix elements defined as

$\sigma_{kl}(t)=\frac{1}{2}\expectationvalue{R_{k}R_{l}+R_{l}R_{k}}-\expectationvalue{R_{k}}\expectationvalue{R_{l}}.$
(24)

Using the master equation in Eq. (II), it is possible to show that the
covariance matrix fulfills the equation (see, for example, Ref. Nicacio  _et
al._ (2016)),

$\dot{\sigma}=D\sigma+\sigma D^{T}+T,$ (25)

where we have defined the noise matrix,

$T=\frac{1}{2}\bigoplus_{i=1}^{N}\sum_{\alpha=1}^{N_{B}}g^{2}_{i,\alpha}(2n_{\alpha}+1)\left(\begin{array}[]{cc}\frac{1}{\omega_{i}(t)}&0\\\
0&\omega_{i}(t)\end{array}\right).$ (26)

|
---|---
|
Figure 2: The heat flows, (a) and (b), and power, (c) and (d) for the one-
oscillator system found by solving the master equation from Eq. (II)
numerically (superscript $n$) versus the “limit cycle” solutions (superscript
lc). We compare slow, (a) and (c), and fast driving, (b) and (d). The initial
conditions used are given in Eq. (III.1). Parameters:
$\delta\omega=\omega_{0}/2,\;g_{c}=g_{h}=\sqrt{\omega_{0}/2},\;\beta_{c}=10/\omega_{0},\;\beta_{h}=5/\omega_{0},\;\Omega_{c}=\omega_{0}/5,\;\Omega_{h}=\omega_{0}/2$.

##  III One-Oscillator System

In this section, we showcase the results for the simplest example of one
system oscillator, with driving frequency,

$\omega(t)=\omega_{0}+\delta\omega\sin(\theta t),$ (27)

in contact with two baths at different temperatures ($N=1,N_{B}=2$). For
clarity of notation, we assume bath 1 to be the cold bath and denote its
quantities with the subscript $c$, for example $n_{1}=n_{c}$ and
$g_{1,1}=g_{c}$. Similarly, we assume bath 2 to be the hot bath and denote its
quantities with the subscript $h$, for example, $n_{2}=n_{h}$ and
$g_{1,2}=g_{h}$.

In this case, the master equation takes the simple form,

$\displaystyle\dot{\rho}_{S}(t)$ $\displaystyle=$
$\displaystyle-\mathbbm{i}\commutator{H_{S}(t)}{\rho_{S}(t)}$ (28)
$\displaystyle+$
$\displaystyle\gamma\bar{n}\left(a_{1}^{\dagger}\rho_{S}a_{1}-\frac{1}{2}\left\\{a_{1}a_{1}^{\dagger},\rho_{S}\right\\}\right)$
$\displaystyle+$
$\displaystyle\gamma(\bar{n}+1)\left(a_{1}\rho_{S}a_{1}^{\dagger}-\frac{1}{2}\left\\{a_{1}^{\dagger}a_{1},\rho_{S}\right\\}\right),$

where $\gamma=g^{2}_{c}+g^{2}_{h}$ is the total dissipation rate and

$\bar{n}=\frac{g^{2}_{c}n_{c}+g^{2}_{h}n_{h}}{g^{2}_{c}+g^{2}_{h}}$ (29)

is the weighted average thermal occupation. Equation (28) shows that the
interaction of a single oscillator with two thermal baths is equivalent to
that of one oscillator in contact with an effective bath with a weighted
average temperature. An analogous result holds in the case of a larger number
of baths. While this effective picture is true for the system dynamics, the
individual heat currents $\dot{Q}_{\alpha}$ depend on the specific
temperatures and couplings of the baths and are given by Eq. (12).

Figure 3: Possible functionings of a single oscillator system where the bar
represents the average over a period of oscillation. (a) Engine, (b)
refrigerator, and (c) accelerator. Below each scheme the conditions for the
operating regime are shown.

Similarly, using the formalism introduced in Sec. II, we can obtain the
equivalent dynamical equation of a system made of one driven QHO in terms of
the covariance matrix [see Eq. (25)],

${\sigma}(t)=\begin{pmatrix}{\sigma}_{xx}(t)&{\sigma}_{xp}(t)\\\
{\sigma}_{xp}(t)&{\sigma}_{pp}(t)\end{pmatrix},$ (30)

where

$\begin{split}\dot{\sigma}_{xx}(t)&=-\gamma\sigma_{xx}(t)+\frac{\gamma(2\bar{n}+1)}{2\omega(t)}+2\sigma_{xp}(t),\\\
\dot{\sigma}_{xp}(t)&=-\gamma\sigma_{xp}(t)+\sigma_{pp}(t)-\sigma_{xx}(t)\omega^{2}(t),\\\
\dot{\sigma}_{pp}(t)&=-\gamma\sigma_{pp}(t)+\frac{\gamma(2\bar{n}+1)\omega(t)}{2}-2\sigma_{xp}(t)\omega^{2}(t).\end{split}$
(31)

From Eq. (12), the heat flows can be calculated as

$\dot{Q}_{\alpha}^{(1)}(t)\\!=\\!-\frac{g_{\alpha}^{2}\Omega_{\alpha}}{2\omega(t)}\\!\left[\sigma_{pp}(t)\\!+\\!\omega^{2}(t)\sigma_{xx}(t)\\!-\\!\omega(t)(2n_{\alpha}+1)\right],$
(32)

where $\alpha=c,h$ and the power is

$\displaystyle\dot{W}^{(1)}(t)$ $\displaystyle=$
$\displaystyle-\frac{\gamma}{2}\sigma_{pp}(t)+\omega(t)\dot{\omega}(t)\sigma_{xx}(t)$
(33) $\displaystyle+$
$\displaystyle\frac{\gamma\omega(t)}{2}[2\bar{n}+1-\omega(t)\sigma_{xx}(t)]$
$\displaystyle-$ $\displaystyle\dot{Q}_{c}^{(1)}(t)-\dot{Q}_{h}^{(1)}(t),$

containing a term explicitly depending on the derivative $\dot{\omega}(t)$ of
the oscillator’s frequency.

In the rest of this section, we are going to show a few examples of the
thermodynamics of one oscillator.

###  III.1 Slow driving

Let us first look at the slow driving case where the driving frequency is
changing slowly compared to the typical timescales of the system,
$\theta_{1}\ll\\{\gamma,\omega\\}$. In this case, we can set
$\dot{\sigma}_{S}(t)=0$, assuming the system to instantaneously follow the
steady state imposed by the instantaneous $\omega(t)$, and solve for
$\sigma_{xx}(t),\sigma_{xp}(t)$, and $\sigma_{pp}(t)$. We find

$\displaystyle\sigma_{xx}(t)$ $\displaystyle=$
$\displaystyle\frac{2\bar{n}+1}{2\omega(t)},$ $\displaystyle\sigma_{xp}(t)$
$\displaystyle=$ $\displaystyle 0,$ $\displaystyle\sigma_{pp}(t)$
$\displaystyle=$ $\displaystyle\omega(t)\frac{2\bar{n}+1}{2}.$ (34)

This state corresponds to a thermal state of a QHO with frequency $\omega(t)$
in equilibrium with an environment with thermal occupation $\bar{n}$. Plugging
these solutions into Eqs. (32) and (33) for the heat flows and power, we get

$\displaystyle\dot{Q}_{c}^{(1)}(t)$ $\displaystyle=$
$\displaystyle\frac{\Omega_{c}g_{c}^{2}g_{h}^{2}(n_{c}-n_{h})}{g_{c}^{2}+g_{h}^{2}},$
$\displaystyle\dot{Q}_{h}^{(1)}(t)$ $\displaystyle=$
$\displaystyle\frac{\Omega_{h}g_{c}^{2}g_{h}^{2}(n_{h}-n_{c})}{g_{c}^{2}+g_{h}^{2}},$
$\displaystyle\dot{W}^{(1)}(t)$ $\displaystyle=$
$\displaystyle\frac{g_{c}^{2}g_{h}^{2}(\Omega_{c}-\Omega_{h})(n_{h}-n_{c})}{g_{c}^{2}+g_{h}^{2}}+\frac{2\bar{n}+1}{2}\dot{\omega}(t).$

Notice that under the assumption of slow driving, the heat currents are
constant. This is due to the fact that instantaneously the system is in the
steady state corresponding to the instantaneous value of the frequency
$\omega(t)$. Note that the second term in $\dot{W}(t)$ is proportional to
$\dot{\omega}(t)$ and therefore, since $\omega(t)$ is periodic, it does not
contribute to the average value.

Figure 2 compares the approximate “limit cycle” solutions (III.1) with the
exact numerical results from Eqs. (32) and (33). Choosing a modulation
frequency $\theta$ small compared to $\omega$ and $\gamma$, as in the top and
bottom left panels, we see that after a brief transient period that depends on
the initial conditions, the heat flows and power still oscillate slightly, but
tend to the limit cycle solutions.

In Fig. 2, we have assumed the system to be initially in a Gibbs state in
thermal equilibrium at the cold inverse temperature $\beta_{c}$ and initial
frequency $\omega(0)$, so that we have the initial conditions,

$\displaystyle\sigma_{xx}(0)$ $\displaystyle=$
$\displaystyle\frac{2n_{0}+1}{2\omega(0)},$ $\displaystyle\sigma_{xp}(0)$
$\displaystyle=$ $\displaystyle 0,$ $\displaystyle\sigma_{pp}(0)$
$\displaystyle=$ $\displaystyle\frac{\omega(0)(2n_{0}+1)}{2},$ (36)

where

$n_{0}=\frac{1}{e^{\beta_{c}\omega(0)}-1}.$ (37)

The slow driving solutions of Eq. (III.1) allow us to interpret the
functioning of the single oscillator system as a thermal device. Depending on
the parameters $n_{c},n_{h},\Omega_{c}$, and $\Omega_{h}$, the device can
function as an engine, a refrigerator or an accelerator. In the following we
briefly illustrate these functionings.

####  III.1.1 Engine

To create an engine we need heat to flow from the hot bath, through the system
and to the cold bath and in doing this, the system should have a certain power
output. In order to achieve this functioning we require that
$\overline{\dot{Q}^{(1)}_{c}}<0,\;\overline{\dot{Q}^{(1)}_{h}}>0$, and
$\overline{\dot{W}^{(1)}}<0$, where the bar indicates the average over a
period of oscillation. The corresponding operating mode is shown in Fig. 3.
From Eq. (III.1), this happens when $n_{c}<n_{h}$, i.e., when
$\beta_{h}\Omega_{h}<\beta_{c}\Omega_{c}$, and $\Omega_{c}<\Omega_{h}$.

####  III.1.2 Refrigerator

For a refrigerator, we instead need a power input that enables heat to flow
from the cold bath, through the system and to the hot bath. This functioning
on average over a period of oscillation requires
$\overline{\dot{Q}^{(1)}_{c}}>0,\;\overline{\dot{Q}^{(1)}_{h}}<0$, and
$\overline{\dot{W}^{(1)}}>0$ which, from Eq. (III.1), happens when
$n_{h}<n_{c}$, i.e., when $\beta_{h}\Omega_{h}>\beta_{c}\Omega_{c}$, and
$\Omega_{c}<\Omega_{h}$ as shown in Fig. 3.

####  III.1.3 Accelerator

An accelerator occurs when, like an engine, we have heat flowing from the hot
bath to the system and from the system to the cold bath, i.e.,
$\overline{\dot{Q}^{(1)}_{c}}<0$ and $\overline{\dot{Q}^{(1)}_{h}}>0$, the
process being sped up, on average, by inputting power over a period of
oscillation, $\overline{\dot{W}^{(1)}}>0$, as shown in Fig. 3. To achieve this
regime in Eq. (III.1), we require that $n_{c}<n_{h}$, meaning
$\beta_{c}\Omega_{c}>\beta_{h}\Omega_{h}$, and $\Omega_{h}<\Omega_{c}$.

During a cycle, the signs of the thermodynamic quantities may change. For
instance, while the average power over a cycle may be negative, causing the
system to operate as an engine overall, there could be instants during the
cycle when the power becomes positive so that the device temporarily switches
to an accelerator.

Figure 4: Plot of the power output against the average efficiency for the one-
oscillator system, both quantities averaged over a period. The dashed line
corresponds to the maximum power which occurs at the Curzon-Ahlborn condition
when $\Omega_{h}=\sqrt{\frac{\beta_{c}}{\beta_{h}}}\Omega_{c}$. Other
parameters:
$g_{c}=g_{h}=\sqrt{\omega_{0}/2},\;\beta_{c}=10/\omega_{0},\;\beta_{h}=5/\omega_{0},\;\Omega_{c}=\omega_{0}/10$.
Figure 5:  COP of a refrigerator for the one-oscillator system, its average
over a period, and the Otto COP value plotted against time. Parameters:
$\delta\omega=\omega_{0}/2,\;\theta/\omega_{0}=\pi/20,\;g_{c}=g_{h}=\sqrt{\omega_{0}/2},\;\beta_{c}=10/\omega_{0},\;\beta_{h}=5/\omega_{0},\;\Omega_{c}=\omega_{0}/10,\;\Omega_{h}=2\omega_{0}$.

####  III.1.4 Efficiency and coefficient of performance

The instantaneous efficiency of the engine can be defined as

$\eta^{(1)}(t)=-\frac{\dot{W}^{(1)}(t)}{\dot{Q}^{(1)}_{h}(t)},$ (38)

with the condition $\dot{W}^{(1)}(t)<0$. Plugging in the limit cycle
solutions, we get

$\eta^{(1)}(t)=\eta_{\rm Otto}+\delta\eta^{(1)}(t),$ (39)

where we have defined the efficiency of an Otto cycle operating between
frequencies $\Omega_{c}$ and $\Omega_{h}$ as $\eta_{\rm
Otto}=1-\Omega_{c}/\Omega_{h}$ and the time-dependent correction,

$\delta\eta^{(1)}(t)=\frac{\dot{\omega}(t)(2\bar{n}+1)}{2\Omega_{h}(n_{h}-n_{c})}\left(\frac{1}{g_{c}^{2}}+\frac{1}{g_{h}^{2}}\right).$
(40)

Since $\omega(t)$ is a periodic function, the average of $\dot{\omega}$ over a
period is zero and we obtain that the average efficiency over a period reduces
to the Otto efficiency, i.e., $\overline{\eta^{(1)}}=\eta_{\rm Otto}$ since
$\overline{\delta\eta^{(1)}}=0$.

For an engine, we are often interested in a regime in which its power output
is maximum and not when its efficiency is the highest. Since the average
efficiency $\overline{\eta^{(1)}}=\eta_{\rm Otto}$ is equal to the Otto value,
we expect the average efficiency at maximum power within these assumptions to
be equal to the Curzon-Ahlborn value Curzon and Ahlborn (1975),

$\eta_{CA}=1-\sqrt{\frac{\beta_{h}}{\beta_{c}}}.$ (41)

This is true for a system of one QHO and occurs when
$\Omega_{h}=\sqrt{\frac{\beta_{c}}{\beta_{h}}}\Omega_{c}$. For one possible
choice of parameters this is illustrated in Fig. 4.

Figure 6:  Plots of average heat flows and power for the one-oscillator system
as a function of $\Omega_{h}$. (i) Accelerator, (ii) engine, (iii)
refrigerator, and (iv) dissipator. (a) Slow driving:
$\theta/\omega_{0}=\pi/200$. (b) Increased driving: $\theta/\omega_{0}=\pi/4$.
(c) Fast driving: $\theta/\omega_{0}=\pi/2$. The initial conditions are given
in Eq. (III.1). Parameters:
$\delta\omega=\omega_{0}/2,\;g_{c}=g_{h}=\sqrt{\omega_{0}}/2,\;\beta_{c}=20/\omega_{0},\;\beta_{h}=10/\omega_{0},\;\Omega_{c}=\omega_{0}/20$.

In the refrigerator regime, the functioning of the thermal device can be
assessed by the coefficient of performance (COP) of the refrigerator defined
as

$COP^{(1)}(t)=\frac{\dot{Q}^{(1)}_{c}(t)}{\dot{W}^{(1)}(t)},$ (42)

with $\dot{Q}^{(1)}_{c}(t)>0$. Plugging in the limit cycle solutions, we get

$\frac{1}{COP^{(1)}(t)}=\frac{1}{COP_{\rm Otto}}+\frac{1}{\delta
COP^{(1)}(t)},$ (43)

where we have defined $COP_{\rm Otto}=\Omega_{c}/(\Omega_{h}-\Omega_{c})$ and

$\frac{1}{\delta
COP^{(1)}(t)}=-\frac{\Omega_{h}}{\Omega_{c}}\delta\eta^{(1)}(t).$ (44)

Since $\delta COP^{(1)}(t)$ oscillates, between positive and negative values,
the instantaneous COP can be larger than the Otto value $COP_{\rm Otto}$. Even
more, due to the non-linear relation in which $\delta COP^{(1)}(t)$ enters the
expression of the COP, the average $\overline{COP^{(1)}}$ can be larger than
$COP_{\rm Otto}$ as shown in Fig. 5. On the other hand, the ratio of the
average of the heat flow to the cold bath and the average of power is equal to
the Otto value $COP_{\rm Otto}$. Note that in Fig. 5 the parameters have been
chosen so that the power remains positive throughout each cycle allowing the
system to operate as a refrigerator at all times. These results show that,
although on average (when taking the ratio of the averages of heat flow and
power independently) the refrigerator with modulated driving performs as
efficiently as an Otto refrigerator, the instantaneous coefficient of
performance can in fact surpass it.

###  III.2 Fast driving

We now analyze what happens in the fast regime when the driving modulation
frequency $\theta$ is comparable to the typical timescale of the system as
dictated by $\gamma$ or $\omega_{0}$. In this regime, we can no longer assume
$\dot{\sigma}_{S}(t)=0$, so we must solve the master equation directly by
numerical integration. We then use the resulting values of the covariance
matrix to obtain the relevant thermodynamic quantities.

After an initial transient period, we find that the heat currents always tend
to some limit cycle in which the functions become periodic even though not
with a simple sinusoidal dependence, as seen in the top and bottom-right panel
of Fig. 2. This is the result of the interplay of the oscillator’s free
evolution, the modulation of its frequency, and the dissipation induced by the
baths. We integrate over a period to find the average values of heat flows and
power, which are plotted in Fig. 6 as a function of $\Omega_{h}$ for
increasing values of $\theta$.

For slow driving ($\theta/\omega_{0}=\pi/200$) there is a value of
$\Omega_{h}$ at which the heat flows and power all go to zero. This is the
Carnot point and occurs when
$\frac{\beta_{c}}{\beta_{h}}=\frac{\Omega_{h}}{\Omega_{c}}$. We find that the
average dissipated power increases with faster driving while the cooling power
decreases, meaning that the COP for the refrigerator decreases. Also, with
faster driving the power remains positive for all values of $\Omega_{c}$ and
$\Omega_{h}$ so it is no longer possible for the system to operate as an
engine. In particular, for larger $\theta$ we do find that the power remains
positive while both heat flows become negative. This regime corresponds to a
dissipator: the increasingly fast driving results in a large injection of
power into the system and is transformed into heat which is dissipated in both
the cold and hot reservoirs. A dissipator may be used to dissipate heat
generated in a mechanical or electronic device. Alternatively, a dissipator
can be used to warm another system, hence the name heater. The maxima of the
dissipated heats $|\dot{Q}^{(1)}_{c}|$ and $|\dot{Q}^{(1)}_{h}|$ correspond to
the maximum of the injected power $\dot{W}^{(1)}$ and arise for a “resonant”
value of $\theta$ for which the conversion work to heat induced by the driving
is more effective. This is seen in Fig. 7. Notice that as the process is
dissipative, there is no thermodynamic constraint limiting the amount of heat
produced.

Figure 7: Average power over a period for the one-oscillator system as a
function of the modulation frequency $\theta$. There is a particular value of
$\theta$ ($\theta/\omega_{0}\approx 2.12$, indicated by a vertical dashed
line) for which the average power injected over a cycle reaches a maximum.
Parameters:
$\delta\omega=\omega_{0}/2,\;g_{c}=g_{h}=\sqrt{\omega_{0}}/2,\;\beta_{c}=20/\omega_{0},\;\beta_{h}=10/\omega_{0},\;\Omega_{c}=\omega_{0}/20,\;\Omega_{h}=3\omega_{0}/20$.

##  IV Two Oscillators

We now extend our study to a system of two driven QHO’s ($N=2$), with driving
frequencies $\omega_{1}(t)$ and $\omega_{2}(t)$, respectively, described by
the Hamiltonian of Eq. (1) with $\lambda_{12}=\lambda$. We assume each QHO to
be coupled to only one reservoir ($N_{B}=2$): oscillator 1 (2) is coupled to a
cold (hot) reservoir characterized by frequency $\Omega_{c}\;(\Omega_{h})$ and
thermal occupation $n_{c}\;(n_{h})$. It follows that in Eq. (7),
$g_{1,2}=g_{2,1}=0$, while we set $g_{c}=g_{1,1}$ and $g_{h}=g_{2,2}$. Our
system can now be described by a 4 by 4 covariance matrix defined in Eq. (24)
with matrix entries

${\sigma}_{S}(t)=\begin{pmatrix}{\sigma}_{x_{1}x_{1}}(t)&{\sigma}_{x_{1}p_{1}}(t)&{\sigma}_{x_{1}x_{2}}(t)&{\sigma}_{x_{1}p_{2}}(t)\\\
{\sigma}_{x_{1}p_{1}}(t)&{\sigma}_{p_{1}p_{1}}(t)&{\sigma}_{x_{2}p_{1}}(t)&{\sigma}_{p_{1}p_{2}}(t)\\\
{\sigma}_{x_{1}x_{2}}(t)&{\sigma}_{x_{2}p_{1}}(t)&{\sigma}_{x_{2}x_{2}}(t)&{\sigma}_{x_{2}p_{2}}(t)\\\
{\sigma}_{x_{1}p_{2}}(t)&{\sigma}_{p_{1}p_{2}}(t)&{\sigma}_{x_{2}p_{2}}(t)&{\sigma}_{p_{2}p_{2}}(t)\end{pmatrix}.$
(45)

Just as in Sec. III, we can find the heat currents

$\displaystyle\dot{Q}^{(2)}_{c}$ $\displaystyle=$
$\displaystyle-\frac{g_{c}^{2}\Omega_{c}}{2\omega_{1}}\left[\sigma_{p_{1}p_{1}}-\omega_{1}(1+2n_{c})+\sigma_{x_{1}x_{1}}\omega^{2}_{1}\right],$
(46) $\displaystyle\dot{Q}^{(2)}_{h}$ $\displaystyle=$
$\displaystyle-\frac{g_{h}^{2}\Omega_{h}}{2\omega_{2}}\left[(\sigma_{p_{2}p_{2}}-\omega_{2}(1+2n_{h})+\sigma_{x_{2}x_{2}}\omega^{2}_{2}\right],$
(47)

which are extensions of Eq. (32). We have omitted the explicit time-
dependence. The expression for the power, corresponding to Eq. (14), in terms
of the covariance matrix entries, is quite lengthy and reported in Appendix A.
As a check of the calculations, we have verified the first law of
thermodynamics expressed in Eq. (15).

Figure 8:  COP of a refrigerator made of two QHO’s in the slow driving regime
plotted as a function of time (solid line) compared with its average over a
period (short-dashed) and the Otto value (long-dashed). Parameters:
$\delta\omega=\omega_{0}/2,\;\theta/\omega_{0}=\pi/20,\;g_{c}=g_{h}=\sqrt{\omega_{0}/2},\;\lambda=\sqrt{2}\omega_{0},\;\beta_{c}=10/\omega_{0},\;\beta_{h}=5/\omega_{0},\;\Omega_{c}=\omega_{0}/10,\;\Omega_{h}=2\omega_{0}$.

###  IV.1 Slow driving

In the slow driving case, solving $\dot{\sigma}_{S}(t)=0$ we can find
analytically the solution for the time-dependence of the covariance matrix in
the limit cycle. The explicit form is not reported as it is quite lengthy.
Using this expression, we obtain the heat currents,

$\begin{split}\dot{Q}^{(2)}_{c}&=\frac{4\Omega_{c}\lambda^{2}g_{c}^{2}g_{h}^{2}\left(g_{c}^{2}+g_{h}^{2}\right)\left(n_{c}-n_{h}\right)}{\left(g_{c}^{2}+g_{h}^{2}\right)^{2}\left(4\lambda^{2}+g_{c}^{2}g_{h}^{2}\right)+4g_{c}^{2}g_{h}^{2}[\omega_{1}-\omega_{2}]^{2}},\\\
\dot{Q}^{(2)}_{h}&=-\frac{4\Omega_{h}\lambda^{2}g_{c}^{2}g_{h}^{2}\left(g_{c}^{2}+g_{h}^{2}\right)\left(n_{c}-n_{h}\right)}{\left(g_{c}^{2}+g_{h}^{2}\right)^{2}\left(4\lambda^{2}+g_{c}^{2}g_{h}^{2}\right)+4g_{c}^{2}g_{h}^{2}[\omega_{1}-\omega_{2}]^{2}}.\end{split}$
(48)

Again, the expression for the power is quite lengthy and reported in Appendix
A.

Notice that if $\omega_{1}(t)\neq\omega_{2}(t)$, due to the term proportional
to $[\omega_{1}(t)-\omega_{2}(t)]^{2}$ in the denominator, the values of
$\dot{Q}^{(2)}_{c}(t)$ and $\dot{Q}^{(2)}_{h}(t)$ [and indeed also of
$\dot{W}^{(2)}(t)$] always have a smaller magnitude than when we have
$\omega_{1}(t)=\omega_{2}(t)$. Therefore, from now on, we will assume that
both QHO’s have the same frequency $\omega_{1}(t)=\omega_{2}(t)=\omega(t)$
(and hence $\theta_{1}=\theta_{2}=\theta$). Under these assumptions, the
expressions of heat currents simplify to

$\displaystyle\dot{Q}^{(2)}_{c}=\frac{\Omega_{c}g_{c}^{2}g_{h}^{2}\left(n_{c}-n_{h}\right)}{\left(g_{c}^{2}+g_{h}^{2}\right)\left(1+\frac{g_{c}^{2}g_{h}^{2}}{4\lambda^{2}}\right)}=\frac{\dot{Q}_{c}^{(1)}}{1+\frac{g_{c}^{2}g_{h}^{2}}{4\lambda^{2}}},$
(49)
$\displaystyle\dot{Q}_{h}^{(2)}=-\frac{\Omega_{h}g_{c}^{2}g_{h}^{2}\left(n_{c}-n_{h}\right)}{\left(g_{c}^{2}+g_{h}^{2}\right)\left(1+\frac{g_{c}^{2}g_{h}^{2}}{4\lambda^{2}}\right)}=\frac{\dot{Q}_{h}^{(1)}}{1+\frac{g_{c}^{2}g_{h}^{2}}{4\lambda^{2}}},$
(50)

where in the last equations we have used the expression for the single
oscillator heat currents of Eq. (III.1). When the working substance is made of
more than two linearly coupled resonant oscillators, with the same inter-
oscillator coupling, the heat flows and the non-oscillating part of the power
coincide with the case with only two oscillators as in Eqs. (49) and (50), see
Ref. Asadian  _et al._ (2013).

Figure 9: Plots of the heat flows (a) and power (b) for the two-oscillator
system with fast driving (solid lines) compared to their values in the slow
driving regime (dashed lines). Both oscillators have the same initial
conditions as in Eq. (III.1). Parameters:
$\delta\omega=\omega_{0}/2,\;\theta/\omega_{0}=\pi,\;g_{c}=g_{h}=\sqrt{\omega_{0}/2},\;\lambda=\sqrt{2}\omega_{0},\;\beta_{c}=10/\omega_{0},\;\beta_{h}=5/\omega_{0},\;\Omega_{c}=\omega_{0}/10,\;\Omega_{h}=3\omega_{0}/10$.

Figure 10: Plots of average heat flows and power for the two-oscillator system
as a function of $\Omega_{h}$. (i) Accelerator, (ii) engine, (iii)
refrigerator, and (iv) dissipator. (a) Slow driving:
$\theta/\omega_{0}=\pi/200$. (b) Increased driving: $\theta/\omega_{0}=\pi/4$.
(c) Fast driving: $\theta/\omega_{0}=\pi/2$. The initial conditions are given
in Eq. (III.1). Parameters:
$\delta\omega=\omega_{0}/2,\;g_{c}=g_{h}=\sqrt{\omega_{0}}/2,\;\lambda=\omega_{0}/\sqrt{2},\;\beta_{c}=20/\omega_{0},\;\beta_{h}=10/\omega_{0},\;\Omega_{c}=\omega_{0}/20$.

As the coupling strength, $\lambda$, between the two system oscillators
increases, the heat flows and power tend to the same expressions for the one-
oscillator system,

$\displaystyle\lim_{\lambda\rightarrow\infty}\dot{Q}^{(2)}_{c}(t)$
$\displaystyle=$
$\displaystyle\frac{\Omega_{c}g_{c}^{2}g_{h}^{2}(n_{c}-n_{h})}{g_{c}^{2}+g_{h}^{2}}=\dot{Q}_{c}^{(1)}(t),$
(51) $\displaystyle\lim_{\lambda\rightarrow\infty}\dot{Q}^{(2)}_{h}(t)$
$\displaystyle=$
$\displaystyle\frac{\Omega_{h}g_{c}^{2}g_{h}^{2}(n_{h}-n_{c})}{g_{c}^{2}+g_{h}^{2}}=\dot{Q}_{h}^{(1)}(t).$
(52)

This result can be understood because in this regime, the two oscillators are
so strongly coupled that only the center of mass effectively interacts with
the two baths and is affected by the driving. The relative motion mode, on the
other hand, stiffens, being characterized by a very large effective frequency
which increases with $\lambda$. In this limit, the expression for the power
also reduces to the corresponding expression for one oscillator,

$\displaystyle\lim_{\lambda\rightarrow\infty}\dot{W}^{(2)}(t)$
$\displaystyle=$
$\displaystyle\frac{g_{c}^{2}g_{h}^{2}(\Omega_{c}-\Omega_{h})(n_{h}-n_{c})}{(g_{c}^{2}+g_{h}^{2})}$
(53) $\displaystyle+$
$\displaystyle\frac{\left[g_{c}^{2}(2n_{c}+1)+g_{h}^{2}(2n_{h}+1)\right]\dot{\omega}(t)}{2(g_{c}^{2}+g_{h}^{2})}$
$\displaystyle=$ $\displaystyle\dot{W}^{(1)}(t).$

Under the assumption $\omega_{1}(t)=\omega_{2}(t)$, using the approach
reported in Appendix A, it is easy to see that the conditions required for the
system to operate on average as an engine, a refrigerator or an accelerator
are the same as those in Fig. 3 for any value of $\lambda$.

One may use Eq. (38) to show that the efficiency for this system of two
quantum harmonic oscillators is

$\eta^{(2)}(t)=\eta_{\text{Otto}}+\delta\eta^{(2)}(t),$ (54)

where

$\delta\eta^{(2)}(t)=\frac{\gamma(1+n_{c}+n_{h})}{4\lambda^{2}\Omega_{h}(n_{c}-n_{h})}\dot{\omega}(t)+2\delta\eta^{(1)}(t).$
(55)

Since $\delta\eta^{(1)}(t)$, and hence $\delta\eta^{(2)}(t)$, is proportional
to $\dot{\omega}(t)$, the average efficiency over a period reduces to
$\overline{\eta^{(2)}}=\eta_{\text{Otto}}$, as in the case of one QHO.

The COP of a refrigerator may be found from Eq. (42) satisfying

$\frac{1}{COP^{(2)}(t)}=\frac{1}{COP_{\text{Otto}}}+\frac{1}{\delta
COP^{(2)}(t)},$ (56)

where

$\frac{1}{\delta
COP^{(2)}(t)}=-\frac{\Omega_{h}}{\Omega_{c}}\delta\eta^{(2)}(t).$ (57)

In Fig. 8 the parameters defining the conditions of the baths, driving of the
system, and couplings of the system to the bath are the same as in Fig. 5,
with the addition of the coupling between the two system oscillators,
$\lambda$. As in the case of Fig. 5, here the system operates as a
refrigerator at all times. Comparing these two plots we see that the
instantaneous COP reaches a higher peak in Fig. 8 and consequently a higher
average over a cycle.

###  IV.2 Fast driving

Let us now pass to the case of fast driving where the condition
$\dot{\sigma}_{S}=0$ does not hold. We continue to assume the resonant
condition $\omega_{1}(t)=\omega_{2}(t)=\omega(t)$. Solving the master equation
numerically we find that, as with the one QHO system, the heat flows and power
reach a limit cycle after a brief transient period, as shown in Fig. 9. Plots
of the average heat flows and power over a period and their corresponding
regimes are shown in Fig. 10. They look very similar to the corresponding
graphs for the system of one QHO in Fig. 6. However, note that for a two-QHO
system with fast driving ($\theta/\omega_{0}=\pi/2$) the power is larger while
the system is functioning as a dissipator.

##  V Squeezed Bath

In this section, we consider the extension of our model in which one of the
environments is composed of squeezed thermal ancillas. These can be prepared
by applying the single-mode squeezing operator defined as

$S(\xi)=\exp(\frac{1}{2}\xi^{*}b_{\alpha,j}^{2}-\frac{1}{2}\xi{b_{\alpha,j}^{\dagger}}^{2}),$
(58)

where $\xi=re^{i\phi}$ is a complex parameter which depends on two real
values: the squeezing magnitude $r\geq 0$ and the squeezing phase $\phi$. For
simplicity, we will look at the case where $\phi=0$, corresponding to the
squeezing of the position-like quadrature and anti-squeezing of the momentum-
like quadrature. As we will see shortly, squeezing effectively raises the
temperature of the bath, so for convenience we will only apply squeezing to
the hotter bath labeled 2. This leads to the following average occupation in
the thermal squeezed state:

$\displaystyle n_{h}^{(\text{eff})}$ $\displaystyle=$ $\displaystyle{\rm
Tr}\left[b^{\dagger}_{2,j}b_{2,j}S(\xi)\rho_{B,2}S^{\dagger}(\xi)\right]$ (59)
$\displaystyle=$ $\displaystyle\frac{(1+2n_{h})\cosh 2r-1}{2}.$

Using the fact that $\cosh 2r\geq 1$ we find $n_{h}^{(\text{eff})}\geq n_{h}$
with the equality obtained only for $r=0$.

For a system of two oscillators, we follow the same method as in Sec.
II—setting $\omega_{1}(t)=\omega_{2}(t)=\omega(t)$—to get the master equation
and energy flows. All the expressions for the heat currents and power are
formally identical to those used in Sec. IV with the replacement $n_{h}\to
n_{h}^{(\text{eff})}$.

As a consequence, at the thermodynamic level, the only modification that the
squeezing introduces is a higher temperature, and thus energy, in the
corresponding bath. This is connected to the known fact that a thermal state
minimizes the energy for fixed entropy. The squeezed thermal state is
characterized by the same entropy but higher energy corresponding to an
effective larger temperature.

To find the effective temperature we can use

$n_{h}^{(\text{eff})}=\Big{(}e^{\beta_{h}^{(\text{eff})}\Omega_{h}}-1\Big{)}^{-1},$
(60)

where $n_{h}^{(\text{eff})}$ is defined in Eq. (59). Solving for
$\beta_{h}^{(\text{eff})}$ we get

$\beta_{h}^{(\text{eff})}=\frac{1}{\Omega_{h}}\log[\frac{\tanh^{2}(r)+\exp(\beta_{h}\Omega_{h})}{1+\tanh^{2}(r)\exp(\beta_{h}\Omega_{h})}\Bigg{]},$
(61)

which coincides with the value found in Refs. Huang  _et al._ (2012); Correa
_et al._ (2015); Alicki and Gelbwaser-Klimovsky (2015); Latune  _et al._
(2019).

Figure 11: Plot of Otto, Carnot and effective Carnot efficiency, for squeezing
parameter
$r=0.3,\;\beta_{c}=10/\omega_{0},\;\beta_{h}=5/\omega_{0},\;\Omega_{c}=\omega_{0}$.

The fact that the squeezed reservoir has a higher effective temperature
modifies the ranges of the operating modes as a function of the frequencies
illustrated in Fig. 3. Let us start by looking at the case of an engine. The
condition for this kind of functioning now becomes

$\beta^{(\text{eff})}_{h}\Omega_{h}<\beta_{c}\Omega_{c}\quad\textrm{and}\quad\Omega_{c}<\Omega_{h},$
(62)

while the average efficiency for slow driving is still the Otto efficiency

${\overline{\eta^{(2)}}}=1-\frac{\Omega_{c}}{\Omega_{h}}=\eta_{\rm Otto}.$
(63)

Combining the conditions given in Eq. (62) we get

$\frac{\beta_{h}^{(\text{eff})}}{\beta_{c}}<\frac{\Omega_{c}}{\Omega_{h}}<1.$
(64)

By definition $\beta^{(\text{eff})}_{h}<\beta_{h}$, so squeezing decreases the
lower bound on $\Omega_{c}/\Omega_{h}$, therefore allowing a larger average
efficiency $\overline{\eta^{(2)}}$ to be attained. This, as also discussed in
Refs. Singh and Müstecaplıoğlu (2020) and Abah and Lutz (2014); Correa  _et
al._ (2015); Alicki and Gelbwaser-Klimovsky (2015); Manzano  _et al._ (2016);
Klaers  _et al._ (2017); Manzano (2018); Latune  _et al._ (2019); Huang  _et
al._ (2012), is not in contradiction with the Carnot bound which is derived
assuming equilibrium baths. Moreover, it does not violate the second law of
thermodynamics provided that one accounts for the preparation of squeezed
ancillas with additional resources.

In this spirit, it is fairer to define a modified Carnot efficiency as

$\eta_{C}^{(\text{eff})}=1-\frac{\beta_{h}^{(\text{eff})}}{\beta_{c}},$ (65)

such that $\overline{\eta^{(2)}}\leq\eta_{C}^{(\text{eff})}$. Figure 11 plots
the Otto efficiency obtained for the squeezed reservoir as well as the Carnot,
and the effective Carnot efficiency, as a function of $\Omega_{h}$ for a
squeezing parameter $r=0.3$. As expected, the efficiency $\eta_{\rm Otto}$ may
become larger than $\eta_{C}$ but never larger than $\eta_{C}^{(\text{eff})}$.

Now, let us move to the case of a refrigerator. The condition for this kind of
functioning now becomes

$\beta_{c}\Omega_{c}<\beta^{(\text{eff})}_{h}\Omega_{h}\quad\textrm{and}\quad\Omega_{c}<\Omega_{h},$
(66)

while the coefficient of performance for slow driving is still given by Eq.
(43). In order to have the largest possible COP, we want
$\Omega_{c}/\Omega_{h}$ to be as close to $1$ as possible. However, combining
the conditions in Eq. (66) gives

$\frac{\Omega_{c}}{\Omega_{h}}<\frac{\beta_{h}^{(\text{eff})}}{\beta_{c}}<1,$
(67)

meaning the upper bound on $\Omega_{c}/\Omega_{h}$ is decreased by squeezing
and therefore so is the COP. Summarizing, while squeezing improves the
performance of the two-oscillator system operating as an engine, it degrades
it when operating as a refrigerator.

The squeezed bath also allows the creation of steady-state entanglement
between the two oscillators that would otherwise be impossible to find in our
model with thermal baths. The absence of entanglement in the steady state for
thermal baths can be analytically proven in the slow driving limit under
symmetry conditions, e.g. equal system’s frequencies and equal couplings to
the respective baths (see Appendix B). We have also verified numerically that
there is no entanglement when these conditions are relaxed. This conclusion
refers to our specific model and does not preclude the existence of
entanglement in the steady state of quantum harmonic oscillators described by
a different model.

Figure 12 shows an example of this squeezing-induced entanglement as measured
by the logarithmic negativity (see Appendix B) for oscillators with constant
frequencies (no modulation, i.e., $\theta=0)$. Such entanglement is quite
fragile and requires the baths to be at very low temperatures. Notice that the
value of the logarithmic negativity increases from zero as we introduce
squeezing until a maximum is reached. From this point, increasing the
squeezing parameter, $r$, only decreases the logarithmic negativity until it
disappears for larger $r$. Choosing the value of $r$ which maximizes the
logarithmic negativity, Fig. 12 also shows that fast driving the frequency of
the oscillators can lead, with appropriately optimized parameters, to higher
instantaneous values of the logarithmic negativity compared to the slow
driving case with constant frequencies.

Summing up: for the thermodynamic quantities, the presence of squeezing only
raises the effective temperature of the squeezed bath and does not change
qualitatively the overall behavior; for the quantum correlations, instead,
squeezing has a pivotal role in the creation of entanglement.

Figure 12: (a) Steady-state logarithmic negativity $\mathcal{E}_{N}$ of the
two system’s oscillators with constant frequencies as a function of the
squeezing parameter $r$ of the hot environment. (b) Logarithmic negativity
$\mathcal{E}_{N}$ (solid line) as a function of time for $r=1.3$ when we drive
the two system oscillators. The dashed line is the maximum value attained with
constant frequencies. Parameters for both plots:
$\lambda=1,\beta_{c}=81.5/\lambda,\beta_{h}=40.7/\lambda,g_{c}=0.495\sqrt{\lambda},g_{h}=0.612\sqrt{\lambda},\omega_{0,1}=\Omega_{c}=4.07\lambda,\omega_{0,2}=\Omega_{h}=0.244\lambda,\delta\omega_{0,1}=\omega_{0,1}/10,\delta\omega_{0,2}=\omega_{0,2}/10,\theta=\pi\lambda$.

##  VI Conclusions

We have investigated the dissipative evolution of a network of driven QHO’s
and studied the energy flows in and out of the system in terms of heat
currents and power. A system composed of both one and two driven QHO’s coupled
to two heat baths of different temperatures may operate as an engine, a
refrigerator, or an accelerator when the system is parametrically driven
slowly. The efficiency of such an engine oscillates periodically in time
around the Otto efficiency. While the instantaneous coefficient of performance
of a refrigerator also oscillates, its average over a cycle can be found to be
above the corresponding value of an Otto refrigerator. A higher maximum value
and average over a cycle of the instantaneous COP has been found when we have
added a second driven oscillator to the system.

As we have increased the driving speed in both the one- and two-oscillator
systems, we have found that there is an increase in the power injected into
the system and it has been found that there is some resonant driving frequency
for which a maximum is reached.

We have also shown that applying squeezing to the hot bath effectively raises
its temperature and leads to a larger effective Carnot efficiency. We have
found that squeezing is a necessary requirement for the sustained generation
of entanglement within a system of two driven dissipative QHO’s. Controlled
driving of the system’s oscillators can increase the degree of entanglement.

Our findings can be readily verified in experimental platforms described by
linearly coupled quantum harmonic oscillators, e.g., trapped ions, atomic
ensembles, or photonic systems. The extension of our model to a large number
of oscillators and the corresponding scaling of thermodynamic quantities will
be the subject of further studies.

###### Acknowledgements.

G.D.C. acknowledges support by the UK EPSRC EP/S02994X/1. H.L. acknowledges
support from the UK EPSRC 2278075. N.P. acknowledges the “Université de
Franche- Comté” for financial support throughout the mobility grant “MOBILITÉ
INTERNATIONALE DES DOCTORANTS 2019” for his visit at the Queen’s University
Belfast (QUB) and thanks the QTeQ group and CTAMOP at QUB for their kind
hospitality on this occasion. B.B. acknowledges support by the French
“Investissements d’Avenir” program, Project ISITE-BFC (Contract No.
ANR-15-IDEX- 03).

## Author Declarations

### Conflict of Interest

The authors have no conflicts to disclose.

##  Appendix A Expression of the power for two driven oscillators

The expression for the work done on two oscillators in contact with two baths
is

$\displaystyle\dot{W}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{\lambda\left(\omega_{1}-\omega_{2}\right)}{\sqrt{\omega_{1}\omega_{2}}}\left(\omega_{1}\sigma_{x_{1}p_{2}}-\omega_{2}\sigma_{x_{2}p_{1}}\right)$
(68) $\displaystyle+$
$\displaystyle\frac{g_{c}^{2}}{2}\left[\omega_{1}\left(1+2n_{c}-\omega_{1}\sigma_{x_{1}x_{1}}\right)-\sigma_{p_{1}p_{1}}\right]$
$\displaystyle+$
$\displaystyle\frac{g_{h}^{2}}{2}\left[\omega_{2}\left(1+2n_{h}-\omega_{2}\sigma_{x_{2}x_{2}}\right)-\sigma_{p_{2}p_{2}}\right]$
$\displaystyle+$
$\displaystyle\omega_{1}\dot{\omega}_{1}\sigma_{x_{1}x_{1}}+\omega_{2}\dot{\omega}_{2}\sigma_{x_{2}x_{2}}$
$\displaystyle-$ $\displaystyle\dot{Q}_{c}^{(2)}-\dot{Q}_{h}^{(2)},$

where for simplicity we have removed the explicit time-dependence of various
quantities. In the slow driving limit ($\dot{\sigma}_{S}=0$), the above
expression becomes

$\displaystyle\dot{W}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{(2n_{c}+1)\dot{\omega}_{1}+(2n_{h}+1)\dot{\omega}_{2}}{2}+$
(69) $\displaystyle+$
$\displaystyle\frac{4\lambda^{2}(g_{c}^{2}+g_{h}^{2})(\dot{\omega}_{1}g_{c}^{2}-\dot{\omega}_{2}g_{h}^{2})(n_{1}-n_{2})}{\left(g_{c}^{2}+g_{h}^{2}\right)^{2}\left(4\lambda^{2}+g_{c}^{2}g_{h}^{2}\right)+4g_{c}^{2}g_{h}^{2}[\omega_{1}-\omega_{2}]^{2}}$
$\displaystyle-$ $\displaystyle\dot{Q}_{c}^{(2)}-\dot{Q}_{h}^{(2)}.$

The expression for the power under the assumption of
$\omega_{1}(t)=\omega_{2}(t)=\omega(t)$ simplifies to

$\displaystyle\dot{W}^{(2)}=\frac{\left\\{\frac{g_{c}^{2}g_{h}^{2}}{\lambda^{2}}\left(g_{c}^{2}+g_{h}^{2}\right)\left(1+n_{c}+n_{h}\right)+4\left[g_{c}^{2}(1+2n_{c})+g_{h}^{2}(1+2n_{h})\right]\right\\}\dot{\omega}-4g_{c}^{2}g_{h}^{2}\left(n_{c}-n_{h}\right)\left(\Omega_{c}-\Omega_{h}\right)}{\left(g_{c}^{2}+g_{h}^{2}\right)\left(4+\frac{g_{c}^{2}g_{h}^{2}}{\lambda^{2}}\right)}.$
(70)

##  Appendix B Entanglement of Gaussian states

In this section, we review the logarithmic negativity as an entanglement
quantifier for two quantum harmonic oscillators Vidal and Werner (2002);
Plenio (2005). We follow the method detailed, for instance, in Ref. Adesso and
Illuminati (2007), and redefine our position and momentum operators as

$\tilde{x}_{\alpha}=\sqrt{2}x_{\alpha},\;\;\;\;\;\tilde{p}_{\alpha}=\sqrt{2}p_{\alpha},$
(71)

where $\alpha=1,\;2$. We define the vector of the system’s quadratures

$\tilde{R}=\Big{\\{}\tilde{x}_{1},\;\tilde{p}_{1},\;\tilde{x}_{2},\;\tilde{p}_{2}\Big{\\}},$
(72)

the matrix

$\tilde{A}_{S}=\frac{\omega(t)}{4}\mathbb{I}_{4},$ (73)

proportional to the $4\times 4$ identity matrix $\mathbb{I}_{4}$ and the
system’s covariance matrix

$\tilde{\sigma}_{ij}=\frac{1}{2}\expectationvalue{\tilde{R}_{i}\tilde{R}_{j}+\tilde{R}_{j}\tilde{R}_{i}}-\langle\tilde{R}_{i}\rangle\langle\tilde{R}_{j}\rangle.$
(74)

Let us rewrite the system’s covariance matrix as

$\tilde{\sigma}_{S}(t)=\begin{pmatrix}\boldsymbol{A}&\boldsymbol{C}\\\
\boldsymbol{C}^{T}&\boldsymbol{B}\end{pmatrix},$ (75)

where

$\boldsymbol{A}=\begin{pmatrix}\tilde{\sigma}_{\tilde{x}_{1}\tilde{x}_{1}}(t)&\tilde{\sigma}_{\tilde{x}_{1}\tilde{p}_{1}}(t)\\\
\tilde{\sigma}_{\tilde{x}_{1}\tilde{p}_{1}}(t)&\tilde{\sigma}_{\tilde{p}_{1}\tilde{p}_{1}}(t)\end{pmatrix},$
(76)
$\boldsymbol{B}=\begin{pmatrix}\tilde{\sigma}_{\tilde{x}_{2}\tilde{x}_{2}}(t)&\tilde{\sigma}_{\tilde{x}_{2}\tilde{p}_{2}}(t)\\\
\tilde{\sigma}_{\tilde{x}_{2}\tilde{p}_{2}}(t)&\tilde{\sigma}_{\tilde{p}_{2}\tilde{p}_{2}}(t)\end{pmatrix},$
(77)
$\boldsymbol{C}=\begin{pmatrix}\tilde{\sigma}_{\tilde{x}_{1}\tilde{x}_{2}}(t)&\tilde{\sigma}_{\tilde{x}_{1}\tilde{p}_{2}}(t)\\\
\tilde{\sigma}_{\tilde{x}_{2}\tilde{p}_{1}}(t)&\tilde{\sigma}_{\tilde{p}_{1}{p}_{2}}(t)\end{pmatrix}.$
(78)

We now define

$\displaystyle I_{1}=\det[\boldsymbol{A}],\;\;\;\;\;$ $\displaystyle
I_{2}=\det[\boldsymbol{B}],$ $\displaystyle
I_{3}=\det[\boldsymbol{C}],\;\;\;\;\;$ $\displaystyle
I_{4}=\det[\tilde{\sigma}_{S}(t)].$ (79)

In order to compute the logarithmic negativity, we need to calculate the
symplectic eigenvalues of the partially-transposed covariance matrix which,
after some algebra, read

$\tilde{\nu}_{\pm}=\sqrt{\frac{1}{2}\left(\tilde{\Lambda}\pm\sqrt{\tilde{\Lambda}^{2}-4I_{4}}\right)},$
(80)

where

$\tilde{\Lambda}=I_{1}+I_{2}-2I_{3}.$ (81)

The logarithmic negativity is then defined in terms of the smallest symplectic
eigenvalue $\tilde{\nu}_{-}$ as

$\mathcal{E}_{N}=\max(0,-\ln\tilde{\nu}_{-}).$ (82)

The expression of the logarithmic negativity is generally complex but can be
simplified in special cases. For instance, applying the limiting cycle
solutions, without squeezing and with the additional conditions
$g_{c}=g_{h}=g$ and $\omega_{1}(t)=\omega_{2}(t)=\omega(t)$, we get

$\displaystyle\tilde{\nu}_{-}$ $\displaystyle=$
$\displaystyle\frac{1}{g^{4}+4\lambda^{2}}\Bigg{[}g^{8}\Big{(}1+2n_{c}(1+n_{c})+2n_{h}(1+n_{h})\Big{)}+4g^{4}\lambda^{2}\Big{(}2+n_{c}(4+n_{c})+n_{h}(4+n_{h})+6n_{c}n_{h})\Big{)}+16\lambda^{4}(1+n_{c}+n_{h})^{2}$
$\displaystyle-$ $\displaystyle
2\Bigg{[}g^{8}(n_{c}-n_{h})^{2}\Big{(}g^{8}(1+n_{c}+n_{h})^{2}+4g^{4}\lambda^{2}(2+n_{c}(4+n_{c})+n_{h}(4+n_{h})+6n_{c}n_{h})+16\lambda^{4}(1+n_{c}+n_{h})^{2}\Big{)}\Bigg{]}^{1/2}\Bigg{]}^{1/2}.$

For entanglement to be present, we must have $\tilde{\nu}_{-}<1$. We find that

$\lim_{\lambda\rightarrow\pm\infty}\tilde{\nu}_{-}=1+n_{c}+n_{h}\geq 1,$ (84)

as $n_{c},n_{h}\geq 0$. The values of $\lambda$ that give $\tilde{\nu}_{-}=1$
are

$\displaystyle\lambda$ $\displaystyle=$
$\displaystyle\pm\mathbbm{i}g^{2}\sqrt{\frac{n_{c}(1+n_{h})}{(n_{c}+n_{h})(2+n_{c}+n_{h})}},$
$\displaystyle\lambda$ $\displaystyle=$
$\displaystyle\pm\mathbbm{i}g^{2}\sqrt{\frac{n_{h}(1+n_{c})}{(n_{c}+n_{h})(2+n_{c}+n_{h})}}.$
(85)

Since the coupling constant $\lambda\in\mathbb{R}$, we cannot have
entanglement with these conditions as there is no $\lambda$ for which
$\tilde{\nu}_{-}<1$.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

## References

  * Kosloff and Levy (2014) R. Kosloff and A. Levy, Annu. Rev. Phys. Chem. 65, 365 (2014).
  * Deffner and Campbell (2019) S. Deffner and S. Campbell,  _Quantum Thermodynamics: An Introduction to the Thermodynamics of Quantum Information_ (Morgan & Claypool, 2019).
  * Mitchison (2019) M. T. Mitchison, Contemp. Phys. 60, 164 (2019).
  * Roßnagel  _et al._ (2016) J. Roßnagel, S. T. Dawkins, K. N. Tolazzi, O. Abah, E. Lutz, F. Schmidt-Kaler, and K. Singer, Science 352, 325 (2016).
  * Maslennikov  _et al._ (2019) G. Maslennikov, S. Ding, R. Hablützel, J. Gan, A. Roulet, S. Nimmrichter, J. Dai, V. Scarani, and D. Matsukevich, Nat. Commun. 10, 202 (2019).
  * Peterson  _et al._ (2019) J. P. S. Peterson, T. B. Batalhão, M. Herrera, A. M. Souza, R. S. Sarthour, I. S. Oliveira, and R. M. Serra, Phys. Rev. Lett. 123, 240601 (2019).
  * Gluza  _et al._ (2021) M. Gluza _et al._ , PRX Quantum 2, 030310 (2021).
  * Kosloff and Rezek (2017) R. Kosloff and Y. Rezek, Entropy 19, 136 (2017).
  * Singh and Müstecaplıoğlu (2020) V. Singh and O. E. Müstecaplıoğlu, Phys. Rev. E 102, 062123 (2020).
  * Reid  _et al._ (2017) B. Reid, S. Pigeon, M. Antezza, and G. De Chiara, Europhys. Lett. 120, 60006 (2017).
  * Deffner (2018) S. Deffner, Entropy 20, 875 (2018).
  * Serafini  _et al._ (2020) A. Serafini, M. Lostaglio, S. Longden, U. Shackerley-Bennett, C.-Y. Hsieh, and G. Adesso, Phys. Rev. Lett. 124, 010602 (2020).
  * Quan  _et al._ (2007) H. T. Quan, Y.-x. Liu, C. P. Sun, and F. Nori, Phys. Rev. E 76, 031105 (2007).
  * Hewgill  _et al._ (2018) A. Hewgill, A. Ferraro, and G. De Chiara, Phys. Rev. A 98, 042102 (2018).
  * Solfanelli  _et al._ (2020) A. Solfanelli, M. Falsetti, and M. Campisi, Phys. Rev. B 101, 054513 (2020).
  * Quan  _et al._ (2005) H. T. Quan, P. Zhang, and C. P. Sun, Phys. Rev. E 72, 056110 (2005).
  * Campisi and Fazio (2016) M. Campisi and R. Fazio, Nat. Commun. 7, 1 (2016).
  * Niedenzu and Kurizki (2018) W. Niedenzu and G. Kurizki, New J. Phys. 20, 113038 (2018).
  * Piccione  _et al._ (2021) N. Piccione, G. De Chiara, and B. Bellomo, Phys. Rev. A 103, 032211 (2021).
  * Uzdin  _et al._ (2015) R. Uzdin, A. Levy, and R. Kosloff, Phys. Rev. X 5, 031044 (2015).
  * Martinez and Paz (2013) E. A. Martinez and J. P. Paz, Phys. Rev. Lett. 110, 130406 (2013).
  * De Chiara  _et al._ (2018) G. De Chiara, G. Landi, A. Hewgill, B. Reid, A. Ferraro, A. J. Roncaglia, and M. Antezza, New J. Phys. 20, 113024 (2018).
  * Hewgill  _et al._ (2021) A. Hewgill, G. De Chiara, and A. Imparato, Phys. Rev. Res. 3, 013165 (2021).
  * Campbell and Vacchini (2021) S. Campbell and B. Vacchini, Europhys. Lett. 133, 60001 (2021).
  * Ciccarello  _et al._ (2021) F. Ciccarello, S. Lorenzo, V. Giovannetti, and G. M. Palma, (2021), arXiv:2106.11974 .
  * Albash  _et al._ (2012) T. Albash, S. Boixo, D. A. Lidar, and P. Zanardi, New J. Phys. 14, 123016 (2012).
  * Scali  _et al._ (2021) S. Scali, J. Anders, and L. A. Correa, Quantum 5, 451 (2021).
  * Braunstein and van Loock (2005) S. L. Braunstein and P. van Loock, Rev. Mod. Phys. 77, 513 (2005).
  * Nicacio  _et al._ (2016) F. Nicacio, M. Paternostro, and A. Ferraro, Phys. Rev. A 94, 052129 (2016).
  * Curzon and Ahlborn (1975) F. L. Curzon and B. Ahlborn, Am. J. Phys. 43, 22 (1975).
  * Asadian  _et al._ (2013) A. Asadian, D. Manzano, M. Tiersch, and H. J. Briegel, Phys. Rev. E 87, 012109 (2013).
  * Huang  _et al._ (2012) X. L. Huang, T. Wang, and X. X. Yi, Phys. Rev. E 86, 051105 (2012).
  * Correa  _et al._ (2015) L. A. Correa, J. P. Palao, D. Alonso, and G. Adesso, Sci. Rep. 4, 03949 (2015).
  * Alicki and Gelbwaser-Klimovsky (2015) R. Alicki and D. Gelbwaser-Klimovsky, New J. Phys. 17, 115012 (2015).
  * Latune  _et al._ (2019) C. Latune, I. Sinayskiy, and F. Petruccione, Sci. Rep. 9, 1 (2019).
  * Abah and Lutz (2014) O. Abah and E. Lutz, Europhys. Lett. 106, 20001 (2014).
  * Manzano  _et al._ (2016) G. Manzano, F. Galve, R. Zambrini, and J. M. R. Parrondo, Phys. Rev. E 93, 052120 (2016).
  * Klaers  _et al._ (2017) J. Klaers, S. Faelt, A. Imamoglu, and E. Togan, Phys. Rev. X 7, 031044 (2017).
  * Manzano (2018) G. Manzano, Phys. Rev. E 98, 042123 (2018).
  * Vidal and Werner (2002) G. Vidal and R. F. Werner, Phys. Rev. A 65, 032314 (2002).
  * Plenio (2005) M. B. Plenio, Phys. Rev. Lett. 95, 090503 (2005).
  * Adesso and Illuminati (2007) G. Adesso and F. Illuminati, J. Phys. A 40, 7821 (2007).