Nonlinear Fluctuating Hydrodynamics for the Classical XXZ Spin Chain

TL;DR

This paper applies nonlinear fluctuating hydrodynamics to the classical XXZ spin chain, analyzing equilibrium correlations and verifying predictions through molecular dynamics simulations, revealing regimes with KPZ scaling, diffusive behavior, and ballistic transport.

Contribution

It provides the first detailed molecular dynamics verification of NFH predictions for the classical XXZ spin chain across different temperature regimes.

Findings

Observation of KPZ scaling in sound peaks at intermediate temperatures

Identification of diffusive peaks at high temperatures

Detection of ballistic transport in the integrable model

Abstract

Using the framework of nonlinear fluctuating hydrodynamics (NFH), we examine equilibrium spatio-temporal correlations in classical ferromagnetic spin chains with nearest neighbor interactions. In particular, we consider the classical XXZ-Heisenberg spin chain (also known as Lattice Landau Lifshitz or LLL model) evolving deterministically and chaotically via Hamiltonian dynamics, for which energy and $z$-magnetization are the only locally conserved fields. For the easy-plane case, this system has a low-temperature regime in which the difference between neighboring spin's angular orientations in the XY plane is an \textit{almost conserved} field. According to the predictions of NFH, the dynamic correlations in this regime exhibit a heat peak and propagating sound peaks, all with anomalous broadening. We present a detailed molecular dynamics test of these predictions and find a reasonably accurate verification. We find that, in a suitable intermediate temperature regime, the system shows two sound peaks with Kardar-Parisi-Zhang (KPZ) scaling and a heat peak where the expected anomalous broadening is less clear. In high temperature regimes of both easy plane and easy axis case of LLL, our numerics show clear diffusive spin and energy peaks and absence of any sound modes, as one would expect. We also simulate an integrable version of the XXZ-model, for which the ballistic component instead moves with a broad range of speeds rather than being concentrated in narrower peaks around the sound speed.

Tables

Table 1. Table 1: G 𝐺 G matrices ( β = 10 , h = 0.3 , ν = 0 formulae-sequence 𝛽 10 formulae-sequence ℎ 0.3 𝜈 0 \beta=10,~{}h=0.3,~{}\nu=0 )

$G$ matrix	MD simulation	NFH
$G^{-}$
$G^0$
$G^{+}$
speed of sound $c$	0.865 (maxima of the sound peaks) 0.859 (normal mode transformation)	0.85217

Table 2. Table 2: G matrices ( β = 10 , h = 0 , ν = 0 formulae-sequence 𝛽 10 formulae-sequence ℎ 0 𝜈 0 \beta=10,~{}h=0,~{}\nu=0 )

$G$ matrix	MD simulation	NFH
$G^{-}$
$G^0$
$G^{+}$
speed of sound $c$	$\begin{array}{c}0.931\text{ (maxima of the sound peaks)}\\ 0.930\text{ (normal mode transformation)}\end{array}$	0.92837

Table 3. Table 3: Summary of the transport properties observed in this work.

Model and the parameter regime	Observation
(Non-integrable LLL chain)
Easy plane at high T	Diffusive spin and energy peaks, no sound modes
Easy plane at low T and $h=0$	Diffusive sound modes∗ and Levy-3/2 heat mode∗
Easy plane at low T and small h	KPZ sound modes and Levy-5/3 heat mode∗
Easy axis at high T	Diffusive spin and energy peaks, no sound modes
(Integrable LLL chain)
Easy plane at high T	Ballistic spin and energy peaks

Equations

$$H=-\sum_{j\in\mathbb{Z}}\big{(}S_{j}^{x}S_{j+1}^{x}+S_{j}^{y}S_{j+1}^{y}+\Delta S_{j}^{z}S_{j+1}^{z}\big{)},$$

$$\frac{d}{dt}{\vec{S}}_j=\{{\vec{S}}_j,H\}={\vec{S}}_j\times {\vec{B}}_j,{\vec{B}}_j=-{\nabla }_{{\vec{S}}_j}H,$$

$$S_j^x=f(s_j)\cos {\varphi }_j,S_j^y=f(s_j)\sin {\varphi }_j,S_j^z=s_j,$$

$$H=-\sum_{j\in\mathbb{Z}}\Big{(}f(s_{j})f(s_{j+1})\cos(\phi_{j+1}-\phi_{j})+\Delta s_{j}s_{j+1}\Big{)}.$$

$$\begin{array}{rl}& \frac{d}{dt}{\varphi }_j=-f^{\prime }(s_j)f(s_{j+1})\cos ({\varphi }_{j+1}-{\varphi }_j)-f(s_{j-1})f^{\prime }(s_j)\cos ({\varphi }_j-{\varphi }_{j-1})\\ & -\Delta (s_{j-1}+s_{j+1}),\end{array}$$

$$\frac{d}{dt}{s}_j=f(s_j)f(s_{j+1})\sin ({\varphi }_{j+1}-{\varphi }_j)-f(s_{j-1})f(s_j)\sin ({\varphi }_j-{\varphi }_{j-1}).$$

$$e_j=-f(s_j)f(s_{j+1})\cos (r_j)-\Delta s_j{s}_{j+1},$$

$${\mathcal{J}}_j^{\mathsf{s}}$$

$${\mathcal{J}}_j^{\mathsf{e}}$$

$$+\Delta f(s_j)f(s_{j+1})s_{j-1}\sin (r_j)+\Delta f(s_j)f(s_{j-1})s_{j+1}\sin (r_{j-1}).$$

$$Z_{N}(\beta,h)^{-1}\exp\Big{[}-\beta\Big{(}H-h\sum_{j}s_{j}\Big{)}\Big{]}\prod_{j}dr_{j}ds_{j},\quad\beta\geq 0,\,\,h\in\mathbb{R},$$

