Semiclassical dynamics of spin density waves
Gia-Wei Chern Kipton Barros Zhentao Wang Hidemaro Suwa Cristian D. Batista







Department of Physics, University of Virginia, Charlottesville, VA 22904, USA



Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA



Department of Physics and Astronomy, The University of Tennessee, Knoxville, TN 37996, USA



Department of Physics and Astronomy, The University of Tennessee, Knoxville, TN 37996, USA




Department of Physics and Astronomy, The University of Tennessee, Knoxville, TN 37996, USA
Quantum Condensed Matter Division and Shull-Wollan Center, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA





We present a theoretical framework for equilibrium and nonequilibrium dynamical simulation of quantum states with spin-density-wave (SDW) order. Within a semiclassical adiabatic approximation that retains electron degrees of freedom, we demonstrate that the SDW order parameter obeys a generalized Landau-Lifshitz equation. With the aid of an enhanced kernel polynomial method, our linear-scaling quantum Landau-Lifshitz dynamics (QLLD) method enables dynamical SDW simulations with N ≃ 10^5 lattice sites. Our real-space formulation can be used to compute dynamical responses, such as dynamical structure factor, of complex and even inhomogeneous  SDW configurations at zero or finite temperatures. Applying the QLLD to study the relaxation of a noncoplanar topological SDW under the excitation of a short pulse, we further demonstrate the crucial role of spatial correlations and fluctuations in the SDW dynamics.


Semiclassical  dynamics of spin density waves
    Cristian D. Batista
    December 30, 2023
=============================================



Quantum states with unusual broken symmetries have long fascinated physicists because of their fundamental importance and potential technological applications. Of particular interest is the regular spatial modulation of electron spin
known as spin-density wave (SDW) state <cit.>. SDWs are ubiquitous in strongly correlated systems and play a crucial role in several intriguing many-body phenomena. For example, the SDW state is proximate to the superconducting phase in several unconventional superconductors, including cuprates and iron pnictides. Indeed, non-Fermi liquid behavior is usually observed in the vicinity of a SDW phase transition <cit.>. Moreover, conduction electrons propagating in a noncoplanar spin texture acquire a nontrivial Berry phase and exhibit unusual transport and topological properties <cit.>.  Consequently, metallic SDW with complex spin structures, such as spirals or skyrmion crystals, offers a novel route to control the charge degrees of freedom through manipulation of spins and vice versa <cit.>.  

While analytical techniques have yielded much insight about itinerant magnetism and SDW states <cit.>, numerical methods continue to provide valuable benchmarks and shed light on controversial issues. Among the various numerical tools <cit.>, quantum Monte Carlo (QMC) simulations provide numerically exact solutions to strongly correlated models <cit.>. However, one severe restriction of most QMC methods is the infamous sign-problem. Powerful alternative  approaches that avoid the sign-problem include dynamical mean-field theory (DMFT) <cit.> and density-matrix renormalization group (DMRG) <cit.>.  Significant developments have also been made in their nonequilibrium extension such as time-dependent (TD) DMFT <cit.> and TD-DMRG <cit.>. Both methods, however, are still very limited in their treatment of complex mesoscopic structures.


In this paper, we present a different numerical approach to SDW dynamics, emphasizing the ability to simulate large-scale lattices and complex SDW orders that often occur in highly frustrated systems. Our starting point is a semi-classical treatment of Hubbard-like models, which  neglects quantum fluctuations, but retains the spatial fluctuations of the SDW field. In a way, this approach is the complement of DMFT, which includes quantum fluctuations at the expense of neglecting spatial correlations. A systematic approach is then developed to reintroduce quantum dynamics to the SDW order parameter. We show that in the leading adiabatic approximation, the SDW dynamics is described by a generalized Landau-Lifshitz (LL) equation in which the effective forces acting on the spins are generated from itinerant electrons. Our numerical scheme can be viewed as a quantum LL dynamics (QLLD), in which the electronic degrees of freedom are integrated out at each time step. By supplementing the LL equation with Ginzburg-Landau type relaxation and stochastic terms,  our QLLD method can be used to simulate SDW dynamics both near and far-from equilibrium.



§ SPIN-FERMION HAMILTONIAN FOR EQUILIBRIUM SDW PHASES



We consider the one-band Hubbard model with an on-site repulsion U > 0,

    ℋ = - ∑_ij, α t_ij c^†_i, α c^ _j, α + U ∑_i n_i,↑ n_i ↓.

After performing the Hubbard-Stratonovich (HS) transformation <cit.>, we obtain the following spin-fermion Hamiltonian

    ℋ_ SDW = - ∑_ij, α t_ij c^†_i, α c^ _j, α - 2 U ∑_i 𝐦_i ·𝐬_i + U ∑_i |𝐦_i|^2,

where 𝐬_i = 1/2 c^†_i,ασ_αβ c^ _i, β is the spin operator of conduction electrons and σ is a vector of the Pauli matrices. The local HS or auxiliary field 𝐦_i is a classical O(3) vector in ℝ^3.

Here we have set ħ = 1.  Importantly, since ℋ_ SDW describes non-interacting electrons coupled to a magnetic background, the fermionic degrees of freedom can be integrated out either in Monte Carlo or dynamical simulations to be described below.

This HS Hamiltonian is typically the starting point for determinant QMC (DQMC) simulations <cit.>. In an alternative approach, one assumes static HS variables; then the above Hamiltonian resembles the so-called spin-fermion model <cit.> and can be simulated with Markov chain Monte Carlo assuming classical “spins” 𝐦_i, while electrons are treated quantum mechanically. Applying this method to the cubic-lattice Hubbard model, the obtained Néel temperature agrees remarkably well with those from DQMC simulations <cit.>. It is worth noting that while Eq. (<ref>) is similar to the Hartree-Fock treatment of the Hubbard model, retaining spatial fluctuations of the local HS fields {𝐦_i} in this static (in imaginary time) HS-field formalism goes beyond the usual mean-field method. For instance, this approach captures the critical fluctuations, and consequently the correct universality class, of any continuous thermodynamic transition into a magnetically ordered state.




Instead of Markov-Chain Monte Carlo, here we employ the stochastic Ginzburg-Landau (GL) relaxation dynamics <cit.> to sample the equilibrium SDW configurations within the static HS-field approximation,

    d𝐦_i/dt = - γ∂⟨ℋ_ SDW⟩/∂𝐦_i   + ξ_i(t).

Here γ is a damping constant, and ξ_i is a δ-correlated fluctuating force satisfying ⟨ξ_i(t) ⟩ = 0  and ⟨ξ_i^μ(t) ξ_j^ν(t') ⟩ = 2 γ  k_B T δ_ijδ_μνδ(t - t') for vector components μ and ν. This equation is similar to the over-damped Langevin dynamics used in Ref. <cit.> for the Kondo-lattice model. We note that the magnitude |𝐦_i| of the O(3) vector 𝐦_i is not fixed in this over-damped dynamics. A fictitious inertial mass term can be added to the above dynamics  to improve the efficiency of the simulation. Unlike conventional GL  simulations, where the force is given by the derivative of a phenomenological energy functional, here the force is computed by solving the equilibrium electron liquid of ℋ_ SDW at each time-step. Or equivalently, the effective energy functional ℰ_ eff({𝐦_i}) = ⟨ℋ_ SDW⟩ is obtained by integrating out electrons on the fly. Our method is thus similar to the so-called quantum molecular dynamics (MD) simulations, in which the inter-atomic force is computed by solving the quantum electron Hamiltonian <cit.>, instead of being derived from a phenomenological classical potential. Drawing on this analogy with MD simulations, our approach can be viewed as a quantum GL method for SDW. Interestingly, the quantum MD method in conjunction with the functional integral theory has already been employed to obtain complex magnetic orderings in itinerant magnet compounds in the past <cit.>.


The GL method is particularly powerful when combined with our recently developed kernel polynomial method (KPM) and its gradient transformation in which forces acting on all spins can be efficiently computed <cit.>; see Ref. <cit.> for more details. The resulting linear-scaling KPM-GL method allows us to simulate large lattices with N ≃ 10^5 to 10^6 sites. We apply the KPM-GL simulations to the triangular-lattice Hubbard model as a benchmark. The phase diagram shown in Fig. <ref> agrees very well with those obtained by holon-doublon mean-field <cit.> and the rotational-invariant slave-boson (SB) <cit.> calculations. Note that the SB method <cit.> describes both large and small U regimes, and shows quantitative agreement with QMC  over a wide range of interaction and doping for the square-lattice Hubbard model <cit.>.



At large U, our simulation finds the expected 120^∘-order, which is the ground state of the Heisenberg Hamiltonian arising from the strong-coupling limit of the half-filled Hubbard model. In fact, the semi-classical Hamiltonian Eq. (<ref>) reduces to the classical Heisenberg spin model in the U/t ≫ 1 limit. To see this, we note that any ground state for t_ij=0 contains exactly one electron in each site with its spin locally aligned with the SDW field 𝐦_i. The magnitude of the SDW field freezes at |𝐦_i|=1/2 in this U →∞ limit because amplitude fluctuations have an energy cost proportional to U.  In analogy with the large U limit of the original Hubbard model,  the ground state manifold is massively degenerate because each  𝐦_i can point in an arbitrary direction. 
The degeneracy is removed to second order in t_ij. The low-energy effective Hamiltonian is obtained by considering virtual  electron hopping processes between two neighboring sites i and j.  The Pauli exclusion principle dictates that an electron at site-i can hop to the neighboring site only when its spin is anti-aligned with that of the local moment at site j. Consequently, the effective hopping constant between the two neighboring sites is t^ eff_ij = t_ij⟨χ_i | χ_j ⟩ = t_ijsin(θ_ij/2), where |χ_i⟩ are local electron spinor eigenstate and θ_ij is the angle between the two local moments. At second-order, the energy gain through the virtual electron hopping produces the effective interaction

    ℋ_ij = -2 t_ij^2 sin^2(θ_ij/2)/U = 4 t_ij^2/U(𝐦_i ·𝐦_j - 1/4).

