Spin diffusion from an inhomogeneous quench in an integrable system

Generalized hydrodynamics predicts universal ballistic transport in integrable lattice systems when prepared in generic inhomogeneous initial states. However, the ballistic contribution to transport can vanish in systems with additional discrete symmetries. Here we perform large scale numerical simulations of spin dynamics in the anisotropic Heisenberg XXZ spin 1/2 chain starting from an inhomogeneous mixed initial state which is symmetric with respect to a combination of spin reversal and spatial reflection. In the isotropic and easy-axis regimes we find non-ballistic spin transport which we analyse in detail in terms of scaling exponents of the transported magnetization and scaling profiles of the spin density. While in the easy-axis regime we find accurate evidence of normal diffusion, the spin transport in the isotropic case is clearly super-diffusive, with the scaling exponent very close to 2/3, but with universal scaling dynamics which obeys the diffusion equation in nonlinearly scaled time.

I ntegrable models, such as the classical Kepler problem, harmonic oscillators, the planar Ising problem and so on, form cornerstones of our understanding of nature. Their equilibrium physics is usually well understood, even for the most complicated among integrable models, for example, the ones solvable by the Bethe ansatz 1 . Non-equilibrium physics of quantum systems on the other hand is much less understood 2 , particularly when going beyond the simplest integrability of quadratic models. This theoretical gap is becoming even more apparent with the advancement of experimental methods that are offering us analogue simulation of models beyond the capability of our best theoretical and numerical methods 3,4 .
Non-equilibrium dynamics of integrable quantum systems is thus one of the main current focuses of both theoretical and experimental condensed matter physics 5 . A macroscopic number of conservation laws existing in such systems 6 provide a variety of ways to break ergodicity, manifesting, for instance, in equilibration processes to non-thermal states or ballistic hightemperature transport of conserved quantities, such as energy, magnetization or charge. A naive classical reasoning might be that, because integrable systems are distinguished by constants of motion that force the dynamics to be simple and almost periodic (for example, orbits winding up the torus), one should expect to see ballistic transport. We shall demonstrate that this picture, while being correct for trivially integrable noninteracting models, such as harmonic oscillator chains 7 , can in fact be wrong for an interacting quantum integrable model.
Recently, a generalization of hydrodynamics has been put forward 8,9 which successfully predicts ballistic currents and scaled density profiles of integrable interacting systems quenched from inhomogeneous initial states [10][11][12][13][14][15] , which is a convenient method to study relaxation and non-equilibrium transport. In this protocol, the system is prepared in the state where the left and the right part, for xo0 and x40, respectively, are in different equilibrium states, and then, at t ¼ 0, let to evolve with a homogeneous interacting Hamiltonian. However, when ballistic transport is prohibited due to generic symmetries, such as is the case for spin transport in the anisotropic Heisenberg spin chain in the easy-axis (Ising) regime, this theory makes no prediction.
In extended interacting integrable system a macroscopic number of local conservation laws exists, in number proportional to the number of degrees of freedom, which can be exploited to develop generalized hydrodynamics 8 . This theory for typical inhomogeneous initial states predicts ballistic scaling f(x ¼ x/t) of densities and currents of conserved quantities, such as energy, charge or magnetization. However, in systems with parity Z 2 ð Þ symmetries, such as particle-hole exchange (or spin reversal), and for observables that are odd under the parity and initial states that are symmetric under the combined parity and spatial reflection x-À x, the ballistic contribution to transport can vanish. In fact, vanishing ballistic transport channel can then be related to the absence of local or quasi-local conserved charges with odd parity 6,16 . This means that the transported conserved quantity at x ¼ 0 grows slower than linear with t.
Here we propose a conjecture, based on large scale simulations, that a quench from an inhomogeneous initial state will in such cases generically result in diffusive spin dynamics. We demonstrate our results on the anisotropic Heisenberg chain (XXZ model). However, we stress that the XXZ model goes beyond being a mere toy model-it has been instrumental in the development of quantum integrability 17,18 and describes interaction in real spin chain materials 19 . Remarkably, in the case of isotropic Heisenberg interaction, spin relaxation is super-diffusive but with universal scaling dynamics which obey the standard diffusion equation in nonlinearly scaled time. Our results thus reveal a surprising property of an important integrable model as well as pose a challenge to theories which at present are unable to account for our observations. Because the parity symmetry is ubiquitous, our set-up should be widely applicable, for instance, we predict a similar physics in the one-dimensional (1D) Hubbard model.

