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

TL;DR

This paper develops a path integral-based classical spin model to efficiently simulate quantum spin systems at finite temperatures, capturing quantum effects as effective anisotropies and reproducing thermal expectation values.

Contribution

It introduces a novel path integral formulation for spin systems that simplifies quantum simulations by mapping them onto effective classical models.

Findings

Effective classical Hamiltonians describe quantum fluctuations.

The model reproduces quantum thermal expectation values.

Facilitates analysis of magnetic ordering at finite temperatures.

Abstract

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.

Equations

$$\{|s,m⟩\},m\in [[-s,s]],$$

$$\left\{|2,-2⟩,|2,-1⟩,|2,0⟩,|2,1⟩,|2,2⟩\right\}.$$

$$|\psi ⟩=c_{-2}|2,-2⟩+c_{-1}|2,-1⟩+c_0|2,0⟩+c_1|2,1⟩+c_2|2,2⟩$$

$$|c_{-2}{|}^2+|c_{-1}{|}^2+|c_0{|}^2+|c_1{|}^2+|c_2{|}^2=1,$$

$$\begin{array}[]{l}\displaystyle{\cal Z}=\int\,d^{2}c_{-2}d^{2}c_{-1}d^{2}c_{0}d^{2}c_{1}d^{2}c_{2}\\ \displaystyle\delta(|c_{-2}|^{2}+|c_{-1}|^{2}+|c_{0}|^{2}+|c_{1}|^{2}+|c_{2}|^{2}-1).\end{array}$$

$$\begin{array}[]{l}\displaystyle{\cal Z}=\int d\psi\langle\psi|e^{-\beta H}|\psi\rangle=\\ \displaystyle\int\,d^{2}c_{-2}d^{2}c_{-1}d^{2}c_{0}d^{2}c_{1}d^{2}c_{2}\\ \displaystyle\delta(|c_{-2}|^{2}+|c_{-1}|^{2}+|c_{0}|^{2}+|c_{1}|^{2}+|c_{2}|^{2}-1)\,e^{-\beta H(c)},\end{array}$$

$$|\psi ⟩=a|\uparrow ⟩+b|\downarrow ⟩,$$

$$\begin{array}{rl}\mathcal{Z}\equiv \mathrm{Tr}(e^{-\beta \hat{\mathcal{H}}})& =\sum _{m=-s}^s⟨s,m|e^{-\beta \hat{\mathcal{H}}[{\hat{S}}_z]}|s,m⟩\\ & =\sum _{m=-s}^s{e}^{-\beta \hat{\mathcal{H}}[m]}\end{array}$$

$$|p⟩\equiv |s,s-p⟩,$$

$$|z⟩\equiv (1+|z{|}^2{)}^{-s}{e}^{z{\hat{S}}_{-}/\hslash }|0⟩$$

$${\hat{S}}_{-}|p⟩$$

$${\hat{S}}_{+}|p⟩$$

$${\hat{S}}_z|p⟩$$

$$|z⟩\equiv {\left(1+|z{|}^2\right)}^{-s}\sum _{p=0}^{2s}{\left(\begin{array}{c}2s\\ p\end{array}\right)}^{1/2}{z}^p|p⟩,$$

$$\mathcal{Z}=\int d\mu (z)⟨z|e^{-\beta \hat{\mathcal{H}}}|z⟩$$

$$d\mu (z)=\frac{2s+1}{\pi }\frac{dz}{{\left(1+|z{|}^2\right)}^2}.$$

$$⟨z|{\hat{S}}_z|z⟩$$

$$⟨z|{\hat{S}}_z^2|z⟩$$

$$⟨z|{\hat{S}}_z^k|z⟩={\left(\hslash s\frac{1-|z{|}^2}{1+|z{|}^2}\right)}^k+\text{corrections}.$$

$$⟨z|{\hat{S}}_z^k|z⟩\approx ⟨z|{\hat{S}}_z{|z⟩}^k,$$

