Quantum Mpemba Effect from Non-Normal Dynamics

1. Introduction

The thermal relaxation of physical systems is commonly assumed to follow an intuitive rule: states closer to equilibrium should relax faster than those further away. Rooted in classical thermodynamics, this logic suggests that, for instance, a hotter object should cool more slowly than a colder one when coupled to a thermal bath. However, this expectation can break down. One striking example is the Mpemba effect (ME), where a system initially at a higher temperature cools faster than a similar system at a lower temperature [1,2,3,4,5,6,7,8,9,10,11,12]. While its classical version remains debated and lacks a universally accepted explanation [4,12], the Mpemba effect has sparked renewed interest in the context of quantum systems [13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46], where it can be explored in cleaner, more controllable settings [46].

In the quantum domain, Mpemba-like phenomena have been studied in both open and closed systems. Broadly, the effect refers to cases where a state initially further from equilibrium relaxes faster than one closer to it. Several mechanisms underlying this behavior have been identified. The strong ME arises when certain initial states have vanishing overlap with the slowest-decaying eigenmodes of the Liouvillian, leading to enhanced asymptotic decay [14,16,17,19,21,36,46]. Another pathway involves initial system-environment correlations, such as entanglement, which effectively alter the open system’s initial conditions, even under Markovian assumptions [43]. Non-Markovian dynamics offer yet another route, where memory effects can produce anomalous relaxation [9,45], including cases where a state reaches equilibrium within the environment’s memory time– a phenomenon referred to as the extreme quantum Mpemba effect [45]. Beyond open quantum systems, Mpemba-like phenomena have also been observed in closed quantum systems evolving unitarily after a quantum quench [18,23,24,26,27,35]. In such cases, anomalous relaxation behavior has been identified in the growth of entanglement entropy, where a system starting further from equilibrium can approach a thermal-like state faster than one closer to equilibrium. These mechanisms illustrate that the ME in quantum systems can emerge from a variety of physical mechanisms.

In this work, we introduce a new and conceptually distinct mechanism that can give rise to a Mpemba-like effect in fully Markovian quantum dynamics, without relying on long-time decay rate hierarchies, non-Markovianity, or entanglement. The key ingredient is the non-normality of the Liouvillian superoperator governing the system’s evolution. Non-normal operators –those that do not commute with their adjoint – have non-orthogonal eigenmodes, leading to transient interference effects between decaying components [47,48,49]. As a consequence, certain initial conditions can exhibit anomalous early-time behavior, including slow early relaxation, transient freezing dynamics or delayed relaxation [48,49,50,51,52,53,54,55,56,57], even when the long-time asymptotic decay is exponential with the decay rate of the slowest decaying mode. We refer to this phenomenon as the non-normal quantum ME. To illustrate this effect, we study a model system drawn from waveguide quantum electrodynamics (QED), in which quantum emitters (modeled as two-level systems) couple to the photonic modes of a one-dimensional waveguide [58,59,60,61]. Such setups facilitate long-range coherent and incoherent interactions mediated by long-lived propagating photons and naturally give rise to non-normality in the Liouvillian. We analyze the relaxation dynamics of different initial states and identify regimes where certain configurations closer to equilibrium, evolving under a memoryless, Lindblad-type dynamics. exhibit slower short-time relaxation due to destructive modal interference, leading to a pronounced ME. Our results enrich the taxonomy of quantum ME by identifying transient interference – rather than asymptotic behavior – as a key source of anomalous relaxation. This opens new directions for controlling nonequilibrium processes in engineered quantum systems, particularly in photonic or circuit QED platforms where non-normality can be designed and exploited for quantum state control or thermodynamic manipulation.

2. Mpemba Effect from Non-Normal Dynamics in Markovian Open Quantum Systems

Let us consider a Markovian open quantum system whose state at time t is described by the reduced density operator $\rho \left(t\right)$. Its time evolution is governed by the Liouvillian superoperator $\mathcal{L}$ in the standard Lindblad form [62], $(d\rho /dt)=\mathcal{L}\rho \left(t\right)$, where

$$\mathcal{L}=-i[H,\rho ]+\sum _k\left(L_k\rho L_k^{†}-{\displaystyle \frac{1}{2}}\left\{L_k^{†}{L}_k,\rho \right\}\right).$$