Results
The set-up. The Hamiltonian of the XXZ chain of n sites reads where D is the anisotropy parameter and s x , ½H; M¼0. We are going to study the spin transport satisfying the continuity equation ds ð3Þ The existence of spin-reversal parity S¼ Q x s ð1Þ x , ½H; S¼0 and odd current j x S ¼ À Sj x , implies an absence of ballistic transport channels based on local conserved charges 20 . We are going to simulate the time evolution of an initial inhomogeneous state composed of two halves with opposite magnetizations.
To this end we choose a product initial state described by a density operator r, where the parameter mA[ À 1, 1] determines the initial magnetization, being hs ð3Þ x!0;o0 i¼ AE 1 2 m. Each of the initial halves can be thought of as being in equilibrium state $ e AE h P x s ð3Þ x at very high temperature and finite magnetization. We are therefore studying high-energy non-equilibrium physics of the model. While the initial state is pure for |m| ¼ 1 (a fully polarized domain wall), evolution of which has been studied in the past 21 , the choice of a mixed state offers several important advantages: it is generic and not plagued by the speciality of m ¼ 1 at D41 for which the dynamics freezes due to the proximity to a gapped eigenstate 22 , and it is, for small m, better suited for numerical simulations. This allows us to study significantly longer timescales as compared to existing literature and infer the scaling functions. We also mention that such an initial state can be thought of as representing an ensemble of pure states with randomized angle j on the Bloch sphere (Methods section).
Scaling exponents. We focus our efforts on DZ1 where there are no analytic results known for the magnetization transport, and the method 8,9 only predicts vanishing ballistic contribution. Two representative examples of a time evolved state r(t), namely the spin and current profiles sðx; tÞ¼ trrðtÞs ð3Þ x , j(x, t) ¼ trr(t)j x , are shown in Fig. 1. To obtain the exact type of transport we shall quantitatively study equilibration of magnetization, in particular the scaling of spin and current profiles as well as the transferred magnetization between the two halves, whose asymptotic scaling power a characterizes the transport type, where j(0, t) is the current at the half-cut. For a ¼ 1/2 the transport is diffusive, for 1/2oao1 it is called super-diffusive, and finally, a ¼ 1 corresponds to ballistic transport. We note that the transport type is connected to current-current correlation function via Green-Kubo linear response theory. In case of diffusive transport, the spin density satisfies the diffusion equation. This notion of diffusion does not necessarily correspond to De Gennes phenomenological theory of spin diffusion which, under much stronger assumptions, in 1D implies 1/ ffiffi t p dependence of local spin density autocorrelation function 23,24 .
We evolved the initial state r(0) (3) up to long times (of order tE160) and set large enough n so that there was no significant finite size effects. From the data we then infer the exponent a using equation (4), see Fig. 2a,b for representative plots. Dependence of the exponent a on D is summarized in Fig. 2c. While the transport is found to be ballistic for Do1, expectedly so for the integrable system, also known rigorously 16 , at DZ1 we find rather clear non-ballistic relaxation. In particular, at D ¼ 1 it is super-diffusive while for D41 the transport is diffusive, observed in driven steady-state setting 25,26 as well as in the Hamiltonian one 24,27-30 . At D ¼ 1 we also observe small dependence of a on m. While for small m, that is, small deviations from an infinite temperature state rB1, the exponent is close to 2/3, closer to pure state m ¼ 1 it appears to be closer to E3/5 (we note that a different numerical procedure is used in the two regimes, see Methods).
Scaling functions. The scaling of the transferred magnetization unequivocally shows a surprising non-ballistic transport in an integrable system which, however, has been observed and discussed before in related contexts, namely within local quench and linear response theory 24,27-30 and boundary driven Lindblad approach 25,26 . But here we can do still more. In Fig. 3 we demonstrate that the spin profiles can be described by a function of a single-scaling variable x/t a -profiles at large times collapse to a single curve. In addition, the profiles of current and magnetization are proportional to each other at different times (Fig. 3c,d), therefore validating Fick's law j ¼ À Drs where the behaviour of the diffusion constant D with respect to the anisotropy D is shown in the inset of Fig. 2c. This comes as no x Þ (a,b) and current (c,d) profile j(x, t) ¼ tr(r(t)j x ) for the isotropic point D ¼ 1 (a,c), and D ¼ 2 (b,d), following an inhomogeneous quench. One can see that the spreading is much faster for D ¼ 1, in both cases though it is slower than ballistic. Dashed green curves guide the eye towards scaling xBt 2/3 in a, and xBt 1/2 in (b). Data are shown for n ¼ 320 and small initial polarization m ¼ p/1,800.
surprise in the diffusive regime D41 where the scaling function of the magnetization (Fig. 3b) is simply the error function sðx; tÞ¼ À m 2 erf x= ffiffiffiffiffiffiffi ffi 4Dt p À Á . However, the same can not be said for the isotropic point D ¼ 1. Proportionality between the magnetization gradient and the current profile (Fig. 3c), this time with a time-dependent ratio D ' K 3 t 1=3 , suggests a diffusion equation in a scaled time @sðx; tÞ @t which again yields error function profile with a different scaling variable sðx; tÞ¼ À m 2 erf K À 1=2 x=t 2=3 À Á with K ¼ 2.33±0.03. In Fig. 3a we compare numerical profiles with the error function, again finding good agreement within accuracy of our simulations. Therefore, the scaling function is, in both cases, D ¼ 1 and D41, the error function, the difference being only in the scaling variable which is x/t 2/3 at the super-diffusive isotropic point. This result is surprising, as anomalous diffusion is usually associated with Levy processes and hence long (non-Gaussian) tails in the profiles. Here it seems it all amounts to a nonlinear rescaling of time. Theoretical explanation of this effect is urgent.
Entanglement entropy and simulation complexity. Finally, we mention a numerical observation that explains why we can simulate dynamics to such long times, and is an interesting property on its own. We use a time-dependent density matrix renormalization group method (tDMRG), see Methods.
The efficiency of tDMRG depends on the entanglement entropy, that is, for pure state evolution on the Von Neumann entropy S¼ À tr r A ln r A ½ of the reduced state r A ¼ tr A |CihC|, whereas for mixed states evolution on an analogous operator space entanglement entropy S # , (ref. 31) of a vectorised density operator r. When starting with a typical product initial state both entropies typically grow linearly with time, regardless of the system being integrable or not 32,33 , causing exponentially fast growth of complexity and with it a failure of these numerical methods. In our case though, see Fig. 4, entropies grow much slower, namely in a power-law fashion with b being o1. The most efficient simulations have been possible with density operators for small m where the exponent b is typically between 0.3 and 0.5.

