Quantum Mpemba effect for operators in open systems

Abstract

The quantum Mpemba effect concerns with anomalous relaxation of quantum states that evolves either under unitary or non-unitary dynamics. In the context of open quantum systems, while most studies focus on quantum states evolving under completely positive trace-presing dynamics described by the Gorini–Kossakowski–Sudarshan–Lindblad (GKSL) master equation, we demonstrate that an analogous effect can arise at the level of operators. In particular, we show that operators that evolves under the adjoint Liouvillian – despite not being a trace-preserving map – can still exhibit a genuine Mpemba effect. We derive general conditions under which this phenomenon can occur and validate our predictions for three different open quantum setups. Our results broaden the scope of the Mpemba effect in quantum systems and provide a framework for controlling the relaxation of physically relevant observables.

Introduction.– Understanding the relaxation dynamics of quantum systems has long been a subject of fundamental research interest [1, 2, 3, 4, 5, 6]. The way a quantum system relaxes to its asymptotic steady-state, not only provides insight into the underlying physical mechanisms, but also plays a pivotal role in shaping the emerging field of quantum technology [7, 8, 9, 10, 11]. Despite significant progress in this direction, the study of non-equilibrium quantum dynamics continues to be an active area of research, particularly in the context of open quantum systems [12, 13] where the system of interest continuously interacts with the surrounding.

Recent theoretical studies in both isolated [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26] and open quantum systems [27, 28, 29, 30, 31, 32, 27, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50] have shown that a system initially farther from equilibrium can relax faster than one that is closer to it. This counterintuitive behavior is the quantum analogue of the Mpemba effect, originally brought to attention by Erasto B. Mpemba, who observed that hot water can freeze faster than cold [51]. Since then, the effect has been extensively studied in classical systems [52, 53, 54, 55, 56, 57, 58], with various mechanisms proposed to explain its origin.

In the quantum regime, the Mpemba effect has attracted growing interest, with studies demonstrating its emergence in both isolated [14, 15, 16, 17] and open systems [27, 28, 29]. In isolated systems, it emerges in the context of symmetry considerations: for two initial pure states that break a conserved $U(1)$ symmetry of an integrable Hamiltonian to different extents, the state with higher symmetry breaking can exhibit faster symmetry restoration at the subsystem level [14]. On the other hand, in open systems, the effect is controlled by the spectrum of the dynamical generator and the overlap with slowly decaying modes [27]. An initial state farther from the steady state, but orthogonal to the slowest decay mode (SDM), can relax faster than a closer state with finite overlap with the SDM [27].

Existing studies so-far have largely focused on the relaxation of quantum states [14, 27]. However, often experiments probe the relaxation of specific physical properties [59, 60, 61], represented by self-adjoint operators, rather than the full quantum state. This naturally raises the question of whether the Mpemba effect can manifest at the level of observables. Recent work [62] has taken initial steps in this direction by formulating the problem in operator space. However, a key challenge remains: the steady state of the operator dynamics generally depends on the choice of the initial operator [62]. This dependence complicates a direct extension of the Mpemba effect from states to observables and motivates the need for a framework in which the notion of steady-state relaxation is independent of the initial operator.

In this work, we develop a framework to characterize the Mpemba effect directly at the level of observables. By formulating the dynamics in operator space, we identify suitable operations under which physical observables can exhibit accelerated relaxation analogous to the traditional quantum state Mpemba effect. In particular, we address a key limitation of previous approaches, namely the dependence of the steady state on the initial operator [62], and provide a formulation in which both the operator and the transformed operator relax to the same steady state. Within this framework, we show that the relaxation speed of an observable is governed by its overlap with its slowest decay modes, and that appropriate transformations can bypass the overlap [27], leading to an exponential speed-up in the relaxation and faster convergence towards the steady state value. We further provide conditions on distance measure to observe Mpemba effect for operators. Our results establish that the Mpemba effect is not restricted to quantum states, but can manifest in experimentally relevant observables as well.

Recipe for observing genuine Mpemba effect for operators.– The time evolution of a quantum state interacting with the external environment is described by a quantum master equation derived under the standard Born-Markov and secular approximations [12, 13, 63]. Under these approximations, the reduced dynamics of the system is governed by the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) master equation [12], which provides the most general form of a Markovian, completely positive and trace-preserving quantum evolution. The GKSL equation is given by

$$\!\!\frac{d\rho}{dt}\!=\!\mathcal{L}\rho\equiv-i[H,\rho]\!+\!\sum_{k}\gamma_{k}\big(L_{k}\rho L_{k}^{\dagger}-\frac{1}{2}\{L_{k}^{\dagger}L_{k},\rho\}\big),$$ (1)

where the first term provides the unitary contribution with $H$ being the Hamiltonian of the system. The second term represents the dissipative contribution with $L_{k}$ being the jump operator for the $k$-th channel and $\gamma_{k}$ is the associated rate of dissipation. Since the Liouvillian $\mathcal{L}$ is a non-Hermitian superoperator, its spectral property is characterized by the left and right eigenoperators. These are defined through $\mathcal{L}r_{k}=\lambda_{k}r_{k},~~~\mathcal{L}^{\dagger}l_{k}=\lambda_{k}^{*}l_{k}$ with $\lambda_{k}$ being the eigenvalue of the Liouvillian. The left and right eigen-operators together form the complete basis and obey the bi-orthonormality condition, i.e., $\mathrm{Tr}(l_{i}r_{j})=\delta_{ij}$ [64, 13, 12].

In comparison to the state, the operator evolves in the Heisenberg picture, under the adjoint-Liouvillian $\mathcal{L}^{\dagger}$ [12]. The time evolution of an operator $\mathcal{O}$ is given by $\mathcal{O}(t)=e^{\mathcal{L}^{\dagger}t}\mathcal{O}(0)$, where $\mathcal{L}^{\dagger}$ reads as,

$$\mathcal{L}^{\dagger}\mathcal{O}=i[H,\mathcal{O}]+\sum_{k}\gamma_{k}\big(L_{k}^{\dagger}\mathcal{O}L_{k}-\frac{1}{2}\{L_{k}^{\dagger}L_{k},\mathcal{O}\}\big).$$ (2)

Following the spectral decomposition of $\mathcal{L}^{\dagger}$, the time evolved operator $\mathcal{O}(t)=e^{\mathcal{L}^{\dagger}t}$ can be expressed as,

$$\mathcal{O}(t)=\mathcal{O}_{\mathrm{ss}}+\sum_{i}\mathrm{Tr}\big[r_{i}\mathcal{O}(0)\big]e^{\lambda_{i}t}l_{i},$$ (3)

