Far-from-equilibrium field theory of many-body quantum spin systems: Prethermalization and relaxation of spin spiral states in three dimensions

We study theoretically the far-from-equilibrium relaxation dynamics of spin spiral states in the three dimensional isotropic Heisenberg model. The investigated problem serves as an archetype for understanding quantum dynamics of isolated many-body systems in the vicinity of a spontaneously broken continuous symmetry. We present a field-theoretical formalism that systematically improves on mean-field for describing the real-time quantum dynamics of generic spin-1/2 systems. This is achieved by mapping spins to Majorana fermions followed by a 1/N expansion of the resulting two-particle irreducible (2PI) effective action. Our analysis reveals rich fluctuation-induced relaxation dynamics in the unitary evolution of spin spiral states. In particular, we find the sudden appearance of long-lived prethermalized plateaus with diverging lifetimes as the spiral winding is tuned toward the thermodynamically stable ferro- or antiferromagnetic phases. The emerging prethermalized states are characterized by different bosonic modes being thermally populated at different effective temperatures, and by a hierarchical relaxation process reminiscent of glassy systems. Spin-spin correlators found by solving the non-equilibrium Bethe-Salpeter equation provide further insight into the dynamic formation of correlations, the fate of unstable collective modes, and the emergence of fluctuation-dissipation relations. Our predictions can be verified experimentally using recent realizations of spin spiral states with ultracold atoms in a quantum gas microscope [S. Hild, et al. Phys. Rev. Lett. 113, 147205 (2014)].


I. INTRODUCTION
The equilibration of isolated quantum many-body systems is a fundamental and ubiquitous question in physics. It plays a central role in understanding a broad range of phenomena, including the dynamics of the early universe [1] , the evolution of neutron stars [2], pump-probe experiments in condensed matter systems [3], and the operation of semiconductor devices [4]. The simplest perspective on the problem is to recognize a dichotomy between ergodic and non-ergodic systems. The former exhibit fast relaxation to local equilibrium states occurring at microscopic timescales, followed by a slower relaxation process to global thermal equilibrium described by classical hydrodynamics of a few conserved quantities [5][6][7]. In contrast, non-ergodic systems possess an extensive set of conservation laws that prevent their thermalization [8,9].
In this work, we discuss the emergence of slow dynamics and prethermalization in translationally invariant spin systems that possess continuous symmetries. In higher dimensions, these systems exhibit thermodynamically stable symmetry broken phases with gapless Goldstone modes. For initial states tuned toward the symmetry broken phases, the slow Goldstone modes give rise to long-lived non-thermal states with a hierarchical relaxation dynamics that closely resembles aging in systems with quenched disorder. A non-perturbative treatment and beyond mean-field corrections are both found to be crucial for describing the relaxation process.
Specifically, we study the dynamics of the three dimensional (3D) isotropic Heisenberg model initially prepared in a spiral state, see Fig. 1 (a). The winding of the spin spiral sets the energy density, which is a conserved quantity in the unitary evolution. The spectrum of the Heisenberg model from ferromagnetic (FM) and antiferromagnetic (AFM) states can be traversed by tuning the spiral winding Q from 0 to π. In the vicinity of the FM and AFM states, we find that the system exhibits a slow hierarchical relaxation and can even come to a dynamical arrest, see Fig. 1 (b-c). Surprisingly, the relaxation dynamics is neither compatible with the trivial relaxation to local thermal equilibrium and slow hydrodynamic evolution, since the spin current is not conserved, nor with the linearized dynamics of the collective modes, which predicts exponentially growing out-of-plane instabilities.
The relaxation of the Néel spin spiral state with Q = π has been previously studied in the 1D Heisenberg model [30][31][32][33]. In contrast to the 1D case, the 3D Heisenberg model exhibits a symmetry broken thermal phase, which is central for the phenomenon discussed here. More recently, the dynamics of the Néel state in the Fermi-Hubbard model on an infinite dimensional Bethe lattice has been investigated [34], however, the approach to the steady state could not be studied due to the small effective exchange interaction.
From a technical perspective, our investigation of the non-equilibrium dynamics of spiral states has been enabled by developing a non-perturbative field theoretic formalism applicable to generic spin-1/2 systems for arbitrary initial states and geometries, which we refer to as the "Spin-2PI" formalism. This is achieved using a Majorana fermion representation of spin-1/2 operators [35,36], enlargement of the spin coordination number by a replica-symmetric extension, and ultimately a systematic 1/N fluctuation expansion of the real-time two-particle irreducible (2PI) effective action [37,38].
The recent rapid progress in the phenomenology of far-from-equilibrium quantum dynamics and its broad applications has been enabled by similar non-perturbative functional techniques.
Understanding the emergence of slow dynamics near thermodynamic phase transitions has implications reaching far beyond the domain of condensed matter physics. For instance, studies of non-equilibrium quantum fields in the context of inflation and early universe dynamics have suggested that the slowing down of quantum evolution near phase transitions is a plausible explanation for the large number of light particles and broken symmetries in the observable universe [53]. Given that the experimental verification of theories about early universe phenomena are typically rather indirect, experiments with synthetic many-body systems that allow precise monitoring of real-time dynamics close to phase transitions could play an important role in elucidating the emergence of slow evolution. The dynamics of various interacting spin systems have been already investigated in experiments with synthetic quantum matter, including domain formation in spinor condensates [54,55], the precise measurement of the evolution of spin flips in the ground state of 1D lattice spin systems [56][57][58][59], quantum coherences in long range models [60], and the relaxation dynamics of spin spiral states in 1D and 2D Heisenberg models [32,61]. The experimental observation of the dynamical phenomena discussed here are thus expected to be within close reach. This paper is organized as follows: In Sec. II, we introduce the Spin-2PI formalism, a technique we develop to study the dynamics of interacting spin systems. Complementary technical details are presented in App. A. We discuss the relaxation of spin spiral states in Sec. III. The phenomenon of dynamical slowing down and arrest will be presented Sec. III A, the long-time thermalization in Sec. III B, and the dynamic formation of correlations and instabilities in Sec. III C. We conclude our findings in Sec. IV.

