Anomalous relaxation and the high-temperature structure factor of XXZ spin chains

TL;DR

This paper investigates the high-temperature spin structure factor of XXZ spin chains, revealing anomalous decay behaviors due to bound states and diffusive effects, verified through matrix product operator calculations.

Contribution

It introduces a comprehensive analysis of the high-temperature structure factor in XXZ chains, including diffusive corrections and the role of magnon bound states, with detailed numerical verification.

Findings

Hierarchy of timescales in the gapped phase

Anomalous power-law decay of the structure factor

Crossover from diffusion to superdiffusion at the isotropic point

Abstract

We compute the spin structure factor of XXZ spin chains in the Heisenberg and gapped (Ising) regimes in the high-temperature limit for nonzero magnetization, within the framework of generalized hydrodynamics including diffusive corrections. The structure factor shows a hierarchy of timescales in the gapped phase, owing to $s$-spin magnon bound states (`strings') of various sizes. Although short strings move ballistically, long strings move primarily diffusively as a result of their collisions with short strings. The interplay between these effects gives rise to anomalous power-law decay of the spin structure factor, with continuously varying exponents, at any fixed separation in the late-time limit. We elucidate the crossover to diffusion (in the gapped phase) and to superdiffusion (at the isotropic point) in the half-filling limit. We verify our results via extensive matrix product operator calculations.

Equations

$$H={\sum }_i\left[S_i^x{S}_{i+1}^x+S_i^y{S}_{i+1}^y+\Delta S_i^z{S}_{i+1}^z\right],$$

$${ϵ}_s(q)=k_s{\Delta }^{1-s}\sin (2q),$$

$$C(x,t)\simeq {\sum }_{s\ge 1}{s}^2{f}^s[J_x(k_s{\Delta }^{1-s}t){]}^2,$$

$$C(x,t)\approx {\sum }_{s=1}^{s_{\ast }}{s}^2\frac{f}{t}{\left(f\Delta \right)}^{s-1},$$

$$C(x,t)\sim \frac{f}{t}{\left(\frac{t}{x}\right)}^{1-\frac{|\log f|}{\log \Delta }}{\log }^2\left(\frac{t}{x}\right)\sim t^{-\frac{2|\mu |}{\eta }}{\log }^{2}t,$$

$$s\gtrsim s_0(t)\sim 1+\frac{|\log (ft)|}{2\log \Delta }.$$

$$C(x\lesssim x_0,t)\sim t^{-\frac{1}{2}-\frac{|\mu |}{\eta }}{\log }^{2}t,2\left|\mu \right|<\eta .$$

$$C(x,t)=\sum _{s=1}^{\infty }\int du{\rho }_s^{\mathrm{tot}}(u){\theta }_s(1-{\theta }_s)(m_s^{\mathrm{dr}}{)}^2{\phi }_t[x-v_s(u)t],$$

$${\int }_{-\infty }^{\infty }dxC(x,t)=\frac{1}{4}(1-{\tanh }^2\mu ),$$

$$D=\frac{2\sinh \eta }{9\pi }\sum _{s=1}^{\infty }(1+s)\left[\frac{s+2}{\sinh \eta s}-\frac{s}{\sinh \eta (s+2)}\right],$$

$$C(x,t)\sim \frac{{\mu }^2}{\eta }\frac{\log (t/x)}{x},e^{-\eta /|\mu |}\ll \frac{x}{t}\ll 1,$$

$$C(x,t)=\frac{1}{4}{\phi }_t(x)=\frac{1}{8\sqrt{\pi D(\eta )t}}{\mathrm{e}}^{-\frac{x^2}{4D(\eta )t}},$$

$$C(x,t)={\mu }^2\left[C_{anom.}(x{\mu }^2,t{\mu }^3)+\frac{1}{t}{C}_{reg.}(x/t)\right],$$

$$C(0,t)=t^{-2/3}g(\mu t^{1/3}),\quad g(y)=\left\{\begin{array}[]{lr}\mathrm{const.}&y\ll 1\\ 1/y&y\gg 1\end{array}\right.$$

$$C(x/t)\sim {\mu }^{2}t/x^2,\mu \ll (x/t)\ll 1.$$

Full text

Many experimentally relevant one-dimensional systems are described by approximately integrable models, such as the Hubbard, Heisenberg, and Lieb-Liniger models Kinoshita et al. (2006); Hild et al. (2014); Rigol et al. (2008). The nonequilibrium dynamics of integrable systems, their failure to thermalize, and their possession of an extensive set of conservation laws, have been explored extensively Caux and Essler (2013); Ilievski et al. (2016, 2015); Pozsgay et al. (2014). (In experiments, integrability is approximate, and gives rise to “prethermal” intermediate-time regimes of effectively integrable dynamics Gring et al. (2012); Langen et al. (2016); Tang et al. (2018).) Integrable systems support stable, ballistically propagating quasiparticles even at high temperature. In the simplest cases (e.g., free fermions), these particles carry the same quantum numbers as the microscopic degrees of freedom, and move with a velocity set by the band structure. In interacting integrable models, however, each quasiparticle is dressed by all the others Takahashi (1999). This dressing can lead to remarkable dynamical effects, for instance in the “gapped” phase of the XXZ model considered here: here, even though quasiparticles move ballistically, finite-temperature spin transport is diffusive in the absence of an external field Sachdev and Damle (1997); Damle and Sachdev (1998); Zotos (1999); Damle and Sachdev (2005); Sirker et al. (2009); Karrasch et al. (2014); Ljubotina et al. (2017); Piroli et al. (2017); Bertini and Piroli (2018); De Nardis et al. (2018a); Gopalakrishnan et al. (2018); De Nardis et al. (2018b); Gopalakrishnan and Vasseur (2019); De Nardis et al. (2019).

Recently, a coarse-grained approach to integrable dynamics has been developed; this approach is termed “generalized hydrodynamics” (GHD) Castro-Alvaredo et al. (2016); Bertini et al. (2016), see also Doyon and Yoshimura (2017); Doyon and Spohn (2017); Ilievski and De Nardis (2017a); Bulchandani et al. (2018); Doyon et al. (2018); Piroli et al. (2017); Ilievski and De Nardis (2017b); Collura et al. (2018); Alba and Calabrese (2017); De Luca et al. (2017); Bertini and Piroli (2018); Bulchandani et al. (2017); Bulchandani (2017); Doyon et al. (2017); Caux et al. (2017); Cao et al. (2018); Doyon (2018); Schemmer et al. (2019); Alba et al. (2019). A core insight of GHD is that an integrable system can be mapped to an appropriate classical soliton gas Doyon et al. (2018). Assuming the system is initially in local equilibrium, the velocities of these solitons can be computed using the thermodynamic Bethe ansatz Bonnes et al. (2014); Castro-Alvaredo et al. (2016); Bertini et al. (2016), which is much more tractable than exactly simulating the full dynamics. In the initial formulation of GHD, the dressing of quasiparticles by interactions was treated at the “Euler” level, yielding purely ballistic hydrodynamics; recently, adding Gaussian fluctuations on top of this treatment was shown to give diffusive corrections to hydrodynamics De Nardis et al. (2018a); Gopalakrishnan et al. (2018); De Nardis et al. (2018b); Gopalakrishnan and Vasseur (2019). In the generic case, diffusive corrections occur on top of ballistic transport; however, in many situations the ballistic term is absent, and transport is dominated by normal or anomalous diffusion.

In this work, we show that even when ballistic transport is present, certain aspects of the structure factor exhibit anomalous exponents. We focus on the XXZ spin chain:

[TABLE]

where $S^{\alpha}=\frac{\sigma^{\alpha}}{2}$ are spin-$\frac{1}{2}$ operators. We are concerned with the case of easy-axis anisotropy $|\Delta|\geq 1$ at infinite temperature (so the sign of $\Delta$ is irrelevant). We define $\eta\equiv\cosh^{-1}\Delta$. The model has a conserved magnetization, $\sigma^{z}_{\rm tot}=\sum\nolimits_{i}\sigma^{z}_{i}$; we denote the associated magnetization density as $h=\tanh\mu$, corresponding to a filling $f=(1+h)/2$. At half-filling ($h=0$), ballistic transport is absent because the propagating quasiparticles carry no spin. One can easily see this for magnons in the ferromagnetic phase at low but nonzero temperatures: a magnon propagates as a $\downarrow$ spin through a domain of $\uparrow$ spins, then continues as an $\uparrow$ spin through a domain of $\downarrow$ spins, etc., so on average it does not carry magnetization. This result, which holds generally, was first noted in the low-temperature limit Damle and Sachdev (1998, 2005), and has recently been incorporated into the GHD framework. Since ballistic transport is absent, the dominant transport mechanism is diffusive. The diffusion constant $D$ has been rigorously bounded Ilievski et al. (2018) and computed using GHD De Nardis et al. (2018b); Gopalakrishnan and Vasseur (2019); De Nardis et al. (2019); $D$ diverges in the isotropic ($\Delta=1$) limit, at which superdiffusion takes place, with a time-dependent diffusion constant $D(t)\sim t^{1/3}$; wavepackets appear to spread with a non-Gaussian front that corresponds to the Kardar-Parisi-Zhang (KPZ) universality class Ljubotina et al. (2017, 2019); Kardar et al. (1986); Quastel and Spohn (2015). Away from half-filling, the $\uparrow$ and $\downarrow$ domains do not precisely cancel, so magnons do carry magnetization and ballistic spin transport is present. Beyond these results, the behavior of the dynamical spin structure factor and the optical conductivity in this high-temperature limit are still poorly understood. (There has been extensive work on these quantities at zero temperature Caux and Maillet (2005); Pereira et al. (2006, 2007); Sirker et al. (2008); Klauser et al. (2011) and in the low-energy field theory regime Sachdev and Damle (1997); Damle and Sachdev (2005, 1998); Sirker et al. (2009, 2011), but as we will see the physics is qualitatively different at high temperature.)

The present work addresses these issues, computing the spin structure factor within GHD. We focus on the connected correlation function $C(x,t)\equiv\langle S^{z}_{i+x}(t)S^{z}_{i}(0)\rangle^{c}$ evaluated at infinite temperature with chemical potential $\mu$; everything we discuss will involve large $x,t$ but arbitrary ratios $x/t$. We find that, even away from half filling, the local behavior of autocorrelators (i.e., for $x/t\ll 1$, corresponding to the return probability) evolves with continuously varying exponents that depend on $\Delta$ and $h$ (Fig. 1). There is a phase transition in the $(\Delta,h)$ plane, between ballistic (i.e., $1/t$) and sub-ballistic (i.e., $1/t^{\gamma}$ for $1/2<\gamma<1$) asymptotic behavior. We compute the exponent $\gamma$ as a function of $(h,\Delta)$, and show that for $\Delta>1$ it universally approaches $1/2$ as $h\rightarrow 0$, recovering (and shedding light on) diffusion at half-filling. This coexistence of ballistic and anomalous behavior was recently demonstrated Agrawal et al. (2019) for disordered integrable spin chains Essler et al. (2018); here we show that the same effect occurs in clean strongly interacting systems. At the Heisenberg point, the phase boundary in the $(\Delta,h)$ plane intersects the ballistic-diffusive phase boundary at $h=0$, and in this sense the isotropic Heisenberg point at $h=0$ is a dynamical multicritical point. We write down a scaling form for the structure factor as one approaches this critical point at finite $h$.

Low-filling limit.—Our results have an elementary interpretation in the limit in which $f\ll 1$. Here, $f\sim e^{2\mu}$ with $\mu\to-\infty$. Nevertheless the system is still at infinite temperature. Further, for the present discussion we take $\Delta\gg 1$. Under these conditions we can calculate the structure factor by elementary methods; we only invoke integrability to claim that quasiparticles are in fact stable. In this limit, the quasiparticles are essentially “bare”: an $s$-string is a sequence of $s$ $\uparrow$ spins on top of a $\downarrow$ background. Since an $s$-string can only move at $s$th order in perturbation theory, its velocity is $v_{s}\sim\Delta^{1-s}$. Neglecting dressing, the $s$-strings have free-particle dispersions of the form

[TABLE]

where $k_{s}$ are constants of order unity. In our discussion of this limit, we take $f\sim{\rm e}^{2\mu}\to 0$ and $\Delta\sim{\rm e}^{\eta}\to\infty$, but allow the ratio $\left|\mu\right|/\eta$ to be of order unity. In the following, we will use GHD to generalize our results to arbitrary filling and $\Delta>1$.

If we ignore diffusive corrections, the model is in effect a gas of free $s$-strings, which occur with probability $f^{s}$. All dressing effects are suppressed by factors of $f$ with no compensating factors of $1/\Delta$, so we neglect them. Then a single-particle calculation yields the structure factor, as follows:

[TABLE]

where $J_{x}$ denote Bessel functions of the first kind Fukuhara et al. (2013). A nontrivial contribution arises if a string beginning at $(0,0)$ has propagated to $(x,t)$. To explore the asymptotics of Eq. (3) we approximate the Bessel function as a step function and ignore irrelevant prefactors, $[J_{x}(k_{s}\Delta^{1-s}t)]^{2}\sim\Theta(x-\Delta^{1-s}t)\Delta^{s-1}/t$. Fixing a point $x$, and counting only those $s$-strings that have reached $x$ by the time $t$, we get

[TABLE]

where $s_{*}=1+\log(t/x)/\log\Delta$. There are two cases. When $f\Delta<1$, higher-order strings are too rare to contribute to the correlation function, which is dominated by the $1/t$ tail of the 1-strings. When $f\Delta>1$, the dominant strings at position $x$ are the heaviest strings that have made it there; the sum in Eq. (4) is given by the term of order $s_{*}$. This then gives the asymptotics

[TABLE]

for $2|\mu|<\eta$. The exponent $\gamma=\frac{2\left|\mu\right|}{\eta}$ in Eq. (5) goes to unity as $\eta\gg|\mu|$, suggesting subdiffusion ($\gamma<\frac{1}{2}$) as $\Delta\rightarrow\infty$ at the Eulerian level. This asymptotic decay will occur through a series of jumps Piroli et al. (2017); we assume here that one coarse-grains over long enough time windows to average out these features.

Diffusive corrections.—The asymptotics (5) arises because—at the Eulerian level—long strings are assumed to be effectively stationary for exponentially long times. At finite $f$ this is not, in fact, the case: when a $q$-string and an $s>q$-string collide, the $s$-string picks up a displacement of $2q$ sites Ganahl et al. (2013); Vlijm et al. (2015). Thus, all strings undergo subleading diffusive motion. In the low-density limit it suffices to consider the diffusion of $s$-strings due to collisions with 1-strings, so the diffusion constant scales as $D\sim f$. Because the model is integrable this diffusion takes place in addition to the ballistic motion of $s$-strings. At intermediate times, an $s$-string moves diffusively; however, there is a crossover to ballistic motion at times such that $\sqrt{ft}<\Delta^{1-s}t$, i.e., for strings satisfying

[TABLE]

When $f\Delta<1$, heavy strings are sparse so it does not matter whether they diffuse. When $f\Delta>1$, Euler-scale results remain valid at distances $x\gtrsim x_{0}\equiv\sqrt{ft}$, but the behavior at distances $x\lesssim x_{0}$ is qualitatively modified. Rather than staying immobile, strings with $s>s_{0}(t)$ move diffusively. Thus the autocorrelator decays as

[TABLE]

The slowest that $C(x,t)$ can decay is as $1/\sqrt{t}$, i.e., a subdiffusive decay of the return probability does not occur in this model. Instead, the generic behavior is an anomalous decay with a continuously varying power law between $\frac{1}{2}$ and $1$.

GHD approach at generic filling.— The above argument is elementary but is restricted to the limit of large $\Delta$ and low filling. We now show that our main conclusion – the anomalous decay (7) of the local autocorrelation function – holds generally for all $h$ whenever $\Delta>1$. In the general case, spin transport can still be understood in terms of a hierarchy of strings, but their interactions are now important and their velocity and effective charge are dressed by the collisions with other quasiparticles. These issues can be addressed using GHD: since we are dealing with a linear response problem, we take advantage of the fact that the quasiparticles are in local thermal equilibrium, and evaluate the dressed quasiparticle dispersion and quasiparticle distribution function using data from the thermodynamic Bethe ansatz solution. Then the hydrodynamic expression for the structure factor takes the form Doyon and Spohn (2017); Ilievski and De Nardis (2017b):

[TABLE]

where $u$ parameterizes the rapidity of a quasiparticle; $m^{\mathrm{dr}}_{s}$ is the dressed magnetization of string $s$, $\rho^{\rm tot}(u)$ is the available density of states, $\theta_{s}$ is its occupation number (Fermi factor), and $v_{s}$ is its effective velocity. These quantities have closed-form expressions for generic $\mu$ at infinite temperature Takahashi (1999); Ilievski et al. (2018); sup . Finally, the function $\varphi_{t}(\zeta)$ is the propagator of a string with quantum numbers $(s,u)$ from $(0,0)$ to $(x,t)$. At the Euler level this propagator would simply be a Dirac delta function. In principle the full form of $C(x,t)$, including diffusive corrections and possible nonlinearities, could be ascertained from flea-gas simulations Doyon et al. (2018). Here we are interested in the asymptotic behavior of this quantity. We therefore include the dominant “diagonal” diffusive corrections by broadening $\varphi_{t}(\zeta)$ to a Gaussian with variance $2D_{s}(\eta,u)t$. The diagonal quasiparticle diffusion constant $D_{s}(\eta,u)$ was computed in Refs. De Nardis et al. (2018a); Gopalakrishnan et al. (2018), and can be evaluated numerically. We can check explicitly that our hydrodynamic form for the structure factor (8) is consistent with the exact sum rule

[TABLE]

since the function $\varphi_{t}(\zeta)$ is normalized to unity sup .

Anomalous local relaxation.—Equipped with this GHD expression, we first consider local relaxation, i.e., $C(x,t)$ at fixed large $x$ when $t\rightarrow\infty$. There are two contributions at time $t$, from light strings (whose motion is primarily ballistic) and from heavy strings (which undergo Brownian motion from collisions with light strings). Regardless of $\mu$, the velocities of very heavy strings scale as $v_{s}\simeq e^{-\eta s}$; also, at infinite temperature, their densities scale as $\rho_{s}(u)=\rho^{\rm tot}_{s}(u)\theta_{s}\sim e^{-2|\mu|s}$ De Nardis et al. (2018b). The dressed magnetization of the heavy strings, meanwhile, is the same as the bare magnetization $m^{\rm dr}_{s}\simeq s$. We see that the asymptotics of $v_{s}$ and $\rho_{s}$ are identical to the low-filling limit: for $\eta>2\left|\mu\right|$, the return probability is dominated by the diffusive strings, $s>s_{0}(t)$, where $s_{0}(t)\simeq\frac{1}{2\eta}\log t$. It follows that Eq. (7) applies for all $\mu$ and $\eta>1$. (Note, however, that away from the perturbative limit $\eta\neq\log\Delta$.)

For a fixed $\eta>0$ and $\mu<\eta/2$, this asymptotic scaling sets in on timescales $t\gtrsim e^{2\eta/|\mu|}$; at shorter times (see below) we expect a smaller apparent exponent, since the dominant strings at those times are not yet exponentially suppressed (Fig. 2). Reaching the asymptotic regime on accessible timescales is numerically challenging: to see many heavy strings at short times, we need $\eta\ll 1$, i.e., near the isotropic limit; however, working near the isotropic limit leads to transient superdiffusion Ljubotina et al. (2017); Ilievski et al. (2018) at short times $t\ll\eta^{-3}$ Gopalakrishnan and Vasseur (2019).

Properties near half-filling.—Near half filling, i.e., for $\mu\ll 1$, we can extract more quantitative information about the structure factor. Again, we classify strings as light and heavy at time $t$, depending on whether their spread up to time $t$ is primarily ballistic or diffusive. For light strings, diffusive corrections are a subleading effect (except at the front) so we treat light strings at the Euler level. For heavy strings near half filling, the diffusive broadening constant has the closed-form expression

[TABLE]

which coincides with the spin diffusion constant De Nardis et al. (2018b); Gopalakrishnan and Vasseur (2019); De Nardis et al. (2019). This expression is a sum over contributions from $s$-strings, and only strings with $s\eta\lesssim 1$ contribute. A slight distance $\mu$ away from half filling, the properties of strings with $s\mu\lesssim 1$ are similar to those at half-filling, while strings with $s\mu\gtrsim 1$ have exponentially suppressed density. When $\mu\ll\eta$ the contributions to Eq. (10) that would be modified by finite $\mu$ are already exponentially suppressed in $\eta$, so heavy strings with $\eta^{-1}\ll s$ diffuse with their $\mu=0$ diffusion constant (10) up to exponentially small corrections.

We now discuss the behavior of $C(x,t)$ near half-filling along rays with $x/t\neq 0$. For nonzero $\mu$, at late times, there are two different types of ballistic strings, depending on the size of $s\mu$. When $s\mu\gg 1$, the ballistic strings behave as in the low-filling limit: their density and velocity are both exponentially suppressed, and we recover Eq. (5). However, for lighter strings with $s\mu\lesssim 1$, the density is only suppressed algebraically as $\rho_{s}\simeq 1/s^{3}$ while the velocity is suppressed exponentially $v_{s}\simeq e^{-\eta s}$. Therefore, at a fixed position $x$, the largest string that has made it out to $x$ has index $s_{*}=\log(t/x)/\eta$. The density of such strings is $1/s_{*}^{3}$ while each carries a small dressed magnetization $s_{*}^{2}\mu$. Thus, in this regime, we have

[TABLE]

where the regime of validity of this result is controlled by $\mu$, and is to be understood on a logarithmic scale for $x/t$. Interestingly, the correlator at a fixed position $x$ grows logarithmically with time, as heavier strings carrying more magnetization appear at $x$. At longer times, the structure factor decays anomalously as (7). Exactly at half filling ($\mu=0$), the structure factor simplifies even further. All strings but the heaviest ones $s\to\infty$ become effectively neutral as $m^{\rm dr}_{s}\sim s^{2}\mu$ goes to 0 for $s\ll\mu^{-1}$, so the structure factor reads

[TABLE]

where we have used the sum rule (9) at half filling, and $v_{\infty}=0$. At half-filling, the structure factor is thus given by the heaviest strings De Nardis et al. (2018b), which are moving purely diffusively because of random collisions with lighter strings, with the spin diffusion constant (10).

Matrix product operator calculations.— We test these predictions by computing the structure factor $C(x,t)$ in the Heisenberg picture by time evolving $S^{z}_{i}$ using matrix product operator (MPO) techniques and the time-dependent density matrix renormalization group (tDMRG) White and Feiguin (2004); Vidal (2003); Schollwoeck (2011); Karrasch et al. (2012, 2013). We find that a fixed truncation error $\epsilon=10^{-8}$ is enough to obtain converged results, and we use a fourth order Trotter decomposition with time step $dt=0.1$. Our calculations are stopped when the bond dimension reaches $\chi\sim 2000$. This approach allows us to compute $C(x,t)$ for any filling (or temperature) from a single calculation of $S^{z}_{i}(t)$. Following Ref. Ljubotina et al. (2019) (see also Vasseur et al. (2015); Karrasch (2017) in the context of the Drude weight), we also compute $C(x,t)$ using a linear response quench setup when the system is initially prepared in a non-equilibrium state with chemical potential $\mu+\delta\mu/2$ in the left half of the system, and $\mu-\delta\mu/2$ in the right half. We then compute the density matrix $\rho(t)$ at time $t$ using MPO methods. Working at fixed truncation error, this approach allows us to reach similar time scales $t\sim 20$ to obtained fully converged results. If we work instead with fixed bond dimension MPOs (with bond dimensions $\chi=100,200,300,400$), the quench time-traces deviate from the exact MPO time trace at short times $t\lesssim 20$ for the return probability and appear to oscillate about it, though they give reasonably converged spatial profiles for $x\neq 0$ out to late times sup , as noted in Refs. Ljubotina et al. (2017, 2019).

These results are plotted in Fig. 2. For $\mu=0.5$ and $\eta=1.5$, the structure factor has a clear ballistic front due to magnons, with a broad diffuse feature behind it, as GHD predicts. The middle panel shows the local autocorrelator (i.e., return probability) as a function of $\mu$ at fixed $\eta=0.5$. Its behavior is consistent with a continuously varying power law that goes from approximately $1/2$ (in fact closer to $0.6$ due to the proximity to the isotropic point $\Delta=1$) at $\mu=0$ to nearly one at large $\mu$. The numerically extracted exponent $\gamma$ is much smaller than the asymptotic GHD prediction. To understand this discrepancy we have numerically evaluated the return probability using Eq. (8). As shown in the right panel of Fig. 2, the GHD curves curve downwards at short times, and only converge to their asymptotic slopes at extremely long times, as anticipated above.

Isotropic point.— Finally, we briefly discuss the structure factor at the isotropic point $\eta=0$ ($\Delta=1$). At half-filling, Eq. (10) implies that the diffusion constant diverges with the number of strings as $D\sim s$. Spin transport at half-filling is therefore superdiffusive Ljubotina et al. (2017); Ilievski et al. (2018), with a time-dependent diffusion constant that was argued to scale as $D(t)\sim t^{1/3}$ due to the anomalous behavior of heavy strings Gopalakrishnan and Vasseur (2019) (see also De Nardis et al. (2019)), consistent with numerical results Ljubotina et al. (2017, 2019).

By considering the approach to half filling at finite $\mu$ one can retrieve the dynamical exponent $x\sim t^{2/3}$ by a simple intuitive argument. A typical thermal state has Gaussian spatial fluctuations of its magnetization, so the effective local magnetization fluctuates as $1/\sqrt{L}$ over a distance $L$. On short enough length-scales, these fluctuations dominate over the average $\mu$. The system averages out these fluctuations and “realizes” it is at $\mu\neq 0$ on a length-scale such that $\mu\sim 1/\sqrt{L}$, i.e., the crossover length scales as $L(\mu)\sim 1/\mu^{2}$. Further, as $\mu\rightarrow 0$, magnetization is primarily transported by the heaviest available strings, for which $s_{*}\simeq 1/\mu$ and $v_{s_{*}}\simeq\mu$. The time it takes these strings to travel a distance $L(\mu)$ is given by $t(\mu)\sim 1/\mu^{3}$. The diffusion constant of such strings diverges as $D\sim\mu^{-1}$, which also gives the same scaling $t(\mu)\sim 1/\mu^{3}$. It follows that $L\sim t^{2/3}$. Moreover, the structure factor near half-filling can be written in the scaling form

[TABLE]

where the first term comes from strings with $s^{*}\sim 1/\mu$ and the second from lighter strings. At precisely half-filling the regular part vanishes as $\mu^{2}$, and only the anomalous part survives. The regimes of $C_{\mathrm{anom.}}(\zeta,\xi)$ are as follows. When both $\zeta,\xi\ll 1$, $C_{\mathrm{anom.}}\sim(t\mu^{3})^{-2/3}f(x/t^{2/3})$, where $f$ was numerically found to have the KPZ form Ljubotina et al. (2019). When $\zeta\gg\xi$ the anomalous part vanishes by causality. The late-time return probability $\xi\gg 1,\zeta\ll\xi$ is dominated by the tail of the heaviest common string, i.e., it goes as $1/(\mu t)$. Putting these together we have

[TABLE]

Meanwhile, the ballistic, regular part can be calculated following the logic of Eq. (11), so $C_{\mathrm{reg.}}(y)\sim 1/y^{2}$ for $y\ll 1$, implying that

[TABLE]

As $\mu\rightarrow 0$ spatial fluctuations of the magnetization dominate the dynamics. If we imagine dividing the system into a large number of hydrodynamic cells with magnetization $m(x,t)$, each cell will have a fluctuating diffusion constant $D[m]\sim 1/m$ and ballistic spin transport coefficient $j_{\rm ballistic}[m]\sim v[m]m\sim m^{2}$ set by its instantaneous magnetization (repeating the argument above with $m$ instead of $\mu$ as the cutoff). Combining these contributions into a hydrodynamic equation for $m$ yields a Burgers equation with a diffusion constant that is singular at low density (see also De Nardis et al. (2019)). We expect this to be compatible with KPZ scaling Kardar et al. (1986): Over a distance $\ell$, Gaussian fluctuations in the initial state lead to $m\sim 1/\sqrt{\ell}$, implying a diffusion constant $D\sim\sqrt{\ell}$ and ballistic velocity $v\sim 1/\sqrt{\ell}$, both implying $t(\ell)\sim\ell^{3/2}$. Moreover, the dominant nonlinearities in the Burgers equation involve anomalous high-density regions, for which the diffusion coefficient is well-behaved, so one might conjecture that the KPZ scaling function is also unaffected, as the numerical evidence Ljubotina et al. (2019) suggests. However, developing this nonlinear fluctuating hydrodynamics Mendl and Spohn (2013) for integrable systems is outside the scope of the present work.

Discussion.—In this work we used generalized hydrodynamics and its diffusive corrections to characterize the structure factor of the XXZ model in the easy-axis regime and at the isotropic Heisenberg point. We argued that even at nonzero magnetization, where ballistic transport is present, the local behavior of the autocorrelation function exhibits rich structure due to heavy “string” quasiparticles. In particular the autocorrelation function for $x\ll\sqrt{Dt}$, i.e., the “return probability,” vanishes with an anomalous exponent $\gamma=\min(\frac{1}{2}+\frac{|\mu|}{\eta},1)$ throughout this phase. Generic response functions therefore behave anomalously at fixed $q$ when $\omega\rightarrow 0$. This behavior is consistent with extensive simulations using MPO methods (Fig. 2). At the isotropic point we wrote down a scaling form for the structure factor, and provided an elementary derivation of the dynamical critical exponent Gopalakrishnan and Vasseur (2019); De Nardis et al. (2019). Many possible extensions present themselves, including a systematic derivation of fluctuating hydrodynamics and long-time tails near half-filling, and an understanding of the scaling properties as one approaches $\Delta\rightarrow 1$ from the easy-plane (“gapless”) regime.

Acknowledgments.—The authors thank Vincenzo Alba, Jacopo De Nardis, David Huse, Christoph Karrasch, and Vadim Oganesyan for helpful discussions. This work was supported by NSF Grant No. DMR-1653271 (S.G.), and US Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-SC0019168 (R.V.). The authors are grateful to the KITP, which is supported by the National Science Foundation under Grant No. NSF PHY-1748958, and the Program “The Dynamics of Quantum Information”, where part of this work was performed.