$$⟨z|{\hat{S}}_z|z⟩$$

$$⟨z|{\hat{S}}_z^2|z⟩$$

$${\mathsf{s}}^2\left\{{\left(\frac{1-|z{|}^2}{1+|z{|}^2}\right)}^2+2\frac{1}{s}\frac{|z{|}^2}{(1+|z{|}^2{)}^2}\right\},$$

$$\sum _{k=0}^{\infty }\frac{1}{k!}⟨z|{\hat{S}}_z^k|z⟩\approx \exp \left(\hslash s\frac{1-|z{|}^2}{1+|z{|}^2}\right),$$

$$\hat{\mathcal{H}}=-\frac{g{\mu }_{\mathrm{B}}}{\hslash }{\hat{S}}_z{B}_z$$

$$\exp \left(-\beta \hat{\mathcal{H}}\right)=\sum _{k=0}^{\infty }\frac{1}{k!}{\left(\beta \frac{g{\mu }_{\mathrm{B}}}{\hslash }{\hat{S}}_z{B}_z\right)}^k,$$

$$⟨z|\exp (-\beta \hat{\mathcal{H}})|z⟩$$

$$⟨z|\exp (-\beta \hat{\mathcal{H}})|z⟩\approx \exp \left(\beta g{\mu }_{\mathrm{B}}{B}_{z}s\frac{1-|z{|}^2}{1+|z{|}^2}\right).$$

$$n_x$$

$$n_y$$

$$n_z$$

$$⟨z|\exp (-\beta \hat{\mathcal{H}})|z⟩\approx \exp \left(\beta g{\mu }_{\mathrm{B}}{B}_{z}s{n}_z\right).$$

$${\mathcal{H}}_{\text{eff}}=-g{\mu }_{\mathrm{B}}{B}_{z}s{n}_z=-{\mu }_s\mathrm{B}\cdot \mathrm{S},$$

$$⟨\hat{O}⟩=\frac{\sum _{m=-s}^s⟨s,m|\hat{O}\exp (-\beta \hat{\mathcal{H}})|s,m⟩}{\sum _{m=-s}^s⟨s,m|\exp (-\beta \hat{\mathcal{H}})|s,m⟩},$$