II. THE SPIN-2PI FORMALISM
Consider a generic Hamiltonian describing the pairwise interaction between localized spin degrees of freedom on a given lattice L:Ĥ where V is an arbitrary interaction, j and k denote lattice sites, and {Ŝ α } are spin-1/2 operators. Summation over the repeated spin indices is assumed. HadŜ been The Spin-2PI formalism illustrated. The spins (green arrows) precess about a fluctuating exchange field (uncertain blue arrows). The quantum fluctuations of the exchange field are mediated by a real vector bosons (wiggly lines) and are suppressed by a factor of 1/N , permitting a systematic expansion.
classical angular momentum variables, the Hamiltonian dynamics of the system would be governed by the (nonlinear) Bloch equation: In case of quantum spins, the Bloch equation only describes the evolution of the spin expectation values Ŝ to the extent of which the mean-field Ansatz Ŝ jŜk ≈ Ŝ j Ŝ k is valid. The latter, however, is only justified for lattices with large coordination number, high spin particles, or in the presence of a high-temperature bath. The crucial role of quantum fluctuations in the dynamics of isolated spin-1/2 systems in finite dimensional lattices is beyond the reach of semi-classical methods, and demands a more careful treatment.
Here, we propose a formalism for transcending the mean-field approximation for spin evolution by a systematic inclusion of quantum corrections. This is achieved using functional methods and a variant of the large-N expansion technique. As a first step, we construct an auxiliary model in which each spin is replicated N times, and each bond is promoted to N 2 bonds between the replicas, with equal weight but with an overall scale factor of 1/N . The Hamiltonian of the auxiliary model is written as: The initial state |Ψ 0 is also subsequently promoted to an uncorrelated product in the replica space, The original problem is recovered by setting N = 1. We refer to the sum appearing in the parentheses in Eq. (3) as the exchange field operator,φ j , which plays the role of an effective fluctuating magnetic field with which the spins interact. The described large-N construction effectively increases the coordination number of each spin, z, to N z, thereby suppressing the fluctuations ofφ j according to the law of large numbers,φ j = ϕ c,j + O(1/ √ N z), where ϕ c,j ≡ ϕ j is the mean exchange field. In the limit of infinite N , the exchange field operator becomes effectively classical such that mean-field dynamics of the original modelĤ emerges as the asymptotically exact description of the dynamics in lim N →∞ĤN . For large but finite N , the fluctuations ofφ are small but not negligible, and can be systematically incorporated into the dynamics order by order in 1/N . This program can be carried out within the functional method of two-particle irreducible (2PI) effective actions. Crucially, truncating the expansion at a finite order in 1/N and taking the limit N → 1 yields non-perturbative and conserving approximations for the spin dynamics. We refer to this method as the Spin-2PI formalism, which is illustrated schematically in Fig. 2. In brief, spins precess about a self-consistently determined exchange mean field, and quantum spin fluctuations are mediated by the local and non-local exchange of a real vector boson whose propagator is suppressed by a factor of 1/N .
In the remainder of this section, we briefly outline the field theoretical developments that underlie the Spin-2PI formalism. Complementary technical details are given in App. A. A path integral for the spin-1/2 operators is constructed using a representation invoking Majorana fermions [35,36] The Majorana operators at each site {η 1 j , η 2 j , η 3 j } satisfy the Clifford algebra {η µ j , η ν k } = δ jk δ µν , from which the SU (2) algebra for spins [Ŝ α j ,Ŝ β k ] = iδ jk ε αβγŜ γ j and the Casimir condition S 2 j = 3/4 follow. The latter ensures a faithful spin-1/2 representation, thereby precluding unphysical states and obviating constraint gauge fields, in contrast to the usual representation using Schwinger charged slave particles [62].
Replacing the spin operators inĤ using Eq. (4), the Hamiltonian is mapped to that of a many-body system of Majorana fermions with quartic interactions. The large-N program can be identically followed by replicating the slave Majorana particles and assigning a replica index to each. We proceed by constructing a path integral for the Majorana fermions using fermionic coherent states on the closed time path (CTP) Schwinger-Keldysh contour.
The Lagrangian is given as: The exchange field is introduced by a Hubbard-Stratonovich decoupling of the quartic term using a real vector boson ϕ j on each lattice site. The non-equilibrium exchange mean field ϕ c , exchange field fluctuation propagator D, and the Majorana propagator G are introduced as: The integer variables are shorthand for the bundle of lattice site, contour time, spin and replica index. According to Eq. (4), the local spin expectation value is proportional to the fermion tadpole: We obtain the real-time evolution equations for G, D and ϕ c using the 2PI effective action formalism [37]. The effective action Γ[G, D, ϕ c ] is found by sourcing G, D and ϕ c and performing Legendre transformations: The bare Majorana and exchange propagators are given as G −1 0 (1, 2) = i∂ t1 δ(1, 2) and D −1 0 (1, 2) = N (V −1 ) δ(t 1 , t 2 ), respectively. M [ϕ c ] is the leading order (LO) self-energy (see Eq. A9). The evolution equations follow from making Γ stationary with respect to G, D, and ϕ c , see Eqs. (A7a)-(A7c).
Save for Γ 2 [G, D], the rest of the terms appearing in Γ[ϕ c , G, D] scale as O(N ) and together comprise the leading-order (LO) approximation. The next-to-leadingorder (NLO) corrections and beyond are represented by Γ 2 [G, D] which formally corresponds to the sum of 2PI vacuum diagrams arising from the cubic interaction vertex: The 1/N expansion of Γ int to the next-to-next-to-leading (NNLO) order is diagrammatically given as: (10) We have used the stationarity condition, Eq. (A7c), to omit ϕ c in favor of G in the LO interaction terms. The Feynman diagram rules are given in Sec. A 2.
Truncating the 1/N expansion of Γ at a finite order and setting N = 1 yields systematic improvements of the mean-field spin dynamics. The ensuing approximate theories are self-consistent and non-perturbative by construction, and respect the conservation laws associated to the global symmetries of the microscopic action, such as magnetization and energy. The latter is crucial for the long-time stability of the non-equilibrium dynamics.
We remark that in systems with large spin coordination number z, fluctuations of the exchange field are inherently suppressed and the expansion parameter is more accurately identified with 1/(zN ). Therefore, the large-N expansion of the Spin-2PI effective action in models with z > ∼ 1 is expected to be controlled and rapidly converging, even after taking the limit N → 1. Studies of the O(N ) model show that the most important correction to the mean-field (LO) approximation is captured by the NLO "fluctuation-exchange" diagram, along with negligible quantitative corrections from the subleading terms [63,64].
The Bloch equation is recovered upon truncating Γ at the LO level, see Sec. A 3. Truncations at NLO and beyond give rise to memory effects due to the dynamical fluctuations of the exchange field and result in a twotime Kadanoff-Baym integro-differential equation instead of the mean field Bloch equation, see Eqs. (A13a)-(A14b). Finally, higher order correlators, in particular the spinspin correlator iχ(1, 2) ≡ T C [Ŝ(1)Ŝ(2)] − Ŝ (1) Ŝ (2) , can be reconstructed with the knowledge of G and D by solving the non-equilibrium Bethe-Salpeter integral equation on the Schwinger-Keldysh contour, see App. A 3.