Discussion
Our numerical results can be interpreted as an evidence of normal spin diffusion and spin Fick's law in the easy-axis anisotropic Heisenberg chain (for anisotropy D41), with spin density satisfying the diffusion equation on large scales. Besides the case D ¼ 2 shown here, we provide additional data for D ¼ 1.05, 1.1, 1.3, 1.5 demonstrating a clear convergence of the diffusive scaling exponents a ¼ 1/2 in all massive cases (Supplementary Note 1), and data for massless cases D ¼ 0, 0.5, 0.7, 0.9 which indicate convergence to ballistic exponent a ¼ 1 (Supplementary Note 2). While for generic, non-spin-reversal-symmetric initial states, the dominant contribution to transport is ballistic as determined by generalized hydrodynamics (or generalized 1D Euler's equations) 8,9,[13][14][15] , the next-to-leading term is now clearly predicted to be diffusive, as following from our work. However, a theoretical explanation, or even derivation of a diffusive contribution to transport in an integrable system with a macroscopic number of conservation laws is still pending. Even more surprising is the discovery of anomalous super-diffusive transport in the isotropic case (D ¼ 1) with the scaling exponent equal to or very close to 2/3. While this might suggest a behaviour described by KPZ (Kardar-Parisi-Zhang) universality class, we find that asymptotic spin density profiles obey the nonlinearly scaled diffusion equation and are distinct from the KPZ profiles. One might conjecture that the scaling exponent 2/3 is a consequence of SU(2) symmetry and not the fact that the model there corresponds to the marginal critical point D ¼ 1. This would be consistent with observed anomalous super-diffusive scalings in SU(4) spin ladders in the set-up of driven steady-state Lindblad dynamics 34 where the scaling exponent appears to be a ¼ 3/5. Curiously, all scaling exponents observed in this work (1/1, 1/2, 2/3, 3/5) are ratios of subsequent Fibonacci numbers 35 .

Methods
Numerical procedures. The time evolution is performed by means of the tDMRG algorithm 36,37 . In particular, for small m data (which is mostly reported here) the most efficient was the matrix product density operator version of tDMRG, with which we could reach times of the order tC200 for system size nC2t using bond dimensions 50-200 resulting in relative truncation errors o1%. One the other hand, for mE1 (close to domain wall pure state), the pure state version of tDMRG becomes more efficient as the corresponding entanglement entropy scaling exponents b are smaller. The two approaches appear to complement one another as can be seen in Fig. 2d. Neither approach allows us to observe long times in the intermediate region of m, where the exponents b become closer to 1. In order to simulate the desired density operator by evolving pure states we define a set of initial states where cðm; fÞ # j i is simply the Bloch sphere representation of a 2-level system and the f x are uniform independent random Scaling of density and current profiles with x/t a . In (a,b) we show the scaling of magnetization profiles, (a) for D ¼ 1 using a ¼ 2/3, and (b) for D ¼ 2 and using a ¼ 1/2 (note that the points for different times overlap almost perfectly; the insets show the convergence of the relative root-mean-square difference (in %) between data s(x, t) and scaled erf-profiles (see text) as a function of time). Frames (c,d) show the emergence of Fick's law at late times (shown at t ¼ 160), comparing current profiles (red) to gradients of spin density (blue)-both indistinguishable from Gaussians, for D ¼ 1 in (c) and D ¼ 2 in (d). In all plots the system size is n ¼ 320. numbers in the range [0, 2p). The density matrix is then obtained as an ensemble average over a set of such pure random states rðtÞ¼E CðtÞ j i CðtÞ h j ð Þ . It is clear that an increasingly large set of random states is needed as the magnetization approaches m-0, where the matrix product density operator simulation is favourable anyway.
Data availability. Data are available on request from the authors.