Decay of spin-spin correlations in disordered quantum and classical spin chains

The real-time dynamics of equal-site correlation functions is studied for one-dimensional spin models with quenched disorder. Focusing on infinite temperature, we present a comparison between the dynamics of models with different quantum numbers $S = 1/2, 1, 3/2$, as well as of chains consisting of classical spins. Based on this comparison as well as by analyzing the statistics of energy-level spacings, we show that the putative many-body localization transition is shifted to considerably stronger values of disorder for increasing $S$. In this context, we introduce an effective disorder strength $W_\text{eff}$, which provides a mapping between the dynamics for different spin quantum numbers. For small $W_\text{eff}$, we show that the real-time correlations become essentially independent of $S$, and are moreover very well captured by the dynamics of classical spins. Especially for $S = 3/2$, the agreement between quantum and classical dynamics is remarkably observed even for very strong values of disorder. This behavior also reflects itself in the corresponding spectral functions, which are obtained via a Fourier transform from the time to the frequency domain. As an aside, we also comment on the self-averaging properties of the correlation function at weak and strong disorder. Our work sheds light on the correspondence between quantum and classical dynamics at high temperatures and extends our understanding of the dynamics in disordered spin chains beyond the well-studied case of $S=1/2$.


I. INTRODUCTION
Noninteracting particles in one and two spatial dimensions are localized even for arbitrarily small values of disorder [1,2]. Generalizing this well-understood concept of Anderson localization to the presence of interactions has been a major objective of modern condensed matter physics. Based on pioneering early works [3,4], and a large number of subsequent (mostly numerical) studies (see e.g. [5][6][7][8][9][10][11]), many-body localization (MBL) is nowadays believed to be a generic property in one-dimensional short-range lattice models with sufficiently strong randomness [12,13]. (For another point of view, see [14].) The phenomenology of the localized phase is best understood in terms of an emergent set of local integrals of motion [15,16], which, e.g., explain the slow logarithmic growth of entanglement in time [17,18]. Moreover, localized systems violate the eigenstate thermalization hypothesis [12,19,20], which reflects itself in the area-law entanglement of eigenstates [21], as well as the failure of MBL systems to reach thermal equilibrium at long times.
Progress in the field of many-body localization has also been fostered by the advance of novel experimental platforms [22,23], which can be very well isolated from the environment. In particular, such experiments also allow to tackle new and exciting questions about the existence of MBL, e.g., in dimensions larger than one [24], in the presence of long-range interactions [25], or in systems which are weakly coupled to a thermal bath [26].
More recently, there has also been much effort to put the phenomenon of MBL into a broader perspective and understand the dynamics of models, which might not be strictly localized but nevertheless exhibit anomalously small relaxation rates and very long equilibration times [27]. Examples include, e.g., disordered models with power-law interactions [28], systems with two particle species and large mass ratio [29][30][31][32], as well as Hubbard models where the disorder only acts on the charge degrees of freedom but not on the spin [33,34].
In this context, the present work scrutinizes the dynamics and the fate of MBL for a class of disordered models which have received less attention so far. In particular, while a plethora of works has explored manybody localization in disordered spin-1/2 systems, much less is known about the effect of disorder for larger spin quantum numbers S > 1/2 [35]. Since genuine MBL is believed to be a pure quantum phenomenon, it is an intriguing and open question to what extent localization can occur upon increasing S, where the quantum spins are supposed to become more and more akin to their classical counterparts. In this work, we shed light onto this question by comparing the infinite-temperature dynamics of equal-site correlation functions for disordered classical and quantum spin chains with S = 1/2, 1, 3/2. Based on this comparison as well as by analyzing the statistics of energy-level spacings, we show that the putative many-body localization transition is shifted to considerably stronger values of disorder for increasing S. Furthermore, introducing an effective disorder W eff , we find a mapping between the dynamics for different quantum numbers S. For small W eff , we show that the realtime correlations become essentially independent of S, and are moreover very well captured by the dynamics of classical spins. Especially for S = 3/2, the agreement between quantum and classical dynamics is remarkably observed even for very strong values of disorder. This paper is structured as follows. In Secs. II A, II B, and II C, we introduce the models and observables which are studied in this work, and explain our numerical approach. In Sec. III, we then present our numerical results for the dynamics in clean and disordered quantum and classical spin chains. Finally, we summarize and discuss our findings in Sec. IV.

