Generalized hydrodynamics, quasiparticle diffusion, and anomalous local relaxation in random integrable spin chains

TL;DR

This paper develops a generalized hydrodynamic framework for random integrable spin chains, revealing diffusive quasiparticle spreading and phase transitions between delocalized and quasi-localized states with distinct relaxation behaviors.

Contribution

It introduces a hydrodynamic theory incorporating disorder-induced diffusion in random integrable systems, highlighting novel phase transition phenomena and anomalous relaxation.

Findings

Diffusive broadening of operator fronts in random integrable systems.

Existence of a phase transition between delocalized and quasi-localized phases.

Anomalous power-law decay of autocorrelation functions in the localized phase.

Abstract

We study the nonequilibrium dynamics of random spin chains that remain integrable (i.e., solvable via Bethe ansatz): because of correlations in the disorder, these systems escape localization and feature ballistically spreading quasiparticles. We derive a generalized hydrodynamic theory for dynamics in such random integrable systems, including diffusive corrections due to disorder, and use it to study non-equilibrium energy and spin transport. We show that diffusive corrections to the ballistic propagation of quasiparticles can arise even in noninteracting settings, in sharp contrast with clean integrable systems. This implies that operator fronts broaden diffusively in random integrable systems. By tuning parameters in the disorder distribution, one can drive this model through an unusual phase transition, between a phase where all wavefunctions are delocalized and a phase in which low-energy wavefunctions are quasi-localized (in a sense we specify). Both phases have ballistic transport; however, in the quasi-localized phase, local autocorrelation functions decay with an anomalous power law, and the density of states diverges at low energy.

Figures

Equations

$$H=\sum _{j=1}^{L/2}$$

$$+K_{2j}\left([{\vec{\sigma }}_{2j}.\left({\vec{\sigma }}_{2j-1}\times {\vec{\sigma }}_{2j+1}\right){]}_{{\Delta }_{2j}^{-1}}+{\Delta }_{2j}^{-1}\right)$$

$$+J_{2j}^{(2)}\left({\vec{\sigma }}_{2j-1}.{\vec{\sigma }}_{2j+1}-1\right).$$

$$J_{2j}^{(1)}$$

$$K_{2j}$$

$$P(\xi )=\frac{1}{\sqrt{2\pi W^2}}{e}^{-{\xi }^2/2W^2}.$$

$$P(\xi )=\frac{\varphi }{2}{e}^{-\varphi |\xi |}.$$

$$H=$$

$$\displaystyle+\rm{tanh}(\xi_{2j})\left[\sigma^{y}_{2j-1}\sigma^{z}_{2j}\sigma^{x}_{2j+1}-\sigma^{x}_{2j-1}\sigma^{z}_{2j}\sigma^{y}_{2j+1}\right]\Bigg{]}.$$

$$H=$$

$$+\sum _{j=1}^{L/2}2i\tanh ({\xi }_{2j})\left({\mathrm{c}}_{2j-1}^{†}{\mathrm{c}}_{2j+1}-h.c.\right),$$

$${\partial }_t{q}_n(x,t)+{\partial }_x{j}_n(x,t)=0,$$

$${\partial }_t{\rho }_{j,\lambda }+{\partial }_x(v_{j,\lambda }^{\text{eff}}{\rho }_{j,\lambda })=0,$$

$$v_{j,\lambda }^{\text{eff}}=\frac{e_j^{\prime }(\lambda )}{p_j^{\prime }(\lambda )}=q_j\frac{e_j^{\prime }(\lambda )}{2\pi {\rho }_{j,\lambda }^T},$$

$$e_j(\lambda )=4J{A}_j(\lambda ),$$

$${\rho }_{j,\lambda }^T=q_j\frac{1}{N}\left(\sum _{k=1}^{N/2}{A}_j\left(\lambda +\frac{4}{\pi }{\xi }_{2k}\right)+\frac{N}{2}{A}_j(\lambda )\right),$$

$$A_j(\lambda )=\frac{\pi }{4}\frac{q_j}{\text{cosh}(\pi \lambda /2)}.$$

$$t_x=\Delta x\sum _{i=1}^M\frac{1}{v_i},$$

$$t_x=\Delta x\frac{2\pi }{e^{\prime }}\sum _{i=1}^M{\rho }_i^T.$$

$$\overline{t_x}=x\frac{2\pi \overline{{\rho }^T}}{e^{\prime }},$$

$$\sigma [t_x]=\sqrt{\frac{\Delta x}{2N_{\mathrm{sub}}}}\sqrt{x}\frac{2\pi \sigma \left[A(\lambda +\frac{2}{\eta }\xi )\right]}{e^{\prime }},$$

$$v^{\star }=\frac{e^{\prime }}{2\pi \overline{{\rho }^T}}={\overline{(v^{\text{eff}}{)}^{-1}}}^{-1}.$$

$$P_x(t)$$

$$\Gamma$$

$$P(x,t)\approx \frac{1}{\sqrt{2\pi \Gamma (v^{\star }{)}^{3}t}}{e}^{-(x-v^{\star }t{)}^2/2\Gamma (v^{\star }{)}^{3}t},$$

$$\rho (x,t)=\int \frac{1}{\sqrt{4\pi Dt}}{e}^{-(x-x_0-v^{\star }t{)}^2/4Dt}{\rho }^0(x^0)d{x}^0,$$

$$\boxed{{\partial }_t{\rho }_{j,\lambda }(x,t)+v_{j,\lambda }^{\star }{\partial }_x{\rho }_{j,\lambda }(x,t)=D_{j,\lambda }{\partial }_x^2{\rho }_{j,\lambda }(x,t),}$$

$$ϵ(x,t)$$

$$s_z(x,t)$$

$$D_q=\frac{1}{(1-q)\log L}\overline{\log {\sum }_x|{\psi }_x{|}^{2q}},$$

$${\rho }_{1,\lambda }^T={⟨\frac{\pi }{8\cosh (\pi \lambda /2+2\xi )}⟩}_{\mathrm{dis}}+\frac{\pi }{8\cosh (\pi \lambda /2)},$$

$$\frac{\varphi \pi }{4}{\int }_{-\infty }^{\infty }d\xi \frac{e^{-\varphi |\xi |}}{\cosh (\pi \lambda /2+2\xi )}.$$

