# Non-Markovian Quantum Exceptional Points

Jhen-Dong Lin<EMAIL_ADDRESS>Department of Physics, National Cheng
Kung University, 701 Tainan, Taiwan  Center for Quantum Frontiers of Research
& Technology, NCKU, 70101 Tainan, Taiwan  Po-Chen Kuo  Department of Physics,
National Cheng Kung University, 701 Tainan, Taiwan  Center for Quantum
Frontiers of Research & Technology, NCKU, 70101 Tainan, Taiwan  Neill Lambert
Theoretical Quantum Physics Laboratory, RIKEN Cluster for Pioneering Research,
Wako-shi, Saitama 351-0198, Japan  RIKEN Center for Quantum Computing (RQC),
Wakoshi, Saitama 351-0198, Japan  Adam Miranowicz  Institute of Spintronics
and Quantum Information, Faculty of Physics, Adam Mickiewicz University,
61-614 Poznań, Poland  Franco Nori<EMAIL_ADDRESS>Theoretical Quantum Physics
Laboratory, RIKEN Cluster for Pioneering Research, Wako-shi, Saitama 351-0198,
Japan  RIKEN Center for Quantum Computing (RQC), Wakoshi, Saitama 351-0198,
Japan  Quantum Research Institute, The University of Michigan, Ann Arbor,
48109-1040 Michigan, USA  Yueh-Nan Chen<EMAIL_ADDRESS>Department
of Physics, National Cheng Kung University, 701 Tainan, Taiwan  Center for
Quantum Frontiers of Research & Technology, NCKU, 70101 Tainan, Taiwan
Physics Division, National Center for Theoretical Sciences, Taipei 106319,
Taiwan

###### Abstract

Exceptional points (EPs) are singularities in the spectra of non-Hermitian
operators, where eigenvalues and eigenvectors coalesce. Recently, open quantum
systems have been increasingly explored as EP testbeds due to their natural
non-Hermitian nature. However, existing works mostly focus on the Markovian
limit, leaving a gap in understanding EPs in the non-Markovian regime. In this
work, we address this gap by proposing a theoretical framework based on two
numerically exact descriptions of non-Markovian dynamics: the pseudomode
mapping and the hierarchical equations of motion. The proposed framework
enables conventional spectral analysis for EP identification, establishing
direct links between EPs and dynamic manifestations in open systems, such as
non-exponential decays and enhanced sensitivity to external perturbations. We
unveil pure non-Markovian EPs that are unobservable in the Markovian limit.
Remarkably, the EP aligns with the Markovian-to-non-Markovian transition, and
the EP condition is adjustable by modifying environmental spectral properties.
Moreover, we show that structured environments can elevate EP order, thereby
enhancing the system’s sensitivity. These findings lay a theoretical
foundation and open new avenues for non-Markovian reservoir engineering and
non-Hermitian physics.

##  I Introduction

Spectral singularities for non-Hermitian systems, known as exceptional points
(EPs), have attracted intense research attention over the past decades [1, 2,
3, 4]. These singularities are pivotal in studying open systems, as
environmental noise inherently breaks the Hermiticity. Early studies of EPs
mainly focus on non-Hermitian Hamiltonians (NHHs) [5, 6, 7, 8] that are
suitable for modeling classical and semiclassical systems. In these settings,
(semi-)classical EPs, often termed Hamiltonian EPs (HEPs), are identified by
the convergence of at least two eigenvalues and their corresponding
eigenstates within NHHs. Notably, it has been shown that environmental noises
can induce exotic effects near HEPs, e.g., EP-induced lasing [9, 10, 11],
programmable mode switching [12], and EP-enhanced sensitivity [13, 14, 15, 16,
17].

Recently, investigations of EPs have extended into the full quantum regime
[18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31], where the temporal
evolution of an open quantum system is governed by a Lindblad master equation
or, equivalently, by a Liouvillian superoperator. Unlike the effective NHHs,
Liouvillian superoperators incorporate the concept of quantum jumps into the
dynamical process [32]. In this context, the EPs associated with Liouvillian
superoperators are termed as quantum EPs or Liouvillian EPs (LEPs). It has
been demonstrated that _pure quantum EPs_ exist [18], which are phenomena
without (semi-)classical counterparts.

To date, the exploration for both HEPs and LEPs has largely been confined to
the memoryless Markovian limit, which is only valid in cases of sufficiently
weak system-environment interaction or environments without any structure.
Recently, several works have indicated that EP-like critical behaviors [33,
34, 35, 36] could manifest in the non-Markovian regime. However, whether the
concepts of EPs can be directly applied to the non-Markovian regime remains an
open question. A primary challenge lies within the structure of the non-
Markovian equation of motion for reduced dynamics, which can be generally
expressed as:

$\frac{d\rho_{\text{S}}(t)}{dt}=\int_{0}^{t}d\tau~{}\mathcal{K}(t,\tau)\rho_{\text{S}}(\tau).$
(1)

Here, the non-Markovian effect is encoded in the memory kernel
$\mathcal{K}(t,\tau)$ [32], and $\rho_{\text{S}}(t)$ denotes the open system’s
reduced density operator. The integral-differential nature of this time-non-
local equation complicates the application of traditional spectral analysis
techniques for identifying EPs.

In this work, we aim to address this theoretical gap. The main result lies in
the development of a systematic framework for quantum EPs associated with
generic non-Markovian open systems. The idea is based on applying the
pseudomode equation of motion (PMEOM) [37, 38, 39, 40, 41, 42, 43, 44, 45, 46,
47] and the hierarchical equations of motion (HEOM) [48, 49, 50, 51, 52, 53,
54, 55, 56, 57], which can be used to describe a large class of system-
environment models. In contrast to the memory kernel approach as described in
Eq. (1), these methods encode the non-Markovian effects in auxiliary degrees
of freedom so that both PMEOM and HEOM can be formulated by time-local
differential equations. Specifically, the dynamics is governed by what we call
_extended Liouvillian superoperators_ , enabling us to perform conventional
spectral analysis, identifying the corresponding EPs, and revealing their
impacts on the non-Markovian open systems.

Moreover, we demonstrate the utility through two analytically tractable
examples. The first one is the spin-boson model with a Lorentzian environment.
We identify a _pure non-Markovian EP_ that cannot be observed in the Markovian
limit. Intriguingly, the EP precisely aligns with the Markovian-to-non-
Markovian transition, thereby revealing a close relationship between the onset
of non-Markovian information backflow and non-Hermitian phase transition. In
addition, we show that the EP criterion is tunable by reservoir engineering.
Specifically, we introduce a band gap to the Lorentzian environment,
demonstrating that the corresponding EP condition requires a smaller system-
environment coupling strength compared to the gap-less scenario.

For the second example, we examine linear bosonic systems by using the adjoint
PMEOM within the Heisenberg picture. We show that the dynamics of the modes’
amplitudes can be determined by effective NHHs, enabling us to explore
potential non-Markovian HEPs. As an example, we consider a two-coupled-modes
system, where a second-order HEP emerges in the Markovian limit. We showcase
that in the non-Markovian regime, the HEP can be transformed to the third
order. Consequently, the system becomes more sensitive to external
perturbations with the help of the quantum memory effect. These findings
reveal the intricate interplay between EPs and the memory effect, laying a
theoretical foundation for exploring non-Hermitian physics towards non-
Markovian and non-perturbative regimes.

##  II Results

### General framework for non-Markovian exceptional points

####  II.0.1 Pseudomode equation of motion

Figure 1: Schematic illustration depicting EPs for (a) a generic non-Markovian
open-system model, where the structured environment is captured by the
spectral density function $J(\omega)$. For a given spectral density function
and the corresponding environmental correlation function, the exact non-
Markovian dynamics can either be described by (b) PMEOM or (c) HEOM with the
corresponding extended Liouvillian superoperators:
$\mathcal{L}_{\text{S}+\text{PM}}$ and $\mathcal{L}_{\text{S}+\text{ADO}}$.
(d) The non-Markovian EPs can then be identified by observing the complex
spectrum of these extended Liouvillian superoperators.

In this section, we present a general approach to characterize EPs for open
systems subject to non-Markovian noise. To begin with, we consider an open
quantum system (S) coupled to a bosonic environment (E), where the total
Hamiltonian is expressed by

