# Toward relaxation asymmetry: Heating is faster than cooling

Tan Van Vu<EMAIL_ADDRESS>Current Address: Department of Physics,
Keio University, 3-14-1 Hiyoshi, Kohoku-ku, Yokohama 223-8522, Japan
Department of Information and Communication Engineering, Graduate School of
Information Science and Technology, The University of Tokyo, Tokyo 113-8656,
Japan  Yoshihiko Hasegawa<EMAIL_ADDRESS>Department of
Information and Communication Engineering, Graduate School of Information
Science and Technology, The University of Tokyo, Tokyo 113-8656, Japan

###### Abstract

An asymmetry in thermal relaxation toward equilibrium has been uncovered for
Langevin systems near stable minima [Phys. Rev. Lett. 125, 110602 (2020)]. It
has been shown that, given the same degree of nonequilibrium of the initial
distributions, relaxation from a lower temperature state (heating) is faster
than that from a higher temperature state (cooling). In this study, we
elucidate this relaxation asymmetry for discrete-state Markovian systems
described by the master equation. We rigorously prove that heating is faster
than cooling for arbitrary two-state systems, whereas for systems with more
than two distinct energy levels, the relaxation asymmetry is no longer
universal. Furthermore, for systems whose energy levels degenerate into two
energy states, we find that there exist critical thresholds of the energy gap.
Depending on the magnitude of the energy gap, heating can be faster or slower
than cooling, irrespective of the transition rates between states. Our results
clarify the relaxation asymmetry for discrete-state systems and reveal several
hidden features inherent in thermal relaxation.

##  I Introduction

Systems attached to thermal reservoirs will relax toward a stationary state.
Such thermal relaxation processes are ubiquitous in nature and possess rich
properties from both dynamic and thermodynamic perspectives. One of the
counterintuitive behaviors is the Mpemba effect Mpemba and Osborne (1969),
where cooling a hot system is faster than cooling a cold system. Such
nonmonotonic relaxation phenomena have been observed in various systems
Auerbach (1995); Lasanta  _et al._ (2017); Baity-Jesi  _et al._ (2019); Kumar
and Bechhoefer (2020); Biswas  _et al._ (2020) and theoretically analyzed for
microscopic dynamics Lu and Raz (2017); Klich  _et al._ (2019); Takada  _et
al._ (2021). In addition, it was found that cooling a system before heating it
could lead to exponentially fast relaxation Gal and Raz (2020). From the
perspective of thermodynamics, thermal relaxation processes exhibit universal
relations regarding irreversibility, which is quantified by irreversible
entropy production Landi and Paternostro (2021). Notably, it has been shown
that irreversible entropy production during thermal relaxation is lower-
bounded by information-theoretical Alhambra and Woods (2017); Shiraishi and
Saito (2019); Van Vu and Hasegawa (2021a) and geometrical Van Vu and Hasegawa
(2021b) distances between the initial and final states in both classical and
quantum regimes. These relations imply stronger inequalities than the
conventional second law of thermodynamics and impose geometrical constraints
on the possible relaxation path. Since thermal relaxation is important in
condensed matter Dattagupta (2012) and heat engines Benenti  _et al._ (2017),
deepening our understanding of thermal relaxation would benefit research in
these areas.

Consider preparing a thermal state corresponding to a given temperature via
thermal relaxation; i.e., attaching a system to a single reservoir and
allowing it to relax toward equilibrium. In this setting, the relaxation time
is a quantity of interest and can be approximately estimated via the
convergence rate of the system state toward the equilibrium state Temme  _et
al._ (2010). Given two identical systems initiated in thermal states, one at a
lower and the other at a higher temperature than the given temperature, a
natural question arises: which one relaxes faster? Recently, this question has
been addressed by Lapolla and Godec Lapolla and Godec (2020) for continuous-
state Langevin systems. By considering a pair of thermodynamically equidistant
temperature quenches (which have the same nonequilibrium free energy
difference), they unveiled an unforeseen asymmetry in thermal relaxation;
i.e., relaxation from a lower temperature is faster than that from a higher
temperature. Roughly speaking, it implies that heating up cold objects is
faster than cooling down hot objects. This phenomenon has been proven for
quenches of dynamics near stable minima; however, it is not universal for
generic systems because counterexamples have been constructed using multi-well
potentials Lapolla and Godec (2020). Another recent study Meibohm  _et al._
(2021) has reported that heating can occur faster or slower than cooling even
for anharmonic single-well potentials, and a crossover region emerges if the
quenches are not too far from equilibrium. Nonetheless, the relaxation
asymmetry may be universal for a specific class of multi-well potentials. It
is well known that overdamped diffusion under a multi-well potential with
sufficiently high barriers converges to Markov jump dynamics at long times. By
projecting the Langevin dynamics onto the Markov jump process between basins,
it has been shown that the general asymmetry is preserved in degenerate
potentials with separated time scales Lapolla and Godec (2020).

Since the average energy of a thermal state increases with the temperature,
the relaxation asymmetry allows us to say that, from the energetic
perspective, uphill relaxation is faster than downhill relaxation, which is
counterintuitive to an extent. Moreover, relaxation speed cannot be
characterized solely by thermodynamic quantities such as dissipation or
frenesy Maes (2020). Therefore, it is highly nontrivial that free energy plays
an essential role as a quantifier of nonequilibrium degree in equidistant
temperature quenches.

In this study, we elucidate the relaxation asymmetry for discrete-state
systems modeled by Markov jump processes, thus improving our understanding of
thermal relaxation. First, we prove that heating is faster than cooling in an
arbitrary two-state system, affirming that there is universality in the
relaxation asymmetry. However, we find that it is not the case for generic
systems with at least three distinct energy levels. By analytically
constructing counterexamples, we demonstrate that heating can be faster or
slower than cooling, depending on the transition rates. Nevertheless,
restricted to a particular class of systems, some universal results on the
relaxation asymmetry are obtained. We show that when the energy levels of the
system are two-state degenerate, there exist two critical energy gap
thresholds. Depending on whether the energy gap is larger or smaller than
these thresholds, it can be concluded with certainty that heating is faster or
slower than cooling. These theoretical results are numerically demonstrated
using several discrete-state systems.

##  II Setup

We consider the thermal relaxation process of an open system with $N$ states.
The system is coupled to a thermal reservoir at the inverse temperature
$\beta_{f}=(k_{\rm B}T_{f})^{-1}$, where $k_{\rm B}$ is the Boltzmann
constant. Owing to interaction with the thermal reservoir, stochastic
transitions between states are induced. The dynamics of the system is governed
by the master equation,

$\partial_{t}|p_{t}\rangle=\mathsf{R}|p_{t}\rangle,$ (1)