where $\mathcal{O}_{\mathrm{ss}}=\mathrm{Tr}\big[\rho_{\mathrm{ss}}\mathcal{O}(0)\big]$ and $\mathrm{Tr}\big[r_{i}\mathcal{O}(0)\big]$ is the overlap of the initial operator $\mathcal{O}(0)$ with the right eigenvector $r_{i}$. The time an operator takes to reach its steady state at the long-time scale is decided by its spectral gap $|\mathrm{Re}(\lambda_{1})|$ and the overlap $\mathrm{Tr}(r_{1}\mathcal{O})$. Bypassing this mode will cause an exponential speed-up in the relaxation process of the operator $\mathcal{O}(t)$, and the dynamics at the long time scale will be governed through the real part of the next SDM. In contrast to the states, bypassing the slowest relaxation mode $r_{1}$ through a unitary transformation on the operators $\mathcal{O}(0)$, i.e., $\mathrm{Tr}\big[r_{1}U^{\dagger}\mathcal{O}(0)U\big]=0$ also changes the asymptotic operator [12, 62]. In other words, the steady-state value $\mathcal{O}_{\mathrm{ss}}=\mathrm{Tr}\big[\rho_{\mathrm{ss}}\mathcal{O}(0)\big]$ changes to $\widetilde{\mathcal{O}}_{\mathrm{ss}}=\mathrm{Tr}\big[\rho_{\mathrm{ss}}\widetilde{O}(0)\big]$ where $\widetilde{O}(0)=U^{\dagger}\mathcal{O}(0)U$. Therefore, unlike for the quantum states, a unitarily transformed operator can relax faster but approaches to a different steady-state value and hence complicates the extension of the Mpemba effect from states to observables. Therefore, in order to achieve the same steady-state i.e., $\mathcal{O}_{\mathrm{ss}}=\widetilde{\mathcal{O}}_{\mathrm{ss}}$, one needs to go beyond the unitary transformation protocol. We here adopt a different approach that removes the contribution of the SDM from the operator while ensuring that the steady state remains unchanged. Specifically, we define the transformed operator as

$$\widetilde{\mathcal{O}}=\mathcal{O}-\mathrm{Tr}\big[r_{1}\mathcal{O}\big]l_{1},$$ (4)

where $r_{1}$ ($l_{1}$) denotes the right (left) eigenvector corresponding to the SDM of the operator. It is important to note that the SDM for an operator need not always be the same as the SDM of the Liouvillian. An important consequence of this construction is that the locality properties of the transformed operator $\widetilde{\mathcal{O}}$ depend on the structure of the removed SDM. In particular, if the corresponding left eigenoperator $l_{1}$ is local, then a local initial operator $\mathcal{O}$ remains local after the transformation. However, when the SDM is intrinsically nonlocal, the subtraction term $\mathrm{Tr}[r_{1}\mathcal{O}]\,l_{1}$ introduces nonlocal contribuition into $\widetilde{\mathcal{O}}$. Consequently, even if the original operator $\mathcal{O}$ is strictly local, the transformed operator can become highly nonlocal [for details see [65]]. Interestingly, the construction in Eq. (4) preserves the steady state while enabling faster relaxation. In particular,

$$\widetilde{\mathcal{O}}_{\mathrm{ss}}=\mathrm{Tr}\big[\widetilde{\mathcal{O}}\rho_{\mathrm{ss}}\big]\mathbb{I}=\mathrm{Tr}\big[{\mathcal{O}}\rho_{\mathrm{ss}}\big]\mathbb{I}-\mathrm{Tr}\big[r_{1}\mathcal{O}\big]\mathrm{Tr}\big[l_{1}\rho_{\mathrm{ss}}\big]\mathbb{I}.$$ (5)

Due to the bi-orthonormality condition, $\mathrm{Tr}[l_{1}\rho_{\mathrm{ss}}]=0$ as $\rho_{\mathrm{ss}}=r_{0}$ and as a consequence $\widetilde{\mathcal{O}}_{\mathrm{ss}}={\mathcal{O}}_{\mathrm{ss}}$. Therefore, the steady-state value of the operator before and after transformation remains unchanged. More generally, this procedure can be extended to eliminate multiple consecutive SDMs. In that case we can write the transformed operator as

$$\widetilde{\mathcal{O}}=\mathcal{O}-\sum_{i}\mathrm{Tr}\big[r_{i}\mathcal{O}\big]\,l_{i},$$ (6)

which leads to a further acceleration of the relaxation dynamics. Note that this procedure only guarantees faster convergence of the time-evolved operator $\mathcal{O}(t)$ to its asymptotic value. However, in order to probe the emergence of the Mpemba effect, it is essential to compare the initial distances: specifically, how far the transformed operator $\widetilde{\mathcal{O}}(0)$ lies from its steady-state $\mathcal{O}_{\rm ss}$ relative to the original operator $\mathcal{O}(0)$. This is what we investigate next.

Distance measure and condition for Mpemba effect for operators.– The dynamical map $e^{\mathcal{L}t}$ governing the evolution of a quantum state is completely positive and trace preserving (CPTP) [12]. As a consequence, contractive distance measures such as the trace distance and quantum relative entropy decay monotonically with time [12, 66, 67]. These measures are therefore commonly used to quantify the distance between the time-evolved state $\rho_{t}$ and the steady state $\rho_{\mathrm{ss}}$, and serve as standard tools to identify the standard Mpemba effect [67]. In contrast, operators evolve under the adjoint map $e^{\mathcal{L}^{\dagger}t}$, which is, in general, not trace preserving [12] (except in the special case of unital dynamics) As a result, standard distance measures such as the trace distance [67, 66] do not necessarily exhibit monotonic decay under adjoint evolution, and hence are not suitable for analyzing the Mpemba effect at the level of operators.

Recent study  [62] has shown that the dressed distance, which can be interpreted as a trace distance under the Petz recovery map [68, 69], remains contractive even under a non-trace-preserving dynamics, generated by $e^{\mathcal{L}^{\dagger}t}$. This makes it a standard distance measure for studying the Mpemba effect in the operator framework. The dressed distance is defined as [62]

$$\mathcal{D}_{\mathrm{dd}}(\mathcal{O},\mathcal{O}_{\mathrm{ss}})=\mathrm{Tr}\left(\sqrt{X^{\dagger}X}\right),$$ (7)

where $X=\sqrt{\rho_{\mathrm{ss}}}\,(\mathcal{O}-\mathcal{O}_{\mathrm{ss}})\,\sqrt{\rho_{\mathrm{ss}}}.$ For unital dynamics ($\rho_{\mathrm{ss}}=\mathbb{I}/d$ with $d$ being the dimension of the Hilbert space), the dressed distance reduces to the standard trace distance up to a constant prefactor.

Following the mode-bypassing procedure introduced in Eq. (4), we construct a transformed operator $\widetilde{\mathcal{O}}$ by removing the contribution of its SDM. By construction, $\widetilde{\mathcal{O}}$ relaxes faster towards the steady state. The key question is whether this transformation also increases the initial distance from the steady-state, i.e., $\mathcal{D}_{\mathrm{dd}}(\widetilde{\mathcal{O}},\mathcal{O}_{\mathrm{ss}})>\mathcal{D}_{\mathrm{dd}}(\mathcal{O},\mathcal{O}_{\mathrm{ss}}).$ If this condition is satisfied at $t=0$, the operator can exhibit the Mpemba effect during the relaxation process. To analyze this further, we define

$$A=\sqrt{\rho_{\mathrm{ss}}}\,(\mathcal{O}-\mathcal{O}_{\mathrm{ss}})\sqrt{\rho_{\mathrm{ss}}},\quad B=\sqrt{\rho_{\mathrm{ss}}}\Big[\mathrm{Tr}\big(r_{1}\mathcal{O}\big)l_{1}\Big]\sqrt{\rho_{\mathrm{ss}}},$$ (8)

