Classical spin simulations with a quantum two-spin correction

TL;DR

This paper introduces a quantum correction to classical spin simulations to better model nuclear spin dynamics in solids, especially when the number of interacting neighbors is small, but faces challenges at longer times.

Contribution

The authors propose adding a two-spin quantum correlation term to classical simulations, improving accuracy for small neighbor interactions in spin chains.

Findings

Improved short-time accuracy in spin dynamics simulations.

Classical plus quantum correction method shows limitations at long times.

Method tested on a spin-1/2 chain with promising results at short times.

Abstract

Classical simulations of high-temperature nuclear spin dynamics in solids are known to accurately predict relaxation for spin 1/2 lattices with a large number of interacting neighbors. Once the number of interacting neighbors becomes four or smaller, classical simulations lead to noticeable discrepancies. Here we attempt to improve the performance of the classical simulations by adding a term representing two-spin quantum correlations. The method is tested for a spin-1/2 chain. It exhibits good performance at shorter times, but, at longer times, it is hampered by a singular behavior of the resulting equations of motion.

Figures

Equations

$$\hat{H}$$

$$C(t)=\frac{\mathrm{Tr}\left({\hat{S}}^x(t){\hat{S}}^x\right)}{\mathrm{Tr}\left({\hat{S}}^{x2}\right)},$$

$${\hat{S}}^x(t)=\exp (i\hat{H}t/\hslash ){\hat{S}}^x\exp (-i\hat{H}t/\hslash )$$

$${\hat{\rho }}_n=\mathcal{O}\left(Z^0\right),{\hat{\rho }}_{m,n}^{\mathrm{corr}}=\mathcal{O}\left(1/Z\right),{\hat{\rho }}_{m,n,r}^{\mathrm{corr}}=\mathcal{O}\left(1/Z^2\right),$$

$$S_n^{\alpha }=⟨{\hat{S}}_n^{\alpha }⟩$$

$$⟨\delta {\hat{S}}_m^{\alpha }\delta {\hat{S}}_n^{\beta }\delta {\hat{S}}_r^{\gamma }⟩$$

$$i{\partial }_t{\hat{\rho }}_n$$

$$i{\partial }_t{\hat{\rho }}_{n,m}^{\mathrm{corr}}$$

$${\hat{\mathcal{L}}}_{n,m}\hat{A}\equiv \left[\sum _{\alpha =x,y,z}{J}_{m,n}^{\alpha }{\hat{S}}_m^{\alpha }{\hat{S}}_n^{\alpha },\hat{A}\right],$$

$${\partial }_t{S}_m^{\alpha }$$

$${\partial }_t{M}_{mn}^{\alpha \beta }$$

$${\hat{\rho }}_m(t=0)=|{\mathbf{s}}_m⟩⟨{\mathbf{s}}_m|=\frac{{\hat{1}}_m+2{\mathbf{s}}_m.{\hat{\mathbf{S}}}_m}{2}.$$

$${\hat{1}}_m=2\oint \frac{d{\mathbf{s}}_m}{4\pi }{\hat{\rho }}_m(t=0),$$

$${\hat{S}}_m^x={\hat{S}}_m^x(t=0)=3\oint \frac{d{\mathbf{s}}_m}{4\pi }{\hat{\rho }}_m(t=0)s_m^x.$$

$${\mathbf{S}}_m(t=0)$$

$$\mathrm{Tr}[{\hat{S}}^x(t){\hat{S}}^x]$$

$$M_{mn}^{\alpha \beta }(0)=0.$$

$${\hat{S}}_{ϵ}^x=\frac{{\hat{S}}^x-ϵ{\hat{S}}^z}{\sqrt{1+{ϵ}^2}},{\hat{S}}_{ϵ}^z=\frac{{\hat{S}}^z+ϵ{\hat{S}}^x}{\sqrt{1+{ϵ}^2}}$$

$${\hat{\rho }}_m^{\pm }(t=0)=|{\mathbf{s}}_m^{\pm }{⟩}_m⟨{\mathbf{s}}_m^{\pm }|=\frac{{\hat{1}}_m+2{\mathbf{s}}_m^{\pm }.{\hat{\mathbf{S}}}_m}{2}.$$

$$\mathrm{Tr}[{\hat{S}}^x(t){\hat{S}}^x]$$

$${\partial }_t{S}_m^{\alpha }$$

$${\partial }_t{M}_{mn}^{\alpha \beta }$$