(1)

In the above equation, H is the system Hamiltonian describing the intrinsic unitary evolution of the system whereas $L_k$ are a set of jump operators describing the dissipative part of the dynamics. The Liouvillian superoperator can be equivalently written as

$$\mathcal{L}=-i\left(H_{eff}\rho -\rho H_{eff}^{†}\right)+\sum _k{L}_k\rho L_k^{†}.$$

(2)

where

$$H_{eff}=H-{\displaystyle \frac{i}{2}}\sum _k{L}_k^{†}{L}_k$$

(3)

is an effective non-Hermitian Hamiltonian and the last term on the right hand side of Equation (2), ${\sum }_k{L}_k\rho L_k^{†}$, describes quantum jumps. We assume that the Liouvillian $\mathcal{L}$ is diagonalizable, and let us denote by ${\lambda }_{\alpha }$ the eigenvalues of $\mathcal{L}$, and by ${\rho }_{\alpha }$ and $l_{\alpha }$ the corresponding right and left eigenvectors, respectively. Since any eigenvalue has a nonpositive real part, they can be sorted in descending order, $0={\lambda }_0>\mathrm{Re}\left({\lambda }_1\right)\ge \mathrm{Re}\left({\lambda }_2\right)\ge ...$. We assume that the zero eigenvalue ${\lambda }_0=0$ is not degenerate, indicating that the steady state is unique and corresponds to the state ${\rho }_0$. The Liouvillian gap is defined as $\Delta =-\mathrm{Re}\left({\lambda }_1\right)$. The right and left eigenvectors of $\mathcal{L}$ form an orthogonal basis, namely $\langle l_{\alpha }|{\rho }_{\beta }\rangle ={\delta }_{\alpha ,\beta }$, where for two operators A and B the scalar product $\langle A|B\rangle$ is defined as $\langle A|B\rangle =\mathrm{Tr}\left(A^{†}B\right)$ (see e.g., [53] ). For an initial state ${\rho }_i$, the general evolution of the density operator $\rho \left(t\right)$ can be written as

$$\rho \left(t\right)=exp\left(\mathcal{L}t\right){\rho }_i={\rho }_0+\sum _{\alpha \ge 1}{c}_{\alpha }{\rho }_{\alpha }exp\left({\lambda }_{\alpha }t\right)$$

(4)

where $c_{\alpha }=\langle l_{\alpha }|{\rho }_i\rangle =\mathrm{Tr}\left(l_{\alpha }^{†}{\rho }_i\right)$. To quantify the distance $d\left(t\right)$ of the system from its equilibrium state ${\rho }_0$ during the relaxation process, different measures can be assumed (see e.g. [46,63]), such as the trace distance $d\left(t\right)=d_T\left(t\right)=(1/2)\mathrm{Tr}\sqrt{A^{†}A}$ and the Frobenius (Hilbert-Schmidt) distance $d\left(t\right)=d_F\left(t\right)=\sqrt{\mathrm{Tr}\left(A^{†}A\right)}$, where $A=\rho \left(t\right)-{\rho }_0$. For example, assuming the Frobenius distance $d\left(t\right)=d_F\left(t\right)$, using Equation (4) one has

$$d^2\left(t\right)=\sum _{\alpha ,\beta =1,2,3,...}{c}_{\alpha }^{*}{c}_{\beta }exp\left[({\lambda }_{\alpha }^{*}+{\lambda }_{\beta })t\right]\mathrm{Tr}\left({\rho }_{\alpha }^{†}{\rho }_{\beta }\right).$$

(5)