A. Model
We consider the one-dimensional Heisenberg model with quenched disorder and periodic boundary conditions, described by the Hamiltonian where L denotes the number of lattice sites, J = 1 is the antiferromagnetic exchange constant, and the on-site magnetic fields h l are randomly drawn from a uniform distribution h l ∈ [−W, W ], with W setting the magnitude of disorder. Moreover, the S l = (S x l , S y l , S z l ) are either quantum spin-S operators with S = 1/2, 1, 3/2 or, in the classical case, real three-component vectors of unit length. Note that H conserves the total magnetization, i.e., [ l S z l , H] = 0, for all values of W and S. Note further that for the particular case of S = 1/2, this model is equivalent via a Jordan-Wigner transformation to interacting spinless fermions hopping in a random potential. Moreover, for S = 1/2 and W = 0, H is integrable in terms of the Bethe ansatz.
It is instructive to briefly consider two limiting cases of the model (1), namely (i) the most quantum case S = 1/2 and (ii) the case of classical spins. On the one hand, for S = 1/2, the Hamiltonian (1) has become a canonical model to study the phenomenon of many-body localization. Specifically, the spin-1/2 model is believed to undergo a transition from an ergodic phase into a manybody localized regime above a critical disorder strength W c ≈ 3.5 (see e.g. [6,9]), although also larger values have been suggested [36,37]. On the other hand, also in the case of classical spins, strong disorder has been shown to cause a drastic reduction of transport coefficients and anomalously slow relaxation [38,39]. The occurrence of genuine MBL, however, is not expected for classical spin models. This can be understood for instance from a rareregion argument, where small parts of the chain are only weakly disordered and exhibit chaotic dynamics, which eventually causes thermalization of the full system at sufficiently long time scales [27].
Not least due to the higher computational requirements, much less is known about the dynamics and the putative MBL transition for the Hamiltonian (1) between these two limiting cases, i.e., for larger spin quantum numbers S > 1/2. The present work attempts to bridge this gap and sheds light on the dynamics in disordered spin chains with S = 1, 3/2.

B. Equal-site correlation function
As a central quantity of interest in this paper, we study the dynamics of the disorder-averaged (denoted by the overbar) equal-site correlation function for the disordered spin model (1) at formally infinite temperature. Note that due to periodic boundary conditions, the specific site l to calculate C(t) is arbitrary (upon averaging the bare C(t) over sufficiently many disorder realizations, cf. Sec. II C 3 below).
In the quantum case, C(t) is given by where d = (2S + 1) L is the dimension of the Hilbert space, and S z l (t) = e iHt S z l e −iHt is the time-evolved operator in the Heisenberg picture. Moreover, in the case of classical mechanics, · denotes the ensemble average over classical trajectories, where for every trajectory all the S l are randomly initialized at t = 0 (within the constraint of |S l | = 1), and evolve in time according to the Hamiltonian equations of motioṅ In addition to C(t), we also study the corresponding spectral function C(ω), which follows from a Fourier transform of the real-time data according to and which is relevant to, e.g., the relaxation rate probed in nuclear magnetic resonance experiments [40]. Note that, since the real-time dynamics of C(t) can in practice only be evaluated up to a finite cutoff time t max < ∞, the Fourier transform (5) yields the spectral function C(ω) with a finite frequency resolution δω = π/t max .

Decay of the correlation function
The autocorrelation function C(t) in Eq.
(2) can be interpreted as the "return probability" of a single spin excitation which is created at some lattice site l, and spreads within the system with respect to an infinitetemperature many-body background.
On the one hand, in the ergodic phase, the dynamics of C(t) can be understood as follows. After a quick initial decay starting from C(t) is eventually expected to display some slower hydrodynamic behavior (since the total S z is conserved), which reflects itself in a power-law decay For one-dimensional systems, α = 1/2 would here correspond to conventional diffusive transport. In this context, let us note that for the spin-1/2 version of Eq. (1), there is numerical evidence for the existence of a subdiffusive phase in the regime of low to intermediate disorder below the localization transition [41][42][43][44]. This subdiffusive regime manifests itself, e.g., in an anomalously slow relaxation of the autocorrelation function C(t) with exponent α < 0.5, as well as sublinear buildup of entanglement [45][46][47][48][49]. Eventually, for even longer times, the spin excitation has explored the entire system (assuming a finite chain length L), and C(t) will saturate to a constant long-time value which scales as ∝ 1/L [50].
On the other hand, in the localized regime, the decay of the autocorrelation function is significantly slower. In particular, C(t) saturates to a nonzero value at long times (or to a value > 1/L in a finite system). In this sense, the long-time value of C(t) can be used as an order parameter for the MBL transition [44].
C. Numerical approach

