Now at Data Storage Institute, A*STAR

# Quantum and classical dynamics of a three-mode absorption refrigerator

Stefan Nimmrichter<EMAIL_ADDRESS>Centre for Quantum Technologies, National
University of Singapore, 3 Science Drive 2, Singapore 117543.  Jibo Dai
Centre for Quantum Technologies, National University of Singapore, 3 Science
Drive 2, Singapore 117543.  Alexandre Roulet  Centre for Quantum Technologies,
National University of Singapore, 3 Science Drive 2, Singapore 117543.  Now at
Department of Physics, University of Basel, Klingelbergstrasse 82, CH-4056
Basel, Switzerland  Valerio Scarani  Centre for Quantum Technologies, National
University of Singapore, 3 Science Drive 2, Singapore 117543.  Department of
Physics, National University of Singapore, 2 Science Drive 3, Singapore
117542.

###### Abstract

We study the quantum and classical evolution of a system of three harmonic
modes interacting via a trilinear Hamiltonian. With the modes prepared in
thermal states of different temperatures, this model describes the working
principle of an absorption refrigerator that transfers energy from a cold to a
hot environment at the expense of free energy provided by a high-temperature
work reservoir. Inspired by a recent experimental realization with trapped
ions, we elucidate key features of the coupling Hamiltonian that are relevant
for the refrigerator performance. The coherent system dynamics exhibits rapid
effective equilibration of the mode energies and correlations, as well as a
transient enhancement of the cooling performance at short times. We find that
these features can be fully reproduced in a classical framework.

##  1 Introduction

Historically, thermodynamics started with the goal of trying to understand how
to convert heat into useful mechanical motion. For that purpose, steam engines
have been developed which revolutionized our lives. Other useful thermal
machines, such as refrigerators and heat pumps, followed. Naturally, these
thermal machines are all large macroscopic entities that one uses classical
laws to describe them and study their performances.

Recently, as our ability to control small quantum systems progresses, it has
now become not only interesting, but also crucial to study what happens when
these machines become microscopic where quantum features become important [1,
2]. More specifically, there has been significant interest in determining the
extent to which quantum effects may help surpass classical limits such as the
Carnot efficiency [3, 4, 5, 6, 7, 8], and in how to export classical notions
such as work to the quantum domain [9, 10, 11, 12, 13]. Various quantum
thermal machines have been proposed and studied, with potential realizations
in optomechanical setups, superconducting circuits, atom-cavity systems and
trapped ions [14, 15, 16, 5, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28].

In the context of refrigerators, two quantum approaches stand out: the
smallest possible refrigerator consists of a qutrit [29] or three interacting
qubits [30], but a quantum absorption refrigerator can also be realised with
three interacting harmonic modes [31, 32]. This latter interaction, which has
been demonstrated with trapped ions [33], is the object of this study. We
focus on the unitary dynamics of the refrigerator itself and highlight
phenomena like effective equilibration, as well as the challenge of
identifying genuine quantum features. Elements on the process of
refrigeration, i.e. the transfer of energy from a cold to a hot bath, mediated
by the machine, are given in the Appendix.

Section 2 introduces the model Hamiltonian, the initial states that will be
considered (independent thermal states of different temperatures) and the
dynamical variables whose evolution is studied (occupation number of each
mode). We discuss why we focus on the unitary dynamics. This dynamics is then
solved and discussed in Section 3, first in general, then in the context of
refrigeration. The following two sections are devoted to the two main observed
features of the dynamics of the occupation numbers. In the short-term regime,
there is a cooling enhancement which would be absent if the interaction were
incoherent (Section 4). In the long-term regime, one observes effective
equilibration of the occupation numbers even if there is no dissipative
dynamics (Section 5). Finally, in Section 6 we study the classical model
obtained by replacing the mode operators with conjugate variables in the
Hamiltonian. We show that the dynamics of the occupation numbers exhibits the
same features observed in the quantum formalism.

##  2 Model

The model we study here is a system of three interacting harmonic oscillators,
with the free Hamiltonian

$\hat{H}_{0}=\sum_{i=h,w,c}\hbar\omega_{i}\left(\hat{a}_{i}^{\dagger}\hat{a}_{i}+\frac{1}{2}\right),$
(1)

and the interaction Hamiltonian

$\hat{H}_{1}=\hbar
g(\hat{a}_{h}^{\dagger}\hat{a}_{w}\hat{a}_{c}+\hat{a}_{h}\hat{a}_{w}^{\dagger}\hat{a}_{c}^{\dagger}).$
(2)

It is convenient to work in the interaction picture where

$\hat{H}_{\mathrm{int}}=\hbar
g(\hat{a}_{h}^{\dagger}\hat{a}_{w}\hat{a}_{c}\,e^{i\Delta
t}+\hat{a}_{h}\hat{a}_{w}^{\dagger}\hat{a}_{c}^{\dagger}\,e^{-i\Delta t}),$
(3)

with $\Delta=\omega_{h}-\omega_{w}-\omega_{c}$ and $g$ the coupling constant.
We will focus on the resonant case where $\Delta=0$ in what follows.

This Hamiltonian describes a wide range of physical processes: parametric
amplification, frequency conversion, and second harmonic generation [34, 35,
36]. This work is based on considering the three interacting modes as an
absorption refrigerator [31, 32, 33].

The refrigeration process works as follows. Each of the three modes is in
contact with a thermal bath: the cold bath at $T_{c}$, the hot bath at
$T_{h}>T_{c}$, and the work mode at $T_{w}$. We shall correspondingly refer to
$\hat{a}_{h}$, $\hat{a}_{w}$ and $\hat{a}_{c}$ as the hot mode, work mode, and
cold mode, respectively. The dynamical variables of interest are the
occupation numbers
$\bar{n}_{i}(t)=\mathrm{tr}\\{\rho(t)\hat{a}_{i}^{\dagger}\hat{a}_{i}\\}$ with
$i=h,w$ or $c$. Absorption refrigeration occurs if the interaction can induce
a stationary heat flow from the cold to the hot mode
($\dot{n}_{c}(t)<\dot{n}_{c}(0)$, $\dot{n}_{h}(t)>\dot{n}_{h}(0)$), the work
mode providing a sufficient amount of free energy.

A full description of the refrigeration process takes into account both the
unitary interaction generated by the trilinear Hamiltonian and the dissipative
dynamics due to the interaction of each mode with its bath. This is a
numerically involved task, further complicated by the subtle issue of
formulating a master equation for composite systems coupled to multiple baths
that is consistent with the second law of thermodynamics [37, 38, 39]. Yet, we
are interested in the regime of weak coupling to the baths such that all the
interesting aspects of the dynamics of refrigeration are captured by
restricting the analysis to the unitary dynamics. This is confirmed in the
Appendix, where we provide a comparison of the dynamics with and without the
dissipative coupling. In the rest of this paper, we shall thus focus on the
unitary dynamics, which has been experimentally implemented with trapped ions
[33].

