Terahertz spin dynamics driven by an optical spin-orbit torque

Spin torques are at the heart of spin manipulations in spintronic devices. Here, we examine the existence of an optical spin-orbit torque, a relativistic spin torque originating from the spin-orbit coupling of an oscillating applied field with the spins. We compare the effect of the nonrelativistic Zeeman torque with the relativistic optical spin-orbit torque for ferromagnetic systems excited by a circularly polarised laser pulse. The latter torque depends on the helicity of the light and scales with the intensity, while being inversely proportional to the frequency. Our results show that the optical spin-orbit torque can provide a torque on the spins, which is quantitatively equivalent to the Zeeman torque. Moreover, temperature dependent calculations show that the effect of optical spin-orbit torque decreases with increasing temperature. However, the effect does not vanish in a ferromagnetic system, even above its Curie temperature.


I. INTRODUCTION
Interest in controlling spins by means of circularly polarised pulses has grown immensely due to its potential applications in spin-based memory technologies [1][2][3]. Apart from the heat-assisted spin manipulations, the spins can also be controlled using the inverse Faraday effect (IFE) [4][5][6], the magnetic field of the terahertz pulses [7][8][9][10] and an optical spin-orbit torque (OSOT) that does not impart angular momentum into the spin system [11,12]. To be able to explain such effects theoretically, one has to simulate spin dynamics including several nonrelativistic and relativistic effects that might appear at ultrashort timescales [8][9][10][11][12].
The Landau-Lifshitz-Gilbert (LLG) equation of motion, consisting of precession of a spin moment around an effective field and a transverse relaxation, has been used extensively in the past to simulate such spin dynamics [13][14][15]. However, for a spin system excited by ultrashort laser pulses, a stochastic LLG equation of motion with atomistic resolution is required due to the strong thermal fluctuations in order to study the dynamical processes [16][17][18][19][20]. In these equations, a stochastic field is added to the effective field, in order to quantify the thermal fluctuations in the spin system. Nonetheless, at ultrashort timescales, the exact form of the LLG equation of motion has to be questioned as several other relativistic phenomena can occur. Therefore, we seek for an equation of motion that can capture all the possible interactions during ultrashort laser excitations of a spin system.
In a previous work, starting from the relativistic Dirac Hamiltonian yet including the magnetic exchange interactions, a rigorous derivation of the LLG equation of motion has been provided [21]. To treat the action of a laser pulse and the corresponding interactions, the Dirac-Kohn-Sham equation with external magnetic vector potential was considered [22][23][24]. To this end, a semirelativistic expansion of the Dirac-Kohn-Sham Hamiltonian that includes several nonrelativistic and relativistic spin-laser coupling terms was derived [25]. Having these coupling terms, an extended equation of motion that includes not only the spin precession and damping, but also other relativistic torque terms * ritwik.mondal@uni-konstanz.de was obtained [26]. One of these torque terms is the fieldderivative torque which appears due to the time-dependent field excitation, e.g., in case of THz pulse excitation [27]. Another new torque term is the OSOT, which stems from the spin-orbit interaction of the applied field with the electron spins i.e., it imparts spin-angular momentum of the applied field to the spins. The new equation of motion for the reduced magnetization vector m i (t) including OSOT cast in the form of an LLG equation is We define the gyromagnetic ratio as γ and D, representing the damping parameter, is in general a tensor. For simplification, we consider only a scalar damping parameter α that can be expressed as α = 1 3 Tr(D) [21]. The effective field is represented as B eff i , which is the derivative of the total magnetic energy without the relativistic light-spin interaction with respect to m i as will be specified later on in detail. The additional field B OSOT describes the optomagnetic field induced by the laser pulse due to the OSOT phenomenon.
In this article, we investigate the effect of the OSOT term within atomistic spin dynamics simulations. First, we revisit the derivation of the OSOT from a general spin-orbit coupling Hamiltonian that has been derived within a relativistic formalism [21]. We find that the OSOT depends on the helicity, frequency and intensity of the light pulse. If the intensity is high and frequency is low, we expect the OSOT terms to show the most significant effects. We simulate the spin dynamics for a spin model, representative for bcc Fe, with the OSOT and find that the OSOT can provide significant contributions at the THz regime. Studying the temperature dependence, we find that the OSOT effects is robust against thermal fluctuations, i.e., we observe no significant reduction of the strength of the OSOT even up to the critical temperature.