Quantum dynamics
The numerical evaluation of the quantum expectation value (3) is a challenging task, not least due to the exponential growth of the Hilbert space and the necessity to study long time scales. We here tackle this numerical challenge by using an efficient pure-state approach based on the concept of dynamical quantum typicality (DQT) [51][52][53][54][55][56][57][58][59][60][61]. Within this concept, the equal-site correlation function C(t) in Eq. (3) can be written as the expectation value within a single pure state [35,62,63], where |ψ(0) is prepared according to and the reference pure state |ϕ = d k=1 c k |k is drawn at random from the full Hilbert space according to the unitary invariant Haar measure [56]. In practice, the states |k here denote the product basis of the local S z projections (Ising basis), and the complex coefficients c k are randomly drawn from a Gaussian distribution with zero mean and unit variance. As an aside, let us note that the application of the square root in Eq. (9) becomes straightforward since S z l + S is nonnegative and diagonal in our working basis. Importantly, the statistical error ǫ = ǫ(|ϕ ) in Eq. (8) scales as ǫ ∝ 1/ √ d [51,57], and becomes very small even for moderate system sizes L (especially for S > 1/2). (We demonstrate the smallness of statistical errors further below in Sec. III B 4.) Moreover, thanks to the sparse matrix structure of generic few-body operators, the time evolution of the pure state |ψ(t) can be efficiently generated by iteratively solving the real-time Schrödinger equation, with a small discrete time step δt by means of, e.g., Runge-Kutta (RK) schemes [58,60], Krylov subspace techniques [64], Trotter product formulas [65], or Chebyshev expansions [66,67]. Due to this fact, we can treat system sizes beyond the range of standard exact diagonalization (ED) (here up to L = 14 for S = 3/2, which is equivalent to L = 28 for S = 1/2), and study comparatively long time scales.
Eventually, let us emphasize that DQT is independent of the validity of the ETH and solely relies on the largeness of the Hilbert space. Thus, DQT also allows for accurate calculations in strongly disordered models which undergo a MBL transition [68].
Independent of the typicality relation in Eq. (8), each side of this relation can be interpreted as a certain type of imperfect "echo protocol" (see also [69]). To make this a bit more transparent, the r.h.s. of Eq. (8) can be slightly rewritten as i.e., after evolving the out-of-equilibrium initial state |ψ(0) for some time t, a "perturbation" is applied in form of the operator S z l . Subsequently, the perturbed state is evolved backwards in time and the overlap between the resulting state e iHt S z l e −iHt |ψ(0) and the initial state |ψ(0) is measured.

Classical dynamics
In contrast to the quantum case, the memory requirements for the classical spin chain do not scale exponentially, but only linearly in system size. Therefore, the Hamiltonian equations of motion (4) can be solved for huge systems and long times with significantly less computational resources. Specifically, we here employ a fourth-order RK scheme, where δt = 0.01 is chosen short enough to conserve the total energy and total magnetization to very high accuracy [70,71]. For a proper comparison between quantum and classical dynamics, however, we here present data mostly for system sizes L, which can also be treated quantum mechanically.

Averaging
As it is standard in context of disordered systems, the bare correlations C(t) have to be averaged over N h independent realizations of the random magnetic fields h l in order to obtain reliable results, Accordingly, the results for C(t) presented in this paper have to be understood as the average from N h = 100 − 1000 random realizations, depending on the system size L and the spin quantum number S. (Note that in the case of classical spins, the ensemble average C(t) for a given disorder realization is obtained from approximately 10 3 independent trajectories.) Moreover, in Sec. III B 4, we are concerned with the question whether or not C(t) is a self-averaging quantity. Here, self-averaging refers to the fact that for larger and larger systems, fewer and fewer disorder realizations are needed to correctly represent the whole statistical ensemble. To this end, we introduce the relative variance R(t) of sample-to-sample fluctuations [72], We call C(t) self-averaging if R(t) decreases upon increasing the system size, i.e., R(t) ∝ L −ν with 0 < ν ≤ 1, where ν = 1 would correspond to strong self-averaging.
As another probe of self-averaging, we also consider the log-averaged correlation function C log (t) [73], defined as Naturally, the corresponding log-averaged spectral function C log (ω) follows from Eq. (5) upon replacing C(t) by C log (t). It is easy to see that the log-average C log (t) is in fact equivalent to the geometric mean (whereas C(t) is the "standard" arithmetic mean), which is also the reason why C log (t) is often referred to as typical correlation [41,73].

