The Einstein-de Haas (EdH) effect, where the spin angular momentum of electrons is transferred to the mechanical angular momentum of atoms, was established experimentally in 1915. While a semiclassical explanation of the effect exists, modern electronic structure methods have not yet been applied to model the phenomenon. In this paper, we investigate its microscopic origins by means of a noncollinear tight-binding model of an O_{2} dimer, which includes the effects of spin-orbit coupling, coupling to an external magnetic field, and vector Stoner exchange. By varying an external magnetic field in the presence of spin-orbit coupling, a torque can be generated on the dimer, validating the presence of the EdH effect. The avoided energy level crossings and the rate of change of magnetic field determine the evolution of the spin. We also find that the torque exerted on the nuclei by the electrons in a time-varying *B* field is not only due to the EdH effect. The other contributions arise from field-induced changes in the electronic orbital angular momentum and from the direct action of the Faraday electric field associated with the time-varying magnetic field.

## I. INTRODUCTION

Materials in burning plasma tokamak fusion devices are expected to be exposed to irradiation in the presence of high strength magnetic fields approaching ∼10 T. In addition to stresses that such fields generate in the reactor components,^{1} they may also affect the dislocation microstructure and plasticity of materials^{2} as well as diffusion of radiation defects. Current tokamak designs, including ITER, use austenitic Fe–Cr–Ni steels^{3} that are nonmagnetic on the macroscopic scale while being microscopically antiferromagnetic.^{4} In a demonstration fusion reactor (DEMOnstration Power Station),^{5} blanket modules are expected to be manufactured from ferromagnetic ferritic-martensitic steels, exhibiting superior resistance to radiation damage. The primary component of a ferritic steel, iron, is strongly ferromagnetic in the body-centered cubic (bcc) phase. The observed directional anisotropy of magnetism in iron^{6–8} suggests that atomic magnetic moments are coupled to the lattice, giving rise to a magnetic contribution to electric and thermal resistivity^{9} and to spin-electron-lattice relaxation effects in collision cascades.^{10}

To help develop the theory of the interaction between electrons and a crystal lattice in a magnetic material in a systematic way, here we consider the Einstein-de Haas (EdH) effect, in which the flipping of spins generates a macroscopic torque on the crystal lattice. This was first predicted in 1908 by Richardson and has been demonstrated numerous times since then.^{11–13} In Einstein and de Haas’s^{14} famous experiment, an iron cylinder was suspended by a thin wire inside a solenoid whose axis coincided with that of the cylinder. Varying the magnetic field within the solenoid led to the generation of a measurable torque on the cylinder.

The present-day interpretation of this result is that a change in the external magnetic field changes the magnetization of a ferromagnetic material by realigning the electronic spins. The dimensionless gyromagnetic factor, *g*′, relates the change in magnetization Δ**M** to a change in electronic angular momentum Δ**J** according to

where *e* is the elementary positive charge and *m*_{e} is the mass of an electron. Conservation of angular momentum requires an equal and opposite change in the mechanical rotational angular momentum of the body of the material. Many experiments have used the EdH effect to measure *g*′ for ferromagnetic materials.^{15}

The EdH effect relies on the transfer of angular momentum from the electronic spins to the lattice of nuclei and thus highlights the role of spin-orbit coupling (SOC) in spin-lattice interactions.^{16} Without SOC, the direction of the electronic spin is decoupled from the orientation of the lattice and the EdH effect does not occur. We note that a good description of SOC is also necessary for the calculation of heat transport coefficients in ferromagnetic materials.^{10,17}

While it may not yet be clear how large an effect magnetism has on thermal conductivity of iron and steels, it should be recalled that the effect on mechanical properties is now known to be significant. For example, nonmagnetic density functional theory (DFT) calculations predict that He interstitials in Fe have the octahedral site as their most favorable position, while magnetic DFT calculations predict that the tetrahedral position is lower in energy.^{18} Self-interstitial defects provide another good example: a self-interstitial atom in magnetic bcc iron adopts a ⟨110⟩ dumbbell structure,^{19} which is different from the linear ⟨111⟩ crowdion configuration adopted by self-interstitials in nonmagnetic bcc transition metals.^{20,21} Perhaps even more significantly, magnetic fluctuations give rise to a sequence of *α*−*γ*−*δ* bcc-fcc-bcc phase transitions in iron and steels at elevated temperatures^{22,23} that require computing free energy differences between competing magnetic phases with megaelectron volt accuracy.

The limitations of our current understanding of magnetic materials are evident from the difficulties in predicting ferromagnetic phenomena such as the sudden collapse of the magnetic moment in hcp-Fe under pressure.^{24} The lack of understanding also hinders attempts to construct microscopic models of heat flow in ferromagnetic materials^{9,25} suitable for use in large-scale atomistic or semiclassical simulations. Our work represents a first step toward addressing this issue.

Although the EdH effect is supported by multiple experiments and illustrates an important way in which magnetic degrees of freedom interact with the lattice degrees of freedom in a material, almost no attempt has been made so far to describe it at a fully quantum-mechanical microscopic level. The authors are unaware of any previous work performing atomic-scale simulations of the EdH effect. However, a few existing studies discuss versions of the EdH effect in a quantum mechanical context. For example, in a numerical model of a pair of dysprosium atoms placed in a spherical harmonic trap interacting through a dipole-dipole interaction term, large orbital angular momenta (*l* > 20) were generated by slow variation of an external magnetic field.^{26} A recent study of a single-molecule magnet coupled to a nanomechanical resonator constitutes an experimental realization of the EdH effect for a single molecule.^{27}