$${⟨\frac{\pi }{8\cosh (\pi \lambda /2+2\xi )}⟩}_{\mathrm{dis}}\approx \frac{\varphi \pi e^{-\pi \lambda /2}}{2(\varphi -2)}(1-e^{-(\varphi -2)\pi \lambda /2})+\dots$$

$${⟨\frac{1}{\cosh (\pi \lambda /2+2\xi )}⟩}_{\mathrm{dis}}\sim e^{-\varphi \pi \lambda /4}.$$

Full text

I Introduction

The study of the dynamics of isolated, many-body quantum systems far from thermal equilibrium has attracted a lot of attention recently, fueled by recent experimental developments on ultracold atoms Bloch et al. (2008); Kinoshita et al. (2006); Langen et al. (2015), trapped ions Leibfried et al. (2003); Duan and Monroe (2010); Blatt and Roos (2012), nitrogen-vacancy centers Doherty et al. (2013); Schirhagl et al. (2014) and superconducting qubits Kelly et al. (2015) platforms. Addressing questions about non-equilibrium transport, thermalization and far-from-equilibrium dynamics pose notable challenges for theory as they are not susceptible to the general principles and methods that govern the physics of low-energy, equilibrium systems.

With the notable exception of many-body localized systems Nandkishore and Huse (2015); Vasseur and Moore (2016); Abanin et al. (2018), generic many-body systems are expected to be “chaotic”, and to thermalize under their own dynamics D’Alessio et al. (2016). This process can be understood as the scrambling of quantum information as it becomes non-local and inaccessible to physical, local measurements. After a local equilibration regime, thermalizing systems can be well-described by classical hydrodynamic equations associated with conserved quantities – typically, energy, particle number and momentum. These hydrodynamic equations describe the evolution of the system from local to global equilibrium. Another class of systems that escape thermalization in the traditional sense are quantum integrable systems, including experimentally relevant examples like the Heisenberg antiferromagnet and the Lieb-Liniger Bose gas in one dimension Kinoshita et al. (2006); Langen et al. (2015); Tang et al. (2018). Such systems have stable quasiparticle excitations even at high temperature, and they possess an extensive number of conserved quantities which strongly constrain their dynamics, and prevent them from thermalizing like generic chaotic systems  Calabrese and Cardy (2006); Prosen (2011); Caux and Essler (2013); Fagotti and Essler (2013); Essler and Fagotti (2016); Prosen and Ilievski (2013); Wouters et al. (2014); Ilievski et al. (2015a, b, 2016); Fagotti et al. (2014). However, contrary to many-body localized systems, integrable systems do thermalize in a generalized sense, as they eventually reach a maximum entropy steady state described by a generalized Gibbs ensemble (GGE) Rigol et al. (2007); Caux and Konik (2012); Rigol (2009); Langen et al. (2015). Such steady-states can exhibit non-zero currents, and are commonly referred to as non-equilibrium steady states (NESS) in the literature Bernard and Doyon (2016); Vasseur and Moore (2016), even though they are natural equilibrium states for integrable systems.

A major step in understanding the non-equilibrium dynamics of quantum integrable systems was the formulation of what is now known as “generalized hydrodynamics” (GHD) Castro-Alvaredo et al. (2016); Bertini et al. (2016), which are Euler hydrodynamics equations (0th order hydrodynamics) obtained in the large space-time limit where the system is locally in equilibrium. While the prospect of solving infinitely-many hydrodynamic equations (one for each conserved quantity in the system) originally appeared daunting, GHD can be conveniently formulated in the basis of quasiparticle excitations: in that language, they can be naturally interpreted as describing a semi-classical gas of solitons (quasiparticles) El (2003); El and Kamchatnov (2005); Bulchandani et al. (2018); Doyon et al. (2018). The key ingredient of GHD is the effective group velocity $v^{\rm eff}$ of the quasiparticles Bonnes et al. (2014); Castro-Alvaredo et al. (2016); Bertini et al. (2016) which depends on the density of all the other quasiparticles in the presence of interactions: at the semi-classical level, quasiparticle wavepackets pick up a phase shift when they collide, leading to a Wigner time-delay. This approach was successfully applied to two-reservoir setups Castro-Alvaredo et al. (2016); Bertini et al. (2016) and more generally to locally equilibrated inhomogeneous initial states Bulchandani et al. (2017); Doyon et al. (2017); Piroli et al. (2017), and has helped addressing a number of key questions in the field concerning Drude weights Ilievski and De Nardis (2017a); Bulchandani et al. (2018); Doyon and Spohn (2017); Ilievski and De Nardis (2017b); Collura et al. (2018), external potentials and traps Doyon and Yoshimura (2017); Cao et al. (2018); Caux et al. (2017), correlation functions Doyon (2018), entanglement dynamics Alba and Calabrese (2017), or even large-deviation functions Myers et al. (2018).

The GHD framework was recently generalized to include diffusive effects in interacting integrable models De Nardis et al. (2018a); Gopalakrishnan et al. (2018); De Nardis et al. (2018b); Gopalakrishnan and Vasseur (2018) – corresponding to “1st order” or Navier-Stokes hydrodynamics, with important consequences for the nature of spin transport in XXZ spin chains De Nardis et al. (2018b); Gopalakrishnan and Vasseur (2018). In particular, Ref. De Nardis et al., 2018a provided a general exact expression of the “diffusion matrix” of the quasiparticles using a form factor expansion of the Kubo formula. Intuitively, diffusive corrections can be seen to arise as follows Gopalakrishnan et al. (2018): the effective velocity of a given quasiparticle depends on the density of quasiparticles in a mean-field fashion that ignores fluctuations; reintroducing thermal fluctuations naturally leads to a diffusive broadening of quasiparticle trajectories in a generic GGE state. Two key ingredients are needed for such diffusive corrections to be present: density-dependent velocities and thermal fluctuations (non-zero entropy states). This immediately implies that in non-interacting integrable models where the group velocity is obtained from band theory and is independent of density, there should be no diffusion Spohn (2018). For non-interacting systems, it is thus natural to expect that the lowest order correction to ballistic GHD comes from higher order derivative terms, which lead to $t^{1/3}$ spreading governed by the Airy kernel (see Ref. Fagotti, 2017 and references therein). Clearly, such higher-order corrections are subleading in the presence of diffusive $t^{1/2}$ spreading.