III. RESULTS
Let us now present our numerical results. In Sec. III A, we begin by studying the clean Heisenberg chain (W = 0) for S = 1/2, 1, 3/2 and classical spins, which serves as a convenient starting point for the actual discussion of disordered models in Sec. III B. Note that, in order to simplify the notation, we drop the overbar for the averaged quantities in the following. Thus, whenever we write C(t) or C log (t), we implicitly refer to C(t) or C log (t).

A. Clean model
In Fig. 1 (a), we show the equal-site correlation function C(t) at a fixed system size L = 14, both for classical 0.01 spins as well as quantum spins with S = 1/2, 1, 3/2 at vanishing disorder W = 0. In order to account for the different quantum numbers S, the time axis is renormalized according to [74] t → t S , with S = S(S + 1) .
(Note that S = 1 in the case of classical spins.) Moreover, the curves are normalized by their initial value C(0) = S 2 /3, cf. Eq. (6). In particular, we find that this rescaling yields a convincing data collapse of C(t), and we can identify an intermediate time window 2 t S 10, where C(t) decays as a power law, C(t) ∝ t −α , with exponent α ≈ 2/3 for all curves shown here. On the one hand, for the integrable S = 1/2 chain, this is consistent with other studies [75][76][77][78], which report that spin transport in the isotropic Heisenberg chain is described by the Kardar-Parisi-Zhang universality class [79,80]. On the other hand, for the nonintegrable S = 1, 3/2 models, the situation is less settled. While it has been put forward that superdiffusion might persist despite the absence of integrability [35,81], it has also been argued that α will eventually approach its diffusive value α = 1/2 asymptotically for larger system sizes and longer time scales [82]. In this context, let us note that even for the classical spin chain, the unambiguous detection of diffusion is known to be a delicate problem [83][84][85][86][87].
In Fig. 1 (a), we moreover find that the spin-1/2 curve exhibits some additional oscillations at intermediate times, whereas the curves for larger S = 1, 3/2 are almost indistinguishable from the dynamics of classical spins. We trace these oscillations back to (i) enhanced quantum fluctuations at S = 1/2, and (ii) specific prop- erties of the model under consideration. For additional numerical results of C(t) in other disorder-free models, we refer to Appendix A. Interestingly, breaking the integrability of the S = 1/2 model, e.g., by means of a next-nearest neighbor interaction, not necessarily leads to a further improvement of the data collapse.
In addition to the real-time data, the corresponding spectral functions C(ω) are shown in Fig. 1 (b). Also in this case, we find that a proper rescaling of the energy scale, ω → ω/ S, results in a convincing data collapse. In fact, notable deviations can only be observed at large frequencies ω/ S 2.5, where C(ω) decays exponentially [88], and this decay turns out to be slightly faster for S = 1/2.
The data presented in Fig. 1 already exemplify one important result of the present paper. Namely, there exist models and observables where quantum dynamics (i) becomes almost independent of the specific quantum number S and (ii) can be well captured by the dynamics of classical spins. Admittedly, this fact might not be entirely surprising since we are dealing with dynamics at formally infinite temperature where quantum effects are less pronounced. Nevertheless, it is still a remarkable result that the dynamics of classical spins not only qualitatively, but even quantitatively reproduces quantum dynamics for S = 1, 3/2. As shown below, this agreement between quantum and classical dynamics becomes even more interesting when studying the dynamics in models with finite disorder W > 0. B. Disordered model