The initial state of the system consists in the three modes being prepared in
uncorrelated thermal states,
$\rho(t=0)=\rho_{h}^{\mathrm{th}}\otimes\rho_{w}^{\mathrm{th}}\otimes\rho_{c}^{\mathrm{th}}$
with

$\rho_{i}^{\mathrm{th}}=\left[1-\exp\left(-\frac{\hbar\omega_{i}}{k_{B}T_{i}}\right)\right]\exp\left(-\frac{\hbar\omega_{i}}{k_{B}T_{i}}\hat{a}_{i}^{\dagger}\hat{a}_{i}\right).$
(4)

This corresponds to the situation in which each modes have been kept in
contact with their respective baths for long time, before turning on the
interaction. From here, focusing on the unitary dynamics implies that the
modes are effectively decoupled from their baths during the evolution; this
will capture the full refrigeration dynamics in the limit of slow
thermalization rate.

Before proceeding, let us notice that some states of this family are
stationary with respect to the unitary dynamics, and therefore also to the
complete dynamics too. Indeed $[\rho_{\text{st}},\hat{H}_{\mathrm{int}}]=0$ if

$\frac{1}{\bar{n}_{h}(0)}+1=\left(\frac{1}{\bar{n}_{w}(0)}+1\right)\left(\frac{1}{\bar{n}_{c}(0)}+1\right).$
(5)

##  3 Unitary dynamics

###  3.1 Methods and general features

We will mostly focus on temperatures $T_{i}$ that correspond to comparably
small initial average occupation numbers,
$\bar{n}_{i}(0)=\mathrm{tr}\\{\rho(t=0)\hat{a}_{i}^{\dagger}\hat{a}_{i}\\}$
with $i=h,w$ or $c$. For later comparison with the classical framework, we
will plot the initial average energy of each mode in the diagrams,
$\epsilon_{i}(0)=\hbar\omega_{i}(\bar{n}_{i}(0)+1/2)$, instead of the
occupation number.

At time $t=0$, the interaction Hamiltonian $\hat{H}_{\mathrm{int}}$ is
switched on and the system evolves unitarily according to
$\rho(t)=\hat{U}\rho_{0}\hat{U}^{\dagger}$ with
$\hat{U}=\mathrm{exp}(-i\hat{H}_{\mathrm{int}}t)$. This coherently transfers
populations between the three modes and changes the average energies
$\epsilon_{i}(t)$. However, we note that even the closed system is not
amenable to an exact analytical solution; previous studies focused on either
short-time behavior, resorted to using products of coherent or Fock states as
the initial state, or considered limiting cases of average occupation numbers
much larger than one [34, 40, 35, 36].

For an efficient simulation of the system dynamics, we take advantage of the
fact that the interaction Hamiltonian couples only Fock states of the form
$|n_{h},n_{w},n_{c}\rangle=|n,N-n,M-n\rangle$ with fixed integers $N$ and $M$
and $0\leq n\leq\mathrm{min}\\{N,M\\}$. That is, $\hat{H}_{\rm int}$ is block-
diagonal with respect to finite-dimensional subspaces characterized by the two
conserved quantities $N$ and $M$ and the dimension $d=\min\\{N,M\\}+1$. The
unitary evolution can then be efficiently computed by diagonalizing the
Hamiltonian in each of the contributing subspaces, up to a cutoff for both $N$
and $M$ that truncates the originally infinite dimensional Hilbert space. For
the simulations presented in this paper, the cutoff is chosen to ignore
populations in the initial density matrix smaller than $10^{-4}$.

Fig. 1 shows the typical time evolution of the modes’ energy under the unitary
dynamics for a thermal state of the form (4). In particular, we have chosen
initial average occupations of $\bar{n}_{h}(0)=0.5$, $\bar{n}_{w}(0)=2.5$, and
$\bar{n}_{c}(0)=2.0$, so that the order of magnitude of the contributing $N$
and $M$ are $10^{0}$ or $10^{1}$. It is readily seen that
$\bar{N}=\bar{n}_{h}(t)+\bar{n}_{w}(t)$, and
$\bar{M}=\bar{n}_{h}(t)+\bar{n}_{c}(t)$ are conserved.

Figure 1: Time evolution of the average energy (in units of $\hbar\omega_{i}$)
of the hot, work and cold mode. The initial state is a product thermal state
with average occupations $\bar{n}_{h}(0)=0.5$, $\bar{n}_{w}(0)=2.5$, and
$\bar{n}_{c}(0)=2.0$. Time is measured in units of the inverse coupling
strength $1/g$. The dashed lines show the values associated to the infinite
time-averaged state (6). The initial transient oscillations are magnified in
the inset.

The curves in Fig. 1 exhibit two noteworthy features. Firstly, the evolution
starts with a transient oscillatory stage, where the largest fluctuation away
from the initial value is found. We shall come back to this feature in Section
4. Secondly, the system energies seem to approach some apparent equilibrium
values (dashed lines), around which only small residual oscillations persist.
One can obtain those values from the infinite time-averaged state

$\sigma\equiv\langle\rho(t)\rangle_{\infty}=\lim_{\tau\rightarrow\infty}\frac{1}{\tau}\int_{0}^{\tau}\text{d}t\,\rho(t).$
(6)

The observation that the expectation value of certain dynamical variables
approach their long-time average, in spite of the fact that the system is not
converging to a steady state because the dynamics is unitary, is a phenomenon
known as effective equilibration [41]. We will get back to it in Section 5.

Figure 2: Change in the cold mode energy from the initial thermal average
$\epsilon_{c}(0)=2.5\hbar\omega_{c}$ to the infinite time-averaged value for
varied initial energies of the hot and the work mode. All energies are given
in units of the respective excitation quanta. The thick line indicates the
stationary configurations of the interaction Hamiltonian, where no change
occurs.

Conservation of $N$ and $M$ implies that any increase in the average energy of
the hot mode corresponds to an equal decrease in the work mode and cold mode
energies. Hence it is sufficient to focus on just one of the modes. In Fig. 2,
we plot the change of the average cold mode energy (or occupation number)
between the long-time equilibrium state $\sigma$ and the initial state for
different initial energies of the work and the hot mode. The initial
occupation number in mode $c$ is fixed at $\bar{n}_{c}(0)=2.0$. One sees that,
depending on the values of $\bar{n}_{h}(0)$ and $\bar{n}_{w}(0)$, the
equilibrium energy of mode $c$ can either increase (red), decrease (blue), or
stay unchanged (thick line).