In this paper, we study a class of integrable random spin chains which support diffusive corrections even in the absence of interactions. These spin chains are a special limit of a more general class of random interacting spin chains that remain integrable. In one-dimensional free-particle problems, disorder generically leads to Anderson localization. Though Anderson localized systems are “integrable” in a sense, here we will use the term “integrability” exclusively to refer to Bethe-ansatz solvable systems with stable ballistically propagating quasiparticles. There are examples of integrable models with impurities Vega and Woynarovich (1992); de Vega et al. (1994); Essler et al. (2018) where disorder is correlated in such a way that integrability in this sense is preserved. Such systems were recently shown to exhibit ballistic transport even at strong disorder in Ref. Essler et al., 2018. It is natural to expect such random systems to exhibit diffusive corrections to ballistic transport even without interactions, as quasiparticles scatter off random static impurities and thus undergo biased random walks. These models illustrate that whereas non-interacting systems are often said to be always integrable, integrability for a random system leads to correlated disorder that can allow one-dimensional random systems to escape Anderson localization. Although most states in these models are only weakly affected by the disorder, quasiparticle states near energy $|E|=0$ have properties that are sensitive to the tails of the disorder distribution. We find disorder distributions for which these quasiparticles have vanishing velocities, so that the behavior of local autocorrelation functions is anomalous. We find that the onset of anomalous behavior in the disorder-averaged autocorrelation functions is associated with the onset of a divergence in the density of states at zero energy, as well as a form of quasi-localization of the low-energy wavefunctions that we discuss below. We compute this anomalous relaxation exponent using GHD. There are strong local correlations between the density of states and the quasiparticle velocity, as both are dominated by rare regions; our work shows how GHD can be adapted to incorporate these rare region effects.

In this work, we propose a hydrodynamic theory to describe such random integrable spin chains, including diffusive corrections due to disorder. The plan of this paper is as follows: in Section II, we recall the definition of a family of random integrable spin chains recently studied in Ref. Essler et al., 2018, and we briefly review their thermodynamics. We then formulate a coarse-grained GHD theory for the dynamics of such systems in Section III, with an emphasis on diffusive corrections due to disorder in non-interacting settings (for which disorder is the only possible source of diffusive corrections). This framework is applied to study non-equilibrium spin and energy transport, and the predictions are compared to numerical results in Section IV. Transport is dominated by fast quasiparticles with energies well away from $|E|=0$. We turn in Sec. V to a more careful discussion of states near $|E|=0$. For these states we find a quasi-localization transition; in the quasi-localized regime, low-energy wavefunctions consist of a few local peaks, quasiparticle velocities vanish, and the density of states diverges in the $|E|\rightarrow 0$ limit. Consequences for operator spreading and scrambling are briefly discussed in Section VI, and a discussion and outlooks for future works are gathered in Section VII.

II Random integrable spin chains

In this section, we introduce a family of random integrable spin chains, following closely Ref. Essler et al., 2018, and we briefly review their thermodynamic Bethe ansatz solution.

II.1 Hamiltonian

Let us consider a random XXZ spin-$\frac{1}{2}$ chain $H=\sum_{i}J_{i}\left[\vec{\sigma}_{i}.\vec{\sigma}_{i+1}\right]_{\Delta_{i}}$, where $\left[\vec{\sigma}_{j}.\vec{\sigma}_{k}\right]_{\Delta_{i}}$ is a shorthand notation for $\sigma^{x}_{j}\sigma^{x}_{k}+\sigma^{y}_{j}\sigma^{y}_{k}+\Delta_{i}(\sigma^{z}_{j}\sigma^{z}_{k}-1)$. In the clean (homogeneous) case, this model is integrable, but the introduction of disorder immediately breaks integrability, and leads to a model that is either chaotic or many-body localized Vosk and Altman (2013); Vasseur et al. (2016). However, it is possible to preserve integrability Vega and Woynarovich (1992); de Vega et al. (1994); Essler et al. (2018) by introducing next-to-nearest neighbor interactions, and by carefully choosing the inhomogeneous couplings:

[TABLE]

The first line of the hamiltonian corresponds to an XXZ interaction, while the last line is an isotropic Heinsenberg interaction. The middle line is more unusual, as it involves three spins. The parameters in the hamiltonian are given by

[TABLE]

with $\xi_{2j}$ a random coupling, while $\eta$ is an overall global parameter that parametrizes the interaction strength. For $\xi_{2j}=0$, one recovers the usual XXZ spin chain. A remarkable feature of this model is that it remains integrable for any choice of the inhomogeneous couplings $\xi_{2j}$. Away from the zero-energy limit, the properties of this model are insensitive to details of the disorder distribution. Therefore, except as specified below (i.e., in Sec. V and subsequently), we will take the $\xi_{2j}$ couplings to be random variables drawn from the gaussian distribution

[TABLE]

Later, we will also consider the exponential distribution

[TABLE]

for which a sharp quasi-localization transition exists.

We will be especially interested in the special point $\eta=\pi/2$. For this value of $\eta$, the XXX part of the Hamiltonian is set to zero, leaving behind a random XX model with three spin interactions,

[TABLE]

The above hamiltonian can be diagonalized via Jordan-Wigner transformation reducing it to a free fermion model

[TABLE]

where the $\xi_{j}$’s are random parameters used in (II.1), extended to odd sites via relation, $\xi_{2j-1}=\xi_{2j}$. We use periodic boundary conditions for the fermions for an even number of sites. While non-interacting and disordered, this model was shown to escape Anderson localization in Ref. Essler et al., 2018, and to exhibit ballistic transport of conserved quantities. It is then natural to ask if transport properties of this model can be captured using generalized hydrodynamic equations, properly adapted to deal with the quenched disorder. If so, it is natural to expect disorder to lead to new hydrodynamic effects, such as diffusion. In the following, we will mostly focus on the special point $\eta=\pi/2$ (eqs. (II.1) (5)), though we expect our approach to generalize to any value of $\eta$. This will be convenient as the free fermion representation of this model allows one to simulate numerically the non-equilibrium dynamics of this system easily, and more importantly, the absence of interactions will allow us to isolate the effect of disorder on diffusion.