III. RELAXATION OF SPIN SPIRAL STATES IN THE 3D HEISENBERG MODEL
In this section, we investigate the unitary evolution of the spin spiral state on a 3D cubic lattice, under the isotropic Heisenberg HamiltonianĤ = −J ij Ŝ i ·Ŝ j using the Spin-2PI formalism developed in the previous section. Here, |→ j denotes the x-polarized state on lattice site j. The spiral is prepared in the xy-plane with a winding wavevector Q. We assume ferromagnetic couplings J > 0 for concreteness, even though the sign of the J does not affect the unitary evolution due to the time reversal symmetry of the Heisenberg model. The spiral state |sp(Q) is a simultaneous eigenstate ofŜ a (Q) ≡T aRz (Q a ), a = x, y, z, whereT a andR z (Q a ) denote the translation by one lattice site along the a-axis and rotation by angle Q a about the z-axis, respectively. The translation and rotation symmetries of the isotropic Heisenberg model imply [Ĥ,Ŝ a ] = 0, so that the spiral state |sp(Q) remains a simultaneous eigenstate of S a (Q) at all times in the course of unitary evolution. As a result, the out-of-plane magnetization Ŝ z j (t) vanishes identically, and the spiral magnetic order with the initial winding Q persists at all times. The transverse magnetization, is the only degree of freedom at the level of single spin observables. Also, M ⊥ (k, t) = 0 for k = Q. We remark that even though the magnetization dynamics is significantly constrained at the level of single spin observables by symmetries, arbitrary spin correlations are allowed to form in the course of evolution, including both in-and out-of-plane spin correlations at arbitrary wavevectors. A simplifying aspect of the present problem is that the apparently broken translation symmetry of the spiral state can be restored using an "unwinding" unitary transformationÛ Q ≡ e i j Q·RjŜ z j under which the spiral state transforms into a uniform x-polarized product state |Ψ 0 =Û Q |sp(Q) = j |→ j . The unwinding transformation, however, transforms the Hamiltonian H →H =Û QĤÛ † Q to an anisotropic Heisenberg model with a Dzyaloshinskii-Moriya term: The translation invariance of the initial state in the spiral frame significantly simplifies the structure of the Spin-2PI equations: G and Σ become local in the real space while D depends only on the distance between the sites. These simplifications hold for arbitrary truncations of Γ int . Additionally, the bosonic self-energy Π becomes local in the real space at the NLO truncation. The magnetization is non-vanishing only along the x-direction in the spiral frame due to the symmetry considerations mentioned earlier. The quantities calculated in the spiral frame can be readily transformed to the lab frame using appropriate rotations. In particular, Eq. (7) gives M ⊥ (Q, t) = (1/2) G 23,> (t, t) with G calculated in the spiral frame. We choose the winding to be along the diagonal direction Q = (Q, Q, Q) hereafter and refer to the spiral winding with the single scalar Q ∈ [0, π].
At the LO level, the spin dynamics is governed by the Bloch equation, Eq. (2). The exchange mean field ϕ c is parallel to the local magnetization at all lattice sites in a spiral state, implying the absence of any dynamics. In other words, the spiral states are fixed points of the mean-field dynamics for all windings Q.
Going beyond the LO dynamics and including the exchange field fluctuations by taking into account the NLO corrections, the spiral state exhibits an intriguing fluctuation-induced relaxation dynamics. States with different windings have different energy densities, along with different strength of in-plane and out-of-plane spin fluctuations, and are found to relax in strikingly different ways. As we discuss below, these factors conspire to give rise to a non-trivial hierarchical relaxation scenario for spiral states lying close to thermodynamically stable orders, exhibiting prethermalization [10], and dynamical arrest resembling glassy systems [65].