where $|p_{t}\rangle\coloneqq[p_{1}(t),\dots,p_{N}(t)]^{\top}$ denotes the
probability distribution of the system at time $t$; the matrix
$\mathsf{R}=[R_{mn}]\in\mathbb{R}^{N\times N}$ is time-independent with
$R_{mn}\geq 0$ denoting the transition rate from state $n$ to state $m\,(\neq
n)$, and $\sum_{m}R_{mn}=0$. Without loss of generality, we assume that
$E_{1}\geq\dots\geq E_{N}$, where $E_{n}$ denotes the energy of state $n$. The
transition rates satisfy the detailed balance condition,
$R_{mn}e^{-\beta_{f}E_{n}}=R_{nm}e^{-\beta_{f}E_{m}}$, which is a sufficient
condition such that the system always relaxes to the thermal Gibbs state
$|\pi^{f}\rangle$ after a sufficiently long time, irrespective of the initial
state. Here,

$\pi_{n}^{f}\coloneqq\frac{e^{-\beta_{f}E_{n}}}{Z_{\beta_{f}}}~{}\text{and}~{}Z_{\beta_{f}}\coloneqq\sum_{n=1}^{N}e^{-\beta_{f}E_{n}}.$
(2)

For $m\neq n$, the transition rate $R_{mn}$ can be expressed as

$R_{mn}=\Gamma e^{-\beta_{f}(B_{mn}-E_{n})},$ (3)

where $B_{mn}=B_{nm}$ are the barrier coefficients, and $\Gamma$ is a positive
constant.

Now, let us formulate the problem. We consider relaxation that initiates from
a thermal state $|\pi^{i}\rangle$ associated with the inverse temperature
$\beta_{i}=(k_{\rm B}T_{i})^{-1}$. This can be regarded as a temperature
quench $T_{i}\to T_{f}$ at time $t=0^{-}$. Given a pair of cold and hot
temperatures, $T_{c}$ and $T_{h}$, satisfying $T_{c}<T_{f}<T_{h}$, we
investigate the relaxation speed depending on the quench direction,
$T_{i}=T_{c}\uparrow T_{f}$ (heating) and $T_{i}=T_{h}\downarrow T_{f}$
(cooling). The degree of nonequilibrium (or free energy) of each initial state
is the same,

$D({\pi^{c}}\|{\pi^{f}})=D({\pi^{h}}\|{\pi^{f}}),$ (4)

where the relative entropy between two distributions $|p\rangle$ and
$|q\rangle$ is given by

$D({p}\|{q})\coloneqq\sum_{n}p_{n}\ln(p_{n}/q_{n}).$ (5)

Note that the characters, $c$ and $h$, are associated with the initial
temperatures; therefore, $c$ and $h$ correspond to the heating and cooling
processes, respectively. For convenience, we define $\mathcal{D}({p})\coloneqq
D({p}\|{\pi^{f}})$. We aim to answer the question of which quench direction
fastens the relaxation. To this end, we first need to quantify the relaxation
speed, which can be evaluated by the distance between the system state and the
thermal state. Analogous to Refs. Lu and Raz (2017); Lapolla and Godec (2020),
the relative entropy is used to measure the distance between states. In
thermal relaxation, the relative entropy is closely related to the free energy
and irreversible entropy production Seifert (2012) as

$\displaystyle\mathcal{D}(p_{t})$
$\displaystyle=\beta_{f}[\mathcal{F}(p_{t})-\mathcal{F}_{f}],$ (6)
$\displaystyle\Sigma_{t}$
$\displaystyle=\mathcal{D}(p_{0})-\mathcal{D}(p_{t}),$ (7)

where $\mathcal{F}(p_{t})=\sum_{n}p_{n}(t)[E_{n}+k_{\rm B}T_{f}\ln p_{n}(t)]$
is the free energy of the distribution $|p_{t}\rangle$,
$\mathcal{F}_{i}\coloneqq-k_{\rm B}T_{i}\ln Z_{\beta_{i}}$, and $\Sigma_{t}$
is the irreversible entropy production. Note that
$\mathcal{F}(p_{0^{+}})\neq\mathcal{F}_{i}$ since $\mathcal{F}(p_{0^{+}})$ and
$\mathcal{F}_{i}$ are evaluated with temperatures $T_{f}$ and $T_{i}$,
respectively. Let $|p^{i}_{t}\rangle$ be the time-evolution distribution
corresponding to the initial state $|\pi^{i}\rangle$, then heating is said to
be faster (slower) than cooling if

$\mathcal{D}({p^{c}_{t}})<(>)\mathcal{D}({p^{h}_{t}})$ (8)

in long times. Throughout this study, we compare the relaxation speeds of the
two quenches in the long-time regime, not over the entire evolution time.
Therefore, it is possible that the heating and cooling curves,
$\mathcal{D}(p_{t}^{c})$ and $\mathcal{D}(p_{t}^{h})$, may intersect at some
point. In what follows, we explain in detail how to determine the relaxation
speed in the long-time regime.

Let $0=\lambda_{1}>\lambda_{2}>\dots>\lambda_{N}$ be the eigenvalues of the
transition matrix and $\\{|v_{n}\rangle\\}_{n}$ be the set of corresponding
eigenvectors,

$\mathsf{R}|v_{n}\rangle=\lambda_{n}|v_{n}\rangle.$ (9)

Notably, all eigenvalues are real numbers since matrix $\mathsf{R}$ satisfies
the detailed balance condition Schnakenberg (1976). The eigenvectors
$\\{|v_{n}\rangle\\}$ form a basis for the space $\mathbb{R}^{N}$ with
$|v_{1}\rangle=|\pi^{f}\rangle$, and $|v_{n}\rangle$ is a traceless vector for
all $n\geq 2$ (i.e., $\langle 1|v_{n}\rangle=0$ since $\langle
1|\mathsf{R}=\langle 0|$). Here, $|0\rangle$ ($|1\rangle$) denotes the
$N$-dimensional vector with all zero (one) elements. Therefore, the initial
distribution $|\pi^{i}\rangle$ can be expressed as a linear combination of
$\\{|v_{n}\rangle\\}$ as follows:

$|\pi^{i}\rangle=|\pi^{f}\rangle+\sum_{n=2}^{N}\gamma_{n}^{i}|v_{n}\rangle,$
(10)

where $\gamma_{n}^{i}$’s are real numbers. Consequently, the probability
distribution at time $t$ can be analytically written in the following form:

$|p^{i}_{t}\rangle=|\pi^{f}\rangle+\sum_{n=2}^{N}\gamma_{n}^{i}e^{\lambda_{n}t}|v_{n}\rangle.$
(11)

In the long-time limit, the probability distribution $|p^{i}_{t}\rangle$ can
be approximated up to the second-order term as