This result corresponds to the classical limit of the S=1/2 Heisenberg model, implying that  the semi-classical SDW Hamiltonian (<ref>) correctly captures the classical limit of the half-filled Hubbard model in the strong-coupling regime <cit.>.

As U is decreased, our simulation shows that the 120^∘ order is replaced by an interesting commensurate collinear SDW as the ground state <cit.>. The collinear SDW at this intermediate U is still gapped electronically and exhibits a zigzag structure. By computing the electron density of states (DOS), we find a metal-insulator transition at U_c2 between an incommensurate spiral and a commensurate collinear SDW phase (see Fig. <ref>). Finally, the metallic spiral SDW undergoes a continuous transition at into a paramagnetic state at U_c1.
We note in passing that the metal-insulator transition, obtained with other numerical techniques (e.g. path-integral renormalization group method <cit.>), is entirely within the paramagnetic regime <cit.>, implying the existence of a paramagnetic insulator (or spin liquid) at intermediate U values. This state  cannot be obtained with our semi-classical approach because it is stabilized by strong fluctuations of the HS fields along the imaginary time axis. However, the existence of this phase remains to be settled. Recent variational QMC <cit.> and DMFT <cit.> calculations show that spiral SDW is more favorable than the spin liquid phase. Nonetheless, if a magnetic phase is stabilized by, e.g. applying a magnetic field, we expect the incommensurate and collinear SDWs obtained here to be strong candidates.




§ SEMICLASSICAL DYNAMICS AND LANDAU-LIFSHITZ EQUATION



Having demonstrated that the semiclassical Hamiltonian ℋ_ SDW provides a viable approach to equilibrium SDW phases, a natural question is whether we can use it to study the  SDW dynamics. To this end, we need to reintroduce physical dynamics to the “static” auxiliary SDW field. We first note that the spin-fermion Hamiltonian ℋ_ SDW can also be obtained from a Hartree-Fock decoupling of the interaction term that varies from one site to another. The important difference relative to the HS approach is that the SDW field satisfies the self-consistent condition 𝐦_i = ⟨ s_i ⟩ =  Tr(ρ 𝐬_i), where ρ is the density matrix characterizing the physical electron state. Consequently, the field 𝐦_i belongs to the sphere of radius 1/2  (0 ≤ |𝐦_i| ≤ 1/2). 





To derive the time dependence of the SDW field, we start with the continuity equation associated with the total spin conservation: d 𝐬_i /dt = - ∑_j𝐉_ij, where 𝐬_i = 1/2 c^†_iασ_αβ c^ _iβ is the electron spin, and 𝐉_ij = -it_ij/2σ_αβ (c^†_iα c^ _jβ - c^†_jαc^ _iβ) is the spin current density on bond ⟨ ij ⟩. We then introduce the single-particle density matrix ρ with elements ρ_iα, jβ≡⟨ c^†_jβ c^ _iα⟩. Taking the average of the continuity equation leads to 

    d𝐦_i/dt = -i/2∑_j t_ijσ^ _βα( ρ_iα, jβ - ρ_jα, iβ).

The SDW dynamics is thus related to the time evolution of the density matrix ρ, which obeys the von Neumann equation dρ/dt = i [ρ, H^ eff]. Up to a constant, H^ eff is the effective single-electron Hamiltonian defined as ℋ_ SDW = ∑_iα, jβ H^ eff_iα, jβ c^†_iα c^ _j,β. Using Eq. (<ref>), we obtain the equation of motion  

    dρ_iα, jβ/dt    =     i (t_ik ρ_kα, jβ - ρ_iα, kβ  t_kj)   
       +     i U  (𝐦_i ·σ_αγ ρ_iγ, jβ 
    	 - ρ_iα, jγ σ_γβ·𝐦_j ).

The electron density matrix is partly driven by the time-varying SDW field. Eqs. (<ref>) and (<ref>) comprise  a complete set of coupled ordinary differential equations for the SDW dynamics. An alternative is to dispense of Eq. (<ref>) and substitute 𝐦_i(t) =  ρ_iα, iβ(t) σ_βα in Eq. (<ref>), giving rise to a set of nonlinear differential equations for ρ. 

In general, time dependence of physical quantities in mean-field approaches can be obtained using the Dirac-Frenkel variational principle <cit.>. Our derivation here, on the other hand, is based on the spin-density continuity equation, hence emphasizing the importance of conservation laws in physical dynamics. This physically intuitive approach can be easily generalized to obtain dynamics for other symmetry-breaking phases.

Our approach here is a real-space formulation of the time-dependent Hartree-Fock (TDHF) method <cit.>, similar to the familiar time-dependent Bogoliubov-de Gennes equation for superconductors or Bose condensates <cit.>. Assuming that the order parameters are characterized by well defined momenta, e.g. 𝐦_i = ∑_r 𝐌_r  exp(i 𝐐_r ·𝐫_i), the above equations can be simplified due to the translation invariance. The problem is then reduced to a set of coupled differential equations for density-matrix elements n_αβ(𝐤, t) ≡⟨ c^†_𝐤α c^ _𝐤β⟩ and g^r_αβ(𝐤, t) ≡ (𝐌^ _r ·σ_αβ) ⟨ c^†_𝐤α   c^ _𝐤+𝐐_r, β ⟩ in momentum space. This mean-field approximation of TDHF has recently been applied to the out-of-equilibrium dynamics of BCS-type superconductors <cit.>, and of Néel-type SDW <cit.>. 

Since the order-parameter field 𝐌_r is assumed to be uniform, spatial inhomogeneity and/or fluctuations are ignored in such k-space approaches. Our formulation here does not require the prerequisite knowledge of ordering patterns, and are particularly capable of simulating complex symmetry-breaking phases, inhomogeneous configurations, and disordered phases with  preformed local moments, such as the paramagnetic state in the large U limit.


The SDW dynamics Eq. (<ref>) can be simplified in the large U limit in the so-called adiabatic approximation, which assumes that electrons quickly relax to the ground state of the instantaneous SDW configuration {𝐦_i}. Using second-order perturbation theory, one readily computes the average spin current: ⟨𝐉_ij⟩ = 4 t^2_ij/U 𝐦_i ×𝐦_j. Substituting this into the right-hand side of Eq. (<ref>) gives rise to a Landau-Lifshitz (LL) equation 

    d𝐦_i/dt = -  ∑_j 4 t_ij^2/U 𝐦_i ×𝐦_j,

with an effective torque computed using the Heisenberg exchange of Eq. (<ref>). More details can be found in Ref. <cit.>.

For intermediate and small U/t values, one needs to solve the von Neumann equation. Since the number of independent density-matrix elements is of order 𝒪(N^2) for a lattice of N spins, the computational cost of integrating the von Neumann equation is tremendous for large lattices, e.g. N∼ 10^5. To further simplify the numerical calculation, here we derive the SDW dynamics in a similar adiabatic limit for arbitrary U. Formally, we employ the multiple-time-scale method <cit.> and introduce an adiabaticity parameter ϵ∼ |d𝐦/dt| such that the fast (electronic) and slow (SDW) times are τ = ϵ t and t, respectively. The single-particle Hamiltonian varies with the slow time, i.e. H^ eff({𝐦_i }) = H^ eff(τ). Expanding the density matrix in terms of the adiabaticity parameter: ρ(t) = ρ^(0)(τ) + ϵρ^(1)(t, τ) + ϵ^2 ρ^(2)(t, τ) + ⋯, and plugging it into the von Neumann equation, we obtain [ρ^(0), H^ eff] = 0 and dρ^(ℓ)/dt - i [ρ^(ℓ), H^ eff] = -dρ^(ℓ-1)/dτ for ℓ≥ 1. This provides a systematic approach to obtain the time dependence of the density matrix.


Here we use the leading adiabatic solution ρ^(0) to compute the expectation value of the spin-current density 𝐉_ij, which is the right-hand side of Eq. (<ref>). We first write H^ eff = T + Σ where T_iα, j β = - t_ijδ_αβ is the tight-binding Hamiltonian and Σ_iα, j β = -U δ_ij 𝐦_i ·σ_αβ is the spin-fermion coupling. It is then straightforward to show that Eq. (<ref>) is simply d 𝐦_i/dt = i σ_αβ [ρ, T]_i β, iα/2. Using the adiabatic equation [ρ^(0), H^ eff] = 0, we have [ρ^(0), T] = -[ρ^(0), Σ], which gives  

    d𝐦_i/dt = -i U /4σ_αβ[ (σ^ _βγρ^(0)_iγ, iα - ρ^(0)_iβ,iγσ^ _γα)
    	·𝐦_i ].

The right-hand side of the above equation can be further simplified using the properties of Pauli matrix multiplication: σ^a σ^b = δ_abℐ_2× 2 + i ϵ_abcσ^c, where a, b, c are x, y, z. For example,  σ_αβ(σ_βγ·𝐦_i) ρ^(0)_iγ,iα = n^(0)_i 𝐦_i + i 𝐦_i ×σ^ _αβρ^(0)_iβ, iα, where n^(0)_i = ρ^(0)_iα, iα is the local electron density. After some algebra, we obtain

    d𝐦_i/dt = U/2𝐦_i ×σ^ _βα ρ^(0)_iα, iβ
    	= - 𝐦_i ×∂⟨ℋ_ SDW⟩/∂𝐦_i .

The second equality comes from the fact that ρ^(0) is computed from the equilibrium electron liquid described by ℋ_ SDW. The local  electron spin ⟨𝐬_i ⟩ = 1/2⟨ c^†_iασ^ _αβ c^ _iβ⟩ = 1/2σ^ _βα ρ^(0)_iα, iβ acts as an effective magnetic field and drives the slow dynamics of the SDW field.
Importantly, this equation shows that the adiabatic SDW dynamics is described by the Landau-Lifshitz (LL) equation <cit.> with an effective energy functional ℰ_ eff({𝐦_i }) = ⟨ℋ_ SDW⟩, obtained from the equilibrium electronic state of the instantaneous spin-fermion Hamiltonian.



 §.§ Benchmark with exact diagonalization


