Dynamic mean-field theory for dense spin systems at infinite temperature

A dynamic mean-field theory for spin ensembles (spinDMFT) at infinite temperatures on arbitrary lattices is established. The approach is introduced for an isotropic Heisenberg model with $S = \tfrac12$ and external field. For large coordination numbers, it is shown that the effect of the environment of each spin is captured by a classical time-dependent random mean-field which is normally distributed. Expectation values are calculated by averaging over these mean-fields, i.e., by a path integral over the normal distributions. A self-consistency condition is derived by linking the moments defining the normal distributions to spin autocorrelations. In this framework, we explicitly show how the rotating wave approximation becomes a valid description for increasing magnetic field. We also demonstrate that the approach can easily be extended. Exemplarily, we employ it to reach a quantitative understanding of a dense ensemble of spins with dipolar interaction which are distributed randomly on a plane including static Gaussian noise as well.

A dynamic mean-field theory for spin ensembles (spinDMFT) at infinite temperatures on arbitrary lattices is established. The approach is introduced for an isotropic Heisenberg model with S = 1 2 and external field. For large coordination numbers, it is shown that the effect of the environment of each spin is captured by a classical time-dependent random mean-field which is normally distributed. Expectation values are calculated by averaging over these mean-fields, i.e., by a path integral over the normal distributions. A self-consistency condition is derived by linking the moments defining the normal distributions to spin autocorrelations. In this framework, we explicitly show how the rotating wave approximation becomes a valid description for increasing magnetic field. We also demonstrate that the approach can easily be extended. Exemplarily, we employ it to reach a quantitative understanding of a dense ensemble of spins with dipolar interaction which are distributed randomly on a plane including static Gaussian noise as well.

I. INTRODUCTION
Nuclear magnetic resonance (NMR) has been an extremely important field for a long time. On the one hand, it constitutes a powerful analytical technique in physical chemistry [1][2][3][4] which helps to understand the structure of molecules on all levels from their primary structure to their tertiary structure. One the other hand, it is a technique which has enabled fundamental steps in quantum computing by taking spins S = 1 2 as quantum bits [5]. The latter development illustrates that the dynamics of the spin degree of freedom has gained enormous attention in particular in recent years. Closely related is the rapid development of the field of quantum sensing based on NV-centers in diamond [6][7][8][9][10][11][12][13] which behave similarly to an elementary spin [14].
A key phenomenon in this field is decoherence, i.e., the loss of coherence of a small quantum system in contact with a larger environment, often called bath. A generic approach to small systems in weak contact with a large bath is the theory of open quantum systems [15]. This is a powerful approach if the energy scales of system and bath are very different. If the bath correlations decay much faster than the system's dynamics quantum master equations reliably capture the physics, for example in radiative decay processes. If, however, the separation of energy scales is not given and the back-action of the system onto the bath cannot be neglected a quantitative description is notoriously difficult. In the context of the coherent control of spins, the small quantum system generically is a single spin. The decoherence can result from a fluctuating environment, for instance from stray magnetic fields or from phonons which may be fast. But very often it results from surrounding spins of electronic or nuclear origin. This is the case often relevant in NMR and in sensing by NV-centers. Then the back-action of the considered spin onto its neighboring spins is important and cannot be neglected. While we cannot provide a comprehensive review over all techniques applicable to spin systems, we present a brief overview of the most commonly used techniques. This allows us to highlight differences to the approach we are proposing in this article. For very small bath sizes of only a few spins the resulting problem can be tackled by exact diagonalization (ED). The Chebyshev polynomial expansion technique (CET) allows for substantially larger but still comparably small finite bath sizes [16,17]. If, however, the degrees of freedom of the bath are too numerous then brute force numerical approaches cannot be applied due to the exponential growth of the Hilbert space with the number of bath spins. For certain geometries such as chains and stars density-matrix renormalization [18,19] provides numerical alternatives. But the maximum times which can be reached are limited.
For approximately star-like topologies, like the one of the central spin model, cluster expansions [20,21] and related methods [22], linkedcluster expansions [23,24] and cluster-correlation expansions [25][26][27] are prominent approaches. But, these approaches become cumbersome for lattices with many different bonds. In addition, they represent expansions in time so that the reachable maximum time is limited by the complexity of the tractable maximum clusters. The eminent problem with (occasionally constrained) Monte Carlo sampling methods [28][29][30] is that statistical errors can become very substantial, see in particular Ref. [31] concerning the effects for up to N = 48 nuclear spins. By means of semi-classical or quantum mechanical master equations for the density matrix of the whole system [32][33][34][35], an access to the overall dynamics can be obtained which is typically easy to realize, but potentially suffers from the drawbacks of mean-field approaches if these are not based on small expansion parameters. arXiv:2107.07821v2 [cond-mat.stat-mech] 2 Feb 2022 Hence, alternative techniques are of significant interest. One useful observation is that the dynamics of spins and in particular the effect of a larger number of spins can be captured fairly well by their classical equations of motion [36][37][38][39][40]. This can be understood as an application of the ideas of the truncated Wigner approximation [41,42] whose foundations date back to the idea of Wigner that a part of the quantumness is captured by averaging over distributions of initial conditions [43]. But it is conceptually very difficult to extend this approach systematically to take more and more quantum effects into account. Apart from the smallness of Planck's constant no small parameter is apparent. In the present article we deal with dense spin systems where each spin interacts with a large number of other spins. In the limit where this number of interaction partners becomes infinite we derive a dynamic meanfield theory for the spin dynamics (spinDMFT) at infinite temperature, i.e., for completely disordered spins. As in all mean-field theories, spinDMFT comprises an effective single-site problem and a self-consistency condition. Similar to the case of the established fermionic dynamic mean-field theories [44] the time dependence of the mean-field is a crucial ingredient. It bears similarities to the mean-field approach for the Sherrington-Kirkpatrick quantum model in which spin glass behavior has been established [45][46][47][48]. A dynamic mean-field approach has also been used for ordered magnetic phases [49] for which, however, the couplings between the spins have to be scaled differently. After this introduction, we derive the spinDMFT in Sect. II in consecutive steps for an isotropic Heisenberg model and discuss details of the numerical implementation. Subsequently, we compare the results of spinDMFT for several systems with results obtained by CET and iterated equations of motion [50,51] in Sect. III. Sect. IV is devoted to the application of spinDMFT to two-dimensional spin ensembles in which the spins couple via dipolar interactions. In particular, we can continuously show how the well-known rotating wave approximation (RWA) becomes more and more reliable for increasing external magnetic field. In Sect. V we conclude the article and give an outlook to future directions of research. The appendices provide technical details of the derivation and the numerical implementation of spinDMFT including an analysis of the achievable accuracy in the numerical simulations.