$\displaystyle
H_{\text{tot}}=H_{\text{S}}+H_{\text{E}}+H_{\text{S}\text{E}}\text{~{}with~{}~{}}H_{\text{E}}=\sum_{k}\omega_{k}b^{\dagger}_{k}b_{k},$
(2) $\displaystyle H_{\text{S}\text{E}}=QX,~{}\text{and
}X=\sum_{k}g_{k}(b^{\dagger}_{k}+b_{k}).$

Here, $H_{\text{S}}$ and $H_{\text{E}}$ represent the free Hamiltonians of the
system and its environment, while $H_{\text{S}\text{E}}$ describes their
interactions. Also, $\omega_{k}$ and $g_{k}$ correspond to the frequency and
the coupling strength for the environmental mode $b_{k}$, respectively, and
$Q$ represents an arbitrary system-environment coupling operator. We consider
that the environment is initialized in a Gibbs state at temperature $T$:
$\rho_{\text{E}}=\exp(-\beta H_{\text{E}})/\operatorname{tr}[\exp(-\beta
H_{\text{E}})]$, where $\beta=(k_{\text{B}}T)^{-1}$ with $k_{\text{B}}$
denoting the Boltzmann constant. In this case, the exact dynamics of the open
system’s reduced density matrix (within the interaction picture) can be
written as

$\displaystyle\rho_{\text{S}}(t)=\hat{\mathcal{T}}\exp\left\\{\hat{\mathcal{F}}\Big{[}Q,C(t)\Big{]}\right\\}\rho_{\text{S}}(0).$
(3)

Here, $\hat{\mathcal{T}}$ denotes the time-ordering operator and
$\hat{\mathcal{F}}$ represents the Feynman-Vernon influence functional [32],
which can be expressed by

$\displaystyle\hat{\mathcal{F}}=-\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}~{}Q(t_{1})^{\times}[C_{\mathbb{R}}(t_{1}-t_{2})Q(t_{2})^{\times}$
(4)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+iC_{\mathbb{I}}(t_{1}-t_{2})Q(t_{2})^{\circ}],$

where $Q(t)=\exp(iH_{\text{S}}t)Q\exp(-iH_{\text{S}}t)$, and we adopt the
superoperator notations:

$\displaystyle
Q(t)^{\circ}=\\{Q(t),\bullet\\}\text{~{}and~{}}Q(t)^{\times}=[Q(t),\bullet].$
(5)

In addition, $C_{\mathbb{R}(\mathbb{I})}$ denotes the real (imaginary) part of
the environmental correlation function

$C(t)=\operatorname{tr}[X(t)X(0)\rho_{\text{E}}]$ (6)

with $X(t)=\exp(iH_{\text{E}}t)X\exp(-iH_{\text{E}}t)$. An essential feature
of $\hat{\mathcal{F}}$ is its exclusive dependence on the system-environment
coupling operator $Q$ and the environmental correlation function $C(t)$. The
latter can be expressed by

$C(t)=\int_{0}^{\infty}d\omega~{}\frac{J(\omega)}{\pi}\left[\coth\left(\frac{\beta\omega}{2}\right)\cos(\omega
t)-i\sin(\omega t)\right],$ (7)

where $J(\omega)=\sum_{k}|g_{k}|^{2}\delta(\omega-\omega_{k})$ represents the
coupling spectral density. This property enables the reproduction of the exact
open system dynamics through an auxiliary model involving a small set of
fictitious damping modes, i.e., the pseudomodes, provided that the correlation
function for the artificial model aligns with Eq. (7).

For a broad range of cases, the correlation function can be efficiently
expressed as a finite weighted summation of exponential terms, i.e.,

$C(t)=\sum_{i}\alpha_{i}^{2}\exp(-i\Omega_{i}t-\gamma_{i}|t|/2).$ (8)

With this expression, one can construct the PMEOM:

$\displaystyle\frac{d}{dt}\rho_{\text{S}+\text{PM}}(t)=\mathcal{L}_{\text{S}+\text{PM}}[\rho_{\text{S}+\text{PM}}(t)]$
(9)
$\displaystyle=-i[H_{\text{S}+\text{PM}},\rho_{\text{S}+\text{PM}}(t)]+\sum_{i}\gamma_{i}\mathcal{L}_{a_{i}}[\rho_{\text{S}+\text{PM}}(t)],$
$\displaystyle\text{with~{}~{}}H_{\text{S}+\text{PM}}=H_{\text{S}}+\sum_{i}\Omega_{i}a_{i}^{\dagger}a_{i}+\alpha_{i}Q(a^{\dagger}_{i}+a_{i}),$

where, $\\{a_{i}\\}$ represent the pseudomodes, and we introduce the
dissipator $\mathcal{L}_{a_{i}}[{\bullet}]=a_{i}\bullet
a_{i}^{\dagger}-\\{a_{i}^{\dagger}a_{i},\bullet\\}/2$. Notably, the
environmental influences on the open system are now captured by the
pseudomodes’ frequencies $\Omega_{i}$, damping rates $\gamma_{i}$, and the
system-pseudomode coupling strengths $\alpha_{i}$. The exact dynamics of S can
be obtained by tracing out the pseudomodes (PM), i.e.,
$\rho_{\text{S}}(t)=\operatorname{tr}_{\text{PM}}[\rho_{\text{S}+\text{PM}}(t)]$,
after solving the PMEOM with these pseudomodes initialized in the thermal
states.

One notable benefit of employing the pseudomode model lies in its facilitation
of establishing physical intuitions regarding the EP criteria. For instance,
recent works have suggested that EPs are closely related to critical damping
points for both classical and quantum systems [23]. As we demonstrate below,
it is feasible to pinpoint EPs by balancing the system-pseudomode coupling
strength and the pseudomode damping, leading the system to a critical damping
point.

Furthermore, a systematic procedure for characterizing non-Markovian EPs based
on a standard spectral analysis can be established. The main idea is grounded
in the observation that both the temporal evolution of
$\rho_{\text{S}+\text{PM}}(t)$ and $\rho_{\text{S}}(t)$ are governed by the
spectral properties of the extended Liouvillian superoperator
$\mathcal{L}_{\text{S}+\text{PM}}$. Specifically, assuming that
$\mathcal{L}_{\text{S}+\text{PM}}$ is diagonalizable, we consider its
eigenvalues and the corresponding eigenmatrices:
$\\{\lambda_{i},\hat{\rho}_{\text{S}+\text{PM},i}\\}_{i}$. The dynamics of
$\rho_{\text{S}+\text{PM}}(t)$ can then be expressed by

$\rho_{\text{S}+\text{PM}}(t)=\sum_{i}c_{i}\exp(\lambda_{i}t)\hat{\rho}_{\text{S}+\text{PM},i}.$
(10)

By tracing out these pseudomodes, one can observe that the exact dynamics of
the system reduced state follows a similar expression, replacing these
eigenmatrices $\hat{\rho}_{\text{S}+\text{PM},i}$ with the reduced
eigenmatrices
$\hat{\rho}_{\text{S},i}=\operatorname{tr}_{\text{PM}}(\hat{\rho}_{\text{S}+\text{PM},i})$,
namely,

$\rho_{\text{S}}(t)=\sum_{i}c_{i}\exp(\lambda_{i}t)\hat{\rho}_{\text{S},i}.$
(11)

To describe EPs, one considers a family of parametrized extended Liouvillian
superoperators $\mathcal{L}_{\text{S}+\text{PM}}(\bm{\xi})$, bearing in mind
that $\bm{\xi}$ includes the parameters related to both the system and the
structured environments. Since $\mathcal{L}_{\text{S}+\text{PM}}(\bm{\xi})$ is
generally non-Hermitian, EPs could potentially exist in the parameter space.
For instance, let $\bm{\xi}_{\text{EP}n}$ represent an $n$th-order EP in the
parameter space, where $n$ different eigenvalues and the corresponding
eigenmatrices $\\{\lambda_{i},\hat{\rho}_{\text{S}+\text{PM},i}\\}_{i\in A}$
coalesce into
$\\{\lambda_{\text{EP}},\hat{\rho}_{\text{S}+\text{PM},\lambda_{\text{EP}}}\\}$.
Here, $A$ denotes a set of indices. Due to the coalescence of the
eigenmatrices, the corresponding $n$-dimensional eigensubspace for
$\mathcal{L}_{\text{S}+\text{PM}}(\bm{\xi}_{\text{EP}})$ cannot be
diagonalized. Nevertheless, a Jordan block for the subspace can be constructed
by introducing generalized eigenmatrices
$\\{\hat{\rho}^{(j)}_{\text{S}+\text{PM},\lambda_{\text{EP}}}\\}_{j=0,\cdots,n-1}$
[7], such that the system dynamics can be expressed by