$${\partial }_t(M_{mn}^{\alpha \beta }+S_m^{\alpha }{S}_n^{\beta })$$

$$({\partial }_t^2+{\omega }_x^2)(S_m^x+S_n^x)=0$$

$$\frac{\mathrm{Tr}\left({\hat{S}}^x(t){\hat{S}}^x\right)}{\mathrm{Tr}\left({\hat{S}}^{x2}\right)}$$

$$\mathrm{Tr}[{\hat{S}}^x(t){\hat{S}}^x{]}^{(2)}$$

$${S_m^y}^{(1)}$$

$${S_m^z}^{(1)}$$

$${S_m^x}^{(2)}$$

$$M_{m,m\pm 1}^{yz(1)}$$

$$\frac{\mathrm{Tr}\left({\hat{S}}^x(t){\hat{S}}^x\right)}{\mathrm{Tr}\left({\hat{S}}^{x2}\right)}$$

Full text

1 Introduction

Nuclear spin-spin relaxation in solids has been a subject of active theoretical researchBloch-46-1 ; VanVleck-48 ; Lowe-57 ; Abragam-61 ; Tjon-66 ; Gade-66 ; Borckmans-68 ; Jensen-73 ; Parker-73 ; Jensen-73A ; Engelsberg-74 ; Engelsberg-75 ; Becker-76 ; Jensen-76 ; Lundin-77 ; Shakhmuratov-91 ; Lundin-92 ; Tang-92 ; Jensen-95 ; Sodickson-95 ; Fine-97 ; Fine-00 ; Dobrovitski-06 ; Zhang-07 ; Elsayed-13-thesis ; Elsayed-13 ; Elsayed-15 since the discovery of nuclear magnetic resonance (NMR) Purcell-46 ; Bloch-46-2 . Nuclear spin dynamics is usually considered in the limit of infinite temperature as the nuclear gyromagnetic ratios and, hence, the interaction energies are very small. Despite this simplification, the dynamics is still nontrivial. Making controllable first-principles predictions of NMR relaxation has remained an elusive goal for the practitioners in the field.

Classical spin simulations were shown to be rather accurate in predicting the NMR spin-spin relaxation for the lattices of quantum spins 1/2 with a large number of interacting neighborsJensen-73 ; Jensen-76 ; Lundin-77 ; Tang-92 ; Elsayed-13-thesis ; Elsayed-15 . However, as the number of the interacting neighbors decreases to four or smaller, classical simulations lead to noticeable discrepancies. In a broader context, similarity and contrast between quantum and classical spin dynamics were also investigated. On the one hand, there exists significant experimental and numerical evidence supported by theoretical argumentsBorckmans-68 ; Engelsberg-74 ; Fabricius-97 ; Fine-03 ; Fine-04 ; Fine-05 ; Morgan-08 ; Sorte-11 ; Meier-12 ; Steinigeweg-12 that the asymptotic high-temperature relaxation in ergodic classical and quantum spin systems has the functional form either exp$(-\gamma t)$ or exp$(-\gamma t)$ cos$(\omega t+\phi)$, where $\gamma$, $\omega$ and $\phi$ are some constants. On the other hand, the classical spin dynamics is known to be chaotic deWijn-12 ; deWijn-13 ; Fine-14 ; Elsayed-15A , while the dynamics of spin-1/2 lattices was conjectured in Refs.Fine-14 ; Elsayed-15A to be non-chaotic in the sense of not exhibiting exponential sensitivity to small perturbations.

In this work, we attempt to improve classical simulations of quantum spin dynamics by introducing corrections representing purely quantum correlations between each pair of interacting spins. The correction is introduced on the basis of an expansion in terms of the inverse of the lattice coordination number NS10 ; QKNS14 . This expansion was already used to describe the quantum dynamics of other lattice systems such as the Bose- and Fermi-Hubbard models NS10 ; QKNS14 ; NQS14 ; KNQS14 and the Heisenberg spin model NTZ17 ; NQS14 . An alternative approach to combining classical and quantum simulations was recently introduced by two of us in Ref. greg .

In section 2, we give general formulation of the NMR spin-spin relaxation problem. In section 3, we present the method of large coordination number expansion, from which we derive the dynamical equations for the spin average and correlations. In section 4, we present the results of applying our method to a chain of spins 1/2. Section 5 contains the conclusions.

2 General formulation

We consider a $D$-dimensional hypercubic lattice of spins 1/2 described by the nearest-neighbor interaction Hamiltonian

[TABLE]