A. Relaxation and dynamical arrest of the transverse magnetization
The spiral state for Q = 0 is a fully polarized FM eigenstate of the Heisenberg model and is therefore stationary. The Q = π spiral, on the other hand, corresponds to an uncorrelated Néel state which in three dimensions has a large overlap with the correlated AFM state lying at the upper end of the spectrum of the FM Heisenberg model. As a result, the system is expected to achieve a steady state marked with a finite staggered magnetization after a short course of dephasing dynamics, provided that the generated effective temperature is below the ordering temperature. The evolution of M ⊥ is shown in Fig. 1 (b) for several choices of Q, along with a global surface plot for Q ∈ [0, π] and tJ ∈ [0, 30] in Fig. 1 (c). The stationarity of the FM state (Q = 0) and the rapid settlement of Néel state (Q = π) to a steady state with finite staggered magnetization is observed.
Short-time dephasing dynamics-For all Q, the first stage of dynamics is a short-time relaxation of the form M ⊥ ≈ 1/2 − ν Q t 2 arising from the dephasing between the eigenstates that overlap with the spiral. A straightforward calculation using the short-time expansion The values of ν Q extracted from the numerically obtained M ⊥ are in agreement with the exact result, see Fig. 6. The second stage of relaxation dynamics depends on the winding of spiral and is either directly thermalizing, or exhibits long-lived prethermalized states preceding the true thermalization. We discuss these cases separately.
Spiral states with Q ∼ π/2-Spin spiral states with Q ∼ π/2 have a high energy density with respect to both the FM and the AFM state. Thus, Q ∼ π/2 spiral states overlap with a large number of eigenstates of the Heisenberg model. Such a broad superposition of states lead to fast dephasing which is found to be within a few exchange times. Our results indicate a rapid onset of exponential decay M ⊥ ∼ e −γ Q t with the fastest rate occurring at Q = 0.55(1)π ∼ π/2. Spiral states with Q ∼ 0 and Q ∼ π-A complex multiscale relaxation scenario emerges for spirals with windings tuned to Q ∼ 0 and Q ∼ π, lying close to FM and AFM magnetic orders, respectively. The transverse magnetization exhibits an intermediate plateau for these initial states which appears continuously upon tuning Q, see Fig. 1. The plot of M ⊥ shown in Fig. 1 (b) for Q = 7π/8 displays the intermediate plateau followed by relaxation at later times. As Q is tuned closer toward 0 or π, the lifetime of the plateau increases abruptly and the magnetization comes to a dynamical arrest. We investigate the nature of such long-lived plateaus in more detail in the following sections.

B. Prethermalization vs. Thermalization
Due to the non-integrability of the 3D Heisenberg model, the energy distribution of spin fluctuations is expected to approach a thermal population in the long time limit, according to the eigenstate thermalization hypothesis (ETH) [66][67][68]. We investigate the nature of steady states emerging in the dynamics by calculating the spinspin correlation and response functions, corresponding to the Keldysh (K) and retarded (R) components of the CTP spin-spin correlator χ µν (t, t ), by solving the nonequilibrium Bethe-Salpeter equation (see App. A 3). At thermal equilibrium, these quantities are related via the bosonic fluctuation-dissipation relation (FDR): where T is the effective temperature. Here, ω refers to the Fourier frequency in the time difference t−t in the steady state achieved at long times. Likewise, one can define an effective temperature for the exchange field fluctuations using the bosonic FDR between D K and D R . We refer to the temperatures obtained from χ and D as T spin and T fluct. , respectively. The effective temperatures obtained from the FDR in the steady state are shown in Fig. 3 (a). For all spiral windings Q, we find that FDR is satisfied to an excellent degree for each bosonic mode once the steady state is reached, see Fig. 3 (b). As we discuss below, however, the effective temperature obtained from different modes may disagree with each other. This allows us to distinguish prethermalization from true thermalization.
Thermalization of spiral states with Q ∼ π/2-For a range of spiral wavevectors π/4 < ∼ Q < ∼ 3π/4, the steady state temperatures obtained from all bosonic modes, i.e. The two temperatures are in agreement for spiral windings near Q ∼ π/2, supporting the true thermalization of the system. The temperatures calculated in the prethermalized plateaus Q ∼ 0, π (shaded regions) disagree with each other, and generically differ from the temperature of the true thermal states that emerge at later times. Inset: the temperature kBT as a function of Q (same data as in the main panel) displays a resonance from positive infinite temperature to negative infinite temperature at the classical duality point Q = π/2. (b) The approach of T fluct. to steady state (light to dark) as obtained from fluctuationdissipation relations for Q = π/4 (left) and Q = π (right). The steady state temperatures are shown on the plots.
local and non-local in-and out-of-plane spin and exchange field fluctuations, agree with each other, suggesting the complete thermalization of the system and in accordance with the ETH. Spirals with Q = π/2 flow to an infinite temperature thermal state, which is understood from the duality Q → π − Q, J → −J present in the classical Heisenberg model. This classical duality extends to the quantum Heisenberg model in the high temperature regime. The duality point Q = π/2 further marks the resonance from positive temperatures for Q < π/2 to negative temperatures for Q > π/2, see the inset of Fig. 3 (a). The T < 0 thermal states of the FM Heisenberg model with coupling −|J| corresponds to T > 0 states of the AFM Heisenberg model with coupling +|J|, and vice versa. Negative temperature states naturally arise in isolated systems with bounded energy spectra as legitimate thermal states and occur when the initial energy density lies closer to the upper edge of the energy spectrum.
Prethermalization of spiral states with Q ∼ 0, π-For spiral states with Q ∼ 0, π, where the system develops a prethermal plateau, the effective spin and exchange field fluctuation temperatures disagree, even though the FDR is satisfied well for each mode individually. This finding supports the prethermalized nature of such steady states. It is understood that the temperatures calculated within the prethermal plateau [shown as shaded regions in Fig. 3 (b)] correspond to the effective temperature of individual modes, and not the true thermodynamical temperature. We expect the two temperatures to approach each other at longer times once the system exits the prethermalized plateau and progresses toward a fully thermalized state.
The spiral state with Q = 0 is an exact ground state of the system and FDR yields T = 0 as expected. In contrast, the Q = π state approaches a finite temperature, which is understood by the fact that the uncorrelated Néel state must be "dressed" with spin correlations before the steady state is reached, see inset of Fig. 3 (a). The 3D Heisenberg model exhibits a finite temperature equilibrium phase transition from the disordered paramagnetic phase to the ordered FM or AFM phase, depending on the sign of the exchange coupling J. For the spiral at Q = π, the FDR of the spin fluctuations are not well fulfilled at accessible times while those for exchange field fluctuations are. The temperature extracted from the latter |T fluct. (Q = π)| = 0.82J lies below the the AFM ordering temperature T AFM c = 0.946(1)J. The latter has been obtained from quantum Monte Carlo simulations [69].
The near-thermal distribution of fluctuations in the prethermalized plateaus and the proximity to the thermodynamically stable FM/AFM ordered phases allow us to explain the dynamical arrest: the spiral winding Q sets the energy density of the system and subsequently the effective temperature T (Q) in the steady state. T (Q) approximately dictates the magnitude of spin fluctuations on the top of the spiral states, which we recall are meanfield saddle points for all Q. Depending on Q, T (Q) can either lie below or above the critical transition temperature, T FM c or T AFM c , thereby providing an approximate condition for the stability of the spiral order.

