††thanks: present address

# Multiple quantum Mpemba effect: exceptional points and oscillations

Amit Kumar Chatterjee<EMAIL_ADDRESS>Yukawa Institute for
Theoretical Physics, Kyoto University, Kitashirakawa Oiwakecho, Sakyo-ku,
Kyoto 606-8502, Japan  Department of Physics, Ramakrishna Mission
Vidyamandira, Belur Math, Howrah 711 202, India  Satoshi Takada
<EMAIL_ADDRESS>Department of Mechanical Systems Engineering and
Institute of Engineering, Tokyo University of Agriculture and Technology,
2–24–16 Naka-cho, Koganei, Tokyo 184–8588, Japan  Hisao Hayakawa
<EMAIL_ADDRESS>Yukawa Institute for Theoretical Physics, Kyoto
University, Kitashirakawa Oiwakecho, Sakyo-ku, Kyoto 606-8502, Japan

###### Abstract

We explore the role of exceptional points and complex eigenvalues on the
occurrence of the quantum Mpemba effect. To this end, we study a two-level
driven dissipative system subjected to an oscillatory electric field and
dissipative coupling with the environment. We find that both exceptional
points and complex eigenvalues can lead to multiple quantum Mpemba effect. It
occurs in an observable when time evolved copies corresponding to two
different initial conditions, one initially having higher observable value
compared to the other and both relaxing towards the same steady state,
intersect each other more than once during their relaxation process. Each of
the intersections denotes a quantum Mpemba effect and marks the reversal of
identities between the two copies i.e. the copy with higher observable value
before the intersection becomes the lower valued copy (and vice versa) after
the intersection. Such multiple intersections originate from additional
algebraic time dependence at the exceptional points and due to oscillatory
relaxation in the case of complex eigenvalues. We provide analytical results
for quantum Mpemba effect in the density matrix in presence of coherence.
Depending on the control parameters (drive and dissipation), observables such
as energy, von Neumann entropy, temperature etc. exhibit either single or
multiple quantum Mpemba effect. However, the distance from steady state
measured in terms of the Kullback-Leibler divergence shows only single quantum
Mpemba effect although the corresponding speed gives rise to either single or
multiple quantum Mpemba effect.

##  I Introduction

Memory effect describes the influence of past information on present events
Keim_2019 . Such memory effects often lead to unexpected events. In physical
systems, several forms of interesting memory effects have been observed in
different contexts e.g. pulse-sign memory in charge-density wave conductors
(memory of last input signal direction) Gill_1981 ; Fleming_1986 ; Kaiser
effect in rocks and metals Kaiser_1950 ; Zang_2014 ; Heimisson_2015  and
Mullins effect in rubbers Mullins_1969 ; Dorfmann_2004 ; Diani_2009  (both
bearing memories of largest applied stress); Kovacs effect (memory of bath
temperature that affects the thermalization process) Kovacs_1963 ; Kovacs_1979
; Bertin_2003 ; Prados_2014 ; Militaru_2021 ; aging and rejuvenation in spin
glasses Jonason_1998 ; Rodriguez_2003 ; Baity-Jesi_2018  and structural
glasses Zheng2003 ; Di2011 ; McKenna2017 ; Scalliet2019 , granular materials
Josserand_2000 ; Kumar16 ; Otsuki20 ; Kawasaki20 , dense suspensions
Pradipto20  and active matter Janssen_2017 .

The Mpemba effect (MPE) was originally conceived as a counter-intuitive
thermal phenomenon of faster cooling of hotter liquids, observed long ago
Aristotle81  and lately rediscovered by Mpemba and Osborne Mpemba69 . The MPE
thus bears the imprint of initial memory in a surprising way and constitutes a
separate class of fascinating memory effects. Later, the MPE has been regarded
as one of the typical anomalous relaxations where a system with initially
higher observable (e.g. temperature, energy, entropy, distance function etc.)
value can relax faster than its copy with initially lower observable value,
when both of them decrease towards the same steady state value. In classical
systems, MPE has been found in variety of classical systems including colloids
Kumar20 ; Kumar22 , granular fluids Santos17 ; Torrente19 ; Biswas20 ;
Biswas21 ; Mompo20 ; Megias22b ; Patron23 , optical resonators Keller18 ;
Santos20 ; Patron21 , inertial suspensions Takada21a ; Takada21b , Markovian
models Lu17 ; Klich19 ; Busiello21 ; Lin22 , Langevin system Biswas23a  and
many others Greaney11 ; Baity-Jesi19 ; Gonzalez2021 ; Yang20 ; Yang22 ;
Gonzalez21 ; Chetrite21 ; Holtzman22 ; Megias22a ; Biswas23 ; Pemartin23 .

In spite of considerable amount of works on MPE in the classical domain, only
a few studies have been done on the quantum Mpemba effect (QMPE), such as
quantum Ising model Nava19 , Dicke model Carollo21 , few-level systems
Manikandan21 , V-type system Ivander23 , XXZ spin chain Ares23 , single level
quantum dot Chatterjee_2023 , XY spin chain Murciano23  and free and
interacting integrable systems Rylands23 . In fact, only quite recently the
quantum analogue of the original thermal MPE i.e. QMPE in temperature, has
been theoretically revealed in a quantum dot Chatterjee_2023  and the
significant role of relaxation modes other than the slowest relaxation mode to
create QMPE, has been shown Chatterjee_2023 . However, the semi-classical
approach taken in Ref. Chatterjee_2023  naturally raises a question about the
generality of QMPE in temperature and other observable, in presence of
coherence in purely quantum mechanical systems. This constitutes one of the
motivations of the present work.

Since the phenomenon of QMPE is a comparatively new area of research, it is
important to explore the plausible connections of QMPE to other interesting
features of quantum systems. In the present work, we would like to address
some such connections. In particular, the study of exceptional points (EPs)
and surfaces in non-Hermitian quantum systems and Liouvillian dynamics has
caught a lot of attention in recent days Heiss12 ; Hashimoto15 ; Doppler16 ;
Zhong19 ; McDonald20 ; Soleymani22 ; Ding22 ; Minganti19 ; Arkhipov20 ;
Arkhipov21 ; Muller08 ; Naghiloo19 ; Saha23 ; Roccati23  along with their
applications in quantum optics Miri19 ; Ozdemir19 . The occurrence of EPs is
directly related to eigen-spectrum of a system and hence to its relaxation
process. Therefore a natural curiosity would be to understand the role of EPs
on the possibility of QMPE. In fact, the occurrence of double MPE (i.e. two
intersections between trajectories) at an EP has already been pointed out for
MPE in classical inertial suspensions Takada21a . Interestingly, the
acceleration of relaxation process in open quantum systems by maximizing the
gap between the zero mode and the slowest eigenmode at an EP, has recently
been studied Zhou23 . In general, such speed up is desired for several
application purposes and therefore it has been a subject of immense interest
in both classical and quantum systems Lapolla20 ; Vu21 ; Kochsiek22 ;
Ghimenti22 ; Bao22 ; Zhou23 . Since QMPE allows faster relaxation of initially
higher distant trajectories, it is important to find criteria for QMPE and to
understand how to utilize QMPE to achieve faster approach to steady state.
Another intriguing feature in non-Hermitian quantum mechanical systems is the
appearance of complex eigenvalues Bender2007 ; Moiseyev2011 ; Nakano2011 ;
Hatano2015 ; Li2023 . It is noteworthy that the analysis of QMPE in the
existing examples Nava19 ; Carollo21 ; Manikandan21 ; Ivander23 ;
Chatterjee_2023  is based on purely exponential relaxation processes, neither
the EP nor the oscillations have been considered yet. There is only an
exception in Ref. Ares23  where more than one intersections between
trajectories (two intersections to be precise, termed as revivals) has been
observed in the entanglement entropy in the XXZ spin chain, although not
discussed much and simply explained as a finite-size effect.

In this paper, we explore the following two topics: (i) the effect of EPs and
complex eigenvalues on QMPE by analysing the Quantum Master Equation, and (ii)
the existence of thermal QMPE in a purely quantum mechanical setting. We show
that in the presence of EPs and complex eigenvalues, two time-evolved
trajectories of an observable (e.g., temperature, energy) can intersect each
other more than once before reaching the same final steady state. In the
context of thermal QMPE, the first intersection implies that the initially
hotter quantum system relaxes faster to intersect the initially colder quantum
system and thereby reverses their roles, i.e. hotter becomes colder and vice
versa. This can be followed by one or more such intersections so that the
present hotter system can become colder again and so on. Similar observations
have been obtained in average energy. Such multiple intersections of
trajectories give rise to the remarkable phenomenon of multiple QMPE. To
demonstrate multiple QMPE as a consequence of EP and complex eigenvalues, we
study a two-level driven dissipative quantum system developed in Ref.
Hatano_2019 , which we call Hatano’s model in this paper, as a typical and
minimal model of open quantum systems. It is well known that many complicated
quantum problems are mapped or reduced to two-level systems; needless to
mention they are experimentally realizable (e.g., two internal electronic
energy levels of an atom controlled by laser) and have applications in qubits,
nuclear magnetic resonance, neutrino oscillations, etc. We also find thermal
QMPE in presence of coherence, some cases remarkably with multiple
intersections, following the same definition of temperature used in Ref.
Chatterjee_2023  . Importantly, we provide analytical results for the time
evolution of the density matrix in the two-level driven dissipative system
considered here. Consequently, we derive analytical criteria for QMPE in
energy and in principle, one can analyze exactly (analytically and
numerically) QMPE in any other observable.

The paper is arranged as follows. In Sec. II, we define the details of the
model and the protocol used to study QMPE. The Sec. III describes different
regions in the control parameter plane depending on the types of eigenvalues
of the Lindbladian and the multiple QMPE is defined. In articular, we focus on
regions with exceptional points and complex eigenvalues to discuss their
effects on QMPE in the density matrix and observables such as energy, entropy,
Kullback-Leibler (KL) divergence Kullback51 ; Sagawa22 ; Shiraishi23 ;
Strasberg22 ; Yoshii23  and temperature; by providing analytical results. In
Sec. IV, we discuss QMPE in a region with second-order EPs. Section V
demonstrates the effect of complex eigenvalues on the time evolution of the
density matrix leading to multiple QMPE in the observables. We summarize the
results and discuss possible future directions in Sec. VI. Explicit
expressions of the eigenvectors used in the analyses in Secs. IV and V, are
provided in Appendix A and Appendix B respectively. We discuss QMPE briefly in
region of purely exponential relaxation in Appendix C and at a third-order EP
in Appendix D.

##  II Model and protocol

We consider a two-level open quantum system, subjected to a periodic drive and
dissipatively coupled to the environment. Following Ref. Hatano_2019 , we
hereafter refer to this model as Hatano’s model in the present manuscript. In
this paper, we adopt $\hbar=1$. The two levels of the system, ground state and
excited state, have energies $E_{\mathrm{g}}$ and $E_{\mathrm{e}}$,
respectively. The gap between the two levels is denoted by
$\Delta:=E_{\mathrm{e}}-E_{\mathrm{g}}$. In Hatano’s model, the transfers
between the energy levels are enabled by the external oscillatory electric
field $E=E_{0}\mathrm{cos}(\omega t)$ and the dissipation of the system to the
environment via coupling $\Gamma$. The parameter $\delta:=\Delta-\omega$
represents the detuning of the periodic drive with respect to the system and
$d:=DE_{0}$, where $D$ is the electric dipole moment, denotes the effective
electric field. The tuning parameters in Hatano’s model are $d$, $\Gamma$ and
$\delta$ . However, one of them fixes the unit of energy, so we essentially
have two free parameters $\tilde{d}:=d/\delta$ and
$\tilde{\Gamma}:=\Gamma/\delta$.

Hatano’s model, within Markovian approximation, can be described by the
Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) equation Lindblad76 ; Gorini76 ;
Manzano20 ; Chruscinski17  (see Hatano_2019  for the present context):

$i\frac{d}{dt}|\widehat{\rho}(t)\rangle=\widehat{\mathcal{L}}|\widehat{\rho}(t)\rangle.$
(1)

The original $2\times 2$ density matrix is recast as the column vector
$|\widehat{\rho}(t)\rangle=(\rho_{\mathrm{eg}}(t),\rho_{\mathrm{ge}}(t),\rho_{\mathrm{ee}}(t),\rho_{\mathrm{gg}}(t))^{T}$
in Eq. (1), where $T$ denotes the transpose. The terms $\rho_{\mathrm{eg}}(t)$
and $\rho_{\mathrm{ge}}(t)$ are the off-diagonal density matrix elements,
characterizing the purely quantum mechanical nature of the system. The
elements $\rho_{\mathrm{gg}}(t)$ and $\rho_{\mathrm{ee}}(t)$ correspond to the
ground state and the excited state probabilities, respectively. For
convenience of notations, we denote the density matrix elements as $\rho_{j}$
where $j=1,2,3,4$ correspond to
$\rho_{\mathrm{eg}}(t),\rho_{\mathrm{ge}}(t),\rho_{\mathrm{ee}}(t),\rho_{\mathrm{gg}}(t)$
respectively. The transition between these density matrix elements is
described by the Lindbladian of the form below

$\widehat{\mathcal{L}}=\left(\begin{array}[]{cccc}1-i\tilde{\Gamma}/2&0&-\tilde{d}/2&\tilde{d}/2\\\
0&-1-i\tilde{\Gamma}/2&\tilde{d}/2&-\tilde{d}/2\\\
-\tilde{d}/2&\tilde{d}/2&-i\tilde{\Gamma}&0\\\
\tilde{d}/2&-\tilde{d}/2&i\tilde{\Gamma}&0\\\ \end{array}\right).$ (2)

Let us clarify a few points about the matrix $\widehat{\mathcal{L}}$ at this
point. Note that the dissipation $\Gamma$ makes the system non-Hermitian. More
importantly, each column sum $\widehat{\mathcal{L}}$ does not equal to zero,
which distinguishes Hatano’s model from the classical Markovian models studied
previously Lu17 . A recent work Chatterjee_2023  showing thermal quantum
Mpemba effect also deals with an effective 4-state quantum system. However,
the semi-classical approach in Ref. Chatterjee_2023  considered a diagonal
form of the density matrix, while the present model includes quantum coherence
in the form of off-diagonal density matrix elements $\rho_{\mathrm{eg}}(t)$
and $\rho_{\mathrm{ge}}(t)$. Naturally, it is interesting to ask whether
Hatano’s model, in the presence of coherence, can exhibit QMPE in observables
such as temperature, energy entropy and KL divergence.

Protocol.- We consider two different sets I and II of initial conditions
$\widehat{\rho}^{\mathrm{I}}(t=0)$ and $\widehat{\rho}^{\mathrm{II}}(t=0)$.
Both of them are chosen in the form of steady state distributions
corresponding to the zero eigenvalue of $\widehat{\mathcal{L}}$. The
difference between the two initial conditions lies in the choice of the
control parameters. To elaborate, we choose
$(\tilde{d}_{\mathrm{I}},\tilde{\Gamma}_{\mathrm{I}})$ for initial condition I
and $(\tilde{d}_{\mathrm{II}},\tilde{\Gamma}_{\mathrm{II}})$ for initial
condition II. Both of them are instantaneously quenched to the same set of
values $(\tilde{d},\tilde{\Gamma})$, so that they will eventually reach the
same steady state. Thus, the entire set of control parameters consists of six
parameters
$(\tilde{d}_{\mathrm{I}},\tilde{\Gamma}_{\mathrm{I}},\tilde{d}_{\mathrm{II}},\tilde{\Gamma}_{\mathrm{II}};\tilde{d},\tilde{\Gamma})$,
where the first four are before-quench parameters and the last two are after-
quench parameters. In most parts of the present work, we perform fixed
dissipation protocol
$\tilde{\Gamma}_{\mathrm{I}}=\tilde{\Gamma}_{\mathrm{II}}=\tilde{\Gamma}$
where we only tune the drive to obtain QMPE. However, in some cases we will
use variable dissipation protocol
$\tilde{\Gamma}_{\mathrm{I}}\neq\tilde{\Gamma}_{\mathrm{II}}\neq\tilde{\Gamma}$
where both drive and dissipation are tuned to search for QMPE in the whole
parameter space.

To mathematically formulate the initial conditions before quench, let us
denote the right eigenvectors and left eigenvectors of $\widehat{\mathcal{L}}$
[Eq. (2)] by $|r_{j}\rangle$ and $\langle\ell_{j}|$ ($j=1,2,3,4$),
respectively. The vector $|r_{1}\rangle$ corresponds to the steady probability
distribution and $\langle\ell_{1}|$ is the corresponding left eigenvector.
Below we provide their explicit forms:

$\displaystyle|r_{1}\rangle$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}-\displaystyle\frac{\tilde{d}(2+\mathrm{i}\tilde{\Gamma})}{4+\tilde{d}^{2}+\tilde{\Gamma}^{2}},&-\displaystyle\frac{\tilde{d}(2-\mathrm{i}\tilde{\Gamma})}{4+\tilde{d}^{2}+\tilde{\Gamma}^{2}},&\displaystyle\frac{\tilde{d}^{2}}{4+\tilde{d}^{2}+\tilde{\Gamma}^{2}},&1\\\
\end{array}\right)^{T},$ (4) $\displaystyle\langle\ell_{1}|$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}0,&0,&\displaystyle\frac{4+\tilde{d}^{2}+\tilde{\Gamma}^{2}}{4+2\tilde{d}^{2}+\tilde{\Gamma}^{2}},&\displaystyle\frac{4+\tilde{d}^{2}+\tilde{\Gamma}^{2}}{4+2\tilde{d}^{2}+\tilde{\Gamma}^{2}}\\\
\end{array}\right).$ (6)

The steady state
$|\widehat{\rho}_{\mathrm{ss}}(\tilde{d},\tilde{\Gamma})\rangle$ is simply
related to $|r_{1}\rangle$ as

$\displaystyle|\widehat{\rho}_{\mathrm{ss}}(\tilde{d},\tilde{\Gamma})\rangle=\frac{4+\tilde{d}^{2}+\tilde{\Gamma}^{2}}{4+2\tilde{d}^{2}+\tilde{\Gamma}^{2}}|r_{1}\rangle,$
(7)

where $\tilde{d}$ and $\tilde{\Gamma}$ are drive and dissipation after the
quench, respectively. We construct the two initial conditions I and II as
follows:

$\displaystyle|\widehat{\rho}^{\mathrm{I}}(t=0)\rangle=|\widehat{\rho}_{\mathrm{ss}}(\tilde{d}_{\mathrm{I}},\tilde{\Gamma}_{\mathrm{I}})\rangle,$
(8)
$\displaystyle|\widehat{\rho}^{\mathrm{II}}(t=0)\rangle=|\widehat{\rho}_{\mathrm{ss}}(\tilde{d}_{\mathrm{II}},\tilde{\Gamma}_{\mathrm{II}})\rangle.$
(9)

The coherence terms $\rho_{\mathrm{eg}}$ and $\rho_{\mathrm{ge}}$ being
complex conjugates of each other, we decompose them into real and imaginary
parts as

$\displaystyle\rho_{\mathrm{eg}}=\rho_{\mathrm{re}}+\mathrm{i}\rho_{\mathrm{im}},\hskip
5.69046pt\rho_{\mathrm{ge}}=\rho_{\mathrm{re}}-\mathrm{i}\rho_{\mathrm{im}}.$
(10)

Specifically, in steady state, $\rho_{\mathrm{re}}$ and $\rho_{\mathrm{im}}$
have the forms (using Eqs. (6) and (7)):

$\displaystyle\rho_{\mathrm{re}}=-\frac{2\tilde{d}}{4+2\tilde{d}^{2}+\tilde{\Gamma}^{2}},\hskip
5.69046pt\rho_{\mathrm{im}}=-\frac{\tilde{d}\,\tilde{\Gamma}}{4+2\tilde{d}^{2}+\tilde{\Gamma}^{2}}.$
(11)

##  III Different regimes in $\tilde{d}-\tilde{\Gamma}$ plane and multiple
QMPE

As discussed in Ref. Hatano_2019 , Hatano’s model exhibits a rich spectrum of
eigenvalues. In the plane of control parameters $\tilde{d}-\tilde{\Gamma}$, it
has six different regions ($\mathrm{a}_{1}$), ($\mathrm{a}_{2}$), (b), (c),
(d) and (e) depending on the nature of the eigenvalues of
$\widehat{\mathcal{L}}$ [Eq. (2)]. The regions ($\mathrm{a}_{1}$) and
($\mathrm{a}_{2}$) consist of complex eigenvalues, region ($\mathrm{b}$) with
purely imaginary eigenvalues, regions ($\mathrm{c}$) and ($\mathrm{d}$) with
second-order EPs and region ($\mathrm{e}$) containing a single third-order EP.
At a $n$-th order EP, $n$ eigenvalues become equal and their eigenvectors
coalesce. The regions ($\mathrm{a_{1}}$) and ($\mathrm{a_{2}}$) can be
distinguished from each other depending on the largest nonzero eigenvalue
which is purely imaginary in ($\mathrm{a_{1}}$) but complex in
($\mathrm{a_{2}}$). They can be also differentiated by their locations in the
parameter space of $(\tilde{d},\tilde{\Gamma})$. The region ($\mathrm{a_{1}}$)
is placed in relatively larger $\tilde{d}$ and smaller $\tilde{\Gamma}$,
whereas the region ($\mathrm{a_{2}}$) is placed in comparatively smaller
$\tilde{d}$ and larger $\tilde{\Gamma}$. All the nonzero eigenvalues in region
(b) are purely imaginary and they are unequal from one another. The line (c)
of second-order EPs separate region ($\mathrm{a_{1}}$) from region (b),
whereas the line (d) of second-order EPs demarcate (b) from
($\mathrm{a_{2}}$). The two lines (c) and (d) of second-order EPs meet at a
single point in the $\tilde{d}-\tilde{\Gamma}$ plane that defines the region
(e) which is a third-order EP. See Figs. 3 and 4 in Ref. Hatano_2019  for the
visual descriptions of these regions. We would like to analyze the full
relaxation of the system that contains all three relaxation modes (apart from
the eigenvalue zero that corresponds to the steady state) of the Lindbladian
$\widehat{\mathcal{L}}$ [Eq.(2)]. The presence of multiple relaxation modes,
EPs and complex eigenvalues naturally triggers the possibility of interesting
features of QMPE in Hatano’s model.

Indeed, in the space of control parameters
$(\tilde{d}_{\mathrm{I}},\tilde{\Gamma}_{\mathrm{I}},\tilde{d}_{\mathrm{II}},\tilde{\Gamma}_{\mathrm{II}};\tilde{d},\tilde{\Gamma})$,
we find multiple QMPE where the two trajectories corresponding to some
observable, intersect each other multiple times before reaching the same
steady state. Specifically, we have obtained two or more intersections in a
oscillatory manner, in regions ($\mathrm{a}_{1}$) and ($\mathrm{a}_{2}$) with
complex eigenvalues. Whereas, at most two intersections or double QMPE can be
found in regions (c), (d) and (e) with EP. Interestingly, even number of
intersections (e.g. double QMPE or quadruple QMPE) imply that the initial
identities of the two copies (higher valued copy and lower valued copy of an
observable, more explicitly hotter and colder copies if the observable is
temperature) are restored finally after the last intersection. Therefore, in
such cases the initially hotter or higher valued copy remains hotter or higher
valued after all the intersections have occurred. On the other hand, odd
number of intersections (single QMPE or triple QMPE) indicates the ultimate
reversal of the initial identities after the final intersection. In these
cases the initially hotter or higher valued copy finally becomes colder or
lower valued copy after the completion of all intersections. This
classification has important consequences when one is not only interested in
faster relaxations or intersections at intermediate times, but also interested
in faster reach to final or steady state. It would be intriguing in future to
find connections between such odd or even QMPE and faster approach to steady
state. Notably, if we restrict on monotonic time volution For that case, even
number of QMPE is of no use. Only odd number of QMPE can help to reach faster
to steady state using QMPE. However, more conditions are required like
avoiding overshoot etc. to use odd QMPE to get faster reach to the steady
state. In the next sections, we discuss multiple QMPE in different regimes in
detail.

In this paper, we illustrate QMPEs in region (d) for quench onto the second-
order EP in Sec. IV (with more analytical details in Appendix A) and region
($\mathrm{a_{1}}$) for the quench onto the point with complex eigenvalues in
Sec. V (with more analytical details in Appendix B). We do not discuss QMPE in
regions (c) and ($\mathrm{a_{2}}$) separately, because the analyses of QMPEs
in these regions are similar to those in regions (d) and ($\mathrm{a_{1}}$),
respectively. The QMPEs in other regions are discussed in the appendices.
Appendix C discusses QMPE in region (b) with purely exponential relaxations
and Appendix D contains the description of the QMPE at the third-order EP in
region (e).

##  IV QMPE in region ($\mathrm{d}$): second-order EP

In this section, we consider the possibility of single and multiple QMPE in
region (d) of Hatano’s model. This region contains second order EP where two
of the nonzero eigenvalues of the Lindbladian $\widehat{\mathcal{L}}$ [Eq.
(2)] become same and their eigenvectors coalesce. In the
$(\tilde{d},\tilde{\Gamma})$ parameter plane, the region (d) is defined by the
set of points

$\displaystyle\tilde{\Gamma}=\sqrt{\frac{\tilde{d}^{4}}{2}+10\tilde{d}^{2}-4+\frac{\tilde{d}}{2}(\tilde{d}^{2}-8)^{3/2}}.$
(12)

The eigenvalues of $\widehat{\mathcal{L}}$ are denoted by
$-\mathrm{i}\lambda_{j}$ ($j=1,2,3,4$) where $\lambda_{1}=0$. Specifically, in
region (d), the eigenvalues obey $0<\lambda_{2}=\lambda_{3}<\lambda_{4}$. To
analyze QMPE analytically, we need to diagonalize $\widehat{\mathcal{L}}$. The
diagonalization is usually done by the relation
$\widehat{\mathcal{L}}_{d}=\widehat{L}\widehat{\mathcal{L}}\widehat{R}$, where
$\widehat{\mathcal{L}}_{d}$ is the diagonal form of the Lindbladian and the
matrix $\widehat{R}$ consists of the right eigenvectors $|r_{j}\rangle$ of
$\widehat{\mathcal{L}}$ as its columns and the matrix $\widehat{L}$ contains
the left eigenvectors $\langle\ell_{j}|$ as its rows, $j=1,2,3,4$. More
explicitly,

$\displaystyle\widehat{R}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}|r_{1}\rangle,&|r_{2}\rangle,&|r_{3}\rangle,&|r_{4}\rangle\\\
\end{array}\right),$ (14) $\displaystyle\widehat{L}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}\langle\ell_{1}|,&\langle\ell_{2}|,&\langle\ell_{3}|,&\langle\ell_{4}|\\\
\end{array}\right)^{T}.$ (16)

The explicit expressions for $|r_{k}\rangle$ and $\langle\ell_{k}|$
$(k=1,2,3,4)$ [Eq. (16)] in region (d) are provided in Appendix A.

It is known that, at the second-order EPs, such completely diagonal form
$\widehat{\mathcal{L}}_{d}$ cannot be obtained. Instead, we achieve a Jordan
normal form $\widehat{\mathcal{L}}_{\mathrm{J}}$ that is not completely
diagonal, through the equation below

$\widehat{\mathcal{L}}_{\mathrm{J}}=\widehat{L}\widehat{\mathcal{L}}\widehat{R},$
(17)

where $\widehat{\mathcal{L}}_{\mathrm{J}}$ has the following form

$\displaystyle\widehat{\mathcal{L}}_{\mathrm{J}}=\left(\begin{array}[]{cccc}0&0&0&0\\\
0&-\mathrm{i}\lambda_{2}&1&0\\\ 0&0&-\mathrm{i}\lambda_{2}&0\\\
0&0&0&-\mathrm{i}\lambda_{4}\\\ \end{array}\right).$ (22)

Note the presence of the off-diagonal unit element in Eq. (22). The explicit
expressions for $\lambda_{2}$ and $\lambda_{4}$ are obtained as

$\displaystyle\lambda_{2}$ $\displaystyle=$
$\displaystyle\frac{2\tilde{\Gamma}}{3}+2\mathrm{cos}\left(\frac{2\pi}{3}\right)\left(\frac{\tilde{\Gamma}}{6}\left(1-\frac{\tilde{d}^{2}}{2}+\frac{\tilde{\Gamma}^{2}}{36}\right)\right)^{1/3},$
(23) $\displaystyle\lambda_{4}$ $\displaystyle=$
$\displaystyle\frac{2\tilde{\Gamma}}{3}+2\left(\frac{\tilde{\Gamma}}{6}\left(1-\frac{\tilde{d}^{2}}{2}+\frac{\tilde{\Gamma}^{2}}{36}\right)\right)^{1/3}.$
(24)

We primarily seek the exact analytical solution for the time evolution of
$|\widehat{\rho}(t)\rangle$ with elements $\rho_{j}(t)$ ($j=1,2,3,4$
correspond to
$\rho_{\mathrm{eg}},\rho_{\mathrm{ge}},\rho_{\mathrm{ee}},\rho_{\mathrm{gg}}$
respectively). In principle, once the time evolution of the density matrix is
obtained, one can calculate any other observable of interest (e.g. energy,
entropy, temperature, the Kullback-Leibler (KL) divergence etc.) to explore
QMPE.

###  IV.1 Density matrix elements

Figure 1: The figure shows the intersection times $\tilde{t}_{0}$ and
$\tilde{t}_{-1}$ from Eq. (33), as functions of the control parameters
$\tilde{d}_{\mathrm{I}}$ and $\tilde{d}_{\mathrm{II}}$. The regions with only
one solution $\tilde{t}_{0}$ indicate single QMPE in the ground state
probability $\rho_{\mathrm{gg}}(t)$ whereas the regions with both solutions
imply double QMPE. Parameters used are
$\tilde{d}=4.0,\tilde{\Gamma}_{\mathrm{I}}=\tilde{\Gamma}_{\mathrm{II}}=\tilde{\Gamma}=\sqrt{(568+64\sqrt{2})/2}$
[Eq. (12)].

Using Eqs. (17) and (22) in Eq. (1), we have obtained the following
expressions of the density matrix elements $\rho_{j}(t)$ $(j=1,2,3,4)$:

$\displaystyle\rho_{j}(t)=\sum_{k=1}^{4}e^{-\lambda_{k}t}r_{k,j}a_{k}-\mathrm{i}te^{-\lambda_{2}t}r_{2,j}a_{3},$
(25) $\displaystyle a_{k}=\sum_{n=1}^{4}\ell_{k,n}\rho_{n}(0),$ (26)

where the initial density matrix elements $\rho_{n}(0)$ before quench are
given by Eq. (9) and the elements $\ell_{k,n}$ are provided in Appendix A.
Interestingly, in the presence of the second-order EP in region (d), the
density matrix elements in Eq. (26) have picked up an additional algebraic
term proportional to time $t$. This is remarkably different from the time
evolution in regions (b) which has only purely exponential relaxations [see
Appendix B] and region ($\mathrm{a_{1}}$) that contains damped oscillations
[see Sec. V]. The linear $t$ dependence originates from the off-diagonal
defect in the Jordan normal form in Eq. (22).

To investigate QMPE, we focus on the ground state density matrix element
$\rho_{\mathrm{gg}}(t)$ $(\equiv\rho_{4}(t))$. To detect the QMPE, we would
like to find the intersections of the two copies I and II during their
relaxation, i.e. solution(s) to the equation $\Delta\rho_{\mathrm{gg}}(t)=0$.
We obtain the following form of $\Delta\rho_{\mathrm{gg}}(t)$

$\displaystyle\Delta\rho_{\mathrm{gg}}(t)=-e^{-\lambda_{2}t}\left[\alpha_{1}e^{-(\lambda_{4}-\lambda_{2})t}+t\,\alpha_{2}+\alpha_{3}\right],$
(27)

where the co-efficients $\alpha_{1},\alpha_{2},\alpha_{3}$ are given by

$\displaystyle\alpha_{1}=a^{\mathrm{I}}_{4}-a^{\mathrm{II}}_{4},\,\,\alpha_{2}=-\mathrm{i}(a^{\mathrm{I}}_{3}-a^{\mathrm{II}}_{3}),\,\,\alpha_{3}=a^{\mathrm{I}}_{2}-a^{\mathrm{II}}_{2}.$
(28)

The solutions to such transcendental equation in Eq. (27) are given by the
Lambert $W$ function, where the equation $ye^{y}=x$ has only one solution
$y=W_{0}(x)$ if $x\geqslant 0$ and two solutions $y=W_{0}(x)$ and
$y=W_{-1}(x)$ if $-1/e\leqslant x<0$. Indeed, such properties of Lambert $W$
function imply that there can be more than one solutions to
$\Delta\rho_{\mathrm{gg}}(t)=0$, meaning multiple QMPE. More precisely, at
most two intersection times i.e. double QMPE can be observed in region (d).
Below we provide the criteria for double QMPE and single QMPE in this region:

Figure 2: The figure illustrates the intersection times $\tilde{t}^{E}_{0}$
and $\tilde{t}^{E}_{-1}$ [Eq. (LABEL:eq:d14)] in energy $E(t)$. The parameter
regions with only one solution $\tilde{t}^{E}_{0}$ give single QMPE in energy
whereas regions with both solutions mark double QMPE. In contrast to
intersection times in $\rho_{\mathrm{gg}}(t)$ [Fig. 1], here there are regions
where QMPE is absent in energy. Parameters used are same as Fig. 1.