Let us remark again that, by preparing the initial system in a product of
thermal states and then bringing them into interaction, we mimic the
thermodynamic situation where these three subsystems each have thermalized
with their own baths and are then allowed to exchange heat among each other.
In Appendix A, we explicitly compare this two-step treatment with a
simultaneous interaction-dissipation model based on a heuristic master
equation. The latter contains the interaction Hamiltonian and independent
thermal dissipators for each mode. When the dissipation rates are small
compared to the coupling rate $g$, we find that the time evolution of the
system energies is similar to that of the purely unitary dynamics.
Specifically, the initial transient features are approximately the same, and
the only relevant difference is a slight deviation of the steady-state
energies from the long-time averages predicted by the unitary model. Strong
dissipation rates, on the other hand, would only thwart the three-mode
interaction dynamics and freeze the system state close to the initial thermal
state.

###  3.2 Analogies and differences with refrigeration processes

After presenting the general features of the closed system dynamics governed
by the trilinear Hamiltonian (3), we now focus on the thermodynamic aspects.
As already mentioned, the interaction is capable of driving a quantum
absorption refrigerator, in close analogy to the earlier proposed three-qubit
fridge [30] with its interaction Hamiltonian
$\hat{H}_{\mathrm{int}}=g(|010\rangle\langle 101|+|101\rangle\langle 010|)$.
The cooling performance of the latter has been extensively studied, and we
will point out similarities and crucial differences for the present three-
oscillator case.

Figure 3: Time evolution of the average cold mode energy (in units of
$\hbar\omega_{c}$) for five different initial conditions ranging from a net
heating to a net cooling scenario (top to bottom). The initial hot and cold
mode energies are fixed at $\epsilon_{h}=1.5\hbar\omega_{h}$ and
$\epsilon_{c}=2.5\hbar\omega_{c}$, while the work mode energy starts from
$\epsilon_{w}=1.5\hbar\omega_{w}$ in the top panel and is increased in steps
of $\hbar\omega_{w}$ to the bottom. Note that this still admits a lower
temperature in the cold than in the hot mode, as $\omega_{c}<\omega_{h}$. The
middle panel corresponds to a stationary configuration.

As proposed in theory [31] and backed by observations with normal modes of
trapped ytterbium ions [33], the average energy of the $c$-mode can decrease
or increase, depending on the initial work mode value $\bar{n}_{w}(0)$. In
this sense, we can speak of cooling or heating, as demonstrated in Fig. 3 for
an exemplary choice of initial values, $\bar{n}_{h}(0)=1.0$,
$\bar{n}_{c}(0)=2.0$, and $\bar{n}_{w}(0)$ varying from $1.0$ to $5.0$ in unit
steps. In fact, the contour plot in Fig. 2 can be viewed as a phase diagram,
where the thick solid line separates the heating (red) from the cooling (blue)
regime. For the latter, the initial occupation numbers must satisfy

$\bar{n}_{w}(0)>\bar{n}_{h}(0)\frac{1+\bar{n}_{c}(0)}{\bar{n}_{c}(0)-\bar{n}_{h}(0)}.$
(7)

Since the initial states are thermal, the inequality translates into a
relation between the temperatures,
$T_{i}=\hbar\omega_{i}/k_{B}\text{ln}(1+1/\bar{n}_{i}(0))$. Notice here that
the refrigerator setting $T_{c}<T_{h}$ does not necessarily imply
$\bar{n}_{c}(0)<\bar{n}_{h}(0)$ due to the different eigenfrequencies.

The cooling regime (7) is consistent with the general _virtual qubit_
framework developed recently [42]. In this picture, one thinks of the hot mode
and work mode forming a set of virtual qubits by the levels
$|n_{h},n_{w}\rangle$ and $|n_{h}-1,n_{w}+1\rangle$ at an _effective virtual
temperature_

$T_{v}=\frac{\hbar\omega_{h}-\hbar\omega_{w}}{\hbar\omega_{h}/T_{h}-\hbar\omega_{w}/T_{w}}.$
(8)

The cold mode at initial temperature $T_{c}$ then interacts with the virtual
mode at $T_{v}$. In the regime $0\leq T_{v}<T_{c}$, cooling occurs since net
energy is transferred from the cold mode to the virtual mode at a lower
virtual temperature. At $T_{v}=T_{c}$, which corresponds to the equilibrium
condition (5), no energy exchange takes place between the two modes.

Figure 4: Thermal Fano factor (9) of the infinite time-averaged state reduced
to the cold mode for varied initial conditions as in Fig. 2. A value of zero
indicates a thermal energy distribution, i.e. that the energy variance is that
of a Gibbs state. We find that the variance is consistently above (below)
thermal in the cooling (heating) regime.

Two important remarks are now in order: Firstly and unsurprisingly, the
asymptotic three-mode state (6), which determines the long-time average
energies after interaction, is in general not a tensor product of thermal
states
$\sigma\neq\rho_{h}^{\mathrm{th}}\otimes\rho_{w}^{\mathrm{th}}\otimes\rho_{c}^{\mathrm{th}}$.
Secondly, the reduced single-mode states $\mathrm{tr}_{j\neq
i}\\{\sigma\\}=\rho_{i}^{\mathrm{th}}$ do not exhibit a thermal energy
distribution either. This can be easily illustrated by looking at the second
moments of the local distributions. The deviation of, say, the cold mode
energy variance in the time-averaged state $\sigma$ from that of a Gibbs
distribution can be captured in the Fano factor

$q=\frac{\left<(\hat{a}_{c}^{\dagger}\hat{a}_{c})^{2}\right>-\left(\left<\hat{a}_{c}^{\dagger}\hat{a}_{c}\right>\right)^{2}}{\left<\hat{a}_{c}^{\dagger}\hat{a}_{c}\right>\left(\left<\hat{a}_{c}^{\dagger}\hat{a}_{c}\right>+1\right)}-1.$
(9)

Fig. 4 shows a contour plot of this quantity as a function of initial average
energies as in Fig. 2. One sees that the variance is consistently higher
($q>0$) than that of a thermal state of the same energy in the cooling regime,
and lower in the heating regime. At the same time, we found that the single-
mode entropy is always lower than that of the thermal state, suggesting that
the time-averaged energy distribution is more biased.

Hence, and even though the reduced single-mode states remain diagonal at all
times, we cannot assign temperatures to the asymptotic state $\sigma$ that the
trilinear system approaches in the long-time average. This is in contrast to
the three-qubit scenario, where one can formally associate a temperature to
any diagonal single-qubit state [43, 44]. Here, the three modes are correlated
and explicitly driven out of thermal equilibrium by the trilinear interaction
(which remains the case in the presence of simultaneous weak bath couplings).
This highlights the difference between genuine thermalization in open systems
undergoing thermodynamic cycles and effective equilibration at an athermal
state in quantum systems with non-trivial Hamiltonians [45, 41, 46].