so that the change in the dressed distance can be expressed as $\Delta=\|A-B\|-\|A\|,$ where $\|X\|=\mathrm{Tr}(\sqrt{X^{\dagger}X})$ is the trace norm. Invoking the polar decomposition form for operators [70, 71] $X=U_{X}|X|$ where $|X|=\sqrt{X^{\dagger}X}$, one can show that [for details see [65]]

$$\|X\|=\mathrm{Tr}(U_{X}^{\dagger}X)\geq\mathrm{Re}\big[\mathrm{Tr}(U^{\dagger}X)\big],$$ (9)

where $U$ is an arbitrary unitary operator and the equality only holds for $U=U_{X}$. Applying the polar decomposition form and the inequality in Eq. (9) to $A$ and $B$ in Eq. (8), we obtain the lower bound on $\Delta$ as

$$\|A-B\| =\mathrm{Re}\big[\mathrm{Tr}[U^{\dagger}_{A-B}(A-B)]\big]$$

$$\geq\mathrm{Re}\big[\mathrm{Tr}[U^{\dagger}_{A}(A-B)]\big]=\|A\|-\mathrm{Re}\big[\mathrm{Tr}[U^{\dagger}_{A}(B)]\big]$$

$$\Delta \geq-\mathrm{Re}\big[\mathrm{Tr}[U^{\dagger}_{A}(B)]\big],$$ (10)

where recall that $\Delta=\|A-B\|-\|A\|.$ Similarly, one can obtain an upper bound on $\Delta$ as

$$\|A\| \geq\mathrm{Re}\big(\mathrm{Tr}[U^{\dagger}_{A-B}(A)]\big)$$

$$=\mathrm{Re}\big(\mathrm{Tr}[U^{\dagger}_{A-B}(A-B+B)]\big)$$

$$\geq\|A-B\|+\mathrm{Re}\big(\mathrm{Tr}[U^{\dagger}_{A-B}B]\big)$$

$$\Delta \leq-\mathrm{Re}\big(\mathrm{Tr}[U^{\dagger}_{A-B}B]\big).$$ (11)

From the above two conditions, as given in Eq. (Quantum Mpemba effect for operators in open systems) and Eq. (Quantum Mpemba effect for operators in open systems), we receive the following inequality for $\Delta$

$$-\mathrm{Re}\Big[\mathrm{Tr}(U_{A}^{\dagger}B)\Big]\leq\Delta\leq-\mathrm{Re}\Big[\mathrm{Tr}(U_{A-B}^{\dagger}B)\Big].$$ (12)

Therefore, a sufficient condition for the increase of the dressed distance i.e., $\Delta\geq 0$ is $\mathrm{Re}\,\big[\mathrm{Tr}(U_{A}^{\dagger}B)\big]\leq 0,$ while a necessary condition is given by $\mathrm{Re}\,\big[\mathrm{Tr}(U_{A-B}^{\dagger}B)\big]\leq 0.$ Interestingly, the sufficient condition admits a simple geometric interpretation. The unitary operator $U_{A}$, obtained from the polar decomposition $A=U_{A}|A|$, captures the orientation of the operator $A$ in operator space. This is analogous to the phase factor in the polar representation of a complex number, $z=re^{i\theta}$, where $e^{i\theta}$ determines the orientation of $z$ in the complex plane. The quantity $\mathrm{Tr}(U_{A}^{\dagger}B)$ can be interpreted as the Hilbert-Schmidt inner product [36] $\langle U_{A},B\rangle$ that measures the alignment between $A$ and the contribution $B$ arising from the removed mode. When this quantity is negative, it indicates that $B$ is oppositely oriented to $A$ in operator space. Consequently, subtracting such a contribution (i.e., bypassing these modes) increases the norm $\|A-B\|$ relative to $\|A\|$, leading to an increase in the dressed distance. This provides an intuitive picture: removing components that are anti-aligned with the operator enhances its initial distance from the steady state, thereby enabling the Mpemba effect for operators. In what follows, we illustrate the genuine Mpemba effect for operators for three paradigmatic open quantum setups.

Example 1: Single-qubit evolving under a unital map.– As a first example, we consider a single qubit with Hamiltonian $H=\omega_{0}\sigma_{z}/2$, subjected to dissipation through the jump operator $\sigma_{x}$ with dissipation strength $\gamma$. In this case, the steady state is maximally mixed, $\rho_{\mathrm{ss}}=\mathbb{I}_{2}/2,$ and the corresponding dynamical map is unital. The Liouvillian eigenvalues are given by $0,-\gamma\pm\sqrt{\gamma^{2}-\omega_{0}^{2}}$, and $-2\gamma$. Depending on the relative magnitudes of $\gamma$ and $\omega_{0}$, the spectrum can consist of purely real eigenvalues, a complex-conjugate pair, or an exceptional point of second order. Here, we focus on the regime $\gamma>\omega_{0}$, where all eigenvalues are real.

Figure 1: Presence of the genuine operator Mpemba effect in a single qubit setup evolving under a unital map. (a) Time evolution of the dressed distance for the original operator $\mathcal{O}=b_{0}\mathbb{I}+\vec{b}\cdot\vec{\sigma}$ and the transformed operator $\widetilde{\mathcal{O}}$, obtained by bypassing the SDM. Although $\widetilde{\mathcal{O}}$ is initially farther from the asymptotic value, it relaxes significantly faster due to the removal of the SDM. (b) Time evolution of the corresponding expectation values evaluated in the initial state $\ket{\psi_{0}}=\ket{+}$, demonstrating the accelerated relaxation of the transformed operator $\widetilde{\mathcal{O}}$. The parameters used are $\omega_{0}=0.5$, $\gamma=1$, $b_{0}=1$, and the vector $\vec{b}$ is chosen randomly with values for the components given as $(-0.01954,\,0.52816,\,0.09795)$.

In this regime, the SDM of the Liouvillian $\mathcal{L}$ corresponds to the eigenvalue $\lambda_{1}=-\gamma+\sqrt{\gamma^{2}-\omega_{0}^{2}},$ which allows for selectively removing it through mode-bypassing protocol, as introduced in Eq. (4). In Ref. [62], it was shown that the operator Mpemba effect is absent for unital dynamics when the acceleration protocol is implemented through a unitary transformation. In contrast, we show here that a genuine operator Mpemba effect can emerge under the same unital dynamics through the mode-bypassing procedure.

To illustrate this point, we consider a generic single-qubit operator $\mathcal{O}=b_{0}\mathbb{I}+\vec{b}\cdot\vec{\sigma},$ and construct the transformed operator $\widetilde{\mathcal{O}}$ by eliminating its overlap with the SDM. As a consequence, the transformed operator relaxes on a faster timescale governed by the next slowest eigenmode. This behavior is shown in Fig. (1). Fig 1(a) shows the evolution of the dressed distance, where the transformed operator $\widetilde{\mathcal{O}}$ initially lies farther from the steady state but converges more rapidly than the original operator $\mathcal{O}$. Fig. 1(b) shows the corresponding expectation values of the operator, before and after transformation, and clearly highlighting the accelerated relaxation induced by the mode-bypassing protocol.