II.2 Thermodynamics

The model introduced above admits an exact solution by Bethe Ansatz. In the following, we very briefly review the Thermodynamic Bethe Ansatz (TBA) approach to integrable systems, with the main goal of introducing some notation and language that will be used in the remainder of this paper 111For a thorough review of TBA, we refer the interested reader to Takahashi (1999)..

The solutions of the Bethe equations are expressed in terms of quasi-momenta or rapidities, and as one takes the thermodynamic limit, one introduces a density of allowed quasi-momenta $\rho^{T}(\lambda)$, corresponding to a total density of states. In general, we will also introduce an additional discrete quasiparticle label $j$; for XXZ-like spin chains, $j$ is known as the string index. Each allowed quasi-momentum state can be occupied or not in a given macrostate, so one defines the hole density, $\rho^{h}_{j}$, and the density of occupied states (quasiparticle density) $\rho_{j}$, and we have $\rho^{T}(\lambda)=\rho_{j}(\lambda)+\rho^{h}_{j}(\lambda)$ by definition. The particle density together with hole density fully characterize the state of the system. The expectation value of a conserved quantity in such a state can be written as $Q=\int\rho_{j}(\lambda)q_{j}(\lambda)d\lambda$ where $q(\lambda)$ is the single-particle charge eigenvalue of the string of type $j$ with quasi-momentum $\lambda$.

Thermodynamic equilibrium properties can be computed by writing the entropy and energy in terms of $\rho$ and $\rho^{h}$, and then maximizing the free energy as usual for a given set of Lagrange multipliers corresponding to a GGE. This leads to the so-called Yang-Yang equation Yang and Yang (1969). The Yang-Yang equation together with the Bethe equation are enough to fully determine the quasiparticle densities $\rho_{j}$ and $\rho^{h}_{j}$, corresponding to a given GGE state. In the following, we will promote these variables locally by assuming local (generalized) equilibrium.

III Generalized hydrodynamics approach

We now formulate a hydrodynamic description of such random integrable systems. The evolution of chaotic quantum systems from local to global equilibrium is described by the framework of hydrodynamics. In that regime, one imagines chopping off the system into hydrodynamic cells that are big enough to assume equilibrium within each cell, but very small compared to the total system size. This separation of scales allows one to assume local equilibrium, where Lagrange multipliers like temperature or chemical potential are allowed to depend on position and time.

There is one hydrodynamic equation per conserved quantity in the system – any other information about the system is “scrambled” by the quantum dynamics into non-local entanglement that is not measurable by local observables. For each conserved quantity $Q_{n}=\sum_{x}q_{n}(x)$, we can write a continuity equation

[TABLE]

where we restricted ourselves to one spatial dimension. Assuming local equilibrium then leads to a relation $j_{n}=F_{n}[\{q_{m}\}]$ between the currents $j_{n}$ and the conserved charges $q_{m}$ at a given position $x$: this relation is an equilibrium property and gives rise to Euler hydrodynamic equations that govern ballistic transport properties. More generally, one can perform a gradient expansion of the (expectation value of the) currents in terms of the charges, where contributions to the currents coming from gradient terms $D_{nm}\partial_{x}q_{m}$ correspond to diffusive contributions to hydrodynamics (see e.g. Spohn, 2012). The diffusion constants $D_{nm}$ are not entirely given by equilibrium properties, and have to be determined by other means such as the Kubo formula, or by using kinetic theory calculations. Once the transport coefficients characterizing the relation between currents and charges are known, hydrodynamics provides a simple set of classical, partial differential equations that govern the non-equilibrium dynamics of the system.

III.1 Generalized hydrodynamics

The hydrodynamic framework summarized above is completely general, and it was successfully applied to integrable systems Castro-Alvaredo et al. (2016); Bertini et al. (2016): the resulting framework is now known as generalized hydrodynamics (GHD), as it describes systems in local GGE equilibrium. There, local equilibrium is characterized by the densities $\rho_{j,\lambda}(x,t),\rho^{h}_{j,\lambda}(x,t)$, with the charges $q_{n}(x)=\sum_{j}\int\rho_{j,\lambda}(x)q_{n,j,\lambda}(x)d\lambda$ and currents $j_{i}(x,t)=\int v^{\text{eff}}_{i,\lambda}(x,t)\rho_{i,\lambda}(x,t)d\lambda$ (ignoring gradient (diffusive) corrections). $v_{j}^{\text{eff}}$ is interpreted as a group velocity which is a functional of the quasiparticle density in general. The continuity equations for the conserved charges then imply a continuity equation for the quasiparticle density Castro-Alvaredo et al. (2016); Bertini et al. (2016)

[TABLE]

For a non-interacting system, $v^{\text{eff}}_{j}$ is independent of $\rho_{j}$. In that case, in clean systems, diffusive corrections to (8) are believed to be absent.

We now turn to our specific example (5), which is noninteracting in the fermionic language. Recall that in the noninteracting limit, $\eta=\pi/2$ in Eq.(2). In this limit, there are two strings $j=1,2$, and their group velocity is simply given by Bonnes et al. (2014); Castro-Alvaredo et al. (2016); Bertini et al. (2016)

[TABLE]

where $e_{j}$ is the quasiparticle energy given by

[TABLE]

$p_{j}$ the quasiparticle momentum, and the total density of states $\rho^{T}_{j,\lambda}=\rho_{j,\lambda}+\rho^{h}_{j,\lambda}$ is given Takahashi (1999); Essler et al. (2018) by the Bethe equation (which is particularly simple since the model is non-interacting)

[TABLE]

for a system of size $N$, with the function $A_{j}(\lambda)$ defined as

[TABLE]

In these equations, $q_{1}=1,q_{2}=-1$. In the interacting case, all these quantities would be “dressed” and would become functional of the quasiparticle densities $\rho_{j}$.

As expected since the model is non-interacting, the group velocity does not depend on the density $\rho$ of the other quasiparticles. However, it does depend on the inhomogeneous variables $\xi_{i}$. It is thus clear that some kind of averaging over these random variables needs to be done in order to formulate a hydrodynamic theory of this random quantum spin chain.

