!!!***** SAVE THIS FILE REGULARLY. CHANGES ARE NOT SAVED AUTOMATICALLY ****!!! 1. ABOUT THE DATASET -------------------- Title: Dataset for "Tunable quantum Mpemba effect in long-range interacting systems" Creator(s): Andrew Hallam[1], Matthew Yusuf[1], Aashish A. Clerk[2], Ivar Martin[2,3], Zlatko Papic[1] Organisation(s): 1. University of Leeds. 2. University of Chicago. 3. Argonne National Laboratory Rights-holder(s): Unless otherwise stated, Copyright 2026 University of Leeds Publication Year: 2026 Description: This dataset contains the numerical simulation data underlying all main-text and Supplemental Material figures of the paper "Tunable quantum Mpemba effect in long-range interacting systems". The data comprise quantum many-body dynamics simulations of a one-dimensional long-range interacting spin-1/2 XYZ/XXZ spin chain, obtained using the time-dependent variational principle (TDVP) for both finite and translation-invariant infinite matrix product states (MPS/iMPS), together with ground-state density matrix renormalization group (DMRG) calculations. The data include entanglement entropies, entanglement asymmetry, trace distance, order-parameter dynamics, and Mpemba-time/critical-interaction-range data used to characterize the quantum Mpemba effect and its interplay with spontaneous symmetry breaking as a function of interaction range, magnetic field strength, and initial-state energy density. Full details of the models, methods, and definitions are given in the accompanying publication and its Supplemental Material. Cite as: Hallam, Andrew; Yusuf, Matthew; Clerk, Aashish A.; Martin, Ivar; Papic, Zlatko (2026): Dataset for "Tunable quantum Mpemba effect in long-range interacting systems". University of Leeds. [Dataset] https://doi.org/10.5518/1825 Related publication: A. Hallam, M. Yusuf, A. A. Clerk, I. Martin, and Z. Papic, "Tunable quantum Mpemba effect in long-range interacting systems", accepteed to Physical Review B (2026). [Update this entry with the final journal reference, volume, page/ article number, and DOI once the article is accepted and published.] Contact: py20mzy@leeds.ac.uk 2. TERMS OF USE --------------- Copyright 2026 University of Leeds. Unless otherwise stated, this dataset is licensed under a Creative Commons Attribution 4.0 International Licence: https://creativecommons.org/licenses/by/4.0/. 3. PROJECT AND FUNDING INFORMATION ---------------------------------- Title: Tunable quantum Mpemba effect in long-range interacting systems Dates: [project start-end - please fill in] Funding organisation(s) and grant(s): - Leverhulme Trust Research Leadership Award RL-2019-015 - EPSRC Grant EP/Z533634/1 - EPSRC Grant EP/W524372/1 - UKRI1337 - NSF PHY-2309135 (supporting the Kavli Institute for Theoretical Physics, KITP) - Erwin Schrodinger International Institute for Mathematics and Physics (support for Z. Papic) - US Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division (support for I. Martin) 4. CONTENTS ----------- File listing FIGURE 2 DATA FILE (main text) ================================ Dataset: XYZ_mpemba_tc_hz.mat Figure(s) generated from this file: Figure 2(a)-(f) of the main text. Description: Single MATLAB data file containing all quench-dynamics data (expectation values of the prethermal Hamiltonian, and entanglement/asymmetry entropies) used to produce Fig. 2 of the main text: quench dynamics under the long-range XYZ Hamiltonian [Eq. (1)] for the tilted product initial states [Eq. (9)], at interaction ranges alpha = 1.5 and alpha = 4, and tilt angles theta = pi/8 and pi/4, as a function of magnetic field strength hz and time t. Software used: - Data were generated using the time-dependent variational principle (TDVP) for translation-invariant infinite matrix product states (iMPS), as described in Sec. III of the main text and Sec. II of the Supplemental Material. - File format: MATLAB v5.0 .mat file. - Reading the file: MATLAB (native), or Python via scipy.io.loadmat. No proprietary toolboxes are required. Naming convention: All variables follow the pattern: _alpha__theta_ - in {D_out, SE, SE_A} - see data dictionary below. - in {15, 4} - interaction range alpha = 1.5 or alpha = 4 (Kac-normalized power-law exponent, dimensionless). - in {0125, 025} - tilt angle of the initial product state [Eq. (9)], given as a fraction of pi: 0125 -> theta = 0.125*pi = pi/8 025 -> theta = 0.25*pi = pi/4 This gives 4 parameter combinations (alpha, theta) x 3 quantities = 12 variables in total. Data dictionary: Number of variables: 12 Array axes: each variable is indexed as (field-strength index, time index[, subsystem index]). Variables: D_out_alpha__theta_ Shape: (13, 401), complex128 Description: D(t) = , expectation value of the U(1)-symmetric prethermal effective Hamiltonian [Eq. (3) of the main text], evaluated along the quench. Used (after normalizing by D(t=0)) for Fig. 2(a)-(b). Values are numerically real; residual imaginary parts (~1e-17) are floating-point/TDVP artifacts and can be discarded (take the real part). Units: dimensionless (units of J) Variables: SE_alpha__theta_ Shape: (13, 401), float64 Description: S(rho_A)(t), von Neumann entanglement entropy of the (non-charge-projected) reduced density matrix for a subsystem of NA = 4 spins. Combined with SE_A_alpha__theta_[:,:,3] this gives the entanglement asymmetry Delta S_A(t) = S(rho_A,Q) - S(rho_A) used in Fig. 2(c)-(f). Units: dimensionless (nats) Variables: SE_A_alpha__theta_ Shape: (13, 401, 4), complex128 Description: S(rho_A,Q)(t), von Neumann entropy of the charge-resolved (symmetrized) reduced density matrix [Eq. (6) of the main text], resolved for subsystem sizes NA = 1, 2, 3, 4 along the third axis (index 0 -> NA=1, ..., index 3 -> NA=4). The NA=4 slice, [:,:,3], is used together with SE_alpha__theta_ to construct Delta S_A(t) in Fig. 2(c)-(f). Residual imaginary parts are numerical artifacts and can be discarded. Units: dimensionless (nats) Axis conventions: Time axis (401 points, second index): t in [0, 10] in steps of dt = 0.025. Field-strength axis (13 points, first index): hz in {0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0}, linearly spaced from 0 to 3. Missing data codes: None. All entries are populated numerical values; there are no NaN or sentinel codes in this dataset. Specialized formats / abbreviations: - alpha: power-law interaction decay exponent, Kac-normalized per Eq. (2) of the main text. - theta: tilt angle of the initial product state, Eq. (9) of the main text. - D_hat: prethermal effective Hamiltonian operator, Eq. (3) of the main text. - rho_A, rho_A,Q: reduced and charge-resolved reduced density matrices over subsystem A, Eq. (6) of the main text. - Complex dtype on real-valued quantities (D_out, SE_A) reflects intermediate complex arithmetic in the TDVP/iMPS code; imaginary parts are at the level of numerical noise (~1e-15 to 1e-17) and are not physically meaningful. FIGURE 3 DATA FILES (main text) ================================== Figure(s) generated from these files: Figure 3(a)-(c) of the main text. File listing: SE_A_data_inf.mat Entanglement asymmetry dynamics for the tilted product states [Eq. (9), theta = pi/8, pi/4] at Jz = -0.75, used for Fig. 3(a). Also contains the corresponding dynamics for the tilted Neel states [Eq. (11), "2-site" variables] at the same Jz = -0.75, from which the Jz = -0.75 point of the 2-site tau_M(alpha) curve (stored in mpemba_time_states.mat, row 2) was extracted. mpemba_time_states.mat Mpemba time tau_M(alpha) at Jz = -0.75, for both the 1-site tilted product states [Eq. (9)] and 2-site tilted Neel states [Eq. (11)]. Used for Fig. 3(b), and supplies the Jz = -0.75 point of the 2-site critical line in Fig. 3(c). mpemba_2site_data.mat Entanglement asymmetry dynamics for the tilted Neel states [Eq. (11), theta = pi/8, pi/4] at Jz = 0, -0.25, -0.5, as a function of alpha. Supplies the remaining three points (Jz = 0, -0.25, -0.5) of the 2-site critical line, alpha_M^2-site (pink dotted line), in Fig. 3(c). mpemba_diagram_data_new.mat (Jz, alpha) grid, corresponding Mpemba time tau_M, and the extracted critical lines alpha_M (first crossing) and alpha_M^(2) (last crossing) for the 1-site tilted product states. Used for Fig. 3(c). mpemba_response_data.mat Order-parameter dynamics ||, together with the underlying entanglement-asymmetry data (resolved by subsystem size NA = 1-4), over the full (alpha, Jz) grid, for theta = pi/4 and pi/8. Used to extract the critical line alpha_sigma+(Jz) in Fig. 3(c). Also used, together with SE_A_data_inf.mat, for SM Figs. 5 and 8 (see below). alpha_critical_192_200.csv Ground-state SSB boundary alpha_c(Jz) (central-charge method, DMRG N = 192, 200), plotted as the solid line in Fig. 3(c). alpha_critical_ES_192_200_1.csv Alternative ground-state SSB boundary estimate alpha_c^ES(Jz) (entanglement-spectrum method), defining the shaded uncertainty band together with alpha_c in Fig. 3(c). Combining the 2-site critical line: the pink dotted line alpha_M^2-site(Jz) in Fig. 3(c), computed over Jz in [-0.75, 0], is assembled from TWO files: the Jz = -0.75 point comes from mpemba_time_states.mat (row 2), and the Jz = -0.5, -0.25, 0 points are computed from mpemba_2site_data.mat. Software used: Same as for Fig. 2: TDVP for translation-invariant iMPS (dynamics), and finite DMRG with the ITensor library (ground-state alpha_c, alpha_c^ES boundaries), as described in Sec. IV of the main text and Secs. II-III of the Supplemental Material. The .mat files are MATLAB v5.0 format; the .csv files are plain comma-separated text. Naming convention: - theta0125 / theta025 Tilt angle theta (or phi, for Neel states) = 0.125*pi = pi/8, or 0.25*pi = pi/4. - 1site / 2site Initial-state family: 1-site = tilted ferromagnetic product state [Eq. (9)]; 2-site = tilted Neel state [Eq. (11)]. - Jz0, Jzm025, Jzm05 Fixed Jz value ("m" = minus): Jz = 0, -0.25, -0.5 respectively. - SE_A__theta (in mpemba_response_data.mat) Charge-resolved entropy S(rho_A,Q) for subsystem size NA = n (n = 1-4), at tilt angle theta. - SE_theta Von Neumann entropy S(rho_A) of the actual (non-projected) reduced density matrix, NA = 4. - Sp_theta Magnitude of the order parameter, ||. Data dictionary: File: mpemba_time_states.mat Variable: mpemba_time_states_total Shape: (3, 30), float64 Description: Row 0: alpha grid, 30 values, 1.1-4.0 in steps of 0.1. Row 1: tau_M(alpha) for the 1-site tilted product states. Row 2: tau_M(alpha) for the 2-site tilted Neel states. Both rows at Jz = -0.75. Units: alpha - dimensionless; tau_M - units of 1/J Missing data code: 0.0 = no finite Mpemba crossing observed (tau_M -> infinity; QME does not occur). ---- File: SE_A_data_inf.mat Variables: SE_A_theta0125_1site, SE_A_theta025_1site Shape: (30, 601), complex128 Description: S(rho_A,Q)(t) for the 1-site tilted product state, theta = pi/8 or pi/4, at Jz = -0.75, over the 30-point alpha grid. Residual imaginary parts are numerical artifacts. Units: dimensionless (nats) Variables: SE_A_theta0125_2site, SE_A_theta025_2site Shape: (30, 601), complex128 Description: Same, for the 2-site tilted Neel state (phi = pi/8, pi/4). Units: dimensionless (nats) Time axis: 601 points, dt = 0.025, giving t in [0, 15]. ---- File: mpemba_2site_data.mat Variables: SE_A_Jz_theta (Z = 0, m025, m05; T = 0125, 025) Shape: (30, 1001), complex128 Description: S(rho_A,Q)(t), 2-site tilted Neel state, tilt phi = pi/8 or pi/4, at fixed Jz in {0, -0.25, -0.5}, over the 30-point alpha grid. Units: dimensionless (nats) Time axis: 1001 points, dt = 0.025, giving t in [0, 25]. ---- File: mpemba_diagram_data_new.mat Variable: Jz Shape: (1, 19), float64 Description: Jz grid, -0.9 to 0.9 in steps of 0.1. Units: dimensionless Variable: alpha_1 Shape: (1, 30), float64 Description: alpha grid, 1.1 to 4.0 in steps of 0.1 (same grid as elsewhere). Units: dimensionless Variable: T_mpemba Shape: (30, 19), float64 Description: tau_M(alpha, Jz) for the 1-site tilted product states, indexed as (alpha_1, Jz). Fig. 3(c) plots the color map of tau_M^-1, i.e. the reciprocal of this array. Units: units of 1/J Missing data code: 0.0 = no finite Mpemba crossing observed within the simulated time window (region below alpha_M, where tau_M -> infinity). Variable: alpha_m Shape: (1, 19), float64 Description: Critical line alpha_M(Jz): largest alpha for which tau_M -> infinity (first crossing), the red solid line in Fig. 3(c). Units: dimensionless Variable: alpha_m2 Shape: (1, 19), float64 Description: Critical line alpha_M^(2)(Jz): analogous boundary defined using the LAST crossing of the entanglement-asymmetry ratio through unity, instead of the first. Units: dimensionless Note: this Jz grid (step 0.1, -0.9 to 0.9) does NOT include Jz = -0.75. The special value Jz = -0.75 used elsewhere in the paper (Figs. 2, 3(a), and several SM figures) is stored separately, at finer/dedicated resolution, in SE_A_data_inf.mat. ---- File: mpemba_response_data.mat Variables: SE_A__theta (n = 1, 2, 3, 4; T = 0125, 025) Shape: (30, 19, 1001), complex128 Description: S(rho_A,Q)(t) for the 1-site tilted product state, tilt theta = pi/8 or pi/4, resolved for subsystem size NA = n, over the full (alpha, Jz) grid (alpha_1 x Jz, from mpemba_diagram_data_new.mat) and 1001 time points. Units: dimensionless (nats) Variables: SE_theta0125, SE_theta025 Shape: (30, 19, 1001), float64 Description: S(rho_A)(t), actual (non-projected) von Neumann entropy, NA = 4, same (alpha, Jz, t) grid. Units: dimensionless (nats) Variables: Sp_theta0125, Sp_theta025 Shape: (30, 19, 1001), float64 Description: ||, magnitude of the U(1) order parameter [Fig. 1(b)], same grid. Verified: at t = 0, Sp_theta0125[alpha, Jz, 0] = 0.5*sin(2*pi/8) = 0.3536, and Sp_theta025[alpha, Jz, 0] = 0.5, matching _(t=0) = 0.5*sin(2*theta). Units: dimensionless Axes: first index = alpha (30 pts, 1.1-4.0), second index = Jz (19 pts, -0.9 to 0.9), third index = time (1001 pts, dt = 0.025, t in [0, 25]). alpha_sigma+(Jz) is derived from Sp_theta025 / Sp_theta0125 by locating where the ratio |^theta=pi/4 / ^theta=pi/8| remains above 1, serving as the upper-bound estimate for alpha_M shown in Fig. 3(c). ---- File: alpha_critical_192_200.csv Row 1: Jz, 16 points, -0.875 to 1.0 in steps of 0.125. Units: dimensionless Row 2: 1/alpha_c(Jz) - this row is the RECIPROCAL of alpha_c (values 0.331-0.368; invert to recover alpha_c approx 2.7-3.0, consistent with alpha_c approx 3 quoted in the main text). Units: dimensionless ---- File: alpha_critical_ES_192_200_1.csv Row 1: Jz (values -0.936 to 2.608, irregularly spaced). Units: dimensionless Row 2: alpha_c^ES(Jz) - given directly, no reciprocal needed (values 2.0-2.94). Units: dimensionless Missing data codes: mpemba_time_states.mat and mpemba_diagram_data_new.mat (T_mpemba): 0.0 = no finite Mpemba crossing detected (tau_M -> infinity). All other files: none. Specialized formats / abbreviations: - alpha_c, alpha_c^ES: ground-state SSB phase boundary (central-charge / entanglement-spectrum methods). - alpha_M, alpha_M^(2), alpha_M^2-site, alpha_sigma+: dynamical critical lines bounding the QME region in Fig. 3(c), obtained respectively from the first crossing, last crossing, 2-site states, and order-parameter ratio. - Complex dtype on real-valued entropy quantities reflects intermediate complex arithmetic in the TDVP/iMPS code; imaginary parts are numerical noise, not physically meaningful. SUPPLEMENTAL MATERIAL FIGURES (SM Figs. 1-8) ============================================== The datasets below correspond to figures in the Supplemental Material (SM), as distinct from Figs. 2 and 3 of the main text documented above. Several SM figures reuse variables already documented in the main-text Fig. 3 entry (SE_A_data_inf.mat, mpemba_response_data.mat, mpemba_diagram_data_new.mat) - in those cases the variable-level dictionary is not repeated below; only the additional file (if any) and the SM-specific combination/derivation are described. SM FIGURE 1 DATA FILES (trace distance, SM Sec. I) ----------------------------------------------------- Figure generated from this file: SM Fig. 1. File: tracedist.zip Contents: 51 CSV files, one per (alpha, theta, hz) combination. Naming convention (folder tracedist/, per-file): tracedist_tracedist_new, N=, alpha=, theta=, hz=, subsize=, Jx=, Jy=, Jz=.csv - N = total system size (100 throughout). - ALPHA in {1.5, 4.0} - interaction range alpha. - TDENOM in {4.0, 8.0} - tilt angle given as the denominator of pi, i.e. theta = pi/TDENOM (TDENOM=4.0 -> theta=pi/4; TDENOM=8.0 -> theta=pi/8). Note: this is a DIFFERENT encoding convention from the theta0125/theta025 naming used in the .mat files elsewhere in this dataset - here the number IS the denominator, not a fraction of pi. - HZ = magnetic field strength, nominally 0.0 to 3.0 in steps of 0.25 (13 values). The (alpha=1.5, theta=4.0, hz=2.75) combination is intentionally absent from the dataset - this run was not performed. - NA (subsize) = 4 throughout. - Jx=-0.5, Jy=-1.5, Jz=-0.75 throughout (matches main-text Fig. 2 XYZ couplings). File contents: each CSV has 2 columns and 200 rows, no header row. Column 1: time t. Values run 0.0 to ~9.95 in steps of dt=0.05. Column 2: trace distance Delta_theta(rho_A)(t) [Eq. (1) of the SM], for the single theta value encoded in the filename (NOT a ratio - the ratio plotted in SM Fig. 1 is obtained by dividing the theta=pi/4 file by the theta=pi/8 file, at matching alpha and hz). Units: t - units of 1/J; trace distance - dimensionless (bounded in [0,1]). Software: generated using TDVP for finite MPS, implemented using the ITensor library, max bond dimension chi=100, as described in the caption of SM Fig. 1 and Sec. I. The reference thermal density matrix rho_A^th was obtained via exact diagonalization on a smaller (N=14 spin) system. Missing data codes: none. SM FIGURE 2 DATA FILE (K-term convergence check, SM Sec. II) ---------------------------------------------------------------- Figure generated from this file: SM Fig. 2 (added in response to referee comment on the K=8 exponential approximation to the power-law interaction). File: mpemba_exp.mat Variables: SE_A_6exp_theta0125, SE_A_6exp_theta025, SE_A_8exp_theta0125, SE_A_8exp_theta025, SE_A_10exp_theta0125, SE_A_10exp_theta025 Shape: (15, 401), complex128 Description: entanglement asymmetry Delta S_A(alpha, t) - already the final asymmetry, not the raw charge-resolved entropy S(rho_A,Q) (see naming note below) - for the tilted product states [Eq. (9) of the main text], theta = pi/8 or pi/4, evolved under the XXZ Hamiltonian with Jx=Jy=-1, Jz=-0.75, where the power-law interaction was fitted using K = 6, 8, or 10 exponential terms [Eq. (3) of the SM]. First axis: 15 representative values of alpha, linearly spaced from 1.2 to 4.0 in steps of 0.2. Second axis: 401 time points, dt=0.025 assumed, t in [0,10]. Units: dimensionless (nats) Naming convention: exp = number of exponential terms K used to fit the power-law interaction (K in {6,8,10}); theta0125 / theta025 = tilt angle pi/8 or pi/4. Note on naming: unlike the SE_A_alpha__theta_ variables in XYZ_mpemba_tc_hz.mat (Fig. 2 main text) and the SE_A_theta variables in SE_A_data_inf.mat (Fig. 3 main text), which store the raw charge-resolved entropy S(rho_A,Q) and require pairing with a separate S(rho_A) variable to construct Delta S_A, the SE_A_exp_theta variables in this file already ARE the final entanglement asymmetry Delta S_A(t). Used to construct: (a) the ratio Delta S_A^(theta=pi/4)/ Delta S_A^(theta=pi/8) for K=6 (dotted), K=8 (dashed), and K=10 (solid); (b) the relative error of the K=8 approximation relative to K=10, (Delta S_A^K=10 - Delta S_A^K=8)/ Delta S_A^K=10, as a function of time. Missing data codes: none. SM FIGURE 3 DATA FILES (finite-size scaling, SM Sec. III) --------------------------------------------------------------- Figure generated from these files: SM Fig. 3. File: sscaling.zip Contents: 60 CSV files, one per (N, alpha, theta) combination. Naming convention (per-file, no subfolder): SA_sscaling, N=, alpha=, theta=, hz=0, subsize=4, Jx=-1, Jy=-1, Jz=-0.75.csv - N in {50, 100, 200} - finite system size. - ALPHA in {1.1, 1.4, 1.7, 2.0, 2.3, 2.6, 2.9, 3.2, 3.5, 3.8} - 10 values, step 0.3. - TDENOM in {4.0, 8.0} - tilt angle pi/TDENOM (same convention as tracedist.zip above: TDENOM is the denominator of pi). - hz=0 (placeholder field-strength parameter, not physically relevant here since this is the effective XXZ Hamiltonian [Eq. (3) of the main text], not the full field-driven model). - Jx=Jy=-1, Jz=-0.75 throughout. File contents: each CSV has 2 columns and 200 rows, no header row. Column 1: time t, 0.0 to ~9.95 in steps of dt=0.05. Column 2: entanglement asymmetry Delta S_A^theta(t) for the finite system size N encoded in the filename. Units: t - units of 1/J; Delta S_A - dimensionless (nats). Software: TDVP for finite MPS, implemented using the ITensor library, max bond dimension chi=100, time step dt=0.05, as stated in SM Sec. III. Combining with the N -> infinity curve: SM Fig. 3 also shows the N -> infinity (iMPS) limit for comparison. This curve is the SAME data already documented under SE_A_theta0125_1site and SE_A_theta025_1site in SE_A_data_inf.mat (see the Fig. 3 main-text entry above) - no separate file is needed for the N -> infinity curve. Missing data codes: none. SM FIGURE 4 DATA FILE (SU(2) point, SM Sec. IV) ---------------------------------------------------- Figure generated from this file: SM Fig. 4. File: SE_A_SU2.mat Variable: SE_ratio_SU2 Shape: (30, 1001), float64 Description: ratio of the SU(2) entanglement asymmetry, Delta S_A^(1)/Delta S_A^(2) [using the SU(2) charge operator, Eq. (4) of the SM], for the tilted Neel states [Eq. (11) of the main text] with phi1=pi/4 and phi2=pi/8, evolved under the long-range Heisenberg Hamiltonian (Jx=Jy=Jz=-1). Already the final ratio - no further combination needed. First axis: alpha, 30 values (1.1 to 4.0 in steps of 0.1, same grid as elsewhere - not stored explicitly in this file but consistent with the convention used throughout). Units: dimensionless Variable: T1 Shape: (1, 1001), float64 Description: explicit time axis for SE_ratio_SU2. Values run 0.0 to 25.0 in steps of dt=0.025 (verified directly from the file). Units: 1/J Missing data codes: none. SM FIGURE 5 DATA (subsystem-size dependence of alpha_M, SM Sec. III) -------------------------------------------------------------------------- Figure generated from: SM Fig. 5. This figure does not require a new file. It is derived entirely from variables already documented under mpemba_response_data.mat in the main-text Fig. 3 entry above: SE_A_1_theta0125, SE_A_1_theta025 (NA = 1) SE_A_2_theta0125, SE_A_2_theta025 (NA = 2) SE_A_3_theta0125, SE_A_3_theta025 (NA = 3) SE_A_4_theta0125, SE_A_4_theta025 (NA = 4) each of shape (30, 19, 1001), complex128, indexed by (alpha, Jz, t) as described previously. Derivation: for each subsystem size NA = 1, 2, 3, 4 separately, the entanglement asymmetry ratio Delta S_A^(theta=pi/4)/ Delta S_A^(theta=pi/8) is constructed from the corresponding SE_A__theta slice (combined with SE_theta, also in the same file, for the S(rho_A) subtraction), and the critical value alpha_M(Jz) at which the Mpemba time diverges is extracted exactly as for the NA=4 case in main-text Fig. 3(c). SM Fig. 5 plots alpha_M(Jz) for each of the four subsystem sizes NA = 1-4 on the same axes. Missing data codes: as noted in the Fig. 3 main-text entry, no additional codes beyond what is already documented for this file. SM FIGURE 6 DATA (other initial-state combinations, SM Sec. V) --------------------------------------------------------------------- Figure generated from: SM Fig. 6. This figure does not require a new file. It is derived entirely from the four variables already documented in SE_A_data_inf.mat under the main-text Fig. 3 entry above: SE_A_theta0125_1site, SE_A_theta025_1site (1-site / theta family) SE_A_theta0125_2site, SE_A_theta025_2site (2-site / phi family) Panel-by-panel combination (all at Jz=-0.75, Jx=Jy=-1): (a) phi1=pi/4 vs phi2=pi/8: SE_A_theta025_2site / SE_A_theta0125_2site (b) theta1=pi/4 vs phi2=pi/4: SE_A_theta025_1site / SE_A_theta025_2site (c) theta1=pi/4 vs phi2=pi/8: SE_A_theta025_1site / SE_A_theta0125_2site (d) theta1=pi/8 vs phi2=pi/8: SE_A_theta0125_1site / SE_A_theta0125_2site (e) phi1=pi/4 vs theta2=pi/8: SE_A_theta025_2site / SE_A_theta0125_1site Missing data codes: none (as documented for this file previously). SM FIGURE 7 DATA FILES (Haar-random initial states, SM Sec. VI) ------------------------------------------------------------------------- Figures generated from these files: SM Fig. 7(a)-(d). File: chi2_rand20_data.mat [SM Fig. 7(a)] Variable: SE_A Shape: (25000, 4), complex128 Description: entanglement asymmetry Delta S_A at t=0 (already the final quantity, same naming convention as mpemba_exp.mat above), for the full ensemble of 25,000 Haar-random MPS states (physical dimension d=2, bond dimension chi=2), resolved for subsystem sizes NA=1,2,3,4 along the second axis. SM Fig. 7(a) plots one NA slice (NA=4, i.e. column index 3) against the corresponding normalized energy density. Units: dimensionless (nats) Variable: eps_out Shape: (25000, 1), float64 Description: normalized energy density epsilon [Eq. (10) of the main text] for each of the 25,000 random states. Units: dimensionless Variables: dat__ (E = 1,2,3,4; G = 1,2) Shape: (20, 2), float64 each Description: (epsilon, Delta S_A) coordinate pairs for the 20 selected/highlighted states per group, used to overlay the red markers in SM Fig. 7(a). E indexes the target normalized energy density: E=1 -> epsilon~0.05, E=2 -> epsilon~0.1, E=3 -> epsilon~0.15, E=4 -> epsilon~0.2 (verified directly from the data ranges). G indexes the entanglement-asymmetry group: G=1 -> "high" Delta S_A, G=2 -> "low" Delta S_A. Units: dimensionless Missing data codes: none. ---- File: deltaSA_rand_chi180_1.mat [SM Fig. 7(b)-(d)] Variables: deltaSA_rand__ (E = 1,2,3; G = 1,2) Shape: (30, 321), complex128 Description: entanglement asymmetry Delta S_A(alpha, t), averaged over the 20 selected Haar-random MPS states, for fixed normalized energy density epsilon (E=1 -> ~0.05, E=2 -> ~0.1, E=3 -> ~0.15, matching the same E convention and value ranges verified in chi2_rand20_data.mat above) and group G (G=1 -> "high" Delta S_A group, G=2 -> "low" Delta S_A group). First axis: alpha, 30 values (1.1-4.0 grid assumed, consistent with the rest of the dataset). Second axis: time, 321 points, dt= 0.025 assumed, t in [0,8]. Units: dimensionless (nats) Note: there is no E=4 (epsilon~0.2) pair in this file - only the E=1,2,3 groups were evolved dynamically; the E=4 group appears only in the static scatter plot data (chi2_rand20_data.mat, SM Fig. 7(a)). Bond dimension: chi=180. Used to construct: SM Fig. 7(b): deltaSA_rand_1_1, plotted against 1/t. SM Fig. 7(c): deltaSA_rand_1_1 / deltaSA_rand_1_2 (epsilon~0.05 ratio). SM Fig. 7(d): deltaSA_rand_2_1 / deltaSA_rand_2_2 (epsilon~0.1 ratio). Missing data codes: none. SM FIGURE 8 DATA (infinite-time extrapolation, SM Sec. VII) ------------------------------------------------------------------- Figure generated from: SM Fig. 8(a)-(b). This figure does not require a new file beyond what is already documented above. It is derived from two files, covering two different parts of the figure: (i) SE_A_theta0125_1site and SE_A_theta025_1site in SE_A_data_inf.mat (main-text Fig. 3 entry above) - the raw Delta S_A(alpha, t) dynamics at the specific value Jz=-0.75, over the t in [0,15] window (601 time points, dt=0.025). This is the dedicated dataset for Jz=-0.75, which is used for SM Fig. 8(a) inset (convergence of alpha_M at Jz=-0.75 with increasing simulation time ts) and for SM Fig. 8(b) (the exponential fit and residual value c(alpha) at Jz=-0.75). (ii) SE_A_1_theta through SE_A_4_theta (together with SE_theta) in mpemba_response_data.mat (main-text Fig. 3 entry above) - the raw Delta S_A(alpha, Jz, t) dynamics over the full 19-point Jz grid (-0.9 to 0.9 in steps of 0.1) and t in [0,25] (1001 time points, dt=0.025). This is the dataset used for SM Fig. 8(a) main panel, which plots alpha_M(Jz) for several values of the maximum simulation time ts (e.g. ts=5,10,15), obtained by truncating this dynamics data at each ts and re-running the same crossing-detection procedure used for main-text Fig. 3(c). Note: the 19-point Jz grid in mpemba_response_data.mat / mpemba_diagram_data_new.mat does NOT include Jz=-0.75 (it steps in units of 0.1: ..., -0.8, -0.7, ...). This is why the Jz=-0.75 inset and the general-Jz main panel of SM Fig. 8(a) draw on two different underlying files, as described above. The dashed "alpha_M(ts -> infinity)" reference line in the inset is obtained from the infinite-time extrapolation of the exponential fit in SM Fig. 8(b) (Delta S_A^theta(t) = a*exp(-b*t) + c), not from mpemba_diagram_data_new.mat. Derivation: SM Fig. 8(a) [main panel]: alpha_M(ts), for ts = 5, 10, 15, across the Jz grid in (ii), by truncating the time series at each ts and locating the crossing of the entanglement-asymmetry ratio through unity, as for main-text Fig. 3(c). SM Fig. 8(a) [inset]: convergence of alpha_M(ts) at the specific value Jz=-0.75, using (i), compared to the infinite-time extrapolation value (dashed line) from the exponential fit described next. SM Fig. 8(b): the entanglement asymmetry curves in (i) are fit to Delta S_A^theta(t) = a*exp(-b*t) + c [a, b, c fit parameters], and the residual value c(alpha) is plotted for theta=pi/4 and pi/8; the crossing point of the two c(alpha) curves (alpha=2.2) marks the threshold for the onset of the QME in the t -> infinity limit at Jz=-0.75. Missing data codes: none beyond what is already documented for the constituent files. 5. METHODS ---------- Dynamics data (Figs. 2 and 3 of the main text, and SM Figs. 1-8) were generated using the time-dependent variational principle (TDVP) for translation-invariant infinite matrix product states (iMPS), applied to the long-range XYZ Hamiltonian [Eq. (1) of the main text] and its effective XXZ descendant [Eq. (3)], as well as TDVP for finite MPS at fixed system sizes (used for the trace distance and finite-size scaling checks in SM Figs. 1 and 3). Power-law interactions were approximated as a sum of K = 8 exponential terms (see Sec. II of the Supplemental Material for the fitting procedure and SM Fig. 2 for a K = 6, 8, 10 convergence check). Simulations used bond dimension chi = 128 (iMPS) or chi = 100 (finite MPS) and a time step of dt = 0.025 (iMPS) or dt = 0.05 (finite MPS), except where otherwise noted per file above. Quantities extracted from the resulting MPS/iMPS wavefunctions include: the expectation value of the prethermal effective Hamiltonian D(t); the von Neumann entanglement entropy S(rho_A) of the reduced density matrix over a subsystem of NA spins; the charge-resolved entropy S(rho_A,Q) obtained by projecting the reduced density matrix onto sectors of fixed U(1) (or SU(2), for the Heisenberg point, SM Fig. 4) charge [Eq. (6) of the main text; Eq. (4) of the SM]; the order parameter = (sigma_x + i*sigma_y)/2; and the trace distance Delta(rho_A)(t) = (1/2)||rho_A(t) - rho_A^th||_1 [Eq. (1) of the SM], with the reference thermal density matrix rho_A^th obtained by exact diagonalization on a smaller reference system and the effective temperature fixed self-consistently via Eq. (2) of the SM. Ground-state boundaries (alpha_c, alpha_c^ES) were obtained from finite-size DMRG using the ITensor library, with system sizes N = 192 and 200, at bond dimension chi = 250. alpha_c was determined from the effective central charge c_eff extracted from entanglement entropy scaling, defining the SSB boundary as the point where c_eff deviates from 1 by more than 10%. alpha_c^ES was obtained independently from the entanglement spectrum of the ground state, following the method of Ref. [108] of the main text. The finite-MPS TDVP simulations (trace distance, SM Fig. 1; finite- size scaling, SM Fig. 3) were likewise implemented using the ITensor library. Haar-random initial states (SM Fig. 7) were generated as an ensemble of 25,000 states with physical dimension d=2 and bond dimension chi=2, by drawing a Haar-random unitary U in U(d*chi), taking its first chi columns, and reshaping the resulting d*chi x chi matrix into a chi x d x chi local MPS tensor, repeated periodically to form a translation-invariant infinite state (SM Sec. VI). Twenty states were selected at each of four target normalized energy densities (epsilon = 0.05, 0.1, 0.15, 0.2) with relatively high or low entanglement asymmetry, and their subsequent dynamics evolved via iMPS-TDVP. The Mpemba time tau_M(alpha) (and the corresponding critical lines alpha_M, alpha_M^(2), alpha_M^2-site, alpha_sigma+, and their subsystem-size- and finite-time-resolved variants in SM Figs. 5 and 8) was extracted from the dynamics data by locating the crossing(s) of the entanglement-asymmetry ratio Delta S_A^theta1/Delta S_A^theta2 (or the analogous order-parameter ratio) through unity, as described in Sec. II B of the main text. Full methodological details, including simulation parameters, convergence checks, and finite-size scaling analysis, are given in the main text (Secs. II-IV) and the accompanying Supplemental Material (Secs. I-VII).