C. Instabilities and Correlations
The dynamical stabilization of spirals near the FM and AFM orders, and consequently the appearance of prethermal plateaus, was understood on the basis of thermodynamical arguments in the previous section. However, even though the spiral states are fixed points of the mean-field dynamical equations, they are unstable and have a tendency to form out-of-plane textures as the energy of spiral states can be reduced by an appropriate out-of-plane tilt. Therefore, in a thermodynamical ensemble where arbitrary out-of-plane fluctuations are allowed, these saddle points fail to give rise to symmetry broken states, leaving Q = 0 and Q = π as the only thermodynamically stable orders in the Heisenberg model. Therefore, the present situation must be regarded from the perspective of quantum dynamics, i.e. the unitary evolution of a pure state |sp(Q) rather than the statistical fluctuations in a mixed thermodynamical ensemble. Here, the system remains in a pure state at all times and the magnetic order is confined to the xy spiral plane due to the symmetries discussed at the beginning of Sec. III. It is therefore conceivable that symmetry-protected dynamical constraints allow thermodynamically unstable saddle points to become long-lived states in the course of unitary dynamics. The evolution of the late-time most enhanced mode for Q = π/4, π/2 (left) and Q = 7π/8, π (right). In cases where the system thermalizes, left column, SU (2) symmetry emerges in the long time limit, while it is broken in the prethermal case Q = 7π/8 and for Q = π, right column. In the latter case, the system can exhibit true long-range order provided its effective temperature is below the critical temperature of the equilibrium phase transition and thus be thermal and simultaneously break SU (2) symmetry.
The out-of-plane instability of the spiral state in the Heisenberg model can be studied either by performing a linear response analysis of the Bloch equations, or similarly from the Holstein-Primakoff spin-wave analysis. Either way, the dispersion of out-of-plane spinwaves forming on the top of the spiral is found as ω k = 2 k − ∆ 2 k [32,70], where: Unstable modes arise when ω k assumes imaginary values. Except for Q = 0, π/2, π, one always finds such unstable modes: for Q < π/2, the fastest growing mode is k = (k, k, k) with k = cos −1 [cos 2 (Q/2)] along with a sharp cutoff |k| ≤ Q; for Q > π/2, unstable modes occur for |k − π| ≤ Q, with the fastest mode always being the staggered k = π mode, independently of Q. A simple estimate for the lifetime of the prethermal plateaus is obtained by calculating the time it takes for typical unstable out-of-plane collective mode to grow to O(1). The rationale behind this estimate is that the in-plane order can not coexist with strong enough out-of-plane fluctuations. Expanding around Q = 0, π, we obtain a scaling t ∼ 1/Q 2 for the FM-like and t ∼ 1/(π − Q) for the AFM-like spirals, up to logarithmic corrections. However, the lifetime of plateaus as found from the Spin-2PI formalism exceeds the above estimates; in particular, as Q is tuned closer toward 0 or π, we observe as rapid increase of the lifetime of prethermal plateaus.
The top panels in Fig. 4 show Im[ω k ] for several values of Q, along with the evolution of equal-time out-of-plane spin correlations iχ zz,K k (t, t) = Ŝ z k (t)Ŝ z −k (t) obtained by solving the non-equilibrium Bethe-Salpeter equation in the NLO approximation, bottom panels. Further, the connected part of in-plane correlations iχ +−,K At t = 0, spin correlations are zero in accordance with the initial spiral state |sp(Q) being an uncorrelated product state. The out-of-plane correlations form at times t ∼ J. The most enhanced correlations coincide with the wavevector predicted by the linear response analysis to a good degree. The sharp cutoffs predicted by this analysis are found to be smeared, which is expected due to the mode coupling embedded in our self-consistent approach. The time scale for the formation of correlations is found to be on the order of the dephasing time, reflecting the fact that the short-time dephasing dynamics and formation of correlations are manifestations of the same phenomenon.
For spiral states that thermalize within the numerically achievable time scales, we observe a smooth shift in both in-plane and out-of-plane spin correlations from the initial Q-dependent enhanced modes to either k = 0 or k = π, depending on whether Q < π/2 or Q > π/2, respectively [see Fig. 4(a), (c), and the inset]. Even though the linear response analysis correctly indicates the wavevector of the fastest growing out-of-plane mode, the spin correlations rapidly saturate to their maximum values, as opposed to an unbounded exponential growth.
A similar rapid dynamical regulation of the growth of unstable modes was previously reported in Ref. [22,45] in the context of parametric resonance in the O(N ) model. In the present problem, this phenomenon explains why the lifetime of the plateaus exceeds the estimate obtained from the linear response analysis, and indicates the important role of nonlinear effects and the necessity of nonperturbative treatments.
For spiral states that exhibit long-lived prethermal plateaus, we study the exchange field correlations D, a quantity that is closely related to χ but can be calculated for much longer times with less computational resources. The evolution of D zz,K k (t, t) and D +−,K k (t, t) for Q = π/4 and Q = 7π/8 are shown in Fig. 5 (a). The former corresponds to a spiral state that thermalizes at about 20J −1 , while the latter exhibits a magnetization plateau up to τ M ∼ 20 J −1 , as shown in Fig. 1. For t < ∼ τ M , the most enhanced in-plane mode occurs at k = Q which upon demagnetization smoothly switches to k = π for t > ∼ τ M . The most unstable out-of-plane mode is always at k = π.
As a further check for thermalizing behavior, we study the restoration of the SU (2) symmetry in the exchange field fluctuations D zz,K k (t, t) → (1/2)D +−,K k (t, t), see Fig. 5 (b). For Q = π/4 and Q = π/2, we find that the SU (2) symmetry is restored at longer times (left column), while for Q = 7π/8 and the Néel initial state Q = π, the SU (2) symmetry remains broken at all accessible times (right column). Notably, the out-of-plane fluctuations for Q = 7π/8 is found to be an order of magnitude stronger than the Q = π, in agreement with the previously mentioned existence of an unstable out-of-plane mode for the former state and its absence in the latter.
The magnitude of in-plane fluctuations remain essentially constant in the plateau for Q = 7π/8 (top right) while the out-of-plane fluctuations monotonically increase and reach a maximum at t ∼ 20 J −1 ∼ τ M , precisely when the prethermal magnetization decays. This finding connects the decay of the the spiral to the growth of out-of-plane fluctuations. The time τ M also marks a reversal in the trend of out-of-plane and in-plane correlations. Even though this change indicates a first step toward establishing SU (2)-symmetric correlations, the condition is far from being satisfied at t ∼ τ M and is bound to occur on much longer time scales, indicating a hierarchical relaxation scenario with the relaxation of magnetization preceding the relaxation of correlations.
The appearance of long-lived prethermal states and the hierarchical relaxation is reminiscent of aging dynamics in classical structural glass models with quenched disorder [65] and kinematically constrained models [71]. Similar multi-scale glassy relaxation dynamics has been recently reported in the quench dynamics of fermions in a nearly-integrable 1D model using a different method [15].
According to the discussions presented so far, the intricate relaxation dynamics of the spiral state is deeply rooted in the quantum nature of spins. In order to further study role of quantum correlations, we compare the predictions of the Spin-2PI formalism with the results obtained from the truncated Wigner approximation (dTWA) [72,73], a variant of the semi-classical TWA method [74]. In Fig. 6 we show a comparison of the magnetization obtained from Spin-2PI (solid black lines) and dTWA (dashed blue lines). The two methods generically agree with the analytic short time expansion (red thick lines), with the exception that dTWA does not reproduce the correct short time dynamics for Q = π/4. The two methods predict strikingly different long time dynamics. Even though dTWA exhibits some degree of dynamical slowing down for FM-like and AFM-like spirals, it produces neither the prethermal plateau for Q = 7π/8, nor the finite steady-state magnetization for Q = π. We note that the latter is supported by exact QMC calculations.