III.2 Averaging and coarse-graining

In order to average of the random variables $\xi_{i}$, we go back to the physical picture of hydrodynamics and divide the system into mesoscopic hydrodynamic cells large enough to be in the thermodynamic limit. Let the system length be $L$ and it be divided in $N\gg 1$ sub-cells of size $\Delta x=L/N\gg a$ with $a$ the lattice spacing. For a given disorder realization, the velocity in each hydrodynamic cell is given by (9). These velocities depend on our choice of sub-cell division, but as we will see this dependence drops out of the final result.

Given these velocities we can easily construct the trajectory of a given quasi-particle from its initial position $x=0$. Let us find the time required for a quasiparticle of type $j$ with rapidity $\lambda$, initially at $x=0$ to reach $x=M\Delta x$. It is given by 222We will momentarily omit the string number and rapidity/momentum subscripts to make the equations more readable.,

[TABLE]

with $v_{i}$ the velocity the $i^{th}$ cell. We then have using (9)

[TABLE]

Since $\rho^{T}$ is a sum of random variables, we can use the central limit theorem to deduce that both $\rho^{T}$ and $t_{x}$ are Gaussian distributed (provided the hydrodynamic cells are large enough, and $M\gg 1$). This also shows that the result is largely independent of the distribution chosen for the random parameter $\xi$, as long as the central limit theorem is applicable, as is the case for the distributions considered in this paper.

Thus we have following result: the time taken for a quasiparticle to move over a distance $x$ is Gaussian distributed, with the average time being given by

[TABLE]

where $\overline{\rho^{T}}$ is the disorder average of $\rho^{T}$. This quantity is clearly independent from our choice of hydrodynamic cells – it only depends on $x$, not $M$, or $\Delta x$. The standard deviation reads

[TABLE]

where $N_{\rm sub}$ is the number of lattice sites inside a cell, and $\sigma[A(\lambda+\frac{2}{\eta}\xi)]$ is the standard deviation of the function defined in equation (12). Note that $\sqrt{\frac{\Delta x}{2N_{sub}}}=\sqrt{\frac{a}{2}}$ with $a$ the lattice spacing, implying that the standard deviation of $t_{x}$ is also independent of our choice of hydrodynamic cells. Thus we conclude that the distribution of $t_{x}$ does not depend on the partition of hydrodynamic cells, and is well defined.

We define the average velocity $v^{\star}$ via the relation $v^{\star}\overline{t_{x}}=x$. This yields

[TABLE]

Note that this is not the average of $v^{\text{eff}}$ over disorder.

The probability distribution of the time it took for quasiparticle to move over a distance $x$ thus reads

[TABLE]

This process is called temporal diffusion Boon et al. (2003); Langmann and Moosavi (2018), as it looks like an usual diffusion process where the roles of space and time are exchanged. However, in the hydrodynamic limit, the spreading of the distribution is confined to region $(x/v^{\star}-t)^{2}=\mathcal{O}(\Gamma x)$ or $x=v^{\star}t\left(1+\mathcal{O}(\sqrt{\frac{\Gamma v^{\star}}{t}})\right)$. Thus in the limit $\frac{\Gamma v^{\star}}{t}\ll 1$ we can replace $x$ by $v^{\star}t$ and get

[TABLE]

which corresponds to a biased random walk. A similar temporal diffusion equation recently appeared in the context of energy transport in a random conformal field theory Langmann and Moosavi (2018). In all the numerical results below, we have checked that the difference between the temporal and ordinary diffusion descriptions are negligible in the hydrodynamic limit.

For generic initial condition of the quasiparticles, $\rho^{0}(x,t=0)$, the evolution should thus reads

[TABLE]

with $D=\Gamma(v^{\star})^{3}/2$ since the quasiparticles are non-interacting. Reintroducing the string and rapidity labels, we find that the quasiparticle density satisfies the following hydrodynamic equation

[TABLE]

where $D_{j,\lambda}\equiv\frac{\Gamma_{j,\lambda}(v^{\star}_{j,\lambda})^{3}}{2}$ is a diffusion constant due to the disorder. We emphasize that the transport coefficients $v^{\star}$ and $D$ in this equation do not depend on the details of our coarse graining procedure – in particular, they do not depend on the size of the hydrodynamic cells $\Delta x$ as long as $L\gg\Delta x\gg a$. We emphasize that contrary to diffusive corrections due to thermal fluctuations, the diffusive term in the above equation is diagonal in the quasiparticle basis, and occurs even in the non-interacting case. For interacting random chains, we expect the argument above to carry over for the average velocity (16), but the diffusion matrix should be more complicated as it will have contributions from interactions and disorder. We leave a detailed analysis of the interacting case for future work.

IV Non-equilibrium transport

We now use the hydrodynamic equation derived above to study non-equilibrium energy and spin transport in the random spin chain (5). This will also allow us to benchmark and test the validity of the hydrodynamic approach, and investigate the importance of the diffusive terms due to disorder. The energy and spin densities can be expressed in terms of the quasiparticle densities as

[TABLE]

where $n_{j}$ is given in our case by: $n_{1}=n_{2}=1$.

IV.1 Energy transport

We first discuss energy transport. We consider a lattice of $L=200$ sites, prepared it in a thermal state with $\beta=1$, but with an interval of $32$ sites in the middle of the system prepared in an infinite temperature state ($\beta=0$). The time evolution of the energy density profile can then be straightforwardly extracted from the free fermion representation (II.1). We note that due to next to nearest interactions we have to define energy density on two sites instead of one. We also emphasize that since the initial condition has a mirror symmetry about the middle of the lattice, one expects the evolution to be symmetric as well. However, this is only guaranteed if the sample of random variables $\xi_{i}$ is symmetric about $\xi=0$ — this is due to the presence of imaginary interaction in Hamiltonian (II.1) which breaks this mirror symmetry. To retain the mirror symmetry exactly in the numerical data, we have ensured that the sample taken for the random variables are symmetric.

