On the equation of motion and the constraining field in ab-initio spin dynamics

It is generally accepted that the effective magnetic field acting on a magnetic moment is given by the gradient of the energy with respect to the magnetization. However, in ab-initio spin dynamics within the adiabatic approximation the effective field is also known to be exactly the negative of the constraining field, which acts as a Lagrange multiplier to stabilize an out-of-equilibrium, non-collinear magnetic configuration. We show that for a physical Hamiltonian both these fields are exactly equivalent, while there can be a finite difference for mean-field Hamiltonians and also for density-functional theory (DFT) calculations. This inequality of the constraining field and the energy gradient is highly relevant for both ab-initio spin dynamics and the ab-initio calculation of exchange constants and effective magnetic Hamiltonians. We argue that the effective magnetic field and exchange constants have highest accuracy in DFT when calculated from the energy gradient and not from the constraining field.


I. INTRODUCTION
The magnetization dynamics of both insulating and metallic materials can in many cases be described within the framework of atomistic spin dynamics (see Refs. [1,2] for an overview). This approach is valid when the electronic Hamiltonian can be mapped onto an effective model of localized spins with constant magnetic moment lengths and interaction parameters that are independent of the spin configuration. While this is generally fulfilled by magnetic insulators, these assumptions may not be valid for magnetic metals, for magnets with non-collinear states [3,4], or for systems far away from equilibrium. An example of the latter is the case of ultrafast demagnetization experiments [5], where a laser pulse demagnetizes a magnetic metal on a sub-picosecond time scale. Such an extreme scenario would require a full nonequilibrium treatment of the electrons, for example with time-dependent density-functional theory (TDDFT) [6]. The numerical difficulty and computational expense of such an approach limits the applications to simulation cells with only a few atoms, which casts doubts on the ability to analyze real experiments. Therefore, it would be desirable to have a method that combines atomistic spin dynamics and electronic structure calculations with a reduced numerical complexity: ab-initio spin dynamics.
Antropov et al. proposed exactly such a formalism where the torques on local magnetic moments are directly calculated from the electronic ground state energy within the adiabatic approximation [7,8]. Stocks et al. pointed out that an arbitrary non-collinear magnetic configuration, that may be formed in such simulations, is not a stable ground state of the energy functional of density functional theory (DFT), and therefore constraining fields are needed to enforce the desired magnetization directions [9,10]. The implementation of arbitrary con-straints within DFT has been worked out previously by Dederichs et al. [11]. The use of constraining fields becomes essential for magnetic configurations that deviate strongly from that of the ground state [12].
Stocks et al. came to the conclusion that the effective field obtained from the energy gradient is exactly the negative of the constraining field, as the constraining field has to cancel the effective field [9,10]. However, to the best of our knowledge, there is no example in the literature where the energy gradient and the constraining field are actually compared to confirm this relation.
Here, we present such calculations for the simple case of an iron dimer, where we find that the equivalence of the constraining field and energy gradient is not exact. Motivated by this surprising numerical result, we derive an exact relation between constraining field and energy gradient: the constraining field theorem. This theorem shows that there is an additional term that can spoil the equality of both fields when the Hamiltonian contains mean-field parameters, which also applies to DFT calculations. We argue that the effective field in DFT should be calculated from the energy gradient and not the constraining field. This implies that exchange constants and effective magnetic Hamiltonians should also be derived from the energy gradient and not from the constraining field.

II. ADIABATIC APPROXIMATION
The adiabatic approximation in ab-initio spin dynamics is based on the assumption that the degrees of freedom can be separated into fast and slow components [7,8,13]. The slow degree of freedom is the magnetization direction, while the fast electronic degrees of freedom, including the magnetic moment lengths, are assumed to equili-brate on much shorter time scales. For the description of the magnetization dynamics it is then valid to consider a quasi-equilibrium state where the magnetic moment directions are held fixed by Lagrange multipliers that act as constraining fields on the magnetic moments. The torques on the magnetic moments can then be calculated from this quasi-equilibrium state.
The constrained Hamiltonian is given bŷ whereĤ 0 is the original Hamiltonian and the constraining term,Ĥ enforces a specific quasi-equilibrium state. Here,Ŝ i is the operator of the total spin at site i, the gyromagnetic ratio γ = −g|e|/(2m e ) with electron spin g-factor g ≈ 2, and B con i is the constraining field. As shown in Fig. 1 this field, combined with intrinsic fields that are present inĤ 0 in Eq. (2.1), acts on an atomic magnetic moment, such that the atomic moment and the constraining field are perpendicular. This causes the field to only constrain the directions of the atomic moments, m i , and not the lengths, M i . While Eq. (2.2) remains finite as an operator, its expectation value vanishes,