##  4 Enhanced cooling in the single-shot regime

One important signature of the cooling dynamics that can be seen from Fig. 3
is that the cold mode always overshoots to a lower-energy state at transient
time scales, before it approaches the long-time average value. Hence, if one
can control the interaction and halt the dynamics at the appropriate time, the
cooling performance can be enhanced. A similar, measurement-induced transient
cooling of a single qubit was found in Ref. [47]. In the three-qubit scenario,
this feature termed single-shot cooling was linked to the presence of quantum
coherence [43, 44]. Here, a similar effect can be observed as well [33]. The
difference is that, in the three-qubit case, the system passes through many
transient oscillations before reaching an equilibrium with help of the
simultaneous thermalization with three independent baths. Here, the closed
system alone approaches an effective equilibrium rather rapidly after the
first overshooting oscillation. Irregular oscillations around the asymptotic
state $\sigma$ then prevail, but at comparably small amplitudes. The higher
the initial temperatures, the more negligible these oscillations become.

Let us now examine the role of coherence in this transient effect. Following
the qubit analogy [44], we compare the coherent energy exchange driven by the
trilinear Hamiltonian to incoherent implementations of the same exchange
process. The simplest incoherent model is to assume the excitation of the hot
mode $\hat{L}=\hat{a}_{h}^{\dagger}\hat{a}_{w}\hat{a}_{c}$ and the de-
excitation $\hat{L}^{\dagger}$ occur spontaneously at the same rate $\gamma$.
The corresponding master equation in the interaction frame reads as

$\dot{\rho}=\gamma\left(\hat{L}\rho\hat{L}^{\dagger}+\hat{L}^{\dagger}\rho\hat{L}-\frac{1}{2}\left\\{\hat{L}^{\dagger}\hat{L}+\hat{L}\hat{L}^{\dagger},\rho\right\\}\right),$
(10)

which replaces the von Neumann equation with the trilinear Hamiltonian (2) of
the coherent model. One can easily check that the master equation still
conserves both $N$ and $M$, and that it does not build up coherence between
the combined Fock basis vectors $|n,N-n,M-n\rangle$; initially diagonal states
remain diagonal at all times. The latter facilitates an efficient numerical
implementation.

Alternatively, one can conceive a random unitary model where the coherent
energy exchange driven by $\hat{H}_{\rm int}\propto\hat{L}+\hat{L}^{\dagger}$
is switched on and off at random times with an average exchange rate $\gamma$,

$\dot{\rho}=\frac{\gamma}{2}\left[\hat{L}+\hat{L}^{\dagger},\left[\rho,\hat{L}+\hat{L}^{\dagger}\right]\right].$
(11)

The master equation describes dephasing in the eigenbasis of the interaction
Hamiltonian. The final equilibrium state is given by the fully dephased, or
infinite time-averaged, initial state (6). Note however that the final state
$\sigma$ does contain nondiagonal coherence terms when represented in terms of
the $|n,N-n,M-n\rangle$ basis. At the same time, neither the coherent nor the
two incoherent models generate coherence locally, i.e. nondiagonal elements in
the Fock state representation of the reduced single-mode states.

Figure 5: Comparison of coherent and incoherent time evolution for the average
cold mode energy, starting from initial thermal energies
$\epsilon_{h}(0)=1.5\hbar\omega_{h}$, $\epsilon_{w}(0)=5.5\hbar\omega_{w}$,
and $\epsilon_{c}(0)=2.5\hbar\omega_{c}$ (bottom panel in Fig. 3). The solid
line represents the coherent evolution, which exhibits a characteristic
transient overshoot to below the long-time average for $gt<1$. This behavior
is not predicted by the incoherent models (10) and (11) for the same energy
exchange between the modes, based on spontaneous excitation processes (dashed)
or on phase-incoherent exchange (dotted), respectively. The latter yields the
correct long-time average. We assume exchange rates equal to the coherent
coupling constant $g$.

An exemplary comparison of the models in terms of the time evolution of the
average cold mode energy is shown in Fig. 5, assuming $\gamma=g$ for
simplicity. As expected, neither the fully incoherent model (green dashed
line) nor the dephasing model (orange dotted) reproduce the enhanced cooling
feature of the coherent model (blue solid). Instead, both incoherent models
lead to an exponential decay towards their equilibrium values; the one of the
dephasing model (11) is the value that the coherent evolution oscillates
around and approaches in the long time average.

The time evolution and the equilibrium states differ slightly here, whereas in
the previously studied three-qubit scenario there is no difference between the
two incoherent models. This can be understood if we restrict to initial states
$\rho(0)$ from the two-dimensional subspace $N=M=1$. The trilinear interaction
is then equivalent to that of three qubits, and the two master equations (10)
and (11) are identical.

Finally, we note that the comparison between the coherent and incoherent
exchange process might suggest a quantum advantage in the transient cooling
performance over “classical” implementations. This argument is based on the
notion of “classical” states as a subset of quantum states without coherence
in the relevant basis representation. In Section 6, we provide an alternative
view by comparing the coherent quantum model to its counterpart in the fully
classical framework. It will turn out that the cooling enhancement is not a
genuine quantum effect. Indeed, the fact that the reduced single-mode quantum
states always remain diagonal in the single-mode Fock basis suggests that the
trilinear Hamiltonian only establishes coherence _between_ the modes. Such a
notion of inter-mode coherence, in the sense of a finite interference
capability and well-defined phase relation over a given coherence time, can
also be understood in classical optics [48].

##  5 Effective Equilibration

A striking feature of the studied three-mode interaction, which we have
observed and exploited in the previous sections, is the fast effective
equilibration of the average single-mode energies: even though the time
evolution is strictly reversible, they appear to converge towards the
stationary values of the time-averaged, fully dephased state, with only minor
oscillations around those values for times $t\gg 1/g$. The phenomenon of
effective equilibration has been explored in finite systems of high dimension
[49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 41, 61], and is observed here
for the case of initially thermal states. The effect is most pronounced at
high initial temperatures, when the thermal occupation of the modes extends to
high-dimensional subspaces associated to large values of the conserved
quantities $N$ and $M$. A proper equilibration in the presence of a bath would
lead to the same values, as long as the three-body coupling $g$ remains the
dominant rate (see Appendix A). This underpins the thermodynamical assessment
based on the time-averaged steady state of the three modes as a model
absorption refrigerator system in Sect. 3.2.