Level-spacing distribution
Before studying dynamics for W > 0, a first orientation for the putative MBL transition in the model (1) with S ≥ 1/2 can be obtained from the mean ratio r of adjacent level spacings, where δ n = |E n+1 −E n | with E n being the energy-ordered eigenvalues of H, and the brackets denote both, averaging over approximately 1/3 of the levels in the center of the spectrum, and over sufficiently many independent realizations of the random h l . For an ergodic system, the level statistics is supposed to follow a Wigner-Dyson distribution with r ≈ 0.53, whereas the occurrence of MBL at strong disorder reflects itself by the onset of Poissonian level statistics with r ≈ 0.39 [6]. Our numerical results for r are summarized in Fig.  2, where we restrict ourselves to a fixed system size L = 8 and eigenstates in the zero-magnetization subsector. (Note that ED becomes very costly for S = 3/2 and L > 8.) As can be seen in Fig. 2 (a), the ratio r interpolates between r = 0.53 and r = 0.39 for all values of S = 1/2, 1, 3/2 upon increasing the disorder strength W . (The small dip of r for S = 1/2 and W = 0.5 can be explained due to the vicinity of the integrable point at W = 0.) Importantly, however, we find that the transition of r from Wigner-Dyson towards Poissonian level statistics occurs at larger and larger values of W if the spin quantum number S is increased. This result suggests that a transition to a many-body localized regime in spin chains with S > 1/2, if existent at all, probably requires a stronger and stronger critical disorder W c (compared to W c ≈ 3.5 for S = 1/2). This is an important result of the present paper.
In view of the successful rescaling found in the context of Fig. 1, let us now introduce an effective disorder strength W eff according to, Plotting r versus W eff in Fig. 2 (b), we in fact observe a reasonable data collapse for small W eff ≈ 1 and large W eff ≈ 8. In the intermediate regime W eff ≈ 4, however, this collapse clearly breaks down and r is still generally larger for larger S. Given the facts that (i) r naturally depends on the system size L [6] (see also Appendix B  for additional data), and (ii) the proper analysis of finitesize data around the critical disorder is known to be very delicate [89], it is difficult to estimate whether or not this behavior persists in the thermodynamic limit L → ∞. Nonetheless, the effective disorder W eff in Eq. (17) actually turns out to be very useful when comparing C(t) and C(ω) for different quantum numbers S below.

Dynamics
Let us now proceed to the actual discussion of C(t) and C(ω) in disordered spin chains. In particular, it is instructive to compare the dynamics for different quantum numbers S, while using the same value of the effective disorder W eff for each S. Since the MBL transition in the S = 1/2 chain is well-studied, this model serves as a point of reference and we set the values of the bare disorder W for the S = 1/2 chain to W = 0.5, 1, 2 and W = 5. In view of the putative critical disorder W c ≈ 3.5, we thus study the dynamics of C(t) and C(ω) both in the ergodic and in the localized regime of the spin-1/2 model. Moreover, given these values of W for S = 1/2, we consequently obtain effective disorder values W eff = W/ 3/4 ≈ 0.58, 1.15, 2.31, 5.77. The corresponding bare disorder values W for S = 1, 3/2 then directly follow from W = W eff S(S + 1) and are summarized in Table I. Let us start our numerical analysis with the case of small disorder. To this end, Fig. 3 (a) shows the equalsite correlation function C(t) for W eff ≈ 0. 58 Fig. 3, but now for stronger effective disorder W eff ≈ 2.31, 5.77. Fig. 1 (a), we observe that the curves for S = 1, 3/2 and for classical spins are essentially indistinguishable from each other. Moreover, while the long-time value of the S = 1/2 curve appears to be slightly higher compared to the curves of the other S, we show later in Sec. III B 4 that this finding is just a finite-size effect and becomes less pronounced for increasing L. Remarkably, as shown in Fig. 3 (b), this mapping between classical and quantum dynamics persists for S = 1, 3/2 also at the slightly stronger disorder W eff ≈ 1. 15. In contrast to the previous case in Fig. 3 (a), however, we now find that the S = 1/2 model clearly deviates from the curves for S ≥ 1 and exhibits distinctly slower dynamics, which is consistent with subdiffusive transport in this parameter regime. The corresponding spectral functions C(ω) for W eff ≈ 0.58, 1.15 are presented in Figs. 3 (c) and (d). Also for C(ω), the rescaling ω → ω/ S and W → W eff leads to a convincing data collapse of all curves shown here. Interestingly, the clear deviations of the real-time S = 1/2 data observed in Fig. 3 (b) are much less pronounced in the frequency representation [see Fig. 3 (d)]. Moreover, we generally find that the overall shape of C(ω) is rather similar compared to the disorder-free case in Fig. 1 (b), except for the regime of very small frequencies ω → 0, where C(ω) starts to diverge upon increasing W .
To proceed, numerical results for the real-time correlation C(t) at stronger disorder W eff = 2.31 and W eff = 5.77 are shown in Figs. 4 (a) and (b). In the case of W eff = 2.31, C(t) now decays distinctly slower for all values of S, which might indicate that subdiffusion also occurs for larger spin quantum numbers S = 1, 3/2. Moreover, compared to Fig. 3 (b), we observe that not only the S = 1/2 data but also the S = 1 curve starts to deviate from the results for S = 3/2 and classical spins. Furthermore, in the case of W eff = 5.77 [ Fig. 4 (b)], we find that C(t) becomes practically time-independent for S = 1/2, which signals the onset of many-body localiza-tion. In contrast, for S = 1, 3/2 as well as classical spins, C(t) clearly has a nonzero slope and continues to decay at long times (albeit this decay is very slow for S = 1 compared to S = 3/2 and classical dynamics).
Eventually, Figs. 4 (c) and (d) show the respective spectral functions C(ω) at strong disorder. Compared to the cases of vanishing or small disorder in Figs. 1 and 3, C(ω) now behaves qualitatively different and develops a long tail at large ω. Nevertheless, at least for W eff = 2.31, we find a reasonable data collapse of C(ω) for all values of S. In Fig. 4 (d), we moreover show that the onset of MBL in the spin-1/2 model reflects itself in a concave shape of C(ω) with a distinct peak at ω/ S ≈ 1, consistent with earlier results in Ref. [73] (see also [90]). While such a pronounced feature is absent for S = 3/2 and classical spins, a small double-peak structure around ω/ S = 1 emerges in the case of S = 1.