$\displaystyle(\lambda_{4}-\lambda_{2})\frac{\alpha_{1}}{\alpha_{2}}e^{(\lambda_{4}-\lambda_{2})\frac{\alpha_{3}}{\alpha_{2}}}>\frac{1}{e}:$
$\displaystyle\hskip 2.84544pt\mathrm{no}\,\,\mathrm{QMPE},$ (29)
$\displaystyle(\lambda_{4}-\lambda_{2})\frac{\alpha_{1}}{\alpha_{2}}e^{(\lambda_{4}-\lambda_{2})\frac{\alpha_{3}}{\alpha_{2}}}>\frac{1}{e}\,\&\,\frac{\alpha_{1}}{\alpha_{2}}>0:$
$\displaystyle\hskip 2.84544pt\mathrm{double}\,\,\mathrm{QMPE},$ (30)
$\displaystyle\frac{\alpha_{1}}{\alpha_{2}}<0:$ $\displaystyle\hskip
2.84544pt\mathrm{usual}\,\,\mathrm{QMPE}.$

The two possible intersection times for $\Delta\rho_{\mathrm{gg}}(t)=0$ in the
form of Lambert W functions are given below:

$\displaystyle\tilde{t}_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{(\lambda_{4}-\lambda_{2})}W_{0}\left((\lambda_{4}-\lambda_{2})\frac{\alpha_{1}}{\alpha_{2}}e^{(\lambda_{4}-\lambda_{2})\frac{\alpha_{3}}{\alpha_{2}}}\right)-\frac{\alpha_{3}}{\alpha_{2}},$
(32) $\displaystyle\tilde{t}_{-1}$ $\displaystyle=$
$\displaystyle\frac{1}{(\lambda_{4}-\lambda_{2})}W_{-1}\left((\lambda_{4}-\lambda_{2})\frac{\alpha_{1}}{\alpha_{2}}e^{(\lambda_{4}-\lambda_{2})\frac{\alpha_{3}}{\alpha_{2}}}\right)-\frac{\alpha_{3}}{\alpha_{2}}.\hskip
17.07182pt$ (33)

The asymptotic expansions of the functions $W_{0}(x)$ and $W_{-1}(x)$ can be
found in Ref. Corless96 . Note that Lambert $W$ functions have found
applications in several fields such as viscous flows More06 , crystal growth
Asadian14 , epitaxial film growth Braun03 , phase separation of polymer
mixtures Bot21 , criticality in temporal networks Moran23  etc.

In Fig. 1, we plot the solutions $\tilde{t}_{0}$ and $\tilde{t}_{-1}$ [Eq.
(33)] as functions of the driving strengths before quench. Indeed, we observe
that there are regions with single QMPE where only $\tilde{t}_{0}$ exists and
there are other regions where both solutions prevail leading to double QMPE.
In particular, we observe double QMPE for large values of
$\tilde{d}_{\mathrm{I}}$ and $\tilde{d}_{\mathrm{II}}$. Note that QMPE appears
to be unavoidable in the ground state probability $\rho_{\mathrm{gg}}(t)$ in
Fig. 1 i.e. there is no such region in the
$(\tilde{d}_{\mathrm{I}},\tilde{d}_{\mathrm{II}})$ where there is no QMPE.

###  IV.2 Energy

Here we analyze QMPE in average energy $E(t)=\mathrm{Tr}[\widehat{\rho}(t)H]$
where $H$ is the system Hamiltonian presented as a $2\times 2$ matrix with
elements $H_{11}=1,H_{22}=0,H_{12}=H_{21}=\tilde{d}/2$ Hatano_2019 .
Consequently, $E(t)$ can be expressed as

$\displaystyle
E(t)=1-\rho_{4}(t)+\frac{\tilde{d}}{2}(\rho_{1}(t)+\rho_{2}(t)).$ (34)

Figure 3: The figure shows an example where both ground state probability
[(a)] and energy [(b)] exhibit double QMPE. We have plotted the difference
between two copies I and II in each case. For (a), the second QMPE in
$\rho_{\mathrm{gg}}(t)$ is much weaker in magnitude whereas for (b) both QMPE
in $E(t)$ are of comparative magnitudes Parameters used are
$\tilde{d}=4.0,\tilde{\Gamma}_{\mathrm{I}}=\tilde{\Gamma}_{\mathrm{II}}=\tilde{\Gamma}=\sqrt{(568+64\sqrt{2})/2},\tilde{d}_{\mathrm{I}}=10.0,\tilde{d}_{\mathrm{II}}=12.0$.

Note that the average energy directly captures the effect of coherence due to
the presence of the off-diagonal density matrix elements in Eq. (34). To
explore the QMPE in $E(t)$, we consider the time evolution of the energies
starting from the two different initial conditions I and II and look for their
intersections during the relaxations. More precisely, we would like to find
the solutions to $\Delta E(t)=0$ where $\Delta
E(t)=E_{\mathrm{I}}(t)-E_{\mathrm{II}}(t)$. We obtain the following expression

$\displaystyle\Delta
E(t)=-e^{-\lambda_{2}t}\left[\gamma_{1}e^{-(\lambda_{4}-\lambda_{2})t}+t\,\gamma_{2}+\gamma_{3}\right],$
(35)

where the coefficients $\gamma_{1},\gamma_{2},\gamma_{3}$ are given by

$\displaystyle\gamma_{1}$ $\displaystyle=$
$\displaystyle\alpha_{1}\left(1+\frac{\tilde{d}^{2}}{1+(\Gamma/2-\lambda_{4})^{2}}\right),$
(36) $\displaystyle\gamma_{2}$ $\displaystyle=$
$\displaystyle\alpha_{2}\left(1+\frac{\tilde{d}^{2}}{1+(\Gamma/2-\lambda_{2})^{2}}\right),$
(37) $\displaystyle\gamma_{3}$ $\displaystyle=$
$\displaystyle\alpha_{3}\left(1+\frac{\tilde{d}^{2}}{1+(\Gamma/2-\lambda_{2})^{2}}\right)-\alpha_{2}(\Gamma/2-\lambda_{4}),\,\,\,\,$
(38)

where the elements $\alpha_{1},\alpha_{2},\alpha_{3}$ are given by Eq. (28).
Note that Eq. (38) is a similar kind of transcendental equation as Eq. (27),
only with different coefficients. Therefore, we can analyze the possibility of
QMPE in energy using Lambert $W$ functions as we did for the ground state
probability. The two possible solutions to $\Delta E(t)=0$ are given by

$\displaystyle\tilde{t}^{E}_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{(\lambda_{4}-\lambda_{2})}W_{0}\left((\lambda_{4}-\lambda_{2})\frac{\gamma_{1}}{\gamma_{2}}e^{(\lambda_{4}-\lambda_{2})\frac{\gamma_{3}}{\gamma_{2}}}\right)-\frac{\gamma_{3}}{\gamma_{2}},$
(39) $\displaystyle\tilde{t}^{E}_{-1}$ $\displaystyle=$
$\displaystyle\frac{1}{(\lambda_{4}-\lambda_{2})}W_{-1}\left((\lambda_{4}-\lambda_{2})\frac{\gamma_{1}}{\gamma_{2}}e^{(\lambda_{4}-\lambda_{2})\frac{\gamma_{3}}{\gamma_{2}}}\right)-\frac{\gamma_{3}}{\gamma_{2}},$

where the superscript $E$ denotes energy and is used to differentiate the
solutions from Eq. (33) obtained for ground state probability.

Figure 2 exhibits the solutions in Eq. (LABEL:eq:d14) as functions of the
control parameters $(\tilde{d}_{\mathrm{I}},\tilde{d}_{\mathrm{II}})$. We
observe region with single QMPE (only $\tilde{t}^{E}_{0}$ exists) as well as
region for double QMPE (both $\tilde{t}^{E}_{0}$ and $\tilde{t}^{E}_{-1}$
exist). In addition, we observe regions where both the solutions are absent
i.e. no QMPE. This crucial difference in Fig. 2 compared to Fig. 1 indicates
that QMPE in ground state probability does not necessarily implies QMPE in
average energy. However, Figs. 1 and 2 also show that there are common
parameter values that lead to QMPE both in $\rho_{\mathrm{gg}}(t)$ and $E(t)$.
We present an explicit example in Fig. 3 where both of these quantities
exhibit double QMPE, for the same set of parameters. However, the two
intersection times for $\rho_{\mathrm{gg}}(t)$ and $E(t)$ are different from
each other. The intersection times for $\rho_{\mathrm{gg}}(t)$ are
comparatively far from each other with the second QMPE being much weaker in
magnitude [Fig. 3(a)], whereas the intersection times for $E(t)$ are
comparatively much closer to each other, with both the first and second QMPE
being comparable in magnitude [Fig. 3(b)].

###  IV.3 von Neumann entropy

Next we consider the von Neumann entropy $S_{\mathrm{vN}}(t)$, which is a
function of $\widehat{\rho}(t)$, defined as

$\displaystyle
S_{\mathrm{vN}}(t):=-\mathrm{Tr}[\widehat{\rho}(t)\,\mathrm{ln}(\widehat{\rho}(t))].$
(41)

Figure 4: The figure illustrates double QMPE in von Neumann entropy for the
same set of parameters used in Fig. 3. The inset shows the second intersection
where the magnitude of the second QMPE is much weaker than the first QMPE
shown in the main figure.

The nonlinearity of $S_{\mathrm{vN}}(t)$ as a function of the density matrix
elements poses difficulties to perform analytical treatment of the QMPE in von
Neumann entropy, in contrast to what we could achieve previously for ground
state probability and energy. To observe QMPE in $S_{\mathrm{vN}}(t)$, we
consider the intersection time(s) for the quantity $\Delta
S_{\mathrm{vN}}(t)=S^{\mathrm{I}}_{\mathrm{vN}}-S^{\mathrm{II}}_{\mathrm{vN}}$.
In Fig. 4, we present the behavior of $\Delta S_{\mathrm{vN}}(t)$ for the same
set of parameter values used in Fig. 3. Indeed, we observe that $\Delta
S_{\mathrm{vN}}(t)$ intersects zero twice, implying the existence of double
QMPE in the von Neumann entropy, similar to $\rho_{\mathrm{gg}}(t)$ and
$E(t)$.

###  IV.4 Temperature

Figure 5: The figure shows thermal QMPE. We have plotted the difference in
temperatures between two copies I and II. After the two copies intersect each
other, the initially hotter copy becomes colder leading to QMPE in
temperature. Parameters used are
$\tilde{d}=4.0,\tilde{\Gamma}=\sqrt{(568+64\sqrt{2})/2},\tilde{\Gamma}_{\mathrm{I}}=15.0,\tilde{\Gamma}_{\mathrm{II}}=2.1,\tilde{d}_{\mathrm{I}}=3.0,\tilde{d}_{\mathrm{II}}=22.9$.

We use the definition of temperature $T(t)$ following Ref. Chatterjee_2023  to
examine the possibility of thermal QMPE in Hatano’s model, as given below

$\displaystyle T(t):=\frac{\partial E(t)}{\partial
S_{\mathrm{vN}}(t)}=\left.\frac{\partial E(t)}{\partial
t}\right/\frac{\partial S_{\mathrm{vN}}(t)}{\partial t},$ (42)

where $E(t)$ and $S_{\mathrm{vN}}(t)$ are defined in Eqs. (34) and (41),
respectively. The QMPE in system temperature has been observed recently in
Ref. Chatterjee_2023  using a semi-classical approach. Naturally, it would be
interesting to explore thermal QMPE in Hatano’s model that involves the effect
of off-diagonal density matrix elements. However, we find that there is no
thermal QMPE for the same set of parameters that have produced QMPE for ground
state probability [Fig. 3(a)], energy [Fig. 3(b)] and entropy [Fig. 4]. To
search for thermal QMPE, we consider a variable dissipation protocol that
allows us to use
$\tilde{\Gamma}_{\mathrm{I}}\neq\tilde{\Gamma}_{\mathrm{II}}\neq\tilde{\Gamma}$
(wider parameter space compared to the fixed dissipation protocol
$\tilde{\Gamma}_{\mathrm{I}}=\tilde{\Gamma}_{\mathrm{II}}=\tilde{\Gamma}$).
Indeed, using such protocol, we have obtained the thermal QMPE shown in Fig. 5
where we present the temperature difference $\Delta T(t)$ between two copies I
and II. We observe that $\Delta T(t)$ becomes zero at some intermediate
relaxation time that denotes the intersection of two copies, after which the
initially hotter copy becomes colder, leading to thermal QMPE.

###  IV.5 Kullback-Leibler divergence

The Kullback-Leibler (KL) divergence is a measure of difference from steady
state which is a monotonically decreasing function of time, and it is given by

$\displaystyle
D_{\mathrm{KL}}(t):=\mathrm{Tr}[\widehat{\rho}(t)(\mathrm{ln}(\widehat{\rho}(t))-\mathrm{ln}(\widehat{\rho}_{\mathrm{ss}}))],$
(43)

where $\widehat{\rho}_{\mathrm{ss}}$ is the steady state density matrix.

Figure 6: The figure (a) shows QMPE in KL divergence [Eq. (43)] where the
initially higher distant copy (here I) reaches to steady state faster than the
initially lower distance copy (here II). The figure (b) demonstrates the
relaxation speed [Eq. (44)] of these two copies. The speed also exhibit QMPE,
but the intersection time is different compared to $D-{\mathrm{KL}}(t)$.
Parameters used for both figures are
$\tilde{d}=4.0,\tilde{\Gamma}=\sqrt{(568+64\sqrt{2})/2},\tilde{\Gamma}_{\mathrm{I}}=20.0,\tilde{\Gamma}_{\mathrm{II}}=4.5,\tilde{d}_{\mathrm{I}}=10.0,\tilde{d}_{\mathrm{II}}=2.0$.

Although the KL divergence is not a distance in strict metric sense, however
it is well-accepted as a statistical distance. Such distance function has been
widely used to investigate MPE in classical systems Lu17 ; Kumar20 ; Klich19
since it is suitable to characterize the approach to steady state because of
its monotonic nature. The KL divergence Kullback51  plays fundamental roles in
non-equilibrium thermodynamics Sagawa22 ; Shiraishi23 ; Strasberg22 ; Yoshii23
and has diverse applications in various fields including information geometry
Amari20 , machine learning Burnham02 , game theory Soklakov20  etc.

To obtain QMPE in KL divergence, we use the variable dissipation protocol. In
Fig. 6(a) we present an example where the copy starting at a higher distance
from steady state relaxes faster compared to the other copy with lower initial
distance. The copies I and II intersect each other leading to QMPE in
$D_{\mathrm{KL}}(t)$ and the copy starting at higher initial distance reaches
the steady state faster. We have not observed multiple QMPE in KL divergence.
It would be an interesting future question to understand if the absence of
multiple QMPE in $D_{\mathrm{KL}}(t)$ is its generic feature and if the reason
is its monotonicity.

An important variable of interest in any relaxation process, be it classical
or quantum, is the speed of the relaxation. Since the KL divergence gives the
distance measure from steady state, a relevant relaxation speed can be defined
as

$\displaystyle v_{\mathrm{KL}}(t)=-\frac{\partial D_{\mathrm{KL}}(t)}{\partial
t},$ (44)

where the negative sign is used to complement the fact that the KL divergence
decreases with time. It would be interesting to see if the QMPE in
$D_{\mathrm{KL}}(t)$ is also manifested in the speed $v_{\mathrm{KL}}(t)$. For
the same set of parameters that produced QMPE in KL divergence [Fig. 6(a)], we
investigate QMPE in the speed in Fig. 6(b). Indeed, we observe that the speeds
of the two copies I and II intersect each other, leading to QMPE. The higher
distant copy goes towards steady state with a higher speed and after the
intersection its speed becomes lower. Notably, the intersection times for
$D_{\mathrm{KL}}(t)$ and $v_{\mathrm{KL}}(t)$ are different. In fact, there is
a time window where the initially higher distant copy becomes lower distant
copy but still maintains its higher speed for some time. It would be
intriguing to investigate whether there is some general upper bound on the
relaxation speed $v_{\mathrm{KL}}(t)$.

We end this section by mentioning that the analysis of the QMPE in region (c)
would follow similar manner as in region (d) discussed elaborately in the
present section. This is because region (c) is also a line of second-order EP
just like region (d), the only difference is that the arrangement of
eigenvalues for region (c) obeys $0<\lambda_{2}<\lambda_{3}=\lambda_{4}$.

##  V QMPE in region ($\mathrm{a_{1}}$): oscillations