Figure 6: Features of effective equilibration: (a) Time $gt_{\rm ex}$ of the
transient overshoot extremum in the coherent evolution of the average mode
energies for varied initial conditions as in Fig. 2. The values were extracted
from a numerical simulation in steps of $0.01/g$. The black pixels correspond
to stationary states satisfying (5), where no such extremum is found. (b)
Purity of the reduced single-mode state as a function of time when the unitary
evolution starts from a pure triple Fock state $|100,N-100,M-100\rangle$ with
varying values for $N=M$. Notice the double-logarithmic scale. The dashed
lines represent the fully dephased states. Both the characteristic half-life
time and the long-time averages of the purity scale in proportion to the
inverse effective Hilbert subspace dimension $N+1$. (c) Exemplary energy gap
spectra associated to the trilinear Hamiltonian in the subspaces of $N=10^{4}$
and varying $M$. The histograms are binned in units of $2\hbar g$.

Quantitatively, one may speak of effective equilibration with respect to a set
of observables if the cumulative time averages of their expectation values
over the coherent time evolution quickly converge close to the values
associated to a dephased stationary state [51, 62, 53, 63]. The characteristic
convergence time scale and the residual difference should typically decrease
with the (finite or effective) system dimension [55, 57]. This implies that
recurrences of the initial values, if present, must be short-lived and rare
under unitary evolution.

We observe such behaviour in the present, formally infinite-dimensional,
scenario where the observables of interest are the single-mode energies. Our
numerical studies show that both the characteristic equilibration time and the
residual difference to the expectation values associated to the dephased state
(6) decrease with increasing initial occupation numbers $\bar{n}_{i}$. We
could not observe recurrences in our simulations, some of which extended over
time intervals a few orders of magnitude greater than $1/g$, see e.g. Appendix
A. In Fig. 6(a), we plot the time at which the cold mode energy assumes its
extremal overshoot value, which can serve as an estimate for the
characteristic equilibration time. The times are extracted from unitary time
evolutions with time increments of $0.01/g$ for varying initial conditions as
in Fig. 2. We observe that the values decrease with growing mode temperatures
and effective Hilbert space dimension covered by the initial state, except at
thermal equilibrium (black line).

We note that for pure initial Fock product states $|n,N-n,M-n\rangle$, one
rather observes persistent strong and irregular oscillations of the average
energies instead of the described equilibration. The latter is a consequence
of averaging over the incommensurate sub-spectra of the trilinear Hamiltonian
associated to its different orthogonal subspaces. Nevertheless, signatures of
effective equilibration also exist at the level of quantum states in high-
dimensional subspaces, e.g. by looking at entangling dynamics [64]. As an
example, we show in Fig. 6(b) how the reduced single-mode purity decays as a
function of time for an initial Fock product state with $n=100$ at varying
subspace dimension. The blue, red, and green lines (top to bottom) in the
double-logarithmic plot correspond to $N=M=300$, $3000$, and $30000$,
respectively. The dashed lines represent the reduced single-mode purities of
the respective dephased states. As expected, the trilinear interaction rapidly
entangles the three modes, and we observe a fast decay of single-mode purity
close to the dephased values. Both these values and the characteristic half-
life time decrease roughly in proportion to the inverse subspace dimension
$d=\min\\{N,M\\}+1$.

On the other hand, rapid equilibration in finite high-dimensional (e.g.
interacting many-body) systems has been linked to the non-Poissonian nature
and spread of the underlying energy gap spectrum [65, 66, 67, 41], and it is
assumed to appear in nonintegrable systems whose classical counterpart may
exhibit signatures of chaos [68, 69, 70].

The energy gap spectrum of the quantum model varies strongly with the choice
of the conserved quantities $(N,M)$, if restricted to a single invariant
subspace. At large subspace dimension, non-Poissonian energy gap distributions
are observed, in agreement with the observation of rapid equilibration. Three
exemplary histograms of energy gaps (i.e. differences between subsequent
energy levels) are shown in Fig. 6(c). They correspond to increasing values of
$M$ at fixed $N=10000$ and subspace dimension $d=N+1$, using a bin width of
$2\hbar g$. In all cases, the probability does not drop with the gap size in a
Poisson-like fashion, as associated to non-equilibrating Hamiltonians.

Taken together, these intricate equilibration features distinguish the
resonant three-mode Hamiltonian (3) from other simple exchange coupling
models. The resonant two-mode exchange coupling, for instance, can be
diagonalized explicitly and leads merely to a normal mode splitting with a
single fixed energy gap given by the coupling frequency. The same reversible
dynamics can be seen in the two-dimensional subspaces (for $N=1$ or $M=1$) of
the present system. At high initial excitations, the three-oscillator
interaction model presented here could provide a minimal, accessible, and
practically relevant testbed for studies of effective equilibration.

##  6 Classical analysis

The present quantum refrigerator model of three trilinearly coupled harmonic
oscillators admits direct benchmarking against the predictions of classical
physics. We arrive at the classical version of the system straightforwardly by
replacing the mode operators $\hat{a}_{i}$ and $\hat{a}^{\dagger}_{i}$ with
canonical complex variables $\alpha_{i}$ and $\alpha^{*}_{i}$. The Hamilton
function becomes

$H=\sum_{i=h,w,c}\hbar\omega_{i}|\alpha_{i}|^{2}+\hbar
g(\alpha_{h}\alpha^{*}_{w}\alpha^{*}_{c}+\alpha_{h}^{*}\alpha_{w}\alpha_{c}).$
(12)

The corresponding Hamilton equations of motion can be solved analytically
after switching to the action-angle representation [71]. For this, the complex
amplitudes are expressed in terms of magnitude and phase,
$\alpha_{i}=\sqrt{\iota_{i}}\mathrm{exp}(i\phi_{i})$, where the $\iota_{i}$
represents the energy contained in mode $i$ in units of $\hbar\omega_{i}$. As
in the quantum case, one finds that the sums $I_{1}=\iota_{h}+\iota_{w}$ and
$I_{2}=\iota_{h}+\iota_{c}$ are conserved, as well as the total phase
difference $\Phi=\phi_{w}+\phi_{c}-\phi_{h}$. Given the initial $\iota_{i}(0)$
and $\phi_{i}(0)$ and the additional constant of motion
$L=\iota_{h}\iota_{w}\iota_{c}\mathrm{cos}^{2}(\Phi)$, this leads to the
solution

$\displaystyle\iota_{h}(t)$
$\displaystyle=c+(b-c)\mathrm{sn}^{2}\left(\pm\sqrt{a-c}gt+\theta_{0}\,|\,m\right),$
(13) $\displaystyle\theta_{0}$
$\displaystyle=\text{sn}^{-1}\left(\sqrt{\frac{\iota_{h}(0)-c}{b-c}}\,\Big{|}\,m\right),\quad
m=\frac{b-c}{a-c}.$