II. OPTICAL SPIN-ORBIT TORQUE
The generalized spin-orbit coupling (SOC) Hamiltonian, as derived within the fully relativistic Dirac framework, can arXiv:2008.07308v2 [cond-mat.mes-hall] 15 Feb 2021 be written as [21,23,26] where the electric fields are represented as E tot = E int + E ext and E ext = − ∂A ∂t −∇Φ with (A, Φ) as magnetic vector and scalar potentials. The physical constants have their usual meanings and σ denotes the electron's spin through Pauli spin matrices and p is the electron momentum.
Note that the SOC can occur in several ways, as described in the following: (i) the angular momentum of an electron couples to the spin of the electron -this can be expressed as σ · (E int × p), (ii) the first-order magnetic vector potential of the electromagnetic (EM) field couples to the spins -this can be expressed as σ · (E int × A), (iii) the spin angular momentum of the EM field couples to the spin -this can be expressed as σ · (E ext × A).
In a spherically symmetric potential, the first type SOC can be written as traditional l · s coupling, where l and s represent the orbital and spin angular momentum respectively. It provides explanations to several relativistic effects in magnetism e.g., magnetic Gilbert damping [21,28,29] and many others [30,31]. Such a SOC exists even without an external field. In contrary, the latter two types of SOC depend explicitly on the external field and can be written as The total internal field depends on the intrinsic field E 0 int that even exist without the applied field and the applied field itself. Within the linear response theory we write E int = E 0 int + ζE ext , where ζ defines the coupling strength of the applied EM field relative to the intrinsic one. We thus consider the optical spin-orbit coupling as Using the definition of Zeeman coupling of −µ s σ · B OSOT , the optical spin-orbit coupling can be shown to induce a field (see Appendix A for details) where B 0 (= E 0 /c), η and ω are the envelope of the oscillating Zeeman field, helicity and angular frequency of the elliptically polarised light, respectively. µ s defines the magnetic moment. Note, that due to the OSOT, the induced field points along the direction of the energy flux of the propagating wave. Additionally, note that the field B OSOT is largest for circularly polarised light (η = ±π/2). However, if the coupling strength ζ is not the same for right and left circularly polarised light, their combined effect could lead to nonzero B OSOT for linearly polarised light. The parameter ζ depends on the electron density and absorption of the light that can be different for right and left circularly polarised light, leading to a magnetic circular dichroism (MCD). In the following, we simulate the Zeeman effect from the electromagnetic THz field and the optical spinorbit coupling effects simultaneously, and make a comparison between these two effects in a ferromagnetic system.

III. ATOMISTIC SPIN SIMULATIONS
In order to compute the spin dynamics, it is convinient to transform the implicit form of our equation of motion to the explicit Landau-Lifshitz (LL) form. For a scalar damping parameters, α, the LLG equation (1) can be recast as The effective field, B eff i , is the derivative of total energy with respect to the magnetic moment, B eff (6), hereby consists of the so called field-like (see also Fig. 1) and the weaker damping-like torque which is proportional to the Gilbert damping coefficient α = 0.01 and describes the coupling of the spins to a heat bath.
In the following we consider a spin model for bcc Fe. The total Hamiltonian of the system (without the relativistic spin-light coupling term) H can be expressed as The first term describes the traditional Heisenberg exchange energy, considering the exchange constants J ij up to the third nearest neighbour. The second term represents the uniaxial anisotropy with energy constant d z and the last term is the Zeeman coupling with a time- The ellipticity of the applied pulse is taken into account through B 0 = B0 √ 2 ŷ + e iηẑ , which considers the major and minor axes of the ellipse and τ defines the pulse duration. In the following we use the abbreviations Zeeman field (ZF), Zeeman torque (ZT), and optical spin-orbit field (OSOF).
A striking difference between ZT and OSOT is that the ZT does not depend on the frequency of the pulse and is proportional to the amplitude of the applied field pulse. However, the OSOT depends inversely on the frequency and is proportional to the square of the field amplitude, i.e. the intensity. To quantify the OSOT effects, one can estimate its amplitude by computing the characteristic field constant of the OSOF Considering a central frequency of f = 1 THz, we find that B(2πf, 0) = 157 T. Hence an electromagnetic field amplitude of about 10 T could induce a OSOF of around 1 T, even for the low limit of the coupling strength, ζ → 0. Needless to mention that the OSOT effects might exceed the ZT if ζ is higher. However, the results shown here were computed for the limit ζ → 0. Such ultra-intense THz fields are currently only achievable with linear polarisation, with recent experiments pointing towards achieving circular polarized laser pulses in the terahertz frequency regime [32,33].