We then compare the numerical results from the solution of the hydrodynamic equation (21), combined with formula (22). Fig. 1 show the hydrodynamic prediction for the energy density (solid line) compared to lattice data, for different disorder strengths $W=0$ (clean case), $W=0.6$ and $W=2.0$. We also show the ballistic hydrodynamic prediction ignoring diffusive corrections (dashed line) – formally setting $D=0$ in the hydrodynamic equation. The agreement between the numerics and hydrodynamics is excellent, and we find that diffusive corrections are needed to accurately describe the numerical data, especially at stronger disorder.

IV.2 Spin transport

We also consider spin transport starting from an initial domain wall state $\left|\psi_{0}\right>=\left|\uparrow\uparrow\dots\uparrow\downarrow\dots\downarrow\downarrow\right>$ (see Ref. Collura et al., 2018 and references therein). This initial state has been considered for clean XXZ spin chains in several recent works: it leads to a steady-state with zero entropy, where diffusive terms due to interactions and thermal fluctuations are expected to vanish. The hydrodynamic prediction for the local magnetization and numerical data are compared in Fig. 2. Here also, eq. (21) provides an excellent description of the numerics and our data clearly shows the presence of diffusion in this non-interacting setting where no thermal fluctuations are present Bulchandani and Karrasch (2018). We expect this initial state to be also useful to extend our approach to interacting random spin chains, as it would allow one to isolate diffusive terms due to from disorder (since the diffusion matrix coming from interactions vanishes). Comparing the predictions from hydrodynamics to numerical results in the interacting case would however be very numerically demanding because of the average over disorder, and is left for future work.

V Properties of low-energy quasiparticles

So far, our discussion has focused on macroscopic energy and spin transport, which is dominated by fast quasiparticles. In addition to these typical, fast, quasiparticles, however, these models also have slow quasiparticles at energy $|E|\approx 0$. These are important, e.g., for low-temperature transport, as well as for the behavior of local autocorrelation functions at late times. We discuss the nature of these quasiparticles here. We first explore the properties of wavefunctions, both numerically and analytically, and find that these undergo a quasi-localization transition for the exponential disorder distribution (4). We then apply the thermodynamic Bethe ansatz results in Sec. III to study the asymptotic behavior of the velocity and density of states as $|E|\rightarrow 0$. Our discussion here is confined to the non-interacting model (5).

V.1 Spatial correlations of eigenstates

We first discuss the properties of single-particle eigenstates. Because of the Bethe-ansatz integrability of the model, eigenstates are never localized in the conventional sense, and most are completely delocalized Essler et al. (2018). However, the states closest in energy to zero have anomalous properties, which show up in the hydrodynamic framework as very slow velocities. We now address these properties more directly. Fig. 3 shows that the inverse participation ratio (IPR) $\sum_{x}\left|\psi_{x}\right|^{4}$ rises steeply for states very near zero energy; thus, these are less spread out than the rest of the spectrum, which is fully delocalized. In the rest of this section we will focus on the few relatively localized states near zero energy.

Given the anomalously large IPRs of states near zero energy, it is natural to investigate their spatial structure. It turns out this structure is simple and striking. The upper panel of Fig. 4 plots the spatial profile of each eigenstate against the eigenstate index. States near $|E|=0$ are halfway up the $y$ axis; they are evidently strongly peaked at certain lattice sites, and the peak locations are the same for each such state. In other words, these states exhibit a very strong form of Chalker scaling Chalker and Daniell (1988); Chalker (1990). Further, the peak locations and intensities follow a simple pattern (lower panel of Fig. 4): the amplitude of each peak is $\propto\cosh(\xi_{j})$, where $\xi_{j}$ is the nearest neighbor hopping at the peak. These observations can be qualitatively understood in terms of a simple picture in which the state has a roughly uniform current, and to maintain this uniform current the density piles up near weak links. Thus some form of “hydrodynamics” applies even at the level of individual wavefunctions in individual realizations. As one moves away from zero energy, the states lose intensity first at the strongest maxima (Fig. 4); higher-energy wavefunctions bypass these weak links instead of piling up near them.

We find, numerically, that these anomalous states are supported entirely on even sites (i.e., sites for which the next-nearest neighbor coupling is absent). This is intuitively plausible since no value of $\xi_{j}$ can simultaneously suppress both the nearest-neighbor and next-nearest-neighbor hopping out of an odd site.

V.2 Multifractality and quasi-localization

The eigenfunctions near zero energy are sharply peaked at a few (even) sites, suggesting that they might be critical rather than conventionally delocalized. The nature of these states is highly sensitive to the tails of the disorder distribution. So far we have considered Gaussian distributions of the $\xi_{j}$; numerical extraction of the IPR at the accessible system sizes does not settle whether the wavefunctions are localized or critical. However, the relationship between the wavefunction amplitudes and $\cosh(\xi_{j})$ allows us to address this question semi-analytically. Among $L$ Gaussian random variables, the largest value is likely to be $\exp(-\xi^{2}/W^{2})\simeq 1/L$, so $\xi=W\sqrt{\log L}$. The (unnormalized) wavefunction amplitude at this site is therefore $\cosh(W\sqrt{\log L})\sim\exp(W\sqrt{\log L})$. Meanwhile, the typical value of $\cosh(\xi_{j})$ is of order unity, and occurs $O(L)$ times. Thus the weakest link in a typical sample does not affect its properties overall. We therefore expect that the low-energy wavefunctions in the Gaussian case are asymptotically not multifractal.

This reasoning also suggests that to get multifractal wavefunctions, it suffices to change the disorder distribution from Gaussian to exponential, $P(\xi)\propto\exp(-\phi\left|\xi\right|)$. In this case, repeating the argument in the previous paragraph gives that the largest typical peak has amplitude $L^{1/\phi}$. When $\phi<2$ the wavefunction piles up at the weakest link, and is “quasi-localized” in the sense that its IPR is independent of system size. Numerical simulations clearly show this quasi-localized behavior at small $\phi$, as well as the expected multifractal behavior at larger $\phi$ (Fig. 5(d)). We characterize multifractality through the quantity Evers and Mirlin (2008)

[TABLE]

where $\overline{O}$ refers to the disorder average of the quantity $O$. For each sample, we take the lowest-$|E|$ state.

