Numerical Simulations of a Spin Dynamics Model Based on a Path Integral Approach

Inspired by path integral molecular dynamics, we build a spin model, in terms of spin coherent states, from which we can compute the quantum expectation values of a spin in a constant magnetic field, at finite temperature. This formulation facilitates the description of a discrete quantum spin system in terms of a continuous classical model and recasts the quantum spin effects within the framework of path integrals in a double $1/s$ and $\hbar s$ expansion, where $s$ is the magnitude of the spin. In particular, it allows for a much more direct path to the low- and high-temperature limits of the quantum system and to the definition of effective classical Hamiltonians that describe both thermal and quantum fluctuations. In this formalism, the quantum properties of the spins emerge as an effective anisotropy. We use atomistic spin dynamics to sample the path integral, calculate thermodynamic observables and show that our effective classical models can reproduce the thermal expectation values of the quantum system within temperature ranges relevant for studying magnetic ordering.


INTRODUCTION
Spin models of magnetic materials are usually either quantum or classical in terms of the elementary building blocks on which they are based.In quantum spin models, the spin states belong to the quantum space of states that includes all linear superpositions of the eigenstates of Ŝz and Ŝ2 , and the spin variables are quantum operators.By contrast, in classical spin models, 'spin' is used colloquially and actually refers to the classical magnetic moment, µS, where S is usually of fixed length with dynamics confined to the surface of the Bloch sphere and µ is the spin magnetic moment in Bohr magnetons.
Quantum models allow an accurate description of both thermodynamics and dynamics, which intrinsically include purely quantum effects such as entanglement and quantum fluctuations.However, the size of systems that can be studied is often limited to tens or hundreds of spins due to the large computational cost, as solving quantum problems exactly amounts to diagonalization of larger and larger matrices, and even approximation schemes thereof suffer from scaling issues.Numerical methods, such as quantum Monte Carlo (QMC), allow calculations of very large quantum spin systems (hundreds of thousands of spins) with very high accuracy.However, there is no access to dynamical quantities, as QMC is intrinsically a description of thermodynamics, where time is absent.Other quantum methods which do provide access to real-time dynamics cannot provide results for such large systems.Additionally, fundamental issues also arise, such as the 'sign problem' in the case of antiferromagnets, since the Hubbard-Stratonovich transformation leads to an effective Hamiltonian that is not hermitian although the evolution operator is unitary 1 .
Classical spin models are frequently used to study the dynamics and thermodynamics of magnetic materials, helping to interpret experiments at "high" temperatures, where quantum effects-such as entanglement-can be neglected.The computational cost is relatively low, and the formalism is easy to parallelize, leading to routine simulations of the dynamics of hundreds of thousands or even millions of spins.While these classical models give a good qualitative description of the magnetic dynamics, issues arise at lower temperatures, where the assumption of classical Boltzmann statistics is no longer appropriate.The magnon Debye temperature tends to be very high and of the same order as the magnetic ordering temperature, so the 'low-temperature' regime may cover most of the temperature range of magnetic ordering 2,3 .Recent efforts have been made to introduce ad hoc corrections to classical spin models to produce results that more closely resemble quantum models and to better agree with experimental measurements 2,[4][5][6][7][8] .However, these approaches are incapable of including quantum effects, such as tunneling between macroscopic states or zero-point fluctuations.These quantum effects are becoming relevant at ever larger length scales and higher temperatures, for example, with the measurement of the motion of domain walls induced by quantum domain fluctuations in Cr up to 40K 9 .Thus, what is still lacking is a dynamical quantum model whose accuracy can bridge the gap between a fully quantum simulation of a few atoms and an effective classical model and that enables simulations scalable to the size of spintronic device components of millions of spins.
Here, we describe a way of constructing a bridge between quantum and classical spin models by employing a path integral formalism for spin dynamics.This is inspired by path integral molecular dynamics 10 where the efficiency of classical molecular dynamics is used to calculate quantum properties, by establishing the appropriate evolution equations to move in the phase space of the quantum system and thus sample configurations therein 11 .However, how to take into account spin degrees of freedom and sample the corresponding phase space is by no means obvious.
First attempts to do so 12 , in particular for molecular magnets 13 express the spin degrees of freedom in terms of equivalent, though fictitious, position and momentum variables and using the known molecular dynamics formalism in this guise.Hence, these involve mapping the spin Hamiltonian to a particle Hamiltonian.This makes the interpretation of the results in terms of classical magnetic moments, the actual experimental observable, much less straightforward, and this mapping is difficult to build for more complex spin interactions.However, the real problem which we must overcome is that the space of positions and momenta is flat; while the space spanned by the spin degrees of freedom is curved.
It is this problem that is solved by using the basis of spin coherent states 12 .While spin coherent states have been used in some quantum methods 14 , these methods incur a non-trivial cost, for large systems, as well as not being well-suited for extracting the information on the individual (classical) spin components.We note, however, that spin coherent states have also recently been used in methods to derive/rederive equations of motion for magnetization dynamics 15 .Introducing spin coherent states comes at a price: these states are no longer eigenstates of the quantum Hamiltonian.Nonetheless, as we are interested in studying the crossover from quantum to classical behaviour, it is precisely these spin coherent states that are best suited for the task.
In this Article we therefore consider the simplest nontrivial spin system: a single spin in an external magnetic field, described by the Zeeman Hamiltonian.We develop a formalism which uses the spin coherent states and the operators that act on them to to compute themal expectation values of the quantum system in exactly solvable cases, and compare the results obtained to numerical calculations performed with classical atomistic spin dynamics methods, in presence of a field, which takes into account the quantum properties of the spins, when in contact with a thermal bath.We demonstrate that this formalism can indeed account for the quantum properties of the spin, across a broad range of temperatures, with deviations appearing only at "very low" temperatures, as expected by intuition.We emphasise here that, we are not seeking an exact classical equivalent of the quantum system, rather, we are building an effective classical model whose thermal expectation values reproduce those of its quantum counterpart through a dynamical stochastic path sampling method.Moreover, the scope of this paper is not a fundamental study of path integrals for spin systems nor is it placed in the context of geometric quantisation schemes, although the literature from these fields has proven particularly useful for building our model and will be important in future works [16][17][18] .
The plan of the paper is as follows: In Section I we start from a description of a quantum spin system in terms of the discrete spin states |s, m⟩ which are eigenvectors of Ŝz and Ŝ2 and switch to the continuous spin coherent state basis to show that from the quantum model, we can recover a continuous description which can be rewritten in terms of the classical spin vectors S. We do this in a systematic double expansion in 1/s and ℏs; To justify this we explicitly recover the classical limit from this formalism as a sanity check of our approach.In Section II we consider special cases where results can be computed directly from the partition function.We compute expectation values using the spin coherent states for the classical limit (an exact result), for several orders of corrections to this classical limit (under our systematic approximation scheme) and for the exact quantum solution using the discrete |s, m⟩ basis.These results serve as reference and are compared with the results obtained from the new method developed in the next section.In Section III we begin by deriving an effective classical Hamiltonian from the quantum partition function in both low (Section III A) and high (Section III B) temperature limits.In both cases, the resulting effective classical magnetic system is sampled by computing stochastic paths on the Bloch sphere using finite temperature atomistic spin dynamics simulations.In fact, for the system at hand the path integral is an integral over the manifold of all possible superpositions, i.e. over a complex projective space.Finally, we compare results from classical atomistic spin dynamics simulations to results from our new enhanced atomistic model, whilst using the results obtained directly from the partition function (cf.Section II) as reference.We show that indeed, we are able to recover the correct quantum thermal expectation values from this effective classical model for most of the temperature range where there is a significant difference between the classical limit and the quantum solution.In section IV, we summarize our findings and discuss key issues to address in further work.