IV. APPLICATION TO FERROMAGNETS
In our following simulations, we consider bcc Fe as a ferromagnetic sample. For the zero-temperature simulations it is sufficient to solve the LLG equation for a single spin, due to the homogeneity of the THz pulse excitation. The choice of exchange parameters is only relevant at elevated temperatures. We shine an electromagnetic pulse having 1 ps pulse duration along thex-direction with y and z field components as shown in the Fig. 1. We will mainly discuss the limit d z ≈ 0, since the uniaxial anisotropy does not have a significant impact on the ultrafast spin dynamics (see Appendix B for details).

A. Zero temperature simulations
In particular we compare the effect of ZT and OSOT at zero temperature. As the optical pulse is applied along thê x-direction (k = |k|x), we calculate the change in magnetization for y and z components in Fig. 2(a). For the case of only ZT the induced OSOT remains obviously zero. For a maximum of 10 T applied ZF, the change in magnetization is about 50 % for the y-component and 5 % for the z-component. The reason for the triggered magnetisation dynamics is the following: the ZF from the optical pulse has two components B y and B z . The B y component exerts a torque along −x on the equilibrium spin direction: with equilibrium magnetisation m 0 . On the other hand, the B z component of the EM field does not exert any torque on the Fe spins at first as the spins are initially aligned alongẑ. However, as soon as the above mentioned torque induced some magnetisation ∆m x , the B z component of the EM pulse, one quarter period later, exerts a torque alongŷ direction: Therefore, there exists a superposition between two torques for the ZF. Note that the change in ∆m y magnetization is antisymmetric in the helicity of light, which suggests that the effect is similar to IFE [34]. Further note that the ab initio calculations of the IFE were calculated using a nonrelativistic theory [35,36]. However, present theory is based on a relativistic interaction Hamiltonian. According to Eq. (5), the induced field diverges in the limit ω → 0, which can be explained by the previous theories of IFE [37][38][39], and is in accordance with the ab initio calculations of IFE [35,36,40]. We would also like to mention that these ab initio calculations showed the IFE being asymmetric in helicity for ferromagnets which includes the absorption of light [36].
The OSOF, on the other hand, has only one component B x alongx direction, see Fig. 2 and changes the magnetization by about 50 % as shown in Fig. 2 (b). According to our theory in Eq. (5), the right and left circular polarization will induce an OSOF along the +x and −x directions respectively. Therefore, the right circular pulse exerts a torque alongŷ (viz.ẑ ×x =ŷ) and similarly, the left circular pulse exerts a torque along −ŷ.
The change in magnetization is also antisymmetric in the helicity of the light pulse, similar to ZT effects in ∆m y . Note that we have assumed that the material dependent parameter ζ is zero, which dictates that the antisymmetric behavior is not universal. In fact, the parameter ζ could be different for right and left circular pulse, meaning that the antisymmetry of the plots above would be broken [36]. The effective contributions of the combined ZT and OSOT have been computed in Fig. 2 (c). We note that the magnetization change in m y is opposite in helicity for ZT and OSOT (e.g., see fifth row plots in Figs. 2 (a) and 2 (b)). Therefore, the final effective contribution becomes only about 5 % change in ∆m y . For the z-component of magnetization, the effective ∆m z is negligibly small, even though individual changes are about 5 % due to ZT and OSOT.
To quantitatively understand the ZT and OSOT effects, we compute the field dependent contributions in terms of the maximum change of each magnetization components ∆m i as a function of the applied field B 0 for left circular polarised light, shown in Fig. 3. From a simple scaling argument we would expect the spin excitation ∆m ⊥ = ∆m xx + ∆m yŷ scaling with the magnitude of B EM and B OSOF , respectively. For the EM field, this is simply the amplitude B EM ∝ B 0 , whereas for the OSOF, that is B OSOF ∝ B 2 0 according to Eq. (5). As already described above, Eq. (9), the spin excitation ∆m x for the EM field is induced by the field-like torque from its B y component, thus, Fig. 3 (a) shows a linear scaling with the EM field amplitude B 0 . The ∆m y component in Fig. 3 (b), however, is induced in a second order process,  described in Eq. (10) and thus scales quadratically with the EM field amplitude B 0 , i.e., ∆m y ∝ B 2 0 . A deviation of this scaling law is observed at lower field amplitudes of B 0 100 mT, where the damping-like torque of the EM field with ∆m y ∝ αB 0 is taking over. This effect can also be seen by comparing the torque amplitudes, which differ by a factor of ∆m y /∆m x ≈ α in this regime. The ∆m z component, i.e., the deviation of the equilibrium component then follows from conservation of magnetization amplitude: ∆m z = m 0 − m 0 − ∆m 2 ⊥ . At low excitation the dominant contribution here is the ∆m x component and one can simplify ∆m z ≈ −∆m 2 x /2m 0 ∝ −B 2 0 . In contrast to the EM field, there is only one field-like torque acting along the +ŷ direction due to the OSOF in x direction. This torque can be seen in Fig. 3 (b), and shows the quadratic scaling ∆m y ∝ −B x ∝ −B 2 0 implied by Eq. (11). The ∆m x excitation in Fig. 3 (a) on the other hand is due to the damping-like OSOT and thus following the same power law, but a factor of α smaller compared to the field-like excitation in Fig. 3 (b). Finally, the ∆m z excitation of the OSOT follows again from conservation of magnetization amplitude, leading to ∆m z ≈ −∆m 2 y /2m 0 and implying ∆m z ∝ −B 4 0 as shown in Fig. 3 (c). These scaling observations are once more summarized in Tab. I where we display the scaling exponents obtained by fitting our simulation data. An interesting observation here is that although the characteristic field of the OSOT is on the order of 100 T in this frequency regime, due to the different torque symmetry compared to the ZT, the OSOT is by no means negligible. For the ∆m y excitation we find that the strength of ZT and OSOT are comparable, though opposite in sign, at fields on the order of 10 2 mT-a value that is in much closer reach experimentally. Altogehter, the OSOT effects can thus provide an equivalent torque compared to the Zeeman effect. Therefore, the OSOT effects cannot be neglected when a circularly polarised light acts on a magnetic system even at the weak coupling limit ζ → 0.
Up to now we did not take into account the role of a finite anisotropy. In order to investigate this, we performed   the same calculation depicted in Fig. 2 with the anisotropy d z = 7.66 µeV for Fe included. We found that the main effect of the anisotropy field is the precession of the induced magnetization ∆m ⊥ around the z-axis over time (see Appendix B for details). On the other hand, no direct effect of the anisotropy can be observed on the ultrafast time scales, i.e., on the time scale of the pulse duration.

