Equations of Correlational Magnetodynamics for Ferromagnetic Materials

A new system of equations for correlational magnetodynamics was developed by means of Bogoliubov hierarchy and new approximation for multiparticle distribution functions. The system consists of two equations. One is Landau–Lifshitz–Bloch like equation, and the other describes the evolution of pair correlations. Computational results show that correlational magnetodynamics model match the direct Landau–Lifshitz better than the standard Landau–Lifshitz–Bloch equation.


1.
Modern studies of magnetic materials demands models taking into account the temperature fluctuations correctly. These demands are due to the progress in nanotechnologies (since nanostructures prone to thermal effects) and advances in governing processes on short time scales. Apart from applying external field there are different ways to achieve magnetization flip. It may be done with aid of spin-polarized current [1], laser [2] or even temperature variation [3].
Full scaled numerical modelling of nonequilibrium processes in magnetic materials taking into account the temperature is required to create cells of magnetoresistive memory which can be flipped in pulse mode by a composition of heating the recording layer and applying spin-polarized current or magnetic field [4].
Models for magnetic materials come in different scales. On the one hand, spin is a quantum phenomenon on the other in certain cases it is preferable to use models which operate on macroscopic scale. The effect in question may depend upon defects in crystal structure, shape or the sample size. The characteristic time may vary from an attosecond to years. Under such conditions it is impossible to use a single model. For example, Monte-Carlo methods are employed to study equilibrium states [5,6]. The most widespread class of models is micromagnetic continuous medium based on Landau-Lifshitz-Gilbert equation [7][8][9][10][11][12][13]. At the moment the most precise commonly accepted micromagnetic model that takes into account evolution of the average magnetization is the Landau-Lifshitz-Bloch equation (LLB) [14][15][16][17][18][19].
Deriving a continuous model from discrete one is a hard problem due to intensive and short ranged exchange interaction and temperature fluctuations which should be taken into consideration. This transi-tion appears to be an important fundamental problem by itself. The mean field approximation which is utilized in LLB derivation produces a number of artifacts due to vanished correlation between magnetic moments of closest atoms. Disappearing of exchange energy in the paramagnetic phase and understated relaxation time for a system [20] are the major drawbacks of this approach. It may be a serious restriction on modelling spintronic and magnetic nanoelectronics devices which are working in a pulse mode.
By choosing another approximation for multiparticle distribution functions one can derive the improved LLB equation that takes into account correlations between the nearest particles. In addition this equation is supplemented with another one describing pair correlations (exchange energy) which is similar to equation on total energy density in hydrodynamics.

2.
We start from discrete model of magnetic material which is governed by the system of Landau-Lifshitz stochastic equations, which describe evolution of N magnetic moments , . Magnetic moments are placed into the nodes of a crystal lattice which have coordinates . The system has the form where γ is the gyromagnetic ratio, α is the damping parameter, H eff is the effective magnetic field, W is the total energy of the system, and T is the temperature, are pairwise independent random vectors composed of δ-correlated random variables with the standard normal distribution, is the gradient operator along the magnetic moment , H ex is the exchange field, J ij is the exchange integral (usually nonzero only for the closest neighbors), H an is the anisotropy field for the easy axis or the easy plane anisotropy, K is the anisotropy parameter, is the unit vector in the anisotropy direction, and H dip is the dipole-dipole (magnetostatic) portion of the field. From now on, we use the dimensionless units.
The system of Eqs.
(1) allows one to take into account thermal fluctuations, crystal structure defects, and a number of other subtle physical phenomena. However it appears to be too computationally demanding for modelling full-sized devices. In order to overcome this restriction, it is necessary to pass to continuous medium equations. The most sophisticated part of this transition is to factor in the exchange interaction between closest neighbors correctly. We are going to use results of numerical integration (1) as ab initio discrete results. The verification of constructed continuous medium equations is carried out on the basis of these results.

Let's introduce notation
Here, is the tangent to the surface of the sphere part of the gradient, is the effective field, D is the diffusion coefficient in the magnetic moment space, and is the distribution function for magnetic moments of a system. Later on, we are going to often use Fokker-Planck-like equations (FPE) [21] hence, the notation will significantly abbreviate the expressions.
Let be the N-particle distribution function in magnetic moments space. It is possible to rigorously enough derive [22] the N-particle FPE from (1) Let's introduce the single-particle distribution function as Here, we denoted integration over a unit sphere by . If we integrate (2) with respect to all but one magnetic moments (see the supplementary material), then we will get the system comprised of N single-particle FPE with integral coefficients where is the two-particle distribution function.

4.
In order to make a complete set of equations one needs to approximate as a functional of and . The standard mean field approximation implies that . In this case, Eq. (3) takes the form Let's assume there is available a smooth (on scales larger than lattice step) six-dimensional distribution function which conforms with at the lattice nodes . We denote corresponding to it distribution of mean magnetization by . In the case only for closest neighbors, the exchange field (7) may be approximated with where a is the length of exchange bond, is the Laplace operator in the configuration space, and is the number of nearest neighbors.

The terms in
have an apparent physical meaning. The first one with describes the interaction larger than the infinitesimal volume for non-homogeneous spacial distribution . It also brings on spin waves. The second term is responsible for interactions inside the infinitesimal volume. It produces spontaneous magnetization and the ferromagnetic-paramagnetic phase transition. The dipoledipole interaction field (8) may be written as integral over the volume, linearly depends on , and has well known methods of calculations [23,24].
After that (6) takes the form (9) where combine the terms which is responsible for interactions on scales larger than infinitesimal volume.