In the quasi-localized regime, $D_{q}=0$ for $q>1$. Strictly speaking, this statement is asymptotically true, in the following sense: wavefunctions at any fixed energy away from zero are delocalized, but as one approaches zero energy their IPR approaches a fixed, size-independent value (Fig. 3). Thus, there is no mobility edge. Although the IPR is finite, unlike an Anderson insulator, the quasi-localized phase has $D_{q}>0$ for $q<1$: although the wavefunction has a finite fraction of its weight concentrated in a few peaks, the rest of the weight is spread out evenly rather than falling off exponentially away from the peaks (Fig. 6). To our knowledge this quasi-localization property has not been previously noticed. Finally, we remark that although its IPR is independent of system size, this does not imply that a particular state will be unaffected by adding sites to one end of the system: the wavefunction will remain very sharply peaked at the weakest link, but the location of this link will move as sites are added to the system.

V.3 GHD approach to quasi-localization

We now use the results in Sec. III to consider the behavior of quasiparticles in the $|E|\rightarrow 0$ limit. From Eq. (10) we see that these quasiparticles correspond to large $|\lambda|$: in fact, $e_{j}(\lambda)\sim e^{-\pi|\lambda|/2}$. In what follows we suppress the index $j$ and denote the energy as $E$. The limiting behavior of the velocity and the density of states is sensitive to the tails of the disorder distribution. For the sample-averaged density of states, Eq. (11) yields

[TABLE]

where $\left\langle\dots\right\rangle_{\rm dis}$ denotes the average over disorder. If the disorder is bounded or falls off faster than exponentially, one can safely approximate $\cosh(x)\approx e^{x}/2$ for large enough $\lambda$. The quasiparticle velocity (given by Eq. (9)) therefore remains nonzero in the $\lambda\rightarrow\infty$ limit, so transport is asymptotically ballistic (but with a slower velocity). Likewise, the density of states remains finite. To find the density of states, one notes that the number of states in a rapidity interval $\delta\lambda$ is given by $\rho^{T}_{1,\lambda}\delta\lambda$. These cover an energy window $\delta E=E^{\prime}(\lambda)\delta\lambda=e^{-\pi\lambda/2}\delta\lambda$. Thus, $\rho(E)=\rho^{T}_{1,\lambda}/E^{\prime}(\lambda)$, which remains finite for Gaussian distributed disorder in the $|\lambda|\rightarrow\infty$ limit (although in practice the suppression is quantitatively quite large).

However, for distributions $P(\xi)$ with exponential or slower tails, computing the expectation value (25) is more subtle. We focus on the exponential case $P(\xi)=\frac{1}{2\phi}\exp(-\phi|\xi|)$, which was discussed numerically above. In this case, the expectation value (25) reads

[TABLE]

To get compact expressions for the asymptotics, we approximate $\cosh(\pi\lambda/2+2\xi)\approx\frac{1}{2}e^{|\pi\lambda/2+2\xi|}$. There are two cases. When $\phi>2$, the integral is dominated by small $|\xi|$ and we get

[TABLE]

where $\ldots$ indicates terms that do not become singular in the limit $\phi\rightarrow 2$. Thus the zero-energy density of states and velocity remain finite when $\phi>2$, but respectively diverge and vanish as $\phi\rightarrow 2^{+}$. Note that there are nonanalytic corrections to the density of states, even in this regime: specifically, $|\rho(E)-\rho(0)|\sim E^{\phi-2}$.

In the opposite limit $\phi<2$, Eq. (25) is dominated by $|\xi|\approx\pi\lambda/4$. In this case we have instead

[TABLE]

The velocity then vanishes as $v(\lambda)\sim\exp[-(\pi\lambda/2)(1-\phi/2)]$, or equivalently $v(E)\sim|E|^{1-\phi/2}$. Correspondingly the density of states $\rho(E)$ diverges as

[TABLE]

This behavior of the density of states is borne out numerically (Fig. 7).

Thus, the dispersion relation of elementary excitations within generalized hydrodynamics agrees with our simple counting estimate in Sec. V: a transition occurs when $\phi=2$. When $\phi>2$ the velocity approaches a finite value as $|E|\rightarrow 0$, so quantities such as the local autocorrelation function behave in an asymptotically ballistic fashion. On the other hand when $\phi<2$ the velocity vanishes as $|E|\rightarrow 0$, and the local autocorrelation function will in general be anomalous. It is interesting to note that despite the very local character of the rare low-energy states due to almost disconnected sites, they appear naturally as slow quasiparticles in the hydrodynamics framework.

VI Operator Spreading and Local Autocorrelations

VI.1 Front broadening

We close this paper by discussing the consequences of our results for operator spreading in random integrable systems. The dynamics of operator spreading has attracted a lot of attention in recent years, due to its possible connection to many-body quantum chaos Larkin and Ovchinnikov (1969); Shenker and Stanford (2014); Maldacena et al. (2016); Liu and Suh (2014); Casini et al. (2016); Mezei and Stanford (2017); Brown and Fawzi (2012); Nahum et al. (2017, 2018); von Keyserlingk et al. (2018); Khemani et al. (2018); Rakovszky et al. (2018); Chan et al. (2018); Bertini et al. (2018); Lin and Motrunich (2018); Xu and Swingle (2018). Under unitary dynamics, initially local operators spread in space ballistically (unless the system is many-body localized), with an operator “front” or light-cone that generically broadens diffusively as  $t^{1/2}$ in one dimensional chaotic (non-integrable) quantum systems. This is in sharp contrast with (clean) non-interacting systems have an operator front that broadens as $t^{1/3}$  Fagotti (2017), governed by an Airy kernel. However, interacting integrable systems have been argued to have a front that also broadens as $t^{1/2}$  Gopalakrishnan (2018); Gopalakrishnan et al. (2018), just like chaotic systems. (With the notable exception of non-generic zero entropy initial states like the spin domain wall initial state discussed above Bulchandani and Karrasch (2018).) In integrable systems, the operator front is governed by the fastest quasiparticle in the system Gopalakrishnan et al. (2018): from our hydrodynamic description of random integrable spin chains, it is natural to expect this front to broaden diffusively even in the non-interacting case, since the fastest quasiparticle follows a biased random walk. For these systems, the diffusive broadening of the operator front is due to the local disorder which causes the fastest quasiparticle (on average) to “wiggle” around its average trajectory.

