Stark Many-Body Localization Induced Quantum Mpemba Effect

1. Introduction

Understanding how isolated quantum many-body systems lose locally accessible information about their initial conditions is a central problem in non-equilibrium physics [1,2,2]. In generic interacting settings, the eigenstate thermalization hypothesis [3,4] suggests that after a global quench small subsystems approach stationary reduced states determined by conservation laws. When a global $U\left(1\right)$ charge Q is conserved, this relaxation has a particularly transparent structural signature: the reduced density matrix of a subregion, ${\rho }_A\left(t\right)$, progressively becomes compatible with the charge-sector decomposition set by the restricted charge $Q_A$, so that coherences between different $Q_A$ sectors are suppressed and local properties become consistent with an equilibrium description constrained by energy and charge [5,6].

This structural relaxation is commonly referred to as symmetry restoration [7,8,9,10,11,12,13,14,15] at the subsystem level. To diagnose this, the EA [16] is utilized: a symmetry-sensitive quantity constructed from ${\rho }_A\left(t\right)$ that measures the remaining locally accessible symmetry-breaking information and vanishes when ${\rho }_A\left(t\right)$ becomes symmetric with respect to $Q_A$. Because it is defined solely from the reduced state and directly targets inter-sector coherence, EA provides a experimental and theoretical proof of symmetry restoration in unitary dynamics. Symmetry restoration also exposes a counterintuitive ordering phenomenon. The QME (see recent reviews [17,18]) captures the possibility that, under identical dynamics, an initial state that is farther from the stationary structure can relax faster than one that starts closer. Quantum Mpemba behaviour has been explored in both open [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,47] and closed [16,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70] systems. Here, unless stated otherwise, we use the term, QME, exclusively in the closed-system sense, where the effect is diagnosed through local equilibration under unitary time evolution. Recent works [16,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70] have shown that this closed-system Mpemba phenomenology arises in Hamiltonian many-body dynamics and in symmetric random quantum circuits, where conservation laws are enforced at the level of local gates and the dynamics provides a minimal setting for unitary thermalization. In the symmetry-restoration context, the Mpemba effect corresponds to an inversion of restoration times, meaning that a more strongly symmetry-broken initial condition can restore symmetry faster according to symmetry-resolved diagnostics such as EA.

Thus, the key motivation of this work is to understand whether such symmetry-restoration Mpemba inversions persist when ergodicity is impeded. Localization [71] provides a natural route to non-ergodic unitary dynamics and long-lived memory, and it reshapes relaxation by restricting which parts of Hilbert space are dynamically explored on accessible time scales. This makes symmetry restoration strongly preparation dependent and suggests a natural mechanism for inversions [62]: different initial states can expand the effectively explored, symmetry-allowed Hilbert space at different rates, leading to different rates of sector redistribution and intra-sector mixing in ${\rho }_A\left(t\right)$, and hence to different decay rates of EA.

Many-body localization has already provided evidence that Mpemba-type inversions can occur even in non-ergodic dynamics [62]. This raises an immediate and practically important question: is the phenomenon tied to quenched randomness, or does it survive in non-ergodic settings driven by a different microscopic mechanism? Stark-MBL [72,73,74] offers a complementary route to localization that is induced by a uniform potential gradient rather than disorder. It avoids sample-to-sample fluctuations, reduces the reliance on disorder averaging, and introduces a continuously tunable control parameter, the field strength, that allows one to traverse ergodic and non-ergodic regimes within a single Hamiltonian family. These features make Stark-MBL an ideal platform to interrogate symmetry restoration and its Mpemba inversion in a controlled manner, and to connect the inversion directly to how a potential gradient reorganizes resonant processes and constrains the growth of the dynamically explored Hilbert space.

In this work, we have provided a comprehensive analysis of symmetry restoration and the QME within the context of Stark-MBL. By investigating the dynamics of the Rényi-2 EA in a tilted XXZ chain, we revealed that the Stark-induced breaking of ergodicity leads to a robust and universal Mpemba effect. Our central findings regarding the phase dependence of the QME are summarized below (see table in Figure 1):