Interestingly, in contrast to the unital map with $\sigma_{x}$ dissipator, if we consider the single qubit setup evolving under a Davies map [36, 72] where the $\sigma_{+}$ and $\sigma_{-}$ are the Lindblad dissipators with rates satisfying the detailed balance condition, the operator Mpemba effect is completely absent (for details see Ref. [65]). A similar absence of the Mpemba effect was recently reported in Ref. [62], where accelerated relaxation was induced through a unitary transformation.

Example 2: Non-equilibrium XX-spin chain.–

Figure 2: Results for the Non-equilibrium $XX$-spin chain: Demonstration of the operator Mpemba effect for operators with different degree of locality for the boundary-driven spin chain with $N=4$. (a) Local operator: population imbalance $\mathcal{O}=n_{2}-n_{3}$. The transformed operator $\widetilde{\mathcal{O}}$ exhibits faster convergence to the steady state despite initially larger dressed distance. (b) Dynamics for the corresponding expectation values starting with the initial state $\ket{\psi_{0}}=\ket{1}^{\otimes 4}$, ($\ket{1}$ is the spin down state) showing accelerated relaxation and convergence to the same steady state. (c) Two-site operator: current $I_{1,2}$. The transformed operator $\widetilde{I}_{1,2}$ relaxes exponentially faster than $I_{1,2}$, despite a larger initial dressed distance. (d) Dynamics for the corresponding expectation values starting with the initial state $\ket{\psi_{0}}=\ket{0}^{\otimes 4}$, ($\ket{0}$ is the spin up state), confirming faster relaxation to the same steady state. (e) Non-local operator: $\mathcal{O}=\mathbb{I}_{1}\otimes\sigma_{z}^{(2)}\otimes\sigma_{z}^{(3)}\otimes\sigma_{z}^{(4)}-\sigma_{z}^{(1)}\otimes\sigma_{z}^{(2)}\otimes\sigma_{z}^{(3)}\otimes\mathbb{I}_{4}.$ Although $\widetilde{\mathcal{O}}$ is far from the steady state, it relaxes faster due to vanishing overlap with the five dominant SDMs. The inset shows the zoomed version of the initial distance. (f) Corresponding expectation values for the initial state $\ket{\psi_{0}}=\ket{1}^{\otimes 4}$, showing accelerated convergence. For simulation, we use the parameters $g=0.05$, $\gamma_{1}=0.2$, $\gamma_{N}=0.5$, $\epsilon_{1,2,3,4}=(1,0.8,0.6,0.5)$, $T_{1}=1.5$, $T_{N}=1$, and $\mu_{1}=\mu_{N}=0$.

As a second example, we consider a boundary-driven $XX$ spin chain described by the Hamiltonian $H=\sum_{i=1}^{N}\epsilon_{i}\,\sigma_{+}^{(i)}\sigma_{-}^{(i)}+g\sum_{i=1}^{N-1}\big(\sigma_{+}^{(i)}\sigma_{-}^{(i+1)}+\sigma_{-}^{(i)}\sigma_{+}^{(i+1)}\big)$, where $\epsilon_{i}$ is the onsite energy, $g$ is the nearest-neighbor hopping strength and $N$ is the number of qubits. The dynamics of the setup is governed by Eq. (1), with jump operators corresponding to two boundary sites are given by $L_{i}^{\pm}=\sigma_{\pm}^{(i)}$, $i\in\{1,N\}$, and associated jump rates are $\gamma_{i}f_{i}(\epsilon_{i})$ and $\gamma_{i}[1-f_{i}(\epsilon_{i})]$, respectively. Here $f_{i}(\epsilon_{i})$ is the Fermi distribution function evaluated at the onsite energy $\epsilon_{i}$ with temperature $T_{i}$ and chemical potential $\mu_{i}$, which are fixed at different values at the two edges. As a result, the system is driven out-of-equilibrium and the setup possesses a unique current carrying steady-state. Within this setup, we investigate the emergence of the genuine Mpemba effect for the operators using the mode-bypassing procedure.

We consider three representative operators with an increasing degree of non-locality. Here we present results for $N=4$. First, as a local (single-site) observable, we study the population imbalance between the consecutive sites, defined as $n_{i}-n_{(i+1)}$ with $n_{i}=(\mathbb{I}+\sigma_{z}^{(i)})/2$. Next, as a two-site observable, we consider bond current which is defined as $I_{i,i+1}=-ig\big[\mathbb{I}_{1}\otimes..\otimes\big(\sigma_{-}^{(i)}\otimes\sigma_{+}^{(i+1)}-\sigma_{+}^{(i)}\otimes\sigma_{-}^{(i+1)}\big)\otimes\mathbb{I}_{i+2}...\otimes\mathbb{I}_{N}\big].$ Finally, as a more non-local (three-site) observable, we analyze a string operator of the form $\mathcal{O}=\mathbb{I}_{1}\otimes\sigma_{z}^{(2)}\otimes\sigma_{z}^{(3)}\otimes\sigma_{z}^{(4)}-\sigma_{z}^{(1)}\otimes\sigma_{z}^{(2)}\otimes\sigma_{z}^{(3)}\otimes\mathbb{I}_{4}.$

In Fig. (2), we demonstrate the emergence of the operator Mpemba effect for these different classes of observables. For the local operator corresponding to imbalance, the transformed operator exhibits faster relaxation compared to the original one, despite being initially farther from the steady state in the dressed distance measure, as shown in Fig (2) (a)-(b). A similar behavior is observed for the two-body current operator $I_{1,2}$, where bypassing the SDM leads to a clear exponential speedup in the relaxation dynamics, as clearly demonstrated in Fig (2) (c)-(d). The effect persists even for the non-local three-site string operator. Here, we bypass the five dominant slow decay modes of the initial observable using the multimode bypassing procedure. Although the transformed observable is initially farther from the steady state, it nevertheless reaches the steady state more rapidly than the original observable, thereby exhibiting the operator Mpemba effect, as clearly shown in Fig (2) (e)-(f). In all these cases, we verify the validity of the bound in Eq. (12) and present results in [65].

Example 3: One-dimensional bosonic lattice with bulk dephasing.–

Figure 3: Results for one-dimensional bosonic lattice subjected to dephasing: (a) The transformed operator $\widetilde{n}_{2}(t)$ initially lies farther from the steady state value $n_{2}^{(\mathrm{ss})}$ in terms of dressed distance (for this case it reduces to trace distance, as the map is unital), yet it relaxes more rapidly than the local number density $n_{2}$ at the second site. In contrast, $n_{2}$, although closer to equilibrium at initial times, exhibits slower relaxation. This counterintuitive behavior highlights the presence of the genuine Mpemba effect for operators. (b) The expectation value evaluated with respect to the initial state $\rho_{0}=\frac{1}{L-1}\sum_{i=1}^{L-1}\ket{i}\bra{i}$ further confirms the accelerated relaxation of the transformed operator. The simulations are performed for parameters $L=11$, $J=1$, and $\gamma=5$.