where $\hat{S}_{n}^{i}$ ($i=x,y,z$) are the spin projection operators for site $n$, $J_{mn}^{i}$ are the coupling constants for the $i$th projections of spins $m$ and $n$.

The quantity of interest in the context of NMR is the infinite-temperature autocorrelation function defined as Lowe-57 ; Abragam-61 :

[TABLE]

where $\hat{S}^{x}\equiv\sum_{m}\hat{S}^{x}_{m}$ is the $x$-component of the total spin polarization. Its time dependence is defined in the Heisenberg representation as

[TABLE]

Function $C(t)$ is proportional to the signal of NMR free induction decay.

A previous work Elsayed-15 obtained $C(t)$ by classical simulations for one-, two- and three-dimensional lattices. For two- and three-dimensional lattices the classical results exhibited reasonable agreement either with the direct quantum calculation or with experiment. However, for one-dimensional chain, the agreement was not very good. Here we aim at achieving an improvement in the latter case by taking into account the next order quantum corrections describing the two-spin correlations.

In order to derive the quantum correction, we perform an expansion in terms of the inverse powers of the lattice coordination number $Z=2D$. When applied up to the order $1/Z^{2}$, this method showed a reasonable convergence in the case of the Bose-Hubbard lattice model in one and two dimensionsKNQS14 . It has also been used successfully for the quantum Ising model NQS14 ; NTZ17 . In the present work, the resulting approximation is to be tested for the worst possible case, namely, for a one-dimensional chain (Z=2).

3 Large coordination number expansion

3.1 Overview of the method

For $Z\gg 1$, the model dynamics is to be described using the method developed in NS10 ; QKNS14 , which we now introduce.

The time evolution of the density matrix $\hat{\rho}$ of the whole lattice is given by the von Neumann-Liouville equation $i\hbar\partial_{t}\hat{\rho}=\left[\hat{H},\hat{\rho}\right]$. This density matrix is usually too complex to be analyzed. Instead, the set of reduced density matrices is introduced, $\hat{\rho}_{\cal S}=\,{\rm Tr}_{\not{\cal S}}\hat{\rho}$, which results from tracing out the Hilbert spaces of all sites except a few: ${\cal S}=\{n_{1},n_{2},\dots,n_{i}\}$. If we keep only one site $n$, then the reduced density matrix is a linear operator $\hat{\rho}_{n}$ acting on the Hilbert space of one lattice site $n$; if we keep two sites $n,m$, then $\hat{\rho}_{n,m}$ acts on the Hilbert space of two sites etc. The decomposition $\hat{\rho}_{m,n}=\hat{\rho}_{m,n}^{\rm corr}+\hat{\rho}_{m}\hat{\rho}_{n}\,,$ and $\hat{\rho}_{m,n,r}=\hat{\rho}_{m,n,r}^{\rm corr}+\hat{\rho}_{m,n}^{\rm corr}\hat{\rho}_{r}+\hat{\rho}_{m,r}^{\rm corr}\hat{\rho}_{n}+\hat{\rho}_{n,r}^{\rm corr}\hat{\rho}_{m}+\hat{\rho}_{m}\hat{\rho}_{n}\hat{\rho}_{r}$, etc. allows us to derive an exact hierarchy of interlinked equations for these operators NS10 ; QKNS14 . This hierarchy is the counterpart of the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) chain but with correlations between sites and not between particles.

In order to treat the time dynamics, it was demonstrated quite generally in NS10 ; QKNS14 that, if the initial state of any quantum lattice system is separable without initial correlations between sites, then the correlations at a later time satisfy – at least for a finite period of time – the following hierarchy scaling $\hat{\rho}^{\rm corr}_{\cal S}\sim 1/Z^{|{\cal{S}}|-1}$, which means that the higher-order correlations are suppressed as an inverse power of the coordination number $Z$. More explicitly,

[TABLE]

and so on. Using the spin representation, this hierarchy can be rewritten as

[TABLE]

where, for an operator $\hat{A}$, $\delta\hat{A}=\hat{A}-\langle\hat{A}\rangle$ and $\langle\hat{A}\rangle=\,{\rm Tr}(\hat{A}\hat{\rho})$.

Many methods of quantum field theory use similar expansion techniques such as the $1/{\cal N}$ expansion, where ${\cal N}$ is the number of field components, or the $1/S$ expansion, where $S$ is the quantum spin. The open question is whether this kind of expansion converges accurately towards the exact solution. In our case, we will expand only up to the first order and test the resulting approximation.