• In the ergodic phase, the occurrence of the QME is contingent upon the microscopic details of the initial state. Specifically, it is present for the Tilted Ferromagnetic State (TFS) but absent for the Tilted Néel State (TNS), as the relaxation is governed by the dimension of the dynamically accessible charge sectors.
• In the Stark-MBL phase, the QME becomes a universal phenomenon. It emerges for both TNS and TFS initial states. This universality arises from a Stark-induced inversion of the steady-state ordering, where stronger initial symmetry breaking leads to a more symmetric local steady state due to the restricted effective Hilbert space dimension.

Furthermore, we analytically determined the characteristic timescale of this anomalous relaxation. We found that the Mpemba time $t_M$ is governed by high-order perturbative processes and scales exponentially with the subsystem size, $t_M\propto {(F/J)}^{N_A-1}$. This prediction is quantitatively corroborated by exact diagonalization simulations, confirming that symmetry restoration in the Stark-MBL phase is driven by slow, off-resonant quantum fluctuations rather than diffusive transport.

These results place Mpemba-type anomalies into a tunable localization setting without quenched disorder, offering a clear protocol to observe this effect in state-of-the-art quantum simulators such as trapped ions or superconducting circuits. Broadly, our work paves the way for utilizing tunable potentials to engineer anomalous relaxation timescales, providing new insights into the interplay between ergodicity breaking and symmetry-protected quantum information.

2. Methods

2.1. Model Hamiltonian

We study a one-dimensional spin-$1/2$ XXZ chain in the presence of a linear Stark potential. The Hamiltonian is given by:

$$H=-J\sum _{i=1}^{N-1}\left(S_i^x{S}_{i+1}^x+S_i^y{S}_{i+1}^y\right)+\Delta \sum _{i=1}^{N-1}{S}_i^z{S}_{i+1}^z+\sum _{i=1}^N[Fi+{(i/N)}^2]S_i^z,$$

(1)

where $S_i^{\alpha }$ ($\alpha \in \{x,y,z\}$) are the spin-$1/2$ operators at site i, satisfying the commutation relations $[S_i^{\alpha },S_j^{\beta }]=i{\delta }_{ij}{ϵ}_{\alpha \beta \gamma }{S}_i^{\gamma }$. Here, J sets the energy scale of the hopping set to $J=2$ hereafter, $\Delta$ represents the nearest-neighbor interaction anisotropy, and F denotes the strength of the effective electric field the Stark potential gradient. We employ open boundary conditions to accommodate the linear potential. The quadratic term ${(i/N)}^2$ is included as a smooth boundary regularization to gently lift finite-size near-degeneracies of the perfectly linear tilt and suppress long-lived coherent oscillations, thereby stabilizing the late-time dephasing behavior [72].

The model conserves the total magnetization $Q={\sum }_i{S}_i^z$, exhibiting a global $U\left(1\right)$ symmetry. The interplay between the interaction $\Delta$ and the tilt F governs the transition from an ergodic thermalizing phase to a Stark-MBL phase, where transport is suppressed, and the system breaks ergodicity.

2.2. Rényi-2 Entanglement Asymmetry

To quantify the restoration of the $U\left(1\right)$ symmetry within a local subsystem, we utilize the EA. We partition the system into a subsystem A of size $N_A$ and its complement B. The state of the subsystem is described by the reduced density matrix ${\rho }_A\left(t\right)={\mathrm{Tr}}_B\left[\left|\psi \right(t)\rangle \langle \psi (t\left)\right|\right]$. The global $U\left(1\right)$ symmetry implies the existence of a local charge operator $Q_A={\sum }_{i\in A}{S}_i^z$. If the symmetry is locally broken, $[{\rho }_A\left(t\right),Q_A]\ne 0$. The symmetrized state, which retains only the diagonal blocks in the charge basis, is obtained via the twirling map:

$${\rho }_{A,Q}\left(t\right)=\sum _q{\mathsf{\Pi }}_q{\rho }_A\left(t\right){\mathsf{\Pi }}_q={\displaystyle \frac{1}{2\pi }}{\int }_0^{2\pi }d\varphi e^{-i\varphi Q_A}{\rho }_A\left(t\right)e^{i\varphi Q_A},$$

(2)

where ${\mathsf{\Pi }}_q$ is the projector onto the eigenspace of $Q_A$ with charge q. The EA is defined as the difference between the entanglement entropy of the symmetrized state and the original state. We specifically employ the Rényi-2 EA, $\Delta S_A^{\left(2\right)}$, due to its analytical tractability:

$$\Delta S_A^{\left(2\right)}\left(t\right)=S^{\left(2\right)}\left({\rho }_{A,Q}\left(t\right)\right)-S^{\left(2\right)}\left({\rho }_A\left(t\right)\right),$$

(3)

where the Rényi-2 entropy is $S^{\left(2\right)}\left(\rho \right)=-log\mathrm{Tr}\left({\rho }^2\right)$. Substituting this definition, the asymmetry can be expressed as the logarithm of the ratio of purities:

$$\Delta S_A^{\left(2\right)}\left(t\right)=log\left({\displaystyle \frac{\mathrm{Tr}\left({\rho }_A^2\left(t\right)\right)}{\mathrm{Tr}\left({\rho }_{A,Q}^2\left(t\right)\right)}}\right).$$

(4)

This quantity is strictly non-negative, $\Delta S_A^{\left(2\right)}\left(t\right)\ge 0$, and vanishes if and only if the state is $U\left(1\right)$ symmetric, i.e., $[{\rho }_A,Q_A]=0$.

2.3. Operator Space Expansion

To analyze the dynamical mechanism of symmetry restoration, it is convenient to expand the density matrix in an orthonormal operator basis. The set of operator strings $P={\prod }_{i\in A}{\mathcal{O}}_i$, where ${\mathcal{O}}_i\in \{I_i,S_i^{+},S_i^{-},S_i^z\}$, forms a complete basis for operators on subsystem A with $S_i^{\pm }={\displaystyle \frac{1}{\sqrt{2}}}(S_i^x\pm i{S}_i^y)$. We normalize the basis such that $\mathrm{Tr}\left(P^{†}{P}^{\prime }\right)=2^{N_A}{\delta }_{P,P^{\prime }}$. The purity of the reduced density matrix can then be expanded as:

$$\mathrm{Tr}\left({\rho }_A^2\left(t\right)\right)={\displaystyle \frac{1}{2^{N_A}}}\sum _{P\subset A}{\left|\langle \psi \left(t\right)\left|P\right|\psi \left(t\right)\rangle \right|}^2\equiv {\displaystyle \frac{1}{2^{N_A}}}f\left(t\right),$$

(5)

where $f\left(t\right)$ represents the total weight of all operator strings.

Similarly, the purity of the symmetrized state ${\rho }_{A,Q}$ contains contributions only from charge-conserving strings, those satisfying $[P,Q_A]=0$:

$$\mathrm{Tr}\left({\rho }_{A,Q}^2\left(t\right)\right)={\displaystyle \frac{1}{2^{N_A}}}\sum _{P\subset A:[P,Q_A]=0}{\left|\langle \psi \left(t\right)\left|P\right|\psi \left(t\right)\rangle \right|}^2\equiv {\displaystyle \frac{1}{2^{N_A}}}{f}^c\left(t\right).$$

(6)

Substituting these into Eq. (4), the Rényi-2 EA, then, becomes:

$$\Delta S_A^{\left(2\right)}\left(t\right)=log\left({\displaystyle \frac{f\left(t\right)}{f^c\left(t\right)}}\right).$$

(7)

Here, $f\left(t\right)$ sums over all strings, while $f^c\left(t\right)$ sums only over the diagonal strings that preserve the local charge.

The dynamics of $\Delta S_A^{\left(2\right)}\left(t\right)$, Eq. (7), are thus governed by the relative decay rates of charge-neutral versus charge-non-conserving operator strings.

2.4. Quantum Mpemba Effect

The QME describes an anomalous relaxation phenomenon where a system initially further from equilibrium relaxes faster than one initially closer. We consider a family of initial states parametrized by a tilt angle $\theta$, denoted as $|{\psi }_{\theta }\rangle$. The distance from symmetry at time t is quantified by $\Delta S_A^{\left(2\right)}(t;\theta )$. We assume that the initial asymmetry $\Delta S_A^{\left(2\right)}(0;\theta )$ is monotonically increasing with $\theta \in (0,\frac{\pi }{2}]$, i.e., “hotter” states have larger initial asymmetry.

A QME is defined to occur iff there exist two parameters ${\theta }_h>{\theta }_c$ such that:

$$\Delta S_A^{\left(2\right)}(0;{\theta }_h)>\Delta S_A^{\left(2\right)}(0;{\theta }_c),$$

(8)

but at some later time $t>t_M$ (the Mpemba time):