$\displaystyle\rho_{\text{S}}(t)=\sum_{i\notin
A}c_{i}e^{\lambda_{i}t}\hat{\rho}_{\text{S},i}+e^{\lambda_{\text{EP}}t}\sum_{j=0}^{n-1}\sum_{m=0}^{j}\frac{t^{m}\tilde{c}_{m}}{m!}\hat{\rho}^{(j)}_{\text{S},\lambda_{\text{EP}}},$
(12)

where the reduced generalized eigenmatrices are introduced as

$\hat{\rho}^{(j)}_{\text{S},\lambda_{\text{EP}}}=\operatorname{tr}_{\text{PM}}(\hat{\rho}^{(j)}_{\text{S}+\text{PM},\lambda_{\text{EP}}}).$
(13)

Equation (12) suggests that the polynomial time dependence, a common dynamical
signature of EPs, could be observed in the reduced dynamics of the open
system. In essence, the PMEOM provides an intuitive and direct route to
investigate EPs beyond the Markovian limit, and this approach is compatible
with the conventional spectral analysis by introducing the (generalized)
reduced eigenmatrices.

####  II.0.2 Hierarchical equations of motion

It is worth noting that an alternative approach is to consider the framework
of HEOM. In this context, a set of auxiliary density operators (ADOs) is
introduced to capture the non-Markovian and non-perturbative effects [53, 57,
55, 58, 59, 60]. Similarly, we can define the extended quantum state
$\rho_{\text{S}+\text{ADO}}$ that contains both the system reduced state and
the ADOs. The dynamics of the extended state is governed by the extended
Liouvillian superoperator $\mathcal{L}_{\text{S}+\text{ADO}}$. The system
reduced state can be obtained through a linear operation, specifically
$\rho_{\text{S}}(t)=\mathcal{P}[\rho_{\text{S}+\text{ADO}}(t)]$, where
$\mathcal{P}$ is a superoperator for discarding all ADOs (see Methods).
Therefore, the EPs for non-Markovian open quantum systems can also be
equivalently characterized under the framework of the HEOM by introducing the
corresponding (generalized) reduced eigenmatrices, which are expressed as
$\tilde{\rho}_{\text{S},i}=\mathcal{P}(\tilde{\rho}_{\text{S}+\text{ADO},i})$
and
$\tilde{\rho}^{(j)}_{\text{S},\lambda_{\text{EP}}}=\mathcal{P}(\tilde{\rho}^{(j)}_{\text{S}+\text{ADO},\lambda_{\text{EP}}})$.

### Spin-boson model

Figure 2: (a) Lorentzian $J_{L}(\omega)$ and band gap $J_{q}(\omega)$ spectral
densities centered at the qubit transition frequency $\omega_{0}$. (d) The
effects of these spectral densities can be represented by two pseudomodes: The
Lorentzian spectral density can be described by the upper mode with qubit-mode
coupling strength $\alpha_{1}=\sqrt{\Gamma\Lambda/2}$ and mode’s damping rate
$\gamma_{1}=2\Lambda$, while the band gap is characterized by the lower mode
with the non-Hermitian coupling $\alpha_{2}=i\sqrt{q\Gamma\Lambda/2}$ and
damping rate $\gamma_{2}=2q\Lambda$, where the red dashed line signifies the
unphysical nature of this mode. (b, c) The real and imaginary parts of the
spectrum of the extended Liouvillian for the gapless scenario ($q=0$). Two
exceptional points, EP2 and an EP3 emerge at the coupling strength
$\Gamma=\Lambda/2$. (e, f) The real and imaginary parts of the spectrum of the
extended Liouvillian when $q=1/4$. The EP criterion becomes
$\Gamma=(1-q)\Lambda/2$. The dotted curves in (b, c, e, f) represent the
spectrum of the extended superoperators for HEOM (see Supplementary Note 1 for
detailed derivations).

To gain a more intuitive understanding, we apply the PMEOM to the prototype
spin-boson model, showing that a pure non-Markovian EP can be observed by
adjusting the structural characteristics of the environment’s spectral
density. To simplify our analysis, we consider a zero-temperature environment
and adopt the rotating-wave approximation (RWA), where the system-environment
interaction in Eq. (2) and the system-pseudomode coupling in Eq. (9) can be
expressed by
$\sum_{k}g_{k}(\tilde{Q}b_{k}^{\dagger}+\tilde{Q}^{\dagger}b_{k})$ and
$\sum_{i}\alpha_{i}(\tilde{Q}a_{i}^{\dagger}+\tilde{Q}^{\dagger}a_{i})$,
respectively. Furthermore, we consider a scenario where the spectral density
function $J(\omega)$ is well localized in the vicinity of a high frequency
$\tilde{\omega}\gg 0$, enabling us to effectively approximate the correlation
in Eq. (7) by extending the lower limit of the integral to negative infinity.
Note that the assumptions mentioned above are also considered in the original
proposals of PMEOM [37, 38, 39].

The model involves a qubit representing the open system, where the system
Hamiltonian and system-environment coupling operator are
$H_{\text{S}}=\omega_{0}|e\rangle\langle e|$ and $\tilde{Q}=\sigma_{-}$,
respectively. Here, $\omega_{0}$ denotes the qubit transition frequency
between the ground state $|g\rangle$ and the excited state $|e\rangle$, and
$\sigma_{+}=|e\rangle\langle g|$ ($\sigma_{-}=|g\rangle\langle e|$) represents
the raising (lowering) operator. We consider a Lorentzian spectral density
that is expressed by

$J_{L}(\omega)=\frac{1}{2}\frac{\Gamma\Lambda^{2}}{(\omega-\omega_{0})^{2}+\Lambda^{2}},$
(14)

where $\Gamma$ and $\Lambda$ denote the coupling strength and the spectral
width, respectively. In the interaction picture and under the earlier
specified assumption of extending the spectral density function to negative
frequencies, the environmental correlation function can be expressed by a
single exponential term, i.e., $C(t)=(\Gamma\Lambda/2)\exp(-\Lambda|t|)$.
Therefore, the PMEOM can be constructed by introducing a single pseudomode
with the damping rate $\gamma=2\Lambda$ and the qubit-pseudomode coupling
strength $\alpha=\sqrt{\Gamma\Lambda/2}$. As aforementioned, the pseudomode
representation provides a physical intuition that the EP could be located by
balancing $\gamma$ and $\alpha$ (or $\Gamma$ and $\Lambda$) so that the system
reaches a critical damping point.

The corresponding extended Liouvillian superpoerator
$\mathcal{L}_{\text{S}+\text{PM}}(\Gamma,\Lambda)$ can be described by a
$9\times 9$ non-Hermitian matrix in the single-excitation subspace
(Supplementary Note 1). The spectrum is illustrated in Fig. 2, revealing that
$\Gamma=\Lambda/2$ corresponds to a second-order EP (EP2) and a third-order EP
(EP3). Notably, these EPs are purely non-Markovian, because they are
unobservable in the Markovian limit. Specifically, in such a limit, the
spectral width (and thus the damping rate of the pseudomode) becomes infinite,
$\Lambda\rightarrow\infty$. Therefore, the pseudomode can be adiabatically
eliminated and the dynamics is governed by a qubit-only Markovian master
equation, i.e.,
$\dot{\rho}_{\text{S}}(t)=\Gamma[2\sigma_{-}\rho_{\text{S}}(t)\sigma_{+}-\\{\sigma_{+}\sigma_{-},\rho_{\text{S}}(t)\\}]/2$.
Intuitively, there is only one qubit decay channel without internal tunneling
between the qubit energy levels, thereby EP does not emerge in this scenario
[18].