We first benchmark our semiclassical SDW dynamics, with and without the adiabatic approximation, against the exact diagonalization (ED) calculation of the original Hubbard model. To this end, we apply our formulation to the two-sublattice collinear Néel state that is obtained for the half-filled Hubbard model on a square lattice. Since we only include NN hopping, the Néel ordering is stable for any positive value of U/|t|.  Specifically, as shown in Fig. <ref>, we compute the dynamical structure factor 𝒮(𝐤, ω) of a 4 × 4 Hubbard cluster with periodic boundaries for  U/t=7.33. Details of the ED calculation are described in Ref. <cit.>.
We compare the ED result at T=0 and the semiclassical SDW dynamics at an extremely low temperature 
(classical moments freeze at T=0 <cit.>). We  set temperature at T=10^-4t and verify that  the results do not change upon decreasing the temperature to T=10^-5, indicating that our results  capture the dynamical response of the classical moments in the
T → 0 limit.

For both the real-space TDHF dynamics [Eqs. (<ref>) and (<ref>)] and the adiabatic dynamics [Eq. (<ref>)], SDW states are first generated by means of GL-Langevin simulations described in Sec. <ref>. The obtained spin configurations, which are representative of the canonical ensemble, are used as the initial condition for dynamical simulations. The dynamical structure factor, 𝒮(𝐤, ω), is calculated by applying the space-time Fourier transform to the time evolution of the auxiliary field 𝐦_i(t). 

In the SDW dynamics,  the elastic peak has a finite width for finite duration of the dynamical simulation (not shown in Fig. <ref>(b) and (c)). The area under the elastic peak is proportional to N ⟨𝐦_i ⟩^2, where N is the number of sites. In contrast, the lowest energy peak of the exact result appears at a small but finite frequency arising from quantum fluctuations neglected by the semiclassical treatment (the exact ground state is a singlet state for a finite size system). This quasi-elastic peak becomes the elastic peak of the spontaneously broken symmetry state in the thermodynamic limit. In Fig. <ref>, we normalize the spectral weights obtained from the SDW dynamics so that the total weight of inelastic peaks equals that obtained from the ED excluding the quasi-elastic peak.





 The low-energy spectrum of the original Hubbard model is not well described by a simple effective spin Hamiltonian for U/|t|=7.33 (charge fluctuations can strongly renormalize the spin-wave dispersion). It is then quite remarkable that all the approaches produce a rather flat magnon dispersion for the wave-vectors included in a 4× 4 square lattice. However, the excitation energies in the adiabatic approximation are roughly 25% lower than the exact result [see Fig. <ref>]. This discrepancy is attributed to two factors: the semi-classical treatment of the spin degrees of freedom and the adiabatic approximation. The excitation energies obtained from the real-space TDHF method are  approximately 15% lower than the exact result. We thus conclude that the adiabatic approximation accounts for roughly 10% of the discrepancy, while the semiclassical  treatment accounts for the remaining 15%. In addition, the normalized spectral weights (areas) in both the semiclassical dynamics are different only approximately 30% from the exact result. A much better quantitative agreement is expected for 3D systems, but their solutions are beyond the scope of state of the art ED methods.

The real-space TDHF method captures not only the transverse modes, but also the longitudinal mode arising from charge fluctuations <cit.>. We note however that the method does not capture the longitudinal  spin fluctuations associated with quantum fluctuations of the magnetic moments. Unlike charge fluctuations, quantum (magnetic) fluctuations persist for arbitrarily large-U/t (they arise from fluctuations of the m-field along the imaginary time axis). These longitudinal fluctuations  correspond to  two-magnon excitations in a 1/S expansion <cit.>. Correspondingly, they have an energy of order J ∝ t^2/U for large U/t. In contrast, the longitudinal spin fluctuations arising from charge fluctuations lead  to the high-energy peaks at ω∼ U, which are observed in the real-space TDHF dynamics, as shown in the inset of Fig. <ref>(c), while they are absent in the adiabatic dynamics. Nevertheless, as expected for this value of U/t, the longitudinal mode is well separated from the transverse modes. 




 §.§ 120^∘ SDW order in Hubbard and Anderson-Hubbard model


Our benchmark study shows that both the TDHF and the adiabatic approach provide a reasonable description of the SDW dynamics. It is worth noting that the linearized TDHF equation of motion corresponds to the random phase approximation (RPA) <cit.>. Our real-space formulation of the TDHF thus provides an efficient and universal numerical approach to describe nonlinear dynamics  beyond the RPA level. Moreover, our approach allows for computation of dynamical response functions  at any finite temperature, including the high-temperature regime in which the magnetic moments only exhibit short range correlations. 




Another unique feature offered by our real-space method is the capability of computing the dynamical response of inhomogeneous SDW. To demonstrate this, here we apply the above adiabatic dynamics Eq. (<ref>) to compute 𝒮(𝐤, ω) for  the 120^∘ SDW  depicted in Fig. <ref>(c) for Hubbard model with quenched disorder.  Specifically, we consider the Anderson-Hubbard model by adding an on-site potential disorder ∑_i,α V^ _i  c^†_iα c^ _iα to Eq. (<ref>). This model has served as a canonical platform for investigating the intriguing interplay of localization and correlations. Relevant to our study here is the effect of disorder on long-range SDW order. For Néel-type SDW on a half-filled bipartite lattice, it has been shown that increasing the disorder first closes the electron spectral gap, while the SDW remains finite <cit.>. The disappearance of the SDW order parameter occurs at a larger disorder <cit.>. This result is relevant to the non-equilibrium dynamics to be discussed below. 

We first compute the dynamical structure factor for the pure Hubbard model on a large lattice with N = 120 × 120 sites. The results shown in Fig. <ref>(a) resemble the linear spin wave dispersion of the Heisenberg  model  <cit.>, except for a significantly renormalized lower-energy branch. 
Next we include a Gaussian disorder with zero mean and  standard deviation σ_V = 0.45  t_ nn, which is relatively large yet not strong enough to destroy the SDW order <cit.>.
The 𝒮(𝐤, ω) computed using the adiabatic SDW dynamics is shown in Fig. <ref>(b). While the overall dispersion is similar to that of the SDW in the disorder-free Hubbard model [Fig. <ref>(a)], there are a few notable new features. Firstly, the magnon dispersion is significantly broadened by the quenched disorder, and the middle of the low-energy branch is further renormalized. Interestingly, a rather sharp dispersion remains near the zone center, indicating that these long-wavelength modes are less sensitive to disorder. Secondly, several new modes appear at low energies, especially below the original gap at the M point. Interestingly, similar disorder-induced low-energy modes are also obtained in bi-layer Heisenberg antiferromagnet using the bond-operator method <cit.>. 
A systematic study of the SDW dynamics with disorder will be left for future studies. 



§ NONEQUILIBRIUM DYNAMICS AT FINITE TEMPERATURES



The adiabatic LL equation can also be used to study nonequilibrium SDW phenomena as long as the electron relaxation is much faster than the SDW dynamics. Here we first generalize the adiabatic dynamics to finite temperatures by adding dissipation and fluctuations to Eq. (<ref>). It is worth noting that the adiabatic SDW dynamics preserves the length of local moments |𝐦_i|. Longitudinal spin relaxation and fluctuations thus come from either higher order terms in the adiabatic expansion or other processes beyond the self-consistent field approach. The standard Gilbert damping also preserves the spin length <cit.>. Instead, here we combine the Ginzburg-Landau relaxation discussed in Eq. (<ref>) with the adiabatic dynamics  of Eq. (<ref>) to account for the longitudinal relaxation <cit.>. This procedure gives rise to the following generalized LL equation

    d𝐦_i/dt = - 𝐦_i ×∂⟨ℋ_ SDW⟩/∂𝐦_i  
    	- γ∂⟨ℋ_ SDW⟩/∂𝐦_i   + ξ_i(t).

 ℋ_ SDW is the spin-fermion Hamiltonian defined in Eq. (<ref>), γ is a damping constant, and ξ_i is a δ-correlated fluctuating force satisfying ⟨ξ_i(t) ⟩ = 0 and ⟨ξ^μ_i(t) ξ^ν_j(t') ⟩ = 2 γ  k_B T δ_ijδ_μνδ(t - t'). The damping coefficient and the stochastic terms are chosen such that the dissipation-fluctuation theorem is satisfied and the above LL equation can be used to faithfully sample the equilibrium Boltzmann distribution at finite temperatures <cit.>. 

A microscopic calculation of the damping coefficient γ is beyond the adiabatic approximation. In the real-space TDHF method, relaxation of SDW mainly arises from the Landau damping mechanism, which describes the energy transfer from the collective SDW mode to single-particle excitations <cit.>. Electron-electron scattering, which is not captured by the TDHF, also contributes to the damping of SDW, especially in ultrafast dynamics of metals <cit.>. Moreover, for open systems as in most pump-probe experiments, coupling of electrons to other degrees of freedom <cit.>, such as phonons, also play an important role in the relaxation of SDW dynamics. Here γ is treated as as a phenomenological parameter which we chose to ensure the adiabatic approximation.  

That the SDW field obeys the LL dynamics can be understood intuitively from the fact that the Heisenberg equation of motion for spin operators corresponds to the classical LL equation <cit.>.  Here we give a microscopic derivation  starting from the Hubbard model, which reveals the condition for the validity of the LL dynamics (adiabatic approximation of the von Neumann equation). In fact,  the adiabatic approximation has been widely employed for spin dynamics in the context of time-dependent spin-density-functional theory <cit.>. Our results thus provide a theoretical foundation for the LL dynamics of SDW, and pave the way for systematic improvements beyond the adiabatic approximation. 

It is worth noting that, in contrast to the conventional LL method, the energy functional in our approach is obtained by solving the spin-fermion Hamiltonian ℋ_ SDW at each time-step. In analogy with the quantum MD simulations <cit.>, our numerical scheme can then be viewed as a quantum LL dynamics (QLLD) method. Although solving the electron Hamiltonian on the fly is computationally expensive, large-scale (N ∼ 10^5) QLLD simulations are enabled by our recently developed KPM algorithm with automatic differentiation, such that the “forces” can be computed along with the total energy without extra overhead <cit.>.