$$\Delta S_A^{\left(2\right)}(t;{\theta }_h)<\Delta S_A^{\left(2\right)}(t;{\theta }_c).$$

(9)

Specifically, we investigate the existence of this in the Stark-MBL regime, where non-ergodicity may protect certain operator sectors, leading to non-trivial relaxation hierarchies.

3. Results

3.1. Analytical Results

3.1.1. From the Stark Hamiltonian to an Effective Diagonal Description

We start from Eq. (1) and focus on the strong-tilt regime $F\gg J,\Delta$. We decompose $H=H_0+V$ with

$$\begin{array}{cc}\hfill H_0& =\sum _{i=1}^N\left[Fi+{(i/N)}^2\right]S_i^z+\Delta \sum _{i=1}^{N-1}{S}_i^z{S}_{i+1}^z,\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill V& =-J\sum _{i=1}^{N-1}\left(S_i^x{S}_{i+1}^x+S_i^y{S}_{i+1}^y\right)=-{\displaystyle \frac{J}{2}}\sum _{i=1}^{N-1}\left(S_i^{+}{S}_{i+1}^{-}+S_i^{-}{S}_{i+1}^{+}\right).\hfill \end{array}$$

(11)

In the Dirac picture with respect to $H_0$, each nearest-neighbor flip-flop term acquires a rapidly oscillating phase $e^{\pm i[({\epsilon }_{i+1}-{\epsilon }_i)+\Delta {\delta }_i]t}$, where ${\epsilon }_i=Fi+{(i/N)}^2$ and ${\delta }_i$ is an $O\left(1\right)$ operator-valued shift set by the adjacent $S^z$ configuration. Since ${\epsilon }_{i+1}-{\epsilon }_i\simeq F$ dominates for $F\gg J,\Delta$, these terms are strongly off-resonant and average out on the long-time scales relevant for the late-time envelopes discussed below.

Consequently, to leading order in $J/F$, the dynamics is governed by an effective Hamiltonian obtained by perturbatively eliminating V (equivalently via a Schrieffer–Wolff transformation). In the strong-Stark limit this effective Hamiltonian is (quasi-)diagonal in the computational basis and takes the generic form

$$H_{\mathrm{eff}}=\sum _i{h}_i{S}_i^z+\sum _{i<j}{J}_{ij}{S}_i^z{S}_j^z+\dots ,$$

(12)

where $h_i\approx Fi+{(i/N)}^2$ up to $O(J^2/F)$ corrections, and the longitudinal couplings $J_{ij}$ arise from virtual sequences of the flip-flop processes and decay rapidly with distance. For our purposes the crucial point is that $H_{\mathrm{eff}}$ generates configuration-dependent phases in the $S^z$ basis, which dephase expectation values of operator strings containing $S^{\pm }$ at late times.

3.1.2. Stark Many-Body Localized EA Plateau

We now derive the late-time plateau rules for the operator-string weights entering the Rényi-2 EA (7). We consider the same tilted product initial-state families as in the main text, namely the TFS and the TNS. Let $\left\{\right|0\rangle ,|1\rangle \}$ denote the local eigenstates of computational basis $S_i^z$. The TFS is the uniform tilted product state

$$|{\psi }_{\theta }^{\mathrm{TFS}}\rangle ={\left(cos\left(\frac{\theta }{2}\right)|0\rangle +sin\left(\frac{\theta }{2}\right)|1\rangle \right)}^{\otimes N}.$$

(13)

The TNS differs only by a staggered relative phase between the $|0\rangle$ and $|1\rangle$ components,

$$|{\psi }_{\theta }^{\mathrm{TNS}}\rangle =\underset{i=1}{\overset{N}{⨂}}\left(cos\left(\frac{\theta }{2}\right)|0\rangle +{(-1)}^{i}sin\left(\frac{\theta }{2}\right)|1\rangle \right).$$

(14)

In both cases the computational-basis probabilities are identical and site independent, $|\langle \sigma |{\psi }_{\theta }\rangle {|}^2={\prod }_{i=1}^N{p}_{{\sigma }_i}$ with $p_0={cos}^2(\theta /2)$ and $p_1={sin}^2(\theta /2)$.