I. FROM THE SPIN STATES TO THE SPIN COHERENT STATES
A. Partition function in the discrete spin states basis In molecular dynamics, the dynamical variables of the quantum system take values in a flat space.This makes the application of path integrals using classical positions and momenta relatively straightforward.For spin systems, the dynamical variables, the components of spin, take values in a curved space and can only take discrete values due to the discrete spectrum of the spin Hamiltonian {|s, m⟩} , m ∈ −s, s , where s is the principal quantum number and m labels all different possible states with this given spin s.For example, with s = 2 there are 2s + 1 = 5 eigenstates: However, all possible states of a quantum system of spin s = 2 are linear combinations of these five states, i.e. they are described as The normalization of these states implies that the coefficients satisfy the constraint which defines a point on the unit sphere in ten dimensions, but the property that five phases can be modded out reduces this to a five-dimensional manifold.The real challenge is to sample this space efficiently.
The partition function of this quantum spin system is the volume of this five-dimensional manifold, which is finite: Upon coupling the magnetic moment to a thermal bath, the partition function takes the form with β = 1/(k B T ), where k B = 1.381 × 10 −23 J/K is the Boltzman constant and T is the temperature in Kelvin.From Eq. ( 6) it is not obvious how the dynamical behavior of the quantum system, defined over the full manifold, goes over to that of a classical system, localized on the five states {|2, −2⟩ , |2, −1⟩ , |2, 0⟩ , |2, 1⟩ , |2, 2⟩} , in the "classical limit" and how this can be defined.
This requires a careful discussion of what we mean by a 'quantum' system and its classical limit.On the one hand, we have the discrete basis of the eigenstates of the Hamiltonian, but on the other hand, we have the quantum superposition of states which leads to a continuous manifold of possible quantum states.Here, we emphasize that we are dealing with classical measurements of quantum systems, which means that the outcome of any single measurement can only be an eigenstate of our Hamiltonian-which is labeled by an integer for spin systems.The prototype of this situation is the experiment by Stern and Gerlach 19 , where, even though the possible quantum states of the electron can belong to a superposition, such that a 2 + b 2 = 1, the outcome of the measurement of the experiment is either |↑⟩ or |↓⟩.This is in contrast to a classical measurement of the projection along the z-axis of a classical magnetic moment for which a single measurement could take any value between +µ s and −µ s where µ s is the total magnetic moment.Thus, if our Hamiltonian is a function of Ŝz only, then the partition function corresponding to the classical measurement of said quantum system is given as a sum over the eigenstates of this Hamiltonian, rather than an integral over the quantum manifold of states, This expression can be evaluated, especially for the case of a single spin; however defining, let alone studying its classical limit is by no means obvious.It is to this end that it's useful to introduce the spin coherent states.