III. EQUATION OF MOTION
We consider the equation of motion of the total spin S i at site i, where the expectation values have to be calculated with respect to the ground state ψ of the constrained Hamil-tonianĤ. In general we have for an operatorÔ, where the ground state ψ({e i }) is a function of the prescribed magnetic moment directions e i with the expectation values of the moment directions fulfilling m i = e i . The total torque onŜ i in the ground state of the full Hamiltonian in Eq. (2.1) is zero. This follows from where E 0 is the ground-state energy eigenvalue of H. This implies that for a Hamiltonian with a constraining field one may write We identify which implies that the effective field The constraining field cancels the effective field, as illustrated in Fig. 1. This shows that the correct torque for a given HamiltonianĤ 0 can be obtained from the constraining field. For a classical spin system with localized rigid spins it can be shown that the effective field is exactly [14,15] which would offer an alternative approach of calculating the effective field of atomistic spin-dynamics, B eff i , compared to Eq. (3.7). However, it is not obvious that this should also hold for an itinerant magnet with a nonclassical Hamiltonian that is not a simple function of spin operators, but of the creation and annihilation operators of the itinerant electrons. Below we evaluate the effective field from the two approaches (3.7) and (3.8).

IV. CONSTRAINING FIELD
In the previous section we introduced the constraining field B con i , but we have so far not given a procedure how to calculate this field. We know that the constraining field at site i has to be tuned so that it cancels the intrinsic field of the HamiltonianĤ 0 and is perpendicular to the magnetic moment direction m i . Furthermore, the constraining field has to be chosen such that at each site the moment points along the prescribed moment direction e i , We are going to consider here two methods of calculating the constraining field: the method proposed by Stocks et al. [9,10] and the method by Ma and Dudarev [16]. Stocks et al. provide the following iterative procedure for calculating the constraining field [9,10], where k is the iteration index and B 0 a free parameter that can be tuned for optimal convergence. The first two terms of Eq. (4.2) ensure that only the contribution perpendicular to e i is carried to the next iteration, while the third term adjusts the constraining field by a term proportional to the difference between the output and prescribed moment direction (again only keeping the perpendicular contribution), which aligns the magnetic moment m i closer along the prescribed direction e i . The algorithm Eq. (4.2) can be derived systematically from the method of Lagrange multipliers [17]. The uniqueness of the constraining field follows from the uniqueness of solutions within constrained DFT [18]. Ma and Dudarev derive the constraining field by imposing an energy penalty for misalignments of m i and e i , which leads to the constraining field [16] where λ determines the strength of the energy penalty and convergence is formally reached for λ → ∞ with Both methods look similar, but have different advantages and disadvantages. The first method (4.2) has the advantage of good convergence since the constraining field is adjusted step by step, but requires this additional iterative calculation, which can be done in parallel to a self-consistent calculation. The second method (4.3) has the advantage that it does not introduce an additional iterative calculation and can be simply included in a selfconsistent calculation, but it has the disadvantage that convergence is problematic if λ is too large. This convergence problem can be circumvented by increasing the value of λ in steps, which keeps the energy penalty sufficiently small.  Figure 2. Comparison of the effective magnetic field calculated from the constraining field B con 2,θ (red symbols) and the energy gradient B grad 2,θ (blue symbols) for an iron dimer, where the moment i = 2 is rotated by θ. The inset shows the difference ∆ = −B grad 2,θ − B con 2,θ between the two effective fields.
Within the formalism of constrained DFT the constraining field can be directly included as a Lagrange multiplier in the energy minimization procedure to determine the constrained ground state [11,12,17,19,20].