In the effective diagonal dynamics of Section 3.1.1, operator strings containing at least one $S^{\pm }$ acquire rapidly varying configuration-dependent phases. After long-time averaging, off-diagonal contributions dephase and one obtains a diagonal-probability expression whose magnitude factorizes into single-site contributions. Each site on which P contains $S^{\pm }$ enforces a flip between the paired computational configurations and contributes a factor $2p_0{p}_1$, whereas each site with I or $S^z$ contributes $p_0^2+p_1^2$. It is therefore convenient to define

$$\begin{array}{cc}\hfill a_{\theta }& \equiv p_0^2+p_1^2={cos}^4(\theta /2)+{sin}^4(\theta /2)={\displaystyle \frac{1+{cos}^2\theta }{2}},\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill b_{\theta }& \equiv 2p_0{p}_1=2{cos}^2(\theta /2){sin}^2(\theta /2)={\displaystyle \frac{{sin}^2\theta }{2}}.\hfill \end{array}$$

(16)

Then for any operator string with $n_{\mathrm{chg}}\left(P\right)=\#\left(S^{+}\right)+\#\left(S^{-}\right)\ge 1$, the late-time plateau obeys

$$\overline{{\left|\langle \psi \left(t\right)\left|P\right|\psi \left(t\right)\rangle \right|}^2}\simeq b_{\theta }^{n_{\mathrm{chg}}\left(P\right)}{a}_{\theta }^{N-n_{\mathrm{chg}}\left(P\right)}.$$

(17)

The factor $a_{\theta }^{N-n_{\mathrm{chg}}\left(P\right)}$ includes contributions from sites outside A because the long-time dephasing expression involves full-system probabilities $|\langle \sigma |{\psi }_{\theta }\rangle {|}^4$. For purely diagonal strings $P\in {\{I,S^z\}}^{\otimes N_A}$ one has $[P,H_{\mathrm{eff}}]=0$, so their expectation values are time independent. Summing their squared expectations in the above operator normalization gives the exact identity

$$\sum _{P\in {\{I,S^z\}}^{\otimes N_A}}{\left|\langle \psi \left(t\right)\left|P\right|\psi \left(t\right)\rangle \right|}^2={(1+{cos}^2\theta )}^{N_A}=2^{N_A}{a}_{\theta }^{N_A}.$$

(18)

We now compute the long-time limits of $f\left(t\right)$ and $f^c\left(t\right)$ defined in Section 2. A convenient way to obtain $f(\theta ,\infty )={\sum }_{P\subset A}\overline{{\left|\langle P\left(t\right)\rangle \right|}^2}$ is: (i) apply the plateau rule Eq. (17) to all $4^{N_A}$ strings, and (ii) replace the purely diagonal subset $P\in {\{I,S^z\}}^{\otimes N_A}$ by its exact contribution Eq. (18).

Step (i): summing Eq. (17) over all strings amounts to summing, on each site in A, the four local contributions I, $S^z$, $S^{+}$, $S^{-}$, which yields

$$(a_{\theta }+a_{\theta }+b_{\theta }+b_{\theta })=2(a_{\theta }+b_{\theta })=2,$$

(19)

since $a_{\theta }+b_{\theta }={(p_0+p_1)}^2=1$. Thus the plateau sum over all strings is $2^{N_A}{a}_{\theta }^{N-N_A}$.

Step (ii): within that plateau sum, the diagonal subset contributes $2^{N_A}{a}_{\theta }^N$, while its exact value is $2^{N_A}{a}_{\theta }^{N_A}$. Replacing the diagonal subset therefore gives

$$f(\theta ,\infty )=2^{N_A}\left(a_{\theta }^{N_A}+a_{\theta }^{N-N_A}-a_{\theta }^N\right).$$

(20)

For $f^c(\theta ,\infty )$, by definition only charge-neutral strings with $[P,Q_A]=0$ contribute. All diagonal strings are charge-neutral, hence $f^c(\theta ,\infty )\ge 2^{N_A}{a}_{\theta }^{N_A}$. Additional charge-neutral contributions arise from strings containing equal numbers of $S^{+}$ and $S^{-}$ within A. These form a strict subset of the spin-changing strings and their total contribution is bounded by the total spin-changing weight, giving

$$0\le f^c(\theta ,\infty )-2^{N_A}{a}_{\theta }^{N_A}\le 2^{N_A}\left(a_{\theta }^{N-N_A}-a_{\theta }^N\right).$$

(21)