$$⟨{\mathcal{J}}_j^{\mathsf{s}}{⟩}_{\beta ,h}=0,⟨{\mathcal{J}}_j^{\mathsf{e}}{⟩}_{\beta ,h}=0.$$

$$C_{\mathrm{ss}}(j,t)$$

$$C_{\mathrm{ee}}(j,t)$$

$$H=-\sum _{j\in \mathbb{Z}}h({\vec{S}}_j,{\vec{S}}_{j+1})$$

$$\displaystyle h(\vec{S},\vec{S^{\prime}})=\log\big{|}\cos(\rho S^{z})\cos(\rho S^{\prime z})+(\cot(\rho))^{2}\sin(\rho S^{z})\sin(\rho S^{\prime z})$$

$$\displaystyle\hskip 80.0pt+(\sin(\rho))^{-2}G(S^{z})G(S^{\prime z})(S^{x}S^{\prime x}+S^{y}S^{\prime y})\big{|},$$

$$\displaystyle G(x)=\big{(}1-x^{2}\big{)}^{-\tfrac{1}{2}}\big{(}\cos(2\rho x)-\cos(2\rho)\big{)}^{\frac{1}{2}}$$

$$h(\vec{S},\vec{S^{\prime}})=\log\big{(}1+\vec{S}\cdot\vec{S^{\prime}}\big{)}.$$

$$r_j={\varphi }_{j+1}-{\varphi }_j,\mathrm{mod}2\pi ,$$

$$\frac{d}{dt}{\tilde{r}}_j(t)=g_j(t)-g_{j+1}(t)$$

$${\tilde{r}}_j(t)=r_j(t)+2\pi n_j(t)$$

$$\sum_{j}\big{(}\tilde{r}_{j}(t)-\tilde{r}_{j}(0)\big{)}=\sum_{j}\int_{0}^{t}ds\big{(}g_{j}(s)-g_{j+1}(s)\big{)}=0,$$

$$\sum_{j}\big{(}r_{j}(t)-r_{j}(0)\big{)}=-2\pi\sum_{j}n_{j}(t)~{}.$$

$$e_{\mathrm{g}}(r,s)=-(1-s^2)\cos r-\Delta s^2-hs.$$

$$E_{\mathrm{slip}}=4\left(1-\frac{1}{4}(1-\Delta {)}^{-2}{h}^2\right)>1$$

$$E_{\mathrm{slip}}=2\left(1-\frac{h}{2}-\frac{1}{2}(1-\Delta {)}^{-1}\Delta h\right).$$

$$H_{\mathrm{lt}}=\sum_{j\in\mathbb{Z}}\big{(}f(s_{j})f(s_{j+1})U(r_{j})-\Delta s_{j}s_{j+1}\big{)}.$$

$$\frac{d}{dt}{\varphi }_j={\partial }_{s_j}{H}_{\mathrm{lt}},\frac{d}{dt}{s}_j=-{\partial }_{{\varphi }_j}{H}_{\mathrm{lt}}.$$

$$r_j,s_j,e_j=f_j{f}_{j+1}{U}_j-\Delta s_j{s}_{j+1},$$

$${\mathcal{J}}_j^{\mathsf{r}}=-f_j^{\prime }{f}_{j+1}{U}_j-f_{j-1}{f}_j^{\prime }{U}_{j-1}+\Delta (s_{j-1}+s_{j+1}),$$

$${\mathcal{J}}_j^{\mathsf{s}}=-f_{j-1}{f}_j{U}_{j-1}^{\prime },$$

$$\mathcal{J}_{j}^{\mathsf{e}}=-f_{j-1}f_{j}f^{\prime}_{j}f_{j+1}\big{(}U^{\prime}_{j-1}U_{j}+U_{j-1}U^{\prime}_{j}\big{)}+\,\Delta\big{(}f_{j-1}f_{j}U^{\prime}_{j-1}s_{j+1}+f_{j}f_{j+1}U^{\prime}_{j}s_{j-1}\big{)}.$$

Full text

1. 1 Introduction
2. 2 High-temperature, non-integrable and integrable Lattice Landau-Lifshitz equations
3. 3 Easy plane at low-temperatures
4. 4 Numerical results at low temperature
5. 5 Discrete time LLL dynamics
6. 6 Conclusions
7. 7 Acknowledgements
8. A Coupling coefficients for nonlinear fluctuating hydrodynamics
9. B Low temperature approximation

1 Introduction

For generic classical and quantum spin chains the only conservation law is the energy, perhaps in addition one spin component (or all three). Momentum conservation is destroyed by the underlying lattice and in thermal equilibrium the average currents vanish. One therefore expects that a local perturbation of the thermal state will spread diffusively, a behavior which is actually observed in a large variety of systems, as prototypical examples we refer toaoki2000bulk ; sirker2011conservation . There are obvious exceptions such as integrable spin chains, for which a small perturbation induces a ballistic response. Also, at least classically, at very low temperatures the harmonic approximation generally becomes valid, which then implies ballistic transport over a suitable time scale.

The goal of our paper is to explain that, beyond the standard folklore, there can be a parameter regime, in which the dynamic correlations consist of a “heat” peak at the origin and in addition two sound peaks symmetrically moving to the right and left. These peaks broaden sub-ballistically but faster than diffusion. The theoretical argument is based on nonlinear fluctuating hydrodynamics (NFH), which refers to long wavelength behavior and should thus be equally valid for both classical and quantum chains. For quantum chains, because of numerical limitations, it is difficult to pin down the phenomenon, even less so the precise scaling form. In this contribution we thus restrict our study to classical spin chains of XXZ type. We note that for off-lattice models, where there is momentum conservation (e.g anharmonic oscillator chains), the presence of sound peaks with anomalous scaling has been observed in several recent studies and is well understood within the framework of NFH MendlSpohn2013 ; das2014numerical . The results presented here for the XXZ model bear close resemblance to those seen for the rotor model, where early results in PhysRevLett.84.2144 ; PhysRevLett.84.2381 ; Yang2003 indicated a diffusive to super-diffusive transport, and has recently been understood in the framework of NFH das2014role ; spohn2014fluctuating . It should also be noted that there are no phase transitions (from diffusive to anomalous transport regimes), rather the anomalous scaling is observed over very long transient time-scales.