The ME arises when comparing the relaxation dynamics of two initial states I and II, described by the density operators ${\rho }_i^{\left(I\right)}$ and ${\rho }_i^{\left(II\right)}$. Assuming that state I is initially further away from equilibrium, i.e., $d^{\left(I\right)}\left(0\right)>d^{\left(II\right)}\left(0\right)$, the ME occurs whenever, after an initial transient, one has $d^{\left(I\right)}\left(t\right)<d^{\left(II\right)}\left(t\right)$ for $t>t_f$, $t_f$ being a suitably large time. Since the Liouvillian spectral gap $\Delta$ determines the decay rate of the slowest relaxation mode, it is naturally expected that its inverse, $1/\Delta$, sets the characteristic timescale for relaxation. This expectation holds when the Liouvillian superoperator $\mathcal{L}$ is normal, i.e., when it commutes with its adjoint: $[\mathcal{L},{\mathcal{L}}^{†}]=0$. In this case, the right and left eigenvectors of $\mathcal{L}$ do coincide, and thus right eigenvectors corresponding to distinct eigenvalues are orthogonal, meaning $\langle {\rho }_{\alpha }|{\rho }_{\beta }\rangle =\mathrm{Tr}\left({\rho }_{\alpha }^{†}{\rho }_{\beta }\right)=0$ for $\alpha \ne \beta$. As a result, the interference terms with $\alpha \ne \beta$ on the right-hand side of Equation (5) vanish, yielding

$$d^2\left(t\right)=\sum _{\alpha ,1,2,3,...}{\left|c_{\alpha }\right|}^{2}exp\left(2\mathrm{Re}\left({\lambda }_{\alpha }\right)t\right).$$

(6)

This result clearly demonstrates that the inverse of the spectral gap, $1/\Delta =1/\left|\mathrm{Re}\right({\lambda }_1\left)\right|$, indeed sets an upper bound on the relaxation time toward equilibrium. More precisely, the actual relaxation time is determined by the largest value of $1/\left|\mathrm{Re}\right({\lambda }_{\alpha }\left)\right|$ for which the corresponding spectral coefficient $c_{\alpha }$ is non-zero. When the initial state displays a non-vanishing spectral coefficient $c_1\ne 0$, the relaxation time is given by $1/\Delta$, i.e., it is limited by the spectral gap. When the coefficient $c_1$ vanishes, the relaxation proceeds faster, governed by the next slowest decay rate, $-\mathrm{Re}\left({\lambda }_2\right)$. This is the essence of the strong Mpemba effect (see, e.g., [14,17,21,36,46]): if state I, initially farther from equilibrium, does not excite the slowest decaying mode (i.e., $c_1=0$ for state I), it relaxes faster at long times than state II, which does have a non-zero overlap with this slowest mode (i.e., $c_1\ne 0$ for state II). This situation is schematically illustrated in Figure 1a.

However, when $\mathcal{L}$ is not a normal operator, i.e., $\mathcal{L}{\mathcal{L}}^{†}\ne {\mathcal{L}}^{†}\mathcal{L}$, the inverse of the spectral gap $1/\Delta$does not set an upper limit to the characteristic relaxation time toward equilibrium, which can be extremely longer than $1/\Delta$. Indeed, several studies have emphasized that the actual relaxation time is not necessarily bounded by the Liouvillian gap [53,54,55,56], although the latter does determine the asymptotic decay rate in the long-time limit [14]. This implies that, in the transient regime, the system can exhibit decay rates that are slower than the Liouvillian gap would suggest. Such anomalously slow relaxation, which arises from the interference terms on the right hand side of Equation (5), can give rise to Mpemba-like effects without requiring accelerated thermalization for one of the two initial states. This type of ME, schematically illustrated in Figure 1b, is fundamentally different from the strong ME and may be referred to as the non-normal Mpemba effect, as it originates from the non-normal nature of the Liouvillian superoperator.

3. Non-Normal Mpemba Effect in Waveguide Quantum Electrodynamics

To illustrate the non-normal Mpemba effect within a key quantum platform, we consider the relaxation dynamics of quantum emitters in the context of waveguide quantum electrodynamics (waveguide QED). In these systems, quantum emitters such as atoms, quantum dots, or superconducting qubits are coupled to one-dimensional photonic waveguides, enabling strong and tunable light-matter interactions [58,59,60,61]. Waveguide QED platforms are particularly well-suited for exploring open quantum system dynamics, as they naturally give rise to structured reservoirs, long-range photon-mediated interactions, and non-Markovian effects. These features naturally lead to a non-normal Liouvillian superoperator, where the relaxation dynamics are governed not only by the spectral gap but also by the non-orthogonality and interference of decay modes. This makes waveguide QED an ideal setting for observing the non-normal ME, where anomalously slow relaxation emerges even in the absence of accelerated decay pathways. Waveguide QED systems offer a broad range of physical implementations [60], including atoms coupled to photonic waveguides [64,65], photonic crystals [66], quantum dots embedded in nanophotonic structures [67], and superconducting circuits [68,69,70,71]. These platforms provide fertile ground for investigating collective phenomena and complex thermalization behavior driven by coherent and incoherent photon-mediated interactions over long distances.