As a final example, we consider a one-dimensional (1D) bosonic lattice of $L$ coupled sites with nearest-neighbor tunneling. The Hamiltonian for the lattice is given as $H=\sum_{n=1}^{L-1}J\left(a_{n}^{\dagger}a_{n+1}+\mathrm{h.c.}\right)$. Here, $J$ denotes the tunneling amplitude between adjacent sites, and $a_{n}^{\dagger}$ ($a_{n}$) are the bosonic creation (annihilation) operators at site $n$. The lattice is further subjected to uniform onsite dephasing [73, 74] with dephasing strength $\gamma_{i}=\gamma$, and the jump operator is given as $L_{n}=a_{n}^{\dagger}a_{n}$, $n=1,2,\cdots L$. For this setup, we investigate the presence of genuine operator Mpemba effect by focusing on the single particle sector under strong dephasing limit $\gamma\gg J$ [73] and intermediate dephasing limit $\gamma\sim J$. Here, we discuss the results in the strong-dephasing limit, and the results for moderate dephasing, $\gamma\sim J$, are presented in Ref. [65].

In the strong dephasing limit, the coherences in the site basis are strongly suppressed, leading to the evolution confined entirely in the population sector $p_{n}=\rho_{n,n}(t)$. The population dynamics can be described by a Pauli master equation [12, 73, 75] $\frac{d\mathbb{P}}{dt}=\mathcal{W}\,\mathbb{P}$ where $\mathbb{P}=(p_{1},p_{2},\cdots,p_{L})^{T}$. Since the rate matrix $\mathcal{W}$ is Hermitian, its left and right eigenmatrices are identical. These eigenmodes govern the relaxation dynamics both for the operators and states and play the central role in the Mpemba effect. For a lattice of size $L$, the SDM of $\mathcal{W}$ takes the form $r_{1}=\sum_{n=1}^{L}\Phi_{1}(n)\ket{n}\bra{n},$ [73] where $\Phi_{1}(n)=\sqrt{\frac{2}{L}}\cos\!\left[\frac{\pi(2n-1)}{2L}\right]$, and $n$ labels the lattice sites. To reveal the operator Mpemba effect, we analyze the dynamics of the local observable $n_{i}=a_{i}^{\dagger}a_{i}$ and its transformed counterpart $\widetilde{n}_{i}$, constructed via the mode-bypassing protocol, as given in Eq. (4).

As shown in Fig. (3)(a), the bare observable $n_{2}$ exhibits slow relaxation towards its steady-state value $n_{2}^{(\mathrm{ss})}$, governed by its finite overlap with the SDM $r_{1}$. In contrast, the transformed operator $\widetilde{n}_{2}$ is explicitly engineered to eliminate this overlap. As a result, it bypasses the SDM and relaxes significantly faster, despite being initially farther from the steady state. This effect is demonstrated in Fig. (3)(b), where $\langle\widetilde{n}_{2}(t)\rangle$ displays a clear exponential speed-up compared to $\langle n_{2}(t)\rangle$. This provides direct evidence of a genuine operator Mpemba effect driven by mode bypassing. Crucially, the transformation preserves the steady state, ensuring that both observables relax to the same asymptotic value while following markedly different dynamics. We further show in [65] that bypassing multiple slow decay modes leads to an even stronger enhancement of the relaxation dynamics, resulting in a further reduction of the relaxation timescale.

Summary.– We investigate the emergence of a genuine Mpemba effect for operators in Markovian open quantum systems. Unlike quantum states, operators evolve under the adjoint Liouvillian and therefore do not follow a CPTP map. Exploiting this distinction, we introduce a systematic procedure to accelerate operator relaxation by selectively removing the slowest decaying mode, while preserving convergence to the same steady-state value. Furthermore, we derive a quantitative bound on the change in the dressed distance induced by this mode-bypassing transformation. Our results are particularly relevant in scenarios where the rapid equilibration of specific observables is more important than the relaxation of the full quantum state. Our work thus establishes that the quantum Mpemba effect is not restricted to the relaxation of quantum states, but can also occur directly at the level of physical observables. This broadens the conceptual scope of anomalous relaxation phenomena in quantum systems and opens promising future directions that includes extending the present analysis for time-dependent systems and for systems evolving under non-Markovian dynamics.

Acknowledgements.– PB acknowledges the University Grants Commission (UGC), Government of India, for the research fellowship (Ref No.- 231620073714). BKA acknowledges CRG Grant No. CRG/2023/003377 from Science and Engineering Research Board (SERB), Government of India.

Supplemental Material: Quantum Mpemba effect for operators in open systems

I Nature of transformed operators under Mode Bypassing

In this section, we discuss how the mode-bypassing protocol can transform an initially local operator into a nonlocal one. The relaxation dynamics of an operator is governed by its decomposition in the eigenmodes of the underlying dynamical generator, which in this case is the adjoint Liouvillian $\mathcal{L}^{\dagger}$. Any operator $\mathcal{O}$ with a finite overlap with its slowest decay mode (SDM) decides its relaxation timescale in the long-time limit. Therefore, controlling this overlap provides a direct route to modify the relaxation behavior. In particular, the mode-bypassing procedure introduced in Eq. (4) suppresses the contribution of the SDM, thereby accelerating the relaxation of the transformed operator. However, this acceleration comes at a cost. If one starts from a local operator, the application of the mode-bypassing transformation generally produces a nonlocal operator. This originates from the fact that the right and left eigenvectors of the Liouvillian $\mathcal{L}$ are, in general, nonlocal when expressed in the Pauli operator basis. Consequently, eliminating the overlap with the SDM inevitably introduces nonlocal components into the transformed operator.

To illustrate this explicitly, we consider Example 2 from the main text, namely the boundary-driven XX spin chain. Considering the initial operator to be the local number operator at the second site, $\mathcal{O}=(\mathbb{I}+\sigma_{z}^{(2)})/2$, we find that applying the mode-bypassing procedure to eliminate its overlap with the SDM and produces a transformed operator that is manifestly nonlocal in Pauli basis. As shown in Fig. S1, although the original operator has support only on a single site, the transformed operator develops finite contributions from both local and increasingly nonlocal operator sectors. In particular, the weight of multi-site operator components, which we define as

$$|c|=\mathrm{Tr}(\widetilde{\mathcal{O}}\,\sigma_{i}^{\otimes 4})/\mathrm{Tr}[r_{1}\mathcal{O}]$$ (S1)

clearly demonstrates that the elimination of the SDM generates extended correlations across the system. Here $\sigma_{i}\in[\mathbb{I}_{2},\sigma_{x}\sigma_{y},\sigma_{z}]$.

It is important to emphasize that the behavior of the transformed operator depends on the structure of the SDM itself. For example, as shown in Example 3 in the main text, under the strong dephasing limit ($\gamma\gg J)$ the SDM is local in the number basis. Consequently, bypassing the SDM in that case does not introduce non-locality into the transformed operator, as a result accelerated or transformed operator remains local.