5.
Multiplying Eq. (9) by and integrating it with respect to (see supplementary materials) afterwards, one can obtain LLB equation [14], which describes evolution of : where means tensor multiplication, is the identity matrix, and is the factor introduced by Garanin in order to take into account mean field fluctuations and get the right critical temperature value [25]. To calculate the integral coefficients , , and , which depend on high-order moments of single-particle distribution function, it is necessary to approximate , e.g., as 1 (11) where vector is an approximation parameter, , , and is the Langevin function. The coefficients , , and may be estimated numerically and fitted analytically (see [26] and supplementary materials). In the equilibrium case, the LLB equation (10) conforms with Curie-Weiss theory with the accuracy of approximation (11). 1 This approximation allows to get qualitative results but is not accurate enough for some cases.
In order to take into account correlations between nearest neighbors, we approximate the twoparticle distribution function as (12) where the parameter describes correlation (indirect included) between closest magnetic moments and , the power is required to comply with the condition . The simplest way to calculate is to solve the equation . Such an approximation may be applied to any two-particle distribution function, but we need only one for the closest particles. Under the condition , approximation (12) approaches mean field approximation while . In the opposite extreme case the exponent (12) tends to Dirac delta while . Despite the presence of the power ρ, the dimensionality of the approximation (12) stays the same as for the two-particle distribution function due to .
The measure of pairwise correlations is defined by Hence, the exchange energy per particle equals .
Let's consider an infinitesimal volume. If the mean magnetization is set there are different levels of pairwise correlations still possible. Each corresponds to distinct configuration of magnetic moments system (Fig. 1). When the pairwise correlations (which is the maximum possible value). The other extreme case is the mean field approximation, in which we get while the pairwise correlations has the lowest possible value . The real world system should be somewhere in between those cases. Hypothetically the ferromagnetic material may have and , which corresponds for the most part to the antiparallel orientation of the closest neighbors.

7.
As mentioned in the preceding part, the mean field approximation works adequately for long-range interactions. Hence, the approximation (12) should be taken into account in estimation of the interaction instead of (9). Note that the exchange field inside the infinitesimal volume becomes antidiffusion with the coefficient in the magnetic moment space. The spontaneous magnetization is governed by the competition between antidiffusion and diffusion caused by the temperature fluctuations. The same competition is also responsible for the phase transition point position.
In the end applying (12) instead of (10) results in LLB-like equation in the form where it is convenient to regard the coefficient as a function of two parameters and (Fig. 2).
In the mean field approximation, Fig. 2 turns into a curve on which . The equilibrium condition is ; therefore, . It corresponds to Curie-Weiss theory which yields the critical temperature greater (Fig. 3) than one obtained from discrete particle modelling (1). In fact, it is due to underrated value and overvalued function . To fix this problem in the mean field approximation, one has to introduce a factor . At the same time, if we take into account the pairwise correlations and assumes , then the value of the function for a given decreases. It leads to a closer match of the discrete modelling results (1).
Here, the summation is performed over all neighbors of the jth particle except for the ith particle. Equations (13) and (14) combine into the system of correlational magnetodynamics equations (CMD).

8.
Three-particle distribution functions which depend on the lattice are required to calculate coefficients . We examine ferromagnetic materials with primitive (SC), bcc, and fcc crystal lattices. See Table 1 and Fig. 4.
We are going to use the notation for three-particle distribution functional in the xcc lattice, which has the kth particle on the sth coordination sphere.
First of all, we study which have all three particles on the line, i.e., , , and . If we neglect the indirect correlations between ith and kth particles, we can approximate such functions with . That is basically a generalization of mean field approximation for the three-particle distribution functions. The corresponding coefficient (see supplementary materials) looks like Another extremal case is the fully symmetric distribution function η  Fig. 4. Primitive, bcc, and fcc lattices and the types of three-particle distribution functions . Indices indicate the coordination sphere which accommodates the third particle. For a primitive cubic lattice with we consider symmetric distribution for particles placed at the vertices of the square On its basis it is possible to obtain which depends on three independent parameters , , and . There are only two constraints, namely, on and . We assume that bonds in diagonal terms of are caused by indirect correlations and deviate from bonds on external edges by the value , i.e., .
Functions present the most challenges. The particles in them are placed at the vertices of the irregular octahedron. That rules out the approach we used to approximate . Comparison with discrete modelling results (1) reveals that approximation produces acceptable results. As a result (14) takes the form (15) The integral coefficients Ψ, Λ, , , , just as are functions of and . In contrast the coefficient also depends on temperature . All the coefficients were estimated numerically on the cluster computer K60 at KIAM RAS [27]. Analytical fits were found for coefficients , Ψ, Λ, and 2 (see the supplementary material). The rest , , and are defined on mesh. 3 All the coefficients are available as the part of aiwlib library [28].
9. The results of the discrete numerical modelling (1) were compared to modelling of the standard LLB model (10) and the CMD (13) modelling for various values of temperature T, external fields , and anisotropy K [27]. A portion of these results are shown of Fig. 5. The discrete particle model in all cases had the computation of the integral coefficients for these materials is significantly more computationally demanding then form ferromagnetic materials.
It is also worth to take a shot at adding quadratic terms to in the exponent of the approximation for the single-particle distribution function (11). Such correction will result in dramatic sophistication of the CMD system since it will be necessary to supplement with additional equations on the second order moments for the distribution function. From the other point of view it may decrease the incoherence with the discrete modelling results for ferromagnetic materials.