We next apply our stochastic QLLD method to investigate the time evolution of a topological SDW on the triangular lattice that arises as a weak-coupling instability at filling fraction n = 3/4 <cit.>. The combination of a van Hove singularity and perfect Fermi surface nesting at this filling fraction gives rise to a magnetic susceptibility that diverges as χ(𝐪) ∝log^2| 𝐪 - 𝐐_η|, where 𝐐_η (η = 1,2,3) are the three nesting wavevectors <cit.>. The system thus tends to develop a triple-𝐐 SDW characterized by three vector order parameters: 𝐦_i = Δ_1 e^i𝐐_1 ·𝐫_i + Δ_2 e^i𝐐_2 ·𝐫_i + Δ_3 e^i𝐐_3 ·𝐫_i. Note that the phase factors e^i 𝐐_η·𝐫_i = ± 1. In general, there are four distinct local moments, leading to a quadrupled magnetic unit cell.
The SDW instability of triangular-lattice Hubbard model at n=3/4 filling is similar to that of the half-filled Hubbard model on square lattice. However, unlike the simple Néel order in the later case, there are several possible triple-𝐐 SDWs <cit.>. 

At the lowest temperatures, the magnetic ordering consists of a non-coplanar SDW with |Δ_1|=|Δ_2|=|Δ_3| and Δ_1 ⊥Δ_2 ⊥Δ_3 <cit.>. This SDW is also called a tetrahedral or all-out order as spins in the unit cell point to the four corners of a regular tetrahedron <cit.>; see Fig. <ref>(a). Moreover, as the spins on each triangular plaquette are non-coplanar, the resulting nonzero scalar spin chirality 𝐦_i ·𝐦_j ×𝐦_k = ± 4Δ^3 also breaks the parity symmetry. Consequently, the tetrahedral SDW is also characterized by a discrete Z_2 chirality order parameter. More importantly, electrons propagating in this non-coplanar SDW acquire a nonzero Berry phase, which is equivalent to a uniform magnetic field. Since the Fermi surface is gapped out by the SDW, the resulting electron state exhibits a spontaneous quantum Hall effect with transverse conductance σ_xy = ± e^2/h <cit.>.



Motivated by  a recent pump-probe experiment on the ultrafast SDW dynamics in chromium <cit.>, we perform simulations of this topological SDW subject to a short heat pulse. For simplicity, we assume that the effect of the pump pulse is to inject energy to the electron system, which quickly equilibrates to a state characterized by temperature T_e. This is consistent with our adiabatic approximation for the SDW dynamics. The time dependence of the effective electron temperature is governed by the rate equation C dT_e/dt = -G (T_e - T_L) + Q(t) <cit.>, where C is the heat-capacity of the electron liquid, G is the coupling to the lattice, T_L is the lattice temperature, and Q(t) ∝exp[-(t - t_p)^2 / w^2] is the heat source due to the pump pulse. We further assume that T_L ≈ 0 throughout the relaxation process. The resultant T_e(t) curve is shown in Fig. <ref>(b).

T_e(t)  is then used for the stochastic noise ξ(t) in our QLLD simulations of Eq. (<ref>). We use the parameters, U = 3, damping γ = 0.1, G/C = 0.02, t_p = 15, and w = 5, in units of the NN hopping t_ nn. The lattice size is N=120^2. Fig. <ref>(c) shows the evolution of the magnetic order parameter at the nesting wavevectors: ℳ = √(|Δ_1|^2 + |Δ_2|^2 + |Δ_3|^2) normalized to its maximum. We also estimate the time dependence of the electron spectral gap ε_ gap from the instantaneous DOS [see Fig. <ref>(c)]. Interestingly, as the temperature rises, the decline of ℳ is rather slow compared with the closing of the energy gap. In fact, the SDW order parameters Δ_η remain finite throughout the process, while the gap closes quickly after the photoexcitation (at t ≈ 12). In equilibrium the SDW order parameters disappear along with the gap above the transition temperature <cit.>, implying that the  photoexcited SDW is in a highly non-equilibrium transient state. As the system relaxes, the gap reopens at a later time [see Fig. <ref>(c)]. 

This picture is further supported by our calculation of instantaneous longitudinal and transverse  conductivities  shown in Fig. <ref>(d). Here we use KPM to compute the Kubo-Bastin formula for  the conductivities <cit.>. The error bars  are estimated from five independent simulations. The electrons exhibit a negligible longitudinal conductivity σ_xx≈ 0 and a quantized Hall conductivity σ_xy = e^2/h in the gapped regimes, as expected for this topological SDW. On the other hand, the longitudinal conductivity increases significantly during the period of vanishing gap, while the transverse conductance decreases  and exhibits small oscillations in the vicinity of the gap-closing transitions.

The closing and subsequent re-opening of the SDW gap  have been reported in recent pump-probe experiment on chromium <cit.>. The ultrafast SDW dynamics seem to be well described by a model that assumes  a thermalized electron gas. However, the closing of the gap is assumed to be always accompanied by the disappearance of the SDW order parameter in Ref. nicholson16, which is not necessarily the case.  As demonstrated in our simulations, an out-of-equilibrium electron state might be gapless while the spin density remains modulated. Indeed, similar pump-probe experiments on the charge density wave (CDW) have revealed a fast collapsing of electronic gap in the time-resolved photoemission spectroscopy <cit.>, and a reduced, yet finite, modulation of charge density inferred from core-level X-ray photoemission <cit.> during the nonequilibrium melting process. Numerical simulations taking into account coupling to the lattice distortion showed that the CDW order parameter can indeed be partially decoupled from the spectral gap dynamics <cit.>. However, it should be noted that the lattice degrees of freedom  introduce a new time scale, in addition to that of the hot electron relaxation. The transient metallic SDW observed in our simulations is probably due to a different mechanism.


To understand the origin of this nonequilibrium metallic SDW, we first note that the electronic gap of this topological SDW arises from the scalar spin chirality <cit.>. Indeed, the electronic gap vanishes for collinear or coplanar triple-𝐐 SDWs <cit.>. This observation leads us to investigate the temporal and spatial fluctuations of the scalar chirality. To this end, we introduce the normalized scalar spin chirality: χ_ =  χ_ijk = (𝐦_i  ·𝐦_j  ×𝐦_k) /  | 𝐦 |^3 for individual triangular plaquettes (here the overline indicates average over all triangles). The time dependence of the (spatial) average and the standard deviation of the scalar chirality χ_ are shown in Fig. <ref>(c). Interestingly, the average chirality remains finite and of the same sign, indicating that the chiral symmetry is still broken in this transient SDW. On the other hand, as shown in Fig. <ref>(c), the standard deviation σ_χ_ increases significantly with T.  In fact, the transient gapless regime coincides roughly with the period when σ_χ_ > χ_, implying that the vanishing gap is due to thermally induced spatial fluctuations of χ_ijk.





Our scenario is confirmed by the spatial distribution of the normalized plaquette chirality at the initial stage (t = 5) and the gapless regime (t  = 40), shown in  Fig. <ref>(a) and (b), respectively. While the chirality is relatively uniform initially (χ_∼ 1), noticeable inhomogeneity develops at later times [see Fig. <ref>(b)]. Histograms of the plaquette chirality h( χ_) and the corresponding electron DOS at various simulation times are shown in Fig. <ref>(c) and (d), respectively. The chirality distribution becomes asymmetric and very broad during the period of vanishing gap. This transient gapless SDW is similar to the disorder-induced metallic antiferromagnetic state observed in the Anderson-Hubbard model <cit.>. As discussed above, the plaquette scalar chirality acts as local magnetic field and it is known that strongly disordered magnetic flux destroys the quantum Hall effect, in agreement with our simulations. Our results thus underscore the importance of thermal fluctuations and spatial inhomogeneity for the nonequilibrium dynamics of SDW, which have been overlooked in most dynamical studies of correlated systems.






§ SUMMARY AND OUTLOOK


We have developed a new theoretical framework for the semiclassical dynamics of SDW in Hubbard-like models. Based on a real-space time-dependent Hartree-Fock (TDHF) method applied to symmetry-breaking phases, our approach provides a Hamiltonian formulation for the SDW dynamics. The time evolution of the SDW field is coupled to the von Neumann equation that describes the dynamics of single-electron density matrix. The formulation correctly reduces to the Holstein-Primarkoff dynamics of magnons (linear spin waves) in the large-U limit at half-filling. We further show that an adiabatic approximation of the von Neumann equation gives rise to a quantum Landau-Lifshitz dynamics (QLLD) for the SDW order parameter. Importantly, the energy functional of the LL equation is computed from an effective spin-fermion Hamiltonian that is obtained from a Hubbard-Stratonovich transformation of the original Hubbard model.  







Our benchmark study of the Néel order on a half-filled Hubbard cluster showed that the semiclassical SDW dynamics agrees reasonably well with the exact diagonalization calculation. We apply our QLLD simulations to compute the dynamical structure factor of a 120^∘ SDW at intermediate values of U/t on the triangular lattice. While the overall spectrum resembles that obtained using linear spin-wave theory for the large U/t limit of the Hubbard model (S=1/2 Heisenberg model),  charge fluctuations produce a significant renormalization of the low-energy branch. We note that quantum fluctuations, not included in our approach,  can also produce a significant renormalization of the spin-wave spectrum of frustrated 2D models <cit.>, like the one considered here. However, renormalization due to quantum fluctuations is much smaller in 3D models, whose dynamical structure factor is typically well described by semiclassical approaches.
Importantly, our real-space approach allows us to include the effects of spatial inhomogeneities of the SDW on large lattices. 

We have demonstrated this unique capability by computing the dynamical structure factor of the same 120^∘ SDW on an Anderson-Hubbard model with disordered on-site potentials. Other than significant broadening of the magnon dispersion, our result shows that the disorder induces many low-energy modes, especially at the boundary of the Brillouin zone. 