Figure S1: Plot to show the nonlocal behaviour of the transformed operator generated by the mode-bypassing procedure in the non-equilibrium $XX$-spin chain model [Example 2 of the main text], starting from a local operator $\mathcal{O}=(\mathbb{I}+\sigma_{z}^{(2)})/2$. The figure shows the distribution of the magnitudes of the Pauli decomposition coefficients, $|c|$, as defined in Eq. (S1) of the transformed operator, grouped according to their Pauli weight (from 1-site to 4-site sectors). Each panel corresponds to a fixed locality sector and displays the coefficients associated with Pauli strings belonging to that sector. Although the initial operator is strictly local, the transformed operator develops finite contributions from increasingly nonlocal sectors after removing the slowest decay mode. The coefficients are normalized by the overlap of the initial operator with the right eigenvector $r_{1}$ of $\mathcal{L}$, corresponding to $\lambda_{1}$ and only terms satisfying $|c|>10^{-4}$ are retained. The parameters used are $g=0.05$, $\gamma_{1}=0.2$, $\gamma_{N}=0.5$, $\epsilon_{1,2,3,4}=(1,0.8,0.6,0.5)$, $T_{1}=1.5$, $T_{N}=1$, and $\mu_{1}=\mu_{N}=0$.

I.1 Experimental strategies to measure the transformed operator

In this subsection, we discuss possible experimental strategies to measure the expectation value of the transformed operators which are generally non-local. In particular, we introduce two different approaches that allow one to access this quantity in practice.

Method 1.– Recall that, the transformed operator is defined as $\widetilde{\mathcal{O}}=\mathcal{O}-\mathrm{Tr}(r_{1}\mathcal{O})\,l_{1},$ where $r_{1}$ and $l_{1}$ denote the right and left eigenoperators corresponding to the SDM of the operator $\mathcal{O}$. This transformation can be interpreted as projecting $\mathcal{O}$ out of the slowest decay mode (SDM) using the projector $\Pi=\mathbb{I}-\ket{l_{1}}\bra{r_{1}}.$, where $\mathbb{I}$ is the Identity operator. As a result, the transformed operator $\widetilde{\mathcal{O}}$ relaxes with a faster timescale $\tau_{2}=1/|\mathrm{Re}(\lambda_{2})|$, where $\lambda_{2}$ is the next slowest decay mode of the operator $\mathcal{O}$. The expectation value of the transformed operator evolves as

$$\mathrm{Tr}\big(\widetilde{\mathcal{O}}(t)\rho_{0}\big)=\mathrm{Tr}(\rho_{t}\mathcal{O})-\mathrm{Tr}(r_{1}\mathcal{O})\,\mathrm{Tr}(l_{1}\rho_{0})\,e^{-\lambda_{1}t}.$$ (S2)

This expression provides a practical route for experimental implementation. Instead of directly measuring the transformed operator $\widetilde{\mathcal{O}}$, one can measure the expectation value of the original operator $\mathcal{O}$ and subtract the known contribution from the SDM [46].

Method 2.– We next discuss an ancilla-based technique to measure the expectation value of the transformed operator $\widetilde{\mathcal{O}}(t)$. We introduce an ancilla qubit prepared in the superposition state $\ket{\psi^{\mathrm{anc}}}=\frac{\ket{0}+\ket{1}}{\sqrt{2}},$ with corresponding density matrix, given as $\rho_{A}=\ket{\psi^{\mathrm{anc}}}\bra{\psi^{\mathrm{anc}}}$. Just before the measurement $(t^{-})$, the total system-ancilla state is given as a product state i.e., $\rho_{\mathrm{tot}}(t)=\rho_{A}\otimes\rho_{s}(t).$ One can then introduce an interaction between the system and the ancilla governed by the Hamiltonian $H_{\mathrm{int}}=g\,\sigma_{z}^{\mathrm{anc}}\otimes\widetilde{\mathcal{O}},$ applied for a duration $\tau$. Here $g$ is the interaction strength. The corresponding unitary evolution is given by

$$U=e^{-ig\tau\sigma_{z}\otimes\widetilde{\mathcal{O}}}=\ket{0}\bra{0}\otimes e^{-ig\tau\widetilde{\mathcal{O}}}+\ket{1}\bra{1}\otimes e^{ig\tau\widetilde{\mathcal{O}}}.$$ (S3)

After the interaction, the total state becomes $\rho_{\mathrm{tot}}(t^{+})=U\big(\rho_{A}\otimes\rho_{s}(t)\big)U^{\dagger}.$ The reduced state of the ancilla is obtained by tracing out the system i.e., $\rho_{A}(t^{+})=\mathrm{Tr}_{s}\big[\rho_{\mathrm{tot}}(t^{+})\big].$ The off-diagonal element of the ancilla is given by

$$\bra{0}\rho_{A}(t^{+})\ket{1}=\frac{1}{2}\mathrm{Tr}_{s}\big[\rho_{s}(t)\,e^{-2ig\tau\widetilde{\mathcal{O}}}\big].$$ (S4)

In the weak-coupling and short interaction time limit, i.e., when $g\tau\ll 1$, one can write $e^{-2ig\tau\widetilde{\mathcal{O}}}=\mathbb{I}-2ig\tau\widetilde{\mathcal{O}}+\mathcal{O}((g\tau)^{2}).$ Using this expansion, one finds that the expectation value of the ancilla observable $\sigma_{y}$ is directly proportional to the expectation value of the transformed operator

$$\langle\sigma_{y}(t+\tau)\rangle=2g\tau\,\langle\widetilde{\mathcal{O}}(t)\rangle.$$ (S5)

This relation provides an operational protocol to measure the expectation value of the transformed operator $\widetilde{\mathcal{O}}$ indirectly via the ancilla [76, 77], avoiding the need for direct measurement of a generic non-local operator.

II Proof for the bound in Eq. (9)

In this section, we aim to show that for any arbitrary unitary operator $U$,

$$\|X\|=\mathrm{Tr}(U_{X}^{\dagger}X)\geq\mathrm{Re}\,\big[\mathrm{Tr}(U^{\dagger}X)\big],$$ (S6)

with equality holds if and only if $U=U_{X}$. This condition was used in the main text to derive the lower and upper bounds for the dressed distance, as given in Eq. (12). Note that $X=U_{X}|X|$ is the polar decomposition for $X$, and $|X|=\sqrt{X^{\dagger}X}\geq 0$. Following the polar decomposition, we can write

$$\|X\|=\mathrm{Tr}|X|=\mathrm{Tr}(U_{X}^{\dagger}X).$$ (S7)

For an arbitrary unitary operator $U$, we define $V=U^{\dagger}U_{X}$, which is also a operator. As a result, we can write

$$\mathrm{Re}\,\big[\mathrm{Tr}(U^{\dagger}X)\big]=\mathrm{Re}\,\big[\mathrm{Tr}(U^{\dagger}U_{X}|X|)\big]=\mathrm{Re}\,\big[\mathrm{Tr}(V|X|)\big].$$ (S8)

Let us now consider the spectral decomposition of $|X|$ and write $|X|=\sum_{i}s_{i}\ket{i}\bra{i}$, with $s_{i}\geq 0$. Then, we can write

$$\mathrm{Re}\,\big[\mathrm{Tr}(V|X|)\big]=\sum_{i}s_{i}\,\mathrm{Re}\,\big[\bra{i}V\ket{i}\big].$$ (S9)

Employing the Cauchy–Schwarz inequality, we write

