Non-thermal melting of Neel order in the Hubbard model

We study the unitary time evolution of antiferromagnetic order in the Hubbard model after a quench starting from the perfect N\'eel state. In this setup, which is well suited for experiments with cold atoms, one can distinguish fundamentally different pathways for melting of long-range order at weak and strong interaction. In the Mott insulating regime, melting of long-range order occurs due to the ultra-fast transfer of energy from charge excitations to the spin background, while local magnetic moments and their exchange coupling persist during the process. The latter can be demonstrated by a local spin-precession experiment. At weak interaction, local moments decay along with the long-range order. The dynamics is governed by residual quasiparticles, which are reflected in oscillations of the off-diagonal components of the momentum distribution. Such oscillations provide an alternative route to study the prethermalization phenomenon and its influence on the dynamics away from the integrable (noninteracting) limit. The Hubbard model is solved within nonequilibrium dynamical mean-field theory, using the density matrix-renormalization group as an impurity solver.


I. INTRODUCTION
Ultra-fast pump-probe experiments on condensed matter systems and experiments with cold gases in optical lattices have opened the intriguing possibility of controlling transitions between complex phases on microscopic timescales. This has motivated intensive theoretical efforts to understand fundamental aspects of the dynamics in interacting many-body systems, and leads to predictions in marked contrast to the naive expectation that interactions imply rapid thermalization [1]: Integrable systems can keep memory of the initial state for all times and relax to a generalized Gibbs ensemble [2,3], but also away from integrability thermalization can be delayed by prethermalizaton [4][5][6][7], and one can identify regimes of different dynamical behavior which are clearly separated by non-thermal critical points [8][9][10][11][12][13].
Of particular interest with respect to complex phases in condensed matter is the dynamics of symmetry brokenstates [14][15][16]. While the relevant relaxation mechanisms after a perturbation are hard to disentangle in a solid, cold atoms in optical lattices provide a versatile platform to investigate isolated quantum systems in ideal situations. The preparation of thermodynamic long-range ordered phases in cold atoms is still a challenge [17], but advanced techniques for lattice design have made it possible to prepare an ordered state on a lattice of isolated sites, * Electronic address: martin.eckstein@mpsd.cfel.de and to probe its dynamics after tunneling between the sites is switched on [18][19][20][21]. In the following we consider such a setup for the Fermi-Hubbard model, a paradigm model for emergent long-range order in condensed matter systems. We will simulate the time evolution starting from a classical Néel state in which neighboring lattice sites of a bipartite lattice are occupied with particles of opposite spin.
In general, one can anticipate fundamentally different pathways for melting of long-range antiferromagnetic order in the weakly and strongly interacting Hubbard model: For strong interaction, long-range order arises from antiferromagnetically coupled local moments, which emerge when charge fluctuations are frozen. Magnetic order could thus possibly melt via the destruction of the local moments themselves, through a reduction of the effective exchange interaction [22] (while moments persist), or, along a quasi-thermal pathway, by the transfer of energy from excited quasiparticles (hot electrons) to spins. The latter mechanism is intensively studied in the context of photo-carrier relaxation in high-T c cuprates [23][24][25][26][27][28], where the investigation of the spin-charge interaction challenges the limits for the time-resolution in state of the art pump-probe experiments [29][30][31]. For weak interaction, on the other hand, quasiparticle states may be important to understand relaxation processes. In the paramagnetic phase the conservation of the quasiparticle momentum occupations imposes constraints on the dynamics which can lead to prethermalization [4,5,9,[32][33][34]. Prethermalizaton, which was recently observed in a one-dimensional Bose gas [7], has been suggested to be a universal feature of near-integrable systems [6], but previous predictions for the Hubbard model rely on a discontinuity of the momentum distribution which is absent at nonzero temperature and thus experimentally hard to observe. Here we show that the symmetry-broken initial state provides an alternative perspective to investigate this physics and its breakdown far from integrability.
Quenches from a Néel state have been explored in quantum spin models [10,[35][36][37], also as a way to prepare ordered states in the Hubbard model [38], but a pure spin model cannot describe the relevant dynamics of charge excitations and local moments. The Hubbard model has been studied in one dimension using the density-matrix renormalization group (DMRG) [39]. For the dynamics of lattice fermion models in more than one dimension, nonequilibrium dynamical mean-field theory (DMFT) [40] is the most promising approach. Quenches within the antiferromagnetic phase of the Hubbard model at strong-coupling [24] are in line with the "quasi-thermal" pathway discussed above. The regime of intermediate interactions, where the notion of local moments becomes ambiguous, or to weak coupling, where prethermalization may be expected, has been elusive so far. Previous numerical solutions of the DMFT equations were based on the self-consistent strong-coupling expansion [41] or weak-coupling impurity solvers [13,42,43], which both fail at intermediate coupling, while weak-coupling quantum Monte Carlo studies [9,41] are most efficient for noninteracting initial states and restricted to short times. In this work we overcome these limitations using a recently developed Hamiltonian based formulation for the impurity model of nonequilibrium DMFT [44], which has opened the possibility to use wave-function based techniques to solve the DMFT equations [45,46]. Here we use DMRG as an impurity solver [46], which allows us to reach sufficiently long times in the evolution to address the above issues.

II. MODEL AND METHODS
Throughout this work we consider the single-band Hubbard model at half-filling, with nearest-neighbor hopping J and on-site Coulomb repulsion U . The Hamiltonian is given by where c † iσ (c iσ ) are electron creation (annihilation) operators for lattice site i and spin σ, and n iσ = c † iσ c iσ . The model is solved using nonequilibrium DMFT [40], for a Bethe lattice in the limit of infinite coordination number Z and hopping J = J * / √ Z, where the approach becomes exact [47]. The energy unit is set by J * = 1, and time is measured in inverse energy, i.e, the free density of states is given by D(ǫ) = √ 4 − ǫ 2 /(2π). To simulate the quench, we choose a time-dependent hopping J * (t) = 0 for t ≤ 0 and J * (t) = 1 for t > 0. For t ≤ 0, the system therefore consists of a set of isolated lattice sites, which are prepared in a classical Néel state, where A and B are sub-lattices of the bipartite Bethe lattice.
In DMFT, the lattice model is mapped to a set of impurity problems, one for each inequivalent lattice site j = A, B, with a time-dependent hybridization function ∆ jσ (t, t ′ ). (In this expression, time arguments lie on the Keldysh contour; see Ref. [40] for a detailed description of nonequilibrium DMFT and the Keldysh formalism.) For the Bethe lattice, the latter is determined self- is the local Green's function. To solve the impurity model with a non time-translationally invariant hybridization function, we derive an equivalent representation in terms of a timedependent Anderson impurity Hamiltonian [44] with up to L = 24 bath orbitals, from which the time-dependent Green's functions are computed using a Krylov timepropagation for matrix product states [46]. The Hamiltonian representation of the DMFT impurity model is exact for small times, but an increasing number of bath sites is needed to reach longer times [48]. We verify the convergence of the solution with the bath size L. Up to L = 12, the results have also been cross checked with a Krylov time-propagation in the full Hilbert space. For further details of the numerical solution, see Appendix A.
III. RESULTS Figure 1 shows the time evolution of the antiferromagnetic order parameter M (t) and the double occupation d(t) = n ↑ (t)n ↓ (t) after the quench, for various values of the Coulomb interaction. In order to account for the trivial reduction of the local spin expectation value by virtual charge fluctuations, we define M (t) as the staggered order M stagg ≡ n A↑ (t) − n A↓ (t) = n B↓ (t) − n B↑ (t) , normalized by the probability P 1 for a site to be singly occupied, P 1 (t) = 1 − 2d(t). To test whether the system thermalizes after the quench, we compare to an equilibrium state at the same internal energy (which is zero for the Néel state). The corresponding effective temperature T eff (Fig. 1c) lies above the Néel temperature T Néel for all values of U [49]. This implies a paramagnetic state after thermalization. While M (t) indeed continues to decay throughout the simulated time interval, the double occupancy saturates to a non-thermal value for U 4 (arrows in Fig. 1b point at the thermalized value d(T eff )), in agreement with earlier studies on the lifetime of doublons in the paramagnetic Mott regime [51][52][53]. Already at a first glance, the relaxation of M (t) and d(t) therefore suggests different mechanisms for small and large values of U , with a rapid and oscillatory decay of M (t),  [54] as an impurity solver). All thermalized states are in the paramagnetic phase, the corresponding inverse temperature 1/T eff is plotted in c). and a long-lived non-thermal state, respectively. In the following we will analyze the two regimes in more detail.