In NS10 ; QKNS14 , an exact set of hierarchy equations has been derived for these density matrices. Up to the first order, these are written as:

[TABLE]

where we define the Liouville operators in terms of the commutators:

[TABLE]

and $\,\widehat{\cal L}^{S}_{m,n}=\,\widehat{\cal L}_{m,n}+\,\widehat{\cal L}_{n,m}$. Note that the trace in the first line of Eq.(7) is carried out on the Hilbert space associated to the site $n$ without a summation over the index $n$. This set forms the basis of the $1/Z$-expansion up to the first order. The term $\hat{\rho}^{\rm corr}_{n,m}$ in Eq.(7) describes quantum pair correlations beyond those determined by single-spin density matrices $\hat{\rho}_{m}$ NQS14 .

3.2 General set of equations including the two-spin quantum correlations

The averages defined in Eqs.(3.1) are used into the Eqs(6-7) in order to determine their dynamics evolution. Assuming $J_{mn}^{\alpha}=J_{nm}^{\alpha}$, we arrive at the following set of expectation values for the single spin and pair operators:

[TABLE]

[TABLE]

where $\epsilon^{\alpha\gamma\delta}$ is the anti-symmetric Levi-Civita tensor.

Equations (9) are identical to the equation of motion for the classical spins, when the correlation terms $M^{\beta\gamma}_{mr}$ are removed. Together, the self-consistent system (9,10) describes the dynamics of classical spin and the quantum correlations between distant spin pairs.

We note that Eqs.(9,10) preserve the conservation laws such as the total energy average, the total spin projection average on the z-axis and the average of the total spin squared (if conserved by ${\hat{H}}$). This is achieved by keeping in Eq.(10) the terms that do not have summation over index $r$ and hence scale as $1/Z^{2}$. We also show in appendix in 6.1 that Eqs.(9,10) lead to the exact result for the case of two spins, which is an indicator of the quantitative promise of the simulation schemeFine-97 .

Equation (9) without correlation terms $M^{\beta\gamma}_{mr}$ preserves the length of the individual spin, follows classical Hamiltonian dynamics and, therefore, does not display any dynamical instability. Once the correlation terms are included, the individual spin lengths are no longer conserved, the dynamics loses the Hamiltonian character, and, as a result, dynamical instabilities eventually set in, accompanied by negative eigenvalues of the reduced density matrices for individual spins. For large $Z$, we expect that the correlation terms $M^{\beta\gamma}_{mr}$ remain small for an extended period of time, which, in turn, delays the onset of the above instabilities. We further expect that the above instabilities are delayed more, if the higher-order terms of the $1/Z$ expansion are included. This subject, however, is beyond the scope of the present work.

3.3 Ensembles of initial conditions

3.3.1 Random spin ensembles

We define quite generally any one-site reduced matrix density associated with the spin direction ${\bf s}_{m}$ of a quantum state:

[TABLE]

We set the normalization to $|{\bf s}_{m}|^{2}=1$ in order to have the density matrix representing a pure state. We then define the ensemble of density matrices by adopting the uniform probability distribution of all possible orientations of vector ${\bf s}_{m}$. This operation is accomplished using the identity:

[TABLE]

where we define the integration over every spin direction. Using these notations, the spin operators along the $x$ axis can be rewritten as:

[TABLE]

This definition implies the following initial expectation value for one realization of initial conditions:

[TABLE]

As a consequence, the autocorrelation function (2) can be rewritten as:

[TABLE]

where $L$ is the number of spins. For the correlation terms, we use the initial conditions

[TABLE]

Thus, the functions $S^{x}_{m}(t)$ become functionals of the set $\{{\bf s}_{m}\}$, which can be obtained by solving the system (9) and (10) with the initial conditions given by Eq.(14, 16). For each set of initial conditions ${\bf s}_{m}$, we determine $S^{x}_{m}(t)$ and then carry out the integral over all initial configurations. In comparison to the classical description, the inclusion of the correlation terms ensures the exact recovery of the short time dynamics, more precisely the second order moment term $C(t)=1-\frac{(J^{z}-J^{y})^{2}}{2}\frac{t^{2}}{2!}+...$ derived in appendix 6.2.

