Non-local emergent hydrodynamics in a long-range quantum spin system
NNon-local emergent hydrodynamics in a long-range quantum spin system
Alexander Schuckert,
1, 2, ∗ Izabella Lovas,
1, 2, † and Michael Knap
1, 2, ‡ Department of Physics and Institute for Advanced Study,Technical University of Munich, 85748 Garching, Germany Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4, D-80799 M¨unchen (Dated: September 5, 2019)Generic short-range interacting quantum systems with a conserved quantity exhibit universaldiffusive transport at late times. We show how this universality is replaced by a more generaltransport process in the presence of long-range couplings that decay algebraically with distanceas r − α . While diffusion is recovered for large exponents α > .
5, longer-ranged couplings with0 . < α ≤ . t − / (2 α − when 0 . < α ≤ .
5. We also extract the associated generalizeddiffusion constant, and demonstrate that it follows the prediction of classical L´evy flights; quantummany-body effects manifest themselves in an overall time scale depending only weakly on α . Ourfindings can be readily verified with current trapped ion experiments. In quantum many-body systems, macroscopic inhomo-geneities in a conserved quantity must be transportedacross the whole system to reach an equilibrium state.As this is in general a slow process compared to local de-phasing, essentially classical hydrodynamics is expectedto emerge at late times in the absence of long-lived quasi-particle excitations [1–8]. The universality of this effec-tive classical description may be understood from thecentral limit theorem: in the regime of incoherent trans-port, short range interactions lead to an effective ran-dom walk with a finite variance of step sizes, leading to aGaussian distribution at late times. This universality isonly broken when quantum coherence is retained, such asin integrable models [9–13] or in the vicinity of a many-body localized phase, where rare region effects give riseto subdiffusive transport [14–20].In this work, we show how this universal diffusivetransport in short range interacting systems is replacedby a more general, non-local effective hydrodynami-cal description in systems with algebraically decayinglong-range interactions. We use semi-analytical non-equilibrium quantum field theory calculations (referredto as spin-2PI below) and a discrete truncated Wignerapproximation (dTWA) to show that in a long rangeinteracting XY spin chain, spin transport at infinitetemperature effectively obeys a classical master equa-tion with long-range , algebraically decaying transitionamplitudes. This effective description can be reformu-lated as a classical random walk with infinite variance ofstep sizes, giving rise to a generalized central limit the-orem and to a late-time description in terms of classicalL´evy flights [21], an example for super diffusive anoma-lous transport. As a result, we demonstrate that thefull spatio-temporal shape of the correlation function C ( j, t ) = (cid:104) ˆ S zj ( t ) ˆ S z (0) (cid:105) , and in particular, the exponent FIG. 1.
Hydrodynamic tails in the spin autocorrelator. (a) For long range coupling exponents α > .
5, autocorrela-tions decay algebraically at late times with an exponent thatdepends on α . By contrast for α ≤ . β α of the hydrodynamic tailobtained from two different approaches (symbols) agree withthe predictions from classical L´evy flights in the thermody-namic limit (dashed curve). Deviations at large α are due tofinite time corrections to scaling which can also be understoodfrom L´evy flights. of the hydrodynamic tail in the autocorrelation function C ( j = 0 , t ), depends strongly on the long-range exponent a r X i v : . [ c ond - m a t . qu a n t - g a s ] S e p α . While for α > . / (2 α −
1) for 0 . < α ≤ . C ( j, t ) possesses a self-similar be-havior, with the scaling function covering all stable sym-metric distributions as a function of α , smoothly cross-ing over from a Gaussian at α = 1 . α = 1 to an even more sharply peaked function as α → .
5. We also extract the generalized diffusion coef-ficient D α from the scaling functions, and explain its α dependence by L´evy flights; quantum effects are incorpo-rated in a many-body time scale depending only weaklyon α . For α ≤ . Model.–
We study the long range interacting quantumXY chain with open boundary conditions, given by theHamiltonianˆ H = − L/ (cid:88) i (cid:54) = j = − L/ J N L,α | i − j | α (cid:16) ˆ S xi ˆ S xj + ˆ S yi ˆ S yj (cid:17) . (1)Here, ˆ S α = ˆ σ α denotes spin- operators given in termsof Pauli matrices, L is the (odd) length of the chain [23],and we set (cid:126) = 1. The interaction strength J is rescaledwith the factor N L,α = (cid:113)(cid:80) j (cid:54) =0 | j | − α in order to re-move the L and α dependence of the time scale associ-ated with the perturbative short time dynamics of thecentral spin at i = 0. The above model shows chaotic(Wigner-Dyson) level statistics for the whole range of α considered here (0 . ≤ α ≤
2) and is an effective de-scription of the long range transverse field Ising modelfor large fields [24, 25]. In particular, it conserves the to-tal S z magnetization, with product states in the S z basisevolving radically differently depending on the complex-ity of the corresponding magnetization sector. For justa few spin flips on top of the completely polarized state,the dynamics can be exactly solved and are described interms of ballistically propagating spin waves, with a di-verging group velocity at α = 1 [24, 26, 27] related to thealgebraic leakage of the Lieb-Robinson bound [28–31].In contrast, here we show that the exponentially largeHilbert space sector for an extensive number of spin flipsgives rise to rich transport phenomena, driven by thelong-range nature of the interactions. Effective stochastic description of long-rangetransport.–
As the model in Eq. 1 is equivalent to long-range hopping hard core bosons, we conjecture the effec-tive classical equation of motion for the transported local density f j ( t ), in our case (cid:68) ˆ S zj ( t ) (cid:69) + , to be of the form[32] ∂ t f j ( t ) = (cid:88) i (cid:54) = j ( W i → j f i − W j → i f j ) . (2)Here, the transition rate W i → j is determined by the mi-croscopic transport processes present in the Hamiltonian,in our case the long-range hopping of spins. More specif-ically, from Fermi’s golden rule the transition rate for aflip flop process between spins i and j is proportional to |(cid:104)↑ i ↓ j | ˆ H | ↓ i ↑ j (cid:105)| , and hence we phenomenologically set W i → j = W j → i = λ | i − j | α , (3)where λ − is a characteristic time scale determined bythe full many-body Hamiltonian.Starting from an initial state with a single excitationin the center of the chain, the solution of this Masterequation is given by [33] f j ( t ) ≈ (cid:40) ( D α t ) − / G (cid:16) | j | ( D α t ) / (cid:17) for α > . D α t ) − β α F α (cid:16) | j | ( D α t ) βα (cid:17) for 0 . < α ≤ . G ( y ) = exp( − y / / √ π denotes the Gaussian distribu-tion, indicating normal diffusion for α > . D α ∝ λ . For 0 . < α ≤ . G ( y ) isreplaced by the family of stable, symmetric distributions F α ( y ), given by F α ( y ) = 14 (cid:90) d k π exp( −| k | /β α ) exp( iyk ) , (5)with the constant prefactor D α = λc α constituting ageneralized diffusion coefficient [34]. We find c α = − − α ) sin( πα ) from the classical Master equation,with Γ denoting the gamma function [33]. Furthermore,the exponent of the hydrodynamic tail β α is given by β α = 12 α − . The Fourier transform in Eq. (5) only leads to elemen-tary functions for α = 3 / α = 1, resulting in aGaussian and a Lorentzian distribution, respectively [35].The scaling functions F α ( y ) are the fixed point distribu-tions in the generalized central limit theorem [36] of i.i.d.random variables with heavy tailed distributions. Impor-tantly, F α ( y ) has diverging variance for α < .
5, unde-fined mean for α ≤
1, and displays heavy tails ∼ | y | − α .The classical Master equation hence predicts a cross-overfrom diffusive ( α ≤ .
5) over ballistic ( α = 1) to super-ballistic (0 . < α <
1) transport.
Quantum dynamics from spin-2PI and dTWA.–
In the following, we demonstrate the emergence of theseeffective classical dynamics in the quantum dynamics ofthe Hamiltonian (1), by studying the unequal time cor-relation function C ( j, t ) := Tr (cid:104) ˆ S zj ( t ) ˆ S z (0) (cid:105) | j=0 (cid:105) = |↑(cid:105) . (6)Here, we perform the trace over product states in the S z basis, restricted to the Hilbert space sector with (cid:80) i S zi = , such that (cid:104) S zi ( t = 0) (cid:105) = δ ,i for all spins i . Thisway, we probe the transport of a single spin excitationmoving in an infinite temperature bath with vanishingtotal magnetization.We employ two complementary, approximate meth-ods to study the dynamics at long times and for largesystem sizes, in a regime that is challenging to accessby numerically exact methods [37]. Schwinger bosonspin-2PI [38, 39], a non-equilibrium quantum field the-ory method, employs an expansion in the inverse co-ordination number 1 /z to reduce the many-body prob-lem to solving an integro-differential equation that scalesalgebraically in system size. As the effective coordina-tion number is large in a long-range interacting system,we expect this approximation to be valid for small α .The discrete truncated Wigner approximation evolvesthe classical equations of motion, while introducing quan-tum fluctuations by sampling initial states from theWigner distribution [40–43] and was shown to be par-ticularly well suited for studying long-range interactingsystems [42, 44]. In both methods, we evaluate C ( j, t ) bystarting from random initial product states in the S z ba-sis and then averaging over sufficiently many such initialstates [45]. If not stated otherwise, all our results havebeen converged with respect to system size, for which weemployed chains with 201 −
601 sites.We study two distinct regimes in the dynamics. Aperturbative short time regime characterized by initialdephasing is followed by the emergent effective classicallong-range transport described by the Master equation.
Perturbative short time dynamics.–
At shorttimes, second-order perturbation theory yieldsTr (cid:104) ˆ S zj ( t ) ˆ S z (0) (cid:105) ≈ (cid:26) (1 − J t ) for j = 0 (cid:16) Jt N L,α (cid:17) | j | α for j (cid:54) = 0 . (7)Physically, in this regime each spin is precessing in theeffective magnetic field created by all other spins. Theautocorrelation function is independent of α and L dueto our choice of the normalization factor N L,α , ensuringthat the typical magnetic field at the center of the chainremains of the order of J . The spatial correlation func-tion at a fixed time inherits the algebraic behavior of theinteraction strength, falling off as | j | − α between spinsof distance j , which is reproduced by both dTWA (notshown) and spin-2PI, see Fig. 2. Hydrodynamic tails.–
The scaling form from classi-cal L´evy flights in Eq. 4 implies the presence of a hydro-dynamic tail in the autocorrelation function C ( j = 0 , t ) FIG. 2.
Short time dynamics.
We compare spin-2PI re-sults with second order perturbation theory, Eq. 7. (a) Thecollapse of the autocorrelator for different exponents α showsthat the short time evolution is independent of α and L whenthe Hamiltonian is rescaled with N α,L . (b) The un-equal-time correlation function for α ∈ { . , , . , } (from top tobottom), shows algebraic tails that are entirely captured bysecond order perturbation theory. We used a moving averageover 5 −
10 lattice sites to smoothen the results. with exponent β α = 1 / (2 α − / α → . β α , these can however be explained by a subtlefinite-time effect also present in classical L´evy flights [33].For α < . L < N L →∞ ,α = ∞ , for α ≤ .
5. On even longer time scales,the discretized Fourier transform underlying the deriva-tion of Eq. 4 is dominated by the smallest wavenumberin finite chains, and the hydrodynamic tail is replacedby an exponential convergence towards the equilibriumvalue 0 . /L with a rate ∼ (1 /L ) α − . Self-similar time evolution of correlations.–
InFig. 3 we show the spreading of t β α C ( j, t ) for two val-ues of α . While for α = 2 a diffusive cone is visible, thespreading for α = 1 is ballistic as expected from the Mas-ter equation. The scaling collapse of these data showsgood agreement with classical L´evy flights, Eq. 4, at latetimes. Interestingly, we find heavy tails even for α ≥ . α > . α = 1 . α (cid:38) α increases. Thesepeaks are remnants of the integrable point at α = ∞ .Such behaviour is not present in the spin-2PI data as this FIG. 3.
Emergent self-similar time evolution.
The correlation function C ( j, t ) obtained from spin-2PI for chains of lengths L = 201 ( α = 2, subfigures (a-d)), L = 301 ( α = 1, subfigures (e-h)). (a,e) C ( j, t ) multiplied by t / (2 α − to account for theoverall decay expected from L´evy flights shows a diffusive cone for α = 2, whereas for α = 1 a ballistic light-cone emerges.The contour lines for α = 1 , t / (2 α − C ( j, t ) = 0 . , − , respectively. (b,f) Rescaling of linearly spacedtime slices for 23 ≤ Jt ≤
84 ( α = 1) and 42 ≤ Jt ≤
226 ( α = 2) (lines become darker as time increases) for the same dataas in (a,e) agrees well with the scaling function expected from classical L´evy flights, Eq. 4. The only fitting parameter is thegeneralized diffusion coefficient. (d,h) Rescaled time slices (2 ≤ Jt ≤
28) on a double-logarithmic scale reveal for α = 1 theheavy tail ∼ y − expected from L´evy flights (Eq. 4), where the dashed-dotted line is the same fit as in (b). The tail ∼ y − (thick black line) for α = 2 (8 ≤ Jt ≤
85) is a finite time effect also present in classical L´evy flights. (c,g) Unscaled data. method is not able to capture integrable behaviour [39].
Generalized diffusion constant.–
The only free pa-rameter of our effective classical description is the gen-eralized diffusion coefficient D α , which we obtain fromthe fits to the scaling function [47]. In Fig. 4 we showthat the leading α dependence of D α can be explainedby D α ∼ c α for α < . λ − , FIG. 4.
Generalized diffusion constant.
The α depen-dence of the diffusion constant obtained from fits with thescaling function of L´evy flights, Eq. 4. The qualitative behav-ior follows the L´evy flight prediction D α ∼ c α for α < . constituting the quantum many-body time scale, dependsonly weakly on α . As expected from their differing ap-proximate treatment of the quantum fluctuations in thesystem, we find slight differences between the values of λ determined by spin-2PI and dTWA, λ ≈ .
25 and λ dTWA ≈ .
15. For α > . Conclusions.–
In this paper, we have shown that spintransport at high temperatures in long-range interactingXY-chains is well described by L´evy flights for long-rangeinteraction exponents 0 . < α ≤ .
5, effectively realizinga random walk with infinite variance of step sizes. In par-ticular, we have shown that the scaling function of theunequal time spin correlation function covers the stablesymmetric distributions, in accordance with the gener-alized central limit theorem. While the system relaxesinstantly for α < .
5, standard diffusion was recoveredfor α > .
5, with heavy tails from finite time correc-tions surviving until extremely long times. We demon-strated the non-trivial dependence of the generalized dif-fusion coefficient D α on α , and found that it is capturedby classical L´evy flights, with the quantum many bodytime scale being approximately independent of α . Whilewe only studied one-dimensional systems, we expect thisphenomenon to generalize straightforwardly to d > d/ < α < d/ d/ (2 α − d ) [33].The long-range transport process found here can beexperimentally studied in current trapped ion experi-ments [49], which can reach the required time scales [50–53]. The effective infinite temperature states can also berealized by sampling over random product states whichare then evolved in time. Acknowledgments.–
We thank Rainer Blatt, EleanorCrane, Iliya Esin, Philipp Hauke, Christine Maier, AsierPi˜neiro Orioli, Tibor Rakovszky for insightful discus-sions and the Nanosystems Initiative Munich (NIM)funded by the German Excellence Initiative for ac-cess to their computational resources. We acknowl-edge support from the Max Planck Gesellschaft (MPG)through the International Max Planck Research Schoolfor Quantum Science and Technology (IMPRS-QST),the Technical University of Munich - Institute for Ad-vanced Study, funded by the German Excellence Initia-tive, the European Union FP7 under grant agreement291763, the Deutsche Forschungsgemeinschaft (DFG,German Research Foundation) under Germany’s Excel-lence Strategy–EXC-2111–390814868, the European Re-search Council (ERC) under starting grant ConsQuan-Dyn and the European Union’s Horizon 2020 researchand innovation program grant agreement no. 77153, fromDFG grant No. KN1254/1-1, and DFG TRR80 (ProjectF8). ∗ [email protected] † [email protected] ‡ [email protected][1] S. Mukerjee, V. Oganesyan, and D. Huse, Phys. Rev. B , 035113 (2006).[2] J. Lux, J. M¨uller, A. Mitra, and A. Rosch, PhysicalReview A , 053608 (2014).[3] A. Bohrdt, C. B. Mendl, M. Endres, and M. Knap, NewJournal of Physics , 063001 (2017).[4] M. Medenjak, K. Klobas, and T. Prosen, Physical Re-view Letters , 110603 (2017).[5] E. Leviatan, F. Pollmann, J. H. Bardarson, D. A. Huse,and E. Altman, arXiv:1702.08894 (2017).[6] V. Khemani, A. Vishwanath, and D. A. Huse, Phys. Rev.X , 031057 (2018).[7] T. Rakovszky, F. Pollmann, and C. W. von Keyserlingk,Phys. Rev. X , 031058 (2018).[8] D. E. Parker, X. Cao, A. Avdoshkin, T. Scaffidi, andE. Altman, arXiv:1812.08657 (2018).[9] S. Gopalakrishnan and R. Vasseur, Physical Review Let-ters , 127202 (2019).[10] M. Ljubotina, L. Zadnik, and T. Prosen, Phys. Rev. Lett. , 150605 (2019).[11] M. Ljubotina, M. ˇZnidariˇc, and T. Prosen, Phys. Rev.Lett. , 210602 (2019).[12] P. Prelovˇsek, J. Bonˇca, and M. Mierzejewski, Phys. Rev.B , 125119 (2018).[13] M. ´Sroda, P. Prelovˇsek, and M. Mierzejewski, Phys. Rev.B , 121110 (2019).[14] K. Agarwal, S. Gopalakrishnan, M. Knap, M. M¨uller,and E. Demler, Physical Review Letters , 160401(2015).[15] Y. Bar Lev, G. Cohen, and D. R. Reichman, Phys. Rev.Lett. , 100601 (2015).[16] R. Vosk, D. A. Huse, and E. Altman, Phys. Rev. X ,031032 (2015).[17] A. C. Potter, R. Vasseur, and S. A. Parameswaran, Phys.Rev. X , 031033 (2015).[18] P. Bordia, H. L¨uschen, S. Scherg, S. Gopalakrishnan,M. Knap, U. Schneider, and I. Bloch, Physical ReviewX , 041047 (2017).[19] S. Gopalakrishnan, K. R. Islam, and M. Knap, PhysicalReview Letters , 046601 (2017).[20] K. Agarwal, E. Altman, E. Demler, S. Gopalakrishnan,D. A. Huse, and M. Knap, Annalen der Physik ,1600326 (2017).[21] V. Zaburdaev, S. Denisov, and J. Klafter, Reviews ofModern Physics , 483 (2015).[22] R. Bachelard and M. Kastner, Phys. Rev. Lett. ,170603 (2013).[23] We always assume integer divisions when we write L .[24] P. Jurcevic, B. P. Lanyon, P. Hauke, C. Hempel, P. Zoller,R. Blatt, and C. F. Roos, Nature , 202 (2014).[25] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess,P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, NaturePhysics , 907 (2016).[26] P. Hauke and L. Tagliacozzo, Physical Review Letters , 20720 (2013).[27] P. Richerme, Z.-X. Gong, A. Lee, C. Senko, J. Smith,M. Foss-Feig, S. Michalakis, A. V. Gorshkov, andC. Monroe, Nature , 198 (2014).[28] E. H. Lieb and D. W. Robinson, Communications inMathematical Physics , 251 (1972).[29] M. B. Hastings and T. Koma, Communications in Math-ematical Physics , 781 (2006).[30] J. Eisert, M. Van Den Worm, S. R. Manmana, andM. Kastner, Physical Review Letters , 260401 (2013).[31] M. C. Tran, A. Y. Guo, Y. Su, J. R. Garrison, Z. El-dredge, M. Foss-Feig, A. M. Childs, and A. V. Gorshkov,Phys. Rev. X , 031006 (2019).[32] A Master equation in terms of a local probability densitymay be obtained by normalizing the local density.[33] See supplementary material.[34] The prefactor 1/4 accounts for the normalization of thecorrelation function, C ( j, t = 0) = δ ,j /
4. Furthermore,reinstating a lattice spacing a , the units of D α dependon α , in particular it is a velocity for α = 1.[35] Other closed form solutions exist [54], for example for α = 1 .
25 in terms of hypergeometric functions.[36] B. Gnedenko and A. N. Kolmogorov,
Limit Distributionsfor Sums of Independent Random Variables (Addison-Wesley, 1954).[37] B. Kloss and Y. Bar Lev, Phys. Rev. A , 032114 (2019).[38] M. Babadi, E. Demler, and M. Knap, Phys. Rev. X ,041005 (2015). [39] A. Schuckert, A. Pi˜neiro Orioli, and J. Berges, Phys.Rev. B , 224304 (2018).[40] W. K. Wootters, Annals of Physics , 1 (1987).[41] A. Polkovnikov, Annals of Physics , 1790 (2010).[42] J. Schachenmayer, A. Pikovski, and A. M. Rey, Phys.Rev. X , 011022 (2015).[43] L. Pucci, A. Roy, and M. Kastner, Phys. Rev. B ,174302 (2016).[44] A. Pi˜neiro Orioli, A. Signoles, H. Wildhagen, G. G¨unter,J. Berges, S. Whitlock, and M. Weidem¨uller, Phys. Rev.Lett. , 063601 (2018).[45] In spin-2PI simulations we used 4 −
16 different initialstates, while in dTWA the averaging over initial statesis performed in parallel with the Monte Carlo averagingover the Wigner distribution; here we typically use ∼ samples.[46] L. Zarfaty, A. Peletskyi, E. Barkai, and S. Denisov,arXiv:1908.02053 (2019).[47] We checked that for α > . D α = lim t →∞ ∂ t ( (cid:80) j j C ( j, t )),gives the same result, i.e. the algebraic tails do not alterthe result. [48] The spurious divergence of c α in the limit α = 1 . N L,α , D ( α ) ∼ N L,α → ∞ for α → . , 207901(2004).[50] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch,C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F.Roos, Science , 260 (2019).[51] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis,P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, andC. Monroe, Nature , 601 (2017).[52] N. Friis, O. Marty, C. Maier, C. Hempel, M. Holz¨apfel,P. Jurcevic, M. B. Plenio, M. Huber, C. Roos, R. Blatt,and B. Lanyon, Phys. Rev. X , 021012 (2018).[53] C. Maier, T. Brydges, P. Jurcevic, N. Trautmann,C. Hempel, B. P. Lanyon, P. Hauke, R. Blatt, and C. F.Roos, Phys. Rev. Lett. , 050501 (2019).[54] W. H. Lee, Continuous and Discrete Properties ofStochastic Processes , Ph.D. thesis, University of Notting-ham (2010).
SUPPLEMENTARY MATERIALScaling functions from the classical master equation
Here, we derive the scaling collapse given by Eqs. (4) and (5) from the master equation Eq. (2), ∂ t f j = (cid:88) i (cid:54) = j W ij ( f i − f j ) , (8)with transition rates W ij = λ | i − j | α . (9)By taking the Fourier transform of this equation, we arrive at ∂ t f ( k ) = [ W ( k ) − W ( k = 0)] f ( k ) , (10)where f ( k ) = L/ (cid:88) j = − L/ e − ikj f j , (11)and W ( k ) − W ( k = 0) = λ − (cid:88) j = − L/ + L/ (cid:88) j =1 (cid:0) e − ikj − (cid:1) / | j | α ≈ λ (cid:90) L/ dx (cos kx − /x α , (12)with k = 2 πn/L , n = − L/ , ..., L/ k (cid:28)
1. For 0 . < α < . k (cid:28) W ( k ) − W ( k = 0) ≈ λ (cid:90) L/ dx (cos kx − /x α − λ (cid:90) dx (cos kx − /x α ≈ λ (cid:90) ∞ dx (cos kx − /x α + λ k (cid:90) dx x − α = (cid:18) − c α | k | α − + k − α (cid:19) λ. (13)Here c α = − (cid:90) ∞ dz (cos z − /z α = − − α ) sin( απ ) , with Γ denoting the gamma function. Note that c α > . < α < . . < α < .
5, the first term in Eq. (13), ∼ | k | α − , will dominate the long time behavior, leading to ∂ t f ( k ) ≈ − λc α | k | α − ⇒ f ( k, t ) = f ( k, e − λc α | k | α − t . (14)In particular, taking an initial state with a single localized excitation, f j ( t = 0) = δ j, / f ( k, ≡ / / S z ) = 1 /
4, we arrive at the scaling ansatz f j ( t ) ≈ (cid:90) dk π exp (cid:0) ikj − λ t c α | k | α − (cid:1) = ( λc α t ) − / (2 α − F α (cid:18) | j | ( λc α t ) / (2 α − (cid:19) , (15)with F α ( y ) given by Eq. (5). Diffusion for α > . .– While we used α < . α > .
5. This can be shown by evaluating the following integral exactly, W ( k ) − W ( k = 0) ≈ λ | k | α − (cid:90) ∞| k | dz (cos z − /z α , and expanding the resulting expression around k = 0. Noting that for α > . | k | α − term is subdominant, wearrive at f α> . j ( t ) ≈ (cid:90) dk π exp (cid:18) ikj − λ α − k t (cid:19) = ( D α t ) − / G (cid:32) | j | ( D α t ) / (cid:33) , (16)reproducing diffusive behaviour for α > . D α = λ/ (2 α −
3) and a Gaussian G ( y ) =exp( − y / / √ π . We note that in contrast to the regime α < .
5, here the dTWA results for the quantum many-body time scale λ show a strong α dependence, with D α increasing as a function of α due to the approach to theintegrable point at α → ∞ . Exponential late time decay of the autocorrelation function.–
For finite system sizes L , approximating thediscrete Fourier sums by integrals eventually breaks down at very long times. In this regime the time evolution willbe dominated by the two smallest non-zero wave-numbers, k = ± π/L , leading to an exponential decay f j ( t ) ≈ L f ( k = 0 , t ) + (cid:88) k = ± π/L e ik j f ( k , t ) = 14 L (cid:104) πj/L ) e − λtc α (2 π/L ) α − (cid:105) (17)for the case of α ≤ .
5. The exponent of this exponential decay is hence expected to scale with the system size as ∼ (1 /L ) α − . This prediction is in agreement with our spin-2PI numerical results. Moreover, this result can be usedto extract the diffusion coefficient D α = λc α from finite size data. Corrections to scaling
For finite times, the two terms in Eq. 13 compete, leading to corrections to the leading-order scaling ansatz, Eq. 4.As we show below, these corrections survive until algebraically long times for α > .
5, and they add a logarithmiccorrection to the scaling ansatz at the threshold value α = 1 .
5, explaining all major deviations from scaling (4) weobserve in our numerical data.
Diffusive regime, α > . .– Taking into account both terms in Eq. 13, then rescaling as k → k √ D α t and j → y = j/ √ D α t where D α = λ/ (2 α − t → ∞ , wearrive at (cid:112) D α t f j ( t ) ≈ (cid:90) dk π exp (cid:0) iky − k (cid:1) (cid:16) − ( D α t ) − α +3 / c α (2 α − | k | α − (cid:17) = G ( y ) − ( D α t ) − α +3 / c α (2 α − α )16 π F (cid:20) α, , − y (cid:21) , (18)with F [ · , · , · ] denoting the Kummer confluent hypergeometric function. Most importantly, the latter exhibits heavytails ∼ y − α for y → ∞ , reproducing the finite time data for α = 2 in Fig. 3 of the main text.We show this behaviour more explicitly in Fig. 5, where we compare our numerical results to a scaling functioninvolving a single fit parameter D , (cid:112) D t f j ( t ) ≈ G ( y ) − √ D t F (cid:20) α, , − y (cid:21) , (19)following from Eq. 18 using lim α → sin( απ )Γ(1 − α ) = π/
12. Furthermore, according to Eq. 18 the approach toGaussian scaling is algebraically slow with exponent 1 . − α → α → .
5. As we show in the following, at thisspecial value α = 1 . Crossover point α = 1 . .– In the limit α → .
5, both prefactors in Eq. (13), c α and 1 / (3 − α ), diverge with theirdifference remaining finite, − c α + 1 / (3 − α ) → γ − / ≈ − .
92, with γ denoting the Euler-Mascheroni constant,resulting in a Gaussian leading order term, ∂ t f ( k ) ≈ − . λ k . However, an additional logarithmic correctionfrom lim α → . c α ( | k | α − − k ) = 0 . k log k also contributes. Following the derivation in Ref. [46], we rescale k as k → k (cid:112) λt Ω( λt ) /
2, with Ω( λt ) a function to be determined, and get (cid:112) λt Ω( λt ) / f j ( t ) ≈ (cid:90) dk π exp( ik ˜ y ) exp (cid:18) k Ω( λt ) ln (cid:18) γ − λt Ω( λt ) (cid:19) + k Ω( λt ) ln( k ) (cid:19) with scaling variable ˜ y = j/ (cid:112) λt Ω( λt ) /
2. The function Ω( λt ) is chosen in such a way that the first term in theexponent is equal to − k , reproducing the leading order Gaussian behavior [46]. This leads toΩ( λt ) = (cid:12)(cid:12)(cid:12)(cid:12) W − (cid:18) − γ − λt (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) , (20) . . . . . . . j/ √ Jt − − − − − − √ J t C ( j , t ) LOLO + NLO j/ p Jt Ω( λt ) / − − − − − p J t Ω ( λ t ) / C ( j , t ) α = 1 . α = 2 LOLO + NLO
FIG. 5.
Corrections to scaling.
The heavy tails found in the scaling functions for α ≥ . /λ ), with its numerical value being approximatelyequal in the fit to the (scaling) functions obtained from the LO and LO+NLO in a simultaneous k → t → ∞ expansion. with W − the secondary branch of the Lambert W-function. As discussed in [46], this gives Ω( λt ) ≈ ln( λt ) ∼ ln t , for t → ∞ , yielding a logarithmic correction to the scaling ansatz. Finally, we expand the resulting expression for largetime Ω( λt ) ∼ ln( λt ) (cid:29) (cid:112) λt Ω( λt ) / f j ( t ) ≈ (cid:90) dk π exp( ik ˜ y ) exp( − k ) (cid:18) k Ω( λt ) ln( k ) (cid:19) = G (˜ y ) (cid:18) λt ) (cid:18)(cid:0) − y (cid:1) ( − γ + ln(4)) + 2 exp (cid:0) ˜ y / (cid:1) F (1 , , (cid:20) , , − ˜ y (cid:21)(cid:19)(cid:19) , (21)with the superscript (1 , ,
0) denoting the derivative with respect to the first argument. This expression shows thelogarithmically slow convergence towards the Gaussian scaling function for α = 1 .
5, as well as a persistent logarithmiccorrection to the scaling form, with scaling variable ˜ y = j/ (cid:112) λt Ω( λt ) . Furthermore, for finite t , the above functionexhibits a heavy tail (cid:112) λt Ω( λt ) f j,α =1 . ( t ) ∼ y − and matches our 2PI results as shown in Fig. 5, using the singlefitting parameter λ . As times were note large enough in the simulations to be in the regime where Ω( λt ) ≈ ln( t ), weused the full expression in Eq. 20 for Ω( λt ) to fit the unrescaled f j ( t ) at a fixed time t . Superdiffusive regime, α < . .– While there are no qualitative corrections to the scaling function in this regime,the term ∼ k in Eq. 13 leads to a correction to the exponent of the hydrodynamic tail as α (cid:37) .
5. When evaluatingthe Fourier transform numerically with the full expression in Eq. (13) for α (cid:46) .
5, we still find an approximatehydrodynamic tail with a modified exponent reproducing the finite-time dTWA and spin-2PI results more closelythan the ’bare’ expression β α and hence accounting for the slight deviations between the numerical results and β α mentioned in the main text. For example, we get β α =1 . ≈ .
57 from this procedure, in agreement with 2PI( β α =1 . ≈ . ± .
02) and dTWA ( β dTWA α =1 . ≈ . ± . Classical master equation in dimension d > In the following, we extend the results of the classical Master equation to spins at locations r i in d dimensions. TheMaster equation (8) then reads ∂ t f r j = (cid:88) i (cid:54) = j W | r i − r j | ( f r i − f r j ) , with W | r i − r j | = λ | r i − r j | α . (22)Fourier transforming again diagonalizes the differential equation, yielding f ( | k | , t ) = f ( | k | ,
0) exp { ( W ( | k | ) − W ( | k | = ) t ) } . (23) Two spatial dimensions d=2.–
Denoting k ≡ | k | in the following, we evaluate the Fourier transform of thetransition amplitudes in continuous space with both an IR (system size L ) and UV (lattice spacing a = 1) cutoff,yielding W ( k ) − W ( k = 0) = λ (cid:90) L d r r (cid:90) π d θ (cid:16) e − ikr cos( θ ) − (cid:17) r α (24)= λ π (cid:90) L d r r − α ( J ( kr ) − , (25)with J ( kr ) denoting the zeroth order Bessel function of the first kind. For α < L → ∞ , hence we expect the dynamics to be described by the infinite ranged mean field modelin that regime. Concentrating on α ≥
1, we can remove the IR cutoff and arrive at W ( k ) − W ( k = 0) ≈ πλ (cid:90) ∞ d r r − α [ J ( kr ) −
1] (26) ≈ πλ (cid:20) k α − (cid:90) ∞ d x x − α ( J ( x ) −
1) + k (cid:90) r − α d r (cid:21) , (27) W ( k ) − W ( k = 0) ≈ λ (cid:20) − c α k α − + k π − α ) (cid:21) , (28)0with c α = 2 π (cid:90) ∞ d x x − α (1 − J ( x )) = − − α π Γ(1 − α )Γ( α ) for 1 < α < . (29)We see that a superdiffusive solution is obtained for 1 < α ≤
2, where the first term ∼ k α − dominates. Neglectingall other terms, we hence arrive at the scaling ansatz for a localized excitation f ( k, t = 0) = f r ( t ) = ( λc α t ) − α − F Dα (cid:32) | r | ( λc α t ) α − (cid:33) , (30)with F Dα ( y ) = 18 π (cid:90) ∞ d k k J ( ky ) e − k α − . (31) Three spatial dimensions d=3.–
We similarly get W ( k ) − W ( k = 0) = 2 πλ (cid:90) L d r r (cid:90) π d θ sin θ (cid:16) e − ikr cos( θ ) − (cid:17) r α (32)= 4 πλ (cid:90) L d r r − α (cid:18) sin ( kr ) kr − (cid:19) , (33)where we get an IR divergence and hence expect mean-field behavior for α < /
2. Considering only α ≥ /
2, we set L → ∞ and get for small k W ( k ) − W ( k = 0) ≈ λ (cid:20) − c α k α − + k π − α ) (cid:21) , (34)with c α = − π sin( πα )Γ(2 − α ) for 1 . < α < . . (35)We see immediately, that now superdiffusive behavior is seen for 1 . < α < . f r ( t ) = ( λc α t ) − α − F Dα (cid:32) | r | ( λc α t ) α − (cid:33) , (36)with F Dα ( y ) = 18 π y (cid:90) ∞ d k k sin( ky ) e − k α − ..