In this section, we explore QMPE in region ($\mathrm{a_{1}}$) of Hatano’s
model which contains complex eigenvalues with both non-zero real and imaginary
parts. The eigenvalues are non-degenerate and they are denoted by
$0,\lambda_{2},\lambda_{3}=\lambda_{\mathrm{re}}+\mathrm{i}\lambda_{\mathrm{im}},\lambda_{4}=\lambda_{\mathrm{re}}-\mathrm{i}\lambda_{\mathrm{im}}$,
where the subscripts $\mathrm{re}$ and $\mathrm{im}$ correspond to real and
imaginary parts, respectively. In addition, the eigenvalues satisfy the
relation $0<\lambda_{2}<\lambda_{\mathrm{re}}$ meaning $\lambda_{2}$ is the
slowest relaxation mode.

Figure 7: The figure shows multiple QMPE in the ground state density matrix
element, where both copies I and II exhibit oscillatory relaxation and
intersect each other multiple times. The inset shows the difference between
the two copies intersect zero multiple times, confirming the presence of
multiple QMPE. Parameters used are
$\tilde{d}=2.5,\tilde{\Gamma}_{\mathrm{I}}=\tilde{\Gamma}_{\mathrm{II}}=\tilde{\Gamma}=0.5,\tilde{d}_{\mathrm{I}}=2.1,\tilde{d}_{\mathrm{II}}=0.51$.

In contrast to the nontrivial approximate diagonalization [Eqs. (17)-(22)] for
second-order EP in Sec. IV, here the complete diagonalization can be obtained
through the usual method

$\displaystyle\widehat{\mathcal{L}}_{d}=\widehat{L}\widehat{\mathcal{L}}\widehat{R},$
(45)

where $\widehat{\mathcal{L}}_{d}$ is the diagonal form of the Lindbladian
$\widehat{\mathcal{L}}$ [Eq. (2)] in region ($\mathrm{a_{1}}$). The explicit
forms of the eigenvalues, right eigenvectors $|r_{k}\rangle$ and left
eigenvectors $\langle\ell_{k}|$ [($k=2,3,4$), which appear in Eq. (16)], are
given in Appendix B.

###  V.1 Density matrix elements

The the individual density matrix elements $\rho_{j}(t)$ ($j=1,2,3,4$) are
given by

$\displaystyle\rho_{j}(t)=\sum_{k=1}^{4}e^{-\lambda_{k}t}r_{k,j}a_{k},\hskip
5.69046pta_{k}=\sum_{n=1}^{4}\ell_{k,n}\rho_{n}(0).$ (46)

To discuss QMPE in the density matrix, let us focus on the ground state
probability i.e. $\rho_{\mathrm{gg}}(t)$. Of course, there is no algebraic
time dependence in the relaxation observed in Sec. IV. Rather, the interesting
part of the relaxation in region $(\mathrm{a_{1}})$ originates from the
complex nature of the eigenvalues and the corresponding coefficients $a_{k}$
in Eq. (46). Specifically, we have found that $a_{3}$ and $a_{4}$ are complex
conjugates of each other, they take the forms
$a_{3}=a_{\mathrm{re}}+\mathrm{i}a_{\mathrm{im}}$ and
$a_{4}=a_{\mathrm{re}}-\mathrm{i}a_{\mathrm{im}}$. Using such expressions, we
obtain the following formula for $\rho_{\mathrm{gg}}(t)$:

$\displaystyle\rho_{\mathrm{gg}}(t)=$
$\displaystyle\frac{4+\tilde{d}^{2}+\tilde{\Gamma}^{2}}{4+2\tilde{d}^{2}+\tilde{\Gamma}^{2}}+e^{-\lambda_{2}t}a_{2}$
(48)
$\displaystyle+2e^{-\lambda_{\mathrm{re}}t}\sqrt{a^{2}_{\mathrm{re}}+a^{2}_{\mathrm{im}}}\,\mathrm{cos}(\lambda_{\mathrm{im}}t-\phi),$

where $\phi=\mathrm{tan}^{-1}(a_{\mathrm{im}}/a_{\mathrm{re}})$.
Interestingly, the Eq. (48) contains an oscillatory term
$\mathrm{cos}(\lambda_{\mathrm{im}}t-\phi)$ and such an oscillatory relaxation
differs from the usual exponential relaxation studied in classical MPE Kumar20
; Klich19 . To investigate the effect of such oscillations in the QMPE, we
consider relaxations from two different initial conditions I and II, and
search for one (single QMPE) or more (multiple QMPE) intersection between the
time evolved trajectories. The detection of such intersections is made from
the solutions of the equation
$\Delta\rho_{\mathrm{gg}}(t)=\rho^{\mathrm{I}}_{\mathrm{gg}}(t)-\rho^{\mathrm{II}}_{\mathrm{gg}}(t)=0$.
In region ($\mathrm{a_{1}}$), the difference $\Delta\rho_{\mathrm{gg}}(t)$ has
the form

$\displaystyle\Delta$
$\displaystyle\rho_{\mathrm{gg}}(t)=2e^{-\lambda_{2}t}\sqrt{(\Delta
a_{\mathrm{re}})^{2}+(\Delta a_{\mathrm{im}})^{2}}\times$
$\displaystyle\left[\frac{\Delta a_{2}}{2\sqrt{(\Delta
a_{\mathrm{re}})^{2}+(\Delta
a_{\mathrm{im}})^{2}}}+e^{-(\lambda_{\mathrm{re}}-\lambda_{2})t}\mathrm{cos}(\lambda_{\mathrm{im}}t-\theta)\right],$

where $\theta=\mathrm{tan}^{-1}(\Delta a_{\mathrm{im}}/\Delta
a_{\mathrm{re}})$ with $\Delta
a_{\mathrm{re}}=a^{\mathrm{I}}_{\mathrm{re}}-a^{\mathrm{II}}_{\mathrm{re}}$
and $\Delta
a_{\mathrm{im}}=a^{\mathrm{I}}_{\mathrm{im}}-a^{\mathrm{II}}_{\mathrm{im}}$.
Due to the presence of the oscillatory cos(.) function in Eq. (LABEL:eq:a15),
it is possible that $\Delta\rho_{\mathrm{gg}}(t)$ can become zero at several
points of time leading to multiple QMPE. A necessary and sufficient condition
for QMPE in region ($\mathrm{a_{1}}$) is

$\displaystyle\mathrm{min}[e^{-(\lambda_{\mathrm{re}}-\lambda_{2})t}$
$\displaystyle\mathrm{cos}(\lambda_{\mathrm{im}}t-\theta)]<\frac{\Delta
a_{2}}{2\sqrt{(\Delta a_{\mathrm{re}})^{2}+(\Delta a_{\mathrm{im}})^{2}}}$
(52)
$\displaystyle<\mathrm{max}[e^{-(\lambda_{\mathrm{re}}-\lambda_{2})t}\mathrm{cos}(\lambda_{\mathrm{im}}t-\theta)],$

where $\Delta a_{2}$ and $\mathrm{cos}(\lambda_{\mathrm{im}}t-\theta)$ must
have opposite signs. Here $\mathrm{min}[f]$ and $\mathrm{max}[f]$ correspond
to minimum and maximum values of the function $f$. Note that
$-1<\mathrm{min}[e^{-(\lambda_{\mathrm{re}}-\lambda_{2})t}\mathrm{cos}(\lambda_{\mathrm{im}}t-\theta)]$
and
$\mathrm{max}[e^{-(\lambda_{\mathrm{re}}-\lambda_{2})t}\mathrm{cos}(\lambda_{\mathrm{im}}t-\theta)]<1$.

In Fig. 7 we present an example of multiple QMPE in the ground state density
matrix element $\rho_{\mathrm{gg}}(t)$, where we have chosen parameters that
satisfy the QMPE criteria Eq. (52) for region ($\mathrm{a_{1}}$). We observe
that the two copies I and II relax in oscillatory manner and intersect each
other multiple times. The inset explicitly shows that the difference between
the two copies $\Delta\rho_{\mathrm{gg}}(t)$ intersects zero multiple times,
confirming the multiple QMPE. Comparing Fig. 3(a) [region (d)] with the inset
of Fig. 7 [region ($\mathrm{a_{1}}$)], we understand that the number of
intersections in region ($\mathrm{a_{1}}$) is not restricted to two as in
region (d). Therefore, in region ($\mathrm{a_{1}}$) one can have many more
intersections and higher order QMPE compared to the maximum double QMPE in the
region (d). Since the amplitudes of the oscillations can be controlled by
tuning the initial driving strengths $\tilde{d}_{\mathrm{I}}$ and
$\tilde{d}_{\mathrm{II}}$, complex eigenvalues seem to be useful to obtain
multiple intersections and thereby multiple QMPE. It would be important in the
future to find the exact number of solutions to
$\Delta\rho_{\mathrm{gg}}(t)=0$ either analytically or numerically.

###  V.2 Energy

The formal expression for average energy $E(t)$ is given by Eq. (34). In the
previous Sec. IV, we have already shown that QMPE in ground state probability
does not necessarily imply QMPE in energy, in general. However, it is possible
to find control parameters that can give rise to QMPE in both observables,
although the intersection time(s) are generally different. Following the same
path, we would like to see if $E(t)$ shows multiple QMPE for the same set of
parameters used to achieve multiple QMPE for $\rho_{\mathrm{gg}}(t)$ in Fig.
7. Indeed, Fig. 8 exhibits multiple intersections between time evolved copies
of energies starting from I and II, both relaxing towards the same final
energy. The inset of Fig. 8 alternatively demonstrates the multiple QMPE
through multiple intersections of the energy difference $\Delta
E(t)=E_{\mathrm{I}}(t)-E_{\mathrm{II}}(t)$ with zero, remarkably different
from Fig. 3(b) with a maximum of two such intersections in region (d).

Figure 8: The figure illustrates multiple QMPE in average energy where two
copies I and II relaxes with oscillations and intersect each other multiple
times showing multiple QMPE. The inset shows the energy difference between the
two copies intersect zero multiple times, leading to multiple QMPE in energy.
Parameters used are same as in Fig. 7.

###  V.3 von Neumann entropy

Next we consider the von Neumann entropy $S_{\mathrm{vN}}(t)$ [Eq. (41)]. To
observe QMPE in $S_{\mathrm{vN}}(t)$, we consider the intersection times for
the entropy difference between two copies $\Delta
S_{\mathrm{vN}}(t)=S^{\mathrm{I}}_{\mathrm{vN}}-S^{\mathrm{II}}_{\mathrm{vN}}$.
In Fig. 9, we present the behavior of $\Delta S_{\mathrm{vN}}(t)$ for the same
set of parameter values used in Figs. 7 and 8. Indeed, we observe that $\Delta
S_{\mathrm{vN}}(t)$ intersects zero multiple times showing the existence of
multiple QMPE in the von Neumann entropy.

Figure 9: The figure shows the temporal behavior of the difference between von
Neumann entropies of two copies I and II, in region ($\mathrm{a_{1}}$). The
entity $\Delta S_{\mathrm{vN}}(t)$ intersect zero multiple times indicating
the occurrence of multiple QMPE. Parameters used are same as in Figs. 7.

###  V.4 Temperature

Here, we study the thermal QMPE in region ($\mathrm{a_{1}}$) using the system
temperature defined in Eq. (42). In Fig. 10, we present the time evolution of
the temperatures of two copies I and II, one initially hotter than the other.
It is fascinating that the temperatures of both the copies relax in
oscillatory manner towards the same steady state temperature and during the
relaxations, they intersect each other multiple times, leading to multiple
thermal QMPE. This, in turn, implies that the hotter and colder copies reverse
their identities multiple times (after each intersection the previous hotter
copy turns into colder copy and remains so until the next intersection). Such
oscillatory multiple thermal QMPE is observed here for the first time in
quantum systems (as far as our knowledge is concerned), although multiple
thermal MPE in classical inertial suspensions has already been reported
Takada21a . The physical interpretation and utility of multiple thermal QMPE
have to be explored in future in more details. Note that we have used the same
parameter set for multiple thermal QMPE as in Figs. 7-9. Similar to Ref.
Chatterjee_2023 , here also we observe negative temperature Landau80 ;
Gauthier19 ; Johnstone19 ; Medley11 ; Ali20  in the early stage of relaxation.
This is due to non-monotonic natures of energy and entropy as functions of
time. However, it is satisfying that even if we neglect such initial negative
temperatures, we observe multiple QMPE at later times of the relaxation
process where the temperatures of both the copies remain positive.

Figure 10: The figure shows multiple thermal QMPE. Two copies, one initially
hotter than the other, relax in oscillatory manner and intersect each other
multiple times. After each intersection, the hotter and colder copies reverse
their identities and these identities are preserved only until the next
intersection. Parameters used are same as in Fig. 7.

###  V.5 Kullback-Leibler divergence

In this subsection we explore the QMPE in KL divergence defined in Eq. (43).
Since we have been able to obtain multiple QMPE in region ($\mathrm{a_{1}}$)
in all other observables considered here, it is natural to ask about the
possibility of multiple QMPE in $D_{\mathrm{KL}}(t)$ as well.

Figure 11: The figure (a) shows single intersection between the KL divergence
of the two copies I and II and thereby exhibit single QMPE instead of multiple
QMPE. Parameters used are same as in Fig. 7 only except
$\tilde{d}_{\mathrm{I}}=1.0$. For these same parameters, the figure (b)
compares the time evolution of the speeds of the two copies. Unlike the
distance, the speeds intersect each other multiple times leading to multiple
QMPE.

However, we have found that the KL divergence does not show any intersection
between two copies I and II for the same parameter values used in Figs. 7-10.
Therefore, we have to use different values of the initial drives (still within
the fixed dissipation protocol) to achieve QMPE in $D_{\mathrm{KL}}(t)$,
presented in Fig. 11(a). Intriguingly, we have only observed single
intersection between the two time evolved copies, resulting in single QMPE in
KL divergence. This means that once the initially higher distant copy
intersects the initially lower distant copy producing QMPE, thereafter the new
lower distant copy (that was initially higher distant) remains closer to the
steady state for the rest of the relaxation process. In fact, for all other
parameter sets that we have checked, we could not find multiple QMPE in
$D_{\mathrm{KL}}(t)$. The proof of the possibility or impossibility of
multiple QMPE in KL divergence in presence of complex eigenvalues remains
unsolved and could be an important question to investigate later.

In this scenario, it would be interesting to examine the corresponding
behaviors of the speeds $v_{\mathrm{KL}}(t)$ [Eq. (44)] of the two copies I
and II. In Fig. 11(b), we present the time evolution of speeds for the same
set of parameters that gave rise to single QMPE in $D_{\mathrm{KL}}(t)$. It is
fascinating to observe that the speeds, unlike the distances, intersect each
other multiple times producing multiple QMPE. The first intersection between
the speeds of the two copies [Fig. 11(b)] shows a time lag in comparison to
the intersection between distances of the same copies [Fig. 11(a)]. Such a
time lag between the QMPE in KL divergence and the speed is in agreement with
the observation in region (d) [Fig. 6]. However, contrary to single QMPE in
$v_{\mathrm{KL}}(t)$ in region (d), we have obtained multiple QMPE in the
speed in region ($\mathrm{a_{1}}$). The contrasting behavior of single QMPE in
$D_{\mathrm{KL}}(t)$ and multiple QMPE in $v_{\mathrm{KL}}(t)$ imply that the
initially higher distant copy that becomes lower distant after the
intersection in distances, can remain at lower distance even if its speed
sometimes become smaller than the speed of the current higher distant copy.
This is interesting because the system almost behaves like an intelligent and
efficient runner in a race who knows how to adjust the speed (i.e. sometimes
running faster and sometimes running slower than the people behind) and still
remains in the first position to finish the race.

##  VI Summary

We have studied QMPE in a two-level quantum system subjected to an oscillatory
electric field and dissipative coupling with the environment. The Lindbladian
of this two-level driven dissipative system, termed as Hatano’s model, has
rich variety of eigen-spectrum, especially EPs and complex eigenvalues which
are generic features of non-Hermitian quantum systems. We have investigated
the roles of complex eigenvalues and EPs on the QMPE. The exact analytical
formulae obtained for time evolved density matrix elements enable us to
propose analytical criteria for the QMPE in observables. Interestingly, at the
EPs, the density matrix elements and consequently observables (e.g. energy)
attain additional algebraic time dependent relaxations apart from the usual
exponential relaxations. Such additional algebraic time dependence at the EPs
leads to double QMPE (twice intersections between trajectories starting at
different initial conditions and relaxing towards the same steady state). We
have procured the exact expressions for intersection times of trajectories,
that mark the occurrence of QMPE and we provide analytical criteria that
demarcate control parameters’ regimes with single or double or no QMPE.
Whereas the multiple QMPE at EPs in the present model is restricted to maximum
double QMPE i.e. two intersections, the multiple QMPE at the complex
eigenvalues can result in even more number of intersections in a oscillatory
manner. Notably, we have established the existence of thermal QMPE (single
QMPE at EPs and multiple QMPE at complex eigenvalues), in presence of
coherence in Hatano’s model. It is fascinating that, similar to energy and
entropy, temperature also exhibits multiple QMPE where the initially hotter
and colder copies intersect each other multiple times. Such multiple reversals
of identities (i.e. hotter becoming colder and vice versa) result in multiple
thermal QMPE. In spite of achieving multiple QMPE in several observables
(ground state probability, energy, entropy, temperature) using either
additional algebraic time dependence at EPs or oscillations at complex
eigenvalues, surprisingly the KL divergence (a distance measure from steady
state) shows single QMPE only. Overall, our theoretical analysis exhibits that
EPs and complex eigenvalues can be used as systematic routes to generate
multiple QMPE. However, we should mention that the presence of EPs or
oscillations are not necessary for multiple QMPE. In this connection, we have
shown that multiple QMPE can be achieved even with purely exponential
relaxations, by putting suitable constraints on the gap statistics of the
eigen-spectrum (Appendix C).