Another important application of our QLLD method is the study of SDW-related non-equilibrium phenomena. Here we generalize the LL dynamics by including a Langevin-type damping and the corresponding stochastic noise to account for longitudinal relaxation and fluctuation. We then apply the generalized QLLD scheme to study the evolution of a topological SDW subject to a heat pulse, similar to the situation in the pump-probe setup. Our simulation shows an intriguing transient non-equilibrium SDW on the triangular lattice. While the SDW order parameter decreases with rising electron temperature, it remains finite even when the electronic gap is closed. The gap reopens at a later time as the system relaxes. Since the electronic gap in this topological SDW originates from the noncoplanar spin configuration, we show that the vanishing gap is due to strong spatial fluctuations of the scalar spin chirality, a quantity measuring the non-coplanarity of plaquette spins. Our real-space QLLD simulations thus underscores again the importance of spatial inhomogeneity and thermal fluctuation of the SDW dynamics.

The theoretical framework and numerical method developed in this paper can be easily generalized to study the dynamics of other symmetry-breaking phases, notably charge-density wave and superconductivity. Compared with other phenomenological method (e.g. time-dependent Ginzburg-Landau simulation for superconductors), keeping the electron degrees of freedom allows us to also look into the instantaneous electronic structure during the evolution of the order-parameter field. 

Our efficient semiclassical approach can also be feasibly integrated with first-principles method such as density functional theory (DFT). Here we note the analogy with molecular dynamics (MD) simulations. While classical MD simulations use phenomenological inter-atomic potentials, the quantum MD method computes the forces by solving the electron Hamiltonian on the fly. The quantum MD methods have proven a powerful tool in many branches of physical sciences. Our method can be viewed as the quantum version of the Landau-Lifshitz dynamics. We envision that our QLLD method combined with DFT calculation will provide a novel new approach to SDW dynamics in realistic materials.

It is worth pointing out that the numerical method presented here is complementary to DMFT. Both approaches are not restricted by the sign-problem that plagues the QMC methods. Conventional DMFT ignores spatial correlations from the outset and focuses on quantum effects or fluctuations along the imaginary time axis. Our semiclassical approach emphasizes the large-scale simulations in order to fully take into account the spatial  correlations and fluctuations of the magnetic order parameter. In developing this method, we are partly motivated by several recent studies emphasizing the important role of emergent nano-scale structures in the functionality of strongly correlated materials <cit.>.  Taking advantage of recent developments of efficient electronic structure method, such as KPM, our scheme is to progressively include the quantum corrections at each time step of the dynamical simulations. Our work here has laid the groundwork for systematic improvement beyond the adiabatic or TDHF approximation, which will be left for future studies.




The authors would like to thank A. Chubukov for useful discussions regarding analytical approaches to SDW dynamics. Work at LANL (K. B.) was carried out under the auspices of the U.S. DOE Contract No. DE-AC52-06NA25396 through the LDRD program. Part of the software and algorithm development  (G.W. Chern and C. D.  Batista) was supported by the Center for Materials Theory   as a part of the Computational Materials Science  (CMS) program, funded by the  U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division.



99

fawcett88 E. Fawcett, Spin-density-wave antiferromagnetism in chromium, Rev. Mod. Phys. 60, 209 (1988).

gruner94 G. Grüner, The dynamics of spin-density waves, Rev. Mod. Phys. 66, 1 (1994).

ritz13 R. Ritz, M. Halder, M. Wagner, C. Franz, A. Bauer, C. Pfleiderer, Formation of a topological non-Fermi liquid in MnSi, Nature 497, 231 (2013).

nagaosa10 N. Nagaosa, J. Sinova, S. Onoda, A. H. MacDonald, and N. P. Ong, Anomalous Hall effect, Rev. Mod. Phys. 82, 1539 (2010).

roessler06 U.K. Roessler, A.N. Bogdanov, C. Pfleiderer, Spontaneous skyrmion ground states in magnetic metals, Nature 442, 797 (2006).

nagaosa13 N. Nagaosa and Y. Tokura, Topological properties and dynamics of magnetic skyrmions, Nature Nanotech. 8, 899 (2013).

moryia85 T. Moriya, Spin Fluctuations in Itinerant Electron Magnetism, (Springer, Berlin 1985) and references therein.

brando16 M. Brando, D. Belitz, F. M. Grosche, T. R. Kirkpatrick, Metallic quantum ferromagnets, Rev. Mod. Phys. 88, 025006 (2016).

simon_collab J. P. F. LeBlanc, et al. (Simons Collaboration on the Many-Electron Problem), Solutions of the two-dimensional Hubbard model: benchmarks and results from a wide range of numerical algorithms, Phys. Rev. X 5, 041041 (2015).

blankenbecler81 R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Monte Carlo calculations of coupled boson-fermion systems. I, Phys. Rev. D 24, 2278 (1981).

hirsch85 J. E. Hirsch, Two-dimensional Hubbard model: Numerical simulation study, Phys. Rev. B 31, 4403 (1985).

white89 S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E. Gubernatis, and R. T. Scalettar, Numerical study of the two-dimensional Hubbard model, Phys. Rev. B 40, 506 (1989).

georges96 A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).

kotliar06 G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Electronic structure calculations with dynamical mean-field theory, Rev. Mod. Phys. 78, 865 (2006).

white92 S. R. White,  Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).

schollwoeck04 U. Schollwoeck,  The density-matrix renormalization group, Rev. Mod. Phys. 77, 259 (2005).



freericks06 J. K. Freericks, V. M. Turkowski, and V. Zlatić,  Nonequilibrium Dynamical Mean-Field Theory, Phys. Rev. Lett. 97, 266408 (2006).

aoki14 H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner,  Nonequilibrium dynamical mean-field theory and its applications, Rev. Mod. Phys. 86, 779 (2014).



white04 S. R. White and A. E. Feiguin,  Real-Time Evolution Using the Density Matrix Renormalization Group, Phys. Rev. Lett. 93, 076401 (2004).

daley04 A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal,  Time-dependent density-matrix renormalization-group using adaptive effective Hilbert spaces, J. Stat. Mech. P04005 (2004).











schulz90 H. J. Schulz, Effective action for strongly correlated fermions from functional integrals, Phys. Rev. Lett. 65, 2462 (1990).

mukherjee14 A. Mukherjee, N. D. Patel, S. Dong, S. Johnston, A. Moreo, and E. Dagotto, Testing the Monte Carlo mean field approximation in the one-band Hubbard model, Phys. Rev. B 90, 205133 (2014).

buhler00 C. Buhler, S. Yunoki, and A. Moreo, Magnetic domains and stripes in a spin-fermion model for cuprates, Phys. Rev. Lett. 84, 2690 (2000).

alvarez02 G. Alvarez, M. Mayr, and E. Dagotto, Phase diagram of a model for diluted magnetic semiconductors beyond mean-field approximations, Phys. Rev. Lett. 89, 277202 (2002).

hohenberg77 P. C. Hohenberg and B. I. Halperin, Theory of dynamic critical phenomena, Rev. Mod. Phys. 49 435 (1977).

bray94 A. J. Bray, Theory of phase ordering kinetics, Adv. Phys. 43, 357 (1994).

barros13 K. Barros and Y. Kato, Efficient Langevin simulation of coupled classical fields and fermions, Phys. Rev. B 88, 235101 (2013).


marx00 D. Marx and J. Hutter, Ab initio molecular dynamics: theory and implementation, in Modern Methods and Algorithms of Quantum Chemistry, 2nd ed., edited by J. Grotendorst (John von Neumann Institute for Computing Julich, Germany, 2000).


kakehashi98 Y. Kakehashi, S. Akbar, and N. Kimura, Molecular-dynamics approach to itinerant magnetism with complex magnetic structures, Phys. Rev. B 57, 8354 (1998).

supp1 Details of KPM and its gradient transformation are discussed in the supplemental material.

krishnamurthy90 H. R. Krishnamurthy, C. Jayaprakash, S. Sarker, and W. Wenzel, Mott-Hubbard metal-insulator transition in nonbipartite lattices, Phys. Rev. Lett. 64, 950 (1990).

capone01 M. Capone, L. Capriotti, F. Becca, S. Caprara, Mott metal-insulator transition in the half-filled Hubbard model on the triangular lattice, Phys. Rev. B 63, 085104 (2001).

kotliar86 G. Kotliar and A. E. Ruckenstein, New functional integral approach to strongly correlated Fermi systems: The Gutzwiller approximation as a saddle point, Phys. Rev. Lett. 57, 1362 (1986).

lilly90 L. Lilly, A. Muramatsu, and W. Hanke, Slave-boson mean field versus quantum Monte Carlo results for the Hubbard model, Phys. Rev. Lett. 65, 1379 (1990).

fazekas99 P. Fazekas, Lecture notes on electron correlation and magnetism, Chapter 5 (World Scientific, Singapore, 1999).


supp2 A formal derivation of the large-U effective Hamiltonian and dynamics can be found in the supplemental material.

yoshioka09 T. Yoshioka, A. Koga, N. Kawakami, Quantum phase transitions in the Hubbard model on a triangular lattice, Phys. Rev. Lett. 103, 036401 (2009).

sahebsara08 P. Sahebsara and D. Sénéchal, Hubbard model on the triangular lattice: spiral order and spin liquid, Phys. Rev. Lett. 100, 136402 (2008).


tocchio13 L. F. Tocchio, H. Feldner, F. Becca, R. Valentí, C. Gros, Spin-liquid versus spiral-order phases in the anisotropic triangular lattice, Phys. Rev. B 87, 035143 (2013).

goto16 S. Goto, S. Kurihara, and D. Yamamoto, Incommensurate spiral magnetic order on anisotropic triangular lattice: Dynamical mean-field study in a spin-rotating frame, Phys. Rev. B 94, 245145 (2016)


dirac30 P. A. M. Dirac, Note on Exchange Phenomena in the Thomas Atom, Proc. Cambridge Philos. Soc. 26, 376 (1930).

frenkel34 J. Frenkel, Wave Mechanics: Advanced General Theory (Clarendon, Oxford, 1934).

mclachlan64 A. D. McLachlan and M. A. Ball,  Time-Dependent Hartree-Fock Theory for Molecules, Rev. Mod. Phys. 36, 844 (1964).