$$|\langle{i}|j\rangle|^{2}\leq|\langle i|i\rangle|^{2}|\langle j|j\rangle|^{2}=1,$$ (S10)

where $\ket{j}=V\ket{i}$ and as $V$ is a unitary operator, it does not change the norm of the vector $\ket{i}$. Now, for any complex number $z$, $\mathrm{Re}\big[z\big]\leq|z|$, we can write

$$\mathrm{Re}\,\big[\bra{i}V\ket{i}\big]\leq 1.$$ (S11)

Therefore,

$$\|X\|=\mathrm{Tr}|X|=\sum_{i}s_{i}\geq\mathrm{Re}\,\big[\mathrm{Tr}(U^{\dagger}X)\big],$$ (S12)

with equality holds if and only if $\mathrm{Re}\,\big[\bra{i}V\ket{i}\big]=1$ for all $i$, which implies $\bra{i}V\ket{i}=1$ for all $i$. This is possible only when $V$ is an identity operator, and hence $U=U_{X}$.

III Enhanced Relaxation for Operators via Multiple Mode-Bypassing

Figure S2: Enhanced relaxation via successive mode-bypassing. (a) Dressed distance dynamics displaying the relaxation dynamics for three cases: the original operator, taken as the local site density $(n_{2})$ (blue dashed), the transformed operator obtained by bypassing the slowest decay mode (SDM1) (green solid), exhibiting the Mpemba effect (see inset), and the further transformed operator obtained by bypassing the two slowest modes (SDM1 and SDM2) (red dotted), demonstrating enhanced relaxation. (b) Corresponding dynamics for operator expectation values evaluated with respect to the initial state $\rho_{0}=\frac{1}{L-1}\sum_{i=1}^{L-1}|i\rangle\langle i|$, showing progressively faster convergence of the local density $n_{2}$ to the same steady state. For the numerics, we consider the parameters $L=11$, $J=1$, and the dephasing strength $\gamma=5$.

In the main text, we established that the accelerated convergence of the transformed operator $\widetilde{\mathcal{O}}$ to the steady state arises by eliminating the contribution of the SDM of the original operator $\mathcal{O}$. This procedure can be systematically extended to further accelerate the dynamics by successively bypassing multiple slow decay modes. To illustrate this explicitly, we revisit Example 3 of the main text, which is a one-dimensional bosonic setup subjected to dephasing. We showed in the main text that, starting from the population density of a local site $(n_{i})$ as the initial operator, we bypass the contribution of its SDM and observe the emergence of the Mpemba effect for operators. One can however proceed further by bypassing other slow decaying modes, which can lead to a further acceleration of the relaxation towards the asymptotic value.

In Fig. S2(a) and (b), we demonstrate the enhanced accelerated relaxation where we bypass not only the slowest decay mode, which we call SDM1, but also the next slow decay mode, SDM2. In Fig.  S2(a), we plot the dynamics of the dressed distance for the local density operator. We observe a significant enhancement in relaxation upon bypassing the two slowest decay modes (SDM1 and SDM2), as shown by the red dotted curve, compared to the relaxation of the original operator (blue dashed line), and the transformed operator with only the slowest decay mode (SDM1) bypassed (green solid line). The associated faster convergence in the expectation value of the operator to the same steady state is observed in Fig. S2(b).

IV Absence of operator Mpemba effect in single qubit evolving under Davies map

