# Spin diffusion from an inhomogeneous quench in an integrable system

## Abstract

Generalised 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 spin 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 magnetisation 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 , but with universal scaling dynamics which obeys the diffusion equation in nonlinearly scaled time.

## I Introduction

Integrable models, such as the classical Kepler problem, harmonic oscillators, the planar Ising problem, etc., form cornerstones of our understanding of nature. Their equilibrium physics is usually well understood, even for the most complicated among integrable models, e.g. the ones solvable by the Bethe ansatz takahashi (). Nonequilibrium physics of quantum systems on the other hand is much less understood RMP (), 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 analog simulation of models beyond the capability of our best theoretical and numerical methods Bloch (); exper ().

Nonequilibrium dynamics of integrable quantum systems is thus one of the main current focuses of both theoretical and experimental condensed matter physics JSTAT (). A macroscopic number of conservation laws existing in such systems qloc_review () provide a variety of ways to break ergodicity, manifesting, for instance, in equilibration processes to non-thermal states or ballistic high-temperature transport of conserved quantities, such as energy, magnetisation 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 (e.g., 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 Lieb (), can in fact be wrong for an interacting quantum integrable model.

Recently, a generalisation of hydrodynamics has been put forward doyon (); bertini () which successfully predicts ballistic currents and scaled density profiles of integrable interacting systems quenched from inhomogeneous initial states ruelle (); bernard (); bhaseen (); enej (); moore (); doyonspohn (), which is a convenient method to study relaxation and nonequilibrium transport. In this protocol, the system is prepared in the state where the left and the right part, for and respectively, are in different equilibrium states, and then, at , 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 generalised hydrodynamics doyon (). This theory for typical inhomogeneous initial states predicts ballistic scaling of densities and currents of conserved quantities, such as energy, charge or magnetisation. However, in systems with parity () symmetries, such as particle-hole exchange (or spin reversal), observables that are odd under the parity and initial states that are symmetric under the combined parity and spatial reflection , 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 prosen (); qloc_review (). This means that the transported conserved quantity at grows slower than linear with .

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 ( model). However, we stress that the model goes beyond being a mere toy model – it has been instrumental in the development of quantum integrability Bethe (); Baxter () and describes interaction in real spin chain materials Hlubek (). 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 setup should be widely applicable, for instance, we predict a similar physics in the one-dimensional Hubbard model.

## Results

The setup.– The Hamiltonian of the chain of sites reads

(1) |

where is the anisotropy parameter and are the spin operators, with Cartesian component , expressed in term of Pauli matrices (we use units ). The Hamiltonian preserves the total magnetisation, , . We are going to study the spin transport satisfying the continuity equation with the current

(2) |

The existence of spin-reversal parity , , and odd current , implies an absence of ballistic transport channels based on local conserved charges zotos (). We are going to simulate the time evolution of an initial inhomogeneous state composed of two halves with opposite magnetisations.

To this end we choose a product initial state described by a density operator ,

(3) |

where the parameter determines the initial magnetisation, being . Each of the initial halves can be thought of as being in equilibrium state at very high temperature and finite magnetisation. We are therefore studying high-energy nonequilibrium physics of the model. While the initial state is pure for (a fully polarised domain wall), evolution of which has been studied in the past Gobert (), the choice of a mixed state offers several important advantages: it is generic and not plagued by the speciality of at for which the dynamics freezes due to the proximity to a gapped eigenstate Gochev (), and it is, for small , 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 randomised angle on the Bloch sphere (see Methods).

Scaling exponents.– We focus our efforts on where there are no analytic results known for the magnetisation transport, and the method doyon (); bertini () only predicts vanishing ballistic contribution. Two representative examples of a time evolved state , namely the spin and current profiles , , are shown in Fig. 1.

In order to obtain the exact type of transport we shall quantitatively study equilibration of magnetisation, in particular the scaling of spin and current profiles as well as the transferred magnetisation between the two halves, whose asymptotic scaling power characterizes the transport type,

(4) |

where is the current at the half-cut. For the transport is diffusive, for it is called superdiffusive, and finally, 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 one-dimension implies dependence of local spin density autocorrelation function affleck11 (); mccoy ().

We evolved the initial state (3) up to long times (of order ) and set large enough so that there was no significant finite size effects. From the data we then infer the exponent using Eq.(4), see Fig.2(a) and (b) for representative plots. Dependence of the exponent on is summarised in Fig. 2(c). While the transport is found to be ballistic for , expectedly so for the integrable system, also known rigorously prosen (), at we find rather clear non-ballistic relaxation. In particular, at it is superdiffusive while for the transport is diffusive, observed in driven steady-state setting jstat09 (); prl11 () as well as in the Hamiltonian one Robin09 (); FHM (); affleck11 (); robin14 (); robin17 (). At we also observe small dependence of on . While for small , i.e., small deviations from an infinite temperature state , the exponent is close to , closer to pure state it appears to be closer to (we note that a different numerical procedure is used in the two regimes, see Methods).

Scaling functions.– The scaling of the transferred magnetisation 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 Robin09 (); FHM (); affleck11 (); robin14 (); robin17 () and boundary driven Lindblad approach jstat09 (); prl11 (). But here we can do more still. In Fig. 3 we demonstrate that the spin profiles can be described by a function of a single scaling variable – profiles at large times collapse to a single curve. In addition, the profiles of current and magnetisation are proportional to each other at different times (Fig. 3(c,d)), therefore validating Fick’s law where the behaviour of the diffusion constant with respect to the anisotropy is shown in the inset of Fig. 2(c). This comes as no surprise in the diffusive regime where the scaling function of the magnetisation (Fig. 3(b)) is simply the error function . However, the same can not be said for the isotropic point . Proportionality between the magnetisation gradient and the current profile (Fig.3(c)), this time with a time-dependent ratio , suggests a diffusion equation in a scaled time

(5) |

which again yields error function profile with a different scaling variable with . In Fig.3(a) 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, and , the error function, the difference being only in the scaling variable which is at the superdiffusive 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.– Lastly, 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 renormalisation group method (tDMRG), see Methods. The efficiency of tDMRG depends on the entanglement entropy, i.e., for pure state evolution on the Von Neumann entropy of the reduced state , whereas for mixed states evolution on an analogous operator space entanglement entropy ProsenPizorn () of a vectorised density operator . When starting with a typical product initial state both entropies typically grow linearly with time, regardless of the system being integrable or not CalabreseCardy (); Chiara06 (), 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

(6) |

with being less than . The most efficient simulations have been possible with density operators for small where the exponent is typically between and .

## Ii 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 ), with spin-density satisfying the diffusion equation on large scales. Besides the case shown here, we provide additional data for demonstrating a clear convergence of the diffusive scaling exponents in all massive cases (Supplementary Note 1), and data for massless cases which indicate convergence to ballistic exponent (Supplementary Note 2). While for generic, non-spin-reversal-symmetric initial states, the dominant contribution to transport is ballistic as determined by generalised hydrodynamics (or generalised one-dimensional Euler’s equations) doyon (); bertini (); enej (); moore (); doyonspohn (), 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 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-ballistic transport in the isotropic case () with the scaling exponent equal to or very close to . While this might suggest a behaviour described by KPZ (Kardar-Parisi-Zhang) universality class, we find that asymptotic spin density profiles obey the non-linearly scaled diffusion equation and are distinct from the KPZ profiles. One might conjecture that the scaling exponent is a consequence of symmetry and not the fact that the model there corresponds to the marginal critical point . This would be consistent with observed anomalous super-diffusive scalings in spin ladders in the setup of driven steady state Lindblad dynamics spinladder () where the scaling exponent appears to be . Curiously, all scaling exponents observed in this work (, , , ) are ratios of subsequent Fibonacci numbers popkov ().

## Iii Methods

Numerical procedures.– The time evolution is performed by means of the tDMRG algorithm vidal (); schollwock (). In particular, for small 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 for system size using bond dimensions resulting in relative truncation errors less than . One the other hand, for (close to domain wall pure state), the pure state version of tDMRG becomes more efficient as the corresponding entanglement entropy scaling exponents are smaller. The two approaches appear to complement one another as can be seen in Fig.2(d). Neither approach allows us to observe long times in the intermediate region of , where the exponents become closer to .

In order to simulate the desired density operator by evolving pure states we define a set of initial states

(7) |

where is simply the Bloch sphere representation of a 2-level system and the are uniform independent random numbers in the range . The density matrix is then obtained as an ensemble average over a set of such pure random states . It is clear that an increasingly large set of random states is needed as the magnetisation approaches , where the matrix product density operator simulation is favourable anyway.

## Iv Acknowledgements

The authors would like to acknowledge support from ERC grant OMNES, as well as grants J1-7279, N1-0025, and P1-0044 of the Slovenian Research Agency.

## Appendix A Supplementary Note 1 - Diffusive Regime

We direct the readers attention to Supplementary Figure 5, where the local-time exponent is shown as a function of time for several values of the anisotropy in the massive regime, namely . The asymptotic scaling exponent converges on the accessible time scale to – expect for the case where the convergence has not yet been reached though the extrapolated asymptotic value is likely the same – is distinctly different from the exponent for the isotropic case .

## Appendix B Supplementary Note 2 - Ballistic Regime

Spin transport is known to be ballistic in the massless regime, . Here we verify that the chosen inhomogeneous quench does indeed reproduce this result. In Supplementary Figure 6 we demonstrate the convergence of local-time scaling exponent for a few cases to an asymptotic ballistic exponent , and again compare it to the isotropic case .

In Supplementary Figure 7 we also show the scaling of spin and current density profiles which clearly exhibit expected ballistic behaviour doyon (); bertini (). In the non-interacting case we find also excellent agreement with analytic solutions antal ().

### References

- M. Takahashi, Thermodynamics of One-Dimensional Solvable Models, Cambridge University Press (1999).
- A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Colloquium: Nonequilibrium dynamics of closed interacting quantum systems, Rev. Mod. Phys. 83, 863-883 (2011).
- I. Bloch, J. Dalibard and S. Nascimbene, Quantum simulations with ultracold quantum gases, Nature Physics 8, 267-276 (2012).
- J.-Y. Choi, et. al., Exploring the many-body localization transition in two dimensions, Science 352, 1547-1552 (2016).
- P. Calabrese, F. H. L. Essler and G. Mussardo, Special Issue on ‘Quantum Integrability in Out of Equilibrium Systems’, J. Stat. Mech. (2016), 064001 (2016).
- See, e.g. a review: E. Ilievski, M. Medenjak, T. Prosen and L. Zadnik, Quasilocal charges in integrable lattice systems, J. Stat. Mech. (2016), 064008 (2016).
- Z. Rieder, J. L. Lebowitz, and E. Lieb, Properties of a harmonic crystal in a stationary nonequilibrium state, J. Math. Phys. 8, 1073-1078 (1967).
- O. A. Castro-Alvaredo, B. Doyon and T. Yoshimura, Emergent hydrodynamics in integrable quantum systems out of equilibrium, Phys. Rev. X 6, 041065 (2016).
- B. Bertini, M. Collura, J. De Nardis, M. Fagotti, Transport in out-of-equilibrium XXZ chains: exact profiles of charges and currents, Phys. Rev. Lett. 117, 207201 (2016).
- D. Ruelle, Natural Nonequilibrium States in Quantum Statistical Mechanics, J. Stat. Phys. 98, 57-75 (2000).
- D. Bernard and B. Doyon, Non-Equilibrium Steady-States in Conformal Field Theory, Ann. Inst. Henri Poincaré 16, 113-161 (2015).
- J. Bhaseen, B. Doyon, A. Lucas and K. Schalm, Energy flow in quantum critical systems far from equilibrium, Nature Physics 11, 509-514 (2015).
- E. Ilievski and J. De Nardis, On the Microscopic Origin of Ideal Conductivity, Preprint at http://arXiv.org/abs/1702.02930 (2017).
- V. B. Bulchandani, R. Vasseur, C. Karrasch and J. E. Moore, Bethe-Boltzmann Hydrodynamics and Spin Transport in the XXZ Chain, Preprint at http://arXiv.org/abs/1702.06146 (2017).
- B. Doyon, H. Spohn and T. Yoshimura, A geometric viewpoint on generalized hydrodynamics, Preprint at http://arXiv.org/abs/1704.04409 (2017).
- T. Prosen, Open XXZ Spin Chain: Nonequilibrium Steady State and a Strict Bound on Ballistic Transport, Phys. Rev. Lett. 106, 217206 (2011).
- H. Bethe, On the Theory of Metals, I. Eigenvalues and Eignefunctions of a Linear Chain of Atoms, Zeits. Physik 74, 205-226 (1931).
- R. J. Baxter, Exactly Solved Models in Statistical Mechanics, (1982).
- N. Hlubek, X. Zotos, S. Singh, R. Saint-Martin, A. Revcolevschi, B. Büchner, and C. Hess, Spinon heat transport and spin–phonon interaction in the spin-1/2 Heisenberg chain cuprates and , J. Stat. Mech. 2012, P03006 (2012).
- X. Zotos, F. Naef, and P. Prelovšek, Transport and conservation laws, Phys. Rev. B 55, 11029-11032 (1997).
- D. Gobert, C. Kollath, U. Schollwöck, and G. Schütz, Real-time dynamics in spin-1/2 chains with adaptive time-dependent density matrix renormalization group, Phys. Rev. E 71, 036102 (2005).
- I. G. Gochev, Contribution to the theory of plane domain walls in a ferromagnet, Sov. Phys. JETP 58, 115-119 (1983).
- K. Fabricius and B. M. McCoy, Spin diffusion and the spin-1/2 XXZ chain at from exact diagonalization, Phys. Rev. B 57, 8340-8347 (1998).
- J. Sirker, R. G. Pereira, and I. Affleck, Conservation laws, integrability, and transport in one-dimensional quantum systems, Phys. Rev. B 83, 035115 (2011).
- T. Prosen and M. Žnidarič, Matrix product simulation of non-equilibrium steady states of quantum spin chains, J. Stat. Mech. 2009, P02035 (2009).
- M. Žnidarič, Spin transport in a one-dimensional anisotropic Heisenberg model, Phys. Rev. Lett. 106, 220601 (2011).
- R. Steinigeweg and J. Gemmer, Density dynamics in translationally invariant spin-1/2 chains at high temperatures: A current-autocorrelation approach to finite time and length scales, Phys. Rev. B 80, 184402 (2009).
- C. Karrasch, J. E. Moore, and F. Heidrich-Meisner, Real-time and real-space spin and energy dynamics in one-dimensional spin-1/2 systems induced by local quantum quenches at finite temperatures, Phys. Rev. B 89, 075139 (2014).
- R. Steinigeweg, J. Gemmer, W. Brenig, Spin-current autocorrelations from single pure-state propagation, Phys. Rev. Lett. 112, 120601 (2014).
- R. Steinigeweg, F. Jin, D. Schmidtke, H. De Raedt, K. Michielsen, and J. Gemmer, Real-time broadening of nonequilibrium density profiles and the role of the specific initial-state realization, Phys. Rev. B 95, 035155 (2017).
- T. Prosen and I. Pižorn, Operator space entanglement entropy in a transverse Ising chain, Phys. Rev. A 76, 032316 (2007).
- P. Calabrese and J. Cardy, Evolution of entanglement entropy in one-dimensional systems, J. Stat. Mech. 2005, P04010 (2005).
- G. De Chiara, S. Montangero, P. Calabrese,and R. Fazio, Entanglement entropy dynamics of Heisenberg chains, J. Stat. Mech. 2006, P03001 (2006).
- M. Žnidarič, Magnetization transport in spin ladders and next-nearest-neighbor chains, Phys. Rev. B 88, 205135 (2013).
- V. Popkov, A. Schadschneider, J. Schmidt, G. M. Schütz, Fibonacci family of dynamical universality classes, Proc. Natl. Acad. Sci. U.S.A. 112, 12645-12650 (2015).
- G. Vidal, Efficient simulation of one-dimensional quantum many-body systems, Phys. Rev. Lett. 93, 040502 (2004).
- U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Ann. Phys. 326, 96-192 (2011).
- T. Antal, Z. Racz, A. Rakos, G. M. Schütz, Transport in the XX chain at zero temperature: Emergence of flat magnetization profiles, Phys. Rev. E 59, 4912-4918 (1999).