For our argument we require that the anisotropy parameter $\Delta=|J_{z}/J_{x}|=|J_{z}/J_{y}|$, where $J_{z}$ is the nearest neighbor coupling between $z$-components of spin and $J_{x}=J_{y}$ between the $x$ or $y$-components respectively, satisfies $\Delta<1$. Then, at high temperatures, energy and the $z$-component of the spin diffuse. However at low temperatures the spin motion is confined to a plane orthogonal to the $z$-axis (the easy-plane) and phase differences between neighboring spins are small. To achieve differences of order $\pi$ ( in other words “phase slips” or, equivalently, “umklapp”) is an activated process and is thus strongly thermally suppressed in a low temperature regime. In this regime, the phase differences are an almost conserved field, so there is a broad range of time scales where this conservation law dictates the hydrodynamics. Under such conditions, nonlinear fluctuating hydrodynamics can be applied. The theory predicts the dynamical correlations contain a central non-propagating heat peak and left- and right-moving sound peaks. The three peaks broaden as nontrivial powers of time according to characteristic explicitly known scaling functions SpohnAHC2014 ; Spohn2016book .

To illustrate the difference between sound peaks and ballistic broadening, we also simulate the integrable spin chain of Fadeev and Takhtajan faddeev2007hamiltonian . The infinite number of conservation laws then leads to a structured scaling function which scales self-similarly as $\sim t^{-1}f(x/t)$.

There has been previous work on non-integrable classical models (and their KPZ connection) such as the Fermi-Pasta-Ulam chain MendlSpohn2013 ; das2014numerical , the discrete nonlinear Schrödinger equation PhysRevA.88.021603 ; PhysRevA.92.043612 ; MendlSpohnNLS2015 , coupled rotors das2014role ; spohn2014fluctuating ; Spohn2016book , and one-dimensional hard-point systems PhysRevE.90.012147 . However, to our knowledge the present paper is the first exploration of NFH in classical spin chains.

2 High-temperature, non-integrable and integrable Lattice Landau-Lifshitz equations

LLL equations — We consider spins of unit length on the one-dimensional lattice, $\vec{S}_{j}=(S_{j}^{x},S_{j}^{y},S_{j}^{z})$ with $|\vec{S}_{j}|=1$, $j\in\mathbb{Z}$. The standard LLL interaction is quadratic in the spins with Hamiltonian

[TABLE]

$\Delta$ the asymmetry parameter, $\Delta\geq 0$. The LLL equations of motion then read

[TABLE]

where the Poisson bracket between two functions, $g_{1},g_{2}$, of the spin variables is defined by \{g_{1},g_{2}\}=\sum_{j}\epsilon_{\alpha\beta\gamma}\big{(}\partial g_{1}/\partial{S_{j}^{\alpha}}\big{)}\big{(}\partial g_{2}/\partial{S_{j}^{\beta}}\big{)}S_{j}^{\gamma} with the usual summation convention. Clearly $|\vec{S}_{j}(t)|=1$ for all times. The Hamiltonian character of the dynamics can be seen also by introducing the position-like angular variable $\phi_{j}\in S^{1}$ and the conjugate canonical momentum-like variable $s_{j}\in[-1,1]$ defined through

[TABLE]

where $f(x)=(1-x^{2})^{1/2}$. Indeed, one checks that $\{s_{i},\phi_{j}\}=\delta_{ij}$, $\{\phi_{i},\phi_{j}\}=0$, $\{s_{i},s_{j}\}=0$. In these variables the hamiltonian (2.1) reads

[TABLE]

Thus at low energies the phases tend to align, while $\Delta$ sets the interactions between the $z$ components. The isotropic model corresponds to $\Delta=1$, easy-plane to $\Delta<1$, and easy-axis to $\Delta>1$. In the new variables the equations of motion become

[TABLE]

For $|s_{j}(t)|<1$ the angle $\phi_{j}(t)$ is defined modulo $2\pi$. For $s_{j}(t)=\pm 1$ the angle $\phi_{j}(t)$ is ill-defined. However, in our context, trajectories where a spin precisely hits either the north or south pole of the unit sphere ($|s_{j}(t)|=1$) have measure zero and therefore can be ignored.

The $s$-$s$ interaction could have additional contributions. One example is the ionic potential $\sum_{j}s_{j}^{2}$ Zagorodny2004 . Many of our results are valid in greater generality, but we explore the simplest case (2).

The LLL dynamics has two locally conserved fields, namely the $z$-component of the spin, $s_{j}$, and the energy

[TABLE]

where for convenience we have introduced the phase difference $r_{j}=\phi_{j+1}-\phi_{j}$. From the equations of motion (2) one deduces the form of the spin current, $\mathcal{J}_{j}^{\mathsf{s}}$, and energy current, $\mathcal{J}_{j}^{\mathsf{e}}$, as

[TABLE]

High-temperature diffusive regime — We first assume that there are no further conservation laws. The equilibrium Gibbs measures are then given by the two-parameter family

[TABLE]

for a chain of length $N$ spins with periodic boundary conditions and partition function $Z_{N}(\beta,h)$. Infinite volume equilibrium averages will be denoted by $\langle\cdot\rangle_{\beta,h}$. We note that the Hamiltonian is even and the currents are odd in $r_{j}$. Hence

[TABLE]

But then also their derivatives with respect to $\beta,h$ vanish and, using the definition of the Drude weight given in MendlSpohnCurrent , one concludes that both Drude weights are zero. Thus the conventional expectation is to have a diffusive spreading of the equilibrium time-correlations for spin and energy. Since they have opposite signature under time reversal, the cross diffusion coefficient should vanish.