Here, $\mathrm{sn}(u\,|\,m)$ denotes the Jacobi elliptic function111We adopt
the usual convention for the elliptic integrals [72], which differs from the
one in [71]. and $a>b>c$ are the three ordered solutions of the equations:
$a+b+c=I_{1}+I_{2}$, $ab+bc+ca=I_{1}I_{2}$, and $abc=L$. The correct sign in
the argument is given by that of $\dot{\iota}_{h}(0)\propto-\sin\Phi(0)$.

The solution (13) is periodic with period $2K(m)$, where $K$ denotes a
complete elliptic integral of the first kind. The time average of the energy
contained in each mode is obtained by integrating (13) over one period,
yielding a compact expression in terms of elliptic integrals of first and
second kind

$\langle\iota_{h}\rangle_{t}=c+\frac{b-c}{m}\left[1-\frac{E(m)}{K(m)}\right].$
(14)

This is the classical counterpart of the quantum long time average, which is
the mean energy associated to a completely dephased state.

For a comparison to the quantum simulation with thermal initial states, we
evaluate the classical time evolution in a Monte-Carlo simulation by drawing
initial angles $\phi_{i}(0)$ according to a uniform distribution and
$\iota_{i}(0)$ according to an exponential distribution with mean values
$\bar{\iota}_{i}=k_{B}T_{i}/\hbar\omega_{i}$. However, the different energy
statistics leaves room for a certain ambiguity in the matching of quantum and
classical initial conditions. One can either match the initial temperatures
$T_{i}$, arguing that the quantum and classical predictions shall be based on
the same physical boundary conditions. Then the initial mean energies
$\hbar\omega_{i}\bar{\iota}_{i}$ and
$\epsilon_{i}=\hbar\omega_{i}(\bar{n}_{i}+1/2)$ will differ slightly. Or one
matches $\bar{\iota}_{i}=\bar{n}_{i}+1/2$, which implies a better matching of
the predicted energy trajectories at the cost of slightly different reservoir
temperatures.

Figure 7: Comparison between the quantum and classical predictions for the
time evolution of the average cold mode energy, starting from the initial
conditions of Fig. 5. The classical results are obtained by a Monte-Carlo
simulation of the solution (13) associated to the Hamilton function (12),
averaging over random samples of initial conditions. The two curves each
correspond to $10^{7}$ trajectories drawn from a thermal distribution using
either the same average mode energies (I, green) or the same initial
temperatures (II, orange) as in the quantum case. The dotted lines indicate
the infinite time averages.

In Fig. 7, we compare both classical options to the quantum prediction for an
exemplary initial configuration in the cooling regime. The classical
simulations were averaged over $10^{7}$ trajectories each. The relative
difference to the quantum result decreases at higher initial temperatures
where the continuous classical distribution of energies approximates well the
discrete quantum statistics.

As expected, the classical result with matching initial energies (I) starts at
the same point and remains closer to quantum case than the result with
matching temperatures (II). Most importantly, the classical model predicts the
same dynamical features as the quantum model, namely the short-time cooling
enhancement and the effective equilibration at long times. The latter can
again be attributed to thermal averaging: Even though single classical
trajectories (13) describe periodic orbits, their periods depend on the
initial energies and relative phases between the modes, and so they vary
broadly over the thermal distribution of initial conditions. The averaged time
evolution is not periodic and does not exhibit recurrences even for very long
interaction times $gt\gg 1$.

In Section 4, we have seen that the transient enhancement known as single-shot
cooling is linked to the presence of coherence between eigenstates of the
interaction Hamiltonian in the quantum framework. Now the quantum-classical
comparison shows that this feature is not relying on genuinely quantum
coherence, but is rather a generic feature of the particular three-mode
exchange interaction. Indeed, it appears as well in the classical framework
where coherence between harmonic modes exists in the form of a fixed phase
relation between their oscillations. The transient cooling enhancement can be
suppressed by subjecting the classical phase coordinate $\Phi(t)$ to external
noise, just as the feature will disappear in the quantum case if dephasing is
present.

##  7 Conclusion

In this article, we have provided an in-depth study of a quantum absorption
refrigerator model based on three resonantly coupled harmonic oscillator
modes, as recently realized in an ion-trap experiment [33]. A closed-system
analysis already exhibits the key features of the earlier studied analoguous
three-qubit refrigerator model, in particular, the description of heating and
cooling regimes in terms of virtual temperatures, and the enhanced cooling
performance at transient time scales.

On the other hand, we have found important differences and new insights in the
more complex three-mode system. Even in the absence of simultaneous bath
coupling, the initially thermal energy distributions of the cold, the hot, and
the work modes equilibrate effectively around athermal distributions
corresponding to a state that is fully dephased in the eigenbasis of the
interaction Hamiltonian. This state exhibits residual coherences in the
combined three-mode Fock representation, and it differs from the steady state
of a fully incoherent implementation based on spontaneous excitation and de-
excitation.

Most strikingly, we found that all essential features of the model can be
explained by entirely classical means, with only minor quantitative
differences that become irrelevant with growing temperatures. This includes,
in particular, the enhanced transient cooling feature that had previously been
attributed to quantum coherence in the three-qubit scenario. With the
possibility of a one-to-one comparison between quantum and classical
predictions at hand, the three-mode system may in fact serve as an ideal
testbed to elucidate the role of quantum resources in thermodynamics.

We acknowledge clarifying discussions with Nicolas Brunner, Christian Gogolin,
Alioscia Hamma, Sandu Popescu, Anthony Short, and Paul Skrzypczyk. This
research is supported by the Ministry of Education, Singapore, and by the
National Research Foundation, Prime Ministers Office, Singapore, under the
Research Centres of Excellence Programme and through the Competitive Research
Programme (Award No. NRF-CRP12-2013-03).

##  Appendix A Influence of simultaneous bath coupling

Here we investigate the performance of the refrigerator when the three
harmonic modes are simultaneously interacting via the trilinear Hamiltonian
and each coupled to their individual thermal reservoirs. For that, we employ
the master equation

$\frac{\text{d}\rho(t)}{\text{d}t}=-i[\hat{H}_{\text{int}},\rho(t)]+\sum_{j=h,w,c}\mathcal{L}_{j}\rho(t),$
(A.1)

where $\mathcal{L}_{j}$ is the superoperator describing the effect of the bath
coupled with mode $j$. We explicitly assume that each mode couples to its own
bath locally and separately, i.e. we employ the usual Lindblad dissipators,

$\mathcal{L}_{j}\rho=\kappa_{j}(1+\bar{n}_{j})\mathcal{D}[\hat{a}_{j}]\rho+\kappa_{j}\bar{n}_{j}\mathcal{D}[\hat{a}_{j}^{\dagger}]\rho,$
(A.2)