In this section, we show that for an arbitrary operator in the single-qubit Hilbert space, evolving under the Davies map, there can not be Mpemba effect for operators. Any general operator in the single-qubit Hilbert space can be written as $\mathcal{O}=b_{0}\,\mathbb{I}+\vec{b}\cdot\vec{\sigma}$, where $\vec{b}\in\mathbb{R}^{3}$ is a real three-dimensional vector and $\vec{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ is the vector of Pauli matrices. The single-qubit system is described by the Hamiltonian: $H=\omega_{0}\,\sigma_{z}/2$, where $\sigma_{z}$ is the Pauli-Z operator and $\omega_{0}$ is the energy gap between the two levels. We choose the jump operators as $\sigma_{+}=\Ket{1}\Bra{0}$ and $\sigma_{-}=\sigma_{+}^{\dagger}$ with corresponding rates $\gamma n_{b}$ and $\gamma(1+n_{b})$, respectively. Here $n_{b}=1/(\exp(\hbar\omega_{0}/k_{B}T)-1)$ is the Bose-Einstein occupancy factor at bath temperature $T$, and $k_{B}$ is the Boltzmann constant. In this case, the eigenvalues $\lambda_{1}$ and $\lambda_{2}$ of the adjoint Liouvillian form a complex-conjugate pair ($\lambda_{2}=\lambda_{1}^{\star}$) corresponding to the slowest decay mode (with equal real parts). The associated right eigenmatrices $r_{1}$ and $r_{2}$ are complex conjugates of each other and are purely off-diagonal, a characteristic feature of the Davies map structure. Therefore, to suppress the contribution of these slowest decaying modes and thereby achieving accelerated relaxation in the operator dynamics, one can perform a mode-bypassing procedure on $\mathcal{O}$ such that it becomes diagonal. The transformed operator then takes the form $\widetilde{\mathcal{O}}=\begin{pmatrix}b_{0}+b_{z}&0\\ 0&b_{0}-b_{z}\end{pmatrix}$ where $b=\sqrt{b_{x}^{2}+b_{y}^{2}+b_{z}^{2}}$. As a result, the operator will quickly relax to the steady state. To investigate how the initial dressed distance relative to the steady state changes after the mode bypassing operation, we compute the dressed distance for the two cases. In our setup, the steady state is a detailed-balance thermal state, given by $\rho_{\mathrm{ss}}=e^{-\beta H}/\mathcal{Z}$, where $H=\frac{\omega_{0}}{2}\sigma_{z}$, $\beta=1/(k_{B}T)$ and the partition function $\mathcal{Z}=2\cosh{(\beta\omega_{0}/2)}$. The initial dressed distance can be written as $\mathcal{D}_{\mathrm{dd}}(\mathcal{O},\mathcal{O}_{\mathrm{ss}})=\sum_{i=1}^{2}|\lambda_{i}^{\prime}|=2\sqrt{b_{x}^{2}+b_{y}^{2}+(2b_{z}/\mathcal{Z})^{2}}/\mathcal{Z},$ where $\lambda_{i}^{\prime}$ are the eigenvalues of the operator $\rho_{\mathrm{ss}}^{1/2}(\mathcal{O}-\mathcal{O}_{\mathrm{ss}})\rho_{\mathrm{ss}}^{1/2}$. After the transformation, the distance between the transformed operator and the steady state becomes $\mathcal{D}_{\mathrm{dd}}(\widetilde{\mathcal{O}},\mathcal{O}_{\mathrm{ss}})=\sum_{i=1}^{2}|\widetilde{\lambda}_{i}^{\prime}|=4b_{z}/\mathcal{Z}^{2},$ which is always smaller than the dressed distance before the mode bypassing operation. Here, $\widetilde{\lambda}_{i}^{\prime}$ denote the eigenvalues of $\rho_{\mathrm{ss}}^{1/2}(\widetilde{\mathcal{O}}-\mathcal{O}_{\mathrm{ss}})\rho_{\mathrm{ss}}^{1/2}$. As a result, no Mpemba-like effect for operators can be observed in the single-qubit case with the chosen dissipators. It is interesting to note that, a similar absence of the Mpemba effect was observed in Ref. [62] where the accelerated relaxation for the operator was obtained via unitary transformation.

V Operator Mpemba effect in the intermediate dephasing regime for bosonic lattice

In this section, we analyze the dynamics of the one-dimensional bosonic lattice with bulk dephasing, introduced in Example 3 of the main text. We will discuss the regime in which the dephasing strength is comparable to the tunneling amplitude, $\gamma\sim J$ which implies that both coherent tunneling and dissipative processes compete with each other.

Unlike the strong dephasing limit, $\gamma\gg J$ discussed in the main text, the evolution cannot be restricted to the population sector alone in the regime $\gamma\sim J$. Instead, coherences in the site basis remain relevant and are coupled to the populations throughout the dynamics. As a result, the effective classical description in terms of the Pauli master equation is no longer valid, and the full Liouvillian dynamics must be considered.

To investigate the possibility of operator Mpemba effect in $\gamma\sim J$ regime, we consider the dynamics of the local observable $n_{i}=a_{i}^{\dagger}a_{i}$ and its mode-bypassed counterpart $\widetilde{n}_{i}$, constructed using Eq. (4). In Fig (S3)(a), we demonstrate the occurrence of the Mpemba effect for the operator where crossing is observed in the dynamics of dressed distance (it is proportional to trace distance, as the map here is unital). The corresponding expectation value of the transformed operator captures the faster relaxation (see inset). In Fig. (S3) (b), we plot the quantity $\Delta=\mathcal{D}_{\mathrm{dd}}(\widetilde{\mathcal{O}},\mathcal{O}_{\mathrm{ss}})-\mathcal{D}_{\mathrm{dd}}(\mathcal{O},\mathcal{O}_{\mathrm{ss}}),$ which is the difference in the distance value after and before the mode bypassing, as a function of $\gamma$. There is an overall enhancement in the distance measure in $\gamma\sim J$ regimes due to the coherence contribution in the dressed distance compared to the $\gamma\gg J$ region, where coherence is strongly suppressed.

Figure S3: Results for one-dimensional bosonic lattice subjected to dephasing in the $\gamma\sim J$ regime: (a) Plot for the dynamics of dressed distance (in this case it reduces to trace distance, as the map is unital) for the transformed operator $\widetilde{n}_{2}(t)$ which initially lies farther from the steady-state value $n_{2}^{(\mathrm{ss})}$, yet it relaxes more rapidly than the local number density of the second site $n_{2}$. As a result, genuine Mpemba effect for operators can be observed even in the $\gamma\sim J$ regime. The inset shows the accelerated convergence of the expectation value of the transformed operator compared to the original operator. (b) Plot for the quantity $\Delta=\mathcal{D}_{\mathrm{dd}}(\widetilde{\mathcal{O}},\mathcal{O}_{\mathrm{ss}})-\mathcal{D}_{\mathrm{dd}}(\mathcal{O},\mathcal{O}_{\mathrm{ss}}),$ which is the change in the distance after and before the mode bypassing, as a function of the dephasing strength $\gamma$. In the regime $\gamma\sim J$, both population and coherence in the site-basis remain coupled. Due to the coherence contribution, there is an overall enhancement in the distance measure after the mode-bypassing. As $\gamma$ increases, the quantity $\Delta$ decreases due to the negligible contribution from coherence. For the plots, we use the parameters $L=11$, $J=1$, and for (a) $\gamma=1$.

Figure S4: Numerical verification of the bounds (a) for the non-equilibrium XX spin chain model [Example 2 of the main text], and (b) for the bosonic lattice with dephasing [Example 3 of the main text]. For (a) we choose the initial operator $\mathcal{O}=\sum_{\alpha_{i}\in\{x,y,z\}}\bigotimes_{i=1}^{4}\sigma_{i}^{\alpha_{i}},$ and construct the transformed operator $\widetilde{\mathcal{O}}$ using the mode-bypassing procedure, where we bypass the SDM. We vary certain parameters of the setup randomly. For all such random realizations, as can be seen, the bounds are always satisfied. The numerical parameters are chosen as $g=0.05$, $\gamma_{1}=0.2$, $\gamma_{N}=0.5$, $\epsilon_{1,2,3,4}=(1,0.8,0.6,0.5)$, $T_{1}=1$, and $\mu_{1}=\mu_{N}=0$, while $T_{N}$ is varied randomly in the range $T_{N}\in[1.1,5]$. For (b), we take the population operator of the second site, $\mathcal{O}=n_{2}$, as the initial observable and construct the transformed operator by eliminating the SDM following the mode-bypassing procedure. The bounds are once again satisfied for all the random realizations. In this case, we use $J=1$ and vary the dephasing strength randomly in the range $\gamma\in[1,5]$.

VI Numerical Verification of the Bounds in Eq. (12)

In the main text, we provided a lower and upper bound on the quantity $\Delta=\mathcal{D}_{\rm dd}(\widetilde{\mathcal{O}},\mathcal{O}_{\rm ss})-\mathcal{D}_{\rm dd}(\mathcal{O},\mathcal{O}_{\rm ss}),$ which quantifies the change in the distance value after and before the mode-bypassing procedure. It satisfies the bounds $\Delta_{\rm lower}\leq\Delta\leq\Delta_{\rm upper}$, where $\Delta_{\rm lower}=-\mathrm{Re}\,\big[\mathrm{Tr}(U_{A}^{\dagger}B)\big]$ and $\Delta_{\rm upper}=-\mathrm{Re}\,\big[\mathrm{Tr}(U_{A-B}^{\dagger}B)\big]$, as given in Eq. (12). Recall that, $A=\mathcal{O}-\mathcal{O}_{\rm ss},B=\mathrm{Tr}(l_{1}\mathcal{O})\,r_{1},$ with $r_{1}$ and $l_{1}$ denoting the right and left slow-mode eigenoperators, respectively. In this section, we demonstrate the validity of both the lower and the upper bounds for example 2 and example 3, considered in the main text. In Fig. S4 (a) and (b) we demonstrate the validity of both lower and upper bounds for the non-equilibrium XX spin chain [Example 2 of the main text] and for the bosonic lattice with dephasing [Example 3 of the main text]. We plot results for 1000 random realizations by randomly varying temperature $T_{N}$ for (a) and dephasing strength $\gamma$ for (b). In Fig. S4 (a) and (b), we consider the initial operators $\mathcal{O}=\sum_{\alpha_{i}\in\{x,y,z\}}\bigotimes_{i=1}^{4}\sigma_{i}^{\alpha_{i}},$ and $\mathcal{O}=n_{2}$, respectively, both of which possess a nonzero overlap with the SDM of the Liouvillian $\mathcal{L}$. Using the mode-bypassing procedure, we construct the transformed operator $\widetilde{\mathcal{O}}$. To verify the bounds, we plot $\Delta-\Delta_{\rm lower}$ and $\Delta-\Delta_{\rm upper}$ for 1000 random realizations and observe the validity of Eq. (12) in all cases.