$|p^{i}_{t}\rangle\simeq|\pi^{f}\rangle+\gamma_{2}^{i}e^{\lambda_{2}t}|v_{2}\rangle.$
(12)

Thus, the relaxation speed can be quantified via the value of
$|\gamma_{2}^{i}|$ Lu and Raz (2017). Accordingly, heating is faster (slower)
than cooling if $|\gamma_{2}^{c}|<(>)|\gamma_{2}^{h}|$ (see Appendix A for
proof).

A closed form of $\gamma_{2}^{i}$ can be obtained analytically Klich  _et al._
(2019). The transition rate matrix $\mathsf{R}$ can be transformed to a
symmetric matrix
$\overline{\mathsf{R}}=[\overline{R}_{mn}]\in\mathbb{R}^{N\times N}$ as
follows:

$\overline{\mathsf{R}}=\mathsf{F}^{1/2}\mathsf{R}\mathsf{F}^{-1/2},$ (13)

where $\mathsf{F}=[F_{mn}]\in\mathbb{R}^{N\times N}$ with
$F_{mn}=e^{\beta_{f}E_{n}}\delta_{mn}$. The elements of matrix
$\overline{\mathsf{R}}$ can be explicitly written in terms of the elements of
matrix $\mathsf{R}$ as

$\overline{R}_{mn}=e^{\beta_{f}(E_{m}-E_{n})/2}R_{mn}.$ (14)

Notably, matrix $\overline{\mathsf{R}}$ has the same eigenvalues as
$\mathsf{R}$, and its eigenvectors $\\{|f_{n}\rangle\\}$ are related to those
of $\mathsf{R}$ as $|f_{n}\rangle=\mathsf{F}^{1/2}|v_{n}\rangle$. Moreover,
these eigenvectors are mutually orthogonal, $\langle
f_{m}|f_{n}\rangle=\langle v_{n}|\mathsf{F}|v_{m}\rangle\delta_{mn}$.
Multiplying $\langle f_{2}|\mathsf{F}^{1/2}$ on both sides of Eq. (10), we can
show that $\gamma_{2}^{i}$ is proportional to the inner product between the
initial distribution $|\pi^{i}\rangle$ and the vector
$|\overline{f}_{2}\rangle$, given by

$\gamma_{2}^{i}=\frac{\langle f_{2}|\mathsf{F}^{1/2}|\pi^{i}\rangle}{\langle
f_{2}|f_{2}\rangle}=\frac{\langle\overline{f}_{2}|\pi^{i}\rangle}{\langle
f_{2}|f_{2}\rangle},$ (15)

where $|\overline{f}_{2}\rangle\coloneqq\mathsf{F}^{1/2}|f_{2}\rangle$. Note
that $\langle\overline{f}_{2}|\pi^{f}\rangle=0$ since $|f_{2}\rangle$ and
$|f_{1}\rangle=\mathsf{F}^{1/2}|\pi^{f}\rangle$ are orthogonal. Because the
sign of $\gamma_{n}^{i}$ can be absorbed by changing the eigenvectors
$|v_{n}\rangle\to-|v_{n}\rangle$, hereafter, we assume $\gamma_{2}^{h}\leq 0$,
which implies that $\langle\overline{f}_{2}|\pi^{h}\rangle\leq 0$.

##  III Results

Given the above setup, we now present our main results on the relaxation
asymmetry, including numerical illustrations and proofs.

###  III.1 Two-state systems

First, we consider two-state systems, for which universality regarding the
relaxation asymmetry can be achieved.

######  Result 1.

For two-state systems, heating is faster than cooling.

Result 1 affirmatively validates that the relaxation asymmetry is universal in
two-state systems. Even for two-state systems, it is highly nontrivial that
heating is faster than cooling. For discrete-state Markovian dynamics, the
speed of state transformation is constrained by time-antisymmetric dissipation
and time-symmetric frenesy (or dynamical activity) Shiraishi  _et al._ (2018);
Vo  _et al._ (2020). Relaxation from a thermal state at a higher temperature
has a higher dynamical activity; more precisely, the average number of jumps
over all the stochastic trajectories in the hot two-state system is always
greater than that in the cold two-state system. Consequently, one may
intuitively expect that cooling, which has higher dynamical activity, is
faster than heating. However, the result is counterintuitive, implying that
the dynamical activity alone cannot account for this relaxation asymmetry.

It is worth discussing the compatibility of Result 1 with the counterexamples
in Ref. Lapolla and Godec (2020), in which cooling can be faster than heating
for some double-well potentials. In the counterexamples therein, there is a
time-scale separation: a local equilibration arises prior to the terminal
exponential relaxation, which is effectively a two-state Markov jump process.
Although heating occurs faster than cooling in the early stages, the asymmetry
is inverted in the later stages. It has been observed that the breaking of the
asymmetry is intimately related to the configuration between the intra-well
and inter-well entropies Lapolla and Godec (2020). Therefore, the apparent
discrepancy between Result 1 and the terminal relaxation trend in the
counterexamples is due to the intra-well entropic contribution, which is
neglected in the Markov jump dynamics considered here.

To illustrate the above result, we use a two-state system [see Fig. 1(a)] with
the following transition matrix:

$\mathsf{R}=\begin{pmatrix}-1&e^{-\beta_{f}\varepsilon}\\\
1&-e^{-\beta_{f}\varepsilon}\end{pmatrix},$ (16)

where $\varepsilon=E_{1}-E_{2}>0$ is the energy gap. We vary the value of
$\varepsilon$ while fixing the inverse temperatures as follows: $\beta_{f}=1$,
$\beta_{h}=0.1$, and $\beta_{c}$ is uniquely determined via the condition in
Eq. (4). The time variation of ratio
$\mathcal{R}_{t}\coloneqq\mathcal{D}(p_{t}^{c})/\mathcal{D}(p_{t}^{h})$ is
plotted as a function of time $t$ in Fig. 1(b). Note that
$\mathcal{R}_{t}<(>)1$ implies that heating occurs faster (slower) than
cooling at time $t$. As shown, ratio $\mathcal{R}_{t}$ is always smaller than
$1$, and there is no intersection between the heating and cooling curves at
any finite time. Therefore, it is numerically verified that heating is always
faster than cooling. Although Result 1 only indicates the long-time behavior,
numerical evidence suggests that heating is faster than cooling for the entire
evolution time.

Figure 1: (a) Schematic of the two-state system. (b) Numerical illustration of
Result 1 in the two-state system. The energy gap, $\varepsilon$, is varied
from $0.2$ to $1.0$. The solid lines depict ratio $\mathcal{R}_{t}$ as a
function of time $t$ for each value of $\varepsilon$.

###### Proof of Result 1.