An important and generic future goal in QMPE research would be to understand
its utility in speeding up quantum processes. In fact, a recent work on
accelerating relaxation in open quantum systems by engineering EPs Zhou23
hints towards the connection of such acceleration to QMPE. Our present work
analyzes in detail how EPs lead to QMPE and provide criteria on control
parameters to observe the QMPE. It is important to find out the deeper
connections between the QMPE and faster relaxations at EPs and for generic
open quantum systems. In this connection, an intriguing question is to find
how the quantum speed limit may affect the degree of faster relaxation that
can be achieved through the QMPE. Our work on QMPE by controlling the
amplitudes of oscillatory electric fields hints towards possible experimental
realizations of the QMPE in qubit, keeping in mind the recent advancement in
qubit related experiment Gyenis21 . The analysis of QMPE in the present study
paves path for future investigations of QMPE in quantum many-body systems.

Acknowledgements.- We thank Naomichi Hatano, Sosuke Ito, Ryo Hanai, Arnab Pal,
Naruo Ohga and Kohei Yoshimura for useful discussions. This work is partially
supported by the Grants-in-Aid for Scientific Research (Grant No. 21H01006 and
No. 20K14428). A.K.C. gratefully acknowledges postdoctoral fellowship from the
YITP. The numerical calculations have been done on Yukawa-21 at the YITP.

##  Appendix A Eigenvectors of $\mathcal{L}$ in region (d)

In this Appendix, we give the explicit expressions for the right eigenvectors
$|r_{k}\rangle$ and left eigenvectors $\langle\ell_{k}|$ of the Lindbladian
[Eq. (2)] in region (d) [Sec. IV], where $k=2,3,4$. The right and left
eigenvectors corresponding to the zero eigenvalue i.e. the steady state. have
already been provided in Eq. (6). Below are the expressions for the right
eigenvectors

$\displaystyle|r_{2}\rangle$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}-\displaystyle\frac{\tilde{d}}{1-\mathrm{i}(\tilde{\Gamma}/2-\lambda_{2})},&-\displaystyle\frac{\tilde{d}}{1+\mathrm{i}(\tilde{\Gamma}/2-\lambda_{2})},&-1,&1\end{array}\right)^{T},$
(2) $\displaystyle|r_{3}\rangle$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}\displaystyle\frac{1-\mathrm{i}(\tilde{\Gamma}/2-\lambda_{4})}{\tilde{d}},&\displaystyle\frac{-1-\mathrm{i}(\tilde{\Gamma}/2-\lambda_{4})}{\tilde{d}},&0,&0\end{array}\right)^{T},$
(4) $\displaystyle|r_{4}\rangle$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}-\displaystyle\frac{\tilde{d}}{1-\mathrm{i}(\tilde{\Gamma}/2-\lambda_{4})},&-\displaystyle\frac{\tilde{d}}{1+\mathrm{i}(\tilde{\Gamma}/2-\lambda_{4})},&-1,&1\end{array}\right),$
(6)

where the superscript $T$ denotes transpose. The left eigenvectors are
obtained to be

$\displaystyle\langle\ell_{2}|=\left(\begin{array}[]{cccc}-\ell^{\mathrm{re}}_{21}-\mathrm{i}\ell^{\mathrm{im}}_{21},&-\ell^{\mathrm{re}}_{21}+\mathrm{i}\ell^{\mathrm{im}}_{21},&\ell_{23},&\ell_{24}\end{array}\right),$
(8)
$\displaystyle\langle\ell_{3}|=\left(\begin{array}[]{cccc}\ell^{\mathrm{re}}_{31}+\mathrm{i}\ell^{\mathrm{im}}_{31},&-\ell^{\mathrm{re}}_{31}+\mathrm{i}\ell^{\mathrm{im}}_{31},&\ell_{33},&\ell_{34}\end{array}\right),$
(10)
$\displaystyle\langle\ell_{4}|=\left(\begin{array}[]{cccc}\ell^{\mathrm{re}}_{21}+\mathrm{i}\ell^{\mathrm{im}}_{21},&\ell^{\mathrm{re}}_{21}-\mathrm{i}\ell^{\mathrm{im}}_{21},&\ell_{43},&\ell_{44}\end{array}\right),\hskip
8.5359pt$ (12)

where the elements $l_{kn}$ are explicitly given by

$\displaystyle\ell^{\mathrm{re}}_{21}$ $\displaystyle=$
$\displaystyle-\frac{(4+(\tilde{\Gamma}-2\lambda_{2})^{2})(4+(\tilde{\Gamma}-2\lambda_{4})^{2})}{4\tilde{d}(\lambda_{2}-\lambda_{4})\eta},$
(13) $\displaystyle\ell^{\mathrm{im}}_{21}$ $\displaystyle=$
$\displaystyle-\frac{(4+(\tilde{\Gamma}-2\lambda_{2})^{2})(\tilde{\Gamma}-2\lambda_{4})(4+(\tilde{\Gamma}-2\lambda_{4})^{2})}{8\tilde{d}(\lambda_{2}-\lambda_{4})\eta},$
(14) $\displaystyle\ell_{23}$ $\displaystyle=$
$\displaystyle\frac{(4+(\tilde{\Gamma}-2\lambda_{2})^{2})\left(\delta_{1}-2\tilde{\Gamma}\lambda_{4}(20+\tilde{\Gamma}^{2}-4\lambda^{2}_{4})-4\lambda^{2}_{4}(\tilde{\Gamma}^{2}-12)\right)}{4(4+2\tilde{d}^{2}+\tilde{\Gamma}^{2})(\lambda_{2}-\lambda_{4})\eta},$
(15) $\displaystyle\ell_{24}$ $\displaystyle=$
$\displaystyle-\frac{(4+(\tilde{\Gamma}-2\lambda_{2})^{2})\left(\delta_{1}+2\tilde{\Gamma}\lambda_{4}(4-3\tilde{\Gamma}^{2}-4\lambda^{2}_{4})-4\lambda^{2}_{4}(4-3\tilde{\Gamma}^{2})\right)}{4(4+2\tilde{d}^{2}+\tilde{\Gamma}^{2})(\lambda_{2}-\lambda_{4})\eta},$
(16) $\displaystyle\ell^{\mathrm{re}}_{31}$ $\displaystyle=$ $\displaystyle
4\tilde{d}(\lambda_{2}+\lambda_{4}-\tilde{\Gamma})\eta^{-1},$ (17)
$\displaystyle\ell^{\mathrm{re}}_{32}$ $\displaystyle=$
$\displaystyle\tilde{d}((\tilde{\Gamma}-2\lambda_{2})(\tilde{\Gamma}-2\lambda_{4})-4)\eta^{-1},$
(18) $\displaystyle\ell_{33}$ $\displaystyle=$
$\displaystyle\frac{4\mathrm{i}\tilde{d}^{2}(4+2\tilde{d}^{2}+\tilde{\Gamma}^{2}+4\lambda_{2}\lambda_{4})}{(4+2\tilde{d}^{2}+\tilde{\Gamma}^{2})\eta},$
(19) $\displaystyle\ell_{34}$ $\displaystyle=$
$\displaystyle-\frac{4\mathrm{i}\tilde{d}^{2}(4+2\tilde{d}^{2}+\tilde{\Gamma}^{2}-4\lambda_{2}\lambda_{4})}{(4+2\tilde{d}^{2}+\tilde{\Gamma}^{2})\eta},$
(20) $\displaystyle\ell_{43}$ $\displaystyle=$
$\displaystyle-\frac{(4+(\tilde{\Gamma}-2\lambda_{4})^{2})\left(\delta_{2}-\delta_{3}-32\lambda_{2}(\tilde{\Gamma}-\lambda_{4})\right)}{4(\lambda_{2}-\lambda_{4})(4+2\tilde{d}^{2}+\tilde{\Gamma}^{2})\eta},$
(21) $\displaystyle\ell_{44}$ $\displaystyle=$
$\displaystyle\frac{(4+(\tilde{\Gamma}-2\lambda_{4})^{2})\left(\delta_{2}+\delta_{3}+4\tilde{\Gamma}\lambda_{2}(4-\tilde{\Gamma}^{2}+2\tilde{\Gamma}\lambda_{4})\right)}{4(\lambda_{2}-\lambda_{4})(4+2\tilde{d}^{2}+\tilde{\Gamma}^{2})\eta},$
(22)

where the eigenvalues $\lambda_{2}$ and $\lambda_{4}$ are given by Eq. (24)
and the parameters $\eta,\delta_{1},\delta_{2},\delta_{3}$ used in Eq. (22)
have the following forms

$\displaystyle\eta$ $\displaystyle=$
$\displaystyle\left(\tilde{\Gamma}^{3}-2\tilde{\Gamma}^{2}(\lambda_{2}+2\lambda_{4})+4\tilde{\Gamma}(-3+2\lambda_{2}\lambda_{4}+\lambda^{2}_{4})+8(\lambda_{2}+2\lambda_{4}-\lambda_{2}\lambda^{2}_{4})\right),$
(23) $\displaystyle\delta_{1}$ $\displaystyle=$
$\displaystyle-16+\tilde{\Gamma}^{4}+2\tilde{d}^{2}((\tilde{\Gamma}-2\lambda_{4})^{2}-4),$
(24) $\displaystyle\delta_{2}$ $\displaystyle=$
$\displaystyle-16+\tilde{\Gamma}^{4}-8\tilde{\Gamma}\lambda_{4}-2\tilde{\Gamma}^{3}\lambda_{4}+2\tilde{d}^{2}\left((\tilde{\Gamma}-2\lambda_{2})(\tilde{\Gamma}-2\lambda_{4})-4\right),$
(25) $\displaystyle\delta_{3}$ $\displaystyle=$ $\displaystyle
4\lambda^{2}_{2}\left(\tilde{\Gamma}^{2}-2\tilde{\Gamma}\lambda_{4}-4\right).$
(26)

##  Appendix B Eigenvectors of $\mathcal{L}$ in region ($\mathrm{a_{1}}$)

In this Appendix, we give the explicit expressions for the right eigenvectors
$|r_{k}\rangle$ and left eigenvectors $\langle\ell_{k}|$ of the Lindbladian
[Eq. (2)] in region ($\mathrm{a_{1}}$) [Sec. V], where $k=2,3,4$. The right
and left eigenvectors corresponding to the zero eigenvalue i.e. steady state.
have already been provided in Eq. (6). Below are the expressions for the right
eigenvectors

$\displaystyle|r_{2}\rangle$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}-\displaystyle\frac{\tilde{d}}{1-\mathrm{i}\left(\frac{\tilde{\Gamma}}{2}-\lambda_{2}\right)},&-\displaystyle\frac{\tilde{d}}{1+\mathrm{i}\left(\frac{\tilde{\Gamma}}{2}-\lambda_{2}\right)},&-1,&1\end{array}\right)^{T},$
(28) $\displaystyle|r_{3}\rangle$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}-\displaystyle\frac{\tilde{d}}{1-\lambda_{\mathrm{im}}-\mathrm{i}\left(\frac{\tilde{\Gamma}}{2}-\lambda_{\mathrm{re}}\right)},&-\displaystyle\frac{\tilde{d}}{1+\lambda_{\mathrm{im}}+\mathrm{i}\left(\frac{\tilde{\Gamma}}{2}-\lambda_{\mathrm{re}}\right)},&-1,&1\end{array}\right)^{T},$
(30) $\displaystyle|r_{4}\rangle$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}-\displaystyle\frac{\tilde{d}}{1+\lambda_{\mathrm{im}}-\mathrm{i}\left(\frac{\tilde{\Gamma}}{2}-\lambda_{\mathrm{re}}\right)},&-\displaystyle\frac{\tilde{d}}{1-\lambda_{\mathrm{im}}+\mathrm{i}\left(\frac{\tilde{\Gamma}}{2}-\lambda_{\mathrm{re}}\right)},&-1,&1\end{array}\right)^{T}.$
(32)

The left eigenvectors are given by

$\displaystyle\langle\ell_{2}|=\left(\begin{array}[]{cccc}\ell^{\mathrm{re}}_{21}+\mathrm{i}\ell^{\mathrm{im}}_{21},&\ell^{\mathrm{re}}_{21}-\mathrm{i}\ell^{\mathrm{im}}_{21},&\ell_{23},&\ell_{24}\end{array}\right),$
(34)
$\displaystyle\langle\ell_{3}|=\left(\begin{array}[]{cccc}\ell^{\mathrm{re}}_{31}+\mathrm{i}\ell^{\mathrm{im}}_{31},&\ell^{\mathrm{re}}_{32}+\mathrm{i}\ell^{\mathrm{im}}_{32},&\ell^{\mathrm{re}}_{33}+\mathrm{i}\ell^{\mathrm{im}}_{33},&\ell^{\mathrm{re}}_{34}+\mathrm{i}\ell^{\mathrm{im}}_{34}\end{array}\right),$
(36)
$\displaystyle\langle\ell_{4}|=\left(\begin{array}[]{cccc}\ell^{\mathrm{re}}_{32}-\mathrm{i}\ell^{\mathrm{im}}_{32},&\ell^{\mathrm{re}}_{31}-\mathrm{i}\ell^{\mathrm{im}}_{31},&\ell^{\mathrm{re}}_{33}-\mathrm{i}\ell^{\mathrm{im}}_{33},&\ell^{\mathrm{re}}_{34}-\mathrm{i}\ell^{\mathrm{im}}_{34}\end{array}\right),\hskip
8.5359pt$ (38)

where

$\displaystyle\ell^{\mathrm{re}}_{21}$ $\displaystyle=$
$\displaystyle-\frac{(4+(\tilde{\Gamma}-2\lambda_{2})^{2})(4(\lambda^{2}_{\mathrm{im}}-1)+(\tilde{\Gamma}-2\lambda_{\mathrm{re}})^{2})}{32\tilde{d}\left(\lambda^{2}_{\mathrm{im}}+(\lambda_{2}-\lambda_{\mathrm{re}})^{2}\right)},$
(39) $\displaystyle\ell^{\mathrm{im}}_{21}$ $\displaystyle=$
$\displaystyle-\frac{(4+(\tilde{\Gamma}-2\lambda_{2})^{2})(\tilde{\Gamma}-2\lambda_{\mathrm{re}})}{8\tilde{d}\left(\lambda^{2}_{\mathrm{im}}+(\lambda_{2}-\lambda_{\mathrm{re}})^{2}\right)},$
(40) $\displaystyle\ell_{23}$ $\displaystyle=$
$\displaystyle-\frac{(4+(\tilde{\Gamma}-2\lambda_{2})^{2})(2\tilde{d}^{2}+\tilde{\Gamma}^{2}+4(1+\lambda^{2}_{\mathrm{im}}+\lambda^{2}_{\mathrm{re}}))}{8(4+2\tilde{d}^{2}+\tilde{\Gamma}^{2})\left(\lambda^{2}_{\mathrm{im}}+(\lambda_{2}-\lambda_{\mathrm{re}})^{2}\right)},$
(41) $\displaystyle\ell_{24}$ $\displaystyle=$
$\displaystyle\frac{(4+(\tilde{\Gamma}-2\lambda_{2})^{2})(2\tilde{d}^{2}+\tilde{\Gamma}^{2}-4(1+\lambda^{2}_{\mathrm{im}}+\lambda^{2}_{\mathrm{re}}))}{8(4+2\tilde{d}^{2}+\tilde{\Gamma}^{2})\left(\lambda^{2}_{\mathrm{im}}+(\lambda_{2}-\lambda_{\mathrm{re}})^{2}\right)},$
(42) $\displaystyle\ell^{\mathrm{re}}_{31}$ $\displaystyle=$
$\displaystyle\frac{(\xi_{1}+\xi_{2})(4(1+\lambda_{\mathrm{im}})^{2}+(\tilde{\Gamma}-2\lambda_{\mathrm{re}})^{2})}{64\tilde{d}\lambda_{\mathrm{im}}\left(\lambda^{2}_{\mathrm{im}}+(\lambda_{2}-\lambda_{\mathrm{re}})^{2}\right)},$
(43) $\displaystyle\ell^{\mathrm{re}}_{32}$ $\displaystyle=$
$\displaystyle\frac{(\xi_{1}-\xi_{2})(4(\lambda_{\mathrm{im}}-1)^{2}+(\tilde{\Gamma}-2\lambda_{\mathrm{re}})^{2})}{64\tilde{d}\lambda_{\mathrm{im}}\left(\lambda^{2}_{\mathrm{im}}+(\lambda_{2}-\lambda_{\mathrm{re}})^{2}\right)},$
(44) $\displaystyle\ell^{\mathrm{im}}_{31}$ $\displaystyle=$
$\displaystyle\frac{\left(\xi_{3}-2\tilde{\Gamma}\lambda_{\mathrm{im}}(\lambda_{\mathrm{im}}-2)+4\lambda_{2}(\lambda_{\mathrm{im}}-1)^{2}\right)(4(1+\lambda_{\mathrm{im}})^{2}+(\tilde{\Gamma}-2\lambda_{\mathrm{re}})^{2})}{64\tilde{d}\lambda_{\mathrm{im}}\left(\lambda^{2}_{\mathrm{im}}+(\lambda_{2}-\lambda_{\mathrm{re}})^{2}\right)},$
(45) $\displaystyle\ell^{\mathrm{im}}_{32}$ $\displaystyle=$
$\displaystyle\frac{\left(\xi_{3}-2\tilde{\Gamma}\lambda_{\mathrm{im}}(\lambda_{\mathrm{im}}+2)+4\lambda_{2}(\lambda_{\mathrm{im}}+1)^{2}\right)(4(-1+\lambda_{\mathrm{im}})^{2}+(\tilde{\Gamma}-2\lambda_{\mathrm{re}})^{2})}{64\tilde{d}\lambda_{\mathrm{im}}\left(\lambda^{2}_{\mathrm{im}}+(\lambda_{2}-\lambda_{\mathrm{re}})^{2}\right)},$
(46)