We consider the dissipative dynamics of a set of N point-like quantum emitters, described as two-level (spin) systems, coupled to a one-dimensional bosonic waveguide, as schematically illustrated in Figure 1c. In the Born-Markov approximation and neglecting retardation effects, the evolution of the reduced density matrix $\rho \left(t\right)$ of the quantum emitters is governed by the Lindblad master equation [58,59,62]

$${\displaystyle \frac{d\rho }{dt}}=\mathcal{L}\rho =-i(H_{eff}\rho -\rho H_{eff}^{†})+{\gamma }_R{L}_R\rho L_R^{†}+{\gamma }_L{L}_L\rho L_L^{†}$$

(7)

where

$$\begin{array}{ccc}\hfill H_{eff}& =& \sum _{n=1}^N\left({\omega }_n-{\omega }_0-i{\displaystyle \frac{{\gamma }_L+{\gamma }_R}{2}}\right){\sigma }_n^{†}{\sigma }_n\hfill \\ & -& i{\gamma }_L\sum _{n<l}{\sigma }_n^{†}{\sigma }_{l}exp\left[ik(x_l-x_n)\right]-i{\gamma }_R\sum _{n>l}{\sigma }_n^{†}{\sigma }_{l}exp\left[ik(x_n-x_l)\right]\hfill \end{array}$$

(8)

is the effective non-Hermitian Hamiltonian and

$$L_L=\sum _{n=1}^N{\sigma }_{n}exp\left(ik{x}_n\right),L_R=\sum _{n=1}^N{\sigma }_{n}exp(-ik{x}_n)$$

(9)

are the jump operators associated to decay into the left- and right-propagating modes of the waveguide. In the above equations, $x_n$ is the spatial position of the n-th quantum emitter in the one-dimensional waveguide, with $x_1<x_2<...<x_{N-1}<x_N$ [Figure 1c], ${\omega }_n$ is its resonance frequency, close to the frequency ${\omega }_0=kc$ of the excited left/right propagating photonic modes in the waveguide, ${\sigma }_n^{†}={|e\rangle }_n\langle g|$ and ${\sigma }_n={|g\rangle }_n\langle e|$ are the atom transition operators with ground state $|g\rangle$ and excited state $|e\rangle$, and ${\gamma }_{L,R}$ are the emission rates of a single quantum emitter into the right (R) and left (L) propagating modes of the waveguide. A chiral waveguide corresponds to ${\gamma }_R\ne {\gamma }_L$[58]. The Hilbert space of the N quantum emitters is spanned by the $2^N$ states $|q_1\rangle \otimes |q_2\rangle \otimes ...|q_N\rangle$, where $|q_n\rangle$ is either $|g\rangle$ or $|e\rangle$, corresponding to the n-th atom being in the ground or excited state. In this basis, the density operator $\rho \left(t\right)$ is described by a $2^N\times 2^N$ matrix, and the superoperator $\mathcal{L}$ by a $2^{2N}\times 2^{2N}$ matrix. We assume that there are not atom-photon bound states, so that any initial atomic excitations irreversibly decay into the ground states, corresponding to the unique stationary state ${\rho }_0=|0\rangle \langle 0|$, where $|0\rangle \equiv |g\rangle \otimes |g\rangle \otimes ...\otimes |g\rangle$ is the pure state with all N atoms in the ground state.

To unveil the emergence of the non-normal ME, let us restrict our attention considering the states corresponding to up to a single excitation in the system [60], so that the Hilbert space is greatly reduced and comprises solely the $(N+1)$ states $|0\rangle$ (the ground equilibrium state) and $|\alpha \rangle \equiv {\sigma }_{\alpha }^{†}|0\rangle$ (the $\alpha$-th atom is in the excited state, while all other atoms are in the ground state, with $\alpha =1,2,...,N$). In this case, the density matrix $\rho \left(t\right)$ is a $(N+1)\times (N+1)$ matrix, with elements