A. Model and Definitions
For concreteness, we consider an isotropic Heisenberg model for an ensemble of spins with S = 1 2 at infinite temperature. The spins are subjected to a static and homogeneous magnetic field aligned in the z direction so that the Hamiltonian reads as Here and henceforth we use bold symbols for quantum mechanical operators and three-dimensional vectors are indicated by the arrow above the symbol. The properties of the underlying spin lattice are encoded in the couplings J ij = J ji . It is useful to introduce the operators of the local environments of each spin Using them the Hamiltonian can be expressed as The prefactor 1 2 occurs here to avoid double counting of the couplings. For later purposes, we also introduce the moments of the coupling constants J m,i := j . ( Note that we do not restrict the model to periodic lattices, but include arbitrary clusters. Both coordination numbers assess the number of spins that constitute the environment of site i. Considering only constant nearestneighbor (NN) interactions both numbers z i and z i equal the number of nearest neighbors z NN,i = z i = z i which is the usual coordination number. A common property of mean-field approaches is that they become exact in the limit z → ∞ [44]. Therefore, 1/z serves as a control parameter allowing us to systematically neglect terms in non-leading order in 1/z. Hence, we examine several quantities with respect to their scaling with the effective coordination numbers. We will show that the spinDMFT becomes exact in the limit of infinite z i and z i . As a consequence, the approach is not optimum for lowdimensional systems. However, since we consider effective coordination numbers instead of the standard one, not only the dimensionality but also the overall behavior of the coupling constants matters. Obviously, long-range couplings will increase the effective coordination numbers at given, fixed dimension. In Sect. IV, we demonstrate that in case of dipolar couplings, i.e., for weakly decreasing couplings with the distance, our approach is successful even in two dimensions. We establish the dynamic mean-field theory for spins (spinDMFT) for the introduced model. This is done in four steps: (i) We replace the local-environment operators by classical time-dependent random local mean-fields.
(ii) We argue that the dynamics of the local mean-field at site i does not depend on the dynamics of the single spin at site i.
(iii) We show that the mean-fields are normally distributed.
(iv) The defining moments of the normal distributions are linked to spin autocorrelations yielding a closed set of self-consistency equations.
Since we consider infinite temperature, the density operator is given by 1/d, where d is the dimension of the Hilbert space. Thus, any quantum expectation values are given by and, consequently, any correlation by where we have set = 1. In the next section, we undertake the first and the second step, (i) and (ii).
B. From the spin ensemble to an effective single site 1.
Step (i) We justify the substitution of the local-environment operators V i (t) by classical fields. For this to hold, it is crucial that the spin ensemble is dense so that quantum fluctuations of the environment are negligible relative to the classical dynamics. The argument is adapted from Ref. 19 and runs as follows. We consider the Frobenius norm of an operator defined by and apply it to for each component of the local-environment operator.
To assess the role of the coordination numbers, we assume that the J 2,i are of the same order of magnitude at every site. They set the relevant energy scale which one should think of being constant when scaling the coordination numbers, i.e., the individual couplings scale roughly like For a classical variable, any commutator would vanish. Hence, we study the commutators of the localenvironment operators and compare their norm to the one of the V i themselves for α = β; for α = β the commutator vanishes. Clearly, for diverging effective coordination number z i → ∞ the commutator vanishes relative to the norm of the operator. Hence its quantumness becomes negligible and the local-environment operators can be replaced by classical mean-fields V i → V i . Note that this is a very common phenomenon in quantum mechanics. Quantities which represent the collective properties of a large number of constituents behave classically. We stress, however, that this argument does not imply that the classical field is static. Hence, we avoid this oversimplification and take the mean-fields as classical, but time-dependent and dynamic. A potential correlation between V i and S i is not ruled out at this stage.

Step (ii)
Here, the aim is to show that the dynamics of the individual spin at site i does not influence the dynamics of V i in the limit of z i → ∞. The basic idea is simple: a single spin contributes only negligibly to the large sum defining V i . But it is not straightforward to cast this idea into a formal argument. What we want to show is that the dynamics of S i does not influence the dynamics of V i , i.e., that no back-action needs to be taken into account. Indeed, we show in Appendix A that the correlation between the spin at site i and at an adjacent site j scales like 1/z for the special case of a Bethe lattice with NN interaction, where z = z i = z i ∀i holds. Hence the correlation between the spin at site i and its local environment V i scales like J ∝ 1/ √ z and becomes negligible for z → ∞. The number of spins in V i scales like z compensating the factor 1/z from the correlations. We stress that this conclusion is subtle. It is valid if the dynamics of V i is induced by a process of order z 0 because the relative error then is indeed of order z −1/2 . But if there is no process inducing a dynamics of order z 0 this does not hold true. Indeed, the central spin model (CSM) provides an instructive example. In this model, a central spin is coupled to a large number of bath spins, but the bath spins are not coupled among themselves wherein C i are arbitrary coupling constants and P denotes the so-called Overhauser field. This looks like a perfect scenario for replacing the P by a classical Overhauser field P with a given dynamics imposed on the central spin. Yet, this approach fails [19,38,52] because the Overhauser field has no dynamics if the central spin is taken out and treated separately. This reasoning shows us that it is essential at this stage to deal with a dense spin system where each spin interacts with many others so that excluding one of them from the dynamics hardly changes the dynamics of all others. This is the case if the coordination number at each site is large. We point out that this is clearly not the case in the CSM where excluding the central spin brings all dynamics to a halt and where the coordination number z for each bath spin is only 1. On the basis of the above arguments, we replace each local-environment field V i in the Hamiltonian (3) by an a priori given dynamic mean-field V i (t) so that the spin dynamics is given by These mean-field Hamiltonians, labeled by the subscript "mf", only contain linear spin-operator terms so that the spin dynamics at a given site is decoupled from the one at other sites once the mean-fields V i (t) are given as function of time. In order to emphasize that not only single values V i (t) are meant, but the whole time dependence we introduce the shorthand V i for it. These time series are drawn from a so far unknown probability functional p V 1 , ..., V N which we will determine in the next section. In conclusion, the original 2 N -dimensional ensemble is mapped to N two-dimensional quantum impurity systems each capturing a single spin subjected to a timedependent mean-field and the external field.

C. Distribution of the mean-fields Vi
Here we carry out step (iii) and step (iv) from the list in Sect. II A.

Step (iii)
The central argument is again that two spins at site i and j are only weakly correlated. This is difficult to show for an arbitrary cluster with arbitrarily linked spins. But for the Bethe lattice with NN interaction we demonstrate in Appendix A that the correlation S α i (t)S β j (0) scales like z −||i−j|| where ||i − j|| is the number of NN links needed to go from site i to j. Moreover, we show that the correlation V α i (t)V β j (0) is suppressed at least like 1/z for i = j. Thus, we conclude that the time series of the local mean fields V i are independent at different sites This allows us to compute any local expectation value by where . . . (sts) Vi stands for the expectation value for a given single time series V i . This contribution is weighted by the probability p i V i to reach the total average. The for the initial condition U i (t 0 , t 0 ) = 1. For future use, it is worth mentioning that the unitary evolution operator U i (t, t 0 ) only depends on the mean-field time series between t 0 and t and not on all times. Hence, the computation of a time-dependent expectation value only requires to average over mean-field time series within the relevant time interval. For the case in (13), it turns out that no averaging over time series is necessary at all. Since we assume that the system is completely disordered at t 0 , i.e., the density matrix is proportional to the identity, the unitary evolution cancels out and one arrives at the second line (13b). Averaging actually does not matter here. For correlations, the evolution in time does matter and hence does the averaging over the time series. We consider which clearly only depends on the time interval [t 1 , t 2 ]. Note that temporal homogeneity is not given for a single time series, but it holds on average so that The next, important conclusion is that each local environment V i is the sum of a large number of essentially independent spins (2). This number becomes infinite for diverging coordination number so that the central limit theorem applies and we conclude that the V i are normally distributed. This means that we need only two moments, the first and second, to determine the distribution. This brings us to the fourth and final step.

Step (iv)
We establish self-consistency conditions which link the first and second moments of the normal distribution to quantum expectation values and correlations. For the first moment, it is straightforward to see from Eq. (13) that it vanishes for all sites i, all components α ∈ {x, y, z}, and all times t because we start from the disordered, T = ∞ case where the expectation values of all spin operators vanish. Hence, we conclude that the distribution p i ( V) is a normal distribution with vanishing first moments This is the first self-consistency condition which is easy to fulfill. The spinDMFT can also be extended to include non-vanishing first moments and even time-dependent moments, but this is not the scope of the present article. For the second moments, we consider where the second line results from the fact that the spinspin correlations between different sites vanish in the limit of infinite coordination number. The third line, finally, results from (15) on average. Self-consistency requires that the second moment computed above equals the correlations of the local mean-fields, i.e., This closes the set of self-consistency conditions. If the second moments only depend on the time difference t 1 − t 2 the resulting two-time spin expectation values only depend on t 1 − t 2 . Hence solutions homogeneous in time exist. Whether they are the only conceivable solutions is an additional question which we do not study in this paper and leave for future research. At this stage, we observe the interesting feature that the resulting mean-field theory is the same that one would obtain for classical spins of the same average length. This is so since the effective single-site problem in Eq. (11) only contains the spin operator linearly. According to the Ehrenfest theorem the quantum mechanical expectation values behave identical to classical variables. We conclude that the classical and quantum mechanical spin system converge to the same spinDMFT for infinite coordination number. We emphasize, however, that for spins larger than 1/2 non-linear local terms may arise, for instance from quadrupolar couplings [53]. Then there is a difference between the quantum and the classical spinDMFT. Henceforth, we use the term 'autocorrelation' to denote the local spin-spin correlation S α i (t 1 )S β i (t 2 ) . Later, when numerical results are presented, we will also distinguish between diagonal autocorrelations (α = β) and cross autocorrelations (α = β). The message of Eq. (15) is that one can compute the autocorrelations at each site if one knows the moments V α i (t)V β i (0) defining the normal distribution of V i . In return, Eq. (19) tells us that the knowledge of the autocorrelations of the spins linked to site i yields the second moments of V i . It is to be expected that this closed set of self-consistency equations can be solved iteratively and our numerical results confirm that this is true. Numerical aspects will be discussed in the next section. In the present paper, we do not intend to use spinDMFT for problems with spatial dependence even though the general formalism derived so far allows for such spatial dependencies. But the concomitant numerical task is quite demanding. Our goal here is first to introduce the approach of spinDMFT and to illustrate its performance. To this end, we opt to consider homogeneous spin ensembles where each site is equivalent to every other site. Certainly, this is the case for periodic lattices but it can also hold for dense random spin ensembles where each spin is interacting on average with the same number of spins and with the same interaction strength. Then, all autocorrelations are the same and hence all second moments of the local mean-fields. Then the self-consistency condition (19) simplifies considerably because the autocorrelations on the right hand side can be taken out of the sum. The site-independent second moments read as Since all sites are equivalent, no site indices need to be denoted. Interestingly, the only energy constant governing the spin dynamics aside from the external magnetic field is the root-mean-square J 2 of the couplings.

D. Numerical implementation
Our aim is to implement a numerical procedure by which the mean-field moments determined by the self-consistent equations can be evaluated. The basic idea is to start with some arbitrary initial function and to converge iteratively to the solution. In each iteration step, one computes the autocorrelations for a single spin via the path integral (15b) and subsequently the mean-field moments via the self-consistent equations (20). This scheme is illustrated in Fig. 1.  The computation of the path integral constitutes the numerical challenge. First, we need to discretize the time so that the number of second moments becomes finite. Henceforth, we set t 0 = 0 and choose an equidistant discretization [t 0 = 0, ...t L ], i.e., t l = lδt. The numerical error resulting from the discretization is discussed and analyzed in appendix B 2. We obtain a [3(L + 1) × 3(L + 1)]dimensional covariance matrix M of which the matrix elements are No site label i occurs because we treat a homogeneous spin ensemble so that the covariance matrix is the same at each site. But a possible generalization to a spatial dependence is obvious. If M is known the corresponding normal distribution reads as To be specific, the vector-matrix-vector product in the argument of the exponential function stands for Second, we use a Monte-Carlo method to carry out the path integral: we draw a time series V from the distribution function, compute the expectation value for each time series and finally calculate the arithmetic mean of the results. This is done sufficiently often to achieve a small enough statistical error which is studied in detail in App. B 1. The strategy to determine the mean-field moments is set up as follows: 1) Choose arbitrary functions for the initial second moments of the mean-fields in the studied time interval.
3) Draw a large number of time series for the meanfield according to the distribution (22). 6) Evaluate the iterated mean-field moments from the self-consistency conditions (20) and return to step 2) or stop if convergence of the second order meanfield moments is achieved within a given tolerance.
For step 2), it is convenient to set up the covariance matrix in blocks depending on the spin components α, β ∈ {x, y, z}: Spin symmetries of the system can easily be exploited to reduce the numerical effort. For instance, for zero magnetic field any block with α = β vanishes, so that the covariance matrix becomes block-diagonal. Furthermore, we stress that M is symmetric. This is actually required for a covariance matrix. Here, it results from the physics at infinite temperature: the quantum expectation values and hence the mean-field moments are symmetric due to the cyclic invariance of the trace. Another crucial property of covariance matrices is their positive semidefiniteness. In App. B 5 we show that M is automatically positive definite because it results from the quantum expectation values of Hermitian operators. In App. B 4 we explain how including time-translation invariance in the algorithm reduces the numerical effort further. Another algorithmic issue is the sampling procedure in step 3). Since the covariance matrix is generally nondiagonal, the mean-fields at different times cannot be drawn independently of each other. Hence, it is indicated to first change basis such that M is diagonal in the new basis. In this basis, for each vector component an independent random variable can be drawn from a one-dimensional normal distribution. Subsequently, we transform back into the original basis obtaining the desired autocorrelation in time. We recommend the following strategy: a) Diagonalize the symmetric, non-negative covariance matrix by the orthogonal transformation O b) Sample a 3(L + 1)-dimensional vector R of uncorrelated Gaussian random numbers in the diagonal basis. Each component has a zero average and a variance given by the corresponding eigenvalue of M , i.e., the corresponding diagonal element of D.
c) Transform the random vector R to the original basis The diagonalization needs to be performed only once in each iteration step since all drawn time series belong to the same covariance matrix. In contrast, steps b) and c) have to be performed for each drawn time series. To compute the time evolution operator, or propagators, in step 4) numerically we split it into a product of propagators over the short time interval between consecutive t l , i.e., over These propagators can be computed efficiently by commutator-free exponential time propagation (CFET) [54]. Since we do not have information about H (mf) at times between two consecutive t l any integral can only be approximated by trapezoidal rule. Therefore, the error of each propagator is at best of order δt 3 so that CFETs of orders larger than two appear pointless. From our numerical experience, we recommend a second-order and an optimized fourth-order CFET [54] where In step 6), one requires an exit condition to decide when a sufficiently converged result has been found. A possible choice is to compute the deviation between the results of current and previous iterations and compare it to a chosen tolerance threshold. If the deviation falls below the tolerance threshold, the iteration is stopped. We discuss the definition of the deviation and the choice of the tolerance in App. B 3. In general, when we graphically show numerical results of spinDMFT, we choose the numerical parameters such that the resulting errors are not larger than the thickness of the lines, if not explicitly discussed otherwise. As mentioned before, Appendix B provides a closer insight into the error sources. In the following sections, we examine the validity of spinDMFT by comparing its results to the ones of established numerical techniques.