V. NUMERICAL RESULTS FOR A DIMER
To investigate the relation between the constraining field and the effective field obtained from the energy gradient, we performed DFT calculations for an Fe dimer system with VASP [21,22]. The constraining field implementation in VASP is based on the method by Ma and Dudarev [16], as described above. See Appendix C for more details.
We start with both magnetic moments of the two iron atoms aligned along the z axis and rotate one of the magnetic moments by an angle θ in the xz plane, which gives an expression for the moment of the second atom m x 2 = sin θ, m y 2 = 0, m z 2 = cos θ. The effective field calculated from the energy gradient is given by Since e i is a unit vector with a fixed length, the gradient has to be defined as where θ i and φ i are the polar and azimuthal angles in spherical coordinates with their corresponding unit vectorsθ i andφ i . The θ component of Eq. (5.2) is therefore In Fig. 2 we show the θ component of the constraining field acting on the rotated spin and we compare this field to what one obtains from the energy gradient. The two calculations, Eqs. (3.7) and (5.2), are plotted in Fig. 2 as a function of θ. From the figure one concludes that the two fields are similar, but not exactly identical. It is also possible to discern that the difference between them becomes bigger the further away one is from the equilibrium configuration.
For comparison, we performed analogous calculations with a mean-field tight-binding model (see Appendix A for details), which, as can be seen in Fig. 3, show similar results. There we also show the fieldB grad i which is calculated without constraining fields, with the constraints only implemented approximately by imposing local quantization axes [23] (see Appendix A). This approximate method underestimates in our case the gradient field by about 25%, even in the limit θ → 0. This implies an underestimate by 25% of the exchange parameter J of the dimer in the absence of constraining fields, see Appendix D. The widely used LKAG formalism [24,25] for the calculation of exchange parameters does not take constraining fields into account [26], which could potentially result in similar inaccuracies [27].
In the following, we analyze the origin of the difference between the constraining and energy gradient fields and we argue why it matters to be aware of this difference.

VI. CONSTRAINING FIELD THEOREM
In this section we derive the constraining field theorem, the main result of this paper, which relates the constraining field and the energy gradient field. We discuss under which circumstances these fields are equivalent and the implications this has for the correct choice of the effective field in spin-dynamics simulations.

A. Derivation of the theorem
To relate the energy gradient field B grad i to the constraining field B con i , we wish to evaluate ∇ ei Ĥ 0 +Ĥ con . Here, and in the following, an expression is used, where the first term on the right-hand side accounts for the gradient of the operatorÔ and the second accounts for the gradient (or rotation) of the wavefunction ψ used to calculate the expectation value. We obtain the Hellmann-Feynman theorem [28,29], where the term that accounts for the wavefunction variation vanishes due to the extremity of the ground state energy, 3) The derivative of the constraining part can be expressed as Here one should note that the eigenstate that minimizes the expression in Eq. (6.3), does not necessarily imply that a variation over the wavefunction vanishes when one considers onlyĤ con , i.e. ∇ ψ ei Ĥ con is non-zero. Similarly, both terms need to be considered when treating onlyĤ 0 in the variation, Using the relationship in Eq. (6.3) leads to Since the moment directions for j = i are held constant and M i · B con i = 0, we find  From Eqs. (5.2), (6.6) and (6.7) we arrive at our main result, the constraining field theorem: For a physical Hamiltonian without any approximations there is no explicit dependence on the moment directions, ∇ eiĤ0 = 0, which directly implies This is analogous to the case of an effective spin-Hamiltonian, Eq. (3.8), where the effective field is also given by the energy gradient.

B. Mean-field theories
In a mean-field treatment, the HamiltonianĤ 0 is symbolized withĤ mf , where the parameters of the Hamiltonian in general depend on the directions {e i }. In this case, the term ∇ eiĤmf can be finite, and it follows from Eq. (6.8) that the relation B grad i = −B con i is not exact. Figure 4 shows that this term gives precisely the difference between constraining and energy gradient fields for a mean-field tight-binding calculation (see Appendix A for details), supporting the constraining field theorem (6.8). The difference between the two fields is determined by how strongly the mean-field parameters depend on the moment directions {e i } and by how strong the correlation effects are that are represented by the mean-field contribution toĤ mf . IfĤ mf is derived from a physical Hamiltonian and includes the operator-independent energy contributions arising from the mean-field decoupling, this leads to a cancellation such that ∇ eiĤmf = 0 and B grad i = −B con i , as we demonstrate in Appendix B. Such operator-independent energy contributions are not included in the tight-binding calculations shown in Figs. 3 and 4.