In finite N with $\theta \in (0,\frac{\pi }{2}]$, so that $a_{\theta }<1$ and $N_A/N<1/2$, this correction is parametrically small compared with $2^{N_A}{a}_{\theta }^{N_A}$. To leading order we therefore take

$$f^c(\theta ,\infty )\simeq 2^{N_A}{a}_{\theta }^{N_A}.$$

(22)

Thus, the steady-state plateau prediction $\Delta S_A^{\left(2\right)}(\infty )=log{\displaystyle \frac{f(\theta ,\infty )}{f^c(\theta ,\infty )}}\simeq log\left[1+a_{\theta }^{N-2N_A}-a_{\theta }^{N-N_A}\right]$. When $a_{\theta }^{N-2N_A}\ll 1$, we obtain

$$\begin{array}{c}\hfill \Delta S_A^{\left(2\right)}(\infty )\simeq a_{\theta }^{N-2N_A}\left(1-a_{\theta }^{N_A}\right).\end{array}$$

(23)

Note that the above plateau rule depends on the initial state only through the computational-basis probabilities $|\langle \sigma |{\psi }_{\theta }\rangle {|}^2$. Since TFS and TNS differ only by site-dependent phases, they share the same $|\langle \sigma |{\psi }_{\theta }\rangle {|}^2$ and hence yield identical plateau predictions for $f(\infty )$, $f^c(\infty )$, and $\Delta S_A^{\left(2\right)}(\infty )$.

3.1.3. Emergent QME in Stark Many-Body Localization and Mpemba-Time Scaling

We now elucidate the physical mechanism underlying the emergent QME and its Rényi-2 EA QME time scaling properties.

First, we address the question of symmetry restoration. In the Stark-MBL regime, the system is non-ergodic and does not relax to a thermal Gibbs ensemble. However, this does not imply that local subsystems retain all information about the initial symmetry breaking. As shown by Eqs. (20) and (22), in the thermodynamic limit $N\to \infty$ with fixed subsystem size $N_A$, the steady-state EA is controlled by the term $a_{\theta }^{N-2N_A}$. Since $a_{\theta }<1$ for any $\theta >0$, this term vanishes as $N\to \infty$. Consequently, $\Delta S_A^{\left(2\right)}(\infty )\to 0$. This indicates that even in the absence of global transport, the dephasing induced by the Stark potential is sufficient to destroy off-diagonal coherences between charge sectors locally, thereby restoring the $U\left(1\right)$ symmetry at the subsystem level.

Second, we determine the existence of the QME and its physical origin. The Stark-MBL regime is governed by a fundamentally different mechanism than the ergodic regime, leading to distinct behaviors for the TNS [also see Figure 1]:

• Ergodic Regime: Symmetry restoration is governed by the dimension of the dynamically accessible Hilbert space sectors.

  -

  For the TNS, the untilted state ($\theta =0$) resides in the half-filling sector ($Q=0$), which has the maximal Hilbert space dimension and thus the fastest thermalization rate. Increasing the tilt $\theta$ spreads the wavefunction into charge sectors away from half-filling, which have smaller dimensions and slower relaxation rates. Since the state with larger initial asymmetry relaxes slower, the EA curves do not cross, and QME is absent for the TNS.

  -

  For the TFS, the situation is reversed. The un-tilted state resides at the edge of the spectrum. Increasing $\theta$ moves the state toward the central, high-dimension sectors, accelerating relaxation. Thus, QME is present.

• Stark-MBL Regime: Relaxation is governed by dephasing in the effective diagonal basis, and the existence of QME is determined by the ordering of the steady-state values.

  -

  As derived in Eq. (23), the steady-state plateau $\Delta S_A^{\left(2\right)}(\infty )$ scales with $a_{\theta }^{N-2N_A}$. Since $a_{\theta }=(1+{cos}^2\theta )/2$ is a monotonically decreasing function of $\theta$, a larger initial tilt $\theta$ (higher initial asymmetry) leads to a lower steady-state asymmetry.

  -

  This inversion of ordering, high initial value leading to low final value, mathematically guarantees that the EA evolution curves for different $\theta$ must cross at intermediate times. Consequently, QME becomes universal in the Stark-MBL regime, appearing for both TNS and TFS regardless of their ergodic behavior.