IV. CONCLUSIONS AND OUTLOOK
We formulated a non-perturbative and conserving field theoretic technique for describing the far-fromequilibrium quantum dynamics of strongly interacting spin-1/2 systems for arbitrary lattices and initial states. Referred to as the Spin-2PI formalism, this method systematically improves upon the mean-field description by including quantum fluctuations by means of an asymptotic 1/N expansion, which is controlled in models with intrinsically large lattice coordination number.
We utilized the Spin-2PI technique to study far-fromequilibrium phenomena in spin systems with continuous symmetries. Specifically, we explored the relaxation dy-namics of spin spiral states in the 3D Heisenberg model, treating the spiral winding Q as a tuning parameter. Going beyond the trivial mean-field (LO) dynamics by including the NLO correction, we found the spiral states with different windings to relax in remarkably different ways. In particular, spiral states resembling FM and AFM ordered states, corresponding to Q ∼ 0 and π respectively, get trapped for long times in non-thermal states, i.e. "false vacuums" whose lifetime diverge as the windings are tuned to Q = 0 or π. In contrast, the spiral states far from Q = 0, π relax rapidly.
We calculated the effective temperature of spin and exchange field fluctuations from the fluctuation-dissipation relation. For Q ∼ π/2 spiral states, all modes reach a single temperature, supporting full thermalization in accordance with the eigenstate thermalization hypothesis. In contrast, the different bosonic modes of prethermalizing spirals settle at different temperatures.
We investigated the dynamical formation of correlations and found that the collective modes predicted to be unstable from a linear response analysis, are self-regularized at rather short time scales, demonstrating the importance of the nonlinear effects and non-perturbative treatments. The growth of out-of-plane fluctuations cause the eventual decay of the prethermal states. The restoration of SU (2) symmetry occurs much later after the decay of magnetization, suggesting a hierarchical relaxation reminiscent of coarsening and aging in classical glassy systems. Our results can be tested readily in ultracold atoms experiments with twocomponent Mott insulators in 3D optical lattices, such as a 3D extension of the experiments in Refs. [32,61].
This work can be extended in several directions. A straightforward extension is to investigate the relaxation of spiral states in anisotropic models, or in lower dimensions. Another immediately accessible direction is studying spin systems with long-range interactions, as realized for instance with Rydberg atoms, polar molecules, or trapped ions, and their instability toward dynamic crystallization. In addition, the evolution of disordered states can be studied as well as the formation of topological defects in quenches to the ordered phase, corresponding to the instantaneous limit of the quantum Kibble-Zurek mechanism [75,76]. For the 3D Heisenberg model with SU (2) symmetry, topologically stable hedgehogs [77] are expected to form with universal scaling laws. Other extensions include coupling the system to a bath, taking into account NNLO corrections to analyze the robustness of the NLO results, and comparing with other systematic expansions such as 1/D-expansion in D-dimensional lattices and 1/S-expansion in large-S spin systems. In this Appendix, we provide supplementary material for the Spin-2PI formalism along with a brief account of the numerical methods. The covered material includes the explicit derivation of the approximate dynamical equations from the NLO truncated 2PI effective action and the reconstruction of real-time spin-spin correlators from the Bethe-Salpeter equation.