Magnetic DFT calculations can be expensive to perform and have only rarely been used to study the coupled dynamics of atoms and spins. By contrast, magnetic tight binding (TB) provides a simpler and quicker approach capable of describing much of the same physics. The anomalous iron–chromium mixing enthalpy was explained simply using a fixed-moment model based on TB.^{28} Magnetic TB was also used to show that, when magnetic defects are present, the energy surfaces of spin-lattice systems can be highly complex, including multiple separate regions with different spin structures.^{29}

In this paper, we study the EdH effect for an O_{2} dimer subject to a time-varying magnetic field. We use a noncollinear TB model that includes coupling to the external magnetic *B* field, SOC, and vector Stoner exchange.^{30} Forces on the atoms are calculated using the Hellman-Feynman theorem.^{31,32} Once the existence of the EdH effect has been established, we investigate the effect of varying the SOC strength. If SOC is neglected, the component of the spin operator parallel to the applied magnetic field **B** commutes with the Hamiltonian and $S\u22c5B^$ is a constant of the motion. For small values of the SOC strength, the EdH torque is correspondingly small, as expected. We find, however, that the EdH effect is not the only torque-generating mechanism experienced by the nuclei due to the electrons under a time-varying *B* field: the field also couples to the electron orbits, and directly to the nuclei via the Faraday effect. We note that the magnetic TB model used here is extensible to larger systems, including solids.

## II. THEORY

Figure 1 shows a schematic diagram of the O_{2} dimer studied in this paper. The dimer axis $\zeta ^$ is held fixed along the $z^$ direction. The internuclear separation is *R*_{nuc}. The magnetic field is aligned along $x^$ and is spatially homogeneous, although its strength may vary in time. The interaction torque exerted on the nuclei by the electrons is given by

where atom 1 is below the *xy*-plane, atom 2 is above it, *F*_{1} is the force exerted by the electrons on nucleus 1, and *F*_{2} is the force exerted by the electrons on nucleus 2.

To create a TB model of O_{2}, it is helpful to know some of its electronic properties. Of the 16 electrons in an O_{2} molecule, 4 fill the 1s states, 4 fill the 2s states, and 8 occupy bonding and antibonding molecular orbitals (MOs) made from the 2*p* atomic orbitals (AOs). A diagram of the electronic structure of the 2*p* levels of the molecule is shown in Fig. 2.

The TB basis functions used here are atomic (or more precisely, atomic-like) *p*_{x}, *p*_{y}, and *p*_{z} orbitals, with separate orbitals for up and down spins. The six basis functions on each atom are denoted as

The basis set does not include the 1*s* and 2*s* orbitals below the 2*p* shell or any orbitals above it, so the dimer Hamiltonian is a 12 × 12 Hermitian matrix and the 12 molecular orbitals (MO) are occupied by 8 electrons. The MOs *ϕ*_{n} are obtained by diagonalizing the Hamiltonian matrix and can be written as linear combinations of the basis functions *χ*_{ασ} with expansion coefficients *d*_{nασ},

where *α* runs over the spatial AOs on both atoms and *σ* runs over spin states.

Consider first a dimer Hamiltonian containing on-site and hopping matrix elements only, without SOC, exchange interactions, or an applied magnetic field. The 12 2*p* MOs are labeled as follows:

where the *σ*, *π*_{x}, and *π*_{y} states are the bonding combinations of the *p*_{z}, *p*_{x}, and *p*_{y} AOs, respectively, with *σ*^{*}, $\pi x*$, and $\pi y*$ being the equivalent antibonding states. As shown in Fig. 2, the 12 MOs have only four different energies. The highest occupied shell of molecular spin-orbitals, the *π*^{*} shell, is occupied by 2 electrons.

Although in a purely mean field picture, the 4 *π*^{*} MOs, $\pi x,\u2191*$, $\pi x,\u2193*$, $\pi y,\u2191*$, and $\pi y,\u2193*$, are degenerate, it is known that the ground state of O_{2} is a triplet, with the spins of the two *π*^{*} electrons aligned. The alignment is caused by the exchange interaction, which splits the four degenerate *π*^{*} spin-orbitals into a threefold degenerate spin-1 triplet and a spin-0 singlet. Thus, a realistic model of ground state O_{2} must take the exchange interaction into account.

### A. Time evolution

The initial molecular orbitals *ϕ*_{n} obtained by diagonalizing the Hamiltonian matrix at the beginning of a simulation evolve into time-dependent molecular orbitals *ψ*_{n}(*t*) according to the time-dependent Schrödinger equation, *iℏ∂*_{t}*ψ*_{n}(*t*) = *H*(*t*)*ψ*_{n}(*t*), subject to the initial condition *ψ*_{n}(*t* = 0) = *ϕ*_{n}. The Hamiltonian *H*(*t*) depends on time if the applied *B* field depends on time, so *ψ*_{n}(*t*) is not in general an exact eigenfunction of *H*(*t*), when *t* > 0. The set of time-evolved MOs does, however, remain orthonormal. The time-dependent expansion coefficients $dn\alpha \sigma $(*t*) are defined by