B. Finite temperature simulations
For the finite temperature simulations, the exchange parameters of Pajda et al. [41] are used, where the first two nearest neighbors are strongly ferromagnetically, however, the third nearest neighbor is weakly antiferromagnetically coupled. For these simulations, we also use the uniaxial anisotropy of d z = 7.66 µeV along the z-axis, in order to align the magnetization. Calculating the temperature dependence of the OSOT for ferromagnetic Fe, we use a simulation grid of 48 3 spins and add a stochastic field to the effective field in Eq. (6), in order to treat the thermal fluctuations [17]. The calculated Curie temperature of this system is T C = 1368 K and thus slightly higher than the true value of Fe. However, since the Curie temperature is the only temperature scale in our simulations, we can basically treat it as a free scaling parameter.
We compare in Fig. 4 the OSOT effects at T = 0 and T = 0.73 × T C (0 K and 1000 K). We find in Fig. 5 that though the OSOT effect in ∆m i appears to decrease with increasing temperature, this reduction is only related to the thermal reduction of magnetic order m eq (T ) at finite temperature. In other words the rotation angle of the normalized magnetization is not sensitive to the temperature. To illustrate this further, we have systematically calculated the temperature dependence of the OSOT by calculating the difference between maximum spin excitation for right and left circular pulses in Fig. 5. For each temperature, we performed ten simulations for each circular pulses, and took the average to determine max [∆m R,y ] − max [∆m L,y ] as a function of temperature, which ensures that the thermal fluctuations are minimised. Far away from T C the spin excitation amplitude is proportional to the equilibrium magnetization curve. Only in the close vicinity of the critical temperature, we find an increase of the net spin excitation, relative to the equilibrium magnetisation. This is simply due to the large amplitude of both ZF and OSOF which induce a transient magnetic order. Therefore, OSOT effects should be observed for bcc Fe even at elevated temperatures i.e., at realistic conditions for ultra-intense spin excitations.

