Introduction. High-temperature transport in quantum magnets is usually assumed to be incoherent and diffusive. In this context, the discovery of superdiffusive spin transport with dynamical exponent z = 3/2 in isotropic integrable quantum spin chains has been especially surprising [1][2][3][4][5][6][7]. While the full picture of superdiffusion in these systems is still coming into focus, a rapidly growing body of work has elucidated several of its features (see [8][9][10][11] for recent reviews). Using recent developments from the theory of generalized hydrodynamics (GHD) [12][13][14][15][16], spin superdiffusion has been argued to emerge from the dynamics of large, semiclassical solitons of Goldstone-like nature [4,17,18] that are generic in integrable models with continuous non-Abelian symmetries [19]. Numerical studies of these models show that the dynamical spin structure factor follows the Kardar-Parisi-Zhang [20] scaling form to high accuracy [7,[21][22][23][24][25]. However, describing higher-order fluctuations seems to require a more complicated fluctuating hydrodynamic theory [26][27][28][29].

Superdiffusive spin transport was also observed in recent inelastic neutron scattering [25], cold atom quantum microscopy [7], and superconducting qubit [29] experiments. Given that experiments are never described by perfectly integrable systems, it is natural to ask why superdiffusion appears to be so robust to integrability-breaking perturbations. In nearly integrable systems, the short-time dynamics are integrable, but at sufficiently long times, the dynamics become chaotic with diffusive transport properties [30][31][32][33]. The nature of this crossover from superdiffusive to diffusive spin transport in nearly integrable isotropic spin chains remains poorly understood. Perturbative arguments indicate that the symmetry of the integrability-breaking perturbation is crucial: for anisotropic perturbations (those that break spinrotation symmetry), the crossover timescales follow generic Fermi's golden rule (FGR) predictions and scale with the square of the perturbation strength [34]. In contrast, the timescales that characterize the crossover to diffusion associated with applying isotropic perturbations are known to be much slower [34][35][36][37], but these timescales have yet to be precisely characterized. This crossover is so slow that it remains controversial whether transport in the long-time limit is indeed diffusive [38][39][40]. The primary obstacles to numerically observing the crossover to diffusion for isotropic perturbations have been the computational expense of simulating quantum systems and the robustness of superdiffusion with respect to isotropic perturbations, which together make accessing the crossover timescale difficult [34,41]. Even for classical spin chains simulated using standard numerical techniques like Runge-Kutta methods, superdiffusive transport in models like the Ishimori spin chain [42] survives for all numerically accessible timescales in the presence of isotropic perturbations [35][36][37].

In this Letter we employ a recently developed classical integrability-preserving discrete-time integration scheme [22,43]. Thanks to the speedup afforded by this novel integration scheme, we are able to simulate times up to t max = 2 16 for systems of N = 200 000 spins, which is sufficient to observe the long timescale associated with the superdiffusiveto-diffusive crossover for isotropic perturbations. For both noisy and Floquet isotropic perturbations of strength λ, we observe a crossover to conventional diffusion t ∼ λ -α that is compatible with α = 6. For noisy chains this slow timescale can be understood in terms of a kinetic theory that characterizes the lifetimes of the large solitons responsible for spin transport.

Isotropic spin chains. For concreteness we focus on integrable spin chains with Hamiltonian H 0 , which are invariant under spin-rotation symmetry, perturbed by some spinrotation-invariant integrability-breaking term V of strength λ:

Here λ ∈ [0, 1] is a parameter interpolating between integrable (λ = 0) and purely nonintegrable (λ = 1) dynamics. We expect that our analysis will hold for any integrable spin chain, classical or quantum, that is invariant under some continuous non-Abelian symmetry. However, since diffusion emerges on very long timescales, we are constrained to simulating models for which one can reliably study very late times. Thus, first we consider classical spin chains: Each site n along the classical spin chain hosts an O(3) vector of unit norm S n = (S x n , S y n , S z n ) with canonical Poisson brackets {S i n , S j m } = i jk S k n δ nm , subject to nearest-neighbor ferromagnetic interactions. Second, to avoid the accumulation of truncation errors in schemes like Runge-Kutta, we will simulate a discretetime, Trotterized form of the dynamics. To this end, we will exploit the existence of Trotterizations of classical spin chains that preserve integrability [22,27,43]. The primary obstacle to previous efforts to study this superdiffusive-todiffusive crossover was the long timescales associated with the crossover, which are unreachable by both quantum and classical continuous-time integration schemes. Thanks to the recent development of these classical Trotterizations, we are now able to access the timescale relevant to the crossover.

In the remainder of this work, the integrable dynamics will be generated by an integrable Trotterization [22,27,43] of the classical Ishimori Hamiltonian H 0 =n ln(1 + S n • S n+1 ) [42]. The isotropic integrability-breaking perturbation will be a (again Trotterized) Heisenberg interaction V = n g n (t ) S n • S n+1 . We will consider two cases: (i) perturbations that are noisy, where g n (t ) will be taken to be random variables uncorrelated in both space and time, g n (t )g m (t ) = δ nm δ(tt ), where O denotes the average of O over noise, and (ii) perturbations that are time periodic, for which g n (t ) = 1.

For further details about the specific numerical implementation here, see Supplemental Material [44].

Spin transport. Numerical, theoretical, and experimental studies have shown that spin transport in the integrable limit λ = 0 is known to be superdiffusive, while energy transport is purely ballistic [1][2][3][4][5][6][7]25]. For simplicity, we will focus on the infinite-temperature regime where spins are initialized at random, although we expect for our conclusions carry over to arbitrary finite temperatures. In order to characterize spin transport, we will focus on fluctuations of spin transfer across a given bond in a background equilibrium state. We define spin transfer as the time-integrated spin current across the link between sites n and n + 1, Q(t ) = t 0 dt j n (t ), where j n (t ) is the spin current between sites n and n + 1, defined by the continuity equation ∂ t S z n + j n+1j n = 0. We invoke translation invariance and average this quantity over each bond n to improve statistics; we note that this form of translation invariance still holds for the case of the noisy perturbation due to the lack of any space or time dependence of the noise. Since there is no net spin transfer in a background equilibrium state, the thermal average Q(t ) = 0 vanishes. Fluctuations in spin transfer can be used to characterize spin transport through the relation

