Stark Many-Body Localization-Induced Quantum Mpemba Effect

TL;DR

This paper explores a quantum phenomenon where systems far from equilibrium relax faster, and finds it becomes a universal feature in certain disordered systems.

Contribution

The study reveals the quantum Mpemba effect becomes universal in Stark many-body localized systems.

Findings

The quantum Mpemba effect transitions from initial-state dependence to universality in Stark-MBL systems.

Mpemba time scales exponentially with subsystem size due to high-order off-resonant processes.

Stark-induced relaxation channels constrain the effective Hilbert space dimension.

Abstract

The quantum Mpemba effect (QME) describes the counterintuitive phenomenon where a system initially further from equilibrium relaxes faster than one closer to it. Specifically, the QME associated with symmetry restoration has been extensively investigated across integrable, ergodic, and disordered localized systems. However, its fate in disorder-free ergodicity-breaking settings, such as the Stark many-body localized (Stark-MBL) phase, remains an open question. Here, we explore the dynamics of local U(1) symmetry restoration in a Stark-MBL XXZ spin-12 chain, using the Rényi-2 entanglement asymmetry (EA) as a probe. Using an analytical operator-string expansion supported by numerical simulations, we demonstrate that the QME transitions from an initial-state-dependent anomaly in the ergodic phase to a universal feature in the Stark-MBL regime. Moreover, the Mpemba time scales exponentially with the subsystem size, even in the absence of global transport, and is governed by high-order off-resonant processes. We attribute this robust inversion to a Stark-induced hierarchy of relaxation channels that fundamentally constrains the effective Hilbert space dimension. The findings pave the way for utilizing tunable potentials to engineer and control anomalous relaxation timescales in quantum technologies without reliance on quenched disorder.

Figures

(a,b) EA dynamics with N=14 and NA=3. The initial state is chosen as TNS. The numerical curves have been smoothed by convolution with a Gaussian, ω(n)=e−(n/σ)2, with σ=5. The shaded curves represent the original numerical data before smoothing. Inset in (a,b): Charge-sector weights pA(qA) with N=14 and NA=4, for the initial TNS. (c) Steady-state EA for F=15>Fc and different subsystem sizes NA, with the initial state prepared as the TNS. Discrete points show the steady-state EA averaged over 2000 samples, while solid lines denote the analytical steady-state values. Inset in (c): analytical initial EA at t=0. (d) QME time tM as a function of subsystem size NA for different F at N=14, starting from the initial TNS. Discrete symbols (with error bars) show the numerically extracted tM, while dashed lines denote the theoretical scaling curves obtained by fitting Equation (24) to the numerical data.

(a,b) EA dynamics with N=14 and NA=3. The initial state is chosen as TFS. The numerical curves have been smoothed by convolution with a Gaussian, ω(n)=e−(n/σ)2, with σ=5. The shaded curves represent the original numerical data before smoothing. Inset in (a,b): Charge-sector weights pA(qA) with N=14 and NA=4, for the initial TFS. (c) Steady-state EA versus θ for F=15>Fc and different subsystem sizes NA, with the initial state prepared as the TFS. Discrete points show the steady-state EA averaged over late-time 2000 samples, while solid lines denote the analytical steady-state values. Inset in (c): analytical initial EA at t = 0. (d) QME time tM as a function of subsystem size NA for different F at N = 14, starting from the initial TFS. Discrete symbols (with error bars) show the numerically extracted tM, while dashed lines denote the theoretical scaling curves obtained by fitting Equation (24) to the numerical data.