To confirm, we performed molecular dynamics simulations at inverse temperature $\beta=1$. In Fig. (1) we show numerical results for the spin and energy correlations defined by

[TABLE]

where $\langle\dots\rangle_{\beta,h}^{\mathrm{c}}$ denotes the connected correlation defined as $\langle Q_{j}(t)Q_{0}(0)\rangle^{\mathrm{c}}_{\rm eq}:=\langle(Q_{j}(t)-\langle Q_{0}\rangle_{\rm{eq}})(Q_{0}(0)-\langle Q_{0}\rangle_{\rm{eq}})\rangle_{\rm{eq}}$.

These simulations are for system size $N=512$ and were run up to time $t=1024$ [simulation details are given later in Sec. 4]. We see that spin and energy autocorrelations indeed show diffusive behavior, while there are no cross correlations.

The novelty of our contribution is to establish that at lower temperatures the dynamical properties change dramatically through the appearance of ballistic sound propagation. But before embarking on that discussion we use the opportunity to illustrate that equilibrium time-correlations for an integrable spin-chain are dominated by a broad ballistic spreading, in contrast to the sub-ballistic broadening of the two sound peaks in the case of a nonintegrable chain. The form of dynamical correlations in classical integrable models have been discussed in several earlier work including zhao06 ; PhysRevE.94.062130 for the Toda chain and in prosen2013macroscopic for a particular integrable spin chain, which we will now discuss.

Integrable LLL model — Faddeev and Takhtajan faddeev2007hamiltonian discovered an integrable version of the LLL model, which still has nearest neighbor coupling but is no longer quadratic. Their Hamiltonian is given by

[TABLE]

with local energy

[TABLE]

with $\rho\geq 0$. The hamiltonian (2.12) seems to be the only known integrable classical spin chain. Easy plane corresponds to $\rho>0$, while in the limit $\rho\to 0$ one recovers the isotropic interaction

[TABLE]

The infinitely extended Faddeev-Takhtajan spin chain has a countable number of locally conserved fields, which are constructed by successive differentiations of the $R$-matrix, see faddeev2007hamiltonian .

Here we focus on the $z$-component of the magnetization, $s_{j}$, and the local energy, $e_{j}$, as the first items in the list of conservation laws. We display their time-correlations at $\rho=1$ with inverse temperature $\beta=1$ and magnetic field $h=0$, see Fig. (2). In the scaling plot Fig. (3), we see that the energy and spin correlations show good ballistic scaling already at short times. In prosen2013macroscopic simulations of the spin current correlations are reported at parameters $\rho=1,\beta=0.25,h=0$.

Without losing integrability, the Hamiltonian (2) can be analytically continued to purely imaginary $\rho$, which amounts to replace the trigonometric functions by their hyperbolic cousin faddeev2007hamiltonian . Physically this corresponds to easy-axis regime. We refer to the discussions in prosen2013macroscopic ; AvijitFT , also reporting on parameters with diffusive spreading.

3 Easy plane at low-temperatures

Low-temperature effective hamiltonian — We now return to the nonintegrable XXZ spin chain. As the temperature is lowered, if we impose easy-axis anisotropy, $\Delta>1$, the spreading of spin and energy correlations is still diffusive, as confirmed by further molecular dynamics simulations. In Fig. (4) we show the easy axis data. Apparently diffusion still holds at $\Delta=1$, although with slow convergence AvijitChaos . However, for easy-plane, $\Delta<1$, phase differences become locally almost conserved, which drastically changes the dynamical behavior, as will now be explained in detail. The equations of motion are given by (2).

Phase differences are defined through

[TABLE]

where we choose coordinates such that $r_{j}\in[-\pi,\pi]$. $r_{j}=0$ corresponds to $\phi_{j}=\phi_{j+1}$ which is the minimum of the cosine-potential. The dynamics of the $r_{j}(t)$’s has the following generic structure. One starts from their extended version with $\tilde{r}_{j}\in\mathbb{R}$ governed by

[TABLE]

for a given collection of smooth functions $\{g_{j}(t),j\in\mathbb{Z}\}$. The restriction to the unit circle $S^{1}$ is achieved by setting

[TABLE]

with $n_{j}(t)$ the integer winding number for bond $j$; we may choose $n_{j}(t=0)=0$, while $r_{j}(t)$ is a smooth function on $S^{1}$. By construction the $\tilde{r}_{j}$’s have the form of a local conservation law, which in integrated version reads

[TABLE]

for a system with periodic boundary conditions. On the other hand,

[TABLE]

Hence $r_{j}(t)$ is locally conserved only if the total winding number remains constant: $\sum_{j}n_{j}(t)=0$. The dynamical events $r_{j}(t)=\pm\pi$ are the “phase slip” processes where the winding number changes. Thus the field of phase differences $r_{j}$ is locally conserved only until the first phase slip event. However, as we show in Fig. (5), at low temperatures most phase slip events come in closely-spaced pairs with no change in the total winding number, so the coarse-grained dynamics actually respects this conservation law until one has unpaired phase slip events.

As illustrated in Fig. (5a) and (5b), for high temperatures there are lots of phase slips. However for easy-plane and low enough temperatures the $z$-component is approximately constant and the phase difference is trapped by the cosine-potential. Then phase slips are very much suppressed, see Fig. (5c) and (5d). Therefore, there is an emergence of an approximately conserved quantity. In other words, phase slips are thermally activated process which, in a low-temperature regime, can safely be ignored on the time scales reached by our simulation and, in fact, much longer. Phase differences are thus approximately conserved in this regime.

Let us first attempt a rough estimate for the presence of a third conservation law. We assume constant $r_{j}=r$ and $s_{j}=s$. Adding also an external field, $h$, the energy of this configuration equals

[TABLE]