These considerations can be generalized for a normalization $|{\bf s}_{m}|$ chosen arbitrarily. When $|{\bf s}_{m}|$ is less than one, the reduced density matrix of a given spin becomes mixed but, when $|{\bf s}_{m}|$ is greater than one, it is difficult to interpret the resulting density matrix physically, because one of its eigenvalues is negative, while the other is greater than 1. Yet, such states may be used in the context of an ensemble average by a straightforward generalization of (15) with the help of a renormalization. Such an approach can be justified by noting that, when a classical-spin limit is taken, the simulated system is no longer a lattice of spins 1/2 — rather it can be thought of as a lattice of large quantum spins. Therefore, one should concentrate on the average spin polarizations rather than on the density matrices. We explore such a renormalization of $|{\bf s}_{m}|$ as follows.

The ensemble based on the “natural” normalization of initial conditions $|{\bf s}_{m}|=1$ imposes the individual quantum spin fluctuations $\langle\delta{\hat{\bf S}}_{m}^{2}\rangle=3/4-{\bf S}_{m}^{2}=1/2$ associated with a pure state. We also consider the value $|{\bf s}_{m}|=\sqrt{3}$, which implies no initial quantum fluctuations $\langle\delta{\hat{\bf S}}_{m}^{2}\rangle=0$, in a close correspondence with classical simulations. The normalization $|{\bf s}_{m}|=\sqrt{3}$ exactly reproduces the second moment already at the level of the classical equations as shown in a previous work Elsayed-15 , thereby leading to a qualitatively correct behavior of $C(t)$ for a spin-1/2 chain. It was also shown analytically in Ref.Lundin-77 , that, in the limit of the infinite number of interacting neighbours, correlation functions $C(t)$ computed classically and quantum-mechanically become identical — consequence of the fact that commutators for quantum spins and Poisson brackets for classical spins have essentially the same structuredeWijn-13 ; greg and, as a result, lead to the same expressions for the time derivatives of $C(t)$.

When correlation terms $M^{\beta\gamma}_{mr}$ are included, the second moment of $C(t)$ (second derivative at $t=0$) is always reproduced exactly irrespective of the initial normalization of $|{\bf s}_{m}|$. The idea behind the use of the normalization $|{\bf s}_{m}|=\sqrt{3}$ in the simulations including correlations $M^{\beta\gamma}_{mr}$ is to reproduce the second moment of $C(t)$ already without $M^{\beta\gamma}_{mr}$, so that the growth of $M^{\beta\gamma}_{mr}(t)$ is additionally delayed.

3.3.2 The z-basis ensemble

In addition to the random ensemble defined in the preceding subsection, we also consider an anisotropic ensemble with initial spin polarizations oriented approximately along the $z$-direction. Such an ensemble might be more adequate for Hamiltonians commuting with the $z$-projection of the total spin, especially when $J^{z}_{m,n}$ is the largest coupling constant. The initial orientations of ${\bf s}_{m}$ cannot simply be all along the $z$-direction, because, if they were, then each realization of such initial conditions would correspond to an unstable fixed point of the dynamics governed by Eqs.(9,10). For this reason, we rotate the z-axis about the y-axis by a vanishing angle $\epsilon\rightarrow 0$, so that the x-projections of spins acquire small numerically tractable values:

[TABLE]

We define the ensemble by the initial set of values ${\bf s}_{m}^{\pm}=\pm(\epsilon,0,1)/\sqrt{1+\epsilon^{2}}$, which leads to the following set of reduced density matrices:

[TABLE]

Using the commutation relation $[\hat{S}^{z},\hat{H}]=0$, we then express the autocorrelation as:

[TABLE]

4 Tests for one-dimensional chains

In our numerical tests, we restrict ourselves to a periodic spin-1/2 chain with $L=12$ sites. We use the values $J^{z}=0.82$, $J^{x}=J^{y}=-0.41$ representing the ratio of the coupling constants between nearest-neighbor spins typical of NMR settings Elsayed-15 ; VanVleck-48 ; Lowe-57 ; Abragam-61 . In such a case, the $z$-projection of the total spin commutes with the Hamiltonian. The reference plots representing fully quantum dynamics are obtained from an exact diagonalization.

Our simulations were subject to the dynamical instabilities mentioned earlier. We found that some initial conditions lead to stable solutions for longer times, while others create instabilities faster. In the results presented in Figs. 1-4, a tiny fraction of completely diverging solutions was not included in the averaging.