Let us proceed with a deeper analysis of the EP dynamical signatures.
According to Figs. 2(b) and 2(c), the EP condition pinpoints a real-to-complex
transition in the spectrum. In other words, the EP condition corresponds to a
critical damping point, separating the overdamped (pure decay) and underdamped
(oscillatory) regimes. By investigating the generalized eigenmatrices for the
extended superoperator at the EP condition, i.e.,
$\mathcal{L}_{\text{S}+\text{PM}}(\Gamma=\Lambda/2)$, we conclude that the
qubit coherence and population dynamics are respectively determined by the EP2
and EP3, resulting in observing first-order and second-order polynomial time
dependencies:

$\displaystyle\langle e|\rho_{\text{S}}(t)|g\rangle=\frac{1}{2}(\Lambda
t+2)e^{-\frac{1}{2}\Lambda t}\langle e|\rho_{\text{S}}(0)|g\rangle,$ (15)
$\displaystyle\langle
e|\rho_{\text{S}}(t)|e\rangle=\frac{1}{4}(\Lambda^{2}t^{2}+4\Lambda
t+4)e^{-\Lambda t}\langle e|\rho_{\text{S}}(0)|e\rangle.$

Notably, by utilizing the Breuer-Lain-Pillo [61] and the Riva-Huelga-Plenio
[62] non-Markovian measures to the qubit dynamics, we find that the EP
condition also aligns with the onset of non-Markovianity (see Supplementary
Note 2). This reveals an unexpected connection between the Markovian-to-non-
Markovian and non-Hermitian EP phase transitions.

Furthermore, the transition from the overdamped to underdamped regimes opens
up the possibility of observing enhanced sensitivity to external perturbations
in the vicinity of the EP. For instance, consider a small perturbation in the
coupling strength, represented as $\Gamma\rightarrow\Lambda(1+\epsilon)/2$,
with $\epsilon>0$. This perturbation induces splittings in the imaginary parts
of the eigenvalues, resulting in oscillatory dynamics. Notably, these
splittings scale as $\sqrt{\epsilon}$, indicating the sensitivity enhancement
in the characteristic frequency of the oscillatory behavior. For this model,
both the qubit excited state population and coherence vanish periodically. We
choose, for simplicity, the first vanishing time $t_{\text{vanish}}$ of the
qubit coherence to capture the system’s oscillations. One can find
$t_{\text{vanish}}^{-1}\approx\Lambda\sqrt{\epsilon}/2+O\left(\epsilon\right)$,
which implies that $t_{\text{vanish}}$ is sensitive to the external
perturbation when the system is prepared at the EP.

Up to this point, we only tune the environmental parameters without changing
the shape of the Lorentzian profile, and we have shown that such an adjustment
does not modify the structure of the PMEOM. Let us further consider a more
complex scenario, where an additional parameter can drastically modify the
spectral shape. To this end, we introduce a band-gap structure to the
Lorentzian background [32], which is modeled as

$\displaystyle
J_{q}(\omega)=J_{L}(\omega)-\frac{1}{2}\frac{\Gamma(q\Lambda)^{2}}{(\omega-\omega_{0})^{2}+(q\Lambda)^{2}},$
(16)

where the band gap is located at the frequency $\omega_{0}$, such that
$J_{q}(\omega_{0})=0$ [32] [see also Fig. 2(a)], and the parameter $q\in(0,1]$
is used to control the relative width for the band gap. The environmental
correlation function is now expressed by two exponential terms:

$C(t)=\frac{\Lambda\Gamma}{2}\exp(-\Lambda|t|)-\frac{q\Lambda\Gamma}{2}\exp(-q\Lambda|t|).$
(17)

Accordingly, the PMEOM for the exact dynamics involves two pseudomodes
($a_{1}$ and $a_{2}$), each characterized by the qubit-pseudomode coupling
strengths and pseudomode damping rates [see also Fig. 2(d)]:

$\\{\alpha_{1}=\sqrt{\Gamma\Lambda/2},\alpha_{2}=i\sqrt{q}\alpha_{1},\gamma_{1}=2\Lambda,\gamma_{2}=q\gamma_{1}\\}.$
(18)

The spectrum of the corresponding extended Liouvillian superoperator is
presented in Figs. 2(e) and 2(f), suggesting that the EP criterion becomes

$\Gamma=(1-q)\Lambda/2.$ (19)

Consequently, including a band gap into the environmental spectral profile
leads to a displacement of the EP in the parameter space. In this case, the
system requires a smaller coupling strength to reach the EP criterion compared
to the case of the previous gapless Lorentzian environment.

We remark that $\alpha_{2}$ is purely imaginary, so it does not directly
correspond to a physical mode. This characteristic causes the
$H_{\text{S}+\text{PM}}$ to become non-Hermitian [42]. Nevertheless, the
reduced dynamics of the qubit is still exact and well-behaved as long as the
PMEOM is consistent with the given $C(t)$ [42]. Furthermore, in Fig. 2 and
Supplementary Note 1, we demonstrate that the spectrum of the extended
Liouvillian operators for the HEOM is consistent with that of PMOEM.

### Linear bosonic systems with adjoint pseudomode-equation of motion

The PMEOM in Eq. (9) offers significant flexibility by allowing both the
system Hamiltonian $H_{\text{S}}$ and the coupling operator $Q$ to be
arbitrary. This utility facilitates the extension of our analysis to a variety
of open quantum systems. As an illustrative application, we now explore linear
bosonic systems, which have garnered significant interest in the study of both
LEPs and HEPs [12, 19, 21, 22, 25]. Specifically, we examine a general
Hamiltonian with $M$ coupled modes, expressed as

$H_{\text{S}}=\sum_{k=1}^{M}\Omega_{k}c^{\dagger}_{k}c_{k}+\sum_{j<k}\chi_{j,k}(c^{\dagger}_{j}c_{k}+c_{k}^{\dagger}c_{j}),$
(20)

where $\Omega_{k}$ denotes the frequency associated with the mode $c_{k}$, and
$\chi_{j,k}$ represents the coherent coupling between the modes $j$ and $k$.
In addition, we consider that the $M$th mode further couples to the
environment with the coupling operator $Q=c_{M}$. In the Heisenberg picture,
the dynamics is governed by the following adjoint PMEOM:

$\displaystyle\frac{d}{dt}O_{\text{S}+\text{PM}}(t)=i[H_{\text{S}+\text{PM}},O_{\text{S}+\text{PM}}(t)]$
(21)
$\displaystyle+\sum_{i}\frac{\gamma_{i}}{2}(2a^{\dagger}_{i}O_{\text{S}+\text{PM}}(t)a_{i}-\\{a_{i}^{\dagger}a_{i},O_{\text{S}+\text{PM}}(t)\\}).$

To facilitate the analysis, we define a vector containing the annihilation
operators for both the system modes and the pseudomodes:

$\mathbf{v}=\big{(}c_{1},\cdots,c_{M},a_{1},\cdots,a_{i},\cdots\big{)}^{T}.$
(22)

Focusing on the amplitudes of these modes, we denote

$\displaystyle\langle\mathbf{v(t)}\rangle=(\langle
c_{1}(t)\rangle,\cdots,\langle c_{M}(t)\rangle,\langle
a_{1}(t)\rangle,\cdots,\langle a_{i}(t)\rangle,\cdots)^{T},$ (23)
$\displaystyle\text{with~{}~{}}\langle
a(t)\rangle=\operatorname{tr}[a(t)\rho_{\text{S}+\text{PM}}(0)].$

Through the adjoint PMEOM, the dynamics for the amplitudes is determined by
the effective NHH

$H_{\text{eff},\text{S}+\text{PM}}=H_{\text{S}+\text{PM}}-i\sum_{i}\gamma_{i}a^{\dagger}_{i}a_{i}/2=\mathbf{v}^{\dagger}\mathbf{H}_{\text{eff},\text{S}+\text{PM}}\mathbf{v},$
(24)

where $\mathbf{H}_{\text{eff},\text{S}+\text{PM}}$ is a matrix representation
of the NHH. Specifically, one obtains

$\langle\dot{\mathbf{v}}(t)\rangle=i\mathbf{H}_{\text{eff},\text{S}+\text{PM}}\langle\mathbf{v}(t)\rangle,$
(25)

indicating that potential EPs can now be encoded in
$\mathbf{H}_{\text{eff},\text{S}+\text{PM}}$.