Its minimum is located at $r=0,s=h/[2(1-\Delta)]$ and we require $|h|<2(1-\Delta)$ to ensure that the minimum lies inside $\{|s|<1\}$. Let us compute the energy required for a phase slip event caused by motion of a single spin. For the case where $s$ remains fixed but $r$ changes to $\pi$, the energy barrier $E_{\rm slip}$ can be easily computed and leads to

[TABLE]

On the other hand for the case where the angle $r$ remains fixed and the spin moves to the north pole (given by $s=1$), the energy barrier is

[TABLE]

We expect phase slips to occur at a rate $\sim e^{-\beta\Delta E_{\rm slip}}$ where $\beta$ is the inverse temperature, and so $\beta E_{\rm slip}>1$ could be a rough criterion for small number of phase-slips and a new approximate conservation law. In Fig. (5) we show space-time plots showing phase-slip events seen in simulations at (a) high temperature $(\beta=1)$, zero magnetic field, (b) high temperature $(\beta=1)$, finite magnetic field $(h=0.3)$, (c) low temperature $(\beta=5)$, zero magnetic field and (d) low temperature $(\beta=5)$, finite magnetic field $(h=0.3)$ with parameters $\Delta=1/2,N=64$. In Fig. (6), we show the dependence of average total number of observed phase slips on the inverse temperature $\beta$ up to $t=512$ seen in the low temperature simulations (with parameters $\Delta=1/2$, $h=0$). We see that $E_{\rm slip}$ takes the values $4$ and $2$ from the two energy estimates mentioned above. From the simulations we find the expected activated dependence form, however with a barrier $\Delta E\approx 1.73$. To understand such behavior we have to investigate in more detail the actual process of phase slips. Let us consider the spin chain with zero external magnetic field and in the easy-plane regime $(h=0,\Delta<1)$. In the ground state all spins are aligned and lying on the $xy$ plane. We are interested in finding the minimum-energy spin configuration in which there is a phase slip. What we find at the lowest temperatures where such phase slip events do happen in our simulations is that the phase slip event is centered on a single spin $j$ that moves out of the $xy$ plane and is very close to either $\theta_{j}=0$ or $\theta_{j}=\pi$. The nearby spins $i$ move to near the configuration of their angles $\theta_{i}$ that minimizes the energy given the special orientation of the one “central” spin $j$, with all spins except the central one oriented near the same $\phi_{i}$. When we do this minimization for $\Delta=0.5$ and $h=0$, the resulting minimum energy of the phase slip event is $E_{\rm slip}\cong 1.73$ ongoing . As expected, the measured density of phase slip events shows a thermally activated dependence on temperature as $\sim\exp(-\beta E_{\rm slip})$, as shown in Fig. (6). These phase slip events are produced at this low density by the chaotic equilibrium dynamics of the LLL chain.

To work with an almost conserved field is somewhat vague and it is more convenient to modify the dynamics such that $r_{j}(t)$ is strictly locally conserved. Of course, this is a valid approximation only in a regime with a very low density of phase slips. We rewrite our Hamiltonian in a slightly more general form as

[TABLE]

If one would set $U(x)=-\cos(x)$, included to have motion on $S^{1}$, then $H_{\mathrm{lt}}=H$. To suppress phase slips entirely we modify $U(x)$ by adding infinitely high potential barriers at $x=\pm\pi$. Then up to the first phase slip event, the dynamics generated by $H$ agrees with the dynamics generated by $H_{\mathrm{lt}}$, but they differ later on. Actually, for what we want to show the precise shape of the potential barriers plays no role, as long as phase slips are forbidden.

Nonlinear fluctuating hydrodynamics (NFH) in the low temperature regime — As before, the hamiltonian equations of motions for (3.9) are

[TABLE]

The conserved fields are phase difference, spin, and energy,

[TABLE]

where for later convenience we introduced the shorthands $f_{j}=f(s_{j})$ and $U_{j}=U(r_{j})$. From the equations of motion we obtain for the current of the phase difference,

[TABLE]

for the spin current,

[TABLE]

and for the energy current,

[TABLE]

The grand canonical ensemble of $H_{\mathrm{lt}}$ for a finite system with $N$ lattice sites and periodic boundary conditions is given by

[TABLE]

with the normalizing partition function

[TABLE]

where $\nu$ is the “chemical potential” for the additional conserved field $r_{j}$. Infinite volume averages with respect to (3.15) are denoted by $\langle\cdot\rangle_{\nu,h,\beta}$. The canonical free energy is defined as

[TABLE]

The infinite volume equilibrium averages of $r_{j}$, $s_{j}$, $e_{j}$ are

[TABLE]

independent of $j$ because of translation invariance. By convexity of $F$, these relations define the inverse mapping $(\mathsf{r},\mathsf{s},\mathsf{e})\mapsto(\nu(\mathsf{r},\mathsf{s},\mathsf{e}),h(\mathsf{r},\mathsf{s},\mathsf{e}),\beta(\mathsf{r},\mathsf{s},\mathsf{e}))$.

Under the constraints (3.8) the LLL equilibrium time-correlations of $s_{j},e_{j}$ should be well approximated by the same time-correlations as computed from the dynamics governed by $H_{\mathrm{lt}}$. But $H_{\mathrm{lt}}$ is just one particular anharmonic chain and, as explained in MendlSpohnNLS2015 , the time-correlations can be predicted from nonlinear fluctuating hydrodynamics. We do not repeat here the details, but merely point out that in normal mode representation one arrives at a three-component fluctuating field, \vec{\phi}(x,t)=\big{(}\phi_{-1}(x,t),\phi_{0}(x,t),\phi_{1}(x,t)\big{)}. The Euler currents have to be expanded to second order, which in approximation then leads to the coupled Langevin equations

[TABLE]