$${\rho }_{0,0}=\langle 0\left|\rho \right|0\rangle ,{\rho }_{0,\alpha }=\langle 0\left|\rho \right|\alpha \rangle ={\rho }_{\alpha ,0}^{*},{\rho }_{\alpha ,\beta }=\langle \alpha \left|\rho \right|\beta \rangle ={\rho }_{\beta ,\alpha }^{*}$$

(10)

($\alpha ,\beta =1,2,...,N$). The explicit evolution equations of the density matrix elements, as derived from the master Equation (7), are given in the Appendix A. A special yet broad subclass of initial states relevant to analytically describe the non-normal ME is provided by the mixed states

$${\rho }_i=\rho (t=0)=p|{\psi }_0\rangle \langle {\psi }_0|+(1-p)|0\rangle \langle 0|$$

(11)

where $|{\psi }_0\rangle ={\sum }_{\alpha =1}^N{c}_{\alpha }^{\left(0\right)}|\alpha \rangle$, ${\sum }_{\alpha =1}^N{\left|c_{\alpha }^{\left(0\right)}\right|}^2=1$ and $0\le p\le 1$. For such a subclass of initial states of the system, the evolved density operator $\rho \left(t\right)$, as obtained from Equations (A3)–(A6) given in the Appendix A, can be written in the form

$$\rho \left(t\right)=p|\psi \left(t\right)\rangle \langle \psi \left(t\right)|+{\rho }_{0,0}\left(t\right)|0\rangle \langle 0|.$$

(12)

In the above equation, we have set

$$|\psi \left(t\right)\rangle =exp(-i{H}_{eff}t)|{\psi }_0\rangle =\sum _{\alpha =1}^N{c}_{\alpha }\left(t\right)|\alpha \rangle ,{\rho }_{0,0}\left(t\right)=1-p\sum _{\alpha }{\left|c_{\alpha }\left(t\right)\right|}^2.$$

(13)

with $c_{\alpha }(t=0)=c_{\alpha }^{\left(0\right)}$. Basically, the state $\left|\psi \right(t)\rangle$ evolves under the effect of the effective non-Hermitian Hamiltonian solely (quantum jumps do not affect its evolution), and the amplitudes $c_{\alpha }\left(t\right)$ are thus found by solving the linear system of coupled equations

$$\begin{array}{ccc}\hfill i{\displaystyle \frac{d{c}_{\alpha }}{dt}}& =& \sum _{\beta =1}^N\langle \alpha |H_{eff}|\beta \rangle c_{\beta }\left(t\right)=\left({\omega }_{\alpha }-{\omega }_0-i{\displaystyle \frac{{\gamma }_L+{\gamma }_R}{2}}\right)c_{\alpha }+\hfill \\ & -& i{\gamma }_L\sum _{\beta >\alpha }exp\left[ik(x_{\beta }-x_{\alpha })\right]c_{\beta }-i{\gamma }_R\sum _{\beta <\alpha }exp\left[ik(x_{\alpha }-x_{\beta })\right]c_{\beta }.\hfill \end{array}$$

(14)

The distance from equilibrium $d\left(t\right)$, as measured by the Frobenius distance $d\left(t\right)=d_F\left(t\right)$, can be readily calculated from Equation (12) and reads

$$d\left(t\right)=\sqrt{2}p{P}_s\left(t\right)$$

(15)

where we has set

$$P_s\left(t\right)=\sum _{\alpha =1}^N{\left|c_{\alpha }\left(t\right)\right|}^2.$$

(16)

The quantity $P_s\left(t\right)$ can be viewed as the survival probability of atomic excitation, i.e., the probability that no photon has been irreversibly emitted into the waveguide at time t[51]. Since we assume that there are not atom-photon bound states, clearly $P_s\left(t\right)\to 0$ as $t\to \infty$, corresponding to $d\left(t\right)\to 0$ as $t\to \infty$. However, the decay behavior of $P_s\left(t\right)$, and thus of $d\left(t\right)$, can display anomalous features related to the non-normal nature of the effective non-Hermitian Hamiltonian $H_{eff}$ (see e.g., [51,52] ). Departure from non-normality of the Hamiltonian $H_{eff}$ can be measured by the Henrici parameter [50,72]