####  II.0.3 Example: two-mode system coupled to a Lorentzian environment

For a concrete example, we examine a two-mode system coupled to a Lorentzian
environment, as described in Eq. (14), assuming that

$\chi_{1,2}=\chi\text{~{}and~{}}\Omega_{1}=\Omega_{2}=\omega_{0}.$ (26)

In the Markovian limit ($\Lambda\rightarrow\infty$), the resulting system-only
effective non-Hermitian Hamiltonian within the rotating frame is given by:

$\mathbf{H}_{\text{eff},\text{S}}=\begin{pmatrix}0&\chi\\\
\chi&i\Gamma\end{pmatrix}.$ (27)

The corresponding eigenvalues are $(i\Gamma\pm\sqrt{4\chi^{2}-\Gamma^{2}})/2$,
indicating the presence of an EP2 if $|\chi|=\Gamma/2$, with the degenerate
eigenvalue $i\Gamma/2$. With a finite width $\Lambda$, the effective
Hamiltonian takes the form:

$\mathbf{H}_{\text{eff},\text{S}+\text{PM}}=\begin{pmatrix}0&\chi&0\\\
\chi&0&\sqrt{\frac{\Gamma\Lambda}{2}}\\\
0&\sqrt{\frac{\Gamma\Lambda}{2}}&i\Lambda\end{pmatrix}.$ (28)

By matching the coefficients of the characteristic polynomial, an EP3 is
identified with the following criteria:

$\left\\{|\chi|=\frac{\Lambda}{3\sqrt{3}},~{}\Gamma=\frac{16\Lambda}{27}\right\\},$
(29)

and the degenerate eigenvalue is $i\Lambda/3$. In other words, the EP can be
transformed from second to third order with the introduction of the structured
environmental characteristics.

This upgrade can lead to a further enhancement in the system sensitivity. For
instance, we introduce a perturbation $\epsilon>0$ to the coupling strength
$\chi\rightarrow\chi(1+\epsilon)$. For the scenario in the exact Markovian
limit, the eigenvalues in the vicinity of the EP2 are

$\left\\{i\frac{\Gamma}{2}-\frac{\Gamma\sqrt{\epsilon}}{\sqrt{2}}+O\left(\epsilon^{\frac{3}{2}}\right),i\frac{\Gamma}{2}+\frac{\Gamma\sqrt{\epsilon}}{\sqrt{2}}+O\left(\epsilon^{\frac{3}{2}}\right)\right\\}.$
(30)

In contrast to this case, for the scenario with a finite spectral width, the
eigenvalues in the vicinity of the EP3 take the form

$\displaystyle\Big{\\{}$ $\displaystyle
i\frac{\Lambda}{3}+x_{1}\Lambda\epsilon^{\frac{1}{3}}+O(\epsilon^{\frac{2}{3}}),i\frac{\Lambda}{3}+x_{2}\Lambda\epsilon^{\frac{1}{3}}+O(\epsilon^{\frac{2}{3}}),$
(31) $\displaystyle
i\frac{\Lambda}{3}+x_{3}\Lambda\epsilon^{\frac{1}{3}}+O(\epsilon^{\frac{2}{3}})\Big{\\}},$

where $x_{1}$, $x_{2}$, and $x_{3}$ are constants. We observe a change from a
square-root bifurcation for the Markovian EP2 to a cubic-root bifurcation for
the non-Markovian EP3 in response to the external perturbation, signifying the
enhancement of the system sensitivity to the external perturbation.

Intuitively, one may anticipate that higher-order EPs may emerge when
considering an environment with a more complicated spectral structure, which
requires more pseudomodes to capture the non-Markovian dynamics and increases
the dimension of the non-Markovian NHHs

##  III discussion

We have presented a general theory on characterizing non-Markovian EP based on
pseudomode mapping and hierarchical equations of motion. This approach is
applied to both the spin-boson model and linear bosonic systems, uncovering
the presence of a purely non-Markovian EP that corresponds to the transition
from Markovian to non-Markovian dynamics. We demonstrate that adjustments to
the EP criteria can be accomplished through the modification of the
environmental spectral profile. Additionally, the incorporation of non-
Markovian effects can increase the order of the EP by effectively increasing
the dimension of both the extended Liouvillian superoperator and NHHs. This
presents an innovative strategy for identifying higher-order EPs.

A direct extension involves exploring more realistic examples that do not rely
on the RWA and restore the detailed-balance condition [42]. In addition,
although this work focuses exclusively on a bosonic environment, the proposed
framework can be directly generalized to the scenarios with arbitrary
combinations of bosonic and fermionic baths [53, 55, 56, 57]. Moreover, beyond
the PMEOM and HEOM, our method of describing non-Markovian EPs via extended
Liouvillian superoperators can also be applied to other pertinent
methodologies, such as the dissipaton-embedded master equation [63, 64] and
reaction-coordinate mapping [65, 66, 67, 68].

Future work involves further generalizing the theory of non-Markovian EPs. An
intriguing possibility involves extending the hybrid-Liouvillian formalism
[20] to the non-Markovian domain by incorporating the postselection of quantum
trajectories. Additionally, it is worthwhile to delve into the potential
applications emerging from the intricate interplay between the (non-)Markovian
exceptional and diabolic points [12, 25] or the exotic topology and geometry
of the parameter space [69, 70, 71, 72, 73]. Such investigations could uncover
new aspects of non-Markovian EPs, enhancing our understanding of open quantum
systems embedded in environments with memory effects.

##  IV Methods

### The extended Liouvillian superoperators for HEOM

Here, we introduce the extended Liouvillian superoperators within the HEOM
formalism, considering both scenarios with and without the rotating wave
approximation (RWA). For the case without RWA, we decompose the correlation
function into exponential terms [as in Eq. (8)]:

$C(\tau)=\sum_{u=\mathbb{R},\mathbb{I}}(\delta_{u,\mathbb{R}}+i\delta_{u,\mathbb{I}})C^{u}(\tau)$
(32)

with $\mathbb{R}$ and $\mathbb{I}$ representing real and imaginary parts, and

$C^{u}(\tau)=\sum_{l}^{l_{\text{max}}}\xi_{l}^{u}\exp(-\chi_{l}^{u}\tau).$
(33)

Through iterative time differentiation of the exact dynamics in Eq. (3), the
HEOM can be expressed as a time-local differential equation within an expanded
space formed by the auxiliary density operators (ADOs): [57, 55, 58, 74, 75]

$\displaystyle\frac{d}{dt}\rho_{\text{S}+\text{ADO}}(t)=\mathcal{L}_{\text{S}+\text{ADO}}\left[\rho_{\text{S}+\text{ADO}}(t)\right].$
(34)

Here, we label the ADOs with a vector $\textbf{j}=[j_{m},\cdots,j_{1}]$,
linking each ADO to a specific exponential term in the correlation function.
The system reduced density operator and these ADOs can then be expressed as
$\rho^{(m)}_{\textbf{j}}(t)$, such that

$\rho_{\text{S}+\text{ADO}}=\begin{bmatrix}\rho^{(0)}_{\textbf{j}}(t)\\\
\rho^{(m>0)}_{\textbf{j}}(t)\\\ \vdots\end{bmatrix},$ (35)

where $m$ indicates the hierarchical level of the ADO with
$\rho_{\textbf{j}}^{(0)}(t)=\rho_{\text{S}}(t)$. Consequently, the system
reduced state can be obtained by performing a projector $\mathcal{P}$ on the
$\rho_{\text{S}+\text{ADO}}$:

$\displaystyle\begin{bmatrix}\rho_{\text{S}}(t)\\\ 0\\\ \vdots\\\
0\end{bmatrix}=\mathcal{P}[\rho_{\text{S}+\text{ADO}}(t)]=\begin{bmatrix}\mathbb{1}&\dots&0\\\
\vdots&\ddots&\vdots\\\
0&\dots&0\end{bmatrix}\begin{bmatrix}\rho^{(0)}_{\textbf{j}}(t)\\\
\rho^{(m>0)}_{\textbf{j}}(t)\\\ \vdots\end{bmatrix},$ (36)

where $\mathbb{1}$ is an identity matrix with the same dimension as the system
reduced density operator. Under this expression of ADOs, the corresponding
extended Liouvillian superoperator of HEOM can be written as [55, 57]