It suffices to prove that $|\gamma_{2}^{c}|<|\gamma_{2}^{h}|=-\gamma_{2}^{h}$.
For $N=2$, an arbitrary probability distribution $|\pi^{i}\rangle$ can be
expressed as a point $(x_{i},1-x_{i})$ in the two-dimensional space. Since
$E_{1}>E_{2}$, all thermal states, $|\pi^{c}\rangle,|\pi^{f}\rangle$, and
$|\pi^{h}\rangle$, lie on the segment with $(0,1)$ and $(1/2,1/2)$ as
endpoints. These probability distributions are geometrically illustrated in
Fig. 1. In the following, we employ a geometrical approach to prove Result 1.
From the conditions, $\langle\overline{f}_{2}|\pi^{f}\rangle=0$ and
$\langle\overline{f}_{2}|\pi^{h}\rangle\leq 0$, we can conclude that
$\langle\overline{f}_{2}|\pi^{c}\rangle\geq 0$ or $\gamma_{2}^{c}\geq 0$.
Thus, it is sufficient to show that

$\gamma_{2}^{c}+\gamma_{2}^{h}<0~{}\text{or}~{}\langle\overline{f}_{2}|(\pi^{c}+\pi^{h})/2\rangle<0,$
(17)

which is equivalent to proving that $(x_{c}+x_{h})/2>x_{f}$. We can rewrite
the equality $D({\pi^{c}}\|{\pi^{f}})=D({\pi^{h}}\|{\pi^{f}})$ as follows:

$S({\pi^{h}})-S({\pi^{c}})=(x_{h}-x_{c})\ln\frac{1-x_{f}}{x_{f}},$ (18)

where $S({p})\coloneqq-\sum_{n}p_{n}\ln p_{n}$ is the Shannon entropy of the
distribution $|p\rangle$. Note that $x_{c}<x_{h}<1/2$ and

$g(x_{i})\coloneqq dS({\pi^{i}})/dx_{i}=\ln[(1-x_{i})/x_{i}]$ (19)

is a strictly convex function over $x_{i}\in[0,1/2]$ (i.e., the second
derivative of $g(x_{i})$ with respect to $x_{i}$ is positive,
$g^{\prime\prime}(x_{i})>0$). Applying the Hermite–Hadamard inequality for
$g(x_{i})$, we obtain the following:

$\frac{1}{x_{h}-x_{c}}\int_{x_{c}}^{x_{h}}g(x_{i})dx_{i}>g\left(\frac{x_{c}+x_{h}}{2}\right),$
(20)

or equivalently,

$S({\pi^{h}})-S({\pi^{c}})>(x_{h}-x_{c})g\left(\frac{x_{c}+x_{h}}{2}\right).$
(21)

Combining Eqs. (18) and (21) results in the following inequality:

$g(x_{f})>g\left(\frac{x_{c}+x_{h}}{2}\right).$ (22)

Since $g(x_{i})$ is a strictly decreasing function, we have
$(x_{c}+x_{h})/2>x_{f}$, which completes the proof. ∎

Figure 2: Geometrical illustration of probability distributions. Two vectors
$|\overline{f}_{2}\rangle$ and $|\pi^{f}\rangle$ are orthogonal,
$\langle\overline{f}_{2}|\pi^{f}\rangle=0$.

###  III.2 Systems with at least three distinct energy levels

Next, we consider more general systems that have at least three distinct
energy levels; i.e., there exist three indices $1\leq i<j<k\leq N$ such that
$E_{i}>E_{j}>E_{k}$. For such systems, we obtain the following result.

######  Result 2.

For systems with more than two distinct energy levels, heating can be faster
or slower than cooling, depending on the transition rates.

Result 2 shows that no universality of the relaxation asymmetry is achieved in
the general case. With an appropriate choice of barrier coefficients, we can
construct a discrete-state system with $|\gamma_{2}^{c}|<|\gamma_{2}^{h}|$ or
$|\gamma_{2}^{c}|>|\gamma_{2}^{h}|$ as desired. This difference between the
continuous- and discrete-state systems can be explained as follows. For
simplicity, we consider a single-particle system described by the overdamped
Langevin equation. In this continuous-state system, the particle tends to
transit to a close place at any instant time. By contrast, in discrete-state
systems, the particle, in principle, can jump to anywhere, provided the
transition rate between these states is positive. This degree of freedom could
lead to a complicated relaxation compared with the continuous-state system. As
shown below, the construction of the transition matrix that determines the
magnitude relation between the heating and cooling rates is somewhat
artificial; therefore, there may not be a realistic system with such a
transition matrix. We anticipate that the universality of the relaxation
asymmetry will be achieved if appropriate constraints are placed on the
transition rates.

Figure 3: (a) Schematic of the three-state system. Numerical illustration of
Result 2 for the cases of (b)(c) faster heating and (d)(e) faster cooling.
(b)(d) The heating and cooling curves, $\mathcal{D}(p_{t}^{c})$ and
$\mathcal{D}(p_{t}^{h})$, are depicted by blue and red solid lines,
respectively. (c)(e) The ratio of
$\mathcal{R}_{t}=\mathcal{D}(p_{t}^{c})/\mathcal{D}(p_{t}^{h})$ for each case
is depicted by the green solid line.

We numerically illustrate Result 2 in a three-state system [see Fig. 3(a)].
The transition rates are determined using $\Gamma=1$ and

$e^{-\beta_{f}B_{mn}}=\frac{\gamma_{mn}}{2\cosh[\beta_{f}(E_{n}-E_{m})]}e^{-\beta_{f}(E_{m}+E_{n})/2},$
(23)

where $\gamma_{12}=\gamma_{21}=0.1$, $\gamma_{13}=\gamma_{31}=10$, and
$\gamma_{23}=\gamma_{32}=0$. The inverse temperatures are fixed as follows:
$\beta_{f}=1$, $\beta_{h}=0.5$, and $\beta_{c}$ is uniquely determined via the
condition in Eq. (4). The energy levels are $E_{1}=2$, $E_{2}=1.5$, and
$E_{3}=\pm 0.5$. The numerical results (i.e., the heating and cooling curves
and the $\mathcal{R}_{t}$ ratio) are plotted in Figs. 3(b) and 3(c) for the
$E_{3}=-0.5$ case and in Figs. 3(d) and 3(e) for the $E_{3}=0.5$ case. As
shown in Fig. 3(c), heating is faster than cooling in the short- and long-time
regimes. However, interestingly, two crossing points exist in the
intermediate-time regime, where cooling temporarily occurs faster than
heating. While it may be difficult to identify these two crossing points in
Fig. 3(b), we can clearly see in Fig. 3(c) that they appear around times
$t=10^{-1}$ and $t=20$. In contrast, in the $E_{3}=0.5$ case, Fig. 3(e) shows
that the $\mathcal{R}_{t}$ ratio is always greater than $1$, implying that
cooling is faster than heating for the entire evolution time. These numerical
results indicate that the relaxation asymmetry is not universal in three-state
systems, which is consistent with the theoretical finding of Result 2. Here,
for simplicity, we have performed numerical calculations with different energy
levels in the two cases. Nevertheless, as shown below, even under the
condition where the energy levels are fixed, it is easy to construct a
transition matrix that yields the desired relaxation trend.