where z is the dynamical exponent, which characterizes transport by relating space and time scaling through t ∼ x z . For diffusive transport we see the dynamical exponent z = 2, whereas we see z = 3/2 in integrable isotropic chains, corresponding to superdiffusive (faster than diffusive) transport with an effective time-dependent diffusion constant D(t ) ∼ t 1/3 . In order to characterize the crossover from z = 3/2 to z = 2 upon applying an integrability-breaking perturbation, we will define the time-dependent dynamical exponent as the logarithmic derivative z(t

As a consistency check, we also characterize spin transport and the crossover to diffusion using the spin autocorrelation function C(t ) = S z n (t )S z n (0) ∼ t -1/z . We extract the dynamical exponent z from C(t ) and observe results compatible with those discussed below (see Supplemental Material [44]). Since transport becomes diffusive at long times with C(t ) ∼ 1/ √ D(λ)t, the autocorrelation function can also be used to compute the diffusion constant D(λ) from lim t→∞ ( √ t C(t )) [39]. Noisy perturbations: Numerics. The discrete-time model used here can be efficiently simulated up to very long times without accumulating Trotter errors, allowing us to reach maximum times of t = 2 16 for systems up to N = 200 000 spins, which is sufficient to observe the crossover to diffusion for noisy isotropic perturbations. The effective time-dependent dynamical exponent z(t ) extracted from the spin transfer shows a clear crossover from z = 3/2 to z = 2 at long times [Fig. 1(a)]. We find that the crossover occurs over a long timescale t ∼ λ -α with an exponent consistent with α = 6. We additionally extract the diffusion constant D(λ) from the autocorrelation function

) and find that it scales as D(λ) ∼ λ -2 (see Supplemental Material [44]). Since matching the diffusive and integrable behaviors at t = t gives D(λ) ∼ t 1/3 ∼ λ -α/3 for t ∼ λ -α , the scaling D(λ) ∼ λ -2 is also consistent with t ∼ λ -6 [44]. (In contrast, a process by which the solitons decay in an s-independent way through golden rule processes yields a crossover time t ∼ λ -2 and a diffusion constant D(λ) ∼ λ -2/3 , consistent with numerical data on anisotropic perturbations [34,44].) Kinetic theory. To understand this unusually long crossover timescale, we turn to a kinetic theory of the quasiparticles responsible for superdiffusion in the integrable case (for relevant recent reviews, see [10,16]). In integrable isotropic spin chains, superdiffusive transport with z = 3/2 is well established in terms of the system's quasiparticle excitations [4,5,18,19,45]. In classical chains, these quasiparticles are infinitely long-lived solitons which remain stable even at high temperatures. Solitons are characterized by their size s, a parameter that is quantized in quantum spin chains and continuous in classical chains. These solitons have velocity v s ∼ 1/s, so large solitons move slowly and occur with density ρ s ∼ 1/s 3 in thermal equilibrium [46]. They also carry spin m s = s in the vacuum, although this net magnetization can be screened in thermal states by a background consisting of other overlapping solitons [46][47][48]. Upon being screened, solitons cease to contribute to spin transport. Integrability results show that solitons are screened at a rate 0 s ∼ s -3 [4,5,49]. This means that at a given time t, small solitons with s s (t ) ∼ t 1/3 are fully screened and carry a net magnetization m s = 0, and so they do not contribute to spin transport. On the other hand, "giant" solitons s s (t ) ∼ t 1/3 remain charged and dominate transport. Combined with the slow velocity v s ∼ 1/s of giant solitons, this screening mechanism is responsible for the superdiffusion exponent z = 3/2 in integrable isotropic chains [4].