Correlation functions on the Schwinger-Keldysh time contour
In the Schwinger-Keldysh formalism, the nonequilibrium dynamics of quantum fields is most elegantly derived from a path-integral defined on the round-trip contour C = C + ∪ C − : The Majorana operators η and real vector boson ϕ are replaced by Grassmann and real vector valued variables in the path-integral, along with an anti-periodic and periodic boundary condition at the contour endpoints, respectively.
The correlation functions defined on the C contour can be thought of 2×2 matrices in the two-dimensional space of the contour branch index. For example, the Majorana 2-point correlator G can be explicitly written as: where the times appearing in the matrix are ordinary times. We have dropped the discrete indices for brevity. The off-diagonal matrix elements are identified with the "lesser" and "greater" explicitly ordered correlators: The diagonal matrix elements are related to each other by the virtue of the unitarity of evolution: which are identified with the usual retarded and advanced response functions, G ++ ≡ G R and G −− ≡ G A . While the lesser and greater correlation functions are independent functions for Dirac (complex) fermions, they are related to each other for Majorana fermions by transposition and negation, as it can be seen from Eq. (A3): In summary, the 2-point correlator of Majorana fermions on the contour is fully specified by a single real-time correlator, e.g. G > (1, 2). It is easily shown that the same decomposition and relations hold for the Majorana self-energy Σ. The correlator of real bosons D and the bosonic self-energy Π admit a similar decomposition, except for the absence of the relative minus sign in the definition of D > and D < : which imply: Similar to Majorana correlators, the 2-point correlator of real bosons on the contour is fully specified by a single real-time correlator, e.g. D > . The same result holds for the bosonic self-energy Π.

Feynman rules for the Spin-2PI formalism
The conventional Feynman diagram rules are used for interpreting the diagrams appearing throughout this work: The integer indices refer to the bundle of lattice site and contour time in the diagrams above. Since the Majorana fermion propagators possess no charge flow direction, one may arbitrarily assign a direction to each line. The overall sign of each diagram, however, must be determined at the end by counting the number of fermion permutations.
The power counting of the large-N extension is performed as follows: (1) each Majorana fermion loop introduces a factor of N resulting from the replica summation, (2) each interaction and boson line introduces a factor of 1/N .
The vacuum diagrams accompany symmetry factors which must be worked out case by case. The self-energy Σ, Π and 4-point vertex Λ (2) diagrams have an extra factor of i and i 2 , respectively.