deumens94 E. Deumens, A. Diz, R. Longo, and Y. Öhrn, Time-dependent theoretical treatments of the dynamics of electrons and nuclei in molecular systems, Rev. Mod. Phys. 66, 917 (1994).

degennes P. G. de Gennes, Superconductivity of Metals and Alloys, (W. A. Benjamin, Inc., New York, 1966).


volkov74 A. F. Volkov and Sh. M. Kogan, Collisionless relaxation of the energy gap in superconductors, Sov. Phys.-JETP 38, 1018 (1975).

barankov04 R. A. Barankov, L. S. Levitov, and B. Z. Spivak, Collective Rabi oscillations and solitons in a time-dependent BCS pairing problem, Phys. Rev. Lett. 93, 160401 (2004).

tsuji13 N. Tsuji, M. Eckstein, and P. Werner, Nonthermal antiferromagnetic order and nonequilibrium criticality in the Hubbard model, Phys. Rev. Lett. 110, 136404 (2013).

sandri13 M. Sandri and M. Fabrizio, Nonequilibrium dynamics in the antiferromagnetic Hubbard model, Phys. Rev. B 88, 165113 (2013).

blinov17 I. V. Blinov, P. Ribeiro, and A. N. Rubtsov, Post-quench dynamics and suppression of thermalization in an open half-filled Hubbard layer, Phys. Rev. B 95, 024309 (2017).

davidson72 See, for example, R. C. Davidson, Methods in Nonlinear Plasma Theory (Academic, New York, 1972), Sec. 1.3.

LL L. D. Landau and L. M. Lifshitz, On the theory of the dispersion of magnetic permeability in ferromagnetic bodies. Physik. Zeits. Sowjetunion 8, 153-169 (1935).

supp3 Details of the ED calculation are discussed in the supplemental material, section C.

note_A In actual computation, floating point inaccuracy produces erratic dynamics even at “zero” temperature.

chubukov92 A. V. Chubukov and D. M. Frenkel, Renormalized perturbation theory of magnetic instabilities in the two-dimensional Hubbard model at small doping, Phys. Rev. B 46, 11884 (1992).

canali93 C. M. Canali and M. Wallin, Spin-spin correlation functions for the square-lattice Heisenberg antiferromagnet at zero temperature, Phys. Rev. B 48, 3264 (1993).



rowe70 D. J. Rowe, Nuclear Collective Motion (Methuen, London, 1970).

pines88 D. Pines and P. Noziéres, The Theory of Quantum Liquids: 1 Normal Fermi Liquids (Addison-Wesley, New York, 1988).

reinhard92 P.-G. Reinhard and Y. K. Gambhir, RPA in wavefunction representation,  Ann. Phys. 1, 598 (1992).


heidarian04 D. Heidarian and N. Trivedi, Inhomogeneous metallic phase in a disordered Mott insulator in two dimensions, Phys. Rev. Lett. 93, 126401 (2004).

lahoud14 E. Lahoud, O. N. Meetei, K. B. Chaska, A. Kanigel, N. Trivedi, Emergence of a novel pseudogap metallic state in a disordered 2D Mott insulator, Phys. Rev. Lett. 112, 206402 (2014).

mourigal13 M. Mourigal, W. T. Fuhrman, A. L. Chernyshev, M. E. Zhitomirsky, Dynamical structure factor of the triangular-lattice antiferromagnet, Phys. Rev. B 88, 094407 (2013).

vojta13 M. Vojta, Excitation Spectra of Disordered Dimer Magnets Near Quantum Criticality, Phys. Rev. Lett. 111, 097202 (2013).


byczuk09 K. Byczuk, W. Hofstetter, and D. Vollhardt, Competition between Anderson localization and antiferromagnetism in correlated lattice fermion systems with disorder, Phys. Rev. Lett. 102, 146403 (2009).


gilbert55 T. L. Gilbert, A Lagrangian formulation of the gyromagnetic equation of the magnetic field, Phys. Rev. 100, 1243 (1955).

solontsov93 A.Z. Solontsov, A.N. Vasil'ev, Spin fluctuation damping of magnons in itinerant electron magnets, Phys. Lett. A 177, 362 (1993).

barankov06 R. A. Barankov and L. S. Levitov, Synchronization in the BCS Pairing Dynamics as a Critical Phenomenon, Phys. Rev. Lett. 96, 230403 (2006).


ma12 P.-W. Ma and S. L. Dudarev, Longitudinal magnetic fluctuations in Langevin spin dynamics, Phys. Rev. B 86, 054416 (2012).

kirilyuk10 A. Kirilyuk, A. V. Kimel, and T. Rasing, Ultrafast optical manipulation of magnetic order, Rev. Mod. Phys. 82, 2731 (2010).

wingreen86 N. S. Wingreen, C. J. Stanton, and J. W. Wilkins, Electron-Electron Scattering in Nondegenerate Semiconductors: Driving the Anisotropic Distribution toward a Displaced Maxwellian, Phys. Rev. Lett. 57, 1084 (1986).

allen87 P. B. Allen, Theory of thermal relaxation of electrons in metals, Phys. Rev. Lett. 59, 1460 (1987).

suarez95 C. Suárez, W. E. Bron, and T. Juhasz, Dynamics and Transport of Electronic Carriers in Thin Gold Films, Phys. Rev. Lett. 75, 4536 (1995).

fatti00 N. Del Fatti, C. Voisin, M. Achermann, S. Tzortzakis, D. Christofilos, and F. Vallée, Nonequilibrium electron dynamics in noble metals, Phys. Rev. B 61, 16956 (2000).

niu99 Q. Niu, X. Wang, L. Kleinman, W.-M. Liu, D. M. C. Nicholson, and G. M. Stocks, Adiabatic dynamics of local spin moments in itinerant magnets, Phys. Rev. Lett. 83, 207 (1999).

capelle01 K. Capelle, G. Vignale, and B. L. Györffy, Spin currents and spin dynamics in time-dependent density-functional theory, Phys. Rev. Lett. 87, 206403 (2001).

qian02 Z. Qian and G. Vignale, Spin dynamics from time-dependent spin-density-functional theory, Phys. Rev. Lett. 88, 056404 (2002).



martin08 I. Martin and C. D. Batista, Itinerant electron-driven chiral magnetic ordering and spontaneous quantum Hall effect in triangular lattice models, Phys. Rev. Lett., 101, 156402 (2008).

nandkishore12 R. Nandkishore, G.-W. Chern, A. V. Chubukov, Itinerant half-metal spin-density-wave state on the hexagonal lattice, Phys. Rev. Lett. 108, 227204 (2012). 

chern12 G.-W. Chern and C. D. Batista, Spontaneous quantum Hall effect via a thermally induced quadratic Fermi point, Phys. Rev. Lett. 109, 156801 (2012).


nicholson16 C. W. Nicholson, C. Monney, R. Carley, B. Frietsch, J. Bowlan, M. Weinelt, and M. Wolf, Ultrafast spin density wave transition in Chromium governed by thermalized electron gas, Phys. Rev. Lett. 117, 136801 (2016).

anisimov74 S.I. Anisimov, B.L. Kapeliovich, T.L. Perel'man, Electron emission from metal surfaces exposed to ultrashort laser pulses, Sov. Phys. JETP 39, 375 (1974).

note1 Our equilibrium QLLD simulations find that both the triple-𝐐 SDW order parameters and the energy gap vanish above a transition temperature T^* ≈ 0.0015.

weisse06 A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, The kernel polynomial method, Rev. Mod. Phys. 78, 275 (2006).

garcia15 J. H. Garcia, L. Covaci, and T. G. Rappoport, Real-space calculation of the conductivity tensor for disordered topological matter, Phys. Rev. Lett. 114, 116602 (2015).

perfetti06 L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen, H. Berger, S. Biermann, P. S. Cornaglia, A. Georges, and M. Wolf, Time evolution of the electronic structure of 1T-TaS_2 through the insulator-metal transition, Phys. Rev. Lett. 97, 067402 (2006).

petersen11 J. C. Petersen, S. Kaiser, N. Dean, A. Simoncig, H. Y. Liu, A. L. Cavalieri, C. Cacho, I. C. E. Turcu, E. Springate, F. Frassetto, L. Poletto, S. S. Dhesi, H. Berger, and A. Cavalleri, Clocking the melting transition of charge and lattice order in 1T-TaS_2 with ultrafast extreme-ultraviolet angle-resolved photoemission spectroscopy, Phys. Rev. Lett. 107, 177402 (2011).

hellmann10 S. Hellmann, M. Beye, C. Sohrt, T. Rohwer, F. Sorgenfrei, H. Redlin, M. Kallne, M. Marczynski-Bhlow, F. Hennies, M. Bauer, A. Fhlisch, L. Kipp, W. Wurth, and K. Rossnagel, Ultrafast melting of a charge-density wave in the Mott insulator 1T-TaS_2, Phys. Rev. Lett. 105, 187401 (2010).

shen14 W. Shen, Y. Ge, A. Y. Liu, H. R. Krishnamurthy, T. P. Devereaux, and J. K. Freericks, Nonequilibrium “Melting” of a charge density wave insulator via an ultrafast laser pulse, Phys. Rev. Lett. 112, 176404 (2014).

chernyshev06 A. L. Chernyshev and M. E. Zhitomirsky, Magnon Decay in Noncollinear Quantum Antiferromagnets, Phys. Rev. Lett. 97, 207202 (2006).

starykh06 O. A. Starykh, A. V. Chubukov, and A. G. Abanov, Flat spin-wave dispersion in a triangular antiferromagnet, Phys. Rev. B 74, 180403(R) (2006).

zhitomirsky13 M. E. Zhitomirsky and A. L. Chernyshev, Rev. Mod. Phys. 85, 219 (2013).

ma16 See, e.g. J. Ma, Y. Kamiya, Tao Hong, H. B. Cao, G. Ehlers, W. Tian, C. D. Batista, Z. L. Dun, H. D. Zhou, M. Matsuda, Static and Dynamical Properties of the Spin-1/2 Equilateral Triangular-Lattice Antiferromagnet Ba_3CoSb_2O_9, Phys. Rev. Lett. 116, 087201 (2016).