where the parameters $\xi_{1},\xi_{2},\xi_{3}$ used in Eq. (46) have the
following forms

$\displaystyle\xi_{1}$ $\displaystyle=$
$\displaystyle\tilde{\Gamma}^{2}\lambda_{\mathrm{im}}-4\tilde{\Gamma}\lambda_{2}\lambda_{\mathrm{im}}+4\lambda_{\mathrm{im}}(\lambda^{2}_{2}-1),$
(47) $\displaystyle\xi_{2}$ $\displaystyle=$ $\displaystyle
4\tilde{\Gamma}(\lambda_{2}-\lambda_{\mathrm{re}})+4(\lambda^{2}_{\mathrm{re}}-\lambda^{2}_{2}+\lambda^{2}_{\mathrm{im}}),$
(48) $\displaystyle\xi_{3}$ $\displaystyle=$
$\displaystyle-\tilde{\Gamma}^{2}\lambda_{2}+2\lambda^{2}_{2}(\tilde{\Gamma}-2\lambda_{\mathrm{re}})-4\lambda_{\mathrm{re}}+\tilde{\Gamma}\lambda_{\mathrm{re}}(\tilde{\Gamma}-2\lambda_{\mathrm{re}})+4\lambda_{2}\lambda^{2}_{\mathrm{re}}.$
(49)

The four other elements
$\ell^{\mathrm{re}}_{33},\ell^{\mathrm{im}}_{33},\ell^{\mathrm{re}}_{34}$ and
$\ell^{\mathrm{im}}_{34}$ in Eq. (38) have too lengthy explicit expressions.
Therefore, we choose to provide their forms in terms of the elements whose
explicit expressions are already given in Eq. (46), as

$\displaystyle\ell^{\mathrm{re}}_{33}$ $\displaystyle=$
$\displaystyle\frac{\tilde{d}\left((\ell^{\mathrm{re}}_{31}+\ell^{\mathrm{re}}_{32})(8\lambda^{2}_{2}-8\tilde{\Gamma}\lambda_{2}-2(4+2\tilde{\Gamma}^{2}+\tilde{d}^{2}))+(\ell^{\mathrm{im}}_{32}-\ell^{\mathrm{im}}_{31})(4\tilde{\Gamma}\lambda^{2}_{2}+4\lambda_{2}(4+\tilde{d}^{2})-\tilde{\Gamma}(4+2\tilde{\Gamma}^{2}+\tilde{d}^{2}))\right)}{(4+2\tilde{\Gamma}^{2}+\tilde{d}^{2})(4+(\tilde{\Gamma}-2\lambda_{2})^{2})},$
(50) $\displaystyle\ell^{\mathrm{im}}_{33}$ $\displaystyle=$
$\displaystyle\frac{\tilde{d}\left((\ell^{\mathrm{im}}_{31}+\ell^{\mathrm{im}}_{32})(8\lambda^{2}_{2}-8\tilde{\Gamma}\lambda_{2}-2(4+2\tilde{\Gamma}^{2}+\tilde{d}^{2}))+(\ell^{\mathrm{re}}_{31}-\ell^{\mathrm{re}}_{32})(4\tilde{\Gamma}\lambda^{2}_{2}+4\lambda_{2}(4+\tilde{d}^{2})-\tilde{\Gamma}(4+2\tilde{\Gamma}^{2}+\tilde{d}^{2}))\right)}{(4+2\tilde{\Gamma}^{2}+\tilde{d}^{2})(4+(\tilde{\Gamma}-2\lambda_{2})^{2})},$
(51) $\displaystyle\ell^{\mathrm{re}}_{34}$ $\displaystyle=$
$\displaystyle\frac{\tilde{d}\left((\ell^{\mathrm{re}}_{31}+\ell^{\mathrm{re}}_{32})(8\lambda^{2}_{2}-8\tilde{\Gamma}\lambda_{2}+2(4+2\tilde{\Gamma}^{2}+\tilde{d}^{2}))+(\ell^{\mathrm{im}}_{32}-\ell^{\mathrm{im}}_{31})(4\tilde{\Gamma}\lambda^{2}_{2}-4\lambda_{2}(\tilde{\Gamma}^{2}+\tilde{d}^{2})+\tilde{\Gamma}(4+2\tilde{\Gamma}^{2}+\tilde{d}^{2}))\right)}{(4+2\tilde{\Gamma}^{2}+\tilde{d}^{2})(4+(\tilde{\Gamma}-2\lambda_{2})^{2})},$
(52) $\displaystyle\ell^{\mathrm{im}}_{34}$ $\displaystyle=$
$\displaystyle\frac{\tilde{d}\left((\ell^{\mathrm{im}}_{31}+\ell^{\mathrm{im}}_{32})(8\lambda^{2}_{2}-8\tilde{\Gamma}\lambda_{2}+2(4+2\tilde{\Gamma}^{2}+\tilde{d}^{2}))+(\ell^{\mathrm{re}}_{31}-\ell^{\mathrm{re}}_{32})(4\tilde{\Gamma}\lambda^{2}_{2}-4\lambda_{2}(\tilde{\Gamma}^{2}+\tilde{d}^{2})+\tilde{\Gamma}(4+2\tilde{\Gamma}^{2}+\tilde{d}^{2}))\right)}{(4+2\tilde{\Gamma}^{2}+\tilde{d}^{2})(4+(\tilde{\Gamma}-2\lambda_{2})^{2})}.$
(53)

The explicit expressions for $\lambda_{2},\lambda_{\mathrm{re}}$ and
$\lambda_{\mathrm{im}}$ in terms of control parameters
$(\tilde{d},\tilde{\Gamma})$ are obtained as

$\displaystyle\lambda_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{12}\left((-\mathrm{i}+\sqrt{3})\chi^{1/3}+8\tilde{\Gamma}-\frac{(\mathrm{i}+\sqrt{3})(12+12\tilde{d}^{2}-\tilde{\Gamma}^{2})}{\chi^{1/3}}\right)$
(54) $\displaystyle\lambda_{\mathrm{re}}$ $\displaystyle=$
$\displaystyle\frac{12\tilde{d}^{2}(\mathrm{i}+\sqrt{3})+(\mathrm{i}-\sqrt{3})\chi^{2/3}+16\tilde{\Gamma}\chi^{1/3}-(\mathrm{i}+\sqrt{3})(\tilde{\Gamma}^{2}-12)}{24\chi^{1/3}},$
(55) $\displaystyle\lambda_{\mathrm{im}}$ $\displaystyle=$
$\displaystyle\frac{12\sqrt{3}\tilde{d}^{2}(\mathrm{i}+\sqrt{3})-\sqrt{3}(\mathrm{i}-\sqrt{3})\chi^{2/3}-\sqrt{3}(\mathrm{i}+\sqrt{3})(\tilde{\Gamma}^{2}-12)}{24\chi^{1/3}},$
(56)

where $\chi$ has the following form

$\chi=\left(\mathrm{i}\tilde{\Gamma}(36-18\tilde{d}^{2}+\tilde{\Gamma}^{2})+6\sqrt{3\tilde{\Gamma}^{2}(\tilde{d}^{4}+20\tilde{d}^{2}-8)-48(1+\tilde{d}^{2})^{3}-3\tilde{\Gamma}^{4}}\right).$
(57)

##  Appendix C QMPE in region ($\mathrm{b}$): purely exponential relaxation

In this Appendix we discuss QMPE in region (b) of Hatano’s model that contain
eigenvalues
$0,-\mathrm{i}\lambda_{2},-\mathrm{i}\lambda_{3},-\mathrm{i}\lambda_{4}$
obeying $0<\lambda_{2}<\lambda_{3}<\lambda_{4}$.

###  C.1 Density matrix elements

The individual density matrix elements $\rho_{j}(t)$ ($j=1,2,3,4$) are given
by

$\displaystyle\rho_{j}(t)=\sum_{k=1}^{4}e^{-\lambda_{k}t}r_{k,j}a_{k},\hskip
5.69046pta_{k}=\sum_{n=1}^{4}\ell_{k,n}\rho_{n}(0).$ (58)

The right eigenvectors $|r_{k}\rangle$ and left eigenvectors
$\langle\ell_{k}|$, are of the form below

$\displaystyle|r_{k}\rangle=\left(\begin{array}[]{cccc}-\frac{\tilde{d}}{1-\mathrm{i}\left(\frac{\tilde{\Gamma}}{2}-\lambda_{k}\right)},&-\frac{\tilde{d}}{1+\mathrm{i}\left(\frac{\tilde{\Gamma}}{2}-\lambda_{k}\right)},&-1,&1\end{array}\right)^{T},\hskip
8.5359pt\langle\ell_{k}|=\left(\begin{array}[]{cccc}\ell^{\mathrm{re}}_{k1}+\mathrm{i}\ell^{\mathrm{im}}_{k1},&\ell^{\mathrm{re}}_{k1}-\mathrm{i}\ell^{\mathrm{im}}_{k1},&\ell_{k3},&\ell_{k4}\end{array}\right),k=2,3,4.$
(61)

The elements $\ell_{kn}$ in Eq. (61) are obtained as

$\displaystyle\ell^{\mathrm{re}}_{k1}$ $\displaystyle=$
$\displaystyle-\frac{(4+(\tilde{\Gamma}-2\lambda_{k})^{2})(-4+\prod_{n=2,n\neq
k}^{4}(\tilde{\Gamma}-2\lambda_{n}))}{32\tilde{d}\prod_{n=2,n\neq
k}^{4}(\lambda_{k}-\lambda_{n})},\hskip
8.5359pt\ell^{\mathrm{im}}_{k1}=-\frac{(4+(\tilde{\Gamma}-2\lambda_{k})^{2})(\tilde{\Gamma}-\sum_{n=2,n\neq
k}^{4}\lambda_{n})}{8\tilde{d}\prod_{n=2,n\neq
k}^{4}(\lambda_{k}-\lambda_{n})},$ (62) $\displaystyle\ell_{k3}$
$\displaystyle=$
$\displaystyle-\frac{(4+(\tilde{\Gamma}-2\lambda_{k})^{2})(4+2\tilde{d}^{2}+\tilde{\Gamma}^{2}+4\prod_{n\neq
k}\lambda_{n})}{8(4+2\tilde{d}^{2}+\tilde{\Gamma}^{2})\prod_{n\neq
k}(\lambda_{k}-\lambda_{n})},\hskip
8.5359pt\ell_{k4}=\frac{(4+(\tilde{\Gamma}-2\lambda_{k})^{2})(4+2\tilde{d}^{2}+\tilde{\Gamma}^{2}-4\prod_{n\neq
k}\lambda_{n})}{8(4+2\tilde{d}^{2}+\tilde{\Gamma}^{2})\prod_{n\neq
k}(\lambda_{k}-\lambda_{n})}.$ (63)

Figure 12: The figure (a) shows the intersection time solutions [Eq. (74)] in
region (b) as functions of the control parameters $\tilde{d}_{\mathrm{I}}$ and
$\tilde{d}_{\mathrm{II}}$. Some of the regimes show finite and positive values
of both $\tilde{t}_{+}$ and $\tilde{t}_{-}$, implying double QMPE. Parameters
used are
$\tilde{d}=6.0,\tilde{\Gamma}_{\mathrm{I}}=\tilde{\Gamma}_{\mathrm{II}}=\tilde{\Gamma}=6\sqrt{17}$
(from Eq. (69)). The figure (b) shows a specific example of double
intersection (double QMPE) in $\rho_{\mathrm{gg}}(t)$. The main figure shows
the first intersection. The inset shows much weaker second intersection.
Parameters additionally used for (b) are
$\tilde{d}_{\mathrm{I}}=16.3,\tilde{d}_{\mathrm{II}}=13.4$.

To explore the possibility of QMPE, we consider the ground state probability
$\rho_{4}(t)\equiv\rho_{\mathrm{gg}}(t)$. The expression for
$\rho_{\mathrm{gg}}(t)$ becomes

$\displaystyle\rho_{\mathrm{gg}}(t)=\frac{4+\tilde{d}^{2}+\tilde{\Gamma}^{2}}{4+2\tilde{d}^{2}+\tilde{\Gamma}^{2}}+\sum_{k=2}^{4}e^{-\lambda_{k}t}\left[2(\ell^{\mathrm{re}}_{k1}\rho_{\mathrm{re}}(0)-\ell^{\mathrm{im}}_{k1}\rho_{\mathrm{im}}(0))+\ell_{k3}+(\ell_{k4}-\ell_{k3})\rho_{\mathrm{gg}}(0)\right].$
(64)

The difference $\Delta\rho_{\mathrm{gg}}(t)=$ has the form

$\displaystyle\Delta\rho_{\mathrm{gg}}(t)$ $\displaystyle=$ $\displaystyle
e^{-\lambda_{2}t}\left[\alpha_{4}x^{\frac{\lambda_{4}-\lambda_{2}}{\lambda_{3}-\lambda_{2}}}+\alpha_{3}x+\alpha_{2}\right],$
(65)

where

$\displaystyle x=e^{-(\lambda_{3}-\lambda_{2})t},\hskip
8.5359pt\alpha_{k}=2(\ell^{\mathrm{re}}_{k1}\Delta\rho_{\mathrm{re}}(0)-\ell^{\mathrm{im}}_{k1}\Delta\rho_{\mathrm{im}}(0))+(\ell_{k4}-\ell_{k3})\Delta\rho_{\mathrm{gg}}(0),$
(66) (67)

where $k=2,3,4,$ and $\Delta y=y^{\mathrm{I}}-y^{\mathrm{II}}$ in Eq. (67). We
want to create the possibility of multiple solutions to
$\Delta\rho_{\mathrm{gg}}(t)(t)=0$. To get $m$ number of intersections, we
have to produce at least an $m$-th order polynomial in Eq. (65). To do so,
clearly we need

$\lambda_{4}-\lambda_{3}=(m-1)(\lambda_{3}-\lambda_{2}).$ (68)

Interestingly, Eq. (68) puts a constraint on the gap statistics of the
eigenvalues of the Lindbladian. Thus the gap statistics of eigenvalues connect
directly to the possibility of multiple QMPE. One can numerically check the
maximum possible value of $m$ in region (b) by estimating
$(\lambda_{4}-\lambda_{3})/(\lambda_{3}-\lambda_{2})$ with the variation of
$(\tilde{d},\tilde{\Gamma})$ in this region.

To consider an example, we study $m=2$ i.e. the double QMPE. The line in
$\tilde{d}-\tilde{\Gamma}$ plane that satisfy this condition is

$\tilde{\Gamma}=3\sqrt{2}\sqrt{\tilde{d}^{2}-2}.$ (69)

Any $(\tilde{d},\tilde{\Gamma})$ on the line Eq. (69) simplifies Eq. (65) such
that the QMPE condition $\Delta\rho_{\mathrm{gg}}(t)=0$ reads as
$\alpha_{4}x^{2}+\alpha_{3}x+\alpha_{2}=0$. The solutions to this equation are

$x_{\pm}=\frac{(-\alpha_{3}\pm\sqrt{\alpha^{2}_{3}-4\alpha_{4}\alpha_{2}})}{2\alpha_{4}}.$
(70)