$$\nu ={\left(\sum _{\alpha ,\beta =1}^N{\left|{\left(H_{eff}\right)}_{\alpha ,\beta }\right|}^2-\sum _{\alpha =1}^N{\left|{\xi }_{\alpha }\right|}^2\right)}^{1/2}$$

(17)

where ${\left(H_{eff}\right)}_{\alpha ,\beta }=\langle \alpha |H_{eff}|\beta \rangle$ and ${\xi }_{\alpha }$ are the N eigenvalues of the matrix${\left(H_{eff}\right)}_{\alpha ,\beta }$. The Henrici parameter $\nu$ vanishes if and only if the matrix ${\left(H_{eff}\right)}_{\alpha ,\beta }$ is normal. Effective non-Hermitian Hamiltonians in chiral waveguides, where $\gamma \_R\ne \gamma \_L$, typically exhibit a high degree of non-normality. However, some level of non-normality can also arise in non-chiral waveguides. As an example, Figure 2a shows the non-normality – quantified by the Henrici parameter $\nu$ – of the effective Hamiltonian $H_{eff}$ as a function of the ratio ${\gamma }_R/{\gamma }_L$, for a system of $N=8$ emitters with identical resonance frequencies ${\omega }_n={\omega }_0$. The emitters are equally spaced along the waveguide with spacing $d=x_{n+1}-x_n$ chosen such that $kd=2\pi$. For this value of $kd$, no atom-photon bound states are present, meaning all eigenvalues ${\xi }_{\alpha }$ have strictly negative imaginary parts – except in the special case of a non-chiral waveguide (${\gamma }_L={\gamma }_R$), where bound states do emerge. Note that the largest value $\nu$ of non-normality is obtained in the chiral cascaded configuration [58]${\gamma }_R/{\gamma }_L=0$, corresponding to unidirectional photon emission, while a normal matrix ($\nu =0$) is found in the limit of a non-chiral waveguide ${\gamma }_R={\gamma }_L$. However, for different values of the phase $kd$ non-normality is obtained for a non-chiral waveguide ${\gamma }_L={\gamma }_R$ as well.

Non-normality is responsible for a variety of transient anomalous relaxation of $P_s\left(t\right)$ toward zero, such a quiescent dynamics, corresponding to an almost non-decaying probability $P_s\left(t\right)\simeq 1$ for a long time interval $0<t<\tau$, after which the decay starts with an abrupt drop of the survival probability [52]. For a given initial condition $\left\{c_{\alpha }^{\left(0\right)}\right\}$, the upper bound for $P_s\left(t\right)$ is given by $P_s\left(t\right)\le {\sigma }_m\left(t\right)$, where ${\sigma }_m\left(t\right)$ is the largest eigenvalue of the matrix ${\mathcal{A}}^{†}\left(t\right)\mathcal{A}\left(t\right)$ with $\mathcal{A}\left(t\right)=exp(-i{H}_{eff}t)$[51]. For a fixed time $t>0$, the largest value ${\sigma }_m\left(t\right)$ is assumed for the initial condition $\left\{c_{\alpha }^{\left(0\right)}\right\}$ of the system which is the eigenvector of ${\mathcal{A}}^{†}\left(t\right)\mathcal{A}\left(t\right)$ with eigenvalue ${\sigma }_m\left(t\right)$. The resilience time $\tau$ can be defined as the largest time $t=\tau$ such that the upper bound ${\sigma }_m$ of the survival probability remains larger than a reference value $P_b\sim 1$ in the entire time interval $(0,\tau )$. An example of resilience dynamics in a chiral waveguide, with $N=8$ emitters, phase $kd=2\pi$ and chiral parameter ${\gamma }_R/{\gamma }_L=0.1$, is shown in Figure 2b. On the other hand, for other initial states, such as when $c_{\alpha }\left(0\right)={\delta }_{\alpha ,1}$ (this condition corresponds to excitation of the atom at the far left edge of the chain), resilience dynamics is not observed and the decay is almost exponential (as ${\gamma }_R/{\gamma }_L\to 0$) with a decay rate ${\gamma }_L$. Such a distinct decay behaviors depending on initial state preparation can yield a pronounced non-normal ME, as illustrated in Figure 3. The figure illustrates the numerically-computed behaviors of the distance $d\left(t\right)$, either the trace distance $d\left(t\right)=d_T\left(t\right)$ [Figure 3a] and Frobenius distance $d=d_F\left(t\right)$ [Figure 3b], corresponding to three different initial states (11) in a highly-chiral waveguide (${\gamma }_L=1$, ${\gamma }_R=0.1$). State I is defined by Equation (11) with $p=1$ and $c_{\alpha }={\delta }_{\alpha ,1}$. This is a pure state with the atom at the far left end of the chain ($\alpha =1$) in the excited state and all other atoms in the ground state. This state displays a near exponential relaxation to equilibrium with a relaxation rate $\simeq {\gamma }_L$; the weak oscillations observed in the decay on a log scale, shown in the insets, arise from weak interference effects originating from a non-vanishing yet small value ${\gamma }_R$.