C. Density functional theory
For the DFT calculations, we have to consider the Kohn-Sham (KS) Hamiltonian [30] (here without spinorbit interaction), where l is the index of the KS quasiparticle with position and momentum operatorsr l andp l and spin operatorŜ l . The effective potential V eff is not only dependent on the position and spin of the quasiparticle, but also depends on the average electron and magnetization densities, n(r) and M(r), and we write where V ext is the Coulomb potential from the ions in the lattice and B ext an external magnetic field. The scalar and magnetic exchange-correlation potentials are given by where E xc is the exchange-correlation energy, which is a functional of the electron and magnetization densities. Since these densities depend on the magnetic moment directions {e i }, we have ∇ eiĤKS = 0 and therefore B con i = −B grad i,KS , according to Eq. (6.8) applied toĤ KS . It is important to note here that the constraining field theorem (6.8) is applied to the KS Hamiltonian and the energy gradient field B grad i,KS is therefore given by the gradient of the KS energy Ĥ KS and not of the total DFT energy E DFT , which contains an additional double counting term E dc [25,30], (6.14) This implies for the energy gradient of the total DFT energy with Eq. (6.8) applied to the KS Hamiltonian, where the last two terms of the first line cancel only partially [31] and ∇ * ei denotes the variation at fixed n(r) and M (r) = |M(r)|. Equation (6.15) is the DFT adaptation of the constraining field theorem (6.8) and explains why the DFT calculations in Fig. 2 show a difference between the constraining and energy gradient fields.
The continuous magnetization density of DFT calculations can be mapped onto atomic magnetic moments, (6.16) where Ω i is the volume associated with the atomic site i, e.g., of the corresponding muffin-tin [9].

VII. CHOICE OF THE EFFECTIVE FIELD
For a HamiltonianĤ 0 with ∇ eiĤ0 = 0 it does not matter if the effective field is calculated from the energy gradient or the constraining field. But what is the right choice for the effective field if that is not the case, and in particular, what choice should one make in calculations based on DFT?
Let us first assume that the DFT energy E DFT reproduces exactly the energies of the physical Hamiltonian H 0 , i.e. for each configuration {e i } The correct effective field is then given by the DFT energy gradient, 2) which is not exactly the same as the negative of the constraining field obtained from the KS Hamiltonian, as shown by Eq. (6.15) and Fig. 2. This implies that the constraining field of the KS HamiltonianĤ KS is not the same as the one of the physical HamiltonianĤ 0 and thereforê H KS does not exactly reproduce the equation of motion, When the constraining field is not equivalent to the energy gradient, it is not exact to construct an effective magnetic Hamiltonian based on the calculation of the constraining field alone. The exchange parameters have to be calculated from the energy gradient [24][25][26].
If we are not considering DFT calculations and we cannot make the assumption (7.1), then the argument above does not apply. The effective field describing the magnetization dynamics of a given HamiltonianĤ 0 is then the negative of the constraining field, as shown in Sec. III.

VIII. SUMMARY
We have shown that the effective magnetic field in the equation of motion within the adiabatic approximation is exactly the negative of the constraining field. For Hamiltonians that do not contain mean-field parameters depending on the moment directions, the effective field derived from the energy gradient is equivalent to the constraining field. We have argued that in the case of DFT the effective field should be calculated from the energy gradient because DFT is designed to reproduce the physically correct energies.
Our results have three important implications: 1) In ab-initio spin dynamics the constraining field alone may be insufficient for calculations of the effective field, which should be obtained from the energy gradient.
2) Therefore exchange constants for an effective magnetic Hamiltonian also should be calculated from energy gradients and not from the constraining fields.
3) Our tight-binding calculations support the notion that an approximate implementation of out-ofequilibrium, non-collinear states without constraining fields can give inaccurate results, even in the vicinity of the ferromagnetic ground state. for the d orbitals and I s = I p = I d /10 for the s and p orbitals [34,35]. The spin operator readŝ with Pauli matrix vector σ. Charge neutrality is enforced by the termĤ counts the number of electrons at site i and n 0 i is the prescribed number of electrons per site. We use in our calculations U lcn = 5 eV [34,35].
BothĤ Stoner andĤ lcn depend on the moment configuration {e i }, explicitly and via the charges n i and moment lengths M i , which leads to a difference between B grad i and −B con i according to Eq. (6.8). This tight-binding model already includes an approximate implementation of the constrained moment directions {e i }, since these directions are favored by the Stoner term (A2). Even without constraining fields, the moments align approximately along those directions. This alignment is not exact, as shown in Fig. 5, which implies that constraining fields are still required for accurate results. Our implementation of the constraining fields for the tight-binding model is based on the method by Stocks et al. [9,10], as given by Eq. (4.2).
For simplicity we consider here only the Heisenberg model as an example, which is parameterized by the exchange constants J ij . The corresponding effective magnetic field is which requires knowledge of the magnetic moment length M i . The magnetic moment length itself may depend on the moment configuration, M i = M i ({e i }), which makes this approach only feasible if the moment length can be assumed to be constant. We had to subtract here the component parallel to e i , which follows from the definition of the gradient (5.3) and the requirement that the effective field is perpendicular to e i . For the dimer, the Heisenberg model depends only on a single exchange constant J, and the effective magnetic field is proportional to J, B eff 2 = J M 2 (e 1 − e 2 (e 2 · e 1 )) . (D6)