Intermediate conclusion
Given our numerical results in Figs. 3 and 4, we conclude that for small values of disorder, the dynamics of C(t) and C(ω) becomes almost independent of the quantum number S when simulated at the same effective disorder strength W eff . In contrast, for stronger values of disorder, this mapping at least partially breaks down. Specifically, while the spin-1/2 model appears to be localized at W eff = 5.77, this does not seem to be the case for S = 1, 3/2 as well as classical spins. In particular, we have found a very good agreement between the dynamics for S = 3/2 and classical spins, at least for the values of W eff and the time scales considered in Figs. 3 and 4.
Generally, the dynamics of C(t) and C(ω) in Figs. 3 and 4 is consistent with the behavior of the level-spacing distribution r discussed in the context of Fig. 2 (b). Namely, for a fixed value of W eff , r is still closer to the chaotic Wigner-Dyson value for a larger value of S. From the combination of Figs. 2-4, one might speculate that MBL eventually occurs also in models with S = 1 and S = 3/2, but the critical disorder strength has to be even stronger than the largest W eff = 5.77 considered by us. (Note that for S = 3/2, this W eff already corresponds to the very large bare disorder W ≈ 11.18, cf. Table I). In view of the small system sizes numerically accessible for S = 1, 3/2, however, we here refrain from a more detailed analysis of the putative onset of MBL in these models.