###### Proof of Result 2.

We prove the result by analytically constructing a transition rate matrix such
that heating is slower than cooling. A transition rate matrix for the opposite
case can also be analogously constructed. First, one can prove that
$|\pi^{c}\rangle,|\pi^{h}\rangle$, and $|\pi^{f}\rangle$ are linearly
independent (see Appendix B). Let $\\{|e_{1}\rangle,|e_{2}\rangle\\}$ be an
orthogonal basis of the space $\mathcal{S}$ spanned by $|\pi^{f}\rangle$ and
$|\pi^{h}\rangle$, and $\mathsf{P}\coloneqq|e_{1}\rangle\langle
e_{1}|+|e_{2}\rangle\langle e_{2}|$ be the projection matrix to the space
$\mathcal{S}$. Define
$|\overline{f}_{2}\rangle\coloneqq|\pi^{c}\rangle-\mathsf{P}|\pi^{c}\rangle$,
then $|\overline{f}_{2}\rangle\neq|0\rangle$ because
$|\pi^{c}\rangle,|\pi^{h}\rangle$, and $|\pi^{f}\rangle$ are linearly
independent. Trivially,
$\langle\overline{f}_{2}|e_{1}\rangle=\langle\overline{f}_{2}|e_{2}\rangle=0$;
thus,
$\langle\overline{f}_{2}|\pi^{f}\rangle=\langle\overline{f}_{2}|\pi^{h}\rangle=0$.
Moreover,
$\langle\overline{f}_{2}|\pi^{c}\rangle=\langle\overline{f}_{2}|\overline{f}_{2}\rangle>0$
since $|\overline{f}_{2}\rangle\neq|0\rangle$.

Next, we construct a transition rate matrix $\mathsf{R}$ that results in
$|\gamma_{2}^{c}|>|\gamma_{2}^{h}|$. Set $B_{mn}=E_{m}+E_{n}$, one can
explicitly calculate that the matrix $\overline{\mathsf{R}}$ has a single zero
eigenvalue associated with the eigenvector
$|f_{1}\rangle=\mathsf{F}^{1/2}|\pi^{f}\rangle$, and the remaining eigenvalues
are all $-Z_{\beta_{f}}$ Klich  _et al._ (2019). Let
$\mathcal{U}=\\{|x\rangle\in\mathbb{R}^{N}~{}|~{}\langle x|f_{1}\rangle=0\\}$
be a subspace orthogonal to $|f_{1}\rangle$. Then, there exists an orthogonal
basis $\\{|f_{2}\rangle,|f_{3}\rangle,\dots,|f_{N}\rangle\\}$ of
$\mathcal{U}$, where
$|f_{2}\rangle\coloneqq\mathsf{F}^{-1/2}|\overline{f}_{2}\rangle$, since
$\langle f_{2}|f_{1}\rangle=0$ (i.e., $|f_{2}\rangle\in\mathcal{U}$).
Obviously, $|f_{n}\rangle~{}(n\geq 2)$ is an eigenvector of
$\overline{\mathsf{R}}$ with the corresponding eigenvalue $-Z_{\beta_{f}}$.
Following the idea in Ref. Klich  _et al._ (2019), we slightly modify
$\overline{\mathsf{R}}$ as follows:

$\overline{\mathsf{R}}\to\overline{\mathsf{R}}+\sum_{n=2}^{N}\epsilon_{n}\frac{|f_{n}\rangle\langle
f_{n}|}{\langle f_{n}|f_{n}\rangle}.$ (24)

Here, $Z_{\beta_{f}}>\epsilon_{2}>\dots>\epsilon_{N}\geq 0$ are small numbers
that ensure the positivity of $\overline{R}_{mn}~{}(m\neq n)$. It is easy to
check that $\overline{\mathsf{R}}|f_{1}\rangle=|0\rangle$ and
$\overline{\mathsf{R}}|f_{n}\rangle=(-Z_{\beta_{f}}+\epsilon_{n})|f_{n}\rangle$
for all $n\geq 2$. Now, the matrix $\overline{\mathsf{R}}$ has $N$ different
eigenvalues, and $|f_{2}\rangle$ is precisely the eigenvector corresponding to
the second-largest eigenvalue $\lambda_{2}=-Z_{\beta_{f}}+\epsilon_{2}$. The
transition rate matrix can be recovered as
$\mathsf{R}=\mathsf{F}^{-1/2}\overline{\mathsf{R}}\mathsf{F}^{1/2}$, and the
detailed balance condition is satisfied due to the symmetry of
$\overline{\mathsf{R}}$. With this construction, the relation between
$|\gamma_{2}^{c}|$ and $|\gamma_{2}^{h}|$ can be clarified as

$|\gamma_{2}^{c}|=\frac{\langle\overline{f}_{2}|\pi^{c}\rangle}{\langle
f_{2}|f_{2}\rangle}>0=\frac{\langle\overline{f}_{2}|\pi^{h}\rangle}{\langle
f_{2}|f_{2}\rangle}=|\gamma_{2}^{h}|,$ (25)

which completes the proof. ∎

###  III.3 Degenerate two-level systems

Last, we consider the remaining case, wherein the energy levels are degenerate
to two energy states. In other words, there exists an index $1\leq n<N$ such
that $E_{1}=\dots=E_{n}>E_{n+1}=\dots=E_{N}$. Such degenerated two-level
systems are seen in atoms Margalit  _et al._ (2013) and have been used to
enhance quantum-annealing performance Watabe  _et al._ (2020) and dissipation-
less heat current Tajima and Funo (2021). For convenience, we define the
energy gap $\Delta E\coloneqq E_{n}-E_{n+1}>0$. Remarkably, we find that,
depending on the magnitude of this energy gap, heating can be faster or slower
than cooling, regardless of the transition rates. Details are summarized in
the following.

######  Result 3.

If $\beta_{h}\Delta E\geq\ln[n/(N-n)]$, then heating is faster than cooling.
Conversely, if $\beta_{c}\Delta E\leq\ln[n/(N-n)]$, then heating is slower
than cooling.