Here $D$ is a constant diffusion matrix and $B\vec{\xi}(x,t)$ is Gaussian white noise, both related through fluctuation-dissipation as $BB^{T}=2D$. This part of the equation is a phenomenological ansatz for the effective noise and dissipation produced by the deterministic chaos. However the sound speed, $c$, and the three symmetric coupling matrices $\vec{G}$ have to be computed from the underlying microscopic model. In particular $\vec{G}$ determines the dynamical universality class. Fortunately, for the LLL chain the magic identity (proven below)

[TABLE]

is available. Using this property the precise form of $\vec{G}$ and its relation to second derivatives of the free energy have been established in MendlSpohnNLS2015 .

We turn to the proof of (3.20).

(i) For the phase difference current we obtain

[TABLE]

(ii) Correspondingly, for the spin current we obtain

[TABLE]

where we used partial integration in the last step.

(iii) For the energy current there are more terms to be considered,

[TABLE]

We conclude that the $H_{\mathrm{lt}}$-currents satisfy (3.20) independently of $U$, provided $U$ diverges sufficiently fast at the two border points so as to have the boundary terms vanish.

When applying the magic identity, the potential $U$ is specified as the cosine with infinitely high barriers added at $x=\pm\pi$. This does not change the free energy. Thus at the end the coupling matrix $\vec{G}$ is given in terms of the LLL free energy. We refer to Appendix A for more details. The chemical potential $\nu$ is not a physical control parameter that we will explore. It is needed in the second order expansion, but will be set to $\nu=0$ afterwards.

4 Numerical results at low temperature

We now present numerical results for the dynamical correlations for two cases with (I) finite magnetic field and (II) zero field, at inverse temperature $\beta=10$. Apart from the spin and energy correlations defined in Eq. (4.1) we now measure

[TABLE]

corresponding to the extra conserved variable. We also compute the normal mode correlations $C_{\sigma,\sigma^{\prime}}$, where $\sigma,\sigma^{\prime}=+1,0,-1$ correspond to the modes $\phi_{+1},\phi_{0},\phi_{-1}$ respectively.

Simulation details: In the simulations we first generated equilibrium configurations corresponding to the canonical distribution, specified by $\beta,h$, using a Metropolis Monte Carlo algorithm. The equilibrated configuration was then evolved with our Hamiltonian dynamics using either a fixed step size and sometimes an adaptive step size 4th order Runge-Kutta (RK4) integrator. This maintains the global conservations (total energy and magnetization) up to an absolute deviation of $\sim 10^{-6}$ and $\sim 10^{-13}$ respectively until the final time of evolution. Individual spin lengths are also conserved up to $\sim 10^{-13}$. Averages were compute over around $10^{6}$ initial conditions.

In all our simulations we fixed $\beta=10$. At much smaller temperatures, phase slips would be even rarer but the dynamics of small fluctuations about the ordered state is expected to be closer to integrable, and the time required to see the asymptotic scaling becomes inaccessible.

Case-I $(\nu=0,~{}h=0.3)$: In this regime we expect to confirm the KPZ sound modes and Levy heat modes as predicted by NFH. We have chosen the following parameters: $\Delta=0.5,~{}\beta=10.0,~{}N=8192$. The $z$-component of the total spin, energy, and individual spin length are conserved up to $10^{-15}$ and $10^{-5}$ and $10^{-10}$ respectively. Phase Slip processes are rare but not completely absent. In Fig. (7), the correlations $C_{\alpha\beta}(x,t)$ are plotted at different times.

In Fig. (8), we plot the sound and heat modes obtained after normal mode transformations and with different scalings. Values of the $G$-matrices are given in Table-1.

As for anharmonic chains, $G^{0}_{00}=0$ which implies that the self-coupling term is absent for the heat mode. Also $G^{\sigma}_{\sigma\sigma}~{}(\sigma=\pm 1)$ are non-zero which in the context of NFH is the crucial property for the KPZ scaling of sound modes and Levy-5/3 scaling of heat mode:

[TABLE]

where $f^{\alpha}_{\text{Levy}}(x)=\text{InverseFourier}\left[e^{-|k|^{\alpha}}\right]\sim\dfrac{1}{\pi}\dfrac{1}{|x|^{\alpha+1}}$ and the universal $f_{\text{KPZ}}$ function is tabulated in prahofer2004exact .

From the results of our simulations, shown in Fig. (8), we find that the sound mode is better described by a KPZ-type scaling while the heat mode appears to be closer to diffusive than the expected Levy-$5/3$. The sound speed from the simulations is $c\approx 0.865$ which is close to the theoretical estimate $c=0.85217$.

Case-II $(\nu=0,~{}h=0)$: In this regime, the $r$- interaction potential is symmetric under reflection. This leads to a distinct universality class. On the basis of NFH the sound modes are expected to be diffusive and the heat mode to be Levy-3/2. We have chosen the following parameters: $\Delta=0.5$, $\beta=10.0,~{}N=8192$. Phase Slip events are absent in this regime. The $z$-component of the spin, energy, and individual spin length are conserved up to $10^{-12}$ and $10^{-5}$ and $10^{-10}$ respectively. In Fig. (9), the correlations $C_{\alpha\beta}(x,t)$ are plotted at different times. In Figs. (10), we plot the sound and heat modes obtained after normal mode transformations and with different scalings. The $G$-matrices are given in Table-2.

We have $G^{0}_{00}=0$ here as well. Unlike the previous case $G^{\sigma}_{\sigma\sigma}~{}(\sigma=\pm 1)$ are zero which gives rise to diffusive sound mode and Levy-3/2 heat mode in NFH:

[TABLE]

As mentioned earlier the expectation from theory is that the sound modes are diffusive while the heat mode is Levy-$3/2$. However, as can be seen from our simulations, i.e. Fig. (10), these scalings are not seen conclusively. The sound speed from the simulations is $c\approx 0.931$ which agrees well with the theoretical estimate $c=0.92837$.

5 Discrete time LLL dynamics