A. Weak-coupling: Residual quasiparticles
For quenches to small U the Hamiltonian is close to the integrable point U = 0. This suggests to study relaxation in terms of the momentum occupation n k (t) = c † k c k , which is conserved at U = 0. For a state with translational symmetry breaking, the single-particle density matrix ρ kk ′ (t) = c † k ′ (t)c k (t) is no longer diagonal in momentum k. (The discussion holds for a general lattice like the Bethe lattice when k denotes the eigenstates of the translationally invariant hopping matrix.) For nearest neighbor hopping on a bipartite lattice, eigenstates come in pairs k,k with single particle energy ǫ k = −ǫk, where the wave functions fork and k differ by a staggered phase ξ i = ± for i ∈ A(B), and ρ kk = 0 if the symmetry between sub-lattices is broken. (On the cubic lattice, k andk = k + (π, π, . . .) are momenta related by the antiferromagnetic nesting vector.) In Fig. 2 we plot the diagonal and off-diagonal components of the single-particle density matrix in terms of the two functions n(ǫ k , t) = c † k (t)c k (t) and π(ǫ k , t) = c † k (t)ck(t) , which depend on k only via ǫ k due to the locality of the self-energy within DMFT. In the thermalized state, π(ǫ) = 0 because the state does not break the sublattice symmetry, while in the localized initial state, n(ǫ, t = 0) = π(ǫ, t = 0) = 1/2. In agreement with the behavior of the double occupancy, n(ǫ, t) does not thermalize at large U (thermalized values π T eff (ǫ) and n T eff (ǫ) are shown by solid lines). For U ≤ 2, however, differences between n(ǫ, t) and n T eff (ǫ) become tiny. This is in stark contrast to the behavior of the paramagnetic system after a quench from U = 0, where prethermalization manifests itself precisely in the difference between n(ǫ, t) and n T eff (ǫ) [5]. Moreover, around U = 3 the relaxation of π(ǫ, t) changes from an oscillatory to a monotonous decay. (In Fig. 2, we plot the real part of π(ǫ, t), the imaginary part shows a similar crossover from oscillatory to non-oscillatory behavior.) This observation may be explained following the perturbative arguments of Refs. [5,6]. To second order in U , the Hamiltonian (1) is unitarily equivalent to a model with a finite residue R k , where incoh. denotes incoherent contributions, i.e., an admixture of particle-hole excitations to higher order in U . Hence the momentum occupation is given by k (the coherent part) is unchanged by the time evolution to second order in U . Back-transforming gives n k (t) = R 4 k n k (0)+incoh., where the incoherent contribution is a smooth function of k. For quenches in the paramagnetic phase, n k (t) thus preserves the initial discontinuity at the Fermi surface, which can be taken as a measure of prethermalization [5]. In the symmetry broken state, however, n k (0) is independent of k, and thus n k (t) does not clearly exhibit the existence of residual quasiparticles. In fact, the numerical results suggest that the incoherent part can accurately be described by a thermal distribution. In contrast, a similar argument for the off-diagonal component shows that π k (t) = R 2 k c † k (t)ck(t) + incoh. = R 4 k e i2ǫ k t π(ǫ, t = 0) + incoh., where we have used the time evolution of the quasiparticle,c k(k) (t) = e ∓iǫ k tc k(k) (0) and R k = Rk. Hence we find that the residual quasiparticle dynamics leading to prethermalization close to the integrable point U = 0 can be studied very conveniently with the symmetry broken initial state in terms of oscillations in the off-diagonal components of the momentum occupation.
Similar to the interaction quench in the paramagnetic phase [9,11], we find that the "prethermalization" regime in which residual quasiparticles dominate the dynamics is limited to small interactions; at large interactions, π(ǫ, t) relaxes to zero monotonously (see the U = 6 data in Fig. 2b) and the distribution becomes flat over the Brillouin zone. Below we will see that the dynamics at large U can be analyzed in terms of well defined localized moments. In contrast to the quench in the paramagnetic phase, the crossover between the weak and strong coupling regimes is relatively smooth and occurs between U = 2 and U = 3: In Fig. 3a, we exemplarily plot Re π(ǫ, t) for fixed ǫ = −1.5 and various U . For U 2, k c k of the momentum occupation (panels a, c, and e) and off-diagonal component Re π(ǫ k ) = c † k ck (panels b, d, and f), where k andk are pairs of single particle states coupled by a staggered potential, plotted for quenches to three different values of U as a function of the energy ǫ k ranging from −2 to 2 in the band of the Bethe lattice. The bold lines indicate momentum distributions obtained in the (paramagnetic) equilibrium state at the same energy.
Amplitudes of the background, b, the oscillations, a, the quasi-particle energy ǫ ′ and the quasi-particle decay rate Γ as a function of U .
the curves can be accurately fit with decaying oscillations f 1 (t) = a exp(−Γt) cos(−2ǫ ′ t + φ), where in agreement with the discussion above the quasi-particle energy ǫ ′ → ǫ, and Γ ∼ U 2 for U → 0 (solid lines in Fig. 3a, fit parameters in Fig. 3b). For U 3, on the other hand, a good fit is a monotonously decaying curve 3, there is a crossover between the two behaviors, as evidenced by the dependence of the amplitudes a and b of the monotonous and the oscillating component on U (Fig. 3b).
Before discussing the strong-coupling regime, we note that off-diagonal momentum distributions can in principle be measured by a modified time-of-flight measurement, if before releasing the cloud, one would switch off the tunneling and the interaction, switch on a staggered potential which is ±∆ on the A and B sub-lattice, respectively, and evolve for a given time t m . Time of flight measures the regular momentum occupation n k = c † k c k after that procedure. In the k,k basis, the staggered potential is given by H = ∆ k (c † k ck + H.c.), so that n k (t + t m ) = n k (t) cos 2 (t m ∆)nk(t) + sin 2 (t m ∆) + sin(2t m ∆)Im π k (t) after propagation in the pure staggered potential from time t to t + t m , and Im π k (t) can be extracted.