Finally, we estimate the scaling of the Mpemba time $t_M$. In the strong Stark limit ($F\gg J$), restoring symmetry within a subsystem of size $N_A$ requires quantum fluctuations that couple different local charge sectors. In the effective Hamiltonian picture, such charge reconfiguration is an off-resonant process. Changing the charge configuration across $N_A$ sites requires virtually occupying intermediate states up the Stark ladder. According to high-order perturbation theory, the effective matrix element for a process spanning $N_A$ sites scales as $J_{\mathrm{eff}}\sim J{(J/F)}^{N_A-1}$. The timescale for symmetry restoration, and consequently the time scale on which the QME occurs, is inversely proportional to this rate:

$$t_M\propto {\displaystyle \frac{1}{J_{\mathrm{eff}}}}\sim {\displaystyle \frac{1}{J}}{\left({\displaystyle \frac{F}{J}}\right)}^{N_A-1}.$$

(24)

This prediction implies that the Mpemba time grows exponentially with the subsystem size $N_A$ and follows a power-law dependence on the Stark field strength F. This confirms that the phenomenon is controlled by the localization length of the system, consistent with the numerical observations in Figure 3(b) and Figure 4(d).

3.2. Numerical Results

3.2.1. Stark Many-Body Localized Diagnostics

We first establish the parameter regime in which the tilted XXZ chain, described by the Hamiltonian in Eq. (1), enters the Stark-MBL phase. To obtain reliable information regarding the energy spectrum and eigenstates, we numerically solve the static and dynamical properties of the system using the Exact Diagonalization method. Figure 2(a) shows the half-chain von Neumann entropy density $\overline{S_{N/2}}/N$ as a function of the Stark field F for several system sizes N. In the small-F regime, $\overline{S_{N/2}}/N$ remains finite and exhibits only a mild size dependence, consistent with the volume-law behavior of the ergodic phase. As F increases, $\overline{S_{N/2}}/N$ decreases substantially and drifts downward with increasing N, indicating a crossover toward the area-law entanglement characteristic of Stark-MBL. The inset of Figure 2(a) demonstrates a finite-size scaling collapse as a function of $(F-F_c)N^{1/\nu }$, yielding a critical field $F_c\simeq 4.0$ (indicated by the vertical green dotted line). Based on this result, we use $F=15>F_c$ as a representative Stark-MBL point in our subsequent analysis.

To further obtain dynamical evidence of non-ergodicity, we examine the charge imbalance $\mathrm{CI}\left(t\right)$ starting from the initial TNS. Figure 2(b) plots the dynamics of $\mathrm{CI}\left(t\right)$ at $N=14$ and $F=15$ for several tilt angles $\theta$. After an early-time transient, $\mathrm{CI}\left(t\right)$ rapidly approaches a non-zero long-time plateau. Crucially, the plateau value remains strongly dependent on the initial angle $\theta$, preserving the ordering of the initial conditions over extremely long time scales (up to $t\sim {10}^{15}$). This persistent memory of the initial staggered pattern provides direct dynamical evidence for the non-ergodicity protected by the Stark-MBL phase.

3.2.2. EA dynamics and the emergence of QME

Having established the phase diagram, we turn to the core objective: investigating the dynamics of the Rényi-2 EA $\Delta S_A^{\left(2\right)}\left(t\right)$ and the resulting QME. To explore how the initial state distribution affects the symmetry restoration process, we first analyze the microscopic composition of the initial states by calculating the charge-sector weights $p_A\left(q_A\right)=\mathrm{Tr}\left({\mathsf{\Pi }}_{q_A}{\rho }_A\right)$, where ${\mathsf{\Pi }}_{q_A}$ is the projection operator onto the eigenspace of the local charge operator $Q_A={\sum }_{j\in A}{S}_j^z$ with eigenvalue $q_A$. To further characterize the spatial memory of the initial configurations, we calculate the charge imbalance (CI), which for our spin-$\frac{1}{2}$ system is defined as

$$\mathrm{CI}\left(t\right)={\displaystyle \frac{1}{N}}\sum _{i=1}^N{(-1)}^i\langle S_i^z\left(t\right)\rangle .$$

(25)