$\displaystyle\mathcal{L}_{\text{S}+\text{ADO}}[\rho^{(m)}_{\textbf{j}}(t)]=\mathcal{L}_{0}[\rho^{(m)}_{\textbf{j}}(t)]-\sum_{r=1}^{m}\chi_{j_{r}}\rho^{(m)}_{\textbf{j}}(t)$
(37)
$\displaystyle~{}~{}~{}~{}-i\sum_{j^{\prime}}\mathcal{A}_{j^{\prime}}[\rho^{(m+1)}_{\textbf{j}^{+}}(t)]-i\sum_{r=1}^{m}\mathcal{B}_{j_{r}}[\rho^{(m-1)}_{\textbf{j}_{r}^{-}}(t)],$

where a multi-index notation is used:
$\textbf{j}^{+}=[j^{\prime},j_{m},\cdots,j_{1}]$, and
$\textbf{j}_{r}^{-}=[j_{m},\cdots,j_{r+1},j_{r-1},\cdots,j_{1}]$, and
$\mathcal{L}_{0}[\cdot]=-iH_{\text{S}}^{\times}$. The system-environment
interaction is encoded in the superoperators $\hat{\mathcal{A}}_{j}$ and
$\hat{\mathcal{B}}_{j}$ that couple the $m$th-level bosonic ADOs to the
$(m+1)$th- and $(m-1)$th-levels, respectively. Their explicit expressions are
given by

$\displaystyle\mathcal{A}_{j}[\cdot]=Q^{\times}\text{~{}and~{}}\mathcal{B}_{j}[\cdot]=\delta_{u,\mathbb{R}}~{}\xi^{\mathbb{R}}_{l}~{}Q^{\times}+i\delta_{u,\mathbb{I}}~{}\xi^{\mathbb{I}}_{l}~{}Q^{\circ}.$
(38)

Let us now consider the case with RWA, where the system-environment
interaction is written as
$H_{\text{SE}}=\sum_{k}g_{k}(\tilde{Q}b_{k}^{\dagger}+\tilde{Q}^{\dagger}b_{k}).$
In this case, we separate the correlation function into the absorption
$(\nu=+)$ and emission $(\nu=-)$ components, namely
$C(t)=\sum_{\nu=\pm}C^{\nu}(t)$, where

$\displaystyle C^{+}(t)=\frac{1}{\pi}\int_{0}^{\infty}d\omega
J(\omega)n^{\textrm{eq}}(\omega)e^{i\omega\tau}$ (39)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}=\sum_{l}^{l_{\text{max}}}\xi_{l}^{\nu=+}\exp(-\chi_{l}^{\nu=+}\tau),$
and $\displaystyle C^{-}(t)=\frac{1}{\pi}\int_{0}^{\infty}d\omega
J(\omega)[n^{\textrm{eq}}(\omega)+1]e^{-i\omega\tau}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}=\sum_{l}^{l_{\text{max}}}\xi_{l}^{\nu=-}\exp(-\chi_{l}^{\nu=-}\tau).$

Here $n^{\textrm{eq}}(\omega)=\\{\exp[\omega/k_{\textrm{B}}T]-1\\}^{-1}$
represents the Bose–Einstein distribution with $k_{\textrm{B}}$ as the
Boltzmann constant and $T$ as the temperature, and $J(\omega)$ represents the
spectral density of the environment.

The Feynman-Vernon influence functional can be derived [55, 57]:

$\displaystyle\hat{\mathcal{F}}[\cdot]=-\sum_{\nu=\pm}\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}~{}\tilde{Q}^{\bar{\nu}}(t_{1})^{\times}$
(40)
$\displaystyle~{}~{}~{}\Big{\\{}C^{\nu}(t_{1}-t_{2})\tilde{Q}^{\nu}(t_{2})[\cdot]+C^{\bar{\nu}\ast}(t_{1}-t_{2})[\cdot]\tilde{Q}^{\nu}(t_{2})\Big{\\}},$

where $\tilde{Q}^{\nu=+}(t)=\tilde{Q}^{\dagger}(t)$ and
$\tilde{Q}^{\nu=-}(t)=\tilde{Q}(t)$. The corresponding extended Liouvillian
superoperator can then be expressed as

$\displaystyle\mathcal{L}_{\text{S}+\text{ADO}}^{\text{RWA}}[\rho^{(m)}_{\textbf{k}}(t)]=\mathcal{L}_{0}[\rho^{(m)}_{\textbf{k}}(t)]-\sum_{r=1}^{m}\chi_{k_{r}}\rho^{(m)}_{\textbf{k}}(t)$
(41)
$\displaystyle-i\sum_{k^{\prime}}\mathcal{A}_{k^{\prime}}[\rho^{(m+1)}_{\textbf{k}^{+}}(t)]-i\sum_{r=1}^{m}\mathcal{B}_{k_{r}}[\rho^{(m-1)}_{\textbf{k}_{r}^{-}}(t)],$

where the superoperators $\mathcal{A}_{k^{\prime}}$ and $\mathcal{B}_{k_{r}}$
are modified as follows

$\displaystyle\mathcal{A}_{k}[\cdot]=\tilde{Q}^{\bar{\nu}\times}~{}\text{and~{}}\mathcal{B}_{k}[\cdot]=\xi^{\nu}_{l}~{}\tilde{Q}^{\nu}[\cdot]-\xi^{\bar{\nu}\ast}_{l}[\cdot]~{}\tilde{Q}^{\nu}.$
(42)

## Acknowledgement

This work is supported by the National Center for Theoretical Sciences and
National Science and Technology Council, Taiwan, Grants No. NSTC
112-2123-M-006-001. N. L. is supported by the RIKEN Incentive Research Program
and by MEXT KAKENHI Grant Numbers JP24H00816, JP24H00820. A.M. was supported
by the Polish National Science Centre (NCN) under the Maestro Grant No.
DEC-2019/34/A/ST2/00081. F.N. is supported in part by: Nippon Telegraph and
Telephone Corporation (NTT) Research, the Japan Science and Technology Agency
(JST) [via the Quantum Leap Flagship Program (Q-LEAP), and the Moonshot R&D
Grant Number JPMJMS2061], the Asian Office of Aerospace Research and
Development (AOARD) (via Grant No. FA2386-20-1-4069), and the Office of Naval
Research (ONR) Global (via Grant No. N62909-23-1-2074).