Soliton decay rates. In the presence of an integrabilitybreaking perturbation of strength λ = 0, solitons acquire a finite lifetime due to both vacuum decay processes and backscattering from collisions. The finite lifetime of a soliton of size s subjected to an integrability-breaking perturbation is characterized by the decay rate s,λ , which we expect to take the form s,λ ∼ n λ α n s -β n , with each term indexed by n corresponding to a different decay process. Note any process with β n 3 will occur slower than the integrable screening process 0 s ∼ s -3 for large solitons and is thus irrelevant. The crossover to diffusion can be understood as a competition between the integrable screening rate 0 s ∼ s -3 and the perturbation-induced decay rate s,λ . The crossover occurs when perturbative contributions to a soliton's decay rate s,λ overpower the screening rate 0 s , that is, the time at which the perturbative decay rate scales like s,λ ∼ s -3 .

The task of understanding the crossover therefore reduces to finding the leading-order term of the decay rate s,λ induced by a noisy isotropic perturbation. In general, analytically evaluating s,λ is very complicated, even perturbatively. We study this question numerically by considering the decay of a single soliton of size s in a vacuum background state, corresponding to an initial state (S x n , S y n ,

and φ n = arctan(tanh n 2s ) + n 2s [50]. This initial state corresponds to an exact right-moving soliton in the Ishimori chain with interaction ln(1 + S n • S n+1 )-left-movers can be obtained by shifting the phase φ n → φ n + π/2. In the context of our discrete-time numerics, this initial state is not technically an exact soliton (even for λ = 0), but it has a very strong overlap with the solitons of the discrete-time model.

In the integrable case (λ = 0), solitons propagate without decaying. Applying a noisy perturbation of strength λ = 0 causes a giant soliton to disintegrate into small solitons that diffuse through backscattering. To quantify this decay, we introduce the inverse participation ratio (IPR) G(t ) = n (S z n (t ) -1) 2 . Spin conservation fixes n (S z n -1) as a time-independent constant of order s; therefore, G(t ) provides a metric to quantify the localization of soliton. For a localized soliton, G(t ) remains a constant of order s, while for a delocalized soliton it decays with system size as 1/L. For large systems and long times, we observe an exponential decay G(t ) ∼ e -s,λ t , which allows us to extract the decay rates s,λ numerically for various soliton sizes s and perturbation strengths λ [Fig. 1(b)]. We find that

in agreement with recent perturbative arguments [Fig. 1(c)] [34]. The long lifetimes ∼s 2 of large solitons can be understood in terms of the Goldstone-like nature of the solitons (3): matrix elements of isotropic perturbations are suppressed as 1/s acting on a soliton of size s [34]. Comparing the decay rates (4) to the integrable screening rates, we find the crossover soliton size s ∼ λ -2 and the crossover timescale

in agreement with the anomalously long timescale observed in our transport data [Fig. 1(a)] and extracted diffusion constant (see Supplemental Material [44]). The argument presented above only relies on the general form of the vacuum decay rate of solitons and on the relative perturbation strength, both of which are generic quantities; therefore, this argument for the crossover timescale is independent of any technicalities of the numerical integration scheme used here and we expect it to hold generically.

While the vacuum decay rates (4) fully explain the crossover timescale t ∼ λ -6 , they also lead to logarithmic corrections to diffusion, since the giant solitons' lifetime ∼s 2 in the vacuum is long enough to lead to anomalous behavior [34]. Whether interactions can induce decay rates ∼s -α with α < 2 that restore diffusion remains an important question.

Floquet perturbations. The previous section considered noisy perturbations V (t ) =n g n (t ) S n • S n+1 , with g n (t ) uncorrelated random variables in both space and time. We now consider non-noisy (Floquet) perturbations with g n (t ) = 1. The crossover associated with a noisy perturbation is primarily driven by the vacuum decay rate s,λ ∼ λ 2 /s 2 ; however, in the absence of noise, isolated solitons remain stable even in the nonintegrable classical Heisenberg chain [51]. The largest contribution to the decay rate must therefore arise from many-body soliton scattering processes; this term should also exist in the decay rate for the noisy perturbations but is masked by the leading-order contribution from the vacuum decay process. Therefore, we expect the crossover timescale associated with Floquet perturbations to be at least as slow as the noisy crossover timescale t ∼ λ -6 , if not slower. This picture is consistent with the recent numerical works on classical spin chains [35][36][37], which were unable to observe a systematic crossover towards diffusion on timescales accessible by standard continuous-time integration schemes like adaptive Runge-Kutta methods.

The discrete-time numerics used here are able to access very long times without finite time-step errors, which allows us to observe timescales inaccessible to Runge-Kutta setups. Even so, we observe a curve collapse for the dynamical exponent z(t ) that is slower than the noisy case but still consistent with t ∼ λ -α scaling with α ≈ 6 [Fig. 2(a)]. For consistency, we extract the diffusion constant D(λ) from the autocorrelation function using the same method as described for the noisy perturbation, which is also compatible with t ∼ λ -6 scaling [Fig. 2(b)]. We note that our data cannot definitively exclude a different exponent in the true scaling limit t → ∞ and λ → 0, and different values of α lead to acceptable collapses as well (see Supplemental Material [44]).

As previously discussed, we expect that the crossover to diffusion is driven by collisions between multiple solitons; however, since the IPR is not a meaningful metric for a finite density of solitons, it is difficult to determine collision decay rates using the same methods as were used to extract the noisy vacuum decay rate. By simulating experiments in which large solitons first interact with a nontrivial background under λ = 0 dynamics, and then filtering out the background by setting λ = 0 and allowing it to "expand" into a vacuum region, we are able to show that the soliton lifetime increases with s, although the data are too noisy to extract a reliable exponent [44]. We note that the crossover scale t ∼ λ -6 is compatible with higher-order perturbative decay rates s,λ ∼ λ 2(n+1) s n-2 for any n. On general grounds we expect the collisions between s-sized solitons at large s and O(1)-sized solitons to scale as 1/s 2 , since a small soliton experiences a large one as a local vacuum rotation. In effect, a thermal background of small solitons acts as noise on the large ones, relating the Floquet and noisy cases. This process would once again lead to logarithmically corrected diffusion. Any higher-order processes with n > 0, if present, would lead to a full crossover to diffusion. Determining the form of the decay rates that emerge from interactions remains an important challenge for future work.

Discussion. We have numerically observed the slow crossover from superdiffusive to diffusive transport for isotropically perturbed integrable classical spin chains, with a crossover timescale t ∼ λ -α with α ≈ 6. Our ability to reach this anomalously long crossover timescale is due to the classical discrete-time algorithm used in this work, in contrast to previous work that studied this question with matrix-product states or classical continuous-time numerical methods. The t ∼ λ -6 scaling holds for both noisy and Floquet perturbations. While the exponent α = 6 is expected perturbatively in the noisy case and has a heuristic justification in the clean case, a full explanation of the crossover to diffusion in either case requires analyzing the decay rates that arise from n-body soliton scattering processes and presents an important challenge for future theoretical work. For example, it is possible that the crossover to diffusion happens in multiple stages, with an intermediate regime of logarithmically corrected diffusion that gets parametrically large in λ. We note that although our numerics are consistent with an exponent α ≈ 6 for Floquet perturbations, we cannot conclusively exclude other values of α in the limit of λ → 0 and hope that future theoretical and numerical advances will pinpoint the exact exponent.

Note added. Recently we became aware of a related work by McRoberts and Moessner [52] reporting a crossover timescale of t ∼ λ -3 in an energy-conserving model and analyzing its temperature dependence. Future work would be needed to understand why the crossovers in both models appear to be dominated by different decay processes on accessible timescales.