In order to characterize operator spreading in our system, we compute a specific out of time ordered commutator (OTOC) Larkin and Ovchinnikov (1969); Shenker and Stanford (2014); Maldacena et al. (2016) given by $C(x,t)=\frac{1}{2}\left\langle[n_{x}(t),n_{0}(0)]^{2}\right\rangle$ in a given thermal state, with $n_{x}=c_{x}^{\dagger}c_{x}$ the particle density in the fermionic language. Fig. 8 shows $C(x,t)$ averaged over 100 disorder realizations for $W=1$ and temperature $T=1$, where a clear ballistic light cone can be observed. We also remark that the OTOC shows significant weight that remains near $x=0$ even after a long time. This behavior has to do with the anomalous low-energy states (Sec. V), and we will return to it below.

We now focus on the operator front (light-cone), and investigate whether it broadens as $t^{1/3}$ as expected for clean non-interacting systems, or as $t^{1/2}$ due to diffusive effects. We start with the clean case (XX model), and show scaling collapses of the OTOC near the front. It is clear from Fig. 9 that the data collapses almost perfectly using a $t^{1/3}$ ansatz. This is consistent with expectations in the clean case, the OTOC shows oscillations characteristic of an Airy Kernel associated with these higher-order corrections to hydrodynamics Fagotti (2017). In the random case, we expect that these oscillations should vanish, and that the front will collapse onto a diffusive form. From Fig. 9, it is clear that there are no characteristic oscillations near the front. However, scaling our data with both $t^{1/2}$ and $t^{1/3}$ leads to equally good collapses, and suggest that such collapses are not a very conclusive way to measure the exponent $\alpha$ of the front broadening $t^{\alpha}$. We also considered collapses from the time evolution following local quenches, with very similar results. Nevertheless, our transport results combined with our hydrodynamic theory strongly suggest that $\alpha=1/2$ even for random non-interacting chains.

VI.2 Slow local relaxation

We now turn briefly to the part of the local operator that remains near its initial position at late times. Since we are considering a noninteracting model, we can equivalently consider the return probability of an initially local wavepacket Crowley et al. (2018). We should distinguish between average and typical behavior: typically the initial site is not a weak link, so anomalous wavefunctions have no support there. However, the site-averaged local autocorrelation function does receive a contribution from rare sites. Once again, we discuss this for the exponential distribution (4) in the quasi-localized phase, $\phi<2$.

We consider the probability that a particle initially localized at site $x$ has moved a distance less than one lattice site at time $t$. This quantity is proportional to the (mean) local autocorrelation function $C_{0}(t)=\overline{\langle S^{z}_{i}(t)S^{z}_{i}(0)\rangle}$. We focus on $i$ even and infinite temperature. In generalized hydrodynamics this can be expressed as

[TABLE]

where at infinite temperature for a free-fermion model, $\rho^{T}\sim\rho$, and $a=1$ is the lattice spacing. Focusing on low-energy quasiparticles, this integral can be written out as

[TABLE]

We now resolve the step function and approximate $\rho_{t}(\lambda,\xi)\simeq e^{-\pi\lambda/2+2\xi}\Theta(\lambda-4\xi/\pi)$, as in Sec. V.3, to rewrite this expression in terms of the double integral

[TABLE]

Higher energy quasiparticles give rise to a ballistic decay $1/t$ that is subleading when $\phi<2$. Thus, throughout the quasi-localized phase the autocorrelation function decays slower than one would expect for a model with ballistic transport. Numerical simulations of the autocorrelation function gives results in very good agreement with this exponent (Fig. 10). We emphasize that in the argument above, it was crucial to disorder-average the full autocorrelation function—separately averaging the velocity and the density of states would yield an incorrect exponent $\phi/(2-\phi)$ in clear disagreement with our numerical results. The local velocity is inversely proportional to the local density of states, and capturing these correlations is essential to deriving the correct anomalous exponent.

We will show elsewhere that anomalous decay of local autocorrelation functions occurs in other models such as XXZ as well; however, the possibility of subdiffusive behavior is specific to the disordered noninteracting models, as generic interacting models will have a finite diffusion constant due to interactions.

VII Discussion

We studied the non-equilibrium dynamics of integrable spin chains with correlated disorder that preserves integrability. Focusing on the non-interacting case, we formulated a (generalized) hydrodynamic theory for such random systems, and describe the emergence of diffusive corrections due to quasiparticles scattering off random impurities. This provides a mechanism for diffusion that is different from the recent theories of diffusive corrections to GHD in clean integrable quantum systems. The predictions from hydrodynamics were compared to numerical results obtained from exact diagonalization of the free fermion problem. Both spin and energy transport can be described very accurately using hydrodynamics, provided diffusive corrections are included. Moreover, we have shown that low-energy quasiparticles are very sensitive to the tails of the disorder distribution, and can become quasi-localized, leading to an anomalous decay of local autocorrelation functions.

We expect our results to generalize naturally to all interacting random integrable systems Vega and Woynarovich (1992); de Vega et al. (1994); Essler et al. (2018). In general, we expect a complicated interplay between diffusive corrections due to disorder, and due to thermal fluctuations and interactions. However, a simpler intermediate setup would be to consider initial states for which thermal fluctuations vanish, such as the spin domain-wall initial state considered above. For such initial states, we expect our predictions to extend naturally to the interacting case, and it would be interesting to compare the hydrodynamic predictions to matrix product state simulations.

Our results also indicate that diffusive broadening of the operator front can occur even in some non-interacting systems, which are clearly non-chaotic. These models could be used as a testbed for future diagnostic tools to distinguish chaotic from integrable systems. These models are also natural from the point of view of integrability breaking: adding integrability-breaking perturbations to a random integrable chain at strong disorder could lead to either thermalization or to many-body localization. The results of Ref. Essler et al., 2018 suggest that the regime of ballistic transport escaping localization might not be as fine-tuned as one could have expected, and it would be interesting to investigate whether hydrodynamics can still describe accurately transport away from the integrable limit.

Acknowledgments.— We thank Vir Bulchandani and Fabian Essler for useful comments on this manuscript. This work was supported by the US Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-SC0019168 (U.A. and R.V.), and by NSF Grant No. DMR-1653271 (S.G.). Most numerical simulations were performed at the Massachusetts Green High Performance Computing Center (MGHPCC).