with
$\mathcal{D}[\hat{L}]\rho=\hat{L}\rho\hat{L}^{\dagger}-\frac{1}{2}\\{\hat{L}^{\dagger}\hat{L},\rho\\}$.
Here $\kappa_{j}$ is the coupling strength of the mode $j$ with its bath and
$\bar{n}_{j}$ is the mean occupation number at the bath temperature $T_{j}$.
In the following, we assume that the bath couplings of the three modes are
equal, $\kappa_{h}=\kappa_{w}=\kappa_{c}=\kappa$.

It was pointed out that using local dissipators (A.2) in a system of inter-
coupled modes might not always yield consistent predictions for transient
quantities such as the energy flow between the subsystems [37, 38, 39].
However, since we are concerned here with the qualitative behavior of local
observables (namely, single-mode energies), local dissipators may be used as
an approximation.

In order to observe the refrigeration effect, the bath coupling $\kappa$
should not be too strong compared to the trilinear coupling $g$. Otherwise the
modes will stay thermalized with their baths and cannot evolve away from the
initial state. The system dynamics will be frozen.

Figure 8: Time evolution of the average cold mode energy for different
coupling rates $\kappa$ between the modes to their thermal baths. The open-
system dynamics is described by the master equation (A.1), which was simulated
numerically using the QuTiP package, whereas $\kappa=0$ denotes the closed
unitary evolution (solid line). We start from thermal initial states at
$\epsilon_{h}(0)=\hbar\omega_{h}$, $\epsilon_{w}(0)=3\hbar\omega_{w}$, and
$\epsilon_{c}(0)=2.5\hbar\omega_{c}$.

For $\kappa\ll g$ on the other hand, the dynamics is not expected to deviate
much from the unitary evolution studied in the main text. Indeed, this is
confirmed in Fig. 8, which compares the evolution of the average cold mode
energy for various thermalization rates to the unitary case. In particular,
the transient oscillations at short times are well reproduced for $\kappa\ll
g$ (green dashed and blue dotted line) and the steady-state value approached
in the long-time limit does not differ significantly from the infinite time
average of the unitary case (black solid line). We ran simulations up to
$gt=10^{4}$ to confirm that deviations as significant as the initial
oscillations do not recur.

Increasing the bath coupling $\kappa$ to the level of $g$ and beyond thwarts
the coherent energy exchange between the modes; the transient oscillations
disappear for $\kappa\geq g$ (red dash-dotted and purple fine-dotted line) and
the mode energy quickly approaches equilibrium close to the initial value.

Concluding, the unitary dynamics studied in the main text approximates the
transient refrigerator dynamics well in the relevant regime of operation,
$\kappa\ll g$, and it gives reasonably good account of the long-time cooling
performance.

