The microscopic Einstein-de Haas effect

TL;DR

This paper models the microscopic Einstein-de Haas effect using a non-collinear tight-binding approach, revealing the roles of spin-orbit coupling, magnetic fields, and electronic orbital contributions in transferring spin to mechanical angular momentum.

Contribution

It introduces a detailed microscopic model of the Einstein-de Haas effect incorporating spin-orbit coupling and external magnetic fields, advancing understanding beyond semi-classical explanations.

Findings

Validation of the EdH effect through generated torque on the dimer.

Identification of additional torque contributions from orbital angular momentum and Faraday electric fields.

Analysis of energy level crossings and magnetic field variation effects.

Abstract

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 semi-classical explanation of the effect exists, modern electronic structure methods have not yet been applied to modelling the phenomenon. In this paper we investigate its microscopic origins by means of a non-collinear tight-binding model of an $\textrm{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. Avoided energy level crossings and the rate of change of magnetic field determine the evolution of the spin. We find also that the torque exerted on the nuclei by the electrons in a time-varying $B$ field is not only due to the EdH effect. 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.

Equations

$$\Delta \mathbf{M}=-g^{\prime }\frac{e}{2m_e}\Delta \mathbf{J},$$

$${\mathbf{\Gamma }}_{\text{int}}=R_{\text{nuc}}(F_{1,y}\hat{\mathbf{x}}-F_{1,x}\hat{\mathbf{y}})=-R_{\text{nuc}}(F_{2,y}\hat{\mathbf{x}}-F_{2,x}\hat{\mathbf{y}}),$$

$$p_{x,\uparrow },p_{x,\downarrow },p_{y,\uparrow },p_{y,\downarrow },p_{z,\uparrow },p_{z,\downarrow }.$$

$${\varphi }_n=\sum _{\alpha \sigma }{d}_{n\alpha \sigma }{\chi }_{\alpha \sigma },$$

$${\sigma }_{\uparrow },{\sigma }_{\downarrow },{\pi }_{x,\uparrow },{\pi }_{x,\downarrow },{\pi }_{y,\uparrow },{\pi }_{y,\downarrow },$$

$${\pi }_{x,\uparrow }^{\ast },{\pi }_{x,\downarrow }^{\ast },{\pi }_{y,\uparrow }^{\ast },{\pi }_{y,\downarrow }^{\ast },{\sigma }_{\uparrow }^{\ast },{\sigma }_{\downarrow }^{\ast },$$

$${\psi }_n(t)=\sum _{\alpha \sigma }{d}_{n\alpha \sigma }(t){\chi }_{\alpha \sigma }$$

$$i\hslash \frac{\partial }{\partial t}{d}_{n\alpha \sigma }(t)=\sum _{{\alpha }^{\prime }{\sigma }^{\prime }}{H}_{\alpha \sigma ,{\alpha }^{\prime }{\sigma }^{\prime }}(t)d_{n{\alpha }^{\prime }{\sigma }^{\prime }}(t).$$

$$i\hslash \frac{\partial D(t)}{\partial t}=H(t)D(t).$$

$$D(t+dt)=\exp \left(\frac{H(t+\frac{1}{2}dt)}{i\hslash }dt\right)D(t).$$

$$d_{n\alpha \sigma }(t+\delta t)=\sum _{{\alpha }^{\prime }{\sigma }^{\prime }}{\left(e^{H(t+\frac{1}{2}\delta t)\delta t/i\hslash }\right)}_{\alpha \sigma ,{\alpha }^{\prime }{\sigma }^{\prime }}{d}_{n{\alpha }^{\prime }{\sigma }^{\prime }}(t).$$

$$d_{n\alpha \sigma }(0)=⟨{\chi }_{\alpha \sigma }||{\varphi }_n⟩.$$

$$H=H_0+H_B+H_{\text{SOC}}+H_{\text{ex}},$$

$$\mathbf{\mu }=-\frac{{\mu }_B}{\hslash }(\mathbf{L}+2\mathbf{S})$$

$$H_{\text{SOC}}=\frac{1}{2m_e^2{c}^2}\frac{1}{r}\frac{dV(t)}{dr}\mathbf{L}\cdot \mathbf{S},$$

$$H_{\text{SOC}}\approx \frac{\xi }{{\hslash }^2}\mathbf{L}\cdot \mathbf{S}.$$

$${\mathbf{m}}_a(t)={\text{tr}}_{\alpha \in A_a}[\rho (t)\mathbf{\sigma }]=\text{tr}[{\rho }_a(t)\mathbf{\sigma }],$$

$$\rho (t)=\sum _{n\text{occ}}|{\psi }_n(t)⟩⟨{\psi }_n(t)|$$

$${\rho }_{{\alpha }^{\prime }{\sigma }^{\prime },\alpha \sigma }(t)=\sum _{n\text{occ}}{d}_{n{\alpha }^{\prime }{\sigma }^{\prime }}^{}(t)d_{n\alpha \sigma }^{\ast }(t)$$

$${\mathbf{m}}_a(t)$$

$$\Delta U=-\frac{1}{4}I\sum _a{\mathbf{m}}_a\cdot {\mathbf{m}}_a,$$

$$\frac{\partial \Delta U}{\partial d_{n\alpha \sigma }^{\ast }}=\sum _{{\alpha }^{\prime }{\sigma }^{\prime }}\Delta H_{\alpha \sigma ,{\alpha }^{\prime }{\sigma }^{\prime }}{d}_{n{\alpha }^{\prime }{\sigma }^{\prime }}.$$

$$\Delta H_{\alpha \sigma ,{\alpha }^{\prime }{\sigma }^{\prime }}=-{\delta }_{\alpha {\alpha }^{\prime }}\frac{1}{2}I{\mathbf{m}}_a\cdot {\mathbf{\sigma }}_{\sigma {\sigma }^{\prime }}.$$

$${\mathbf{m}}_a(t+\frac{1}{2}\delta t)\approx {\mathbf{m}}_a(t)+\frac{1}{2}{\dot{\mathbf{m}}}_a(t)\delta t.$$

$${\mathbf{F}}_a$$

$$\text{expectationvalue}\mathbf{L}=\text{tr}(\rho \mathbf{L}),\text{expectationvalue}\mathbf{S}=\text{tr}(\rho \mathbf{S}),\text{expectationvalue}\mathbf{\mu }=\text{tr}(\rho \mathbf{\mu }).$$

$${\mathbf{\Gamma }}^N$$

$${\mathbf{\Gamma }}^{N,EM}=\sum _{a\in N}{\mathbf{r}}_a\text{crossproduct}{\mathbf{F}}_a^{EM},$$

$${\mathbf{F}}_a^{EM}=q_a({\mathbf{v}}_a\times {\mathbf{B}}_a)+q_a{\mathbf{E}}_a$$

$$\frac{d\text{expectationvalue}\mathbf{J}}{dt}=-{\mathbf{\Gamma }}_{\text{int}}-\mathbf{B}\text{crossproduct}\text{expectationvalue}\mathbf{\mu },$$

$${\tilde{\mathbf{L}}}_{mn}$$

Full text

††preprint: APS/123-QED

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 $\sim$10 T. In addition to stresses that such fields generate in the reactor components Gade et al. [2015], they may also affect the dislocation microstructure and plasticity of materials Alshits et al. [2003] as well as diffusion of radiation defects. Current tokamak designs, including ITER, use austenitic Fe-Cr-Ni steels Kawashima et al. [2001] that are non-magnetic on the macroscopic scale while being microscopically antiferromagnetic Wróbel et al. [2015]. In a demonstration fusion reactor (DEMO) Kohyama et al. [1998], 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 VanVleck [1937]; Daniel et al. [2008]; L. Daniel [2018] suggests that atomic magnetic moments are coupled to the lattice, giving rise to a magnetic contribution to electric and thermal resistivity Pepperhoff and Acet [2001] and to spin-electron-lattice relaxation effects in collision cascades Ma et al. [2012].

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 O. W. Richardson and has been demonstrated numerous times since then Richardson [1908]; Barnett [1915]; Kittel [1949]. In Einstein and de Haas’s Einstein and de Haas [1915] 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^{\prime}$, relates the change in magnetization $\Delta\mathbf{{M}}$ to a change in electronic angular momentum $\Delta\mathbf{{J}}$ according to

[TABLE]

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^{\prime}$ for ferromagnetic materials Scott [1962].

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 Zhang and Hübner [2000]. 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 Ohta and Shimizu [1982]; Ma et al. [2012].

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 favourable position, while magnetic DFT calculations predict that the tetrahedral position is lower in energy Seletskaia et al. [2005]. Self-interstitial defects provide another good example: a self-interstitial atom in magnetic bcc iron adopts a $\langle 110\rangle$ dumbbell structure Fu et al. [2004], which is different from the linear $\langle 111\rangle$ crowdion configuration adopted by self-interstitials in non-magnetic bcc transition metals Derlet et al. [2007]; Nguyen-Manh et al. [2006]. Perhaps even more significantly, magnetic fluctuations give rise to a sequence of $\alpha$-$\gamma$-$\delta$ bcc-fcc-bcc phase transitions in iron and steels at elevated temperature Yu. Lavrentiev et al. [2010]; Ma et al. [2017] that require computing free energy differences between competing magnetic phases with meV 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 Iota et al. [2007]. Lack of understanding also hinders attempts to construct microscopic models of heat flow in ferromagnetic materials Pepperhoff and Acet [2001]; Dudarev [2013] suitable for use in large-scale atomistic or semi-classical simulations. Our work represents a first step towards 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 Górecki and Rza̧żewski [2016]. A recent study of a single-molecule magnet coupled to a nanomechanical resonator constitutes an experimental realisation of the EdH effect for a single molecule Ganzhorn et al. [2016].

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 Paxton and Finnis [2008]. 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 Soin et al. [2011].

In this paper we study the EdH effect for an O2 dimer subject to a time-varying magnetic field. We use a non-collinear TB model that includes coupling to the external magnetic $B$ field, SOC, and vector Stoner exchange Coury et al. [2016]. Forces on the atoms are calculated using the Hellman-Feynman theorem Hellmann [1937]; Feynman [1939]. 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 $\mathbf{{B}}$ commutes with the Hamiltonian and $\mathbf{{S}}\cdot\hat{\mathbf{B}}$ 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.

This paper is structured as follows. Section II describes the method used for the calculations. The results of the simulations are discussed in section III and conclusions are drawn in section IV.

II Theory

Fig. 1 shows a schematic diagram of the O2 dimer studied in this paper. The dimer axis $\hat{\bm{\zeta}}$ is held fixed along the $\hat{\mathbf{z}}$ direction. The internuclear separation is $R_{\textrm{nuc}}$. The magnetic field is aligned along $\hat{\mathbf{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

[TABLE]

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

To create a TB model of O2 it is helpful to know some of its electronic properties. Of the 16 electrons in an O2 molecule, 4 fill the 1s states, 4 fill the 2s states, and 8 occupy bonding and anti-bonding molecular orbitals (MOs) made from the $2p$ atomic orbitals (AOs). A diagram of the electronic structure of the $2p$ 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

[TABLE]

The basis set does not include the $1s$ and $2s$ orbitals below the $2p$ shell or any orbitals above it, so the dimer Hamiltonian is a $12\times 12$ Hermitian matrix and the 12 molecular orbitals (MO) are occupied by 8 electrons. The MOs $\phi_{n}$ are obtained by diagonalizing the Hamiltonian matrix and can be written as linear combinations of the basis functions $\chi_{\alpha\sigma}$ with expansion coefficients $d_{n\alpha\sigma}$:

[TABLE]

where $\alpha$ runs over the spatial AOs on both atoms and $\sigma$ 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 twelve $2p$ MOs are labeled as follows:

[TABLE]

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

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

II.1 Time Evolution

The initial molecular orbitals $\phi_{n}$ obtained by diagonalizing the Hamiltonian matrix at the beginning of a simulation evolve into time-dependent molecular orbitals $\psi_{n}(t)$ according to the time-dependent Schrödinger equation, $i\hbar\partial_{t}\psi_{n}(t)=H(t)\psi_{n}(t)$, subject to the initial condition $\psi_{n}(t=0)=\phi_{n}$. The Hamiltonian $H(t)$ depends on time if the applied $B$ field depends on time, so $\psi_{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 $d_{n\alpha}(t)$ are defined by

[TABLE]

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

[TABLE]

The method used to evolve this equation in time is time-reversible and unitary. Rewriting the equation of motion, Eq. (6), with $(D)_{n\alpha\sigma}=d_{n\alpha\sigma}$ and $(H)_{\alpha^{\prime}\sigma^{\prime},\alpha\sigma}=H_{\alpha^{\prime}\sigma^{\prime},\alpha\sigma}$ as matrices, we obtain

[TABLE]

Letting $dt$ be an infinitesimal positive time interval gives Hatano and Suzuki [2005]

[TABLE]

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

[TABLE]

The initial values of the expansion coefficients of the time-dependent MO $\psi_{n}(t)$ are obtained from the corresponding eigenfunction $\phi_{n}$ of the initial Hamiltonian $H(t=0)$. Since $\psi_{n}(0)=\phi_{n}$, we have:

[TABLE]

The time step used in all simulations is 4 Hartree atomic units, which corresponds to approximately $97\text{\,}\mathrm{a}\mathrm{s}$.

II.2 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 O2 it is also necessary to include exchange interactions. A non-collinear TB method can accommodate all of these requirements.

The Hamiltonian employed is

[TABLE]

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_{\textrm{SOC}}$ is the SOC term, and $H_{\textrm{ex}}$ is the exchange term.

$H_{B}$ is taken from the standard Pauli Hamiltonian and is given by $-\bm{\mu}\cdot\bm{B}(t)$, where

[TABLE]

is the total magnetic moment, $\mathbf{{L}}$ and $\mathbf{{S}}$ are the canonical orbital and spin angular momentum operators in the Coulomb gauge, $\mu_{B}$ is the Bohr magneton, and $\bm{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 Foldy and Wouthuysen [1950], which, for spherical potentials, gives

[TABLE]

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, $\xi$, so that Landau and Lifshitz [1965]

[TABLE]

This is the form of $H_{\textrm{SOC}}$ used in this work. The next subsection describes the Stoner exchange term, $H_{\textrm{ex}}$.

II.3 Vector Stoner Exchange

All of the results reported in this paper have the effects of the vector Stoner exchange term ($H_{\textrm{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 Coury et al. [2016]. 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 $\mathbf{{m}}_{a}(t)$ on atom $a$ $(\in\{1,2\})$ at time $t$ is evaluated as

[TABLE]

where $\bm{\sigma}$ is the vector of Pauli spin matrices, $\textrm{tr}_{\alpha\in A_{a}}$ denotes a trace over the atomic orbitals $\alpha$ belonging to atom $a$ only, and $\rho_{a}(t)\equiv P_{a}\rho(t)P_{a}$ is the projection of the time-dependent one-particle density operator

[TABLE]

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

[TABLE]

and the exchange moment of atom $a$ is

[TABLE]

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

[TABLE]

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 $\Delta H_{\alpha\sigma,\alpha^{\prime}\sigma^{\prime}}$, where

[TABLE]

For $\alpha\in A_{a}$ this gives,

[TABLE]

The presence of the Stoner exchange term complicates the time-evolution algorithm described in Sec. II.1 (and makes it less than perfectly time reversible) because the Hamiltonian at time $t$, $H(\mathbf{{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 $\mathbf{{m}}_{a}(t)$ by a backward difference, $\dot{\mathbf{{m}}}_{a}(t)\approx(\mathbf{{m}}_{a}(t)-\mathbf{{m}}_{a}(t-\delta t))/\delta t$, and Taylor expanding to first order in $\delta t$:

[TABLE]

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

II.4 Model Parameters

The TB model requires several experimental parameters. Reference Tinkham and Strandberg [1955] reports the value of the SOC parameter $\xi$ for an oxygen dimer as $\xi=2.604~{}\textrm{meV}$. The same paper gives an experimental value for the bond length $R_{\textrm{nuc}}$ as $1.21\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. This value is in agreement with ab initio results based on the relativistic Pauli-Breit Hamiltonian Fedorov et al. [2001]. The other tight-binding parameters required to describe an O2 dimer are the Stoner $I$, the on-site energy for $p$ orbitals, $\epsilon_{p}$, and the hopping parameters, $v_{\pi}$ and $v_{\sigma}$. $v_{\pi}$ is the hopping parameter between $p$ orbitals perpendicular to the dimer axis and $v_{\sigma}$ the hopping parameter between $p$ orbitals parallel to the dimer axis. The on-site energy used in this work is $\epsilon_{p}=-16.77$ eV and the hopping parameters are $v_{\pi}=-0.63\frac{\hbar^{2}}{m_{e}R_{\textrm{nuc}}^{2}}$ and $v_{\sigma}=2.22\frac{\hbar^{2}}{m_{e}R_{\textrm{nuc}}^{2}}$ Harrison [2004]. Using the experimental bond length, $R_{\textrm{nuc}}=$ $1.21\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, gives $v_{\pi}=-3.28$ eV and $v_{\sigma}=11.55$ eV. The value of the Stoner $I$ parameter, $I=$0.98\text{,}\mathrm{e}\mathrm{V}$$, is taken from Ref. Schweitzer and Schmidt [2003].

II.5 Calculating Observables

Forces on the atoms are evaluated using the time-dependent equivalent of the Hellman-Feynman theorem Hellmann [1937]; Feynman [1939]. Let $a\in\{1,2\}$ index the two atoms in the dimer and $\mathbf{{R}}_{a}$ denote the position of the nucleus of atom $a$. Using the Hellman-Feynman theorem, we obtain

[TABLE]

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 EM field. The only contributions to the Hamiltonian of Eq. (11) that depend on ${\mathbf{{R}}}_{1}$ and ${\mathbf{{R}}}_{2}$ are the hopping parameters. Using the distance-scaling of the hopping matrix elements discussed in Sec. II.4 and Slater-Koster tables Sutton [1993] for the $p$ orbitals, $\bm{\nabla}_{a}H$ 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,

[TABLE]

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

[TABLE]

The electromagnetic torque $\mathbf{{\Gamma}}^{N,EM}$ can be expressed as

[TABLE]

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

[TABLE]

is the Lorentz force on nucleus $a$. Here $q_{a}$ is the charge of nucleus $a$, $\mathbf{{B}}_{a}$ and $\mathbf{{E}}_{a}=-\bm{\nabla}_{a}\varphi_{a}-\partial_{t}\mathbf{{A}}_{a}$ are the applied magnetic and electric fields at the position of nucleus $a$, $(\varphi_{a}/c,\mathbf{{A}}_{a})=(\varphi(\mathbf{{r}}_{a},t)/c,\mathbf{{A}}(\mathbf{{r}}_{a},t))$ is the electromagnetic four-potential experienced by nucleus $a$, ${\mathbf{{v}}_{a}=\frac{1}{m_{a}}(\mathbf{{p}}_{a}-q_{a}\mathbf{{A}}_{a})}$ is the nuclear velocity, $m_{a}$ is the nuclear mass, and $\mathbf{{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 $\mathbf{{v}}_{a}\times\bm{B}_{a}$ terms in the Lorentz force. The direct EM torque, $\mathbf{{\Gamma}}^{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 $-8.9\times 10^{-4}\hat{\mathbf{x}}$ Hartree atomic units. This is significant, but our simulations last only 1000 atomic units of time, equivalent to approximately $2.4\text{\times}{10}^{-14}\text{\,}\mathrm{s}$, while real experiments enact the change in field over times of the order of ${10}^{-2}\text{\,}\mathrm{s}$ Einstein and de Haas [1915]. 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, $\mathbf{{\Gamma}}_{\textrm{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 build up 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 ($\bm{B}=\bm{0}\,$T), (ii) a constant $B$ field of ${10}^{3}\text{\,}\mathrm{T}$ ($\bm{B}=${10}^{3}$\hat{\bm{x}}\,$T), and (iii) a linear ramp in $B_{x}$ from $-5\text{\times}{10}^{3}\text{\,}\mathrm{T}$ to $5\text{\times}{10}^{3}\text{\,}\mathrm{T}$ over ${10}^{3}\text{\,}\mathrm{H}$artree atomic units of time, described by the equation $\bm{B}(t)=(-$5\text{\times}{10}^{3}$+10t)\hat{\bm{x}}\,$T, where $t$ is measured in atomic units. Such large magnetic fields strengths are used because the magnetic Hamiltonian $H_{B}=-\bm{\mu}\cdot\bm{B}$ associated with a field of 1\text{\,}\mathrm{T}$=$4.25\text{\times}{10}^{-6} atomic units is very small on the atomic energy scale. The decision to reverse the field by applying a linear ramp instead of rotating a $\mathbf{{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, ${\xi=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 $\xi$ than are physically realistic for O2. 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 {\xi=0.4\text{,}\mathrm{e}\mathrm{V}}; and (ii) a large coupling strength of {\xi={10}^{3}\text{,}\mathrm{e}\mathrm{V}}. The value of {\xi=0.4\text{,}\mathrm{e}\mathrm{V}} is large enough to have observable effects but not unreasonable for heavier atoms. (Values of $\xi$ as large as $0.4\text{\,}\mathrm{e}\mathrm{V}$ have been reported for Sr2IrO4 Kim et al. [2008]. In bcc Fe the value of $\xi$ is close to 0.06 eV Desjonquères et al. [2007].) The unphysically large value of {\xi={10}^{3}\text{,}\mathrm{e}\mathrm{V}} 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.).

III.1 Without spin-orbit coupling

The results discussed in this section were all obtained in the absence of SOC, i.e., with ${\xi=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\text{\times}{10}^{3}+10t)}\,T.

The $B_{x}=0\,$T simulation illustrates the behaviour 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 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 O2.

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

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

In appendix A we show that

[TABLE]

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

A simulation of the dimer was also performed in the presence of a constant ${10}^{3}\text{\,}\mathrm{T}$ external magnetic field. The magnetic field points in the $\hat{\mathbf{x}}$ direction and the magnetic field term in the Hamiltonian ($H_{B}$) is minimized with $\bm{\mu}$ parallel to $\bm{B}$, so the spin lies anti-parallel to the $\mathbf{{B}}$ field.

The orbital angular momentum, which is non-zero, 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.

Fig. 3 shows the simulation results for the magnetic field profile $\mathbf{{B}}(t)=B_{x}(t)\hat{\mathbf{x}}=($-5\text{\times}{10}^{3}$+10t)\hat{\mathbf{x}}$. This corresponds to a linear ramp from $-5\text{\times}{10}^{3}\text{\,}\mathrm{T}$ to $5\text{\times}{10}^{3}\text{\,}\mathrm{T}$ over a duration of ${10}^{3}\text{\,}\mathrm{a}$.u. The spin stays constant, pointing in the $\hat{\mathbf{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 $\hat{\mathbf{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.2 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 orbital angular momentum of the electrons. Although the torque is oscillatory, its time average is non-zero. 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\expectationvalue{\mathbf{{S}}}}/{dt}=\mathbf{{0}}$ and there is no spin-orbit coupling.

III.2 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}=-\bm{\mu}\cdot\bm{B}$ in the MO basis. Let $(M^{T})_{n\alpha\sigma}\equiv d_{n\alpha\sigma}$ be the transformation matrix from the AO basis to the time-independent MO basis at $B=0$ T. The expansion coefficients $d_{n\alpha\sigma}$ are as defined in Eq. (4). The zero-field MOs $\phi_{n}$ are given by $\phi_{n}=(M^{T})_{n,\alpha\sigma}\chi_{\alpha\sigma}$, where $\chi_{\alpha\sigma}$ is an AO and repeated suffices are summed. The matrix representations of the angular momentum operators transform as

[TABLE]

where $\bm{L}$ and $\bm{S}$ are matrices in the atomic orbital basis and $\tilde{\bm{L}}$ and $\tilde{\bm{S}}$ are matrices in the $B=0$ MO basis. The $B=0$ Hamiltonian $H_{0}$ is diagonal in the MO basis, so only the Zeeman term $-\bm{\mu}\cdot\bm{B}=-\mu_{x}B_{x}$ mixes MOs when $B_{x}\neq 0$. By quantizing the spins along $\hat{\mathbf{x}}$ rather than $\hat{\mathbf{z}}$, so that $\tilde{S}_{x}$ is diagonal, we ensure that the only off-diagonal contributions to $-\mu_{x}B_{x}\sim(L_{x}+2S_{x})B_{x}$ are those arising from the matrix $\tilde{L}_{x}$.

In this form the $12\times 12$ Hamiltonian matrix becomes block diagonal, consisting of four uncoupled states ($\pi_{x\uparrow}$, $\pi_{x\downarrow}$, $\pi^{*}_{x\uparrow}$ and $\pi^{*}_{x\downarrow}$), and four $2\times 2$ blocks, each of which couples one of the following pairs of MOs:

$\sigma_{\uparrow}$ and $\pi^{*}_{y\uparrow}$, 2. 2.

$\sigma_{\downarrow}$ and $\pi^{*}_{y\downarrow}$, 3. 3.

$\pi_{y\uparrow}$ and $\sigma^{*}_{\uparrow}$, 4. 4.

$\pi_{y\downarrow}$ and $\sigma^{*}_{\downarrow}$.

For example, the $2\times 2$ matrix describing the coupling between $\sigma_{\uparrow}$ and $\pi^{*}_{y\uparrow}$ is

[TABLE]

which has eigenvalues

[TABLE]

The magnitude of the difference between the two eigenvalues is

[TABLE]

This energy difference corresponds to a Rabi angular frequency $\omega_{1}=\Delta\epsilon_{1}/\hbar$ and hence to an oscillation in $\expectationvalue{L_{x}}$. The magnitude of $\Delta\epsilon_{1}$ varies slowly as the applied B field varies.

To check that this explanation accounts for the observed oscillations in $\expectationvalue{L_{x}}$, consider the case when {B_{x}=5\text{\times}{10}^{3}} T, for which $\Delta\epsilon_{1}=0.54552$ Ha and $T=1/f=h/\Delta\epsilon_{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.

III.3 With spin-orbit coupling

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

Fig. 4 shows results from a simulation with ${\xi=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 behaviour shown in Fig. 3. The most noticeable difference is a slight decrease in the magnitude of the initial value of $\expectationvalue{L_{x}}$. Since $\hat{S}_{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). Towards the end of the simulation the spin changes sign, which leads to the change in direction of the torque shown in Fig. 4(b). Thus, unlike the SOC-free case, the spin is able to reverse 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(\mathbf{{m}}_{a}(t),t)$, as a function of time. These ‘instantaneous’ eigenvalues are calculated by diagonalizing $H(\mathbf{{m}}_{a}(t),t)$ non-self-consistently. Because $H(\mathbf{{m}}_{a}(t),t)$ depends on $\mathbf{{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 $\mathbf{{m}}_{a}$ is to be aligned along $\hat{\mathbf{x}}$ since this minimizes the spin contribution to $-\bm{\mu}\cdot\mathbf{{B}}$. As $B_{x}(t)$ rises through zero, the new minimum energy direction for $\mathbf{{m}}_{a}$ becomes $-\hat{\mathbf{x}}$. Thus, if a simulation were performed in which the system were 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, $\mathbf{{m}}_{a}(t)$ is calculated from its previous values and the $\frac{1}{2}I\mathbf{{m}}_{a}\cdot\bm{\sigma}_{\sigma\sigma^{\prime}}$ vector Stoner exchange term in the effective Hamiltonian acts to keep $\mathbf{{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 $\mathbf{{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 towards $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 $\pi$ and $\pi^{\ast}$ molecular orbitals does not occur at the same time as the avoided crossing, but slightly afterwards. 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 $\dot{B}$ is varied and it is found that for low $\dot{B}$, the change in sign of the spin of the time-dependent $\pi$ and $\pi^{\ast}$ molecular orbitals tends toward being simultaneous with the change in sign of the spin of the $\pi$ and $\pi^{\ast}$ instantaneous eigenstates of the time-evolved Hamiltonian.

The simulation results for the case of extremely strong SOC ($\xi=10^{3}$ eV) are shown in Fig. 7. The $(\xi/\hbar^{2})\bm{L}\cdot\bm{S}$ spin-orbit term is now by far the largest in the Hamiltonian. Figure 7(a) shows that $\expectationvalue{\mathbf{{L}}}$ and $\expectationvalue{\mathbf{{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

[TABLE]

where {\xi={10}^{3}} eV. The orbital energies are $\frac{\xi}{2}(j(j+1)-l(l+1)-s(s+1))=\frac{\xi}{2}(j(j+1)-2.75)$, with $j=1/2$ or $3/2$. The six eigenstates comprise a doublet with energy $-{10}^{3}\text{\,}\mathrm{e}\mathrm{V}$, corresponding to the $j=1/2$ states with $m_{j}=\pm 1/2$, and a quadruplet with energy $500\text{\,}\mathrm{e}\mathrm{V}$, 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 Cohen-Tannoudji et al. [1986], it can be shown that the valence electrons residing in the $j=3/2$ states have $\bra{\phi}\mathbf{{S}}\ket{\phi}=\frac{1}{2}\bra{\phi}\mathbf{{L}}\ket{\phi}$ which explains the proportionality between $\expectationvalue{\mathbf{{S}}}$ and $\expectationvalue{\mathbf{{L}}}$ in Fig. 7(a).

Above we noted the link between the torque exerted by the electrons on the nuclei, $\bm{\Gamma}_{\text{int}}$, and the rate of change of electronic angular momentum, $\expectationvalue{\mathbf{{J}}}$. Eq. (39) of appendix A 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\expectationvalue{\mathbf{{J}}}}/dt$ using finite differences of the total electronic angular momentum $\expectationvalue{\mathbf{{J}}}$ and comparing the result with the torque $\bm{\Gamma}_{\text{int}}$ acting on the nuclei.

The torque on the nuclei due to $-{d\expectationvalue{\mathbf{{J}}}}/{dt}$ can be separated into orbital and spin contributions. The orbital contribution, $-{d\expectationvalue{\mathbf{{L}}}}/{dt}$, can be interpreted as in classical physics. The spin contribution, $-{d\expectationvalue{\mathbf{{S}}}}/{dt}$, corresponds to the EdH effect. In the absence of SOC, as was shown in section III.1, the spin direction cannot reverse when the applied field reverses by changing its magnitude along a fixed axis. In this case $\hat{\mathbf{B}}\cdot\bm{\Gamma}_{\textrm{int}}=-\hat{\mathbf{B}}\cdot{d\expectationvalue{\mathbf{{L}}}}/{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 $\hat{\mathbf{B}}\cdot\bm{\Gamma}_{\textrm{int}}=-\hat{\mathbf{B}}\cdot{d\expectationvalue{\mathbf{{L}}}}/{dt}-\hat{\mathbf{B}}\cdot{d\expectationvalue{\mathbf{{S}}}}/{dt}$ contribute appreciably in this case. For extremely strong SOC strengths (Fig. 7), $\expectationvalue{\mathbf{{L}}}$ and $\expectationvalue{\mathbf{{S}}}$ lock together and its average over a short time window does not change significantly as a function of time. The initial magnitude of $\expectationvalue{\mathbf{{S}}}$ is smaller than in the previous simulations, which is due to the $j=3/2$ eigenfunctions of the SOC term having $\bra{\phi}\mathbf{{S}}\ket{\phi}=\frac{1}{2}\bra{\phi}\mathbf{{L}}\ket{\phi}$ as is described above.

IV Conclusions

This paper investigated a non-collinear 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 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\mathbf{{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 $\mathbf{{L}}$ in systems without rotational symmetry, such as the iron clusters investigated in Ref. Niemeyer et al. [2012]? 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.

Acknowledgements.

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 grant 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]. WMCF’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 NA 0002911. We acknowledge support from the Thomas Young Centre under grant TYC-101.

Appendix A 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, $(\varphi_{a}/c,\mathbf{{A}}_{a})=(\varphi(\mathbf{{r}}_{a},t)/c,\mathbf{{A}}(\mathbf{{r}}_{a},t))$ where $c$ is the speed of light. The single-particle Hamiltonian of a particle in an EM field White [2013], generalised to multiple interacting particles is,

[TABLE]

where particle $a$ has mass $m_{a}$, charge $q_{a}$, SOC parameter $\xi_{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 $\bm{p}_{a}-q_{a}\bm{A}_{a}$, we obtain the Lorentz force operator for particle $a$:

[TABLE]

where $\bm{E}_{a}=-\bm{\nabla}_{a}\varphi_{a}-\partial_{t}\mathbf{{A}}_{a}$ is the electric field operator on particle $a$ and $\mathbf{{v}}_{a}=\frac{1}{m_{a}}(\mathbf{{p}}_{a}-q_{a}\mathbf{{A}}_{a})$ 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 (35) 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 $\mathbf{{\Gamma}}^{N,EM}=\sum_{a\in N}\mathbf{{r}}_{a}\crossproduct\mathbf{{F}}_{a}^{EM}$, where $\mathbf{{r}}_{a}$ is the classical position of nucleus $a$ and $\mathbf{{F}}_{a}^{EM}$ is the direct electromagnetic force on nucleus $a$. The total torque exerted on the set $N$ of nuclei in the molecule is given by

[TABLE]

where $\mathbf{{\Gamma}}_{\textrm{int}}=\expectationvalue{\sum_{a\in N,b\in E}\mathbf{{r}}_{a}\crossproduct\mathbf{{F}}_{ab}}$ is the Hellmann-Feynman torque exerted on the nuclei by the set $E$ of electrons, and $\mathbf{{F}}_{ab}=-\bm{\nabla}_{a}V(\lvert\mathbf{{r}}_{a}-\mathbf{{r}}_{b}\rvert)$ is the operator for the Hellman-Feynman force on nucleus $a$ due to electron $b$.

Let $\mathbf{{J}}=\sum_{a\in E}\mathbf{{J}}_{a}$ be the total electronic canonical angular momentum, where $\mathbf{{J}}_{a}=\mathbf{{L}}_{a}+\mathbf{{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 $\mathbf{{J}}$ can be calculated using the Ehrenfest theorem:

[TABLE]

Simplifying the result by assuming that the applied $B$ field is uniform, working in the Coulomb gauge, and neglecting terms of second order in $\mathbf{{A}}$ yields

[TABLE]

and thus,

[TABLE]

where $\hat{\mathbf{B}}$ denotes the unit vector pointing in the direction of the $B$ field.