B. Partition function in the continuous spin coherent state basis
One way to sample the partition function over the quantum space of states, that is particularly useful in studying the crossover to the classical limit, is to recast the system in terms of the so-called spin coherent states 20 .Indeed, not only do the spin coherent states form a continuous basis for the spin system, enabling a mapping onto the continuous description in terms of a unit vector living on a sphere, but it has also been shown that their behavior is close to the classical limit 21 .Thus, they enable us to, on one hand efficiently sample the manifold of quantum states and other other hand to consistenly define the classical limit.The spin coherent states have previously been used to study fundamental aspects such as emerging supersymmetry in spin systems 22 , semiclassical transition probabilities 23 , and energy gap computations within mean-field quantum perturbation theory 24 .
We now proceed by introducing the spin coherent states and showing that the matrix elements of Ŝz can be written as a sum of the classical limit plus corrections.These corrections are essential for including quantum fluctuations into our effective model.We show results for both the purely classical limit of the spin coherent states and how the systematic inclusion of these corrections brings the expectation values closer to their quantum counterparts.
To use the spin coherent states, we work as follows: for a given quantum spin number s, we set where p ∈ {0, 1, . . ., 2s−1, 2s} using the labeling introduced above and we define the spin coherent states |z⟩ , labeled by a complex number z, by the action of the lowering operator 25 , where the 1/ℏ factor is a bookkeeping device needed to keep the exponential dimensionless.Its role in setting the scale of the quantum fluctuations will emerge in what follows.The action of Ŝ+ , Ŝ− and Ŝz on |p⟩ produces The expression in ( 10) is equivalent to which, as we shall see, is more convenient for computing the action of spin operators on the spin coherent states.In this basis, we can write the partition function ( 8) as an integral over the complex label z for the spin coherent states as where the measure must be properly normalized as dµ(z) |z⟩ ⟨z| = 1.In this case C. Crossover from the quantum system to the classical limit To study the quantum system close to the classical limit, we must calculate the matrix elements of Ŝz and its powers on the states |z⟩.The first two powers are In general, it can be shown that the higher-order terms are all of the form The first term is the leading term in the classical limit.If we were to simply approximate we would be discarding all quantum fluctuations.However, it is the systematic inclusion of the quantum fluctuations that we aim to achieve in later Section III.The second term in ( 16) is an example of a correction term, but there is no general, closed expression for the correction terms of increasing order in k.These corrections terms expresses the fact that the manifold of the spin states is curved and is not intrinsically due to the noncommutivity of quantum mechanical operators.Essentially these terms are the difference in the trajectory between states on a flat surface compared to a curved surface; rotations on a classical sphere don't commute.However, in taking quantum states to be all possible superpositions of the basis states, the classical states emerge in the limit ℏ → 0, s → ∞ while keeping the product, s ≡ ℏs fixed.The correction terms are always of the same order in ℏ as the leading term.Thus neglecting these terms does not simply correspond to the semiclassical ℏ expansion and needs to be justified differently.To show this we rewrite equation (16) as which highlights the property that the correction terms, which are sensitive to the curvature of the manifold of spin superpositions, are of higher order in an 1/s expansion; and that the operators, that have a sensible large-spin, i.e. semiclassical, limit are Ŝk z /s k .Indeed, this limit entails taking ℏ → 0, s → ∞ while keeping the product, s ≡ ℏs fixed.It is precisely these corrections that will be refered to in the rest of the text as noncommutative corrections.Indeed, these corrections arise as the spin coherent states are eigenstates of Ŝ− but not of Ŝz , and these operators do not commute.
In addition to this double expansion, we are interested in the dependence of the partition function on β which characterizes the thermal bath with which our quantum system is in equilibrium.To this end we perform a standard high temperature expansion of the partition function, i.e. we rewrite the exponential series e −β Ĥ in powers of β.
Therefore, the corrections to the classical limit we are computing are obtained by a two-fold approximation scheme, both in the noncommutative terms as depicted in (19), and in the high temperature β expansion.
The first term on the right-hand side of (17) (ignoring the noncommutative terms), can be written as an exponential series We now define the Hamiltonian for a single spin (whose electromagnetic properties will be described by its g−factor) in an applied magnetic field that is constant along the zdirection, For the electron, g ≈ 2.002 = |g e | is the absolute value of the electron g-factor, µ B = 9.274 × 10 −23 J/T is the Bohr magneton, ℏ = 1.05457182 × 10 −34 J/K is Planck's constant and B z is the applied magnetic field in Tesla.Choosing a fixed field direction (which can always be taken to be along z) simplifies the calculation by reducing the noncommutativity as we work with the exponential of operators.
To compute the partition function, we again express the exponential as a series and compute the matrix elements ⟨z| exp(−β Ĥ) |z⟩, which, using equation (20), can be approximated by (23) Thus, the matrix elements take the simple form The complex number z (and its conjugate z) can then be mapped onto a unit 2-sphere by defining a unit spin coherent state vector 26 , n, with components and using this we can rewrite the matrix elements (24) as This leads immediately to the definition of the classical Hamiltonian where we identify S = n as the classical spin vector (magnetic moment) with length µ s = sgµ B .Therefore, dropping the noncommutative terms, indeed yields the expected classical limit of this quantum system.We emphasize that all the powers of Ŝk z are needed to recover the classical limit-only the noncommutative terms have been dropped.As we go to the large-spin limit, since the curvature of the sphere is proportional to 1/s 2 , it becomes smaller and smaller, which justifies neglecting these terms.
The vector n defined by the spin coherent states plays the role of the spin unit vector, which is commonly used in classical Heisenberg spin models.Thus, not only does the spin coherent states basis provide us with a continuous (integral) description of the quantum system, but it also yields a straightforward interpretation of the quantum system (described by its states and operators) in terms of the continuous classical system (described by the magnetization vector).We would like to clarify that the convergence of the quantum infinite spin limit towards the classical limit has been rigorously proven long ago, using spin coherent states, in the more general context of the quantum Heisenberg model in the thermodynamic limit 27 , and in more recent work, without using spin coherent states 28 or thermodynamic limit assumptions 29 .However, these approaches are based on constructing lower and upper bounds for the quantum partition function (and/or free energy), which converge to the classical limit in the infinite spin limit and don't aim to build a spin dependent classical approximation, which is our goal in this paper.The cornerstone of our model is the double expansion, on the one hand relative to the curvature of the spin manifold-i.e. the 1/s expansion (which can be understood as a large N or t'Hooft expansion 30 , keeping ℏs fixed), on the oher hand the high temperature expansion in powers of β, for the exponential series (22).It is from the interplay of these two expansions that we obtain the effective classical Hamiltonian, when equilibrium with both baths (of quantum and thermal fluctuations) is assumed.The aim of our approach is to build an efficient numerical method for computing controlled approximations to the thermal expectation values of the quantum system.
We shall now use the partition function in the spin coherent state basis to compute expectation values for the quantum spin Hamiltonian, close to the classical limit, by performing an expansion in increasing orders of β.These serve as a reference to compare to numerical calculations using atomistic spin dynamics in Section III.