Figure 3(a) and Figure 3(b) show the EA dynamics for the TNS at $F=2$ and $F=15$, respectively. At the ergodic point $F=2$ [Figure 3(a)], while larger $\theta$ corresponds to a larger initial EA, the ordering of the curves does not invert during relaxation, and the QME is absent. The inset shows that $p_A\left(q_A\right)$ for the TNS is nearly insensitive to $\theta$. However, deep in the Stark-MBL regime ($F=15$), the EA curves exhibit clears, signaling the emergence of the QME. In contrast, for the TFS at $F=2$ [Figure 4(a)], the QME is already observable, which is highly consistent with the microscopic feature shown in its inset where $p_A\left(q_A\right)$ significantly shifts toward the central sectors as $\theta$ increases.

Then, to verify the accuracy of our aforementioned analytical predictions, we perform a quantitative comparison between the long-time numerical results and the analytical plateau values. In Figure 3(c) (TNS) and Figure 4(c) (TFS), we plot the steady-state EA $\Delta S_A^{\left(2\right)}\left(t_{\mathrm{ss}}\right)$ versus $\theta /\pi$ for $F=15$. The discrete symbols represent the numerical results averaged over 2000 samples, while the solid lines denote the analytical values. The results show that the numerical simulations and analytical predictions are in excellent agreement across all subsystem sizes $N_A$. This provides strong evidence that our analytical framework quantitatively reproduces the long-time dephasing behavior of Stark-MBL systems under symmetry constraints.

3.2.3. Scaling characteristics of the Mpemba time

Finally, we quantify the characteristic timescale of the QME. We define the Mpemba time $t_M$ (corresponding to the first intersection of the EA dynamical envelopes) by extracting the first point of the EA curves and fit the numerical results using the analytically derived scaling relation Eq. (24). As shown in Figure 3(d) and Figure 4(d), the fitting curves (dashed lines) perfectly capture the exponential growth of $t_M$ with respect to both subsystem size $N_A$ and field strength F. This not only further validates the theoretical predictions but also reveals that in the Stark-MBL environment, the QME timescale is dominated by suppressed high-order operator dynamics, where the effective transition rate $J_{\mathrm{eff}}$ decays as ${(J/F)}^{N_A-1}$.

4. Conclusions

To sum up, in this work, we have provided a comprehensive analysis of symmetry restoration and the QME within the context of Stark-MBL. By investigating the dynamics of the Rényi-2 EA in a tilted XXZ chain, we revealed that the Stark-induced breaking of ergodicity leads to a robust and universal Mpemba effect, distinct from its behavior in thermalizing systems.

Our central finding is that while the occurrence of the QME in the ergodic phase is contingent upon the specific choice of the initial state (e.g., present for TFS but absent for TNS), the Stark-MBL regime enforces a universal inversion of relaxation hierarchies. Through an analytical operator-string expansion, we demonstrated that this universality stems from the structure of the steady-state EA. In the strong-tilt limit, the effective Hilbert space dimension accessible to local observables is strictly constrained by the Stark potential, leading to a “hotter-is-faster” relaxation scenario where states with larger initial asymmetry ultimately approach a more symmetric steady state. Furthermore, we established that the timescale of this anomalous relaxation, the Mpemba time $t_M$, is governed by high-order off-resonant processes. We derived a scaling law showing that $t_M$ grows exponentially with the subsystem size and follows a power-law dependence on the field strength. This prediction is fully supported by our numerical simulations, confirming that symmetry restoration in Stark-MBL is driven by slow, perturbative quantum fluctuations rather than diffusive transport. Our results highlight the unique potential of Stark-MBL systems as tunable platforms for exploring non-equilibrium thermodynamics and symmetry dynamics. Unlike disorder-driven MBL, the Stark potential offers a clean, continuously adjustable control parameter to engineer relaxation pathways. This suggests promising avenues for observing the QME in present-day quantum simulators, such as cold atoms [75,76,77,78,79,80,81] in optical lattices or trapped ions [82], where linear potentials can be implemented with high precision.

Finally, our study opens up several specific directions for future theoretical investigation. While we focused on closed quantum systems, examining the stability of the Mpemba effect under dissipative or non-Hermitian dynamics [83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98] would be crucial for understanding its fate in realistic open environments. Additionally, extending the analysis to the framework of nonstabilizerness [99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154] offers a fresh perspective. Utilizing metrics such as the stabilizer Rényi entropy as dynamical probes could quantify the complexity of the relaxation process, potentially revealing a rich interplay between symmetry restoration, localization, and the scrambling of quantum magic.