The test for the random spin ensemble with $|{\bf s}_{m}|=1$ is presented in Fig. 1. It shows a rather noticeable discrepancy with the reference plot. In comparison, the test for the random spin ensemble with $|{\bf s}_{m}|=\sqrt{3}$ shown in Fig. 2 exhibits a much better initial agreement. However, this initial performance was still not better than that of purely classical simulations reported in Ref.Elsayed-15 . At longer times, the individual solutions generated with both $|{\bf s}_{m}|=1$ and $|{\bf s}_{m}|=\sqrt{3}$ ensembles begin exhibiting singular behavior leading to an increasingly poor convergence of the statistical averaging procedure.

The test for the $z$-basis ensemble with $|{\bf s}^{\pm}_{m}|=1$ is shown in Fig. 3. Here the averaging was carried out over all possible $2^{L}=4048$ initial conditions. This discrete summation procedure appears to be more efficient than that for the isotropic random spin ensemble. Yet, for $t>2$ a noticeable discrepancy sets in. We also tested the $z$-basis ensemble with $|{\bf s}^{\pm}_{m}|=\sqrt{3/2}$, which would lead to the correct second moment, without the contribution from the correlation term. The resulting initial performance has improved and, in fact, become arguably better than that of the purely classical simulations of Ref.Elsayed-15 . Yet the agreement with the reference plot at longer times remains unsatisfactory.

5 Summary and outlook

We analyzed the performance of the large coordination number expansion method with the purpose of modeling NMR free induction decays in solids. The method exhibited promising analytical structure. However, despite certain improvements at the initial and intermediate times, the numerical tests indicate that the first-order approximation within this method applied to spin-1/2 chains is still not sufficient to outperform the existing classical calculations. A possible reason for this lies in the sensitivity of the method to the choice of the ensemble of initial conditions. To conclude, the application of the large coordination number expansion in the NMR context requires further investigations including higher-order terms, different ensembles of initial conditions, and, perhaps, different ways of including the quantum corrections to the classical equations of motion.

Acknowledgements

This work was supported by a grant of the Russian Science Foundation (Project No. 17-12-01587).

6 Appendix

6.1 Two spin cases

In the case of two spins $m$ and $n$, we set $Z=1$ and the system of Eqs.(9) and Eqs.(10) simplifies into:

[TABLE]

These two set of equations can be combined into

[TABLE]

The term inside the first time derivative has the same form as the right-hand-side of Eq.(6.1). Taking the second time derivative of (6.1), we eliminate this term to obtain for the $x$ component:

[TABLE]

where $\omega_{x}=\frac{|J_{mn}^{y}-J_{mn}^{z}|}{2}$ is frequency of oscillation with the solution $S^{x}(t)=S^{x}(0)\cos(\omega_{x}t)$. Therefore, up to the first order in the large coordination number expansion, we recover the exact purely oscillatory autocorrelation function $C(t)=\cos(\omega_{x}t)$.

6.2 Second moment of $C(t)$

The initial behavior of $C(t)$ can be determined by expanding it up to the second order in $t$, thereby obtaining the second moment as a coefficient in front of $t^{2}/2$:

[TABLE]

This exact result can be compared with the one obtained from our approximations. For this purpose, we solve Eqs.(9) through a perturbation expansion $S^{\alpha}_{m}(t)={S^{\alpha}_{m}}+{S^{\alpha}_{m}}^{(1)}t+{S^{\alpha}_{m}}^{(2)}t^{2}/2!+\dots$ and $M^{\alpha\beta}_{mn}(t)={M^{\alpha\beta}_{mn}}^{(1)}t+\dots$ and for the autocorrelation function:

[TABLE]

After the expansion, we obtain the relevant set of equations:

[TABLE]

Solving this set and inserting the result into (24), we recover, after setting $S^{i}_{m}=s^{i}_{m}/2$, the exact result (23).

On the other hand, if we neglect the quantum pair correlation terms $M^{yz(1)}_{m,m\pm 1}$ for the random spin ensemble, we obtain for any normalization $|{\bf s}_{m}|$:

[TABLE]

The random spin ensemble with the normalization $|{\bf{s}}_{m}|=\sqrt{3}$ leads to the correct second moment (23) without the quantum terms $M^{\alpha\beta}_{mn}$ but differs from the “natural” normalization $|{\bf{s}}_{m}|=1$ by a factor $\sqrt{3}$. In contrast, for the $z$-basis ensemble, we obtain instead

[TABLE]

Since the coupling values are chosen such that $J^{z}-J^{y}=3J^{z}/2$, the normalization is fixed to $|{{\bf s}}_{m}^{\pm}|=\sqrt{3/2}$ in order to reproduce the exact second moment.