State II is defined by Equation (11) with $p=0.5$ and with $c_{\alpha }$ given by the distribution shown in the lower inset of Figure 2b. This is a mixed state where with 50% probability the system is already in the equilibrium state (all atoms are in the ground state) and with 50% probability the atoms are prepared in the single-excitation entangled state of Figure 2b. This state displays resilience to decay for a time interval $\tau \sim 16$, leading to a pronounced ME when compared to the relaxation of state I. Finally, state III is defined by Equation (11) with $p=0.5$ and $c_{\alpha }={\delta }_{\alpha ,8}$. This is a mixed state where with 50% probability the system is already in the equilibrium state (all atoms are in the ground state) and with 50% probability the atom at the far right end of the chain ($\alpha =N=8$) is in the excited state. Compared to state II, the initial preparation of atoms in state III is clearly much more feasible from an experimental standpoint. Non-normal ME is also observed when comparing the relaxation dynamics of states I and III. Finally, we note that in the limiting case of cascaded non-chiral waveguides, where ${\gamma }_R/{\gamma }_I=0$[58], the effective Hamiltonian $H_{eff}$ becomes non-diagonalizable and exhibits a higher-order exceptional point when ${\omega }_n={\omega }_0$. In this regime, the expansion analysis must include generalized eigenvectors due to the defective nature of both $\mathcal{L}$ and $H_{eff}$ (see, e.g., [51,52]). Nevertheless, non-normality – the fundamental mechanism underlying the non-normal memory effect – remains present. Numerical simulations show a relaxation behavior qualitatively similar to that in Figure 3, except that the weak oscillations in the decay dynamics, stemming from interference between nearby resonances, are completely suppressed. In other words, exceptional points are not necessary for the observation of the non-normal ME; the essential ingredient is non-normality itself.

4. Conclusions

In this work, we have identified and characterized a novel mechanism underlying the quantum Mpemba effect, arising from the non-normality of the Liouvillian superoperator in Markovian open quantum systems. Unlike previously known mechanisms – such as strong Mpemba effects tied to spectral suppression, or those relying on non-Markovian memory or initial correlations – our results demonstrate that transient interference between non-orthogonal decay modes can alone induce anomalous relaxation behavior. In particular, we have shown that certain initial states, despite being closer to thermal equilibrium, may exhibit slower early-time relaxation than states prepared farther away, purely due to non-normal dynamics. We illustrated this phenomenon using a waveguide QED model, a natural platform where directionality and radiative coupling introduce non-normality in a controlled manner. Our analysis reveals clear signatures of the non-normal quantum Mpemba effect, such as transient freezing and delayed relaxation, even in the absence of memory effects or asymptotic spectral degeneracies. These results broaden the theoretical understanding of quantum nonequilibrium phenomena and point to non-normality as a universal and designable resource for controlling thermalization dynamics.

Looking forward, this work opens several promising directions. On the theoretical side, it suggests a need to systematically classify non-normal contributions to quantum relaxation and to develop tools for diagnosing them in complex systems. Experimentally, waveguide QED platforms – both in photonic and superconducting circuit architectures – offer feasible testbeds to probe and exploit non-normal Mpemba behavior. Ultimately, harnessing such effects may provide new knobs for quantum state preparation, information processing, and thermodynamic control in engineered quantum environments.