and satisfy the following discrete equivalent of the time-dependent Schrödinger equation:

The method used to evolve this equation in time is time-reversible and unitary. Rewriting the equation of motion, Eq. (6), with (*D*)_{nασ} = *d*_{nασ} and (*H*)_{α′σ′,ασ} = *H*_{α′σ′,ασ} as matrices, we obtain

Letting *dt* be an infinitesimal positive time interval gives^{33}

Thus, a finite step of the time evolution can be performed approximately by

The initial values of the expansion coefficients of the time-dependent MO *ψ*_{n}(*t*) are obtained from the corresponding eigenfunction *ϕ*_{n} of the initial Hamiltonian *H*(*t* = 0). Since *ψ*_{n}(0) = *ϕ*_{n}, we have

The time step used in all simulations is 4 hartree a.u., which corresponds to approximately 97 as.

### B. The Hamiltonian

Describing the EdH effect requires a Hamiltonian that incorporates (i) the coupling of electrons to a magnetic field, (ii) arbitrary electron spin directions and magnitudes, which may depend on spatial position, and (iii) spin-orbit coupling. To obtain physically meaningful results for O_{2}, it is also necessary to include exchange interactions. A noncollinear TB method can accommodate all of these requirements.

The Hamiltonian employed is

where *H*_{0} is the basic TB dimer Hamiltonian used to obtain Fig. 2, *H*_{B} is the coupling of the electrons to the magnetic field, *H*_{SOC} is the SOC term, and *H*_{ex} is the exchange term.

*H*_{B} is taken from the standard Pauli Hamiltonian and is given by −** μ** ·

**(**

*B**t*), where

is the total magnetic moment, **L** and **S** are the canonical orbital and spin angular momentum operators in the Coulomb gauge, *μ*_{B} is the Bohr magneton, and ** B**(

*t*) is the external magnetic field acting on the molecule.

SOC is a relativistic effect, and the SOC term in the Hamiltonian may be derived from the Dirac equation via the Foldy-Wouthuysen transformation,^{34} which, for spherical potentials, gives

where *m*_{e} is the mass of an electron, *c* is the speed of light, and *V* (*t*) is the potential experienced by the electron due to the atomic nucleus and the other electrons belonging to that atom in the central field approximation. The gradient of the nuclear potential is largest very close to the nucleus, so the spin-orbit term can be assumed to couple atomic orbitals on the same atom only. The radial part of the SOC matrix element between two orbitals in the same shell on a given atom may then be approximated by a constant given in units of energy, *ξ*, so that^{35}

This is the form of *H*_{SOC} used in this work. Subsection II C describes the Stoner exchange term, *H*_{ex}.

### C. Vector Stoner exchange

All of the results reported in this paper have the effects of the vector Stoner exchange term (*H*_{ex}) included. The more familiar collinear form of the Stoner exchange term is inappropriate for use in describing spin dynamics as it breaks rotational symmetry in spin space.^{30} The vectorial exchange interaction is treated at the mean-field (Hartree-Fock) level, leading to a self-consistent independent-electron problem.

The remainder of this section outlines the vector Stoner exchange implementation. The exchange moment vector **m**_{a}(*t*) on atom *a* (∈{1, 2}) at time *t* is evaluated as

where ** σ** is the vector of Pauli spin matrices, $tr\alpha \u2208Aa$ denotes a trace over the atomic orbitals

*α*belonging to atom

*a*only, and

*ρ*

_{a}(

*t*) ≡

*P*

_{a}

*ρ*(

*t*)

*P*

_{a}is the projection of the time-dependent one-particle density operator,

on to the basis of atomic orbitals on atom *a*. We refer to **m**_{a}(*t*) as the exchange moment vector in order to differentiate it from the magnetic moment vector ** μ**. In terms of expansion coefficients [see Eq. (5)], the matrix representation of the density operator in the atomic orbital basis is

and the exchange moment of atom *a* is

From now on, for the sake of notational simplicity, we often suppress the time dependence.

The total energy due to the vector Stoner term is given as

where *I* is the Stoner exchange parameter given in units of energy. The TB representation of the self-consistent exchange potential experienced by molecular orbital *n* is a matrix Δ*H*_{ασ,α′σ′}, where

For *α* ∈ *A*_{a} this gives

The presence of the Stoner exchange term complicates the time-evolution algorithm described in Sec. II A (and makes it less than perfectly time reversible) because the Hamiltonian at time *t*, *H*(**m**_{a}(*t*), *t*), is also a function of the exchange moments at time *t*, which are not calculated in a reversible manner. The evolution of the exchange moments is found by determining the time derivative of the current moment vector **m**_{a}(*t*) by a backward difference, $m\u0307a(t)\u2248(ma(t)\u2212ma(t\u2212\delta t))/\delta t$, and Taylor expanding to the first order in *δt*,

The self-consistency cycle is run only once at the beginning of the calculation when determining the initial set of eigenstates, *ϕ*_{n}. For all subsequent time steps, **m**_{a} is calculated using Eq. (17).