Finite-size analysis and self-averaging
While we have shown numerical result for C(t) and C(ω) in Figs. 3 and 4 for a fixed system size L = 12, let us briefly comment also on the finite-size scaling of these quantities. Since the numerical calculation of larger system sizes for S = 1, 3/2 becomes unfeasible very quickly, we here restrict ourselves to the cases of S = 1/2 and classical spins.
To begin with, Fig. 5 (a) shows the real-time correlation function C(t) for spin-1/2 chains with different system sizes L = 10, 12, . . . , 20 and two different values of the effective disorder W eff ≈ 0.58, 5.77. As a comparison, we also depict corresponding data for classical chains with L = 20. On the one hand, for strong W eff = 5.77, we observe that the spin-1/2 curves are essentially independent of L for the time scales and systems sizes shown here. On the other hand, for W eff ≈ 0.58, one finds that C(t) remains converged up to times t S 10, while finite-size effects occur at later times and the long-time value scales as C(t → ∞) ∝ 1/L, as expected for the ergodic regime. Moreover, comparing the L = 20 curves for S = 1/2 and classical spins at this W eff , we find that their long-time behaviors agree very well with each other. Connecting to our earlier discussion in the context of Fig.  3 (a), this observation confirms that the effective disorder W eff , at least for small disorder, provides a useful mapping between classical and quantum dynamics with different S.
Since it is possible to treat classical spin chains with significantly larger system sizes, Fig. 5 (b) exemplarily shows numerical results for L = 20 and L = 200 up to rather long time scales t ≤ 1000. (Note that this is still far from the maximum L and t values accessible [38,70]). For weak W eff ≈ 1.15 and large L = 200, we observe a pronounced diffusive decay, C(t) ∝ t −1/2 , which persists essentially over the entire time window depicted. Furthermore, for stronger W eff ≈ 5.77, 8.00, one finds that although the dynamics of C(t) is very slow for t 100, the clearly visible decay at longer times t 500 is incompatible with localization. Comparing data for L = 20 and L = 200, we moreover find that finite-size effects remain irrelevant at these large disorder values, even for the long time scales studied in Fig. 5 (b). Our results indicate that while high-temperature spin transport in classical spin chains becomes strongly suppressed upon increasing disorder, a genuine MBL phase is absent in the classical model. These findings for the equal-site correlation function are consistent with earlier results from Ref. [38], which has focused on energy transport instead.
Next, let us discuss the self-averaging properties of C(t). To this end, the relative variance R(t) of sampleto-sample fluctuations [cf. Eq. (13)] is shown in Fig. 6 for spin-1/2 chains with vanishing disorder W = 0, as well as weak and strong disorder W = 0.5, 5. In particular, we restrict ourselves to short time scales t ≤ 10, where C(t) is free of trivial finite-size effects, cf. Fig. 5 (a).
On the one hand, for W = 0, R(t) should in principle be strictly zero since there is no disorder. Note, however, that we calculate C(t) by means of a typicality approach which relies on randomly drawn pure states and comprises a finite statistical error, cf. Eq. (8). The nonzero data for R(t) shown in Fig. 6 (a) can therefore be interpreted as the accuracy of this typicality approximation. Consistent with our discussion in Sec. II C 1, we find that R(t) decreases exponentially upon increasing L for all times shown here. This demonstrates that typicality indeed provides an accurate numerical method to determine C(t) and, in the spirit of Eq. (13), C(t) is "super" self-averaging at W = 0. On the other hand, for finite disorder W > 0, the behavior of R(t) is more complicated [see Figs. 6 (b) and (c)]. For very short times, we find that the disorder apparently has no impact, such that R(t) is dominated by the typicality contribution and decreases exponentially with L. In contrast, for longer times, this exponential scaling clearly breaks down and the curves for different L are very similar to each other. While for W = 0.5 one might still argue that R(t) slowly decreases with increasing L [see also inset in Fig. 6 (b) with linear axis], R(t) appears to be independent of L at strong W = 5 (except for residual statistical errors in our numerics). Consistent with recent results from Ref. [72], this indicates that the equal-site correlation function C(t) is self-averaging for short times t → 0. At longer times, self-averaging is either much weaker (at small disorder), or might break down (at stronger disorder). Note that, since our conclusions are based on finite-size data with L ≤ 20, we cannot rule out that a power-law scaling R(t) ∝ L −ν eventually emerges for larger system sizes (especially if ν is small).
To proceed, we discuss finite-size effects for the spectral function C(ω). In Figs. 7 (a)-(c), C(ω) is shown for three different disorder values W = 0.5, 2, 5 and three different system sizes L = 12, 16, 20. Generally, one observes that C(ω) is essentially free of finite-size effects and the curves for different L coincide very well with each other. The only exception to this finding is the case of weak disorder W = 0.5 and very small ω → 0 [ Fig. 7 (a)], where the curves do not collapse anymore (see also Ref. [73] for similar results and further discussion).
Eventually, we also study the log-averaged spectral function C log (ω). In this context, let us note that Ref. [73] investigated a very similar quantity and reported on a breakdown of self-averaging at large disorder. Specifically, Ref. [73] found that the logarithmic average became strongly dependent on the system size L at large W (while it was L-independent for smaller W ). As can be seen in Figs. 7 (d)-(f), our numerical results for C log (ω) do not confirm this finding. [Note that we have considered the same values of W and L as in the previous discussion of C(ω)]. Namely, we find that the log-averaged spectral function C log (ω) exhibits no notable dependence on L, both for weak as well as strong disorder. Moreover, C log (ω) is very similar to the standard average C(ω) shown in Figs. 7 (a)-(c). We explain the apparent discrepancy between our results and the findings from Ref. [73] by the different ways the log-averages are defined (see Appendix C for more details.)