Full text

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 $\hat{S}_{z}$ and $\hat{\bm{S}}^{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, $\mu{\bm{\mathrm{S}}}$, where ${\bm{\mathrm{S}}}$ is usually of fixed length with dynamics confined to the surface of the Bloch sphere and $\mu$ 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 Ceperley and Alder (1986).

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 Barker and Bauer (2019); Barker et al. (2020). 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 Woo et al. (2015); Bergqvist and Bergman (2018); Evans et al. (2015); Barker and Bauer (2019); Anders et al. (2022); Walsh et al. (2022). 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 Shpyrko et al. (2007). 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 Parrinello and Rahman (1984) 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 Habershon et al. (2013). 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 Runeson and Richardson (2020), in particular for molecular magnets Coronado (2019) 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 Runeson and Richardson (2020). While spin coherent states have been used in some quantum methods Bossion et al. (2022), 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 Zhang and Batista (2021). 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 Kochetov (1998); Cabra et al. (1997); Klauder (1979).

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 $\ket{s,m}$ which are eigenvectors of $\hat{\bm{S}}_{z}$ and $\hat{\bm{S}}^{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 ${\bm{\mathrm{S}}}$. We do this in a systematic double expansion in $1/s$ and $\hbar 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 $\ket{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.1) and high (Section III.2) 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

I.1 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

[TABLE]

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:

[TABLE]

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

[TABLE]

The normalization of these states implies that the coefficients satisfy the constraint

[TABLE]

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:

[TABLE]

Upon coupling the magnetic moment to a thermal bath, the partition function takes the form

[TABLE]

with $\beta=1/(k_{B}T)$, where $k_{B}=1.381\times 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 $\left\{\ket{2,-2},\ket{2,-1},\ket{2,0},\ket{2,1},\ket{2,2}\right\},$ 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 Gerlach and Stern (1922), where, even though the possible quantum states of the electron can belong to a superposition,

[TABLE]

such that $a^{2}+b^{2}=1$, the outcome of the measurement of the experiment is either $\ket{\uparrow}$ or $\ket{\downarrow}$. 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 $+\mu_{s}$ and $-\mu_{s}$ where $\mu_{s}$ is the total magnetic moment. Thus, if our Hamiltonian is a function of $\hat{S}_{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,

[TABLE]

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.

I.2 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 Radcliffe (1971). 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 Lee Loh and Kim (2015). 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 Stone (1989), semiclassical transition probabilities Stone et al. (2000), and energy gap computations within mean-field quantum perturbation theory Koh (2018).

We now proceed by introducing the spin coherent states and showing that the matrix elements of $\hat{S}_{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

[TABLE]

where $p\in\{0,1,\ldots,2s-1,2s\}$ using the labeling introduced above and we define the spin coherent states $\ket{z},$ labeled by a complex number $z,$ by the action of the lowering operator 111Of course one can also define these in terms of the raising operator or any linear combination of these Nemoto (2000), $\hat{S}_{-}=\hat{S}_{x}-i\hat{S}_{y}$, as

[TABLE]

where the $1/\hbar$ 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 $\hat{S}_{+}$, $\hat{S}_{-}$ and $\hat{S}_{z}$ on $\ket{p}$ produces

[TABLE]

The expression in (10) is equivalent to

[TABLE]

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

[TABLE]

where the measure must be properly normalized as $\int d\mu(z)\ket{z}\bra{z}=1$. In this case

[TABLE]

I.3 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 $\hat{S}_{z}$ and its powers on the states $\ket{z}$. The first two powers are

[TABLE]

In general, it can be shown that the higher-order terms are all of the form

[TABLE]

The first term is the leading term in the classical limit. If we were to simply approximate

[TABLE]

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 $\hbar\to 0,s\to\infty$ while keeping the product, ${\sf s}\equiv\hbar s$ fixed. The correction terms are always of the same order in $\hbar$ as the leading term. Thus neglecting these terms does not simply correspond to the semi-classical $\hbar$ expansion and needs to be justified differently. To show this we rewrite equation (16) as

[TABLE]

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. semi-classical, limit are $\hat{S}_{z}^{k}/{\sf s}^{k}$. Indeed, this limit entails taking $\hbar\to 0,s\to\infty$ while keeping the product, ${\sf s}\equiv\hbar 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 $\hat{S}^{-}$ but not of $\hat{S}_{z}$, and these operators do not commute.

In addition to this double expansion, we are interested in the dependence of the partition function on $\beta$ 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^{-\beta\hat{\cal H}}$ in powers of $\beta$.

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 $\beta$ expansion.

The first term on the right-hand side of (17) (ignoring the noncommutative terms), can be written as an exponential series

[TABLE]

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 $z$-direction,

[TABLE]

For the electron, $g\approx 2.002=|g_{e}|$ is the absolute value of the electron $g$-factor, ${\mu_{\mathrm{B}}}=9.274\times 10^{-23}$ J/T is the Bohr magneton, $\hbar=1.05457182\times 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

[TABLE]

and compute the matrix elements $\bra{z}\exp(-\beta\hat{\cal H})\ket{z}$, which, using equation (20), can be approximated by

[TABLE]

Thus, the matrix elements take the simple form

[TABLE]

The complex number $z$ (and its conjugate $\bar{z}$) can then be mapped onto a unit 2-sphere by defining a unit spin coherent state vectorKarchev (2012), ${\bm{\mathrm{n}}}$, with components

[TABLE]

and using this we can rewrite the matrix elements (24) as

[TABLE]

This leads immediately to the definition of the classical Hamiltonian

[TABLE]

where we identify ${\bm{\mathrm{S}}}={\bm{\mathrm{n}}}$ as the classical spin vector (magnetic moment) with length $\mu_{s}=sg{\mu_{\mathrm{B}}}$. Therefore, dropping the noncommutative terms, indeed yields the expected classical limit of this quantum system. We emphasize that all the powers of $\hat{S}_{z}^{k}$ 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 ${\bm{\mathrm{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 limitLieb (1973), and in more recent work, without using spin coherent statesConlon and Solovej (1990) or thermodynamic limit assumptionsMillard and Leff (2003). 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’t Hooft (1974), keeping $\hbar s$ fixed), on the oher hand the high temperature expansion in powers of $\beta$, 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 $\beta$. 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 $\hat{O}$ for the discrete quantum spin system is

[TABLE]

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.

[TABLE]

As mentioned above, the spin coherent states are not eigenstates of $\hat{S}_{z}$, making the exponentiation more subtle. The action of the exponential of $\hat{S}_{z}$ on $\ket{s,m}$ simply yields the exponentiation of the eigenvalue

[TABLE]

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 $\langle\hat{S}_{z}\rangle$ 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 Greiner et al. (2000). The expectation value $\langle\hat{S}_{z}\rangle$ 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 $\hat{S}_{z}\ket{s,m}=\hbar m\ket{s,m}$, gives

[TABLE]

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

[TABLE]

Using these expressions for the discrete quantum model (31) and the classical limit of the spin coherent state (32), we plot the expectation value $\langle\hat{S}_{z}\rangle$ as a function of temperature in Figure 1. Neglecting the terms due to the non-comutativity of $\hat{S}_{z}$ and $\hat{S}_{\pm},$ 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 $\langle S_{z}\rangle$ 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 Anders et al. (2022).

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 $\beta$,

[TABLE]

Higher-order terms beyond $\hat{S}_{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 $\langle\hat{S}_{z}\rangle$ in the spin coherent state basis in increasing orders of the $\beta$ expansion, which includes the terms due to noncommutativity of $\hat{S}_{z}$ to higher orders. The results are shown with dashed lines in Figure 1. ‘1 correction term’ includes noncommutative corrections for $\hat{S}_{z}^{2}$, ‘2 correction terms’ corrections for $\hat{S}_{z}^{3}$ 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 $\beta$. Going to higher orders in $\beta$ 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 [math] K, the approximation as a power series in $\beta$ breaks down and diverges because $\beta$ is the inverse of the temperature. We emphasize, however, that already at first order in $\beta$, 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

III.1 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 ${\cal H}_{\textrm{eff}}$ which approximates $\bra{z}\exp(-\beta\hat{{\cal H}})\ket{z}$ as $\exp(-\beta{\cal H}_{\textrm{eff}})$. By finding such an approximate expression, we recast the quantum system with partition function (8) into an effective classical system with partition function

[TABLE]

where ${\cal H}_{\textrm{eff}}$ yields the same expectation values as for the quantum case and $\tilde{\mu}(z)$ describes a potentially enlarged, higher-dimensional, phase space, as is the case in path integral molecular dynamics approachesDeymier et al. (2016).

We consider the partition function with the first noncommutative correction (16), and seek an expression such that

[TABLE]

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 $\bra{z}\hat{S}_{z}^{k}\ket{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

[TABLE]

which, written in terms of the spin coherent state vector ${\bm{\mathrm{n}}}$, is

[TABLE]

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 ${\cal H}_{\textrm{Qeff}}$.

To calculate the thermal expectation values using this effective Hamiltonian, we use the techniques of atomistic spin dynamics (ASD) Halilov et al. (1998); Chubykalo et al. (2003); Mryasov et al. (2005); Skubic et al. (2008); Evans et al. (2014). This is usually used to model the dynamics of localized spin magnetic moments ${\bm{\mathrm{\mu}}}=\mu_{s}{\bm{\mathrm{S}}}$ where ${\bm{\mathrm{S}}}$ is a unit vector and $\mu_{s}=gs\mu_{B}$ is the size of the spin magnetic moment. The moments interact with a local effective magnetic field ${\bm{\mathrm{B}}}_{\textrm{eff}}$ obtained from a Hamiltonian ${\cal H}_{\textrm{eff}}$ that encodes the different magnetic interactions of the system. Here we will use the normalised vector ${\bm{\mathrm{n}}}$ rather than ${\bm{\mathrm{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

[TABLE]

where $\gamma\equiv\frac{g{\mu_{\mathrm{B}}}}{\hbar}$ is the gyromagnetic ratio in $\text{rad}\cdot\text{s}^{-1}\cdot\text{T}^{-1}$, $\alpha$ is a dimensionless damping parameter, and the effective field ${\bm{\mathrm{B}}}_{\textrm{eff}}$ in Tesla is calculated as

[TABLE]

thus, the field from our effective Hamiltonian (37) is

[TABLE]

where ${\bm{\mathrm{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 ${\bm{\mathrm{B}}}_{\textrm{eff}}\to{\bm{\mathrm{B}}}_{\textrm{eff}}+{\bm{\mathrm{\eta}}}$ 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 methodsParrinello and Rahman (1984) is using molecular dynamics Chen and Kim (2004) to sample the partition function. The stochastic field is defined through the fluctuation dissipation theorem, which in the classical case requires ${\bm{\mathrm{\eta}}}$ to be a white noise with the properties

[TABLE]

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 Barker and Bauer (2019). 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 Thibaudeau and Beaujouan (2012) 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

[TABLE]

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 $\hat{S}_{z}$

[TABLE]

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 high-temperature 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 Kuz’min (2005). 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 noncommutative correction we must again try to derive a $\mathcal{H}_{\mathrm{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.

III.2 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(-\beta\hat{\mathcal{H}})$ (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

[TABLE]

in the $\ket{s,m}$ basis, in which $\hat{S}_{z}$ is diagonal, resulting in

[TABLE]

Using the definition of $\ket{z}$ and the action of $\hat{S}_{z}$ on $\ket{p}$ we find

[TABLE]

for which we need to rewrite the integrand

[TABLE]

as a single exponential of the form $F[\beta,z]\equiv\exp({-\beta{\cal H}_{\textrm{eff}}})$ in order to identify an effective Hamiltonian. Through a series of identities (see appendix C), we can write

[TABLE]

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 $\beta\to 0$. Thus in the high-temperature limit (which we later find to be quite low)

[TABLE]

Mapping to the spin coherent state vector components using $(1-|z|^{2})/(1+|z|^{2})=n_{z}$ and $|z|^{2}/(1+|z|^{2})=(1-n_{z}^{2})/4$, we can write a temperature-dependent effective Hamiltonian:

[TABLE]

From the temperature-dependent Hamiltonian (50) and the definition of the effective field (39), we derive

[TABLE]

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 $\braket{\hat{S}_{z}}$ we should be using an expression of the form of Eq. (29) as

[TABLE]

but instead (see appendix D), the consistent approximation is given by

[TABLE]

We know that in the quantum case given by Eq. (31), $\braket{\hat{S_{z}}}_{\textrm{quantum}}$ goes to $s$ as $\beta\to\infty$. We can show that in the same limit, for Eq. (68), we have

[TABLE]

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 temperature-dependent effective field. This is not a surprise as the space of states is curved. In contrast to the previous section (III.1), 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.1, 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 $\beta$ to compute the expectation values of $\hat{S}_{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 $\hat{S}_{x}$, $\hat{S}_{y}$ and $\hat{S}_{\pm}$. 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 effects 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 integralsZwanziger et al. (1990) 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 Nussle et al. (2023). 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.

Author Contributions

Thomas Nussle: conceptualization, methodology, investigation, software, writing - original draft. Stam Nicolis: methodology, writing - review and editing. Joseph Barker: conceptualization, methodology, software, data curation, writing - review and editing, funding acquisition.

acknowledgments

This work was supported by the Engineering and Physical Sciences Research Council [grant number EP/V037935/1]. JB acknowledges funding from a Royal Society University Research Fellowship. The authors thank A. Sylla, F. Labéy and T. Raujouan for very insightful mathematical discussions, as well as J. Hodrien and A. Coleman from the University of Leeds Research Computing team for their help with optimizing the Python code on which this work is relying.

Appendix A Correspondence of the spin coherent states with the classical limit

Here we show that the observable $\braket{\hat{S}_{z}}$ from the spin coherent states with the commutators neglected (i.e. in the classical limit (32)) is identical to $\braket{S_{z}}$ calculated from the classical Heisenberg model. For a classical Heisenberg spin with Hamiltonian

[TABLE]

where ${\bm{\mathrm{S}}}$ lives on the unit sphere, the partition function is

[TABLE]

for which the expectation value of the $z$-component of ${\bm{\mathrm{S}}}$ is given by

[TABLE]

If the external field is constant along the $z$-direction then we have

[TABLE]

as the integrals over $S_{x}$ and $S_{y}$ in the numerator and denominator cancel each other out. Comparing this to $\braket{\hat{S}_{z}}$ for the spin coherent state (32) and using $n_{z}=(1-|z|^{2})/(1+|z|^{2})$ and $\mu_{s}S_{z}=gs{\mu_{\mathrm{B}}}n_{z}$ we see that (58) and (32) are identical up to a factor of $\hbar$, as the classical spin vector has no units, whereas the quantum expectation value of $\braket{\hat{S}_{z}}$ is in units of $\hbar$.

Appendix B Coarse approximation method

We expand the operator exponential series (22) up to second order in $\beta$

[TABLE]

we can show that by taking

[TABLE]

and expanding the effective classical exponential up to the same order in $\beta$, we get

[TABLE]

This is where our approximation becomes more qualitative than quantitative. Indeed, the fifth and sixth terms on the right-hand side of (61) are not present in (59) even though they are not of higher order in $\beta$, 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 (60)) as the correction (third term on the right-hand side of (59)) comes from the square term in the exponential series. Taking the correction (second term on the right-hand side of (60)) to be negative implies that

[TABLE]

or in terms of the spin coherent state vector

[TABLE]

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.

Appendix C High temperature model exponential form

Starting from (47), we rewrite

[TABLE]

hence (47) can be rewritten as

[TABLE]

Appendix D High-temperature model normalization

We approximate

[TABLE]

as our approximation scheme for the partition function aims to move from a quantum description in terms of states and operators to a classical description

[TABLE]

Within this approximation, we can rewrite

[TABLE]

which is the expression we use for our averages, as it corresponds to the same approximation as the atomistic model, as proven by the exact overlap of both the averages computed from the partition function (53) and the atomistic average over time and the number of realizations (42).

What is of peculiar interest is that the ratio

[TABLE]

which reminds us of the fact that the eigenvalues of $\hat{{\bm{\mathrm{S}}}}^{2}$ are $s(s+1)$ as in

[TABLE]

rather than simply $s^{2}$. Indeed, in the classical limit $s\to\infty$ we recover

[TABLE]

We would like to emphasize that this required normalization factor is identical for both the results of the atomistic simulations (42) and the results from the approximate partition function (53).

The expectation values for $\braket{\hat{S}_{z}}_{\textrm{app}}$ with and without normalization are given in Figure 4, along with the appropriate quantum solution.

This is very important for more general applications of this model as this means that the normalization of the curves does not require an additional fitting parameter of any kind but is rather analytically computable and has a general, closed expression.

Appendix E Higher order correction for the high-temperature model

As mentioned in section III.2 our method can technically carry out this approximation scheme to any order in the noncommutative terms, numerically, without requiring to compute these corrections using pen and paper. But as this relies on a Taylor expansion around the high-temperature limit $\beta\to 0$ there is a limit as to how low in temperature we can provide accurate results. Indeed there is no reason for this high-temperature expansion to converge to the quantum solution for temperatures around [math] K. This is shown in Figure 5.