In this section we discuss a discrete version of the lattice Landau-Lifshitz equations, using an integration algorithm that explicitly preserves the total angular momentum and the total energy of the system to machine precision independent of the time step used in the numerical integration, while allowing a bounded error in the fixed-length constraint on each classical spin (in the anisotropic case with $\Delta\neq 1$, it is of course only the $z$ component of the angular momentum which is thus conserved). As will be apparent from the description below, the smallness of the violation of the fixed-length constraint for each spin is controlled by the size of the time-step used in the numerical integration. Additionally, this error does not grow with the total time over which the system is evolved. Since all our general arguments and theoretical analysis rely heavily on the existence of a conserved energy and angular momentum density, and the fixed-length constraint on each individual classical spin does not play a similar central role in the theoretical analysis, our procedure allows us to use relatively large time-steps while preserving the universality class of the dynamics to machine precision.

Our procedure may be viewed as a modification of the so-called “odd-even dynamics” that has been employed previously in the literature frank1997geometric . In the odd-even decomposition, one splits the time-evolution into two steps, one of which sets in motion all the odd spins, allowing them to precess in the exchange field supplied by the (temporarily) static even spins, while the other step reverses the role of even and odd spins to evolve the even spins. Below we refer to these individual steps as odd and even precessions.

Our modification is suggested by a particular approximation (the Cayley approximation) to the evolution operators that implement these even or odd precessions over a small time $\delta$ in a way that preserves the orthogonal nature of the evolution operator. Consider Eqn. (2.2) with $j$ odd and all $B_{j}$ held fixed (by keeping all even spins fixed). With all $B_{j}$ fixed, the right hand side is a $B_{j}$ dependent rotation of all the odd spins. Let us denote the corresponding linear operator by $R_{o}(B_{e})$ (where $o$ stands for odd and the subscript $e$ on $B$ reminds us that $B$ depends on the current configuration of all even spins, held fixed for the duration of this step). Here $R_{o}$ and $R_{e}$ are antisymmetric matrices that serve as the generators of the corresponding rotations. The odd spins after a time $\delta$ can be obtained by applying the operator $\exp(\delta R_{o})$ to the initial configuration of odd spins.

We begin with the Cayley approximation to individual precessions,

[TABLE]

and similarly for $R_{e}$. With this in hand, we define the operator

[TABLE]

Here, we have explicitly displayed the dependence of $R_{o}$ ($R_{e}$) on $B_{e}$ ($B_{o}$) in the configuration on which the operator acts. Let us schematically denote the configuration of odd spins after the action of the first term by $S^{1/2}_{o}$. The second term is a first order approximation to the precession of even spins in the exchange field provided by $S^{1/2}_{o}$, starting with the initial configuration of even spins, denoted schematically by $S^{0}_{e}$. Now, the first term amounts to precessing odd spins backwards in time (again using the same first order approximation), using the exchange field provided by the configuration $S^{0}_{e}$, but starting with the configuration $S^{1/2}_{o}$ in order to end up finally with the initial configuration $S^{0}_{o}$ of the odd spins.

As a result, it is easy to see that $O_{1}$ explicitly preserves both the total angular momentum (only the $z$ component of the angular momentum if $\Delta\neq 1$) and the total energy of the system, independently of the value of $\delta$. However, $O_{1}$ is not orthogonal. $O_{2}$ defined as

[TABLE]

has the same properties.

Next we note that the product $O_{2}O_{1}$ provides a second order accurate (in $\delta$) approximation to the actual spin dynamics. This approximation preserves the total angular momentum and total energy independently of $\delta$, while preserving the fixed-length constraint on each spin only to second order. Additionally, our numerical tests reveal that the violation of the fixed length constraint does not grow with total time of integration if each time step is implemented by the operator $O_{2}O_{1}$, even for relatively large values of $\delta$.

Therefore, we integrate the lattice Landau Lifshitz equations using $O_{2}O_{1}$ to evolve the system over one time step. As noted earlier, this can be viewed as a modification of the standard odd-even dynamics used earliear.

We display our results from discrete time dynamics with $dt=1$.

Case-III $(\beta=8.0,h=0.3,N=8192):$ In Fig. (11), we show plots for evolution of the correlations $C_{ss}(x,t)$ and $C_{ee}(x,t)$, while in Figs. (12,13) we plot the sound mode and heat modes after appropriate translation and with different scalings. As with the continuous time dynamics, in Fig. (8), we find again that the sound mode is better described by a KPZ-type scaling while the heat mode appears to be closer to diffusive than the expected Levy-$5/3$. The sound speed from the simulations is $c\approx 0.816$ which is close to the theoretical estimate $c=0.809$.

Case-IV $(\beta=8.0,h=0.0,N=8192):$ In Fig. (14), we show plots for evolution of the correlations $C_{ss}(x,t)$ and $C_{ee}(x,t)$, while in Figs. (15,16) we plot the sound mode and heat modes after appropriate translation and with different scalings. In this case the expectation from theory is that the sound modes are diffusive while the heat mode is Levy-$3/2$. However, as can be seen in Figs. (15,16), these scalings are not seen conclusively in the simulations. The sound speed from the simulations is $c\approx 0.831$. Note that, in contrast to the continuous time model, we do not have a simple construction of the stationary equilibrium state. Hence even for the sound speed we do not have a reliable prediction.

6 Conclusions