## References

  * Horodecki and Oppenheim [2013] M. Horodecki and J. Oppenheim, Nat. Commun. 4, 2059 (2013).
  * Skrzypczyk  _et al._ [2014] P. Skrzypczyk, A. J. Short, and S. Popescu, Nat. Commun. 5, 4185 (2014).
  * Scully  _et al._ [2003] M. O. Scully, M. S. Zubairy, G. S. Agarwal, and H. Walther, Science 299, 862 (2003).
  * Dillenschneider and Lutz [2009] R. Dillenschneider and E. Lutz, Europhys. Lett. 88, 50003 (2009).
  * Roßnagel  _et al._ [2014] J. Roßnagel, O. Abah, F. Schmidt-Kaler, K. Singer, and E. Lutz, Phys. Rev. Lett. 112, 030602 (2014).
  * Correa  _et al._ [2014] L. A. Correa, J. P. Palao, D. Alonso, and G. Adesso, Sci. Rep. 4, 3949 (2014).
  * Strasberg  _et al._ [2017] P. Strasberg, G. Schaller, T. Brandes, and M. Esposito, Phys. Rev. X 7, 021003 (2017).
  * Niedenzu  _et al._ [2017] W. Niedenzu, V. Mukherjee, A. Ghosh, A. G. Kofman, and G. Kurizki,  (2017), arXiv:1703.02911 .
  * Allahverdyan  _et al._ [2004] A. E. Allahverdyan, R. Balian, and T. M. Nieuwenhuizen, Europhys. Lett. 67, 565 (2004).
  * Talkner  _et al._ [2007] P. Talkner, E. Lutz, and P. Hänggi, Phys. Rev. E 75, 050102 (2007).
  * Campisi  _et al._ [2011] M. Campisi, P. Hänggi, and P. Talkner, Rev. Mod. Phys. 83, 771 (2011).
  * Roncaglia  _et al._ [2014] A. J. Roncaglia, F. Cerisola, and J. P. Paz, Phys. Rev. Lett. 113, 250601 (2014).
  * Talkner and Hänggi [2016] P. Talkner and P. Hänggi, Phys. Rev. E 93, 022131 (2016).
  * Geva and Kosloff [1996] E. Geva and R. Kosloff, J. Chem. Phys. 104, 7681 (1996).
  * Venturelli  _et al._ [2013] D. Venturelli, R. Fazio, and V. Giovannetti, Phys. Rev. Lett. 110, 256801 (2013).
  * Kosloff and Levy [2014] R. Kosloff and A. Levy, Annu. Rev. Phys. Chem. 65, 365 (2014).
  * Mari  _et al._ [2015] A. Mari, A. Farace, and V. Giovannetti, J. Phys. B 48, 175501 (2015).
  * Hofer  _et al._ [2016a] P. P. Hofer, J.-R. Souquet, and A. A. Clerk, Phys. Rev. B 93, 041418 (2016a).
  * Joulain  _et al._ [2016] K. Joulain, J. Drevillon, Y. Ezzahri, and J. Ordonez-Miranda, Phys. Rev. Lett. 116, 200601 (2016).
  * Mitchison  _et al._ [2016] M. T. Mitchison, M. Huber, J. Prior, M. P. Woods, and M. B. Plenio, Quantum Sci. Technol. 1, 015001 (2016).
  * Hofer  _et al._ [2016b] P. P. Hofer, M. Perarnau-Llobet, J. B. Brask, R. Silva, M. Huber, and N. Brunner, Phys. Rev. B 94, 235420 (2016b).
  * Karimi and Pekola [2016] B. Karimi and J. P. Pekola, Phys. Rev. B 94, 184503 (2016).
  * Roulet  _et al._ [2017] A. Roulet, S. Nimmrichter, J. Arrazola, S. Seah, and V. Scarani, Phys. Rev. E 95, 062131 (2017).
  * Bissbort  _et al._ [2017] U. Bissbort, C. Teo, C. Guo, G. Casati, G. Benenti, and D. Poletti, Phys. Rev. E 95, 062143 (2017).
  * Zhang and Zhang [2017] K. Zhang and W. Zhang, Phys. Rev. A 95, 053870 (2017).
  * Reid  _et al._ [2017] B. Reid, S. Pigeon, M. Antezza, and G. De Chiara,  (2017), arXiv:1708.07435 .
  * Hardal  _et al._ [2017] A. Ü. C. Hardal, N. Aslan, C. M. Wilson, and Ö. E. Müstecaplıoğlu,  (2017), arXiv:1708.01182 .
  * Mu  _et al._ [2017] A. Mu, B. K. Agarwalla, G. Schaller, and D. Segal,  (2017), arXiv:1709.02835 .
  * Palao  _et al._ [2001] J. P. Palao, R. Kosloff, and J. M. Gordon, Phys. Rev. E 64, 056130 (2001).
  * Linden  _et al._ [2010] N. Linden, S. Popescu, and P. Skrzypczyk, Phys. Rev. Lett. 105, 130401 (2010).
  * Levy and Kosloff [2012] A. Levy and R. Kosloff, Phys. Rev. Lett. 108, 070604 (2012).
  * Goold  _et al._ [2016] J. Goold, M. Huber, A. Riera, L. del Rio, and P. Skrzypczyk, J. Phys. A Math. Theor. 49, 143001 (2016).
  * Maslennikov  _et al._ [2017] G. Maslennikov, S. Ding, R. Hablutzel, J. Gan, A. Roulet, S. Nimmrichter, J. Dai, V. Scarani, and D. Matsukevich,  (2017), arXiv:1702.08672 .
  * Walls and Barakat [1970] D. F. Walls and R. Barakat, Phys. Rev. A 1, 446 (1970).
  * [35] G. P. Agrawal and C. L. Mehta, J. Phys. A , 607.
  * Gambini [1977] R. Gambini, Phys. Rev. A 15, 1157 (1977).
  * Levy and Kosloff [2014] A. Levy and R. Kosloff, Europhys. Lett. 107, 20004 (2014).
  * González  _et al._ [2017] J. O. González, L. A. Correa, G. Nocerino, J. P. Palao, D. Alonso, and G. Adesso,  (2017), arXiv:1707.09228 .
  * Hofer  _et al._ [2017] P. P. Hofer, M. Perarnau-Llobet, L. D. M. Miranda, G. Haack, R. Silva, J. B. Brask, and N. Brunner,  (2017), arXiv:1707.09211 .
  * Bonifacio and Preparata [1970] R. Bonifacio and G. Preparata, Phys. Rev. A 2, 336 (1970).
  * Gogolin and Eisert [2016] C. Gogolin and J. Eisert, Rep. Prog. Phys. 79, 56001 (2016).
  * Brunner  _et al._ [2012] N. Brunner, N. Linden, S. Popescu, and P. Skrzypczyk, Phys. Rev. E 85, 051117 (2012).
  * Brask and Brunner [2015] J. B. Brask and N. Brunner, Phys. Rev. E 92, 062101 (2015).
  * Mitchison  _et al._ [2015] M. T. Mitchison, M. P. Woods, J. Prior, and M. Huber, New J. Phys. 17, 115013 (2015).
  * Kulkarni  _et al._ [2012] M. Kulkarni, K. L. Tiwari, and D. Segal, Phys. Rev. B 86, 155424 (2012).
  * Farrelly  _et al._ [2017] T. Farrelly, F. G. S. L. Brandão, and M. Cramer, Phys. Rev. Lett. 118, 140601 (2017).
  * Erez  _et al._ [2008] N. Erez, G. Gordon, M. Nest, and G. Kurizki, Nature 452, 724 (2008).
  * Glauber [2007] R. J. Glauber,  _Quantum Theory of Optical Coherence_ (Wiley-VCH, 2007).
  * Jensen and Shankar [1985] R. V. Jensen and R. Shankar, Phys. Rev. Lett. 54, 1879 (1985).
  * Tasaki [1998] H. Tasaki, Phys. Rev. Lett. 80, 1373 (1998).
  * Reimann [2008] P. Reimann, Phys. Rev. Lett. 101, 190403 (2008).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 (2008).
  * Short [2011] A. J. Short, New J. Phys. 13, 53009 (2011).
  * Ponomarev  _et al._ [2011] A. V. Ponomarev, S. Denisov, and P. Hänggi, Phys. Rev. Lett. 106, 010405 (2011).
  * Short and Farrelly [2012] A. J. Short and T. C. Farrelly, New J. Phys. 14, 13063 (2012).
  * Reimann [2012] P. Reimann, Phys. Scr. 86, 58512 (2012).
  * Malabarba  _et al._ [2014] A. S. L. Malabarba, L. P. García-Pintos, N. Linden, T. C. Farrelly, and A. J. Short, Phys. Rev. E 90, 012121 (2014).
  * Malabarba  _et al._ [2015] A. S. L. Malabarba, N. Linden, and A. J. Short, Phys. Rev. E 92, 062128 (2015).
  * Goldstein  _et al._ [2015] S. Goldstein, T. Hara, and H. Tasaki, New J. Phys. 17, 45002 (2015).
  * Kaufman  _et al._ [2016] A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Science 353, 794 (2016).
  * García-Pintos  _et al._ [2017] L. P. García-Pintos, N. Linden, A. S. L. Malabarba, A. J. Short, and A. Winter, Phys. Rev. X 7, 031027 (2017).
  * Linden  _et al._ [2009] N. Linden, S. Popescu, A. J. Short, and A. Winter, Phys. Rev. E 79, 061103 (2009).
  * Reimann and Kastner [2012] P. Reimann and M. Kastner, New J. Phys. 14, 43020 (2012).
  * Shaffer  _et al._ [2014] D. Shaffer, C. Chamon, A. Hamma, and E. R. Mucciolo, J. Stat. Mech. 2014, P12007 (2014).
  * Sengupta  _et al._ [2004] K. Sengupta, S. Powell, and S. Sachdev, Phys. Rev. A 69, 053616 (2004).
  * Manmana  _et al._ [2007] S. R. Manmana, S. Wessel, R. M. Noack, and A. Muramatsu, Phys. Rev. Lett. 98, 210405 (2007).
  * Luitz  _et al._ [2015] D. J. Luitz, N. Laflorencie, and F. Alet, Phys. Rev. B 91, 081103 (2015).
  * Kollath  _et al._ [2007] C. Kollath, A. M. Läuchli, and E. Altman, Phys. Rev. Lett. 98, 180601 (2007).
  * Rigol [2009] M. Rigol, Phys. Rev. Lett. 103, 100403 (2009).
  * Bañuls  _et al._ [2011] M. C. Bañuls, J. I. Cirac, and M. B. Hastings, Phys. Rev. Lett. 106, 050405 (2011).
  * Bajer and Miranowicz [2001] J. Bajer and A. Miranowicz, J. Opt. B Quantum Semiclass. Opt. 3, 251 (2001).
  * Abramowitz and Stegun [1965] M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables_ (Dover, New York, 1965).