### D. Model parameters

The TB model requires several experimental parameters. Reference 36 reports the value of the SOC parameter *ξ* for an oxygen dimer as *ξ* = 2.604 meV. The same paper gives an experimental value for the bond length *R*_{nuc} as 1.21 Å. This value is in agreement with *ab initio* results based on the relativistic Pauli-Breit Hamiltonian.^{37} The other tight-binding parameters required to describe an O_{2} dimer are the Stoner *I*, the on-site energy for *p* orbitals, *ϵ*_{p}, and the hopping parameters, $v$_{π} and $v$_{σ}. $v$_{π} is the hopping parameter between *p* orbitals perpendicular to the dimer axis, and $v$_{σ} is the hopping parameter between *p* orbitals parallel to the dimer axis. The on-site energy used in this work is *ϵ*_{p} = −16.77 eV and the hopping parameters are $v\pi =\u22120.63\u210f2meRnuc2$ and $v\sigma =2.22\u210f2meRnuc2$.^{38} Using the experimental bond length, *R*_{nuc} = 1.21 Å, gives $v$_{π} = −3.28 eV and $v$_{σ} = 11.55 eV. The value of the Stoner *I* parameter, *I* = 0.98 eV, is taken from Ref. 39.

### E. Calculating observables

Forces on the atoms are evaluated using the time-dependent equivalent of the Hellman-Feynman theorem.^{31,32} Let *a* ∈ {1, 2} index the two atoms in the dimer and **R**_{a} denote the position of the nucleus of atom *a*. Using the Hellman-Feynman theorem, we obtain

for the force on nucleus *a* due to the electrons, where we have used the definition of the density matrix given in Eq. (16). Since our method does not allow the dimer nuclei to move, the nuclei are effectively clamped in position by artificial externally applied forces which oppose any force enacted by the electrons or electromagnetic (EM) field. The only contributions to the Hamiltonian of Eq. (11) that depend on **R**_{1} and **R**_{2} are the hopping parameters. Using the distance-scaling of the hopping matrix elements discussed in Sec. II D and Slater-Koster tables^{40} for the *p* orbitals, $\u2207aH$ is calculated analytically.

All the other expectation values computed in this study are found by taking the trace of the operator multiplied by the density matrix, for example,

As shown in the Appendix, the classical nuclei experience both an internal torque, **Γ**_{int} [defined in Eq. (2)], exerted by the quantum mechanical electrons, and a direct torque, **Γ**^{N,EM}, exerted by the applied classical electromagnetic field. The total torque acting on the nuclei is the sum of these two contributions

The electromagnetic torque **Γ**^{N,EM} can be expressed as

where *N* is the set of nuclei, **r**_{a} is the position of nucleus *a*, and