dagotto05 E. Dagotto, Complexity in Strongly Correlated Electronic Systems, Science 309, 257 (2005).

fradkin15 E. Fradkin, S. A. Kivelson, and J. M. Tranquada, Theory of intertwined orders in high temperature superconductors, Rev. Mod. Phys. 87, 457 (2015).












Supplemental Materials





§ A. KERNEL POLYNOMIAL METHOD AND ITS GRADIENT TRANSFORMATION


The kernel polynomial method (KPM) <cit.> and the technique of automatic differentiation <cit.> are crucial to our implementation of efficient QLLD simulations.  Here we briefly review these numerical techniques. Conventional KPM provides an efficient approach to computing the system free-energy ℱ. However, central to our QLLD simulations is the calculation of the `forces' acting on spins: ∂ℱ / ∂𝐦_i, where the effective energy functional is calculated from the quadratic fermion Hamiltonian: ℱ = ⟨ℋ_ SDW⟩. Specifically, the force is

    ∂ℱ/∂𝐦_i = -U ⟨𝐬_i ⟩ = -U/2ρ_iα, iβ σ_βα.

Computing the force is thus equivalent to evaluating the single-particle density matrix ρ_iα, jβ. We first introduce the single-particle Hamiltonian H_iα, jβ such that the quadratic spin-fermion Hamiltonian is expressed as (up to a term that is independent of fermions)

    ℋ_ SDW = ∑_iα,j β H_iα,jβ c^†_iα c^ _jβ = ∑_IJ H_IJ c^†_I   c^ _J.

Here we have introduced notation I = (i, α), J = (j, β), ⋯ for simplicity. The density matrix is then given by the derivative

    ρ_IJ = ⟨ c^†_J  c^ _I ⟩ = ∂ℱ/∂ H_IJ.

Next we outline the KPM procedure for computing the free energy which is expressed in terms of the DOS as ℱ = ∫ρ(ϵ)   f(ϵ) dϵ, where f(ϵ) = -T log [1 + e^-(ϵ-μ)/T]. KPM begins by approximating the DOS as a Chebyshev polynomial series,

    ρ(ϵ) = 1/π√(1 - ϵ^2)∑_m=0^M-1 (2 - δ_0, m)  μ_m  T_m(ϵ),

where T_m(x) are Chebyshev polynomials, and μ_m are the expansion coefficients. The expansion is valid only when all eigenvalues of H_IJ have magnitude less than one. This can in general be achieved through a simple shifting and rescaling of the Hamiltonian. Moreover, damping coefficients g_m are often introduced to reduce the unwanted artificial Gibbs oscillations. Substituting ρ(ϵ) into the free energy expression gives

    ℱ = ∑_m=0^M-1 C_m μ_m,

where coefficients C_m = (2 - δ_0,m) g_m ∫_-1^1T_m(ϵ) f(ϵ)/π√(1-ϵ^2) dϵ are independent of the Hamiltonian and may be efficiently evaluated using Chebyshev-Gauss quadrature. 

The key step of KPM is to replace computation of the Chebyshev moments μ_m =  Tr  T_m(H) by an ensemble average μ_m = ⟨ T_m(H) ⟩ =  1/R∑_ℓ=1^R r_ℓ^† H r_ℓ over random normalized column vectors r <cit.>. Taking advantage of the recursive relation of Chebyshev polynomials: T_m(H) = 2 H · T_m-1(H) - T_m-2(H), the moments can be evaluated recursively as follows:

    μ_m = r^†·α_m,

where r is a random vectors with complex elements drawn from the uniform distribution |r_I|^2 = 1. The random vectors α_m are given by

    α_m = {[                  r,                                   m = 0;              H · r,                                     m=1; 2 H ·α_m-1 - α_m-2,                                     m>1;                     ].

The above recursion relation also indicates that evaluation of μ_m that are required for computing ℱ only involves matrix-vector products. For sparse matrix H with 𝒪(N) elements, this requires only 𝒪(M N) operations, where M is the number of Chebyshev polynomials. On the other hand, even with the efficient algorithm for ℱ, a naive calculation of the derivatives ∂ℱ/∂ H_IJ based on finite difference approximation is not only inefficient but also inaccurate. The computational cost of finite difference is similar to the KPM-based Monte Carlo method with local updates.

To circumvent this difficulty, we employ the technique of automatic differentiation with reverse accumulation <cit.>. Instead of directly using Eq. (<ref>), the trick is to view ℱ as a function of vectors α_m and write

    ∂ℱ/∂ H_IJ = ∑_m=0^M-1∂ℱ/∂α_m, K∂α_m, K/∂ H_IJ,

Here α_m, K denotes the K-th component of vector α_m, and summation over the repeated index K is assumed. Using Eq. (<ref>), we have 

    ∂α_0, K/∂ H_IJ = 0,     ∂α_1, K/∂ H_IJ = δ_IK α_0, J,     ∂α_m, K/∂ H_IJ = 2 δ_IK α_m-1, J   (m>1)

The expression of ∂ℱ/∂ H_IJ can be simplified by introducing a new set of random vectors:

    β_m ≡∂ℱ/∂α_m+1,

From Eqs. (<ref>) and (<ref>) , we obtain

    ∂ℱ/∂ H_IJ = β_0, I α_0, J + 2 ∑_m=1^M-2β_m, I α_m, J.

Remarkably, the vectors β_m can also be computed recursively. To this end, we note that the recursion relation (<ref>) implies that ℱ depends on α_m through three paths:

    ∂ℱ/∂α_m, K = ∂ℱ/∂μ_m∂μ_m/∂α_m, K
    	+ ∂ℱ/∂α_m+1, L∂α_m+1, L/∂α_m, K
    	+ ∂ℱ/∂α_m+2, L∂α_m+2, L/∂α_m, K.

The various terms above can be straightforwardly calculated:

    ∂ℱ/∂μ_m = C_m,   ∂μ_m/∂α_m, K = r^*_K,   ∂α_m+1, L/∂α_m, K = 2 H_LK,   ∂α_m+2, L/∂α_m, K = -δ_LK.

Consequently, 

    β_m = C_m+1  r^† + 2 β_m+1· H - β_m+2,      (m < M-1).

Restoring the site and spin indices, we obtain the following expression for the density matrix

    ρ_iα, jβ = β_0, iα α_0, jβ + 2 ∑_m=1^M-2β_m, iα α_m, jβ.

As in standard KPM, there are two independent sources of errors in our method <cit.>: the truncation of the Chebyshev series at order M-1, and the stochastic estimation of the moments using finite number R of random vectors.  The performance of the stochastic estimation can be further improved using correlated random vectors based on the probing method <cit.>. 
Most simulations discussed in the main text were done on a 120× 120 triangular lattice. The number of Chebyshev polynomials used in the simulations is in the range of M = 1000 to 2000. The number of correlated random vectors used is R = 64 to 144.



§ B. LARGE-U LIMIT: FORMAL DERIVATION
 

Here we present a formal derivation of the effective SDW Hamiltonian in the large-U limit, which is expected to be equivalent to that of the original Hubbard model. Our first step is to write the spin-fermion Hamiltonian Eq. (<ref>) in a new reference frame, such that the local quantization axis of site i coincides with the direction of the SDW field. Let 𝐦_i = |𝐦_i| (sinθ_i cosϕ_i, sinθ_i sinϕ_i, cosθ_i), the fermionic operators in the new reference frame are

    c̃^†_i, +   =    e^-i ϕ_i/2sin(θ_i/2)   c^†_i, ↓ + e^i ϕ_i /2cos(θ_i/2)  c^†_i, ↑,  
    c̃^†_i, -   =    e^-i ϕ_i/2cos(θ_i/2)   c^†_i, ↓ - e^i ϕ_i /2sin(θ_i/2)  c^†_i, ↑,

The inverse transformation is

    c^†_i, ↑ = e^-i ϕ_i/2cos(θ_i/2) c̃^†_i, + - e^-i ϕ_i /2sin(θ_i/2) c̃^†_i, -, 
    
    	c^†_i, ↓ = e^+i ϕ_i/2sin(θ_i/2)  c̃^†_i, + + e^+ i ϕ_i/2cos(θ_i/2)  c̃^†_i, -.

We separate the spin-fermion Hamiltonian into two parts ℋ_ SDW = 𝒯 + 𝒰. The kinetic hopping term can be re-expressed in the new reference frame as 

    𝒯 = - ∑_⟨ ij ⟩∑_μν = ± t_ij(η^μν_ij c̃^†_iμ c̃^ _j ν +  h.c.),

where

    η^++_ij   =    +cos(ϕ_j - ϕ_i/2) cos(θ_j - θ_i/2) + i sin(ϕ_j - ϕ_i/2) cos(θ_j + θ_i/2), 
    η^+-_ij   =    -cos(ϕ_j - ϕ_i/2) cos(θ_j - θ_i/2) - i sin(ϕ_j - ϕ_i/2) cos(θ_j + θ_i/2),  
    η^-+_ij   =    -cos(ϕ_j - ϕ_i/2) cos(θ_j - θ_i/2) + i sin(ϕ_j - ϕ_i/2) cos(θ_j + θ_i/2),  
    η^–_ij   =    +cos(ϕ_j - ϕ_i/2) cos(θ_j - θ_i/2) - i sin(ϕ_j - ϕ_i/2) cos(θ_j + θ_i/2).

To derive the effective Hamiltonian, we use the standard perturbation approach by treating the hopping 𝒯 as a perturbation to the coupling term

    𝒰 = - U ∑_i |𝐦_i| ( c̃^†_i,+c̃^ _i, + - c̃^†_i,-c̃^ _i, -) + U ∑_i |𝐦_i|^2.

For convenience, we first introduce the resolvent of  𝒰: Ĝ_0(ϵ) = 1/(ϵ - 𝒰 ). The effective Hamiltonian up to second order in t_ij is given by

    ℋ_ eff = 𝒫 𝒯Ĝ_0(ϵ_0) 𝒯𝒫,

where 𝒫 is a projector onto the lowest energy subspace with one electron per site whose spin is parallel to local moment 𝐦_i, and ϵ_0 = - NU/4 is the energy of the degenerate large-U ground state (N is the number of lattice sites). Then, in the new reference frame, the only processes contributing to ℋ_ eff are the spin-flip hoppings which annihilate electrons with spin μ=+ and create electrons in a difference site with spin μ=-. Given that each state of the lowest energy subspace is fully characterized by the field configuration {𝐦_i }, the effective Hamiltonian can be expressed in terms of the the SDW field:

    ℋ_ eff   =    -∑_⟨ ij ⟩t_ij^2/U( |η^-+_ij|^2 + |η^+-_ij|^2 ) 
        =    ∑_⟨ ij ⟩ t_ij^2/U[ 1 - sinθ_i sinθ_j cos(ϕ_i - ϕ_j) - cosθ_i cosθ_j ]  
    	 =  ∑_⟨ ij ⟩4 t_ij^2/U(𝐦_i ·𝐦_j - 1/4)

which is the expected Heisenberg exchange interaction for spin-1/2 in the large U limit <cit.>.

We next derive the dynamics equation in the large-U limit, which is given by the Landau-Lifshitz equation. To this end, we first consider the Heisenberg equation of motion for local spin operator d 𝐬_i / dt = -i [𝐬_i, ℋ_ H]. Here 𝐬_i = 1/2 c^†_iασ^ _αβ c^ _iβ, and ℋ_ H is the Hubbard Hamiltonian. Expressing the on-site U interaction in terms of spin operators, it is given by
    ℋ_ H = - ∑_⟨ ij ⟩, α t_ij(c^†_iα c^ _jα +  h.c.) - 2U/3∑_i 𝐬_i^2 + N_e U/2.

where N_e is the number of electrons. Obviously, the U term of the Hubbard Hamiltonian commutes with the spin operator, and we have d 𝐬_i / dt =-i [𝐬_i, 𝒯]. For example, we consider the z component first. Using commutation relations [s^z_i, c^†_i, ↑] = 1/2 c^†_i, ↑ and [s^z_i, c^†_i, ↓ ] = -1/2 c^†_i,↓, we have

    -i [s^z_i, (c^†_i α c^ _j α + c^†_jα c^ _iα) ] 
    	= -i σ_α/2(c^†_i α c^ _j α - c^†_jα c^ _iα) 
    	= -i/2(c^†_i ασ̂^z_αβ  c^ _j β - c^†_jασ̂^z_αβ  c^ _iβ).

Here σ_α = ± 1 for α = ↑, ↓, respectively. Summing over repeated indices is also implied. The equation for the x and y components of 𝐬_i can be obtained by applying π/2 rotations to this equation. We have

    d𝐬_i/dt = - ∑_j 𝐉_ij,

where we have defined the spin current density operator 𝐉_ij. 

    𝐉_ij = -i σ_αβ/2(c^†_iα c^ _jβ - c^†_jα c^ _iβ).

It is worth noting that Eq. (<ref>) is simply the continuity equation for the spin density. Taking the expectation value with respect to the ground state gives rise to the equation of motion for the SDW field

    d𝐦_i/dt = -∑_j⟨𝐉_ij⟩,

Next we compute the expectation value of 𝐉_ij in the large U limit. In the t / U → 0 limit, obviously ⟨𝐉_ij⟩ = 0. A nonzero contribution comes from the second-order perturbation due to electron hopping. The procedure is similar to what we did to derive the effective Hamiltonian. Specifically, we project the spin current operator into the degenerate low-energy manifold of ℋ_ SDW with the electronic eigenstates corrected up to first order in the perturbation. For the z-component, we have

    ⟨ J^z_ij⟩ = 1/2[ 𝒫 J^z_ij G_0(ϵ_0) 𝒯𝒫 + 𝒫𝒯 G_0(ϵ_0) J^z_ij𝒫].

Once again, it is convenient to work in the new reference frame. For example, the current density operator

    J^z_ij = -i t_ij/2∑_μν = ±(τ^μν_ijc̃^†_i μ c̃^ _jν - (τ^μν_ij)^* c̃^†_jν c^ _iμ).

Here we have introduced

    τ^++_ij   =    +cos(ϕ_j - ϕ_i/2) cos(θ_j + θ_i/2) + i sin(ϕ_j - ϕ_i/2) cos(θ_j - θ_i/2), 
    τ^+-_ij   =    -cos(ϕ_j - ϕ_i/2) cos(θ_j + θ_i/2) - i sin(ϕ_j - ϕ_i/2) cos(θ_j - θ_i/2),  
    τ^-+_ij   =    -cos(ϕ_j - ϕ_i/2) cos(θ_j + θ_i/2) + i sin(ϕ_j - ϕ_i/2) cos(θ_j - θ_i/2),  
    τ^–_ij   =    +cos(ϕ_j - ϕ_i/2) cos(θ_j + θ_i/2) - i sin(ϕ_j - ϕ_i/2) cos(θ_j - θ_i/2).

Using expressions (<ref>) and similar one for the kinetic term Eq. (<ref>), we obtain

    ⟨ J^z_ij⟩ = -i t_ij^2/U{[ τ^-+_ij (η^-+_ij)^* - (τ^-+_ij)^* η^-+_ij]
    	+ [ τ^+-_ij (η^+-_ij)^* - (τ^+-_ij)^* η^+-_ij]
    	 }

Using the definitions for τ^μν_ij and η^μν_ij, it can be shown that

    [ τ^-+_ij (η^-+_ij)^* - (τ^-+_ij)^* η^-+_ij]    =   [ τ^+-_ij (η^+-_ij)^* - (τ^+-_ij)^* η^+-_ij] 
       =    2 i sin(ϕ_j - ϕ_i) [ sin^2 (θ_j + θ_i/2) - sin^2 (θ_j - θ_i/2) ] 
       =    2 i sin(ϕ_j - ϕ_i) sinθ_i sinθ_j

In terms of the SDW field 𝐦_i, the right-hand side of the above equation is 2 i  𝐦̂_i ×𝐦̂_j ·𝐳̂, where 𝐦̂_i is a unit vector along the local moment direction. Using the fact that |𝐦_i| = 1/2 in the large U limit at half filling, this result indicates the following vector identity for the spin current

    ⟨𝐉_ij⟩ = 4 t_ij^2/U𝐦_i ×𝐦_j.

Substituting this into Eq. (<ref>) gives the well known Landau-Lifshitz equation of motion (<ref>) in the main text for the Heisenberg exchange Hamiltonian.



§ C. EXACT DIAGONALIZATION CALCULATION OF DYNAMICAL STRUCTURE FACTOR


Here we provide details of the exact diagonalization (ED) calculation of Hubbard Model on 4×4 square lattice with periodic boundary condition (PBC). For simplicity, we consider the case where SU(2) symmetry is conserved in the model, which leads to 𝒮(q,ω)=3𝒮^zz(q,ω). 

To calculate 𝒮^zz(q,ω) at T=0, we first obtain the ground state |Ψ_0⟩ in the total S_z=0 sector at half-filling, by using the implicitly restarted Arnoldi method provided through the ARPACK libary <cit.>. The dynamical structure factor can be expressed through the fluctuation-dissipation theorem:

    S^zz(q,ω)    =    -2 Imχ_zz(q,ω) 
       =    -2 ⟨Ψ_0  | S_q^z S_-q^z | Ψ_0⟩ Im⟨ϕ_0 | (ω + i η + E_0-H)^-1 | ϕ_0 ⟩,

where | ϕ_0 ⟩≡ S^z(-q)|Ψ_0⟩ / √(⟨Ψ_0  | S_q^z S_-q^z | Ψ_0⟩), and E_0 is the ground state energy: H |Ψ_0 ⟩ = E_0 | Ψ_0 ⟩.

The matrix inverse (z-H)^-1 can be calculated through the Lanczos algorithm <cit.>:
1em
1em

In the new basis { |ϕ_0⟩, | ϕ_1 ⟩, | ϕ_2 ⟩,…}, the Hamiltonian is expressed by a tridiagonal matrix:

    H= 
    [ a_0 b_1    ; b_1 a_1 b_2;     b_2 a_2   ⋱;           ⋱   ⋱ b_n;             b_n a_n ].


With Cramer's rule, the first element of the inverse matrix can be expressed as a continued fraction:

    ⟨ϕ_0 | (z-H)^-1 | ϕ_0 ⟩ =
    [ (z-a_0)-b_1^2/(z-a_1)-b_2^2/(z-a_2)-⋯]^-1,

which leads to

    𝒮^zz(q,ω)     =-2 ⟨Ψ_0  | S_q^z S_-q^z | Ψ_0⟩·Im[ (z-a_0)-b_1^2/(z-a_1)-b_2^2/(z-a_2)-⋯]^-1,

where z≡ω+iη+E_0.


99

weisse06b A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, The kernel polynomial method, Rev. Mod. Phys. 78, 275 (2006).

barros13b K. Barros and Y. Kato, Efficient Langevin simulation of coupled classical fields and fermions, Phys. Rev. B 88, 235101 (2013).

silver94 R. N. Silver and H. Röder, Density of states of mega-dimensional Hamiltonian matrices, Int. J. Mod. Phys. C 5, 735 (1994).

griewank89 A. Griewank, in Mathematical Programming: Recent Developments and Applications, edited by M. Iri and K. Tanabe (Kluwer Academic, Dordrecht, The Netherlands, 1989), pp. 83–108.

tang12 J. M. Tang and Y. Saad, A probing method for computing the diagonal of a matrix inverse, Numer. Linear Algebra Appl. 19, 485 (2012).

fazekas99b P. Fazekas, Lecture notes on electron correlation and magnetism, Chapter 5 (World Scientific, Singapore, 1999).

lehoucq98 R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK Users' Guide: Solution of Large-scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods (Society for Industrial and Applied Mathematics, 1998).

lanczos90 C. Lanczos, An Iteration Method for the Solution of the Eigenvalue Problem of Linear Differential and Integral Operators, J. Res. Nat. Bur. Stand. 45, 255 (1950).

dagotto94 E. Dagotto, Correlated electrons in high-temperature superconductors, Rev. Mod. Phys. 66, 763 (1994).