IV. DISCUSSION
In conclusion, this work has shed light on the correspondence between quantum and classical dynamics, and moreover extended our understanding of the dynamics in disordered spin chains beyond the well-studied case of S = 1/2. Specifically, we have compared the infinitetemperature dynamics of equal-site correlation functions C(t) and their spectral functions C(ω) in disordered quantum and classical spin chains with S = 1/2, 1, 3/2. statistics of energy-level spacings, we have shown that the putative many-body localization transition is shifted to larger and larger values of disorder upon increasing the spin quantum number S. Especially for vanishing or small values of disorder, we found that the dynamics of C(t) and C(ω) becomes almost independent of S and agrees with the classical result, upon introducing an effective disorder W eff . Developing a better understanding and defining proper criteria where such a type of "universality" occurs, promises to be an interesting avenue of future research. This question is not only of conceptual relevance, but is also of practical importance if the dynamics of strongly correlated quantum systems can be obtained from a much less demanding simulation of a suitable classical spin model instead. In this context, it is also interesting to extend the comparison between quantum and classical dynamics to a wider class of observables such as, for instance, the full space-time profiles S z l (t)S z l ′ and their respective structure factors in momentum space.
While the very good agreement between classical and quantum dynamics with S = 3/2 was found to persist even at very large values of W eff (given the system sizes and time scales available), the mapping between different quantum numbers S at least partially breaks down at stronger disorder. Specifically, while C(t) and C(ω) exhibit distinct signatures of MBL in the spin-1/2 model, the dynamics for S = 1, 3/2 appear delocalized even at the strongest W eff considered by us. Clarifying the existence of a MBL transition and determining the exact critical disorder strength in models with S > 1/2 will require further numerical and analytical efforts in the future.

Acknowledgments
We thank M. Serbyn for fruitful discussions. This work has been funded by the Deutsche Forschungsgemeinschaft (DFG) -Grants No. 397067869 (STE 2243/3-1), No. 355031190 -within the DFG Research Unit FOR 2692.
Appendix A: Dynamics in clean models for additional model parameters In addition to the data presented in the context of Fig. 1, let us briefly discuss the dynamics of C(t) also for other disorder-free models. To this end, we extend the Hamiltonian (1) from the main text by considering an additional anisotropy in the z direction and a nextnearest neighbor interaction, i.e., H = J L l=1 H l now reads H l = S x l S x l+1 + S y l S y l+1 + ∆S z l S z l+1 + ∆ ′ S z l S z l+2 . (A1) In Fig. 8 (a), C(t) is shown for ∆ = 1.5 and ∆ ′ = 0. Note that also for this choice of ∆, ∆ ′ , the spin-1/2 model remains integrable. In comparison to Fig. 1, we find that the agreement between quantum dynamics at different S and classical mechanics becomes even better in the anisotropic model. In particular, the oscillations of the S = 1/2 curve are much less pronounced. To proceed, Fig. 8 (b) shows C(t) for a finite nextnearest neighbor interaction ∆ = ∆ ′ = 1, which breaks the integrability. While one might expect that such a breaking of integrability improves the agreement between S = 1/2 and larger S ≥ 1, the results in Fig. 8 (b) do not confirm this expectation. Specifically, we observe a pronounced dip of the S = 1/2 data at t S ≈ 2, which is not present for S = 1, 3/2 or classical spins.
Appendix B: Finite-size effects of the level-spacing ratio r In Fig. 9, we show a finite size scaling of the mean ratio r of adjacent level spacings. As can be seen in Fig. 9 (a) for the case of S = 1/2, the transition from Wigner-Dyson to Poissonian level statistics occurs more abruptly for increasing L (see also [6]). Albeit the accessible system sizes are very small for S = 1 and S = 3/2, a similar behavior can also be found for these values of S in Figs. 9 (b) and (c). Therefore, one can expect that the behavior of r discussed in the context of Fig. 2 where the |n and E n naturally depend on the specific disorder realization. Importantly, note that the logarithm in Eq. (C1) is outside the sum over the full Hilbert space. On the other hand, Ref. [73] considered the average of logarithms of individual matrix elements, ln | m| S z l |n | 2 , which are then binned in respective energy windows. These two quantities are not necessarily the same and, in particular, the correlation function (C1) appears to be insensitive to the breakdown of selfaveraging reported in [73].
Complementary to Fig. 7, we present in Fig. 10 additional results for the log-averaged correlation C log (t) in real time. For the two disorder values W = 1 and W = 3 considered, we find that while the standard average C(t) and the log-average C log (t) are slightly different from each other, their overall behavior is very similar (see also [41]). Moreover, both C(t) and C log (t) are almost independent of L, at least for the system sizes and time scales depicted.