is the Lorentz force on nucleus *a*. Here, *q*_{a} is the charge of nucleus *a*, **B**_{a} and **E**_{a} = −$\u2207a$*φ*_{a} − *∂*_{t}**A**_{a} are the applied magnetic and electric fields at the position of nucleus *a*, respectively, (*φ*_{a}/*c*, **A**_{a}) = (*φ*(**r**_{a}, *t*)/*c*, **A**(**r**_{a}, *t*)) is the electromagnetic four-potential experienced by nucleus *a*, $va=1ma(pa\u2212qaAa)$ is the nuclear velocity, *m*_{a} is the nuclear mass, and **p**_{a} is the canonical momentum of nucleus *a*.

We are interested in the dynamics of the nuclei, as these correspond to the motion of the lattice in a solid body. Since the nuclei are held fixed in position in our simulations, they do not experience a force due to the **v**_{a} × *B*_{a} terms in the Lorentz force. The direct EM torque, **Γ**^{N,EM}, thus arises solely from the circulating electric field produced by the time-varying magnetic field via Faraday’s law. The rate of change of the applied magnetic field in our simulations is 10 T per atomic time unit, implying a direct torque of $\u22128.9\xd710\u22124x^$ hartree a.u. This is significant, but our simulations last only 1000 a.u. of time, equivalent to approximately 2.4 × 10^{−14} s, while real experiments enact the change in field over times of the order of 10^{−2} s.^{14} The Faraday torque in our simulations is thus about 10^{12} times larger than it would be in a real experiment. Section III reports values of the torque on the nuclei due to the electrons, **Γ**_{int}. This quantity includes the effects of the Lorentz and Faraday forces on the electrons, which are built-in to the Hamiltonian in Eq. (11), but omits the torque arising from the direct action of the Faraday electric field on the nuclei. If required, the true mechanical torque on the nuclei can be evaluated using Eq. (24).

## III. RESULTS

To facilitate the interpretation of the physics through a gradual buildup in complexity, the discussion of the results is split into two parts. First, the results without SOC are presented, followed by the results with SOC. Without SOC, we discuss data obtained from simulations with (i) no *B* field (** B** =

**0**T), (ii) a constant

*B*field of 10

^{3}T ($B=103x^\u2009$ T), and (iii) a linear ramp in

*B*

_{x}from −5 × 10

^{3}T to 5 × 10

^{3}T over 10

^{3}hartree a.u. of time, described by the equation $B(t)=(\u22125\xd7103+10t)x^\u2009$ T, where

*t*is measured in atomic units. Such large magnetic field strengths are used because the magnetic Hamiltonian

*H*

_{B}= −

**·**

*μ***associated with a field of 1 T = 4.25 × 10**

*B*^{−6}a.u. is very small on the atomic energy scale. The decision to reverse the field by applying a linear ramp instead of rotating a

**B**vector of constant magnitude was made because the EdH experiment was performed in a magnetic field generated by a solenoid with an oscillatory current. At its center, the solenoid is only capable of producing a magnetic field pointing along its axis.

The experimental value of the spin-orbit coupling strength in oxygen, *ξ* = 2.604 meV, is so small that the simulation results with SOC are not visibly different from those without SOC. Thus, in order to observe the EdH effect clearly, it is necessary to use larger values for *ξ* than are physically realistic for O_{2}. In doing so, more can be learned about the effects of SOC on the dimer’s dynamics. We consider two different SOC strengths: (i) an intermediate coupling strength of *ξ* = 0.4 eV and (ii) a large coupling strength of *ξ* = 10^{3} eV. The value of *ξ* = 0.4 eV is large enough to have observable effects but not unreasonable for heavier atoms. (Values of *ξ* as large as 0.4 eV have been reported for Sr_{2}IrO_{4}.^{41} In bcc Fe, the value of *ξ* is close to 0.06 eV.^{42}) The unphysically large value of *ξ* = 10^{3} eV is chosen in order to investigate the limit in which SOC is the dominant energy scale.

Unless stated otherwise, the results below are expressed in hartree atomic units (a.u.).

### A. Without spin-orbit coupling

The results discussed in Subsection III C were all obtained in the absence of SOC, i.e., with *ξ* = 0 eV. The effects of exchange and the interaction with an external magnetic field are included. The simulations considered have (i) *B*_{x} = 0 T, (ii) *B*_{x} = 10^{3} T, and (iii) *B*_{x}(*t*) = (−5 × 10^{3} + 10*t*) T.

The *B*_{x} = 0 T simulation illustrates the behavior of the electrons in the absence of an external magnetic field. We find that there is some spin in each of the Cartesian directions. This is because, in the absence of SOC, the spin subsystem is decoupled from the geometry of the dimer which thus does not introduce a preferred direction. Furthermore, since *B* = 0 T, the energy is unaffected by the orientation of the spin; thus, the system is completely degenerate with respect to the spin direction. The direction of the spin, which is determined by the self-consistency cycle, is therefore the same as the direction of the initial random guess. This arbitrary dependence of the spin on the random direction of the initial guess is a necessary consequence of choosing a spin direction in an otherwise spherically symmetric system. The magnitude of the spin vector is unity, as expected for the triplet ground state of O_{2}.

The electrons have zero net orbital angular momentum because the *σ* and *π* subshells are filled, leaving one spin-up electron in each of the two *π*^{*} orbitals. The total orbital angular momentum is the sum of contributions from the two singly occupied *π*^{*} orbitals and is thus zero.

The orbital and spin angular momenta **L** and **S** are both constant in time because the Hamiltonian is independent of time and the system is in an energy eigenstate.

In the Appendix, we show that

where **Γ**_{dim} is the torque on the dimer and **J** = **L** + **S** is the total canonical electronic angular momentum in the Coulomb gauge. The torque on the dimer is therefore zero in this case since **B** = **0** T and **L** and **S** are constant.

A simulation of the dimer was also performed in the presence of a constant 10^{3} T external magnetic field. The magnetic field points in the $x^$ direction and the magnetic field term in the Hamiltonian (*H*_{B}) are minimized with ** μ** parallel to

**, so the spin lies antiparallel to the**

*B***B**field.

The orbital angular momentum, which is nonzero, also points opposite to the *B* field and remains constant in time again because the Hamiltonian is independent of time and the system is in an energy eigenstate. As is expected for this static system, Eq. (27) shows that the torque is zero.

Figure 3 shows the simulation results for the magnetic field profile $B(t)=Bx(t)x^=(\u22125\xd7103+10t)x^$. This corresponds to a linear ramp from −5 × 10^{3} T to 5 × 10^{3} T over a duration of 10^{3} a.u. The spin stays constant, pointing in the $x^$ direction throughout. The applied field *B*_{x} is negative at *t* = 0, and the spin is oriented antiparallel to this field in the initial ground state. The spin remains constant as it is initially aligned along $x^$ and *S*_{x} commutes with the Hamiltonian. The spin is unable to flip in response to the reversal of the applied magnetic field.

The evolution of the orbital angular momentum over time has two notable features. The first is an adiabatic effect in which the amount of orbital angular momentum pointing along the *B* field is proportional to minus the field. The second effect is a small sinusoidal variation, barely visible in Fig. 3(a), caused by a Rabi oscillation between eigenstates split by the *B* field. Section III B explains this effect by deriving the energy difference of the splitting and thus the frequency of the oscillations. A torque is exerted on the dimer by the change in the orbital angular momentum of the electrons. Although the torque is oscillatory, its time average is nonzero. The dimer would therefore start to rotate if its nuclei were not clamped in position. This torque is not related to the EdH effect since −*d*⟨**S**⟩/*dt* = **0** and there is no spin-orbit coupling.

### B. Oscillations induced by the time-dependent magnetic field

An analytic expression for the oscillation frequency of the orbital angular momentum can be found by considering the form of *H*_{B} = −** μ** ·

**in the MO basis. Let $(MT)n\alpha \sigma \u2261dn\alpha \sigma $ be the transformation matrix from the AO basis to the time-independent MO basis at**

*B**B*= 0 T. The expansion coefficients

*d*

_{nασ}are as defined in Eq. (4). The zero-field MOs

*ϕ*

_{n}are given by $\varphi n=(MT)n,\alpha \sigma \chi \alpha \sigma $, where

*χ*

_{ασ}is an AO and repeated suffices are summed. The matrix representations of the angular momentum operators transform as

where ** L** and

**are matrices in the atomic orbital basis and $L\u0303$ and $S\u0303$ are matrices in the**

*S**B*= 0 MO basis. The

*B*= 0 Hamiltonian

*H*

_{0}is diagonal in the MO basis, so only the Zeeman term −

**·**

*μ***= −**

*B**μ*

_{x}

*B*

_{x}mixes MOs when

*B*

_{x}≠ 0. By quantizing the spins along $x^$ rather than $z^$ so that $S\u0303x$ is diagonal, we ensure that the only off-diagonal contributions to −

*μ*

_{x}

*B*

_{x}∼ (

*L*

_{x}+ 2

*S*

_{x})

*B*

_{x}are those arising from the matrix $L\u0303x$.

In this form, the 12 × 12 Hamiltonian matrix becomes block diagonal, consisting of four uncoupled states (*π*_{x↑}, *π*_{x↓}, $\pi x\u2191*$, and $\pi x\u2193*$), and four 2 × 2 blocks, each of which couples one of the following pairs of MOs:

*σ*_{↑}and $\pi y\u2191*$,*σ*_{↓}and $\pi y\u2193*$,*π*_{y↑}and $\sigma \u2191*$,*π*_{y↓}and $\sigma \u2193*$.

For example, the 2 × 2 matrix describing the coupling between *σ*_{↑} and $\pi y\u2191*$ is

which has eigenvalues

The magnitude of the difference between the two eigenvalues is

This energy difference corresponds to a Rabi angular frequency *ω*_{1} = Δ*ϵ*_{1}/*ℏ* and hence to an oscillation in ⟨*L*_{x}⟩. The magnitude of Δ*ϵ*_{1} varies slowly as the applied B field varies.

To check that this explanation accounts for the observed oscillations in ⟨*L*_{x}⟩, consider the case when *B*_{x} = 5 × 10^{3} T, for which Δ*ϵ*_{1} = 0.545 52 hartree and *T* = 1/*f* = *h*/Δ*ϵ*_{1} = 11.52 a.u. This agrees well with the oscillation period seen in the simulation results in Fig. 3. The period associated with the other three pairs of coupled MOs is the same.

### C. With spin-orbit coupling

Since the results of the simulations without SOC are practically indistinguishable from results obtained using the small value of *ξ* appropriate for a real oxygen atom, two further simulations were carried out, one for an intermediate value of *ξ* and another for a very large value of *ξ* intended to approach the limit in which SOC is the dominant energy scale.

Figure 4 shows results from a simulation with *ξ* = 0.4 eV and a linear ramp in the *B* field. At the beginning of the simulation, the dynamics look similar to the SOC-free behavior shown in Fig. 3. The most noticeable difference is a slight decrease in the magnitude of the initial value of ⟨*L*_{x}⟩. Since *Ŝ*_{x} no longer commutes with the Hamiltonian when SOC is included, its expectation value can now change with time as demonstrated in Fig. 4(a). Toward the end of the simulation, the spin changes sign, which leads to the change in the direction of the torque shown in Fig. 4(b). Thus, unlike the SOC-free case, the spin is able to reverse the direction in response to the reversal of the applied field. The spin flip takes place near *t* = 800 a.u., well after the time, *t* = 500 a.u., at which the applied field passes through zero.

Figure 5 shows the eigenvalues of the time-evolved Hamiltonian, *H*(**m**_{a}(*t*), *t*), as a function of time. These “instantaneous” eigenvalues are calculated by diagonalizing *H*(**m**_{a}(*t*), *t*) non-self-consistently. Because *H*(**m**_{a}(*t*), *t*) depends on **m**_{a}(*t*), which evolved from previous time steps, the instantaneous eigenvalues and eigenstates are history dependent. Stoner exchange thus introduces a memory into the system. The energy level crossings near *t* = 800 a.u. determine the timing of the spin flip.

It is important to distinguish the instantaneous eigenvalues shown in Fig. 5 from the self-consistent eigenvalues at time *t* (not shown), which are obtained by carrying out a fully self-consistent diagonalization of the Hamiltonian in the presence of a fixed applied magnetic field *B*(*t*). Because the magnetic moments are always at their self-consistent values, the self-consistent eigenvalues are not history dependent. Neither the instantaneous eigenfunctions nor the self-consistent eigenfunctions are the same as the time-evolved molecular orbitals.

For *B*_{x} < 0 T, the lowest energy direction for **m**_{a} is to be aligned along $x^$ since this minimizes the spin contribution to −** μ** ·

**B**. As

*B*

_{x}(

*t*) rises through zero, the new minimum energy direction for

**m**

_{a}becomes $\u2212x^$. Thus, if a simulation were performed in which the system was solved self-consistently at each time step, an energy level crossing would occur at

*t*= 500 a.u., when the magnetic field

*B*

_{x}(

*t*) passed through zero.

In a time-evolved simulation, **m**_{a}(*t*) is calculated from its previous values and the $12Ima\u22c5\sigma \sigma \sigma \u2032$ vector Stoner exchange term in the effective Hamiltonian acts to keep **m**_{a}(*t*) aligned along its current value. A lower energy state with all spins reversed exists, but an energy barrier has to be surmounted to reach it. This causes **m**_{a}(*t*) to maintain its original direction for longer than in the corresponding self-consistent calculation, explaining the delay in the reversal of the spin of the time-evolved Hamiltonian shown in Fig. 5. Additional simulations show that the time at which the avoided crossing occurs moves toward *t* = 500 a.u. as the Stoner exchange parameter *I* is reduced and can be made to occur later by increasing *I*.

Comparing Figs. 5 and 6, we see that the reversal in the sign of the spin of the time-dependent *π* and *π*^{*} molecular orbitals does not occur at the same time as the avoided crossing, but slightly afterward. The adiabatic theorem says that the evolution of a time-dependent MO closely follows the evolution of the corresponding eigenstate of the Hamiltonian if the rate of change of *B* is small enough. The electron therefore remains predominantly in the adiabatically connected eigenstate of the time-evolved Hamiltonian, which changes its spin direction as the avoided crossing is traversed. Thus, the delay in the rotation of the spins of the time-evolved molecular orbitals as they respond to the change in the time-evolved Hamiltonian is due to the simulation being diabatic, which increases the tendency of the time-evolved MO to remain similar to its original instantaneous MO of the time-evolved Hamiltonian during the avoided crossing. This is supported by further simulations in which *Ḃ* is varied and it is found that for low *Ḃ*, the change in sign of the spin of the time-dependent *π* and *π*^{*} molecular orbitals tends toward being simultaneous with the change in sign of the spin of the *π* and *π*^{*} instantaneous eigenstates of the time-evolved Hamiltonian.

The simulation results for the case of extremely strong SOC (*ξ* = 10^{3} eV) are shown in Fig. 7. The (*ξ*/*ℏ*^{2})** L** ·

**spin-orbit term is now by far the largest in the Hamiltonian. Figure 7(a) shows that ⟨**

*S***L**⟩ and ⟨

**S**⟩ are proportional to each other, with a proportionality constant of 1/2.

This result can be understood by considering a Hamiltonian constructed solely from the SOC term. This decouples the two atoms since the SOC interaction acts on-site in our model. The Hamiltonian of one atom takes the form

where *ξ* = 10^{3} eV. The orbital energies are $\xi 2(\u2009j(\u2009j+1)\u2212l(l+1)\u2212s(s+1))=\xi 2(\u2009j(\u2009j+1)\u22122.75)$, with *j* = 1/2 or 3/2. The six eigenstates comprise a doublet with energy −10^{3} eV, corresponding to the *j* = 1/2 states with *m*_{j} = ±1/2, and a quadruplet with energy 500 eV, corresponding to the *j* = 3/2 states with *m*_{j} = 3/2, 1/2, −1/2, and −3/2. In the ground state of an oxygen atom (which contains 4 electrons in our *p*-band TB model), the two *j* = 1/2 orbitals are filled first, leaving the other two electrons to occupy the four available *j* = 3/2 orbitals. Using the Wigner-Eckart theorem,^{43} it can be shown that the valence electrons residing in the *j* = 3/2 states have $\varphi |S|\varphi =12\varphi |L|\varphi $ which explains the proportionality between ⟨**S**⟩ and ⟨**L**⟩ in Fig. 7(a).

Above we noted the link between the torque exerted by the electrons on the nuclei, **Γ**_{int}, and the rate of change of electronic angular momentum, ⟨**J**⟩. Equation (A6) of the Appendix shows that their vector components parallel to the *B* field are related by a minus sign. This link can be verified computationally by approximating −*d*⟨**J**⟩/*dt* using finite differences of the total electronic angular momentum ⟨**J**⟩ and comparing the result with the torque **Γ**_{int} acting on the nuclei.

The torque on the nuclei due to −*d*⟨**J**⟩/*dt* can be separated into orbital and spin contributions. The orbital contribution, −*d*⟨**L**⟩/*dt*, can be interpreted as in classical physics. The spin contribution, −*d*⟨**S**⟩/*dt*, corresponds to the EdH effect. In the absence of SOC, as was shown in Sec. III A, the spin direction cannot reverse when the applied field reverses by changing its magnitude along a fixed axis. In this case $B^\u22c5\Gamma int=\u2212B^\u22c5d\u27e8L\u27e9/dt$, implying that the EdH effect cannot occur without SOC. For intermediate SOC strengths (Fig. 4), we saw that the spin flip influences the direction of the torque. Both contributions to $B^\u22c5\Gamma int=\u2212B^\u22c5d\u27e8L\u27e9/dt\u2212B^\u22c5d\u27e8S\u27e9/dt$ contribute appreciably in this case. For extremely strong SOC strengths (Fig. 7), ⟨**L**⟩ and ⟨**S**⟩ lock together and its average over a short time window does not change significantly as a function of time. The initial magnitude of ⟨**S**⟩ is smaller than in the previous simulations, which is due to the *j* = 3/2 eigenfunctions of the SOC term having $\varphi |S|\varphi =12\varphi |L|\varphi $ as is described above.

## IV. CONCLUSIONS

This paper investigated a noncollinear TB model capable of simulating the EdH effect in a dimer. Since the EdH effect is based on the realignment of spins in a uniaxial *B* field, it cannot occur in the absence of SOC. Based on simulation results showing the torque on the dimer for three very different SOC strengths, we showed how the avoided crossings and the rate of change of magnetic field affect the reversal of electronic spins. We also showed that the EdH effect is not the only source of torque on a dimer in a time-varying *B* field: there is an additional contribution from the change in the electronic orbital angular momentum. The direct action on the nuclear charges of the Faraday electric field associated with the rate of change of the applied magnetic field also exerts a torque on the nuclei. This is very small for experimentally achievable values of *d***B**/*dt* (although significant in our simulations).

Future work on this topic could aim to better understand the mechanism that translates the microscopic EdH effect to the better known macroscopic observations. How is the EdH effect modified by the quenching of **L** in systems without rotational symmetry, such as the iron clusters investigated in Ref. 44? What are the consequences of irreversibility and the loss of energy into other microscopic degrees of freedom? Now that the dimer model has been established, it can be scaled up to treat larger assemblies of atoms and a more diverse range of atom types. Another potentially fruitful area of investigation would be to look for effective classical models capable of emulating the effects of SOC and yet simple enough to be used in large-scale simulations of real materials problems, including studies of radiation damage in ferromagnetic steels.

## ACKNOWLEDGMENTS

This work was supported through a studentship in the Centre for Doctoral Training on Theory and Simulation of Materials at Imperial College London, funded by EPSRC under Grant No. EP/L015579/1. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training program 2019–2020 under Grant Agreement No. 633053 and from the RCUK Energy Programme (Grant No. EP/P012450/1). W.M.C.F.’s visit to the University of Illinois at Urbana-Champaign, where much of his work on this project was done, was supported by the US Department of Energy under Grant No. NA 0002911. We acknowledge support from the Thomas Young Centre under Grant No. TYC-101.

### APPENDIX: RELATING THE DIMER TORQUE TO THE ELECTRONIC ANGULAR MOMENTUM

Below, we show how the total torque on the dimer nuclei is a sum of the direct torque due to the interaction of the nuclei with the electromagnetic (EM) field and the torque due to the Coulomb attraction between the electrons and the nuclei.

We consider a system of *n* particles (nuclei and electrons), all of which have spins and interact with the externally applied EM four-potential, (*φ*_{a}/*c*, **A**_{a}) = (*φ*(**r**_{a}, *t*)/*c*, **A**(**r**_{a}, *t*)), where *c* is the speed of light. The single-particle Hamiltonian of a particle in an EM field,^{45} generalized to multiple interacting particles is

where particle *a* has mass *m*_{a}, charge *q*_{a}, and SOC parameter *ξ*_{a}, and *V* (*r*) is the interaction potential energy for a pair of particles separated by a distance *r*.

By calculating the rate of change of the kinetic momentum *p*_{a} − *q*_{a}*A*_{a}, we obtain the Lorentz force operator for particle *a*,

where *E*_{a} = −$\u2207a$*φ*_{a} − *∂*_{t}**A**_{a} is the electric field operator on particle *a* and $va=1ma(pa\u2212qaAa)$ is the velocity operator for particle *a*. Since we are working in the rest frame of the molecule, the nuclei have no orbital angular momentum and experience no spin-orbit coupling. Equation (A2) was derived quantum mechanically but reduces to the familiar classical Lorentz force law if the nuclei are treated classically, as is the case in this work. The direct electromagnetic torque on the classical nuclei is $\Gamma N,EM=\u2211a\u2208Nra\u2009\xd7\u2009FaEM$, where **r**_{a} is the classical position of nucleus *a* and $FaEM$ is the direct electromagnetic force on nucleus *a*. The total torque exerted on the set *N* of nuclei in the molecule is given by

where **Γ**_{int} = ⟨*∑*_{a∈N,b∈E}**r**_{a} × **F**_{ab}⟩ is the Hellmann-Feynman torque exerted on the nuclei by the set *E* of electrons and **F**_{ab} = −$\u2207a$*V*$$|**r**_{a} − **r**_{b}|$$ is the operator for the Hellman-Feynman force on nucleus *a* due to electron *b*.

Let **J** = ∑_{a∈E} **J**_{a} be the total electronic canonical angular momentum, where **J**_{a} = **L**_{a} + **S**_{a} is the sum of the canonical orbital and spin angular momenta of electron *a*. The rate of change of the expectation value of **J** can be calculated using the Ehrenfest theorem,

Simplifying the result by assuming that the applied *B* field is uniform, working in the Coulomb gauge, and neglecting terms of second order in **A** yield

and thus,

where $B^$ denotes the unit vector pointing in the direction of the *B* field.