Evolution of correlations functions in the
Spin-2PI formalism The transition from the path-integral to the 2PI effective action Γ[G, D, ϕ] was briefly outlined in the main text and is a straightforward generalization of the results of Cornwall, Jackiw, and Tomboulis [37]. Within this formalism, the evolution equations follow from a variational principle, reminiscent of Lagrangian dynamics of classical particles, with the quantum correlators playing the role of generalized coordinates. Going back to Eqs. (8a) and (8b) and making Γ stationary with respect to G, D, and ϕ, we obtain: With the spin, replica, time, and space indices laid out explicitly, the "bare" fermion and boson propagators are written as: respectively. The contour Dirac δ-function is defined as δ C (t 1 , t 2 ) = ±δ(t 1 − t 2 ) with the ± sign corresponding to t 1 , t 2 ∈ C ± , respectively. In Eq. (A7a), M [ϕ c ](1, 2) = −iδ σ1σ2 δ C (t 1 , t + 2 ) δ j1j2 ϕ µ c,j1 (t 1 ) ε µα1α2 the LO interaction effect and describes the coupling of Majorana fermions with the classical spin mean-field ϕ c . According to Eq. (A7c), the latter is instantaneously determined by the Majorana tadpole contracted with a bare interaction line. Thus, we find: which resembles the familiar Hartree self-energy that describes the mean-field effects. We emphasize that the Majorana tadpole is identified with the magnetization in our formalism. Note that the M (1, 2) ∝ δ C (t 1 , t + 2 ) is instantaneous and carries no memory effect. We will later show that truncating the approximation at this level and neglecting the self-energies indeed yields the Bloch equation.
Going beyond the LO approximation, Σ[G, D] and Π[G, D] describe memory effects associated from the spatiotemporal fluctuations of the exchange field. By definition, these self-energies are obtained from the variations of Γ 2 [G, D]: We recall that Γ 2 [G, D] is formally equivalent to the sum of 2PI vacuum diagrams constructed from the interaction vertex L int [η, ϕ] = i 2 ε αβγ ϕ α j η β;σ j η γ;σ j and admits a systematic expansion in 1/N . Here, we truncate the series at the NLO level: where Π µν 0 (1, 2) = i ε µαβ ε νγλ N σ=1 G α;σ,γ;σ j1j2 (t 1 , t 2 ) G β;σ,λ;σ j1j2 (t 1 , t 2 ) is the Majorana bubble. Since D ∼ 1/N and the factor of N resulting from the replica summation in the Majorana bubble, we find Γ NLO 2 ∼ O(1). This must be compared to the LO term in Γ int which is O(N ). The resulting NLO self-energies are given as: Having derived the explicit expressions for the selfenergies, we discuss the derivation of evolution equations as the next step. Our starting point are the coupled Dyson's equations given in Eqs. (A7a) and (A7b). Strictly speaking, Dyson's equations are differential identities on the contour functions. They can be cast into a more useful form by acting them from the left and right hand side by G and D, respectively, resulting in a set of contour integro-differential equations: We have defined shorthand δ(1, 2) ≡ δ j1j2 δ α1α2 δ C (t 1 , t 2 ) and the contour integral Eq. (A13a) and its adjoint Eq. (A13b) are referred to as Kadanoff-Baym (KB) equations. The convolution integrals of selfenergies and correlators manifestly show memory effects, which is a shared feature of beyond mean-field approximations.
The spatial structure of Eqs. (A13a)-(A14b) can be simplified by noting that physical initial states imply initial correlations between pairs of Majorana operators on the same site, i.e. G α1α2 j1j2 (t 0 , t 0 ) ∝ δ j1j2 . Crucially, this property extends to all times in the KB dynamics, independent of the order of truncation in 1/N . To see this, one first establishes that the assumption G j1j2 (t 1 , t 2 ) ∝ δ j1j2 for t 0 ≤ t 1 , t 2 ≤ T implies Σ j1j2 (t 1 , t 2 ) ∝ δ j1j2 in the same domain. The causal structure of Eqs. (A13a)-(A13b) subsequently extends this property to an infinitesimally larger domains, and eventually to all times by induction. Therefore, we can always make the following simplifying substitution in the KB equation: The LO approximation: The KB equations reduce to the mean-field Bloch equation upon truncation at the LO level which amounts to neglecting fluctuation self-energy corrections Σ → 0. In this limit, the KB equations imply: i∂ t1 G α1α2,> jj (t 1 , t 2 ) + iϕ ν c,j (t 1 ) ε να1µ G µα2 jj (t 1 , t 2 ) = 0, −i∂ t2 G α1α2,> jj (t 1 , t 2 ) + iϕ ν c,j (t 2 ) ε νµα2 G α1µ jj (t 1 , t 2 ) = 0.
Subtracting the equations from one another, setting t 2 = t 1 = t and using Eq. (7), we finally obtain: which is the Bloch equation as anticipated.
The NLO approximation: Including self-energy corrections, the time convolutions appearing in the KB equations prohibit us from arriving at a closed equation for the equal-time Green's functions and we inevitably need to solve for the complete unequal time Green's function. For concreteness, we consider the case of spin spirals hereafter. The spatial structure of the KB equations can be significantly simplified by applying the unwinding unitary transformation, either directly on Eqs. (A13a)-(A14b) or on the spin Hamiltonian. Either way, the initial spiral state transforms into an uncorrelated xpolarized FM state |Ψ 0 ≡ j | → j at the expense of an anisotropic interaction (see Eq. 13). The 2-point correlator of Majorana fermions at t = t 0 is easily found as: The exchange field correlator at t = t 0 is not an independent degree of freedom and is determined by G(t 0 , t 0 ), see Eq. (A12). For translationally invariant initial states as such, G and Σ further become independent of the lattice site. Furthermore, D j1j2 depends only on the distance and at the NLO level, Π j1j2 is local as well. The simplified structure of the correlators and self-energies is summarized as follows: The KB equations can be written explicitly in terms of G > and D > using the Langreth rules [78]. We quote the final result, setting N = 1: The self-energies Σ > and Π > are read from Eq. (A12). The last explicit equations are suitable for devising a numerical forward propagation scheme. Starting from G 0 (t 0 , t 0 ), we calculate Σ(t 0 , t 0 ) and Π(t 0 , t 0 ) from Eq. (A12), and D(t 0 , t 0 ) from Eq. (A20a). The casual structure Eqs. (A19a)-(A20b) allows us to propagate {G, Σ, D, Π} in (t 1 , t 2 ) in discrete time steps of size ∆t. This is achieved using a robust predictor-corrector method with guaranteed accuracy to O(∆t 3 ).
The first step in solving the BS equation is to recast it in terms of functions of ordinary times. In comparison to the KB equation, this step is significantly more involved here due to the complex real-time structure of 3-time and 4-time CTP functions and multiple contour integrals. We leave the cumbersome details for a separate publication and solely outline the procedure here. We showed earlier in Sec. A 1 that the 4 real-time matrix elements of 2-time functions such as G and D can be fully specified using a single real-time function, e.g. G > . A similar analysis of Γ α1α1;µ q (t 1 , t1; t 2 ), taking into account symmetries and unitarity of evolution, reveals that the 8 real-time components of 3-time function as such can be fully specified by 3 independent functions. Accordingly, the contour BS equation can be explicitly written as 3 coupled two-dimensional integral equations in ordinary times; the latter is numerically solved by discretizing the integrals using approximate quadratures and solving the resulting linear system.