In this work we have mainly studied the classical XXZ or LLL model and an integrable counterpart of it in easy-plane regimes. We are mostly concerned with the basic conserved quantities — energy and magnetization in both models (integrable and non-integrable) at different regimes of model parameters and external parameters like temperature and magnetic field. For the LLL model, being non-integrable, these two are the only known conserved quantities but at low temperature a third almost conserved quantity emerges, producing sound modes. We have put this into the framework of NFH to get some predictions about the dynamical scaling. We have shown that the non-integrable XXZ chain in the easy-plane regime displays diffusive behaviour at high temperatures, and sound modes with KPZ broadening at low temperatures. The integrable spin chain however demonstrates ballistic behaviour in the easy plane regime. We also find that at high temperature the easy-axis nonintegrable case gives us perfectly diffusive behavior. Our various results are summarized in Table 3. For the heat mode, while NFH predicts anomalous Levy broadening, we have not been able to clearly observe this yet in the simulations. It is possible that one requires larger system sizes and longer times to observe this and further studies are necessary. Similar problems related to non-observation of the expected asymptotic behaviour, due to possible finite size or finite time effects, have been reported in other recent transport studiesZhong12 ; Das2014 ; Wang13 ; Dhar2018 ; Cintio18 , and understanding their origins is an open problem.

7 Acknowledgements

AD would like to thank the support from the grant EDNHS ANR-14-CE25-0011 of the French National Research Agency (ANR) and from Indo-French Centre for the Promotion of Advanced Research (IFCPAR) under project 5604-2. MK gratefully acknowledges the Ramanujan Fellowship SB/S2/RJN-114/2016 from the Science and Engineering Research Board (SERB), Department of Science and Technology, Government of India. MK also acknowledges support the Early Career Research Award, ECR/2018/002085 from the Science and Engineering Research Board (SERB), Department of Science and Technology, Government of India. MK would like to acknowledge support from the project 6004-1 of the Indo-French Centre for the Promotion of Advanced Research (IFCPAR). DAH was supported in part by (USA) DOE grant DE-SC0016244. The numerical simulations were done on Mowgli, Mario and Tetris clusters of ICTS-TIFR and Gaggle and Pride clusters of DTP-TIFR.

Appendix A Coupling coefficients for nonlinear fluctuating hydrodynamics

Given the magic identity (3.20), one can follow MendlSpohnNLS2015 to obtain the couplings of the quadratic nonlinearities of NFH. Our simulations are for $\nu=0$. Restricting to this case, the averages below, denoted by $\langle\cdot\rangle$, $\langle\cdot;\cdot\rangle$, $\langle\!\langle\cdot;\cdot\rangle\!\rangle$, refer to fixed $h,\beta$ at $\nu=0$. One obtains

[TABLE]

and

[TABLE]

Here $c_{\mathrm{s}}$ denotes the speed of sound determined through

[TABLE]

where $\Gamma=\langle\!\langle s_{0};s_{0}\rangle\!\rangle\langle\!\langle e_{0};e_{0}\rangle\!\rangle-\langle\!\langle s_{0};e_{0}\rangle\!\rangle^{2}$ and we have introduced the shorthand notation

[TABLE]

with $\langle f;g\rangle=\langle fg\rangle-\langle f\rangle\langle g\rangle$ denoting the second cumulant. Also T denotes the transpose relative to the anti-diagonal and

[TABLE]

The thermodynamic averages and cumulants can be obtained as appropriate derivatives of the free energy $F(h,\nu,\beta)$ defined in (3.17) with $\nu$-derivative at $\nu=0$,

[TABLE]

As second order Taylor coefficients, the $G$-matrices are symmetric.

The free energy has to be numerically evaluated. An efficient method is to use transfer operator techniques MendlSpohnNLS2015 . Inserting our simulation parameters yields (rounded to 4 digits):

$\Delta=0.5$, $\beta=10$, $h=0.3$

, speed of sound $c_{\mathrm{s}}=0.85217$,

[TABLE]

$\Delta=0.5$, $\beta=10$, $h=0$

, speed of sound $c_{\mathrm{s}}=0.92837$,

[TABLE]

In TABLE I and II these results are compared with data coming from molecular dynamics.

Appendix B Low temperature approximation

It is instructive to work out the prefactors of $G$-matrices (A.1), (A.2) in the harmonic approximation. In principle, also the next order correction could be computed, compare with MendlSpohnNLS2015 . We start from the hamiltonian

[TABLE]

with $s_{j}\in[-1,1]$, phase difference $r_{j}=\phi_{j+1}-\phi_{j}$, and external field $h$.

The minimum of $H$ is assumed for $r_{j}=0$ and $s_{j}=\frac{h}{2(1-\Delta)}$ for all $j$. We expand $H$ in a Taylor series around this minimum and set $t_{j}=s_{j}-\frac{h}{2(1-\Delta)}$. Neglecting boundary terms one obtains

[TABLE]

with the ground state energy $e_{\mathrm{g}}=-1-\frac{h^{2}}{4(1-\Delta)}$, the quadratic contribution

[TABLE]

and the cubic correction

[TABLE]

Here we introduced the shorthands

[TABLE]

We write $H_{0}=\frac{1}{2}\langle x,Ax\rangle$ with $x=(t_{-N/2},\dots,t_{N/2-1},r_{-N/2},\dots,r_{N/2-1})\in\mathbb{R}^{2N}$ and a block-diagonal matrix $A$ consisting of two blocks $A_{t}$ and $A_{r}=c\mathbb{I}_{N}$, with tridiagonal

[TABLE]

The partition function of the quadratic part plus $Ne_{\mathrm{g}}$ is defined as

[TABLE]

For large $\beta$, the integration domain of $s_{j}$ and $r_{j}$ can be approximately extended to $\mathbb{R}$, which yields

[TABLE]

The determinant of the tridiagonal matrix $A_{t}$ is available in closed form HuOConnell1996 , such that

[TABLE]

Inserting these relations results in the free energy

[TABLE]

In the thermodynamic limit $N\to\infty$,

[TABLE]

With this input we work out the free energy derivatives (A) to leading order in $1/\beta$ and obtain

[TABLE]

Therefore the speed of sound is given by

[TABLE]

For $h=0$ and $\Delta=0.5$, we get $c_{s}=1$.

For the $G$ matrices (at $h=0$) we arrive at

[TABLE]

which for our values of parameter ($\Delta=0.5$) gives,

[TABLE]