B. Dynamics of local moments
In a Mott insulator at large U one can expect the existence of well-defined local moments. It is an intriguing question whether these moments persist in the quenched state while the long-range order disappears, and to what extent the crossover in relaxation behavior from weak to strong coupling can be characterized in terms of these local moments. In the following, we propose a simple experiment to distinguish the existence and strength of moments in the quenched state: one spin in the initial Néel state on a given site (the probe site "o") is flipped to the x-direction (see Fig. 4d, inset). Choosing o on the A-sublattice of the Néel state, the initial state (2) of the dynamics is changed to ( In a perfect local moment picture, the spin should then precess in the exchange field of it's neighbors. The inhomogeneous setup with one probe spin can be solved within DMFT, where it corresponds to a modified impurity problem at site o, while the rest of the lattice is unchanged (see App. A). Figure 4a-c shows the local spin expectation values S z , S y , and S x at site o for various values of the interaction. In Fig. 4d we show the trajectory of the spin in the S x -S y plane, starting from S x = 1, S y = 0 at time t = 0. For large U one can indeed observe a precessional motion in the S x -S y plane, as expected for a local moment subject to an exchange field in the z-direction. For U = 0, on the other hand, the spin-dynamics is entirely longitudinal, showing no sign of well-defined local moments. (For U = 0, the dynamics can be solved analytically, yielding S x,o = J 1 (t) 2 /t 2 , while S o,y = 0 for the Bethe lattice at Z = ∞, where J 1 (x) is the first Bessel function (Appendix C)). There is a crossover between the two relaxation regimes.
Although the exchange interaction is in principle not an instantaneous interaction on the timescale of the electronic hopping [22], it is illustrative to quantify the precession dynamics in terms of an effective exchange field. For this purpose we follow Refs. [22,55] and define B eff such that S(t) satisfies the equation of motion d/dt S(t) = B eff × S(t) . We can assume that B eff = B effẑ acts only in the z direction (parallel to the order parameter M stagg on the neighboring sites), and use the parametrization B eff = J ex M stagg to define an effective exchange interaction J ex ; the latter is then given by J ex =φ(t)/|M stagg | where φ(t) = atan(S y (t)/S x (t)) is the angle of the spin in the x-y plane. The resulting value J ex is plotted in Fig. 4e. For large U , J ex shows very good agreement with the perturbative value of the exchange in the Hubbard model, 4J 2 * /U , and is not substantially decreasing with time even for quenches at intermediate interaction (U ≈ 4) where the order parameter quickly decays to zero (see Fig. 1). Equilibrium estimates of the local moment in the intermediate coupling regime (Fig. 4f) furthermore show tendencies of moment formation at elevated temperatures, which may explain why some spin precession occurs even for U = 2. The combination of these results shows that the melting of long-range order proceeds by the "quasi-thermal" pathway discussed in the introduction, i.e., a disordering of exchange-coupled moments, rather than by a change of the exchange interaction or a destruction of the moments.

C. Strong coupling: spin-charge interaction
At large U , a quench within a Mott insulator freezes virtual charge fluctuations, leaving behind a certain density n δ of long-lived mobile carriers [24]. The mechanism for the decay of the antiferromagnetic order is thus expected to be the transfer of energy from excited quasiparticles to the spins, which is currently intensively investigated in condensed matter pump-probe experiments. Although this mechanism is rather well understood in To investigate the decay of long-range order systematically, one has to vary the excitation density. Here we use a quench protocol where in addition to switching on the hopping at time t = 0, the interaction is changed to an intermediate interaction value U i for a short time 0 ≤ t ≤ 0.5, before it is set to the final value U for t > 0.5. (Note that various other protocols, such as an intermediate time-dependent modulation of the hopping, would have the same effect.) Small values U i lead to a larger double occupancy (Fig. 5a), and indeed also a more rapid decay of M (t) (Fig. 5b). We also note that an exponential fit M (t) ∼ ae −t/τ + b would be consistent with a threshold behavior in which M (t) extrapolates to a finite value b for small excitation density (U i close to U i = 8) and to b = 0 for large excitation density, consistent with earlier quench studies based on the non-crossing approximation impurity solver [24], but the times are not sufficient to analyze this long-time behavior in detail.
For a quantitative analysis of the short time behavior, we determine the number n δ of doublons and hole carriers in the quenched state (Fig. 5c inset). Due to virtual charge fluctuations, n δ is not exactly given by an instantaneous expectation value d(t) in the Hubbard model, and we compute n δ from the total weight in the upper Hubbard band (App. B) [56]. For small times, the curves M (t) for various values U i can then be scaled on top of each other by plotting [1 − M (t)]/n δ (Fig. 5c). Such a scaling implies that the number of flipped spins, 1 − M (t), is proportional to the number of carriers. This is consistent with the picture that spin-flips are inserted by mobile carriers, which are initially localized and thus act independently up to times depending on n δ . For large times there is a deviation from the scaling due to the gradual melting of the order parameter.
To further corroborate this picture, we analytically compute the spin-flip rate per carrier in the low density limit from the behavior of a single carrier which is initially localized at a given site 0 in a Ising spin background. Following Ref. [28], we omit the transverse dynamics of the spins, which is on the timescale of J ex and much slower than the hopping, and keep only the z-component of the exchange coupling J ex (t-J z model). The model can then be reduced to a tight-binding model for a single particle on the lattice, with effective Hamiltonian H = −J * / √ Z ij c † i c j + ZJ ex /2 j |j|c † j c j ; the number of flipped spins is simply given by the displacement |j| from the origin, and the second term in the Hamiltonian accounts for the corresponding exchange energy cost, i.e., the particle is bound to the origin by a linear potential due to the "string" of flipped spins left behind [57]. The dotted line in Fig. 5c shows the mean displacement R(t) of the particle in this model, which indeed coincides with the mean number of flipped spins per particle in the numerical DMFT results. As evident from a comparison of the two curves for J ex = 0 and J ex = 0.5 (the perturbative value for the Hubbard model at U = 8, see also Fig. 4e), the effect of J ex becomes important only at longer times (when numerical data already depend on n δ ), because initially the kinetic energy of the carrier is much larger than J ex . In order to measure the effect of J ex on the charge-carrier interaction, one would have to reduce the number of excitations (e.g., by switching on the hopping slowly), which however makes an accurate determination of n δ increasingly difficult.

IV. CONCLUSION
In conclusion, we have studied the short-time relaxation dynamics of the Néel state in the single-band Hubbard model by means of nonequilibrium DMFT, using DMRG to solve the quantum impurity model. We find qualitatively different relaxation behaviors for weak and strong interactions, separated by a crossover around U ≈ 0.6×bandwidth: For strong interaction, local magnetic moments persist while their order is destroyed by spinflips due to the hopping of mobile charges. The latter resembles the femtosecond carrier spin interaction which is relevant for the dynamics of photo-induced states in high-T c cuprates [31]. To demonstrate the persistence of local moments we proposed a spin precession experiment, which could be implemented similar to the proposed measurement of dynamic spin-spin correlation functions in equilibrium [58]. At weak interaction, the dynamics of the Néel state is governed by almost conserved quasiparticles, which are also the origin for prethermalization in nearly integrable systems [4,6,7]. In the symmetrybroken state, the breakdown of these quasiparticles away from integrability leads to a crossover from oscillatory to non-oscillatory relaxation behavior, which can provide a clear experimental signature that does now rely on a quantitative comparison to the thermal equilibrium state.
Our simulations within DMFT are exact in the infinite dimensional limit, and it is thus interesting to compare to recent results for one dimension [39]. Similar to our results, in d = 1 one finds a rapid saturation of the double occupancy and a slower dynamics of the order parameter at large U , but the decay of antiferromagnetic order is of different origin: In large dimensions, the fastest melting processes after the quench take place on the timescale of the hopping due to the strong charge-spin interaction, while the latter is absent in d = 1 so that the dynamics happens on the timescale of the exchange interaction [39]. The quasiparticle physics at weak coupling and in the crossover regime has not been addressed in Ref. [39], but based on the perturbative argument given above the signatures in the off-diagonal components of the momentum distribution should persist also in lower dimensions. (Also in the paramagnetic case, a long-lived jump in the momentum distribution function is found in d = 1 [32,34,59] and d = 2 [33,34].) Quench experiments starting from the Néel state have recently been performed with noninteracting fermions in one dimension [20], and bosons in two-dimensions [21]. Hence this setup should be a feasible approach to study fundamental aspects of the decay of antiferromagnetic long-range order in the paradigmatic Hubbard model. Moreover, on the numerical side our work emphasizes the high potential of DMRG as an impurity solver for future applications of nonequilibrium DMFT, to explore the intermediate coupling regime which is inaccessible by weak or strong coupling perturbation theory.
To simulate the dynamics of a lattice model which is initially in equilibrium at temperature T = 1/β, we adopt the formulation of dynamical mean-field theory within the Keldysh framework (nonequilibrium DMFT), for an L-shaped time contour C which extends from initial time t = 0 to a maximal time t max along the real time axis, back to time 0, and along the imaginary time axis to −iβ. For a general description of the formalism, as well as the notation and definition of contour-ordered functions, we refer to Ref. [40]. In this appendix we summarize the specific setup for the quench from the Néel state, and the solution of the DMFT equations using DMRG.
In DMFT, the lattice model is mapped to a set of impurity problems, one for each inequivalent lattice site j, with time-dependent hybridization functions ∆ jσ (t, t ′ ). The action of the impurity model is given by on the Keldysh contour C, which yields the local contour ordered Green's function G jσ (t, t ′ ) = −iTr[T C e Sj c σ (t)c † σ (t ′ )]/Z. The hybridization function ∆ jσ (t, t ′ ), must be defined self-consistently. For the Bethe lattice, one has [50] ∆ jσ (t, where the sum runs over nearest neighbors of j. In the antiferromagnetic state, all sites on the A and B sublattices are equivalent, respectively. With the additional symmetry G A,σ = G B,−σ only one impurity model must be solved with ∆ σ (t, t ′ ) = J * (t)G −σ (t, t ′ )J * (t ′ ), where we used the scaling J(t) = J * (t)/ √ Z with the coordination number Z. For the initial product state with J * (t) = 0 for t < 0, ∆(t, t ′ ) = 0 if one time argument is on the imaginary branch of C. Furthermore equivalence under a simultaneous spin and particle-hole transformation implies the symmetry To compute the Green's function we follow Ref. [44] and map the impurity model to a time-dependent Anderson Hamiltonian in which the impurity is coupled to L bath orbitals (p = 1, ..., L). The parameters V pσ (t) and ǫ p are determined such that the local Green's functions obtained from (A1) and (A4) are identical. As derived in Ref. [44], for J * (t < 0) = 0, one can choose V pσ (t) = 0 for t < 0, and the mapping condition is satisfied by (assuming L even) where ǫ pσ = 0, and the bath orbitals p = 1...L/2 and p = L/2 + 1...L are initially doubly occupied and empty, respectively. Equations (A5) and (A6) are solved by a Cholesky fit of the real-time matrix ∆(t, t ′ ), which quickly converges for small times with the number of bath orbitals required [48]. Due to the symmetry (A3) we use The impurity site is initially occupied with a spin σ =↑ (for a site on the A sublattice), i.e., the initial state for the impurity model is a product state |Ψ imp,A = c † ↑ L/2 i=1 a † p↑ a † p↓ |0 , and the Green's function is obtained by solving where time evolution is determined by (A4). We use a Krylov time-propagation for matrix product states [46] with up to L = 24 bath orbitals.

Inhomogeneous setup
For the inhomogeneous setup we assume that in the initial state on the lattice the spin at one site o of the lattice is flipped in the x direction. Without loss of generality we assume that o is on the A sub-lattice. From the self-consistency equation (A2) one can see that the hybridization on all other sites differs from the homogeneous case only in order 1/Z, i.e. for Z → ∞ the back-action of the probe site on the rest of the lattice can be neglected. On the probe site we solve an impurity problem with the same (nonequilibrium) hybridization function ∆ A as on all remaining A-sites, i.e., an impurity problem (A4) wit the same parameters V pσ , but with a different initial state, (A10)
In the translationally invariant case (no probe site), we also determine diagonal and off-diagonal components of the momentum occupations n(ǫ, t) and π(ǫ, t), which are obtained from the momentum resolved Green's function (for the definition of k andk, see the main text) .
(A11) (Here and in the following, bold-face quantities denote 2 × 2 matrices and we omit spin indices for simplicity).
The self-energy is local in space but depends on the sublattice and spin; in the k,k representation it thus assumes the (2 × 2) form so that G ǫ is obtained from the lattice Dyson equation where the dispersion in the k,k representation reads ǫ = ǫτ z because ǫ k = −ǫk. The components Σ j of the self-energy (j = A, B) are obtained from the impurity Dyson equation (i∂ t +µ−∆ j −Σ j ) −1 = G j . In praxis, we solve an integral equation Appendix B: Mobile carrier density in the excited state.
In the Mott insulating phase of the Hubbard model, a well-defined measure for the number of doublon or hole carriers is given by the total occupied spectral weight in the upper Hubbard band and the total unoccupied weight in the lower Hubbard band, respectively. The double occupancy, in contrast, depends on virtual charge fluctuations which are nonzero also in the insulating ground state. Specifically, we define the occupied density of states as the partial Fourier transform is the local Green's function, and δ = 1.5 ensures a smooth cutoff (which does not influence the results unless its inverse width is longer than the inverse of the gap). The spectrum N σ (ω, t) is plotted in Fig. 6a for two different times, for the same quench parameters as in Fig. 5 of the main text. The right panel shows the integrated density W σ (t) = ∞ 0 dωN σ (ω, t). While the weight in the upper and lower band differs considerably between majority and minority spin, the integrated weight W σ (t) reflects the doublon density and is thus independent of σ. It is interesting to point out that as a function of time spectral weight is both redistributed between the lower Hubbard bands of the two spin components (which reflects the decay of the Néel order), and within the upper Hubbard band (which reflects the change of the kinetic energy of the doublons), while the total weight in the upper band is roughly constant (see, e.g., Ref. [24]). For the analysis in the main text, we take n δ = (W ↑ (t max ) + W ↓ (t max ))/2. For U = 0 the time evolution of the Néel state on the Bethe lattice can be obtained analytically by solving the Heisenberg equations of motion for the c-operators, which provides a good check for the numerical implementation. For completeness, we provide this solution in the following. We choose site 0 to be the origin of the Bethe lattice, which is on the A sub-lattice without loss of generality. One can map the solution of equations of motion on the Bethe lattice to a one-dimensional semi-infinite chain by introducing operators which are invariant under all permutations of the branches of the Bethe lattice [60], where Z n = i:|i−0|=n is the number of sites on the n-th nearest neighbor shell. Then the action of the Hamiltonian is determined by [H, Hence, eigenvectors for the eigenvalue ǫ satisfy the equation and are thus given by the Chebychev polynomials of second kind [61], φ(ǫ) n = U n (ǫ/2) for −1 ≤ ǫ/2 ≤ 1. The U n can be conveniently written as U n (cos(θ)) = sin[(n + 1)θ] sin(θ) , from which one can also see the orthogonality is given by C 0 (t) = The second to last line is a variable transformation θ → π − θ, and in the last line we have used the integral representation of the Bessel function [61], J n (z) = (−i) n π π 0 dθ cos(nθ)e iz cos(θ) .
The explicit form of the C-operators can be used to obtain local observables where the expectation values are simple initial state values. We start by evaluating the time evolution of the magnetic order, n ↑ − n ↓ , at site 0, in the classical Neel state. For the latter we have and hence (where the summation index has been shifted by one). We can now use Gegenbauers addition theorem for Bessel functions [61] to obtain the final result Ψ Néel |n ↑ (t) − n ↓ (t)|Ψ Néel = J 1 (4t) 2t , which fits the numerics.
Next we compute site 0 expectation values on the probe site. Now the initial state is a superposition We evaluate the cross-spin expectation values Spin-flip expectation values are only non-zero in the initial state at site 0, where we have Hence S + 0 (t) is purely real, so that the dynamics is entirely longitudinal in the S x -S y -plane, S y 0 (t) = 0. (C24)