## References

  * Heiss [2012] W. Heiss, The physics of exceptional points, J. Phys. A Math. Theor. 45, 444016 (2012).
  * El-Ganainy  _et al._ [2018] R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Non-Hermitian physics and $\mathcal{PT}$ symmetry, Nat. Phys. 14, 11 (2018).
  * Miri and Alu [2019] M.-A. Miri and A. Alu, Exceptional points in optics and photonics, Science 363, 7709 (2019).
  * Özdemir  _et al._ [2019] Ş. K. Özdemir, S. Rotter, _et al._ , Parity–time symmetry and exceptional points in photonics, Nat. Mater. 18, 783 (2019).
  * Bender and Boettcher [1998] C. M. Bender and S. Boettcher, Real Spectra in Non-Hermitian Hamiltonians Having $\mathcal{P}\mathcal{T}$ Symmetry, Phys. Rev. Lett. 80, 5243 (1998).
  * Bender [2007] C. M. Bender, Making sense of non-Hermitian Hamiltonians, Rep. Prog. Phys. 70, 947 (2007).
  * Moiseyev [2011] N. Moiseyev,  _Non-Hermitian quantum mechanics_ (Cambridge University Press, 2011).
  * Mostafazadeh [2002] A. Mostafazadeh, Pseudo-Hermiticity versus PT symmetry: the necessary condition for the reality of the spectrum of a non-Hermitian Hamiltonian, J. Math. Phys. 43, 205 (2002).
  * Peng  _et al._ [2014] B. Peng, Ş. Özdemir, S. Rotter, _et al._ , Loss-induced suppression and revival of lasing, Science 346, 328 (2014).
  * Jing  _et al._ [2014] H. Jing _et al._ , $\mathcal{PT}$-Symmetric Phonon Laser, Phys. Rev. Lett. 113, 053604 (2014).
  * Zhang  _et al._ [2018] J. Zhang _et al._ , A phonon laser operating at an exceptional point, Nat. Photon. 12, 479 (2018).
  * Arkhipov  _et al._ [2023] I. I. Arkhipov _et al._ , Dynamically crossing diabolic points while encircling exceptional curves: A programmable symmetric-asymmetric multimode switch, Nat. Commun. 14, 2076 (2023).
  * Chen  _et al._ [2017] W. Chen, Ş. Kaya Özdemir, G. Zhao, J. Wiersig, and L. Yang, Exceptional points enhance sensing in an optical microcavity, Nature 548, 192 (2017).
  * Hodaei  _et al._ [2017] H. Hodaei, A. U. Hassan, S. Wittek, H. Garcia-Gracia, R. El-Ganainy, D. N. Christodoulides, and M. Khajavikhan, Enhanced sensitivity at higher-order exceptional points, Nature 548, 187 (2017).
  * Rechtsman [2017] M. C. Rechtsman, Optical sensing gets exceptional, Nature 548, 161 (2017).
  * Wiersig [2020] J. Wiersig, Review of exceptional point-based sensors, Photon. Res. 8, 1457 (2020).
  * Kuo  _et al._ [2020] P.-C. Kuo, N. Lambert, A. Miranowicz, H.-B. Chen, G.-Y. Chen, Y.-N. Chen, and F. Nori, Collectively induced exceptional points of quantum emitters coupled to nanoparticle surface plasmons, Phys. Rev. A 101, 013814 (2020).
  * Minganti  _et al._ [2019] F. Minganti, A. Miranowicz, R. W. Chhajlany, and F. Nori, Quantum exceptional points of non-Hermitian Hamiltonians and Liouvillians: The effects of quantum jumps, Phys. Rev. A 100, 062131 (2019).
  * Arkhipov  _et al._ [2020a] I. I. Arkhipov, A. Miranowicz, F. Minganti, and F. Nori, Quantum and semiclassical exceptional points of a linear system of coupled cavities with losses and gain within the Scully-Lamb laser theory, Phys. Rev. A 101, 013812 (2020a).
  * Minganti  _et al._ [2020] F. Minganti, A. Miranowicz, R. W. Chhajlany, I. I. Arkhipov, and F. Nori, Hybrid-Liouvillian formalism connecting exceptional points of non-Hermitian Hamiltonians and Liouvillians via postselection of quantum trajectories, Phys. Rev. A 101, 062112 (2020).
  * Arkhipov  _et al._ [2020b] I. I. Arkhipov, A. Miranowicz, F. Minganti, and F. Nori, Liouvillian exceptional points of any order in dissipative linear bosonic systems: Coherence functions and switching between $\mathcal{PT}$ and anti-$\mathcal{PT}$ symmetries, Phys. Rev. A 102, 033715 (2020b).
  * Arkhipov  _et al._ [2021] I. I. Arkhipov, F. Minganti, A. Miranowicz, and F. Nori, Generating high-order quantum exceptional points in synthetic dimensions, Phys. Rev. A 104, 012205 (2021).
  * Khandelwal  _et al._ [2021] S. Khandelwal, N. Brunner, and G. Haack, Signatures of Liouvillian exceptional points in a quantum thermal machine, PRX Quantum 2, 040346 (2021).
  * Zhang  _et al._ [2022] J.-W. Zhang, J.-Q. Zhang, G.-Y. Ding, J.-C. Li, J.-T. Bu, B. Wang, L.-L. Yan, S.-L. Su, L. Chen, F. Nori, _et al._ , Dynamical control of quantum heat engines using exceptional points, Nat. Commun 13, 6225 (2022).
  * Perina Jr  _et al._ [2022] J. Perina Jr, A. Miranowicz, G. Chimczak, and A. Kowalewska-Kudlaszyk, Quantum Liouvillian exceptional and diabolical points for bosonic fields with quadratic Hamiltonians: The Heisenberg-Langevin equation approach, Quantum 6, 883 (2022).
  * Han  _et al._ [2023] P.-R. Han, F. Wu, X.-J. Huang, H.-Z. Wu, C.-L. Zou, W. Yi, M. Zhang, H. Li, K. Xu, D. Zheng, H. Fan, J. Wen, Z.-B. Yang, and S.-B. Zheng, Exceptional Entanglement Phenomena: Non-Hermiticity Meeting Nonclassicality, Phys. Rev. Lett. 131, 260201 (2023).
  * Zhou  _et al._ [2023] Y.-L. Zhou, X.-D. Yu, C.-W. Wu, X.-Q. Li, J. Zhang, W. Li, and P.-X. Chen, Accelerating relaxation through Liouvillian exceptional point, Phys. Rev. Res. 5, 043036 (2023).
  * Bu  _et al._ [2023] J.-T. Bu, J.-Q. Zhang, G.-Y. Ding, _et al._ , Enhancement of Quantum Heat Engine by Encircling a Liouvillian Exceptional Point, Phys. Rev. Lett. 130, 110402 (2023).
  * Khandelwal  _et al._ [2023] S. Khandelwal, W. Chen, K. W. Murch, and G. Haack, Chiral Bell-state transfer via dissipative Liouvillian dynamics, arXiv:2310.11381 (2023).
  * Abo  _et al._ [2024] S. Abo, P. Tulewicz, K. Bartkiewicz, Ş. K. Özdemir, and A. Miranowicz, Liouvillian Exceptional Points of Non-Hermitian Systems via Quantum Process Tomography, arXiv:2401.14993 (2024).
  * Laha  _et al._ [2024] A. Laha, A. Miranowicz, R. K. Varshney, and S. Ghosh, Correlated nonreciprocity around conjugate exceptional points, Phys. Rev. A 109, 033511 (2024).
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, 2002).
  * Cheung  _et al._ [2018] H. F. H. Cheung, Y. S. Patil, and M. Vengalattore, Emergent phases and critical behavior in a non-Markovian open quantum system, Phys. Rev. A 97, 052116 (2018).
  * Garmon  _et al._ [2021] S. Garmon, G. Ordonez, and N. Hatano, Anomalous-order exceptional point and non-Markovian Purcell effect at threshold in one-dimensional continuum systems, Phys. Rev. Res. 3, 033029 (2021).
  * Sergeev  _et al._ [2023] T. Sergeev, A. Zyablovsky, E. Andrianov, and Y. E. Lozovik, Signature of exceptional point phase transition in Hermitian systems, Quantum 7, 982 (2023).
  * Mouloudakis and Lambropoulos [2022] G. Mouloudakis and P. Lambropoulos, Coalescence of non-Markovian dissipation, quantum Zeno effect, and non-Hermitian physics in a simple realistic quantum system, Phys. Rev. A 106, 053709 (2022).
  * Garraway [1997] B. M. Garraway, Nonperturbative decay of an atomic system in a cavity, Phys. Rev. A 55, 2290 (1997).
  * Dalton  _et al._ [2001] B. J. Dalton, S. M. Barnett, and B. M. Garraway, Theory of pseudomodes in quantum optical processes, Phys. Rev. A 64, 053813 (2001).
  * Mazzola  _et al._ [2009] L. Mazzola, S. Maniscalco, J. Piilo, K.-A. Suominen, and B. M. Garraway, Pseudomodes as an effective description of memory: Non-Markovian dynamics of two-state systems in structured reservoirs, Phys. Rev. A 80, 012104 (2009).
  * Pleasance and Petruccione [2021] G. Pleasance and F. Petruccione, Pseudomode description of general open quantum system dynamics: non-perturbative master equation for the spin-boson model, arXiv:2108.05755 (2021).
  * Tamascelli  _et al._ [2018] D. Tamascelli, A. Smirne, S. F. Huelga, and M. B. Plenio, Nonperturbative Treatment of non-Markovian Dynamics of Open Quantum Systems, Phys. Rev. Lett. 120, 030402 (2018).
  * Lambert  _et al._ [2019] N. Lambert, S. Ahmed, M. Cirio, and F. Nori, Modelling the ultra-strongly coupled spin-boson model with unphysical modes, Nat. Commun. 10, 3721 (2019).
  * Pleasance  _et al._ [2020] G. Pleasance, B. M. Garraway, and F. Petruccione, Generalized theory of pseudomodes for exact descriptions of non-Markovian quantum processes, Phys. Rev. Res. 2, 043058 (2020).
  * Cirio  _et al._ [2023a] M. Cirio, N. Lambert, P. Liang, P.-C. Kuo, Y.-N. Chen, P. Menczel, K. Funo, and F. Nori, Pseudofermion method for the exact description of fermionic environments: From single-molecule electronics to the Kondo resonance, Phys. Rev. Res. 5, 033011 (2023a).
  * Luo  _et al._ [2023] S. Luo, N. Lambert, P. Liang, and M. Cirio, Quantum-Classical Decomposition of Gaussian Quantum Environments: A Stochastic Pseudomode Model, PRX Quantum 4, 030316 (2023).
  * Cirio  _et al._ [2023b] M. Cirio, S. Luo, P. Liang, F. Nori, and N. Lambert, Modeling the unphysical pseudomode model with physical ensembles: simulation, mitigation, and restructuring of non-Markovian quantum noise, arXiv:2311.15240 (2023b).
  * Menczel  _et al._ [2024] P. Menczel, K. Funo, M. Cirio, N. Lambert, and F. Nori, Non-Hermitian Pseudomodes for Strongly Coupled Open Quantum Systems: Unravelings, Correlations and Thermodynamics, arXiv:2401.11830 (2024).
  * Jin  _et al._ [2008] J. Jin, X. Zheng, and Y. Yan, Exact dynamics of dissipative electronic systems and quantum transport: Hierarchical equations of motion approach, J. Chem. Phys. 128 (2008).
  * Kreisbeck  _et al._ [2011] C. Kreisbeck, T. Kramer, M. Rodriguez, and B. Hein, High-performance solution of hierarchical equations of motion for studying energy transfer in light-harvesting complexes, Journal of Chemical Theory and Computation 7, 2166 (2011).
  * Ma  _et al._ [2012] J. Ma, Z. Sun, X. Wang, and F. Nori, Entanglement dynamics of two qubits in a common bath, Phys. Rev. A 85, 062323 (2012).
  * Li  _et al._ [2012] Z. H. Li, N. H. Tong, X. Zheng, D. Hou, J. H. Wei, J. Hu, and Y. J. Yan, Hierarchical Liouville-space approach for accurate and universal characterization of quantum impurity systems, Phys. Rev. Lett. 109, 266403 (2012).
  * Dunn  _et al._ [2019] I. S. Dunn, R. Tempelaar, and D. R. Reichman, Removing instabilities in the hierarchical equations of motion: Exact and approximate projection approaches, J. Chem. Phys. 150 (2019).
  * Tanimura [2020] Y. Tanimura, Numerically “exact” approach to open quantum dynamics: The hierarchical equations of motion (HEOM), J. Chem. Phys. 153, 020901 (2020).
  * Ikeda and Scholes [2020] T. Ikeda and G. D. Scholes, Generalization of the hierarchical equations of motion theory for efficient calculations with arbitrary correlation functions, J. Chem. Phys. 152 (2020).
  * Cirio  _et al._ [2022] M. Cirio, P. C. Kuo, Y. N. Chen, F. Nori, and N. Lambert, Canonical derivation of the fermionic influence superoperator, Phys. Rev. B 105, 035121 (2022).
  * Lambert  _et al._ [2023] N. Lambert, T. Raheja, S. Cross, P. Menczel, S. Ahmed, A. Pitchford, D. Burgarth, and F. Nori, QuTiP-BoFiN: A bosonic and fermionic numerical hierarchical-equations-of-motion library with applications in light-harvesting, quantum control, and single-molecule electronics, Phys. Rev. Res. 5, 013181 (2023).
  * Huang  _et al._ [2023] Y.-T. Huang, P.-C. Kuo, N. Lambert, M. Cirio, S. Cross, S.-L. Yang, F. Nori, and Y.-N. Chen, An efficient Julia framework for hierarchical equations of motion in open quantum systems, Commun. Phys. 6, 313 (2023).
  * Kuo  _et al._ [2023] P.-C. Kuo, N. Lambert, M. Cirio, Y.-T. Huang, F. Nori, and Y.-N. Chen, Kondo QED: The Kondo effect and photon trapping in a two-impurity Anderson model ultrastrongly coupled to light, Phys. Rev. Res. 5, 043177 (2023).
  * Nakamura and Tanimura [2021] K. Nakamura and Y. Tanimura, Optical response of laser-driven charge-transfer complex described by Holstein–Hubbard model coupled to heat baths: Hierarchical equations of motion approach, J. Chem. Phys. 155, 064106 (2021).
  * Kato and Tanimura [2016] A. Kato and Y. Tanimura, Quantum heat current under non-perturbative and non-Markovian conditions: Applications to heat machines, J. Chem. Phys 145, 224105 (2016).
  * Breuer  _et al._ [2009] H.-P. Breuer, E.-M. Laine, and J. Piilo, Measure for the Degree of Non-Markovian Behavior of Quantum Processes in Open Systems, Phys. Rev. Lett. 103, 210401 (2009).
  * Rivas  _et al._ [2010] A. Rivas, S. F. Huelga, and M. B. Plenio, Entanglement and Non-Markovianity of Quantum Evolutions, Phys. Rev. Lett. 105, 050403 (2010).
  * Yan [2014] Y. Yan, Theory of open quantum systems with bath of electrons and phonons and spins: Many-dissipaton density matrixes approach, J. Chem. Phys. 140 (2014).
  * Wang and Yan [2022] Y. Wang and Y. Yan, Quantum mechanics of open systems: Dissipaton theories, J. Chem. Phys. 157 (2022).
  * Iles-Smith  _et al._ [2014] J. Iles-Smith, N. Lambert, and A. Nazir, Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems, Phys. Rev. A 90, 032114 (2014).
  * Iles-Smith  _et al._ [2016] J. Iles-Smith, A. G. Dijkstra, N. Lambert, and A. Nazir, Energy transfer in structured and unstructured environments: Master equations beyond the Born-Markov approximations, J. Chem. Phys. 144 (2016).
  * Restrepo  _et al._ [2019] S. Restrepo, S. Böhling, J. Cerrillo, and G. Schaller, Electron pumping in the strong coupling and non-Markovian regime: A reaction coordinate mapping approach, Phys. Rev. B 100, 035109 (2019).
  * Anto-Sztrikacs and Segal [2021] N. Anto-Sztrikacs and D. Segal, Capturing non-Markovian dynamics with the reaction coordinate method, Phys. Rev. A 104, 052617 (2021).
  * Ding  _et al._ [2022] K. Ding, C. Fang, and G. Ma, Non-Hermitian topology and exceptional-point geometries, Nat. Rev. Phys. 4, 745 (2022).
  * Bergholtz  _et al._ [2021] E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Exceptional topology of non-Hermitian systems, Rev. Mod. Phys. 93, 015005 (2021).
  * Ju  _et al._ [2022] C.-Y. Ju, A. Miranowicz, F. Minganti, C.-T. Chan, G.-Y. Chen, and F. Nori, Einstein’s quantum elevator: Hermitization of non-Hermitian Hamiltonians via a generalized vielbein formalism, Phys. Rev. Res. 4, 023070 (2022).
  * Ju  _et al._ [2024a] C.-Y. Ju, A. Miranowicz, Y.-N. Chen, G.-Y. Chen, and F. Nori, Emergent parallel transport and curvature in Hermitian and non-Hermitian quantum mechanics, Quantum 8, 1277 (2024a).
  * Ju  _et al._ [2024b] C.-Y. Ju, F.-H. Huang, and G.-Y. Chen, Event-Horizon-Like Singularities and Quantum Phase Transitions, arXiv:2403.16503 (2024b).
  * Cainelli  _et al._ [2022] M. Cainelli, R. Borrelli, and Y. Tanimura, Effect of mixed Frenkel and charge transfer states in time-gated fluorescence spectra of perylene bisimides H-aggregates: Hierarchical equations of motion approach, J. Chem. Phys 157, 084103 (2022).
  * Ikeda and Tanimura [2017] T. Ikeda and Y. Tanimura, Probing photoisomerization processes by means of multi-dimensional electronic spectroscopy: The multi-state quantum hierarchical Fokker-Planck equation approach, J. Chem. Phys. 147, 014102 (2017).