II. PARTITION FUNCTION AND EXPECTATION VALUES
The expectation value of an operator Ô for the discrete quantum spin system is where the denominator is the partition function (8).In the spin coherent state basis, the expectation value is expressed in terms of integrals, rather than sums, viz.
As mentioned above, the spin coherent states are not eigenstates of Ŝz , making the exponentiation more subtle.The action of the exponential of Ŝz on |s, m⟩ simply yields the exponentiation of the eigenvalue e Ŝz/ℏ |s, m⟩ = e m |s, m⟩ ; but in the spin coherent state basis, we cannot exactly compute the action and must resort to approximations such as the double expansion described in Section in I expansion and the high-and low-temperature expansions.We proceed by calculating the expectation value ⟨ Ŝz ⟩ as a function of temperature with the Zeeman Hamiltonian (21).This is known to be qualitatively different for classical and quantum spin models due to spin quantisation 31 .The expectation value ⟨ Ŝz ⟩ can be identified with the magnetization induced by an external field (in the limit when the exchange interaction can be neglected).
The exact quantum expectation value, calculated from the discrete basis, where the action of Ŝz |s, m⟩ = ℏm |s, m⟩, gives The expectation value in the classical limit is calculated with the spin coherent states using equation ( 29) and the approximation in equation ( 24) which neglects the terms proportional to powers of 1/s, yielding Using these expressions for the discrete quantum model (31) and the classical limit of the spin coherent state (32), we plot the expectation value ⟨ Ŝz ⟩ as a function of temperature in  31).Blue solid line -the classical limit of the spin coherent state basis from Eq. (32).Dashed lines are successive corrections to the partition function to include noncommutative terms such as appear in Eq. ( 16).The applied field is Bz = 1 T for all figures.
Figure 1.Neglecting the terms due to the non-comutativity of Ŝz and Ŝ± , i.e. working to leading order in the 1/s expansion, means the representation by the spin coherent states produces the classical limit (blue solid line), as expected, with an immediate decay of the spin alignment with the external field as soon as the temperature is non-zero.Equation ( 32) is, in fact, identical to the expectation value ⟨S z ⟩ of a classical spin, as is expected from Ehrenfest's theorem-a useful sanity check (see Appendix A).In the quantum case (red solid line) the expectation value remains almost flat-at low temperatures-and displays a slower characteristic decay around the zero temperature value, along with an initial inflection point that is expected on general grounds 7 .
These characteristic differences between quantum and classical models of single spins are well known and well studied.Of practical interest is that we can obtain an intermediate approximation for the quantum thermal expectation values by retaining terms related to the noncommutativity of operators.Indeed, if we wish to include quantum features into the classical model in a rigourous manner, we cannot neglect all the noncommutative terms by simply using the approximation of (18).To do this, the exponential functions in the spin coherent state expectation value (29) must be expanded as a series in β, Higher-order terms beyond Ŝz contain the effects of the noncommutativity of operators, as seen in ( 16), and we now include these terms as we evaluate the expectation value.We calculate ⟨ Ŝz ⟩ in the spin coherent state basis in increasing orders of the β expansion, which includes the terms due to noncommutativity of Ŝz to higher orders.The results are shown with dashed lines in Figure 1.'1 correction term' includes noncommutative corrections for Ŝ2 z , '2 correction terms' corrections for Ŝ3 z and so on.We see that including even the first noncommuting term in this expansion yields a solution that is already significantly different from the classical result and close to the quantum solution at temperatures of the order of 1 K and above.The agreement improves as the temperature increases, as expected for an expansion in powers of β.Going to higher orders in β causes the expectation value to converge more quickly to the quantum solution (Figure 1), thus producing a continuous description of the discrete quantum system, which is one of our main objectives.
For very low temperatures, close to 0 K, the approximation as a power series in β breaks down and diverges because β is the inverse of the temperature.We emphasize, however, that already at first order in β, this semi-classical model accurately captures the salient features of the thermal spin statistics of the quantum system at temperatures of the order of 1 K. Next, we build a numerical sampling technique for this partition function based on classical, atomistic, spin dynamics.