Full text

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]. 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 $[eqn]$ charge Q is conserved, this relaxation has a particularly transparent structural signature: the reduced density matrix of a subregion, $[eqn]$ , progressively becomes compatible with the charge-sector decomposition set by the restricted charge $[eqn]$ , so that coherences between different $[eqn]$ 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 $[eqn]$ that measures the remaining locally accessible symmetry-breaking information and vanishes when $[eqn]$ becomes symmetric with respect to $[eqn]$ . Because it is defined solely in terms of the reduced state and directly targets inter-sector coherence, EA provides an 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] and closed [16,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65] 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,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65] 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 [66] 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 [58]: 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 $[eqn]$ , 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 [58]. This raises an immediate and practically important question: Is the phenomenon tied to quenched randomness, or does it persist in non-ergodic settings driven by a different microscopic mechanism? Stark-MBL [67,68,69] 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 for interrogating symmetry restoration and its Mpemba inversion in a controlled manner, and for connecting 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 provide a comprehensive analysis of symmetry restoration and the QME in the context of Stark-MBL. By investigating the dynamics of the Rényi-2 EA in a tilted XXZ chain, we found that the Stark-induced breaking of ergodicity gives rise to a robust and universal Mpemba effect. Our central findings regarding the phase dependence of the QME are summarized below (see the 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 $[eqn]$ is governed by high-order perturbative processes and scales exponentially with the subsystem size, $[eqn]$ . 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 in a tunable localization setting without quenched disorder, providing a clear protocol for observing 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- $[eqn]$ XXZ chain in the presence of a linear Stark potential. The Hamiltonian is given by:

$$[eqn]$$

where $[eqn]$ ( $[eqn]$ ) are the spin- $[eqn]$ operators at site i, satisfying the commutation relations $[eqn]$ . Here, J sets the energy scale of the hopping set to $[eqn]$ thereafter; $[eqn]$ 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 $[eqn]$ 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 [67].

The model conserves the total magnetization $[eqn]$ , exhibiting a global $[eqn]$ symmetry. The interplay between the interaction $[eqn]$ and the tilt F governs the transition from an ergodic thermalizing phase to a Stark-MBL phase, in which transport is suppressed, and the system breaks ergodicity.

2.2. R Ényi-2 Entanglement Asymmetry

To quantify the restoration of the $[eqn]$ symmetry within a local subsystem, we utilize the EA. We partition the system into a subsystem A of size $[eqn]$ and its complement B. The state of the subsystem is described by the reduced density matrix $[eqn]$ . The global $[eqn]$ symmetry implies the existence of a local charge operator $[eqn]$ . If the symmetry is locally broken, $[eqn]$ . The symmetrized state, which retains only the diagonal blocks in the charge basis, is obtained via the twirling map:

$$[eqn]$$

where $[eqn]$ is the projector onto the eigenspace of $[eqn]$ 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, $[eqn]$ , due to its analytical tractability:

$$[eqn]$$

where the Rényi-2 entropy is $[eqn]$ . Substituting this definition, the asymmetry can be expressed as the logarithm of the ratio of purities:

$$[eqn]$$

This quantity is strictly non-negative, $[eqn]$ , and vanishes if and only if the state is $[eqn]$ symmetric, i.e., $[eqn]$ .

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 $[eqn]$ , where $[eqn]$ , forms a complete basis for operators on subsystem A with $[eqn]$ . We normalize the basis such that $[eqn]$ . The purity of the reduced density matrix can then be expanded as follows:

$$[eqn]$$

where $[eqn]$ represents the total weight of all operator strings.

Similarly, the purity of the symmetrized state $[eqn]$ contains contributions only from charge-conserving strings, those satisfying $[eqn]$ :

$$[eqn]$$

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

$$[eqn]$$

Here, $[eqn]$ sums over all strings, while $[eqn]$ sums only over the diagonal strings that preserve the local charge.

The dynamics of $[eqn]$ , Equation (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 in which a system that is initially farther from equilibrium relaxes faster than one that is initially closer. We consider a family of initial states parametrized by a tilt angle $[eqn]$ , denoted as $[eqn]$ . The distance from symmetry at time t is quantified by $[eqn]$ . We assume that the initial asymmetry $[eqn]$ is monotonically increasing with $[eqn]$ , i.e., “hotter” states have larger initial asymmetry.

A QME is defined to occur iff there exist two parameters $[eqn]$ such that

$$[eqn]$$

but at some later time $[eqn]$ (the Mpemba time):

$$[eqn]$$

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 Equation (1) and focus on the strong-tilt regime $[eqn]$ . We decompose $[eqn]$ with

$$[eqn]$$ $$[eqn]$$

In the Dirac picture with respect to $[eqn]$ , each nearest-neighbor flip-flop term acquires a rapidly oscillating phase $[eqn]$ , where $[eqn]$ and $[eqn]$ is an $[eqn]$ operator-valued shift set by the adjacent $[eqn]$ configuration. Since $[eqn]$ dominates for $[eqn]$ , 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 $[eqn]$ , 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

$$[eqn]$$

where $[eqn]$ up to $[eqn]$ corrections, and the longitudinal couplings $[eqn]$ arise from virtual sequences of the flip-flop processes and decay rapidly with distance. For our purposes, the crucial point is that $[eqn]$ generates configuration-dependent phases in the $[eqn]$ basis, which dephase expectation values of operator strings containing $[eqn]$ 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 $[eqn]$ denote the local eigenstates of the computational basis $[eqn]$ . The TFS is the uniform tilted product state

$$[eqn]$$

The TNS differs only by a staggered relative phase between the $[eqn]$ and $[eqn]$ components,

$$[eqn]$$

In both cases, the computational-basis probabilities are identical and site-independent, $[eqn]$ with $[eqn]$ and $[eqn]$ .

In the effective diagonal dynamics of Section 3.1.1, operator strings containing at least one $[eqn]$ acquire rapidly varying configuration-dependent phases. After long-time averaging, off-diagonal contributions dephase, yielding a diagonal-probability expression whose magnitude factorizes into single-site contributions. Each site on which P contains $[eqn]$ enforces a flip between the paired computational configurations and contributes a factor $[eqn]$ , whereas each site with I or $[eqn]$ contributes $[eqn]$ . It is therefore convenient to define

$$[eqn]$$ $$[eqn]$$

Then for any operator string with $[eqn]$ , the late-time plateau obeys

$$[eqn]$$

The factor $[eqn]$ includes contributions from sites outside A because the long-time dephasing expression involves full-system probabilities $[eqn]$ . For purely diagonal strings $[eqn]$ , one has $[eqn]$ , so their expectation values are time independent. Summing their squared expectations in the above operator normalization gives the exact identity

$$[eqn]$$

We now compute the long-time limits of $[eqn]$ and $[eqn]$ defined in Section 2. A convenient way to obtain $[eqn]$ is to (i) apply the plateau rule Equation (17) to all $[eqn]$ strings and (ii) to replace the purely diagonal subset $[eqn]$ by its exact contribution Equation (18).

Step (i): Summing Equation (17) over all strings amounts to summing, on each site in A, the four local contributions I, $[eqn]$ , $[eqn]$ , $[eqn]$ , which yields

$$[eqn]$$

since $[eqn]$ . Thus, the plateau sum over all strings is $[eqn]$ .

Step (ii): Within that plateau sum, the diagonal subset contributes $[eqn]$ , while its exact value is $[eqn]$ . Replacing the diagonal subset, therefore, gives

$$[eqn]$$

For $[eqn]$ , by definition, only charge-neutral strings with $[eqn]$ contribute. All diagonal strings are charge-neutral, hence $[eqn]$ . Additional charge-neutral contributions arise from strings containing equal numbers of $[eqn]$ and $[eqn]$ 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

$$[eqn]$$

In finite N with $[eqn]$ so that $[eqn]$ and $[eqn]$ , this correction is parametrically small compared with $[eqn]$ . To leading order we therefore take

$$[eqn]$$

Thus, the steady-state plateau prediction $[eqn]$ . When $[eqn]$ , we obtain

$$[eqn]$$

Note that the above plateau rule depends on the initial state only through the computational-basis probabilities $[eqn]$ . Since TFS and TNS differ only by site-dependent phases, they share the same $[eqn]$ and hence yield identical plateau predictions for $[eqn]$ , $[eqn]$ , and $[eqn]$ .

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 Equations (20) and (22), in the thermodynamic limit $[eqn]$ with fixed subsystem size $[eqn]$ , the steady-state EA is controlled by the term $[eqn]$ . Since $[eqn]$ for any $[eqn]$ , this term vanishes as $[eqn]$ . Consequently, $[eqn]$ . 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 $[eqn]$ symmetry at the subsystem level.

Second, we determine whether the QME exists 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 ( $[eqn]$ ) resides in the half-filling sector ( $[eqn]$ ), which has the maximal Hilbert space dimension and thus the fastest thermalization rate. Increasing the tilt $[eqn]$ 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 $[eqn]$ 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 Equation (23), the steady-state plateau $[eqn]$ scales with $[eqn]$ . Since $[eqn]$ is a monotonically decreasing function of $[eqn]$ , a larger initial tilt $[eqn]$ (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 $[eqn]$ 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 $[eqn]$ . In the strong Stark limit ( $[eqn]$ ), restoring symmetry within a subsystem of size $[eqn]$ 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 $[eqn]$ sites requires virtually occupying intermediate states up the Stark ladder. According to high-order perturbation theory, the effective matrix element for a process spanning $[eqn]$ sites scales as $[eqn]$ . The timescale for symmetry restoration, and consequently the timescale on which the QME occurs, is inversely proportional to this rate:

$$[eqn]$$

This prediction implies that the Mpemba time grows exponentially with the subsystem size $[eqn]$ 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 3b and Figure 4d.

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 Equation (1), enters the Stark-MBL phase. To obtain reliable information on the energy spectrum and eigenstates, we numerically solve the static and dynamical properties of the system using the Exact Diagonalization method. Figure 2a shows the half-chain von Neumann entropy density $[eqn]$ as a function of the Stark field F for several system sizes N. In the small-F regime, $[eqn]$ remains finite and exhibits only a mild size dependence, consistent with the volume-law behavior of the ergodic phase. As F increases, $[eqn]$ decreases substantially and drifts downward with increasing N, indicating a crossover toward the area-law entanglement characteristic of Stark-MBL. The inset of Figure 2a demonstrates a finite-size scaling collapse as a function of $[eqn]$ , yielding a critical field $[eqn]$ (indicated by the vertical green dotted line). Based on this result, we use $[eqn]$ as a representative Stark-MBL point in our subsequent analysis.

To further obtain dynamical evidence of non-ergodicity, we examine the charge imbalance $[eqn]$ starting from the initial TNS. Figure 2b plots the dynamics of $[eqn]$ at $[eqn]$ and $[eqn]$ for several tilt angles $[eqn]$ . After an early-time transient, $[eqn]$ rapidly approaches a non-zero long-time plateau. Crucially, the plateau value remains strongly dependent on the initial angle $[eqn]$ , preserving the ordering of the initial conditions over extremely long time scales (up to $[eqn]$ ). 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 $[eqn]$ 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 $[eqn]$ , where $[eqn]$ is the projection operator onto the eigenspace of the local charge operator $[eqn]$ with eigenvalue $[eqn]$ . To further characterize the spatial memory of the initial configurations, we calculate the charge imbalance (CI), which for our spin- $[eqn]$ system is defined as

$$[eqn]$$

Figure 3a,b show the EA dynamics for the TNS at $[eqn]$ and $[eqn]$ , respectively. At the ergodic point $[eqn]$ [Figure 3a], while larger $[eqn]$ 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 $[eqn]$ for the TNS is nearly insensitive to $[eqn]$ . However, deep in the Stark-MBL regime ( $[eqn]$ ), the EA curves exhibit clear signals, showing the emergence of the QME. A broader scan in the $[eqn]$ plane for the TNS is provided in Appendix A; see Figure A1, which further supports the absence of QME at small F and its emergence beyond a characteristic field scale. In contrast, for the TFS at $[eqn]$ [Figure 4a], the QME is already observable, which is highly consistent with the microscopic feature shown in its inset where $[eqn]$ significantly shifts toward the central sectors as $[eqn]$ 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 3c (TNS) and Figure 4c (TFS), we plot the steady-state EA $[eqn]$ versus $[eqn]$ for $[eqn]$ . 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 $[eqn]$ . This provides strong evidence that our analytical framework quantitatively reproduces the long-time dephasing behavior of Stark-MBL systems under symmetry constraints.

(a,b) EA dynamics with N=14 and NA=3. The initial state is chosen as TNS. The numerical curves have been smoothed by convolution with a Gaussian, ω(n)=e−(n/σ)2, with σ=5. The shaded curves represent the original numerical data before smoothing. Inset in (a,b): Charge-sector weights pA(qA) with N=14 and NA=4, for the initial TNS. (c) Steady-state EA for F=15>Fc and different subsystem sizes NA, with the initial state prepared as the TNS. Discrete points show the steady-state EA averaged over 2000 samples, while solid lines denote the analytical steady-state values. Inset in (c): analytical initial EA at t=0. (d) QME time tM as a function of subsystem size NA for different F at N=14, starting from the initial TNS. Discrete symbols (with error bars) show the numerically extracted tM, while dashed lines denote the theoretical scaling curves obtained by fitting Equation (24) to the numerical data.

(a,b) EA dynamics with N=14 and NA=3. The initial state is chosen as TFS. The numerical curves have been smoothed by convolution with a Gaussian, ω(n)=e−(n/σ)2, with σ=5. The shaded curves represent the original numerical data before smoothing. Inset in (a,b): Charge-sector weights pA(qA) with N=14 and NA=4, for the initial TFS. (c) Steady-state EA versus θ for F=15>Fc and different subsystem sizes NA, with the initial state prepared as the TFS. Discrete points show the steady-state EA averaged over late-time 2000 samples, while solid lines denote the analytical steady-state values. Inset in (c): analytical initial EA at t = 0. (d) QME time tM as a function of subsystem size NA for different F at N = 14, starting from the initial TFS. Discrete symbols (with error bars) show the numerically extracted tM, while dashed lines denote the theoretical scaling curves obtained by fitting Equation (24) to the numerical data.

3.2.3. Scaling Characteristics of the Mpemba Time

Finally, we quantify the characteristic timescale of the QME. We define the Mpemba time $[eqn]$ (corresponding to the first intersection of the EA dynamical envelopes) by extracting the first point of the EA curves and fitting the numerical results using the analytically derived scaling relation Equation (24). As shown in Figure 3d and Figure 4d, the fitting curves (dashed lines) perfectly capture the exponential growth of $[eqn]$ with respect to both subsystem size $[eqn]$ 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 $[eqn]$ decays as $[eqn]$ .

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 found 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 depends on the specific choice of the initial state (e.g., present for TFS but absent for TNS), the Stark-MBL regime imposes 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 $[eqn]$ , is governed by high-order off-resonant processes. We derived a scaling law showing that $[eqn]$ 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, which confirm that symmetry restoration in Stark-MBL is driven by slow, perturbative quantum fluctuations rather than by 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 [70,71,72,73,74,75,76] in optical lattices or trapped ions [77], 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 [78,79,80,81,82,83,84,85,86,87,88,89,90,91,92] would be crucial for understanding its fate in realistic open environments. Additionally, extending the analysis to the framework of nonstabilizerness [55,93,94,95,96,97,98,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] 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.