Figure 13: The figures (a) and (b) show nature of QMPE in energy and von
Neumann entropy in region (b), respectively, for the same parameters
exhibiting double QMPE in $\rho_{\mathrm{gg}}(t)$ in Fig. 12(b). We observe
that although energy shows double QMPE (figure (a)), the von Neumann entropy
exhibits single QMPE only (figure (b)).

The criteria for different types of QMPE in ground state probability
$\rho_{\mathrm{gg}}(t)$ in the region (b), are given by:

$\displaystyle
0<x_{+}<1\,\,\&\,\,0<x_{-}<1:\,\,\mathrm{double}\,\,\mathrm{QMPE},$ (71)
$\displaystyle
0<x_{+}<1\,\,\mathrm{or}\,\,0<x_{-}<1:\,\,\mathrm{single}\,\,\mathrm{QMPE},$
(72) $\displaystyle
x_{+}\not\in(0,1)\,\,\&\,\,x_{-}\not\in(0,1):\,\,\mathrm{no}\,\,\mathrm{QMPE}.$
(73)

Then, the intersection times $\tilde{t}_{\pm}$ are obtained as (using Eq.
(67))

$\displaystyle\tilde{t}_{\pm}=-\frac{1}{\lambda_{3}-\lambda_{2}}\mathrm{ln}(x_{\pm}).$
(74)

In Fig. 12(a), we present the intersection time solutions denoting QMPE in
$\rho_{\mathrm{gg}}(t)$ as functions of the control parameters
$\tilde{d}_{\mathrm{I}}$ and $\tilde{d}_{\mathrm{II}}$. Indeed, we observe
that there are parameter regions where both $\tilde{t}_{+}$ and
$\tilde{t}_{-}$ are positive and finite, indicating double QMPE. Other regions
have only one intersection time, either $\tilde{t}_{+}$ and $\tilde{t}_{-}$,
implying single QMPE. To demonstrate double QMPE, the Fig. 12(b) exhibits a
specific example of double intersection in $\rho_{\mathrm{gg}}(t)$. However,
the two copies become close enough so that the second intersection (shown in
inset of Fig. 12(b)) is weak in magnitude and hard to detect.

###  C.2 Energy

We look for solutions to $\Delta E_{\mathrm{s}}(t)=0$, which can give to
single or multiple QMPE depending on the number of positive finite
intersection time solutions. The Fig. 13(a) shows that the average energy
exhibits double QMPE indicated by the twice zero intersections of $\Delta
E_{\mathrm{s}}(t)$, for the same set of parameter values producing double QMPE
in $\rho_{\mathrm{gg}}(t)$ in Fig. 12(b).

Figure 14: The figure (a) shows thermal QMPE in region (b) where the
temperatures of the two copies intersect each other once leading to single
thermal QMPE. Parameters used are
$\tilde{d}=6.0,\tilde{\Gamma}_{\mathrm{I}}=0.5,\tilde{\Gamma}_{\mathrm{II}}=0.4,\tilde{\Gamma}=6\sqrt{17},\tilde{d}_{\mathrm{I}}=3.0,\tilde{d}_{\mathrm{II}}=2.0$.
The figure (b) exhibits single QMPE in KL divergence in region (b). Parameters
used for (b) are
$\tilde{d}=6.0,\tilde{\Gamma}_{\mathrm{I}}=0.5,\tilde{\Gamma}_{\mathrm{II}}=0.4,\tilde{\Gamma}=6\sqrt{17},\tilde{d}_{\mathrm{I}}=3.0,\tilde{d}_{\mathrm{II}}=2.0$.

###  C.3 von-Neumann entropy

Here we ask if the entropy, like density matrix and energy, also exhibit
double QMPE for the same set of parameters used in Figs. 12(b) and 13(a).
However, in Fig. 13(b), we observe that the entropy difference $\Delta
S_{\mathrm{vN}}(t)$ between two copies intersect zero only once, giving rise
to single QMPE instead of double QMPE.

###  C.4 Temperature

In Fig. 14(a), we present the time evolution of temperatures for the two
copies I and II, one initially starting as hotter with respect to the other.
We observe that the temperatures of the two copies intersect each other i.e.
the hotter becomes colder and vice versa. Notably, unlike ground state
probability [Fig. 12(b)], energy [Fig. 13(a)] that can exhibit double QMPE, we
have observed only single QMPE in temperature. If observation of multiple QMPE
is possible or not with purely exponential relaxations, remains an open
question.

###  C.5 Kullback-Leibler divergence

Here we study the possibility of QMPE in KL divergence for the same parameters
used to explore thermal QMPE in Figs. 14(a). The Fig. 14(b) shows QMPE in KL
divergence, however the intersection time corresponds to negative temperatures
in Fig. 14(a).

##  Appendix D QMPE in region $\mathrm{(e)}$: 3rd order EP

The region (e) in the $(\tilde{d},\tilde{\Gamma})$ plane of Hatano’s model is
composed of a single point

$\displaystyle\tilde{d}=2\sqrt{2},\,\,\tilde{\Gamma}=6\sqrt{3}.$ (75)

It is a third order EP where all the non-zero eigenvalues become equal and
their eigenvectors coalesce. Since this is a single point with fixed values of
$\tilde{d}$ and $\tilde{\Gamma}$, the right eigenvectors are explicitly given
by

$\displaystyle|r_{1}\rangle$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}\frac{-1-3\mathrm{i}\sqrt{3}}{15\sqrt{2}},&\frac{-1+3\mathrm{i}\sqrt{3}}{15\sqrt{2}},&\frac{1}{15},&1\\\
\end{array}\right)^{T},\hskip
11.38092pt|r_{2}\rangle=\left(\begin{array}[]{cccc}\frac{-1+\mathrm{i}\sqrt{3}}{\sqrt{2}},&\frac{-1-\mathrm{i}\sqrt{3}}{\sqrt{2}},&-1,&1\\\
\end{array}\right)^{T},$ (78) $\displaystyle|r_{3}\rangle$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}\frac{1+\mathrm{i}\sqrt{3}}{2\sqrt{2}},&\frac{-1+\mathrm{i}\sqrt{3}}{2\sqrt{2}},&0,&0\\\
\end{array}\right)^{T},\hskip
11.38092pt|r_{4}\rangle=\left(\begin{array}[]{cccc}\frac{1}{2\sqrt{2}},&\frac{1}{2\sqrt{2}},&0,&0\\\
\end{array}\right)^{T}.$ (81)

Of course $|r_{1}\rangle$ and $\langle\ell_{1}|$ have same expressions both in
Eqs. (6) and (81). The left eigenvectors are

Figure 15: The figure (a) illustrates the two intersection time solutions
$\tilde{t}_{\pm}$ [Eq. (98)] in region (e) as functions of
$\tilde{d}_{\mathrm{I}}$ and $\tilde{d}_{\mathrm{II}}$. For comparatively high
values of these control parameters, both solutions exist, leading to double
QMPE. Parameters used are
$\tilde{d}=2\sqrt{2},\tilde{\Gamma}=\tilde{\Gamma}_{\mathrm{I}}=\tilde{\Gamma}_{\mathrm{II}}=6\sqrt{3}$.
The figure (b) shows a surface from the three dimensional plot in figure (a).
It demonstrates the intersection times [Eq. (98)] by fixing
$\tilde{d}_{\mathrm{I}}$ and varying $\tilde{d}_{\mathrm{II}}$. We observe
single QMPE as well as double QMPE depending on the parameter choices.
Depending on the sign of $(\tilde{d}_{\mathrm{I}}-\tilde{d}_{\mathrm{II}})$,
the two solutions interchange their branches. Parameters used are
$\tilde{d}=2\sqrt{2},\tilde{\Gamma}=\tilde{\Gamma}_{\mathrm{I}}=\tilde{\Gamma}_{\mathrm{II}}=6\sqrt{3},\tilde{d}_{\mathrm{I}}=8.0$.

$\displaystyle\langle\ell_{1}|$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}0,&0,&\frac{15}{16},&\frac{15}{16}\\\
\end{array}\right),\hskip
8.5359pt\langle\ell_{2}|=\left(\begin{array}[]{cccc}0,&0,&-\frac{15}{16},&\frac{1}{16}\\\
\end{array}\right),$ (84) $\displaystyle\langle\ell_{3}|$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}\sqrt{2},&-\sqrt{2},&\frac{\mathrm{i}9\sqrt{3}}{4},&\frac{\mathrm{i}\sqrt{3}}{4}\\\
\end{array}\right),\hskip
8.5359pt\langle\ell_{4}|=\left(\begin{array}[]{cccc}\sqrt{2}+\mathrm{i}\sqrt{6},&\sqrt{2}-\mathrm{i}\sqrt{6},&5,&1\\\
\end{array}\right).$ (87)

The non-zero eigenvalue of the Lindbladian [Eq. (2)] in region (e) is
$-\mathrm{i}\lambda=-\mathrm{i}4\sqrt{3}$. Importantly, as a consequence of
the 3rd order EP, we can only achieve an approximately diagonal form i.e.
Jordan normal form $\widehat{\mathcal{L}}_{\mathrm{Jn}}$ leading to the
equation

$\displaystyle\widehat{\mathcal{L}}_{\mathrm{Jn}}=\widehat{L}\widehat{\mathcal{L}}\widehat{R},$
(88)

where $\widehat{\mathcal{L}}_{\mathrm{Jn}}$ has the explicit form as below

$\displaystyle\widehat{\mathcal{L}}_{\mathrm{Jn}}=\left(\begin{array}[]{cccc}0&0&0&0\\\
0&-\mathrm{i}4\sqrt{3}&1&0\\\ 0&0&-\mathrm{i}4\sqrt{3}&1\\\
0&0&0&-\mathrm{i}4\sqrt{3}\\\ \end{array}\right).$ (93)

Due to the off-diagonal nonzero defects in Eq. (93), the time evolution of
density matrix elements and other observables pick up additional algebraic
terms in time.

###  D.1 Density matrix elements

In region (e), the density matrix elements have the expressions below

$\displaystyle\rho_{j}(t)$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{4}e^{-\lambda_{k}t}r_{k,j}a_{k}-\mathrm{i}te^{-\lambda
t}(r_{2,j}a_{3}+r_{3,j}a_{4})-\frac{t^{2}}{2}e^{-\lambda t}r_{2,j}a_{4},$ (94)
$\displaystyle a_{k}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{4}\ell_{k,n}\rho_{n}(0).$ (95)

To investigate QMPE, we consider two initial conditions I and II as usual and
calculate the difference in time evolution between these two copies. We have
obtained the following expression for ground state probability difference
$\Delta\rho_{\mathrm{gg}}(t)=\rho^{\mathrm{I}}_{\mathrm{gg}}-\rho^{\mathrm{II}}_{\mathrm{gg}}$
in region (e):

$\displaystyle\Delta\rho_{\mathrm{gg}}(t)=e^{-4\sqrt{3}t}\left[\gamma_{2}t^{2}+\gamma_{1}t+\gamma_{0}\right],$
(96)
$\displaystyle\gamma_{2}=2\rho_{\mathrm{gg}}(0)-\sqrt{2}\left(\rho_{\mathrm{re}}(0)+\sqrt{3}\rho_{\mathrm{im}}(0)\right),\hskip
8.5359pt\gamma_{1}=2\sqrt{2}\rho_{\mathrm{im}}(0)-2\sqrt{3}\rho_{\mathrm{gg}}(0),\hskip
8.5359pt\gamma_{0}=\Delta\rho_{\mathrm{gg}}(0).$ (97)

The possibility of single or multiple QMPE in region (e) is dictated by the
equation $\Delta\rho_{\mathrm{gg}}(t)=0$ and its solutions which are the
intersection times, are

$\displaystyle\tilde{t}_{\pm}=\frac{-\gamma_{1}\pm\sqrt{\gamma^{2}_{1}-4\gamma_{2}\gamma_{0}}}{2\gamma_{2}}.$
(98)

The criteria for different types of QMPE in energy in region (e), are listed
below

$\displaystyle
0<\tilde{t}_{+}<\infty\,\,\&\,\,0<\tilde{t}_{-}<\infty:\,\,\mathrm{double}\,\,\mathrm{QMPE},$
(99) $\displaystyle
0<\tilde{t}_{+}<\infty\,\,\mathrm{or}\,\,0<\tilde{t}_{-}<\infty:\,\,\mathrm{single}\,\,\mathrm{QMPE},$
(100)
$\displaystyle\tilde{t}_{+}\not\in(0,\infty)\,\,\&\,\,\tilde{t}_{-}\not\in(0,\infty):\,\,\mathrm{no}\,\,\mathrm{QMPE}.$
(101)

Figure 16: The figure (a) shows an example of double QMPE in the ground state
density matrix element in region (e). The main figure shows the first
intersection between two copies I and II whereas the inset shows their second
intersection that is weaker in magnitude. The figure (b) shows QMPE in average
energy in region (e). The main figure shows that the energy difference between
two copies I and II intersect zero two times, resulting in double QMPE. The
insets show the intersections in the time evolution of the individual energies
of the two copies. Parameters used are
$\tilde{d}=2\sqrt{2},\tilde{\Gamma}=\tilde{\Gamma}_{\mathrm{I}}=\tilde{\Gamma}_{\mathrm{II}}=6\sqrt{3},\tilde{d}_{\mathrm{I}}=5.0,\tilde{d}_{\mathrm{II}}=7.0$.

In Fig. 15(a), we present existence of $\tilde{t}_{+}$ and $\tilde{t}_{-}$ as
function of the control parameters $\tilde{d}_{\mathrm{I}}$ and
$\tilde{d}_{\mathrm{II}}$ using the fixed dissipation protocol. The solutions’
surfaces in Fig. 15(a) reveal that there are regions where only one of the
solutions $\tilde{t}_{+}$ or $\tilde{t}_{-}$ exist leading to single QMPE, and
in other regions both the solutions exist resulting in double QMPE.
Particularly, lower values of $\tilde{d}_{\mathrm{I}}$ and
$\tilde{d}_{\mathrm{II}}$ correspond to single intersection time solution and
therefore single QMPE. To understand the behaviors of $\tilde{t}_{+,-}$ more
clearly, we fix $\tilde{d}_{\mathrm{I}}$ and display $\tilde{t}_{+,-}$ as
functions of $\tilde{d}_{\mathrm{II}}$ in Fig. 15(b). The Fig. 15(b) shows
that the two branches $\tilde{t}_{+}$ and $\tilde{t}_{-}$ interchange as
$\tilde{d}_{\mathrm{I}}<\tilde{d}_{\mathrm{II}}$. More specifically, for
$\tilde{d}_{\mathrm{I}}>\tilde{d}_{\mathrm{II}}$ we have
$\tilde{t}_{+}>\tilde{t}_{-}$ and for
$\tilde{d}_{\mathrm{I}}<\tilde{d}_{\mathrm{II}}$ we have
$\tilde{t}_{+}<\tilde{t}_{-}$. The underlying reason is that the discriminant
$\sqrt{\gamma^{2}_{1}-4\gamma_{2}\gamma_{0}}$ in Eq. (98) is proportional to
sign of $(\tilde{d}_{\mathrm{I}}-\tilde{d}_{\mathrm{II}})$. Thus the two
intersection times interchange their roles by changing the sign accompanying
the discriminant as we go from
$\tilde{d}_{\mathrm{I}}>\tilde{d}_{\mathrm{II}}$ to
$\tilde{d}_{\mathrm{I}}<\tilde{d}_{\mathrm{II}}$. The Fig. 16(a) demonstrates
an explicit example of double QMPE in ground state probability in region (e).
Like other regions, here also the second intersection (shown in the inset of
Fig. 16(a)) is much feeble in comparison to the first intersection.

###  D.2 Average energy

For the same parameters used to obtain QMPE for $\rho_{\mathrm{gg}}(t)$ in
Fig. 16(a), we would like to examine the possibility of QMPE in average energy
$E(t)$. Indeed, we observe in Fig. 16(b) that energy also exhibits double
QMPE, marked by the twice intersection of zero by the energy difference
$\Delta E(t)$ between the two copies. The insets of Fig. 16(b) shows
separately the two intersections in the time evolution of the energies for the
two copies. The second intersection is much weaker in magnitude compared to
the first intersection just like the case of ground state probability. Note
that the first intersection between the energies of I and II occur after they
overshoot their steady state value which they reach during some early time of
the relaxation process.

###  D.3 von Neumann entropy