V. CONCLUSIONS
To conclude, we incorporated a new torque into the LLG equation that should appear in ultrafast spin dynamics, namely OSOT, and investigated this effect via computer simulations of an atomistic spin model, representative for bcc Fe. The OSOT originates from the spin-orbit coupling of the electron spins to an external EM field. The strength of this OSOT, unlike the first-order ZT, depends on the intensity of the ultrafast light pulse, as well as on the frequency. In addition, the OSOT depends on the ellipticity of the pulse and it provides maximum torque for circularly polarised light pulses. Throughout the simulations presented here, we considered the weak coupling limit ζ → 0. However, the coupling depends on the electronic configurations of the system and could potentially further increase the strength of the OSOT. Although the OSOT is a higher order contribution in the external field, we found that the ZT and the OSOT provide quantitatively equivalent torques on the spins for circularly polarised laser pulses in the magnetization component perpendicular to the k-vector and equilibrium magnetisation m 0 . The effect of the OSOT resembles an already known effect, namely the IFE and it can be considered as a relativistic contributions to the IFE. The temperature dependence study of OSOT shows that the OSOT effect is present at elevated temperatures, even up to the Curie temperature.

ACKNOWLEDGMENTS
We thank László Szunyogh for valuable discussions and acknowledge financial support from the Alexander von Humboldt-Stiftung, Zukunftskolleg at Universität Konstanz via grant No.
Appendix A: Derivation of optical spin-orbit torque Following Eq. (4), the induced field can be written as We use the time dependent field as E ext = R E 0 e i(k·r−ωt) and the amplitude as the elliptically polarised light E 0 = E0 √ 2 ŷ + e iηẑ , with η as the ellipticity of the light. Therefore, the electric field can be written as (when only the timedependent part is taken) E ext = ∂A ∂t ⇒ A = − E ext dt that can be calculated as follows Therefore, the induced field can be taken from Eq. (A1) as In the last step of the calculation, we used the relation E 0 = cB 0 . In our simulations, we apply time-dependent Zeeman fields along y and z-directions and the corresponding induced optical spin-orbit field acts along x-direction.

Appendix B: Effect of anisotropy
Here, we compute the influence of uniaxial magnetic anisotropy on the spin dynamics induced by the ZT and OSOT. Fig. 6 shows the magnetization dynamics taking into account the uniaxial anisotropy for bcc Fe. The main effect of anisotropy can be noticed by comparing with the zero-anisotropy simulations in Fig. 2. A small increase in ∆m x is observed in the case of only ZF or OSOF after the pulse has passed. This is due to the slow precession of the induced ∆m y component which starts to precess around the anisotropy field alongẑ axis. In case of the superposition of ZF and OSOF the excitation of ∆m y is mostly compensated and therefore no ∆m x emerges either. Additionally, we mention that the anisotropy energies do not affect the ∆m y and ∆m z on the ultrafast time scale and the net excitation remains the same irrespective of the anisotropyat least for the typically weak magnetic anisotropies of the 3d ferromagnets.