Note that Result 1 can be derived from Result 3 by setting $N=2$ and $n=1$.
Result 3 indicates that there are two critical thresholds of the energy gap
$\Delta E$. As the energy gap is above or below these thresholds, a universal
conclusion on asymmetry in thermal relaxation can be drawn. It is also highly
nontrivial that the energy gap affects the relaxation speeds of heating and
cooling in this way. When $\Delta E$ is large, the jump from energy state
$E_{n+1}$ to $E_{n}$ is less likely to occur compared with the opposite jump.
Thus, heating is expected to be slower than cooling. However,
counterintuitively, heating is faster than cooling as $\Delta E$ is
sufficiently large. In addition, provided $n\leq N/2$, heating is faster than
cooling, regardless of the value of $\Delta E$. This implies that the number
of excited states also plays a crucial role in determining the relaxation
speed.

Again, we numerically illustrate Result 3 in a three-state system with $n=2$
[see Fig. 4(a)]. The transition rates are analogously defined as in Eq. (23)
with $\gamma_{13}=\gamma_{31}=10$, $\gamma_{23}=\gamma_{32}=0.1$, and
$\gamma_{12}=\gamma_{21}=0$. We consider two cases of the parameter setting:
(i) $\Delta E=(1+\kappa)\beta_{h}^{-1}\ln[n/(N-n)]$,
$(\beta_{f},\beta_{h})=(1,0.5)$ and (ii) $\Delta
E=\kappa\beta_{c}^{-1}\ln[n/(N-n)]$, $(\beta_{f},\beta_{c})=(1,1.5)$. Here,
$\kappa\in(0,1)$ is a tuning parameter. Note that once the above parameters
are provided, the other parameters are uniquely determined. According to
Result 3, cases (i) and (ii) correspond, respectively, to faster heating and
faster cooling. We vary the value of $\kappa$ from $0.1$ to $0.5$ and plot the
numerical results of cases (i) and (ii) in Figs. 4(b) and 4(c), respectively.
Figure 4(b) shows that $\mathcal{R}_{t}<1$ in the long-time regime, implying
that heating is faster than cooling. Interestingly, there are crossing points
between the heating and cooling curves in the intermediate-time regime. In
both the system considered here and the two-state system considered in the
previous section, the energy levels are two-state degenerate. However, the
former shows a more complicated relaxation trend than the latter. For case
(ii), conversely, $\mathcal{R}_{t}$ is always greater than $1$, implying that
cooling always occurs faster than heating. Consequently, all these numerical
results are consistent with Result 3.

Figure 4: (a) Schematic of the three-state system. Numerical illustrations of
Result 3 for (b) case (i) and (c) case (ii). The solid lines depict the ratio,
$\mathcal{R}_{t}=\mathcal{D}(p_{t}^{c})/\mathcal{D}(p_{t}^{h})$, for each
value of $\kappa$. The value of parameter $\kappa$ is varied from $0.1$ to
$0.5$.

###### Proof of Result 3.

We employ the same strategy used in proving Result 1. We prove the former case
first, i.e., $\beta_{h}\Delta E\geq\ln[n/(N-n)]$ leads to a faster heating. It
can be observed that all points lying on the segment $\ell$ with
$|\pi^{c}\rangle$ and $|\pi^{h}\rangle$ as endpoints are thermal states. It is
also evident that $|\pi^{c}\rangle,|\pi^{f}\rangle$, and $|\pi^{h}\rangle$ are
linearly dependent, i.e., there exists a real number $a\in(0,1)$ such that
$|\pi^{f}\rangle=a|\pi^{c}\rangle+(1-a)|\pi^{h}\rangle$. Since
$\langle\overline{f}_{2}|\pi^{f}\rangle=0$ and
$\langle\overline{f}_{2}|\pi^{h}\rangle\leq 0$,
$\langle\overline{f}_{2}|\pi^{c}\rangle\geq 0$ or $\gamma_{2}^{c}\geq 0$
follows from the continuity of the inner product
$\langle\overline{f}_{2}|p\rangle$ for $|p\rangle\in\ell$. Thus, it suffices
to show that $\gamma_{2}^{c}+\gamma_{2}^{h}<0$ or
$\langle\overline{f}_{2}|(\pi^{c}+\pi^{h})/2\rangle<0$, which is equivalent to
proving that $(x_{c}+x_{h})/2>x_{f}$, where $x_{i}=\pi_{n}^{i}$ completely
characterizes the thermal state $|\pi^{i}\rangle$. The condition
$D({\pi^{c}}\|{\pi^{f}})=D({\pi^{h}}\|{\pi^{f}})$ can be rewritten as follows:

$S({\pi^{h}})-S({\pi^{c}})=(x_{h}-x_{c})n\ln\frac{1-nx_{f}}{(N-n)x_{f}}.$ (26)

Note that $x_{c}<x_{h}\leq 1/(2n)$ since $\beta_{h}\Delta E\geq\ln[n/(N-n)]$,
and $g(x_{i})\coloneqq dS({\pi^{i}})/dx_{i}=n\ln[(1-nx_{i})/(N-n)x_{i}]$ is a
strictly convex function over $x_{i}\in[0,1/(2n)]$. Applying the
Hermite–Hadamard inequality for $g(x_{i})$ and following the same steps as in
Eqs. (20) and (21), we obtain $(x_{c}+x_{h})/2>x_{f}$, which proves the former
case.

When $\beta_{c}\Delta E\leq\ln[n/(N-n)]$, one can derive that $1/(2n)\leq
x_{c}<x_{h}$, and $g(x_{i})$ is a strictly concave function over
$x_{i}\in[1/(2n),1]$ (i.e., the second derivative of $g(x_{i})$ with respect
to $x_{i}$ is negative, $g^{\prime\prime}(x_{i})<0$). Applying the
Hermite–Hadamard inequality for $g(x_{i})$, we obtain the following:

$g(x_{f})=\frac{1}{x_{h}-x_{c}}\int_{x_{c}}^{x_{h}}g(x_{i})dx_{i}<g\left(\frac{x_{c}+x_{h}}{2}\right),$
(27)

or $x_{f}>(x_{c}+x_{h})/2$, which implies $|\gamma_{2}^{h}|<|\gamma_{2}^{c}|$.

∎

##  IV Summary and Discussion

In this study, we elucidated the relaxation asymmetry for discrete-state
systems described by Markov jump processes. We proved that the relaxation
asymmetry is universal in two-state systems, but not in generic systems with
more than two distinct energy levels. For systems with two degenerate energy
levels, we obtained some universal results indicating that the asymmetry in
thermal relaxation depends on the energy gap and the number of excited states.