III. EFFECTIVE HAMILTONIAN AND ATOMISTIC SPIN DYNAMICS A. Low-temperature expansion of the matrix elements
Building a classical Hamiltonian dynamics model to emulate a quantum system, expressed in the spin coherent states basis, requires finding an effective classical Hamiltonian H eff which approximates ⟨z| exp(−β Ĥ) |z⟩ as exp(−βH eff ).By finding such an approximate expression, we recast the quantum system with partition function (8) into an effective classical system with partition function where H eff yields the same expectation values as for the quantum case and μ(z) describes a potentially enlarged, higherdimensional, phase space, as is the case in path integral molecular dynamics approaches 32 .
We consider the partition function with the first noncommutative correction (16), and seek an expression such that where the first term on the right-hand side is the classical limit and the second term is the first noncommutative term which appears on the right-hand side of ( 16).We ignore higher-order non commutation terms in ⟨z| Ŝk z |z⟩, beyond k = 2, keeping only the first noncommutative correction.This is the same level of approximation used in '1 correction term' in Fig. 1.
As a first and very coarse approximation (for more details, see appendix B) we take which, written in terms of the spin coherent state vector n, is The apparent non analyticity in these equations ( 36)- (37), is an artifact of our parametrization.The first term is again the purely classical Zeeman Hamiltonian (27).The second term arises due to the quantization of spin and energetically favors the spin to align with the quantization axis (z).It has a form similar to magnetocrystalline anisotropy, but its origin is the quantum behavior of the spin rather than any physical interaction.We will refer to this term as H Qeff .
To calculate the thermal expectation values using this effective Hamiltonian, we use the techniques of atomistic spin dynamics (ASD) [33][34][35][36][37] .This is usually used to model the dynamics of localized spin magnetic moments µ = µ s S where S is a unit vector and µ s = gsµ B is the size of the spin magnetic moment.The moments interact with a local effective magnetic field B eff obtained from a Hamiltonian H eff that encodes the different magnetic interactions of the system.Here we will use the normalised vector n rather than S to emphasize that we are solving the dynamics of the spin coherent state vector rather than making an a priori assumption of classical spin magnetic moments.
Calculations of the thermodynamic quantities of classical spins can be performed with ASD or Monte Carlo calculations, but ASD is trivial to parallelize across large ensembles of spins, allowing efficient calculation as well as the ability to calculate real-time dynamics.The classical spin dynamics is described by the Landau-Lifshitz-Gilbert (LLG) equation of motion where γ ≡ gµB ℏ is the gyromagnetic ratio in rad • s −1 • T −1 , α is a dimensionless damping parameter, and the effective field B eff in Tesla is calculated as thus, the field from our effective Hamiltonian ( 37) is where e z is the unit vector along z.This expression is apparently singular for n z = 1; this singularity simply indicates that the magnetic field doesn't have any effect on a moment that is aligned with it; we realize, indeed, that such an initial condition, which must be treated separately, is very improbable at any finite temperature.
Temperature is included in the formalism by adding a stochastic field B eff → B eff +η that turns the Landau-Lifshitz-Gilbert equation of motion (38) into a Langevin equation.This is where our method gets its path integral name from.We sample the partition function of this system using several stochastic realisations (or paths) on the Bloch sphere to evaluate the properties of the statistical distribution of the spin vector.The analogue in path integral molecular dynamics methods 10 is using molecular dynamics 38 to sample the partition function.The stochastic field is defined through the fluctuation dissipation theorem, which in the classical case requires η to be a white noise with the properties where i, j are Cartesian components.
In our work the quantum nature of the spin is included directly into the effective field without making any assumption of the statistical distribution.
Recently, stochastic fields using the quantum fluctuation dissipation theorem have been used, enforcing a Bose-Einstein statistical distribution for the noise 2 .This assumes that the relevant thermally occupied objects in this case are magnons, which should obey bosonic statistics.
We numerically integrate the LLG equation ( 38) using a symplectic integration scheme 39 with a timestep of 0.05 ps.The expectation values from the numerical method are calculated as averages over time and multiple realizations of the stochastic dynamics where N s is the number of independent spin trajectories and N t is the number of time samples.The average in time is taken after an equilibration period where the system relaxes from the initial state to a thermalized state.The simulations performed here equilibrate within a few nanoseconds; therefore, we started the averaging procedure after an equilibration period of 5 ns.The averaging time is 15 ns and N s = 20.
From the effective Hamiltonian (36), we compute the expectation value for Ŝz and compare these to results from (42).The results for different values of the principal quantum number s = 1/2, 2, 5 are shown in Figure 2. All three models, classical, quantum and the effective Hamiltonian (43) converge to the same values in the hightemperature limit.In figure 2a for s = 1/2 the effective model differs only slightly from classical model and is not close to the quantum model.Only the slope at zero temperature shows any of the quantum behavior with a small inflection point.This is a feature which several effective models have attempted to force artificially on the studied spin systems to reproduce the experimental behavior for magnetization curves 40 .However, our effective classical atomistic model does not impose any assumptions on the system and has no fitting parameters.The additional computational cost of making the classical system more closely resemble its quantum avatar is minimal, requiring only the addition of a field that amounts to an effective anisotropy.
Although this coarse approximation scheme provides results that are closer to the quantum results, there is no way to systematically improve it.For each higher-order noncommu-tative correction we must again try to derive a H eff ad hoc that satisfies equation (34).Therefore, we continue by developing a more systematic method for which computing the thermal expectation values to higher orders of accuracy is straightforward.