III. COMPARISON OF spinDMFT TO OTHER APPROACHES
Before applying the advocated spinDMFT to various physical systems it is advisable to compare results of spinDMFT with results of different well-established methods. Since the main idea of spinDMFT is based on a large number of interaction partners we expect the agreement to become the better the larger the coordination number of the spin ensemble is.

A. Methods for comparison
Below, we use two methods to obtain results for comparison. The first method is the Chebyshev expansion technique [16,17] (CET), the second method is the iterated equations of motion (iEoM) approach [50,51]. The CET is numerically exact up to a systematically controlled error threshold. The iEoM approach is an approximate approach controlled by the number m of iterations performed.
To obtain the time dependence O(t) of an observable using CET we expand the unitary time evolution operator U = e −iHt in terms of Chebyshev polynomials defined recursively by All polynomials T n are defined on the closed interval I = [−1; 1]. To ensure that the energy spectrum of a given Hamiltonian H lies in I we rescale the Hamiltonian according to H → H = (H − b)/a. Then, the Chebyshev polynomials can be used as an orthogonal functional basis. In order to perform the rescaling an estimate of the extremal eigenvalues [55][56][57] of H is needed to obtain a = (E max − E min ) /2 and b = (E max + E min ) /2. Rough estimates in the form of upper (lower) bounds for E max (E min ) are sufficient because the rescaling only has to ensure that the rescaled eigenvalues lie within I. Subsequently, the expanded time evolution operator reads as where the time-dependent coefficients contain the Bessel functions of first kind J n (at). Given an initial state |ψ 0 its time evolution reads as Here, the basis states of the expansion are |φ 0 := |ψ 0 , |φ 1 := H |ψ 0 , and |φ n+1 := 2H |φ n − |φ n−1 .
In the numerical implementation, the infinite series (32) must be terminated at some finite value N c < ∞. The time dependence of the prefactors is essentially determined by the time dependence of the Bessel functions J n (t) [58]. The higher the order n the longer it takes the Bessel function J n (t) to contribute noticeably to the series. Given the cut-off N c of the series the truncation error of the CET series is estimated by Note that the truncation error is not only related to the cut-off N c , but also depends on the maximum time up to which results are calculated as well as on the parameter a which equals half the width of the energy spectrum. The important property of the CET is that N c , required to keep the error low, increases only linearly with the time t up to which one intends to compute the evolution. The second method we employ for comparison is the iEoM approach [50,51] which approximates the time dependence of an operator in the Heisenberg picture. Starting with an arbitrary operator A 1 one expands where all time dependence is incorporated in the complex prefactors h i (t with the Liouville superoperator L. Expanding the result of L(A i ) in terms of the chosen basis {A i } by means of leads to the Liouvillian matrix L, also called dynamic matrix. For a compact notation the time-dependent prefactors h i (t) are combined to a vector h(t) of which the dynamics is obtained by inserting both expansions (34) and (36) in (35) yielding The Liouvillian matrix is most easily computed for an orthonormal operator basis {A i } (ONOB) so that each matrix element is given by As previously argued [50,51], it is crucial to achieve Hermiticity of L to avoid exponentially diverging solutions which are unphysical. The Hermiticity of L is equivalent to the self-adjointness of L which depends on the used operator scalar product. A convenient choice is the Frobenius scalar product Due to the invariance of the trace under cyclic permutations L is indeed self-adjoint and thus L is assured to be Hermitian [50,51]. The ONOB is found iteratively by applying the Liouville superoperator m times which is called the loop order. Starting from a spin operator at a given site, the application of L creates more and more increasingly complicated expressions which are sums of operator monomials, i.e., sums of products of local operators. The number of such monomials is finite for all m, but grows exponentially for increasing m. For a spin S = 1 2 the site local operators are those given in Tab. I. After m iterations monomials involving up to m + 1 sites occur. This means that in these monomials Pauli matrices at up to m + 1 sites can occur. The identity operator is trivial and does not need to be tracked. If A = A 1 is the initial spin operator the initial vector h has the components

B. Observables and symmetries
Since we consider infinite temperature, any expectation value of a single-time observable is actually timeindependent, see Eq. (13). Therefore, the primarily interesting observables are the spin autocorrelations which we denote by There are nine different autocorrelations of the above type due to the choices for α, β ∈ {x, y, z}. However, the symmetries of the Hamiltonian imply a number of relations between them so that only a small number needs to be considered. We briefly discuss the symmetries of the system in the following. The original Hamiltonian (1) is invariant under any spin rotation around the z axis, in particular about the angle π/2 implying As a consequence any correlation between the transversal and longitudinal spin components disappear while the transversal cross autocorrelations g xy and g yx fulfill By means of cyclic permutations in the trace and homogeneity in time we additionally derive The transversal diagonal autocorrelations are equal In case of zero magnetic field, the system is also invariant under time reversal because the Hamiltonian is bilinear in spin operators which implies so that all cross autocorrelations vanish in this case. Furthermore, the diagonal autocorrelations g αα are equal due to complete isotropy of the model.    For short times and zero magnetic field, a Gaussian fit describes g αα very well in the linear plot in Fig. 2. Some deviation is discernible from intermediate times onwards.
To analyze this deviation in more detail the functions are plotted in Fig. 4 on a logarithmic scale vs. t 2 . Interestingly, the diagonal autocorrelations appear to show Gaussian behavior at short and at long times, but with different standard deviations. For longer times, the decay is slowed down.

C. Comparison to results of other approaches
We compare results from the spinDMFT to results from exact diagonalization (ED), iterated equations of motion (iEoM), and Chebyshev expansion technique (CET). ED is a very well-known technique and the latter two approaches have been explained above. A conceptual difficulty lies in the fact that these alternative techniques work best for small and low-dimensional systems while spinDMFT is rather justified in large, high-dimensional systems. But comparing results from spinDMFT to these alternatives is the best option at hand. Note that such comparisons are particularly challenging for spinDMFT. Considering the self-consistency problem (20), we stress that all lattice properties are embodied in a single coupling constant, namely the root-mean-square J 2 .
No site index appears because we deal with homogeneous systems. Time is naturally measured in units of 1/J 2 . First, we consider one-dimensional (1D) spin chains with S = 1 2 . For finite pieces of chains with N sites, periodic boundary conditions (PBC) are taken. The Hamiltonian in the isotropic case reads which entails Figure 5 compares the results from the above mentioned methods to the data obtained from spinDMFT. The ED data are taken from Ref. 59. The results from ED and CET coincide very nicely in the considered time interval. Moreover, no finite-size effects are discernible in this interval. For not too long times, the iEoM result also matches very well. It has the advantage to consider the infinite system by construction. These results almost coincide with the spinDMFT data until roughly t ≈ 3/J 2 . The subsequent deviations can be attributed to the small coordination number of the 1D chain with z = 2 which, obviously, is a challenge for a mean-field approach. The spinDMFT shows quick and rather complete decoherence while the genuine 1D results show weak coherent revivals at t ≈ 5/J 2 and t ≈ 9/J 2 . This is not surprising because the integrable 1D system is strongly constrained in its dynamics due to its extensive number of conserved quantities [60]. Figure 6 compares the results of CET and iEoM in 2D, i.e., for the isotropic Heisenberg model on the square lattice with NN coupling J to the data obtained from spinDMFT. In this case, J 2 = 2J holds. The agreement between CET and iEoM is good up to t ≈ 3.5/J 2 ; then, the effects of finite loop order m kick in. In 2D, it is unfortunately not possible to reach larger values of m. Up to this range, the spinDMFT is in nice agreement with the other approaches as well. What is even more interesting is to see the evolution from 1D to 2D. To this end, we include the CET result in 1D. Clearly, passing from 1D to 2D improves the agreement between the genuine numerical results and spinDMFT. This is exactly what one had to expect in view of the derivation of spinDMFT as mean-field theory which becomes exact for infinite coordination number. Hence, this observation constitutes a good confirmation of the validity of spinDMFT. In order to corroborate the foundation of spinDMFT further we consider the Heisenberg model on complete graphs, i.e., graphs where each site is connected to all other sites Such graphs or clusters are called infinite-range clusters in the physical literature. The Heisenberg model on such graphs is highly symmetric if the coupling is the same for all links. This leads to rather special autocorrelations. In order to avoid features from non-generic high symmetries we consider a random model where the couplings are drawn from a Gaussian distribution. Then, they are normalized, i.e., multiplied by a suitable constant λ j , such that holds for all j. The results for the autocorrelations are averaged over 100 sets of random couplings.  compares the CET results to the data from spinDMFT for various values of N . The symbols display the data extrapolated to N = ∞ by a linear fit in 1/N 3/2 of the data for the three largest values of N . This particular power law fit is chosen in view of the scaling of each J ij ∝ 1/ √ N stemming from the normalization (52). This suggests to use power law fits ∝ 1/N p/2 with some integer p. We found that p = 3 is most robust.  The obvious trend is that the data for finite N appear to converge to the spinDMFT. This observation again justi-fies the systematic derivation of the advocated dynamic mean-field theory. Finally, we mention that the Ising model on complete graphs in the limit N → ∞ has the autocorrelations [61] 4g xx = e − 1 The energy scale is J 2 = √ N − 1J. These results are reproduced by spinDMFT; we refrain from displaying them because the graphs coincide perfectly. Due to the spin anisotropy the self-consistency conditions are changed to On the basis of the above results, we conclude that spinDMFT is a systematically controlled dynamic meanfield approach to disordered spin systems at infinite temperature which becomes exact for infinite coordination number. It is designed to provide quantitative information of the local spin dynamics. It is a valid approximation for finite coordination numbers which nicely captures essential physics such as rapid decoherence, spin anisotropies, and Larmor precession. Due to the required moderate computational resources it is an attractive tool to understand spin dynamics in various setups. We illustrate this last point by applying spinDMFT to a spin ensemble with dipolar interactions.

IV. APPLICATION TO A DIPOLAR SURFACE SPIN ENSEMBLE
In this section, we want to show that spinDMFT can be adapted to models which describe experimental setups or are very close to experimental questions. We illustrate that spinDMFT can be applied to complex physical situations because of its flexibility and, furthermore, that the resulting numerical task is feasible and does not require excessive compute resources. The model which we will address is motivated by intensive studies of localized defect spins of electronic origin with S = 1 2 and a g-factor of g ≈ 2 on diamond surfaces which are observed by NV centers [7,11,13]. These spins were seen and examined in recent studies [12,62,63]. The precise origin of the defect spins is still a matter of debate although recent progress indicates that they are formed by trapped electrons very close to the surface [64]. They appear to be inhomogeneously distributed over the surface [7,13,65]. Driving these spins reduces decoherence in shallow NV centers [66]. The surface spins interact with one another and they are subject to additional, slow noise. The origin of the latter is not yet clarified: nuclear proton spins are candidates [11,67] which agrees with the importance of the precise chemical and morphological conditions at the surface [65]. A second candidate is phonons which also play a role [63]. In addition, NV centers can also measure the dynamics of 13 C nuclear spin baths which are distributed three-dimensionally in the bulk of diamonds [68]. Here, we do not aim at describing one of the above exciting experiments in detail, but to address a generic model comprising the essential features. To this end, we consider a random ensemble of localized electronic spins S = 1 2 on a planar surface interacting by dipolar couplings. Aside from these interactions, the spins are subjected to a global magnetic field as well as to local static magnetic field noise. The latter can be viewed as being generated by slowly fluctuating nuclear spins stemming, for instance, from protons. Figure 8. Sketch of the considered dipolar spin ensemble on a planar surface. The magnetic field B (red) defines the zdirection; the isotropic random local magnetic fields bi (blue) vary in strength and direction from site to site.

A. Model and Definitions
We consider an ensemble of spins S = 1 2 distributed randomly over a planar surface which interact via dipoledipole interaction. Each spin is subjected to an external magnetic field B whose direction defines the z direction which forms an angle ϑ with the normal n of the plane. Additionally, the spins see local random magnetic noise, i.e., local magnetic fields b i with zero average. While the system is inhomogeneous we assume that the random distribution is such that each spin sees the same environment on average, i.e., the system behaves on average like a homogeneous system. This means that the average quantities do not depend on the site. The Hamiltonian is given by where is the dipolar coupling. We apply spinDMFT which relies on average dynamic mean-fields. Since we want to treat the random local magnetic fields b i in addition we have to distinguish three types of averages: (i) the one from spinDMFT which we denote by an overline as before, (ii) the one solely due to the average over the local magnetic fields which we denote by an overline with index 'n' for 'noise', and (iii) the complete average comprising (i) and (ii) which we denote by an overline with index 'c' for 'complete'. For simplicity, we assume that the local magnetic fields are distributed isotropically according to a normal distribution defined by the moments and The latter implies that the local fields are independent from one another. The distribution reads as As mentioned above, we allow for an angle ϑ between the surface normal n and the external magnetic field B.
We introduce the in-plane polar coordinates R ij , ϕ ij to express the distance vectors between site i and j by The complete system is sketched in Fig. 8 including the various introduced quantities.

B. spinDMFT
As motivated in the previous section where we introduced spinDMFT we define local operators describing the environments of the spins where the anisotropies are incorporated in the matrix With their help, the Hamiltonian can be rewritten where a factor 1 2 is introduced to avoid double counting. From the derivation of spinDMFT in Sect. II we know that large coordination numbers provide the justification for this dynamic mean-field theory. Thus, we consider the effective coordination numbers z and z defined in (5)  As expected, the long range of dipolar interactions increases the effective coordination number to considerably larger values compared to the NN coupling, see for instance the triangular lattice. This effect is larger for z than for z because the sums for higher moments converge faster than those for the second or first moment. We point out that the first moment "just" converges like RdR R 3 ∝ 1 R . For randomly distributed spins the issue is more complicated and depends on how the spins are located on the surface. In case of random positions without any restrictions, the effective coordination numbers are fairly small because there is a high probability for each spin to have a single neighbor close to it which dominates z and z . We found for completely random distributions z ≈ 1 − 10. Then, this close partner governs the dynamics and the application of spinDMFT is not well justified. However, considering restrictions for the location of the spins, in particular a minimum distance between the spins, the effective coordination numbers increase substantially to about the values of the triangular lattice which is the closest-packed lattice in two dimensions. We found z ≈ 5−15. We emphasize that such restrictions are very plausible: a minimum distance can result from the surface structure which does not allow the spin-carrying adatoms to be located very close to one another. In addition, a local repulsion between them would also ensure a minimum distance between the spins. We conclude that the lattice is dense enough so that spinDMFT is justified. Next, we replace the local-environment fields by dynamic mean-fields considering the local mean-field model The self-consistency conditions need to be complemented by the effect of the random noise fields b i . Averaging has to be done over the random time series for V i and the random local magnetic fields. We perform this in a single step and thus pass to combined fields and perform a single average over the combined distribution where the modified covariance matrix is given by Then, the noise leads only to an offset in the second moments which is constant in time.
The self-consistency condition of the first moment is still trivial For the second moments we consider the self-consistency and thus the complete self-consistency reads as This equation still constitutes a challenging numerical issue because it amounts up to a self-consistency condition for each spin. But, as stated at the beginning of Sect. IV A, we assume that the system is dense enough to be treated on average as a homogeneous system. Essentially, this means that J 2,i takes the same value at each site i. Then, the site indices can be omitted and we obtain much simpler self-consistency conditions The constants J and χ αβ ργ embody the energy scale and the spin anisotropies. The key idea is to approximate the discrete sums in the self-consistency conditions by integrals assuming a continuous distribution of spins with density n 0 = 1/r 2 min . Of course, this is not exact, but it provides a well-justified quantitative relation between the dipolar interaction in Eq. (57) and the prefactors of the self-consistency condition (73). We replace the sum by the integration Since we are dealing with a system constant in time we study self-consistent solutions which depend only on the time difference t 1 − t 2 . Hence, from now on we only consider correlation functions with t 1 = t and t 2 = 0, In the remainder of this section, we specialize the above general equations to the case of a perpendicular magnetic field, i.e., ϑ = 0, for simplicity. From Eq. (75b) we obtain by straightforward analytic calculation while the other coefficients vanish. For a brief symmetry discussion, we consider the original Hamiltonian (56). Since ϑ = 0, the transversal x and y spin components lie in the plane of the surface. Thus, the dipolar interaction term and the magnetic-field term are invariant under a rotation in spin and real space about the angle π/2 around the z axis. This does not hold true for the noise term ∝ b i . But the noise distribution (60) is isotropic so that on average this term also remains invariant and we have a rotational symmetry of the total system. In particular, this implies g xx (t) = g yy (t) (78a) g xy (t) = −g yx (t) (78b) g xz (t) = g zx (t) = g yz (t) = g zy (t) = 0.
(78c) Summarizing, we obtain the self-consistency equations It is worth mentioning that the noise explicitly appears in these equations because we included it in the meanfield W . The magnetic field, on the other hand, enters the physical problem in the Hamiltonian (66). Another important observation is that the transversal and longitudinal equations (79a) and (79d) show different prefactors. Certainly, this leads to different behavior of the corresponding autocorrelations. For zero magnetic field, where the cross autocorrelations g xy = −g yx vanish due to time reversal symmetry, we expect g xx = g yy to decay slower than g zz , since the transversal mean-field is stronger so that the z components are more strongly precessing. Henceforth, we measure the time in units of 1 J and it makes sense to define a dimensionless noise strength and a dimensionless magnetic field (80b) Fig. 9 shows the numerical results of the self-consistent equations (79) for various noise strengths and zero magnetic field. Considering C = 0.0, we observe a clear difference between the transversal and longitudinal signal. This anisotropy is expected due to the different prefactors in (79a) and (79d) as mentioned before. For C = 1.0, the difference is still present, however, both curves show a similar trend: a local minimum at the beginning followed by a rather slow decay. As we increase C, the anisotropy further diminishes because the isotropic noise contributions in (79a) and (79d) dominate more and more over the dipolar terms. In case of a very large noise strength, the mean-field contributions can be neglected and the remaining dynamics can be solved analytically [19,69] 4g αα stat (t) = Inspecting Fig. 9, the data approach the curve (81) for large values of C and short times. At larger times, i.e., beyond the local minimum, the dynamics due to the dipolar couplings makes itself felt and the signal decays below the analytical plateau. This is not surprising since the analytical consideration only includes static noise neglecting any mean-field dynamics. The corresponding coupling may be weak compared to the noise, but it certainly affects the long-time behavior of the signals. of g xx disappear very early although the signal is still finite. Subsequently, this correlation shows a slow longtime decay without discernible precession. We attribute this behavior to the presence of transversal noise stabilizing g xx = g yy . This noise component is certainly weakened by the longitudinal magnetic field, but it appears to be still strong enough to keep the signal finite for quite a while. Remarkably, the combination of noise and magnetic field causes a slow down of the longitudinal autocorrelation. This behavior is studied in detail in the next section where we use the RWA to tackle the problem for considerably larger magnetic fields. Since the transversal noise vanishes for such large fields, we expect the transversal long-time signal to vanish very quickly.

C. Strong-Field Regime and RWA
In the preceding section, we treated a general external magnetic field of arbitrary strength, weak or strong, but perpendicular to the plane. In this respect, the situation was specific. In the present section, we choose the angle ϑ of the external field with the surface normal in an arbitrary way, but consider a strong field so that the RWA is valid. First, we switch from the laboratory frame to the frame rotating with the Larmor frequency of precession This leads to a time-dependent effective Hamiltonian with oscillating terms. They oscillate the faster the stronger the magnetic field is. The RWA consists in averaging these fast oscillations yielding an effective timeindependent Hamiltonian again. The spinDMFT is applied to this effective time-independent Hamiltonian. This requires to solve a closed set of self-consistency equations capturing the spin dynamics in the rotating frame.
In order to consider the system in the Larmor rotating frame one has to rotate any observables of the lab frame backwards by the unitary evolution operator where the Zeeman term is the last-but-one term in Eq. : The spin operators are given in the rotating frame by The full time evolution in the rotating frame results from the complete time evolution operator Its Schödinger equation is derived by inserting (86) in the original Schödinger equation with the Hamiltonian Clearly, the Zeeman term H Z in H is canceled in this way which was the goal of this transformation. The remaining terms, i.e., the dipole interaction and the noise are rotated by U Z . As a consequence, the Hamiltonian H rot (t) is strongly time-dependent comprising fast oscillating terms such as cos(ω L t) or sin(ω L t) which are averaged in the sense of a Magnus expansion [70] in first order. The neglected terms are smaller by a factor 1/ω L . Put simply, the larger the Larmor frequency ω L is relative to the typical dipolar interaction frequency ω DD = J / , the better the RWA is justified. Thus, the strong-field regime is realized for In this regime, one replaces all fast oscillating terms in the Hamiltonian by their average values, i.e., In our case, we obtain which is again time-independent by construction. Note that the transversal components of the magnetic field noise are eliminated while the longitudinal component remains unchanged. Therefore, one expects growing differences between transversal and longitudinal autocorrelations upon increasing the noise strength σ N .
Next, spinDMFT is applied to the Hamiltonian (91). We only highlight expressions which differ from the previous application. Expressed in local-field operators, the Hamiltonian reads as where D rot is given by Replacing the local-field operators by the corresponding mean-fields leads to the local mean-field Hamiltonian and we again combine noise and mean-field in since both follow from Gaussian distributions. While the first moment is still zero, the second moments obey The coupling J is given by Eq. (75a) and the coefficients χ rot result from Eq. (75b) with D replaced by D rot . Defining the function of the polar angle ϑ we can express the prefactors concisely by Any other coefficient vanishes because D rot is diagonal. We reconsider the symmetries of the underlying system to be able to formulate the minimum set of self-consistency conditions. The original rotating-frame Hamiltonian (91) is invariant under spin rotation around the z axis by construction: all transversal spin components only occur in pairs. The dipolar part is invariant under time reversal. While this does not hold for the noise term for each individual b i , their distribution remains unchanged under time reversal. Spin-rotation symmetry and time-reversal symmetry allow us to conclude This enables us to reduce the general self-consistency conditions (96) to Note that there is a natural anisotropy between the transversal and longitudinal equation again, but this time by a factor of four. Furthermore, the noise only acts in the z-direction. Considering the results of the previous section, we expect even bigger differences between the autocorrelations here. Henceforth, we no longer use the terms "diagonal" and "cross" because no cross autocorrelations appear in the RWA dipole model. If not stated otherwise, we set ϑ to the so-called magic angle ϑ magic := arcsin 2 3 I(ϑ magic ) = 1 2 in our numerical calculations. Figs. 12 and 13 show our numerical findings for the solutions of the self-consistent equations. Analyzing them, we conclude two important facts: (i) Because of the natural anisotropy and the noise, the longitudinal signal decays considerably slower than the transversal signal. Increasing the noise strength amplifies this difference.
(ii) The transversal signal decays very accurately following a Gaussian, while the longitudinal decay is weaker than exponential. Considering the results, the first statement is rather obvious. For C = 0.0, the longitudinal signal already survives significantly longer than the transversal one due to the anisotropic factors in (100). Increasing the noise strength causes the spin to precess more and more quickly and randomly about the z axis. As a result, the transversal spin components experience an even stronger decoherence than before and the corresponding autocorrelations decay faster. In contrast, the z component of the spin is stabilized by the additional rotations about the z axis so that the longitudinal signal relaxes only very slowly. converge to zero. We emphasize that this negative dip is only very small, ∝ 10 −3 , and hence not visible in the provided figures. Still, we ensured that it does not result from numerical inaccuracies. Interestingly, its height decreases upon increasing the noise width. The clearly linear behavior in the plot has a slope r = 2 so that it clearly indicates Gaussian behavior. Hence, we use the fit function to extract the standard deviation σ as function of the noise strength displayed in Fig. 16. This behavior can be understood by an analytical consideration [19,69]. We consider purely static noise neglecting the mean-field contributions in the RWA Hamiltonian (94): The analytical averaging yields the transversal signal which explains the Gaussian behavior of the transversal signal for large values of C. The analytical standard deviation is also depicted in Fig. 16 as blue dashed line. For increasing noise strength it describes the fitted data better and better. This is expected because the larger the static noise is relative to the mean-fields, the better the static noise model captures the dynamics. For small values of C and in particular for C = 0.0, the transversal signal still follows a Gaussian to good accuracy. The mean-field contribution changes the standard deviations; quite unexpectedly the mean-fields appear to counteract the static noise partly reducing the standard deviation. It turns out that this effect is well captured by the fit as can be seen in Fig. 16 with This quantifies the contribution of the mean-fields to the transversal spin dynamics. The behavior of the longitudinal autocorrelations is more complex; Fig. 17 shows that the decay is weaker than exponentially as we stated before. We refer to Sect. IV E for a detailed examination of the behavior.

D. Transition from weak to strong external magnetic field
If we consider the case of a surface perpendicular to the external magnetic field, i.e., ϑ = 0, we can compare the results from Sect. IV B for arbitrarily strong magnetic fields to the results from the previous Sect. IV C based on the RWA. This allows us to study how well the RWA reproduces the exact result. In particular, we can determine above which magnetic fields the RWA is reliable and to which extent. First, we compute the results of the self-consistency problem in the lab frame (exact spinDMFT) (79) and in the Larmor rotating frame using RWA (RWA spinDMFT) (100). Then, we transform the exact spinDMFT results to the Larmor rotating frame via so that a quantitative comparison is possible. justified if the energy scale of the magnetic field is larger than the energy scales of any other interaction in the system, including the static noise.

E. Long-time behavior
Now, we come back to the long-time behavior of the longitudinal autocorrelation. This is an interesting issue because various ideas exist on the origin of the rather slow decay and its functional form [71]. Fig. 17 shows that the autocorrelations do not decay in a Gaussian fashion at all. Such decay would have led to a negative curvature downwards. Instead we discern a positive curvature upwards which implies that the decay is even slower than exponential. The question arises which functional dependencies describe this decay. Considering this, our most successful fitting attempt is a power law according to 4g zz fit: with parameters B and m. We checked also stretched exponentials since these were suggested in Ref. [71]. But we did not achieve satisfactory fits for an ansatz according to 4g zz fit: with parameters A, α, and the exponent ν. The longitudinal results including the power law fits can be seen in Fig. 24. Note that much longer times are not easily accessible for two reasons. First, the numerical effort increases as t 2 for the Monte-Carlo simulation and as t 3 for the diagonalization. Second, the smaller the autocorrelation is the more difficult it becomes to determine it with good relative accuracy in view of the statistical way of computing it. This is also the reason why we cannot go to longer times for C ≈ 0. The exponent m displays a pronounced dependence on the relative noise strength C as depicted in Fig. 25. The dependence m(C) can be described heuristically by This fit works surprisingly well in spite of the divergence for C → 0. Limited by the numerical accuracy, we can hardly say if we actually obtain this divergence or if it can be truncated, e.g., by replacing C → C +C 0 in (109). This issue is associated to the question if the power law behavior solely results from the presence of the noise or if it is a valid feature of spinDMFT. All in all, these results provide evidence that the longitudinal autocorrelations are relatively long-lived. Clearly, their long-time behavior poses an interesting question which calls for further research, both numerical and analytical.

V. CONCLUSIONS
In this paper, we introduced and justified a mean-field theory designed to capture the spin dynamics in disordered dense spin systems. The key idea is that it is not sufficient to introduce a static mean-field but that the mean-field is dynamic itself so that we call it "spin dynamic mean-field theory" (spinDMFT). As usual, this approach becomes exact if each site has an infinite number of interaction partners, i.e., the coordination number becomes infinitely large. Historically, the same limit led to the introduction of the fermionic dynamic mean-field theory [44,72]. For spins, we established that the important correlations are the autocorrelations and that these define the dynamic mean-fields to which each spin is subjected. These mean-fields are normally distributed and the dynamic variances of these normal distributions are given by the  autocorrelations. This constitutes the self-consistency problem which has to be solved for spinDMFT. We showed how this can be done stochastically.
If the effective single-site problem is linear in the spin operators it does not matter whether we consider classical spins averaged over all directions or a quantum spin given that the average length is scaled to be the same. In this sense, the quantum spin system and the classical spin system have the same spinDMFT. For spin-1/2 this has to be the case since locally only linear spin operators can appear. For larger spins, however, higher powers may arise such as anisotropies of various kinds. Then, the quantum spinDMFT and the classical spinDMFT are different.
We gauged the advocated spinDMFT against numerical results for isotropic spin systems obtained by other methods, namely exact diagonalization, iterated equations of motion, and Chebyshev expansion. This can only be done for rather small spin clusters in low dimensions so that the reproduction of the results by spinDMFT is particularly challenging. Nevertheless, encouraging agreement could be established. Subsequently, we applied spinDMFT to a twodimensional ensemble of spins with dipolar interactions including local static noise, i.e., fluctuations of magnetic fields. This underlines that spinDMFT is capable of dealing with anisotropic interactions of long range as well.
We studied the case of an arbitrary external magnetic field perpendicular to the plane of spins and the RWA for a large tilted magnetic field. For zero tilt, we compared both approaches quantitatively. This allowed us to show quantitatively to which extent the RWA is justified and above which magnetic field it yields reliable results. We showed that the transversal autocorrelations behave essentially like Gaussians in time. The longitudinal autocorrelations, however, display a more complex behavior with a rather slow decay towards long times. Evidence for power law behavior is found. This certainly calls for further investigations.
We are confident that spinDMFT can be applied successfully to many more physical systems aside from the ones that we mentioned so far. Ample applications can be found for nuclear magnetic resonance (NMR), electron spin resonance (ESR), quantum information storage and processing based on spins in solid state systems, in particular in nanostructures, and for all phenomena related to spin diffusion in such systems. Conceptually, further issues to be addressed are the treatment of explicitly time-dependent Hamiltonians and spatially inhomogeneous solutions with distributions of mean-fields which vary along the samples. Both extensions are of greatest interest, for example in the coherent control of spin degrees of freedom. A third fascinating issue consists in an extension of spinDMFT to finite temperatures. So far, we derived the approach for disordered ensembles. But in view of related developments for spin glasses [47,48] and the general analogy between real and imaginary timess a dynamic mean-field theory for spins at finite temperatures should also exist. Its development would enhance the applicability of spinDMFT even further. nian with arbitrary couplings J αβ = J βα allowing for spin anisotropy. The factor z − 1 2 is denoted separately to explicitly keep track of the scaling with z. It must be chosen in this way to keep the energy scale of the dynamics, J 2 in Eq. (4), constant for z → ∞. The argument runs as follows. First, we use the Heisenberg equations of motion to set up a system of differential equations for the temporal evolution of the correlation functions. Second, we postulate the scaling of the correlations. Third, we show that the postulated scaling is consistent with the initial conditions and with the differential equations, i.e., the scaling is fulfilled by the initial conditions at t = 0 and by the differential equations. As starting point, we calculate the time derivative of the general pair correlation function, where the Levi-Civita tensor occurs due to the spin algebra. The complexity of the expectation value is increased by the Liouville operator L which consists in the commutation with the Hamiltonian. This leads to an an additional time-dependent spin operator in the expectation value.
The same behavior is observed for higher derivatives: the application of L results in expectation values with incremented or decremented number of spin operators by one. The decrement occurs if the additionally generated spin operator hits an already present spin operator at the same site and of the same component due to (S α ) 2 = 1/4. The quickly growing number of products of spin operators makes the book keeping tedious. A solution consists in considering the correlation function of a general cluster at time t with a single spin at time 0 and site 0 Here, C denotes an arbitrary product of spin operators at different sites where c is a set of lattice sites r ∈ {0, ..., N } and α is a set of components α r ∈ {x, y, z}. The time derivative of such cluster correlations (A3) reads where the sum runs over all clusters C that can be reached from C by one application of the Liouville operator. All factors that are independent of the coordination number, e.g., the couplings J αβ and factors 1/4 r = 0 from products of the same spin operators, are collected in J(C, C ). Hence, this generalized coupling J(C, C ) does not contribute to the scaling which we want to determine. By (A5) we have formally defined the system of differential equations describing the dynamics of all cluster correlation functions. The starting conditions read as It is challenging to keep track of the more and more complex clusters. To this end, we define a measure of the cluster size or complexity and link it to the scaling. Since S = 1 2 , at each site r of a cluster the active operator can be the identity or one of the three spin components. All products with more factors can be reduced to this case. We call a site "occupied" if one of the three spin components is present at this site. Otherwise, we call it "unoccupied". We introduce the measure κ(c) which consists of two components κ(c) := κ 1 (c) + κ 2 (c). (A7) The first term κ 1 (c) measures the overall size of the cluster c. Consider a covering of the minimum number of bonds needed to link all sites in c and the origin 0, see the green bold links in Fig. 26. We stress that this covering is unique, i.e., there is only one such covering due to the properties of the Bethe lattice where two sites are linked by one specific path. There are no loops except self-retracing paths. The second term κ 2 (c) counts the number of empty, unoccupied sites which are touched by this covering. We show below that κ(c) is a lower bound for the minimum number n(C) of applications of the Liouville operator L needed to generate C(c, α) from an initial cluster C 0 ({0}, α 0 ) with a single spin operator S α0 0 at the origin. For the sake of completeness, we mention that in the exceptional cases n = κ(c) + 2 = 2 holds for reaching S ρ 0 from S γ 0 . From the assertion n(C) κ(c) and the starting conditions (A6) we deduce This motivates our central assertion that the cluster correlation functions are scaling with z according to Obviously, this agrees with the starting conditions while all other cluster correlation functions vanish at t = 0.
To validate the claim (A10) for arbitrary times, we show that it is consistent with the equations of motion (A5). A single application of the Liouville operator L to a cluster C generates a sum of multiple clusters C . Since we consider nearest-neighbor interaction, each of these clusters C differs from C only at one link (i, j) where i and j are adjacent. Therefore, it is sufficient to study the possible processes on this link where we denote the subcluster of C or C on this link by C (i,j) and C (i,j) , respectively.
can be categorized in three different types. The relevant commutators read as Interestingly, the last commutator yields zero for S = 1 2 , so that this process does not contribute. The remaining two processes (a) and (b) are analyzed further. Since κ(c) depends on the covering of c, we distinguish whether the considered link (i, j) is part of this covering or not. This leads to two subcases for processes (a) and (b) as shown in Fig. 28. Figure 28. Relevant link processes for a single application of L including the effect on the covering which is shown in green.
As an example, we derive that the first process (a.i) preserves the asserted scaling. The extension of the covering by one link leads to while the number of unoccupied sites remains the same In order to assess the complete scaling one has to count how often a process can occur. We call this its multiplicity m(c, z). The multiplicity of process (a.i) can be bounded from above by counting at maximum z neighbors of each occupied site in c. There are κ 1 (c)−κ 2 (c)+1 of such sites so that holds. The indicators κ i depend on the size of the cluster, but not on the coordination number so that the scaling resulting from the differential equation (A5) is where g γ (C (a.i) , t) denotes the sum of the contributions of the clusters C to the correlation of cluster C via the link process (a.i). Clearly, the asserted scaling of g γ (C, t) is confirmed. Note that the first factor z − 1 2 in (A18a) stems from the overall scaling on the right hand side of Eq. (A5). In conclusion, we showed that the first process is consistent with the claimed scaling (A10). We do not repeat the line of argument for the processes (a.ii), (b.i), and (b.ii). Their effects on the scaling are Table III. Induced scaling of the contributions of the link process to the differential equation for C. The number −1, 0, +1 indicate the increments in κ1 and κ2. The entries of the column of the multiplicity provides the scaling factors in m(c, z). The last column shows the resulting scaling of the contribution of C to C.
summarized in Tab. III. While the processes (a.ii) and (b.i) yield the same scaling as (a.i), the last process (b.ii) is even suppressed by an additional factor z − 1 2 . Hence, we have derived that none of the processes violates the asserted scaling. Since this scaling holds initially, see Eq. (A11), we deduce by continuous induction via the differential equations (A5) that the scaling (A10) holds at all times on the Bethe lattice at infinite temperature. As mentioned in the beginning, we assume that this behavior is generic, i.e., that it applies to general lattices and clusters. This justifies the application of spinDMFT for dense systems by which we mean systems with large effective coordination numbers. A direct corollary applies to the two-time pair autocorrelations S α i (t 1 )S β j (t 2 ) . In this particular case, the cluster reads as so that only the site i belongs to c and Using the above multiplicities in combination with the scaling of the spin-spin correlations (A22) we obtain that the pair correlations scale like so that they are suppressed by at least the factor z −1 . Hence, no correlations between the local-environment fields need to be accounted for. By self-consistency, this extends to the second moments of the local mean-fields.
In a similar fashion, we show that the spin dynamics of S i is uncorrelated to its corresponding local-environment field V i : as well as to any other local-environment field V j (j = i): which both tend to 0 for z → ∞. For the last conclusion, one has to distinguish two cases again, see Fig. 30. First, site i can be nearer to j than to k. Second, site i can be nearer to k than to j. The latter case is dominant with a multiplicity of m ∝ 1 and a scaling of ∝ z − j−1 +1 . Together with the prefactor 1/ √ z this yields the provided scaling.   pectation values which are calculated by path integrals averaged over the distribution of the mean-field time series V. Numerically, we estimate the path integrals using a Monte-Carlo method. The autocorrelations are computed for M time series V and averaged which converges to g αβ (t 1 , t 2 ) for M → ∞. Since the time series are drawn independent of each other the variance of g αβ M (t 1 , t 2 ) is given by the variance of a single time series divided by M where σ 2 (. . . ) denotes the variance of the quantity in the bracket.
The standard deviation σ of a single time series depends on many parameters and cannot be calculated analytically in a simple way. But it is clear that its value is bounded by the maximum value of the autocorrelations, i.e., by 1 4 . Figs. 31, 32, and 33 show generic time dependencies of σ computed for different physical situations. Interestingly, we find that all of them converge to the value 1/(4 √ 3) for t → ∞. The rapidity of this convergence depends on the actual decay of the corresponding autocorrelation. This behavior can be understood by the following argument. The autocorrelations vanish for large t. Hence, the variance equals the quadratic mean lim t→∞ σ 2 g αβ (t, 0) = lim In a next step, we consider the vector-valued signal where R V (t, 0) denotes the orthogonal rotation matrix which describes the rotation of the initial spin vector due to its temporal evolution subjected to the fluctuating time series V (t). Obviously, the square of g β V (t, 0) is constant before and thus also after averaging.
Next, it is plausible and in accord with all previous results that the decoherence is sufficiently strong to have the spin vector R V (t, 0) S(0) lose all information about its initial direction. If it points in z direction at t = 0 it will point into any direction of the unit sphere after sufficiently long time. Therefore, each component α of the rotated spin vector has the same variance for t → ∞ and contributes equally to the squared vector in (B6). This allows us to conclude in perfect agreement with the numerical findings. This enables us to derive the reliable estimate for the statistical error from averaging over M times series. It becomes even exact for t → ∞.

Discretization Error
The goal of this appendix is to estimate the error resulting from the discretization of time, i.e., the error resulting from working with a finite time step δt > 0 instead of taking δt → 0. Certainly, this error will depend on details of the model. We discuss it for both the isotropic Heisenberg model and the dipole model without and with RWA.
We consider the discretization error of the diagonal autocorrelations α = β up to some maximum time t max which we determine such that holds, i.e., the autocorrelation has decreased to 1/25 of its initial value. For simplicity, we focus on the autocorrelation that decays slowest in each physical scenario, since it is most susceptible to the discretization error which accumulates in the course of its temporal evolution. Mostly this is the longitudinal autocorrelation.
To be specific, we compute g αα (t) on the time interval t ∈ [0, t max ] for various step widths The results for g αα (t) for ν < ν max are compared to the 'best' solution for ν max , i.e., the solution from the finest discretization which is used as reference. We define the discretization error where L(ν) = t max /δt(ν) is the number of time steps. We stress that the input data for ∆q(ν) are not exact,  so large that the deviations are of order 2/4 because the fluctuating diagonal autocorrelations g αα are bounded by 1/4. Thus, the relevant regime is at intermediate time steps between these two plateaus. Here the effect of the discretization can be discerned. The double logarithmic plots are consistent with the conclusion ∆q ∝ δt 2 as can be read off by comparing to the straight lines resulting from the quadratic power law. This can be easily understood by the approximation we have to use for the unitary time evolution operators which propagate the system from t to t + δt. The employed second-order CFET and the trapezoidal rule entail a discretization error scaling like δt 3 . But since this error accumulates over time one has to multiply this scaling by the number of step from t = 0 to t = t max so that we obtain ∆q ∝ Lδt 3 ∝ δt 2 .
This explains the observed quadratic scaling of the discretization error with the time step size δt. Aside from the above discussed scaling further conclusions on the influence of the discretization can be drawn. Inspecting Figs. 34 and 35 shows that a finite magnetic field increases ∆q. The reason is that the Larmor precessions need to be resolved. If the step size δt is too long, approaching the Larmor period, sizable discretization errors occur. As a rule of thumb, δt should be at maximum a tenth of the Larmor precession period Since the time discretization error accumulates over time longer lasting correlations imply larger discretization errors. This can be seen for instance in Fig. 36 where the dipole model is studied in RWA. The error ∆q increases with the noise strength because the decay of the longitudinal autocorrelation is slowed down upon increasing C.
Finally, we point out that a small discretization error is desirable. But in view of the efficiency of the total algorithm it does not pay to reduce the discretization error below the statistical error. Hence, the parameters should be set such that ∆q ≈ σ(g) holds.

Termination condition for the iteration
We determine the solution of the self-consistency conditions iteratively. If the algorithm is stable, the autocorrelation functions g αβ (i) (t) of the iteration i converge to the exact results for i → ∞. In practice, it is necessary to define a termination condition to decide when the iterations can be stopped. For this we use This quantity measures the difference between the results of iteration j+1 and j, where j ≥ 1 . The execution of the code is stopped if ∆I αβ (i) falls below a certain threshold.
Since the iteration error itself is limited by the statistical accuracy this threshold cannot be chosen smaller than the statistical error. It needs to be set above the statistical standard deviation to achieve reliable termination. In our numerics, it turned out that ∆I threshold = 2σ g αβ M (t 1 , t 2 ) = is a reasonable choice. It avoids unnecessary iterations while it is strict enough to yield sufficient convergence. Figure 37 depicts the iteration error of g zz for the isotropic Heisenberg model without magnetic field. The iterations started from two different initial diagonal autocorrelations g αα while the cross autocorrelations g αβ with α = β are set to zero. The iteration error decreases very fast in the beginning. At about the fourth iteration it reaches the magnitude of the threshold. Beyond the fourth iteration, statistical fluctuations stemming from the averaging over M time series occur and dominate the iteration error ∆I xx . For optimum computational efficiency, the code should terminate before the statistical fluctuations take over. We emphasize that the "converged" results obtained in this way are independent of the initially chosen autocorrelations. The deviation between the iterated results from the exponential and the Gaussian initial autocorrelation is of the same magnitude as the error threshold (B15). For this reason, we consider the chosen iterative algorithm robust enough to determine physically meaningful solutions.

Time-Translation Invariance
For both the isotropic and the dipolar spin systems we considered a time-independent Hamiltonian so that the systems are invariant under time translation and so are all two-time spin autocorrelations S α (t 1 )S β (t 2 ) = S α (t 1 − t 2 )S β (0) .
Applying spinDMFT, the situation is slightly more complicated. The autocorrelations with respect to a single time series V i are not time-translation invariant because the mean-field is dynamic, i.e., it depends on time.
But the physically meaningful expectation values are obtained by averaging over the distribution of possible mean-fields. In the path integral the mean-field V i can be shifted in time so that the time evolution operators and the distribution functional p[ V i ] are shifted. Timetranslation invariance is ensured if and only if the distribution remains unchanged by this shift. Specifically, this is equivalent to the condition for the second moments of the mean-fields, i.e., they are time-translation invariant themselves. As stated in the main text, the time invariance of the moments and hence of the distribution implies time translation invariance of the spin-spin autocorrelations. So timetranslational self-consistent solutions are possible, but it cannot be excluded that solutions exist which are not time-translationally invariant. Here, we do not consider them because they cannot occur in a quantum system at infinite temperature although time-crystalline behavior cannot be excluded generally for specific situations. In practice, the averaging based on Monte-Carlo sampling allows for small unphysical violations of time translational invariance due to the finite statistical error: the self-consistently computed moments and autocorrelations may not be exactly time-translation invariant, even though the initially inserted moments have this property.
The reduced set of data together with the assumption of time translation invariance provides all required information. Due to the time-translational invariance the covariance matrix has constant matrix elements on all diagonals This procedure has two advantages: (i) time-translation invariance is built-in by construction; (ii) the autocorrelations need only be computed at the time differences ∆t. Hence, the effort of steps 4) and 5) of the numerical procedure is reduced from O(L 2 ) to O(L) where L is the number of time steps.