Notably, the relaxation asymmetry has recently been numerically studied for
few-level open quantum systems described by the Lindblad master equations
Manikandan (2021). When the initial density matrix contains no coherence, the
quantum systems can be described by classical Markov jump processes with the
population distributions. It has been shown that heating is always faster than
cooling for two-level systems, whereas heating can be faster or slower than
cooling for three- and four-level systems. These numerical demonstrations are
consistent with Results 1 and 2, thus affirmatively supporting our theoretical
findings.

Although the relaxation asymmetry is not universal in generic discrete-state
systems, universality may be achieved by imposing some constraints on the
transition rates. Such a question requires further investigation and will be
addressed in future work. It would also be interesting to study the relaxation
asymmetry for non-Markovian systems Yang and Hou (2020). For instance, the
relaxation of some non-Markovian processes, such as a tagged particle in a
single file and the end-to-end distance in a Rouse polymer, has been studied
in Ref. Lapolla and Godec (2020). Another possible direction involves
formulating and investigating the relaxation asymmetry in open quantum systems
Breuer and Petruccione (2002); Nava and Fabrizio (2019) where coherence is
present in the initial state.

###### Acknowledgements.

We thank Keiji Saito for insightful discussions. This work was supported by
Ministry of Education, Culture, Sports, Science and Technology (MEXT) KAKENHI
Grant No. JP19K12153.

##  Appendix A Quantification of relaxation speed

Given
$|\pi^{i}\rangle=|\pi^{f}\rangle+\sum_{n=2}^{N}\gamma_{n}^{i}|v_{n}\rangle$
and $|\gamma_{2}^{c}|<|\gamma_{2}^{h}|$, we will show that
$\mathcal{D}({p^{c}_{t}})<\mathcal{D}({p^{h}_{t}})$ as $t\to\infty$. Note that
$|p^{i}_{t}\rangle=|\pi^{f}\rangle+\sum_{n=2}^{N}\gamma_{n}^{i}e^{\lambda_{n}t}|v_{n}\rangle$.
In the long-time limit, the term
$\sum_{n=2}^{N}\gamma_{n}^{i}e^{\lambda_{n}t}|v_{n}\rangle$ vanishes. Since

$D({p}\|{p+dp})=\sum_{k=1}^{N}\frac{dp_{k}^{2}}{p_{k}}+O(\Delta^{3}),$ (28)

where $\Delta=\sum_{k=1}^{N}|dp_{k}|$, we can approximate

$\mathcal{D}({p^{i}_{t}})=\sum_{k=1}^{N}\frac{\left(\sum_{n=2}^{N}\gamma_{n}^{i}e^{\lambda_{n}t}v_{nk}\right)^{2}}{\pi_{k}^{f}}+O(e^{3\lambda_{2}t}).$
(29)

Here, $|v_{n}\rangle=[v_{n1},\dots,v_{nN}]^{\top}$. Consequently, we have

$\displaystyle\mathcal{D}({p^{h}_{t}})-\mathcal{D}({p^{c}_{t}})$
$\displaystyle=\left[(\gamma_{2}^{h})^{2}-(\gamma_{2}^{c})^{2}\right]\sum_{k=1}^{N}\frac{v_{2k}^{2}}{\pi_{k}^{f}}e^{2\lambda_{2}t}$
(30) $\displaystyle+\sum_{\begin{subarray}{c}m,n\geq 2\\\
\max\\{m,n\\}>2\end{subarray}}a_{mn}e^{(\lambda_{m}+\lambda_{n})t}+O(e^{3\lambda_{2}t}),$

where $a_{mn}$ are constants. The first term on the right-hand side is
positive since $|\gamma_{2}^{c}|<|\gamma_{2}^{h}|$. The remaining terms may be
negative; however, they are negligible compared with the first term in the
long-time limit. Therefore,
$\mathcal{D}({p^{c}_{t}})<\mathcal{D}({p^{h}_{t}})$ as $t\to\infty$.

##  Appendix B Proof of linear independence

To prove the linear independence of $|\pi^{c}\rangle,|\pi^{f}\rangle$, and
$|\pi^{h}\rangle$, it is sufficient to show that the determinant of the
following matrix is negative.

$\mathsf{X}=\begin{pmatrix}e^{-\beta_{c}E_{1}}&e^{-\beta_{c}E_{2}}&e^{-\beta_{c}E_{3}}\\\
e^{-\beta_{f}E_{1}}&e^{-\beta_{f}E_{2}}&e^{-\beta_{f}E_{3}}\\\
e^{-\beta_{h}E_{1}}&e^{-\beta_{h}E_{2}}&e^{-\beta_{h}E_{3}}\end{pmatrix}.$
(31)

Here, $E_{1}>E_{2}>E_{3}$ and $\beta_{c}>\beta_{f}>\beta_{h}$. Without loss of
generality, we assume that $E_{3}=0$ and $\beta_{h}=0$. In this case, the
determinant can be calculated as follows:

$|\mathsf{X}|=(1-e^{-\beta_{c}E_{1}})(1-e^{-\beta_{f}E_{2}})-(1-e^{-\beta_{c}E_{2}})(1-e^{-\beta_{f}E_{1}}).$
(32)

Therefore, $|\mathsf{X}|<0$ is equivalent to

$\frac{1-e^{-\beta_{c}E_{1}}}{1-e^{-\beta_{f}E_{1}}}<\frac{1-e^{-\beta_{c}E_{2}}}{1-e^{-\beta_{f}E_{2}}}.$
(33)

Set $h(x)\coloneqq(1-e^{-\beta_{c}x})/(1-e^{-\beta_{f}x})$, Eq. (33) is
equivalent to $h(E_{1})<h(E_{2})$. We need only prove that $h(x)$ is a
strictly decreasing function over $x>0$. Taking the derivative of $h(x)$ with
respect to $x$, we have

$\frac{dh(x)}{dx}=\frac{e^{(\beta_{f}-\beta_{c})x}[\beta_{c}(e^{\beta_{f}x}-1)-\beta_{f}(e^{\beta_{c}x}-1)]}{(e^{\beta_{f}x}-1)^{2}}.$
(34)

Since

$\displaystyle\beta_{c}(e^{\beta_{f}x}-1)-\beta_{f}(e^{\beta_{c}x}-1)$
$\displaystyle=\sum_{n=1}^{\infty}\frac{x^{n}}{n!}(\beta_{c}\beta_{f}^{n}-\beta_{f}\beta_{c}^{n})$
$\displaystyle=\sum_{n=1}^{\infty}\frac{x^{n}}{n!}\beta_{c}\beta_{f}(\beta_{f}^{n-1}-\beta_{c}^{n-1})$
$\displaystyle<0~{}(\because~{}0<\beta_{f}<\beta_{c}),$

we have $dh(x)/dx<0$, which completes the proof.