B. High-temperature spin coherent states expansion
The effective model in the previous section produced by approximating the integrand of the partition function by an exponential is very coarse but yields part of the quantum corrections and at a very low computational cost.We now improve on this to try to recover a behavior more similar to the expansion of the partition function in Figure 1.We do this by including higher-order noncommutative terms in the expansion of exp(−β Ĥ) (22) in a more systematic way.
To this end, we return to the partition function ( 13) and, similar to the path-integral molecular dynamics approaches, introduce the resolution of unity as At this stage, the expression is still exact and includes all noncommutative corrections to the classical limit and all orders of temperature.We then approximate (48) with a Taylor expansion as β → 0. Thus in the high-temperature limit (which we later find to be quite low) From the temperature-dependent Hamiltonian (50) and the definition of the effective field (39), we derive (51) We use this effective field in numerical atomistic simulations, and sample several stochastic paths of these effective dynamics over the Bloch sphere.We compare the results with the expectation values computed directly from the partition function (43) and the relevant terms, according to the order of the approximation, of the effective Hamiltonian (50).The results are shown in Figure 3.
When we include only the first correction for the effective field, namely the first and second terms on the right-hand side of (50) then, contrary to the previous section (Figure 2), the low-temperature limit is far from both classical and quantum solutions.However, around 1 K, the results become very close to the quantum solution and converge to be almost identical as the temperature increases.
Including higher-order terms (for example, using all the terms in (51)) we see that although at low temperatures the model is initially further away from the quantum solution, the rate of convergence towards the quantum model is much faster than for lower order corrections.For the first correction, once close to the quantum solution, it takes a while before both curves are indistinguishable, and this happens much quicker when including the second term (see the inset of Figure 3).As our approximation is computed to higher orders, the convergence becomes faster.We note that there is no reason why this high-temperature expansion should become valid at much lower temperatures as we go to higher orders.
Another issue that we have to deal with is that these expectation values have to be normalized in order for the atomistic simulations to overlap with the direct computation from the partition function.Indeed, when we compute the expectation value for ⟨ Ŝz ⟩ we should be using an expression of the form of Eq. ( 29) as but instead (see appendix D), the consistent approximation is given by (53) We know that in the quantum case given by Eq. ( 31), ⟨ Ŝz ⟩ quantum goes to s as β → ∞.We can show that in the same limit, for Eq.(D3), we have hence our expectation values need to be normalized by this factor to yield the correct results (see appendix D for more details).In summary, using this approximation scheme, we can compute expectation values for the quantum system from an equivalent classical atomistic simulation where the quantum nature of the system is represented by a temperaturedependent effective field.This is not a surprise as the space of states is curved.In contrast to the previous section (III A), these then need to be properly rescaled.However, we can compute a closed expression for this rescaling factor, which once again depends only on the principal quantum spin number s.Once this step is fulfilled, the results are almost identical to the fully quantum expectation values for high enough temperatures, which are of the order of 1 K for the single spin in a magnetic field studied here.The low-temperature behavior of this scheme is not as well behaved as in Section III A, which is not surprising, as this is a high-temperature expansion (see Appendix E).