Definiteness of the Covariance Matrix
In the self-consistency conditions the two-time correlations are taken as matrix elements of a covariance matrix. The physical justification is given in Sect. II. In addition, we discuss whether this identification is mathematically possible. To this end, the quantum correlations must provide (i) real symmetric matrix elements and the resulting matrix must be (ii) positive semi-definite. In Sect. II D we already showed property (i) at infinite temperature. Property (ii) is equivalent to for arbitrary real coefficients λ α t and t ∈ I where I is an arbitrary set of time instants. By shifting the sums into the expectation value we find α t∈I Since B is a Hermitian operator, its square is a nonnegative operator and thus any of its expectation values are non-negative. This proves the required positive semidefiniteness and applies to the expectation values of the path integral (15) and extends to the averages (B2). There is, however, a subtlety in the implementation described in the previous section. We guaranteed time translation invariance by using For M = ∞ this is relation holds exactly true. But the statistical fluctuations at finite M imply that the covariance matrix computed by Eq. (B20) need not be nonnegative. If we refrained from setting all matrix elements on a diagonal constant, but used the two-time autocorrelations, L would be non-negative. Enforcing time translation invariance breaks the positive semi-definiteness to the extent of the statistical standard deviation. Since we keep this deviation low the computed covariance matrices are very close to being non-negative. Due to the Monte Carlo approach to the averaging all numerical results are subjected to the statistical error so that the statistical inaccuracies of the approximate covariance matrix are no source of additional deviations.
Practically, we diagonalize the approximate covariance matrix and set the eigenvalues which are negative to zero for the next iteration step. We stress that the modulus of these negative eigenvalues is of the order of the statistical error, i.e., small in a systematically controlled way. Therefore, the employed procedure provides consistent physical solutions for the dynamic autocorrelations within the discussed errors.