Figure 17: The figure (a) illustrates QMPE in von Neumann entropy in region
(e). The difference of entropies between two copies I and II intersect zero
twice, leading to double QMPE. The inset shows the second intersection that is
weaker in magnitude with respect to the first intersection shown in the main
figure. Parameters used in figure (a) are same as in Fig. 16. The figure (b)
illustrates an example of thermal QMPE where the intersection between
initially hotter and colder copies occur after the initial stage of negative
temperature. The inset of figure (b) shows the full time evolution including
the initial negative temperatures. Parameters used are
$\tilde{d}=2\sqrt{2},\tilde{\Gamma}=6\sqrt{3},\tilde{\Gamma}_{\mathrm{I}}=14.0,\tilde{\Gamma}_{\mathrm{II}}=15.0,\tilde{d}_{\mathrm{I}}=6.0,\tilde{d}_{\mathrm{II}}=3.0$.

Here we investigate the von Neumann entropy $S_{\mathrm{vN}}(t)$, starting
from two copies I and II. We use same parameters that have produced double
QMPE both in ground state probability [Fig. 16(a)] and average energy [Fig.
16(b)]. In fact, the Fig. 17(a) shows that the entropy also exhibits double
QMPE indicated by twice zero intersections (the second one shown in inset of
Fig. 17(a)) of the entropy difference $\Delta S_{\mathrm{vN}}(t)$ between the
two copies.

###  D.4 Temperature

In this section we study the thermal QMPE in region (e). We observe thermal
QMPE in region (e) by using variable dissipation protocol in Fig. 17(b). There
we observe single intersection and reversal of identities (i.e. hotter
becoming colder and vice versa) occur after the negative temperatures have
happened. So, such temperatures can be considered as well-defined and
therefore this intersection in Fig. 17(b) corresponds to single QMPE in
temperature in region (e).

## References

  * (1) N. C. Keim, J. D. Paulsen, Z. Zeravcic, S. Sastry and S. R. Nagel, Rev. Mod. Phys. 91, 035002 (2019)..
  * (2) J. Gill,  Solid State Commun. 39, 1203 (1981).
  * (3) R. M. Fleming and L. F. Schneemeyer, Phys. Rev. B 28, 6996 (1983).
  * (4) J. Kaiser, An investigation into the occurrence of noises in tensile tests or a study of acoustic phenomena, Ph.D. thesis (PhD thesis, Technical University, Munich, Germany) (1950).
  * (5) A. Zang, B. Oye, P. Jousset, N. Deichmann, R. Gritto, A. McGarr, E. Majer and D, Bruhn, Geothermics 52, 6 (2014).
  * (6) E. R. Heimisson, P. Einarsson, F. Sigmundsson and B. Brandsdóttir, Geophys. Res. Lett., 42, 7958 (2015).
  * (7) L. Mullins,  Rubber Chem. Technol. 42, 339 (1969).
  * (8) A. Dorfmann and R. W. Ogden, Int. J. Solids Struct. 41, 1855 (2004).
  * (9) J. Diani, B. Fayolle, and P. Gilormini , European Polymer Journal 45, 601 (2009).
  * (10) A. J. Kovacs, Adv. Polym. Sci 3, 394 (1963).
  * (11) A. J. Kovacs, J. J. Aklonis, J. M. Hutchinson and A. R. Ramos, J. Polym. Sci. B 17, 1097 (1979).
  * (12) E. M. Bertin, J. P. Bouchaud, J. M. Drouffe, and C. Godrèche, J. Phys. A: Math. Gen. 36, 10701 (2003).
  * (13) A. Prados and E. Trizac, Phys. Rev. Lett. 112, 198001 (2014).
  * (14) A. Militaru, A. Lasanta, M. Frimmer, L. L. Bonilla, L. Novotny, and R. A. Rica, Phys. Rev. Lett. 127, 130603 (2021).
  * (15) K. Jonason, E. Vincent, J. Hammann, J. P. Bouchaud and P. Nordblad, Phys. Rev. Lett. 81, 3243 (1998).
  * (16) G. F. Rodriguez, G. G. Kenning and R. Orbach, Phys. Rev. Lett. 91, 037203 (2003).
  * (17) M. Baity Jesi, Spin Glasses: Criticality and Energy Landscapes(Springer, 2018).
  * (18) Y. Zheng and G. B. McKenna, Macromolecules 36, 2387 (2003).
  * (19) X. Di, K. Z. Win, G. B. McKenna, T. Narita, F. Lequeux, S. R. Pullela and Z. Cheng, Phys. Rev. Lett. 106, 095701 (2011).
  * (20) G. B. McKenna and S. L. Simon, Macromolecules 50, 6833 (2017).
  * (21) C. Scalliet and L. Berthier, Phys. Rev. Lett. 122, 255502 (2019).
  * (22) C. Josserand, A.V. Tkachenko, D.M. Mueth and H.M. Jaeger, Phys. Rev. Lett. 85, 3632 (2000).
  * (23) N. Kumar and S. Luding, Granul. Matter 18, 58 (2016).
  * (24) M. Otsuki and H. Hayakawa, Phys. Rev. E 101, 032905 (2020).
  * (25) T. Kawasaki and K. Miyazaki, arXiv:2003.10716.
  * (26) Pradipto and H. Hayakawa, Soft Matter 16, 945 (2020).
  * (27) L. M. C. Janssen, A. Kaiser and H. Löwen, Sci Rep. 7, 5667 (2017).
  * (28) Aristotle and W. D. Ross, Aristotle’s Metaphysics (Clarendon, Oxford, UK, 1981).
  * (29) E. B. Mpemba and D. G. Osborne, Phys. Educ. 4, 172 (1969).
  * (30) A. Kumar and J. Bechhoefer, Nature 584, 64 (2020).
  * (31) A. Kumar, R. Cheétrite and J. Bechhoefer, Proc. Natl. Acad. Sci. U.S.A. 119, e2118484119 (2022).
  * (32) A. Lasanta, F. Vega Reyes, A. Prados, and A. Santos, Phys. Rev. Lett. 119, 148001 (2017).
  * (33) A. Torrente, M. A. López-Castaño, A. Lasanta, F. Vega Reyes, A. Prados, and A. Santos, Phys. Rev. E 99, 060901(R) (2019).
  * (34) A. Biswas, V. V. Prasad, O. Raz, and R. Rajesh, Phys. Rev. E 102, 012906 (2020).
  * (35) A. Biswas, V. V. Prasad, and R. Rajesh, EPL 136, 46001 (2021).
  * (36) E. Mompó, M. A. López Castaño, A. Torrente, F. Vega Reyes, A. Lasanta, Phys. Fluids 33, 062005 (2021).
  * (37) A. Megías and A. Santos, Front. Phys. 10, 971671 (2022).
  * (38) A. Patrón, B. Sánchez-Rey, C. A. Plata and A. Prados, EPL 143, 61002 (2023).
  * (39) T. Keller, V. Torggler, S. B. Jäger, S. Schütz, H. Ritsch, and G. Morigi, New J. Phys. 20, 025004 (2018).
  * (40) A. Santos and A. Prados, Phys. Fluids 32, 072010 (2020).
  * (41) A. Patrón, B. Sánchez-Rey and A. Prados, Phys. Rev. E 104, 064127 (2021).
  * (42) S. Takada, H. Hayakawa and A. Santos, Phys. Rev. E 103, 032901 (2021) .
  * (43) S. Takada, EPJ Web. Conf. 249, 04001 (2021).
  * (44) Z. Lu and O. Raz, Proc. Natl. Acad. Sci. U.S.A. 114, 5083 (2017).
  * (45) I. Klich, O. Raz, O. Hirschberg, and M. Vucelja, Phys. Rev. X 9, 021060 (2019).
  * (46) D. M. Busiello, D. Gupta, amd A. Maritan, New J. Phys. 23, 103012 (2021).
  * (47) J. Lin, K. Li, J. He, J. Ren and J. Wang, Phys. Rev. E 105, 014104 (2022).
  * (48) A. Biswas, R. Rajesh and A. Pal, J. Chem. Phys. 159, 044120 (2023).
  * (49) R. G. González, N. Khalil and V. Garzó, Phys. Fluids 33, 053301 (2021)..
  * (50) 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. Muñoz 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, Proc. Natl. Acad. Sci. U.S.A. 116, 15350 (2019).
  * (51) P. A. Greaney, G. Lani, G. Cicero, and J. C. Grossman, Metall. Mater. Trans. A 42, 3907 (2011).
  * (52) Z.-Y. Yang and J.-X. Hou, Phys. Rev. E 101, 052106 (2020).
  * (53) Z.-Y. Yang and J.-X. Hou, Phys. Rev. E 105, 014119 (2022).
  * (54) I. González-Adalid Pemrtín, E. Mompó, A. Lasanta, V. Martín-Mayor and J. Salas, Phys. Rev. E 104, 044114 (2021).
  * (55) R. Cheétrite, A. Kumar and J. Bechhoefer, Front. Phys. 9, 654271 (2021).
  * (56) R. Holtzman and O. Raz, Communications Physics 5, 280 (2022).
  * (57) A. Megías, A. Santos, and A. Prados, Phys. Rev. E 105, 054140 (2022).
  * (58) A. Biswas, V. V. Prasad, and R. Rajesh, arXiv:2303.10900 (2023).
  * (59) I. G. Pemartin, E. Mompó, A. Lasanta, V. Martin-Mayor and J. Salas, arXiv.2308.04094.
  * (60) A. Nava and M. Fabrizio, Phys. Rev. B. 100, 125102 (2019).
  * (61) F. Carollo, A. Lasanta, and I. Lesanovsky, Phys. Rev. Lett. 127, 060401 (2021).
  * (62) S. K.Manikandan, Phys. Rev. Res. 3, 043108 (2021).
  * (63) F. Ivander, N. Anto-Sztrikacs, and D. Segal, arXiv:2301.06135 (2023).
  * (64) F. Ares, S. Murciano, and P. Calabrese, Nat. Commun. 14, 2036 (2023).
  * (65) A. K. Chatterjee, S. Takada and H. Hayakawa, Phys. Rev. Lett. 131, 080402 (2023).
  * (66) S. Murciano, F. Ares, I. Klich and P. Calabrese, arXiv:2310.07513.
  * (67) C. Rylands, K. Klobas, F. Ares, P. Calabrese, S. Murciano and B. Bertini, arXiv:2310.04419.
  * (68) W. D. Heiss, J. Phys. A: Math. Theor. 45, 444016 (2012).
  * (69) K. Hashimoto, K. Kanki, H. Hayakawa and T. Petrosky, Prog. Theor. Exp. Phys. 2015, 023A02 (2015).
  * (70) J. Doppler, A. Mailybaev, J. Böhm, U. Kuhl, A. Girschik, F. Libisch, T. J. Milburn, P. Rabl, M. Moyseyev, and S. Rotter, Nature 537, 76–79 (2016).
  * (71) Q. Zhong, J. Ren, M. Khajavikhan, D. N. Christodoulides, S. K. Özdemir, and R. El-Ganainy, Phys. Rev. Lett. 122, 153902 (2019).
  * (72) A. McDonald and A. A. Clerk, Nat. Commun. 11, 5382 (2020).
  * (73) S. Soleymani, Q. Zhong, M. Mokim, S. Rotter, R. El-Ganainy, and S. K. Özdemir, Nat. Commun. 13, 599 (2022).
  * (74) K. Ding, C. Fang, and G. Ma, Nat. Rev. Phys. 4, 745 (2022).
  * (75) F. Minganti, A. Miranowicz, R. W. Chhajlany, and F. Nori, Phys. Rev. A 100, 062131 (2019).
  * (76) I. I. Arkhipov, A. Miranowicz, F. Minganti and F. Nori, Phys. Rev. A 101, 013812 (2020).
  * (77) I. I. Arkhipov, F. Minganti, A. Miranowicz and F. Nori, Phys. Rev. A 104, 012205 (2021).
  * (78) M. Müller and I. rotter, J. Phys. A: Math. Theor. 41, 244018 (2008).
  * (79) M. Naghiloo, M. Abbasi, Y. N. Joglekar and K. W. Murch, Nature Phys. 15, 1232 (2019).
  * (80) M. Saha, B. K. Agarwalla, M. Kulkarni and A. Purkayastha, Phys. Rev. Lett. 130, 187101 (2023).
  * (81) F. Roccati, A. Purkayastha, G. Massimo Palma and F. Ciccarello, Eur. Phys. J. Spec. Top. 232, 1783 (2023).
  * (82) M. A. Miri and A. Alù, Science 363, eaar7709 (2019).
  * (83) S. K. Özdemir, S. Rotter, F. Nori, and L. Yang Nat. Mater. 18, 783–798 (2019).
  * (84) Yan-Li Zhou, Xiao-Die Yu, Chun-Wang Wu, Xie-Qian Li, Jie Zhang, Weibin Li, and Ping-Xing Chen, Phys. Rev. Res. 5, 043036 (2023).
  * (85) A. Lapolla and A. Godec, Phys. Rev. Lett. 125, 110602 (2020).
  * (86) T. Van Vu and Y. Hasegawa, Phys. Rev. Res. 3, 043160 (2021).
  * (87) S. Kochsiek, F. Carollo and I. Lesanovsky, Phys. Rev. A 106, 012207 (2022).
  * (88) F. Ghimenti and F. van Wijland, Phys. Rev. E 105, 054137 (2022).
  * (89) R. Bao and Z. Hou, arXiv:2212.11170.
  * (90) C. M. Bender, Rep. Prog. Phys. 70, 947 (2007).
  * (91) N. Moiseyev, Non-Hermitian Quantum Mechanics. Cambridge University Press; Cambridge, UK: 2011.
  * (92) R. Nakano, N. Hatano and T. Petrosky, Int. J. Theor. Phys. 50, 1134 (2011).
  * (93) N. Hatano and T. Petrosky, AIP Conf. Proc. 1648 200005 (2015).
  * (94) X. Li, Y. Li, and J. Jin, Phys. Rev. A 107, 032219 (2023).
  * (95) N. Hatano, Mol. Phys. 117, 2121 (2019).
  * (96) S. Kullback and R. A. Leibler, Ann. Math. Statist. 22, 79 (1951).
  * (97) P. Strasberg, Quantum Stochastic Thermodynamics: Foundations and Selected Applications (Oxford University Press, Oxford, 2022).
  * (98) T. Sagawa, Entropy, Divergence and Majorization in Classical and Quantum Thermodynamics (Springer Berlin, 2022).
  * (99) N. Shiraishi, An Introduction to Stochastic Thermodynamics: From Basic to Advanced (Springer, 2023).
  * (100) R. Yoshii and H. Hayakawa, Phys. Rev. Res. 5, 033014 (2023).
  * (101) G. Lindblad, Commun. Math. Phys. 48, 119 (1976).
  * (102) D. Manzano, AIP Advances 10, 025106 (2020) .
  * (103) V. Gorini, A. Kossakowski and E. C. G. Sudarshan, J. Math. Phys. 17, 821 (1976).
  * (104) D. Chruscinski and S. Pascazio, Open Sys. Inf. Dyn. 24, 1740001 (2017).
  * (105) R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey and D. E. Knuth, Adv. Comput. Math. 5, 329 (1996).
  * (106) A. A. More, Chem. Eng. Sci. 61, 5515 (2006).
  * (107) M. Asadian, H. Saeedi, M. Yadegari and M. Shojaee, J. Cryst. Growth 396, 61 (2014).
  * (108) A. Braun, A. Wokaun and H. G. Hermanns, Appl. Math. Model. 27, 47 (2003).
  * (109) A. Bot, B. Dewi and P. Venema, ACS Omega. 6, 7862 (2021).
  * (110) J. Moran, M. Romeijnders, P. Le Doussal, F. P. Pijpers, U. Weitzel, D. Panja and J-P. Bouchaud, arXiv:2309.15070.
  * (111) Shun-ichi Amari, Information Geometry and Its Applications, Applied Mathematical Sciences. Vol. 194. (Springer Japan, 2016).
  * (112) K. P. Burnham and D. R. Anderson, Model Selection and Multi-Model Inference, (Springer, 2002).
  * (113) A. N. Soklakov, Entropy 22, 860 (2020).
  * (114) L. D. Landau, and E. M. Lifshitz, Statistical Physics. Part 1, (Pergamon Press, Oxford, 449, 1980).
  * (115) G. Gauthier, M. T. Reeves, X. Yu, A. S. Bradley, M. A. Baker, T. A. Bell, H. Rubinsztein-Dunlop, M. J. Davis, and T. W. Neely, Science. 364, 1264 (2019).
  * (116) S. P. Johnstone, A. J. Groszek, P. T. Starkey, C. J. Billinton, T. P. Simula, and K. Helmerson, Science. 365, 1267 (2019).
  * (117) P. Medley, D. M. Weld, H. Miyake, D. E. Pritchard, and W. Ketterle, Phys. Rev. Lett. 106, 195301 (2011).
  * (118) M. M. Ali, W-M. Huang, and W-M. Zhang, Sci. Rep. 10, 13500 (2020).
  * (119) A. Gyenis, P. S. Mundada, A. Di Paolo, T. M. Hazard, X. You, D. I. Schuster, J. Koch, A. Blais and A. A. Houck, PRX Quantum 2, 010339 (2021).