IV. CONCLUSIONS AND OUTLOOK
In this Article, we have built an effective, classical, dynamical model for quantum spin systems from a path integral approach inspired by path integral molecular dynamics in the simplest case of a single spin of arbitrary size in a constant magnetic field described by a Zeeman Hamiltonian.While path integral models of spin have a long history and have been investigated in fundamental contexts such as supersymmetry or, more closely related to our work for molecular magnets, a systematic approach bridging the gap from small-size fully quantum simulations to large-scale dynamical simulations with quantum features has been lacking.Our work here is the first step towards this direction.
We have started by expressing the partition function for spin systems in the spin coherent state basis to obtain a continuous description in terms of an integral rather than a sum, to make the connection to classical spin dynamics.This allows the use of highly efficient atomistic spin dynamics simulations for quantum spin systems and makes the connection between the quantum system defined by its states and the Hamiltonian operator and classical spin dynamics more explicit.We then proceeded to expand the relevant matrix elements of the partition function in powers of β to compute the expectation values of Ŝz directly from the partition function and from atomistic spin dynamics.Here, we have seen that in this first approximation this could be done very simply and efficiently by adding an anisotropic effective field, which could be directly inferred from the quantum spin number of the system.For small spin values, we have seen that the improvement is quite small but increases with the spin.Of course, spin s = 1/2 represents the most extreme limit of spin quantization.As the magnitude of the spin increases to s = 2 and s = 5 (Fig. 2b,c) the corrections in the effective model take the system closer to the quantum solution.Many magnetic materials of practical relevance have s in the range 3/2 to 7/2 so having an improved quantum description for these larger spin values is already very useful.
We also investigated a different method of approximating the integrand of the partition function by an exponential by allowing the effective Hamiltonian of the system to be explicitly temperature-dependent, yielding a temperature-dependent effective field for describing in this way the quantum nature of the system.This method proved to be more accurate for higher temperatures, above 1 K, than the low-temperature expansion, but with the drawback that the expectation values computed using this method require renormalization.However, this renormalization factor has a closed general expression that depends only on the quantum spin number s of the system.
The next step we aim to investigate is the more general case of a general, time-dependent, magnetic field.This introduces more noncommutativity issues with operators Ŝx , Ŝy and Ŝ± .Beyond this, more complex Hamiltonians including the exchange interaction and magnetocrystalline anisotropy in a quantum fashion will allow the large-scale calculation of the thermodyamics of magnetic materials including quantum ef-fects with a relatively low computational cost.In the present case of a constant magnetic field and for a single spin, we have seen that, conversely to path integral methods for molecular dynamics, we did not need to introduce copies of the spin which interact with itself.We do not expect this to hold in more complex Hamiltonians.
It is important to note that despite being a dynamical sampling method, our method provides accurate results for the evalutation of the thermal expectation values, but is not guaranteed to provide accurate real-time dynamics when quantum fluctations drive the system far from the classical limit.In future studies we aim to explore how the dynamical behaviour changes in this context and we expect some fundamental aspects of spin path integrals 41 which apparently do not arise in our model, to resurface for real-time dynamics, even at higher temperatures.