## References

  * Mpemba and Osborne (1969) E. B. Mpemba and D. G. Osborne, Cool?, Phys. Educ. 4, 172 (1969).
  * Auerbach (1995) D. Auerbach, Supercooling and the Mpemba effect: When hot water freezes quicker than cold, Am. J. Phys. 63, 882 (1995).
  * Lasanta  _et al._ (2017) A. Lasanta, F. Vega Reyes, A. Prados, and A. Santos, When the hotter cools more quickly: Mpemba effect in granular fluids, Phys. Rev. Lett. 119, 148001 (2017).
  * Baity-Jesi  _et al._ (2019) M. Baity-Jesi, E. Calore, A. Cruz, L. A. Fernandez, J. M. Gil-Narvión, A. Gordillo-Guerrero, D. Iñiguez, A. Lasanta, A. Maiorano, E. Marinari, V. Martin-Mayor, J. Moreno-Gordo, A. M. Sudupe, D. Navarro, G. Parisi, S. Perez-Gaviro, F. Ricci-Tersenghi, J. J. Ruiz-Lorenzo, S. F. Schifano, B. Seoane, A. Tarancón, R. Tripiccione, and D. Yllanes, The Mpemba effect in spin glasses is a persistent memory effect, Proc. Natl. Acad. Sci. U.S.A. 116, 15350 (2019).
  * Kumar and Bechhoefer (2020) A. Kumar and J. Bechhoefer, Exponentially faster cooling in a colloidal system, Nature 584, 64 (2020).
  * Biswas  _et al._ (2020) A. Biswas, V. V. Prasad, O. Raz, and R. Rajesh, Mpemba effect in driven granular Maxwell gases, Phys. Rev. E 102, 012906 (2020).
  * Lu and Raz (2017) Z. Lu and O. Raz, Nonequilibrium thermodynamics of the Markovian Mpemba effect and its inverse, Proc. Natl. Acad. Sci. U.S.A. 114, 5083 (2017).
  * Klich  _et al._ (2019) I. Klich, O. Raz, O. Hirschberg, and M. Vucelja, Mpemba index and anomalous relaxation, Phys. Rev. X 9, 021060 (2019).
  * Takada  _et al._ (2021) S. Takada, H. Hayakawa, and A. Santos, Mpemba effect in inertial suspensions, Phys. Rev. E 103, 032901 (2021).
  * Gal and Raz (2020) A. Gal and O. Raz, Precooling strategy allows exponentially faster heating, Phys. Rev. Lett. 124, 060602 (2020).
  * Landi and Paternostro (2021) G. T. Landi and M. Paternostro, Irreversible entropy production: From classical to quantum, Rev. Mod. Phys. 93, 035008 (2021).
  * Alhambra and Woods (2017) A. M. Alhambra and M. P. Woods, Dynamical maps, quantum detailed balance, and the Petz recovery map, Phys. Rev. A 96, 022118 (2017).
  * Shiraishi and Saito (2019) N. Shiraishi and K. Saito, Information-theoretical bound of the irreversibility in thermal relaxation processes, Phys. Rev. Lett. 123, 110603 (2019).
  * Van Vu and Hasegawa (2021a) T. Van Vu and Y. Hasegawa, Lower bound on irreversibility in thermal relaxation of open quantum systems, Phys. Rev. Lett. 127, 190601 (2021a).
  * Van Vu and Hasegawa (2021b) T. Van Vu and Y. Hasegawa, Geometrical bounds of the irreversibility in Markovian systems, Phys. Rev. Lett. 126, 010601 (2021b).
  * Dattagupta (2012) S. Dattagupta,  _Relaxation Phenomena in Condensed Matter Physics_ (Elsevier, New York, 2012).
  * Benenti  _et al._ (2017) G. Benenti, G. Casati, K. Saito, and R. S. Whitney, Fundamental aspects of steady-state conversion of heat to work at the nanoscale, Phys. Rep. 694, 1 (2017).
  * Temme  _et al._ (2010) K. Temme, M. J. Kastoryano, M. B. Ruskai, M. M. Wolf, and F. Verstraete, The $\chi^{2}$-divergence and mixing times of quantum Markov processes, J. Math. Phys. 51, 122201 (2010).
  * Lapolla and Godec (2020) A. Lapolla and A. c. v. Godec, Faster uphill relaxation in thermodynamically equidistant temperature quenches, Phys. Rev. Lett. 125, 110602 (2020).
  * Meibohm  _et al._ (2021) J. Meibohm, D. Forastiere, T. Adeleke-Larodo, and K. Proesmans, Relaxation-speed crossover in anharmonic potentials, Phys. Rev. E 104, L032105 (2021).
  * Maes (2020) C. Maes, Frenesy: Time-symmetric dynamical activity in nonequilibria, Phys. Rep. 850, 1 (2020).
  * Seifert (2012) U. Seifert, Stochastic thermodynamics, fluctuation theorems and molecular machines, Rep. Prog. Phys. 75, 126001 (2012).
  * Schnakenberg (1976) J. Schnakenberg, Network theory of microscopic and macroscopic behavior of master equation systems, Rev. Mod. Phys. 48, 571 (1976).
  * Shiraishi  _et al._ (2018) N. Shiraishi, K. Funo, and K. Saito, Speed limit for classical stochastic processes, Phys. Rev. Lett. 121, 070601 (2018).
  * Vo  _et al._ (2020) V. T. Vo, T. Van Vu, and Y. Hasegawa, Unified approach to classical speed limit and thermodynamic uncertainty relation, Phys. Rev. E 102, 062132 (2020).
  * Margalit  _et al._ (2013) L. Margalit, M. Rosenbluh, and A. D. Wilson-Gordon, Degenerate two-level system in the presence of a transverse magnetic field, Phys. Rev. A 87, 033808 (2013).
  * Watabe  _et al._ (2020) S. Watabe, Y. Seki, and S. Kawabata, Enhancing quantum annealing performance by a degenerate two-level system, Sci. Rep. 10 (2020).
  * Tajima and Funo (2021) H. Tajima and K. Funo, Superconducting-like heat current: Effective cancellation of current-dissipation trade-off by quantum coherence, Phys. Rev. Lett. 127, 190604 (2021).
  * Manikandan (2021) S. K. Manikandan, Equidistant quenches in few level quantum systems, arXiv preprint arXiv:2102.06161 (2021).
  * Yang and Hou (2020) Z.-Y. Yang and J.-X. Hou, Non-Markovian Mpemba effect in mean-field systems, Phys. Rev. E 101, 052106 (2020).
  * Breuer and Petruccione (2002) H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, New York, 2002).
  * Nava and Fabrizio (2019) A. Nava and M. Fabrizio, Lindblad dissipative dynamics in the presence of phase coexistence, Phys. Rev. B 100, 125102 (2019).