DATA ACCESS
Python code and output data to reproduce all results and figures reported in this paper are openly available from the Zenodo repository: Sources for: Numerical Simulations of a Spin Dynamics Model Based on a Path Integral Approach.https:// doi.org/10.5281/zenodo.7692092 42.The repository contains: • Python code to generate analytic equations derived herein.
• Python code to perform enhanced atomistic spin dynamics calculations with the quantum effective fields.
• Python scripts to reproduce all figures.
The software and data are available under the terms of the MIT License.
(B3) This is where our approximation becomes more qualitative than quantitative.Indeed, the fifth and sixth terms on the righthand side of (B3) are not present in (B1) even though they are not of higher order in β, however, we have taken advantage of the freedom of choice for the sign of the extra term in the effective Hamiltonian (second term on the right-hand side of (B2)) as the correction (third term on the right-hand side of (B1)) comes from the square term in the exponential series.Taking the correction (second term on the right-hand side of (B2)) to be negative implies that or in terms of the spin coherent state vector which means that our expectation value remains close to the classical expectation value, especially for lower temperatures where the spin preferentially aligns with the z-axis.Although this constitutes quite a coarse approximation, it is definitely a relevant primer to understand the subtleties of the path integral spin dynamics method.

FIG. 1 .
FIG.1.Expectation value ⟨ Ŝz⟩ for spin s = 1/2 as a function of temperature.Red solid line -the exact quantum solution in the discrete spin basis |s, m⟩ from Eq. (31).Blue solid line -the classical limit of the spin coherent state basis from Eq.(32).Dashed lines are successive corrections to the partition function to include noncommutative terms such as appear in Eq. (16).The applied field is Bz = 1 T for all figures.

FIG. 2 .
FIG. 2. Expectation value forŜz as a function of temperature for the classical limit (green solid curve), quantum solution (red solid curve) and effective model (blue solid curve) from partition function.Equivalent results from enhanced atomistic spin dynamics simulation for classical limit (purple dashed curve) and effective model (orange dashed curve).(a) Top pane s = 1/2, (b) middle pane s = 2 and (c) bottom pane s = 5

QASD 1
FIG. 3. Expectation value forŜz for s = 2 as a function of temperature for classical limit (green solid curve) and quantum solution (red solid curve) and effective model with the first correction (light blue solid curve) and second correction (dark blue solid curve) from partition function.Equivalent results from enhanced atomistic spin dynamics simulation for classical limit (purple dashed curve) and second effective model with first correction (orange dashed curve) and second correction (yellow dashed curve)