Current Induced Damping of Nanosized Quantum Moments in the Presence of Spin-Orbit Interaction

TL;DR

This paper develops a quantum mechanical formalism using Keldysh Green functions to study current-induced damping and magnetization dynamics in nanoscale ferromagnetic systems with spin-orbit interaction, revealing regimes of damping behavior and spin switching.

Contribution

It introduces a quantum approach to analyze current-induced magnetization dynamics, contrasting with classical models, and explores damping and switching phenomena in nanoscale ferromagnets with spin-orbit coupling.

Findings

Finite size effects prevent divergence of Gilbert damping in ballistic regime.

Bias can induce spin switching and demagnetization above a critical threshold.

Bias-induced antidamping efficiency varies significantly between regimes.

Abstract

Motivated by the need to understand current-induced magnetization dynamics at the nanoscale, we have developed a formalism, within the framework of Keldysh Green function approach, to study the current-induced dynamics of a ferromagnetic (FM) nanoisland overlayer on a spin-orbit-coupling (SOC) Rashba plane. In contrast to the commonly employed classical micromagnetic LLG simulations the magnetic moments of the FM are treated {\it quantum mechanically}. We obtain the density matrix of the whole system consisting of conduction electrons entangled with the local magnetic moments and calculate the effective damping rate of the FM. We investigate two opposite limiting regimes of FM dynamics: (1) The precessional regime where the magnetic anisotropy energy (MAE) and precessional frequency are smaller than the exchange interactions, and (2) The local spin-flip regime where the MAE and precessional frequency are comparable to the exchange interactions. In the former case, we show that due to the finite size of the FM domain, the \textquotedblleft Gilbert damping\textquotedblright does not diverge in the ballistic electron transport regime, in sharp contrast to Kambersky's breathing Fermi surface theory for damping in metallic FMs. In the latter case, we show that above a critical bias the excited conduction electrons can switch the local spin moments resulting in demagnetization and reversal of the magnetization. Furthermore, our calculations show that the bias-induced antidamping efficiency in the local spin-flip regime is much higher than that in the rotational excitation regime.

Equations

$$H_{tot}=\sum _{r{r}^{\prime },\sigma {\sigma }^{\prime }}T{r}_{\{s_d\}}\left[\left(1_s{\hat{H}}_{r{r}^{\prime }}^{\sigma {\sigma }^{\prime }}+{\delta }_{r{r}^{\prime }}{\delta }_{\sigma {\sigma }^{\prime }}{1}_s{\mu }_r+{\delta }_{r{r}^{\prime }}{J}_{sd}{\vec{\sigma }}_{\sigma {\sigma }^{\prime }}\cdot {\vec{\mathbf{s}}}_d(r)+{\delta }_{\sigma {\sigma }^{\prime }}{\delta }_{r{r}^{\prime }}{\mathbf{H}}_M\right){\psi }_{\{s_d^{\prime }\}{r}^{\prime }{\sigma }^{\prime }}^{\ast }{\psi }_{\{s_d\}r\sigma }\right].$$

$${\mathbf{H}}_M=$$

$$-\sum _{⟨r,r^{\prime }⟩}\frac{J_{r{r}^{\prime }}^{dd}}{s_d^2}{\vec{\mathbf{s}}}_d(r^{\prime })\cdot {\vec{\mathbf{s}}}_d(r)-\sum _r\frac{J_{sd}}{s_d}{\vec{\mathbf{s}}}_c(r)\cdot {\vec{\mathbf{s}}}_d(r),$$

$${\hat{H}}_{x{x}^{\prime }}^{\sigma {\sigma }^{\prime }}(k_{y}a)=[t\cos (k_{y}a){\delta }_{\sigma {\sigma }^{\prime }}-t_{so}\sin (k_{y}a){\sigma }_{\sigma {\sigma }^{\prime }}^x]{\delta }_{x{x}^{\prime }}$$

$$+t({\delta }_{x,x^{\prime }+1}+{\delta }_{x+1,x^{\prime }}){\delta }_{\sigma {\sigma }^{\prime }}+i{t}_{so}({\delta }_{x,x^{\prime }+1}-{\delta }_{x+1,x^{\prime }}){\sigma }_{\sigma {\sigma }^{\prime }}^y.$$

$$\left(E-i\eta -\hat{\mu }-\hat{H}-{\mathbf{H}}_M-\frac{J_{sd}}{2}\hat{\vec{\sigma }}\cdot {\widehat{\vec{\mathbf{s}}}}_d\right){\hat{\mathbf{G}}}^r(E)=\hat{1},$$

$$\hat{\mathbf{\rho }}=\int \frac{dE}{\pi }{\hat{\mathbf{G}}}^r(E)\eta f(E-\hat{\mu }-{\mathbf{H}}_M){\hat{\mathbf{G}}}^a(E).$$

$${\mathcal{T}}_m$$

$${\mathcal{T}}_m^{\pm }$$

$$\mathcal{T}={\sin }^2(\theta )t_{so}(c_1{t}_{so}\hslash \omega +c_{2}e{V}_{bias}),$$

$$\mathcal{T}=-\frac{\omega }{\eta }\sum _n\int \frac{d{k}_{x}d{k}_{y}d\phi }{(2\pi {)}^3}(v_{n\vec{k}}^{\phi }{)}^2{f}^{\prime }({\epsilon }_{n\vec{k}}(\phi )).$$

$$\alpha \approx \frac{1}{\eta }{t}_{so}^2\left[(k_F^{+}a{)}^2{D}^{+}(E_F)+(k_F^{-}a{)}^2{D}^{-}(E_F)\right](1+{\cos }^2(\gamma )),$$

$$[{\widehat{\vec{\mathbf{s}}}}_{d,r}{]}_{r_1{r}_2}={\vec{s}}_d(r_1){\delta }_{r_1{r}_2}(1-{\delta }_{r_{1}r})1_s+\frac{1}{2}{\delta }_{r_1{r}_2}{\delta }_{r_{1}r}\vec{\mathbf{\tau }},$$

$$\vec{B}(r)=g{\mu }_B{\vec{B}}^{ext}+4\sum _{r^{\prime }}{J}_{r{r}^{\prime }}^{dd}{\vec{s}}_d(r^{\prime })+2J_{sd}{\vec{s}}_c(r).$$

$$\left(E-\hat{\mu }-{\mathbf{H}}_M(r)-\hat{H}-\frac{J_{sd}}{2}\hat{\vec{\sigma }}\cdot {\widehat{\vec{\mathbf{s}}}}_{d,r}\right){\hat{\mathbf{G}}}_r^r(E)=\hat{1}.$$

$$\mp i\frac{\partial }{\partial t}{\mathbf{S}}^{\pm }(t)={\mathbf{h}}^z{\mathbf{S}}^{\pm }(t)-{\mathbf{h}}^{\pm }(t){\mathbf{S}}^z(t)$$

$$-i\frac{\partial }{\partial t}{\mathbf{S}}^z(t)=\frac{1}{2}\left({\mathbf{h}}^{+}(t){\mathbf{S}}^{-}(t)-{\mathbf{h}}^{-}(t){\mathbf{S}}^{+}(t)\right)$$

$${\vec{\mathbf{h}}}_{m{m}^{\prime }}(t)=\frac{1}{\hslash }\sum _r{J}_{sd}{\vec{s}}_c^{m{m}^{\prime }}(r)+g{\mu }_B{\delta }_{m{m}^{\prime }}\vec{B}(t),$$

$$\left(E-i\eta -\hat{H}(k)+n\omega -\frac{n}{2S}{J}_{sd}(k){\sigma }^z\right){\hat{\mathbf{G}}}_{nm}^r(E,k)$$

$$-\frac{\sqrt{S(S+1)-n(n+1)}}{2S}{J}_{sd}(k){\sigma }^{-}{\hat{\mathbf{G}}}_{n+1m}^r(E,k)$$

$$-\frac{\sqrt{S(S+1)-n(n-1)}}{2S}{J}_{sd}(k){\sigma }^{+}{\hat{\mathbf{G}}}_{n-1m}^r(E,k)=\hat{1}{\delta }_{nm}$$

$${\hat{\mathbf{\rho }}}_{nm}=e^{-i(n-m)\omega t}\sum _{p=-S}^S\int \frac{dE}{2\pi }{\hat{\mathbf{G}}}_{np}^{r}2\eta f_{p\hat{\mu }}{\hat{\mathbf{G}}}_{pm}^a$$

$${\hat{g}}_n^{u,r}(E,k)=\frac{1}{E-{\omega }_n-i{\eta }_n-\hat{H}(k)-{\hat{\Sigma }}_n^r(E,k)-\frac{n}{2S}{J}_{sd}(k){\sigma }^z-\frac{(S_n^{-}{)}^2}{4S^2}{J}_{sd}(k){\sigma }^{+}{\hat{g}}_{n-1}^{u,r}(E,k){\sigma }^{-}{J}_{sd}(k)}$$

$${\hat{\Sigma }}_n^{u,<}(E,k)=-\sum _{\alpha }\left(2i{\eta }_n+{\hat{\Sigma }}_{n,\alpha }^r(E,k)-{\hat{\Sigma }}_{n,\alpha }^a(E,k)\right)f_{n\alpha }+\frac{(S_n^{-}{)}^2}{4S^2}{J}_{sd}{\sigma }^{+}{\hat{g}}_{n-1}^{u,r}{\hat{\Sigma }}_{n-1}^{u,<}{\hat{g}}_{n-1}^{u,a}{\sigma }^{-}{J}_{sd}$$

$${\hat{g}}_n^{d,r}(E,k)=\frac{1}{E-{\omega }_n-i{\eta }_n-{\hat{\Sigma }}_n^r(E,k)-\hat{H}(k)-\frac{n}{2S}{J}_{sd}(k){\sigma }^z-\frac{(S_n^{+}{)}^2}{4S^2}{J}_{sd}(k){\sigma }^{-}{\hat{g}}_{n+1}^{u,r}(E,k){\sigma }^{+}{J}_{sd}(k)}$$

$${\hat{\Sigma }}_n^{d,<}(E,k)=-\sum _{\alpha }\left(2i{\eta }_n+{\hat{\Sigma }}_{n,\alpha }^r(E,k)-{\hat{\Sigma }}_{n,\alpha }^a(E,k)\right)f_{n\alpha }+\frac{(S_n^{+}{)}^2}{4S^2}{J}_{sd}{\sigma }^{-}{\hat{g}}_{n+1}^{d,r}{\hat{\Sigma }}_{n+1}^{d,<}{\hat{g}}_{n+1}^{d,a}{\sigma }^{+}{J}_{sd}$$

$${\hat{\mathbf{G}}}_{n,m}^r(E,k)$$

$$=\frac{S_n^{+}}{2S}{\hat{g}}_n^{u,r}(E,k)J_{sd}(k){\sigma }^{-}{\hat{\mathbf{G}}}_{n+1,m}^r(E,k),n\ne m$$

$$=\frac{S_n^{-}}{2S}{\hat{g}}_n^{d,r}(E,k)J_{sd}(k){\sigma }^{+}{\hat{\mathbf{G}}}_{n-1,m}^r(E,k),n\ne m$$

$${\hat{\Sigma }}_n^{r,u}=\frac{(S_n^{-}{)}^2}{4S^2}{J}_{sd}(k){\sigma }^{+}{\hat{g}}_{n-1}^{u,r}(E,k){\sigma }^{-}{J}_{sd}(k)$$

$${\hat{\Sigma }}_n^{r,d}=\frac{(S_n^{+}{)}^2}{4S^2}{J}_{sd}(k){\sigma }^{-}{\hat{g}}_{n+1}^{d,r}(E,k){\sigma }^{+}{J}_{sd}(k)$$

$${\mathcal{T}}_n^{\prime }=\frac{1}{2}({\mathcal{T}}_n^{\prime -}-{\mathcal{T}}_n^{\prime +})=\frac{1}{2}\Im \left(\sum _{k}Tr[\frac{S_n^{-}}{2S}{\sigma }^{+}{J}_{sd}(k){\hat{\mathbf{\rho }}}_{nn+1}(k)-\frac{S_n^{+}}{2S}{\sigma }^{-}{J}_{sd}(k){\hat{\mathbf{\rho }}}_{nn-1}(k)]\right)$$

$${\mathcal{T}}_n^{\prime -/+}=\sum _k\int \frac{dE}{2\pi i}\mathrm{Tr}\left\{\left[{\hat{\Sigma }}_n^{d/u,r}(E,k)-{\hat{\Sigma }}_n^{d/u,a}(E,k)\right]{\hat{\mathbf{G}}}_{nn}^{<}(E,k)+{\hat{\Sigma }}_n^{d/u,<}(E)\left[{\hat{\mathbb{G}}}_{nn}^r(E,k)-{\hat{\mathbb{G}}}_{nn}^a(E,k)\right]\right\},$$

Full text

I Introduction

Understanding the current-induced magnetization switching (CIMS) at the nanoscale is mandatory for the scalability of non-volatile magnetic random access memory (MRAM) of the next-generation miniaturized spintronic devices. However, the local magnetic moments of a nanoisland require quantum mechanical treatment rather than the classical treatment of magnetization commonly employed in micromagnetic simulations, which is the central theme of this work.

The first approach of CIMS employs the spin transfer torque (STT) Slonczewski1996 ; Berger1996 in magnetic tunnel junctions (MTJ) consisting of two ferromagnetic (FM) layers (i.e., a switchable free layer and a fixed layer) separated by an insulating layer, which involves spin-angular-momentum transfer from conduction electrons to local magnetization Ralph2008 ; Theodonis2006 . Although STT has proven very successful and brings the precious benefit of improved scalability, it requires high current densities ($\geq$ 1010 A/cm2) that are uncomfortably high for the MTJ’s involved and hence high power consumption. The second approach involves an in-plane current in a ferromagnet-heavy-metal bilayer where the magnetization switching is through the so-called spin-orbit torque (SOT) for both out-of-plane and in-plane magnetized layers. Gambardella2011 ; Miron2010 ; Liu2012 ; Liu2012b The most attractive feature of the SO-STT method is that the current does not flow through the tunnel barrier, thus offering potentially faster and more efficient magnetization switching compared to the MTJs counterparts.

As in the case of STT, the SO-STT has two components: a field-like and an antidamping component. While the field-like component reorients the equilibrium direction of the FM, the antidamping component provides the energy necessary for the FM dynamics by either enhancing or decreasing the damping rate of the FM depending on the direction of the current relative to the magnetization orientation as well as the structural asymmetry of the material. For sufficiently large bias the SOT can overcome the intrinsic damping of the FM leading to excitation of the magnetization precession. Liu2012b The underlying mechanism of the SOT for both out-of-plane and in-plane magnetized layers remains elusive and is still under debate. It results from either the bulk Spin Hall Effect (SHE) Dyakonov1971 ; Hirch1999 ; Jungwirth2012 ; Sinova2014 , or the interfacial Rashba-type spin-orbit coupling, Bijl2012 ; Kim2012 ; Kurebayashi2014 ; MahfouziPRB2016 or both Freimuth2014 ; Kim2013 ; Xiao2014 .

Motivated by the necessity of scaling down the size of magnetic bits and increasing the switching speed, the objective of this work is to develop a fully quantum mechanical formalism, based on the Keldysh Green function (GF) approach, to study the current-induced local moment dynamics of a bilayer consisting of a FM overlayer on a SOC Rashba plane, shown in Fig. 1.

Unlike the commonly used approaches to investigate the magnetization dynamics of quantum FMs, such as the master equation Chudnovskiy2014 , the scattering Fahnle2011 or quasi-classical Swiebodzinski2010 methods, our formalism allows the study of magnetization dynamics in the presence of nonequilibrium flow of electrons.

We consider two different regimes of FM dynamics: In the first case, which we refer to as the single domain dynamics, the MAE and the precession frequency are smaller than the exchange interactions, and the FM can be described by a single quantum magnetic moment, of a typically large spin, $S$, whose dynamics are governed mainly by the quantized rotational modes of the magnetization. We show that the magnetic degrees of freedom entering the density matrix of the conduction electron-local moment entagled system simply shift the chemical potential of the Fermi-Dirac distribution function by the rotational excitations energies of the FM from its ground state. We also demonstrate that the effective damping rate is simply the net current along the the auxiliary $m$-direction, where $m$= -S, -S+1, $\ldots$, +S, are the eigenvalues of the total $S_{z}$ of the FM. Our results for the change of the damping rate due to the presence of a bias voltage are consistent with the anti-damping SOT of classical magnetic moments, Li2015 ; MahfouziPRB2016 , where due to the Rashba spin momentum locking, the anti-damping SOT, to lowest order in magnetic exchange coupling, is of the form, $\vec{m}\times(\vec{m}\times\hat{y})$, where $\hat{y}$ is an in-plane unit vector normal to the transport direction.

In the adiabatic and ballistic transport regimes due to the finite S value of the nanosize ferromagnet our formalism yields a finite “Gilbert damping”, in sharp contrast to Kambersky’s breathing Fermi surface theory for damping in metallic FMs. Kambersky2007 On the other hand, Costa and Muniz  Costa2015 and Edwards Edwards2016 demonstrated that the problem of divergent Gilbert damping is removed by taking into account the collective excitations. Furthermore, Edwards points out Edwards2016 the necessity of including the effect of long-range Coulomb interaction in calculating damping for large SOC.

In the second case, which corresponds to an independent local moment dynamics, the FM has a large MAE and hence the rotational excitation energy is comparable to the local spin-flip excitation (exchange energy). We investigate the effect of bias on the damping rate of the local spin moments. We show that above a critical bias voltage the flowing conduction electrons can excite (switch) the local spin moments resulting in demagnetization and reversal of the magnetization. Furthermore, we find that, in sharp contrast to the single domain precessional dynamic, the current-induced damping is nonzero for in-plane and out-of-plane directions of the equilibrium magnetization. The bias-induced antidamping efficiency in the local moment switching regime is much higher than that in the single domain precessional dynamics.

The paper is organized as follows. In Sec. II we present the Keldysh formalism for the density matrix of the entagled quantum moment-conduction electron system and the effective dampin/antdamping torque. In Sec. III we present results for the current-induced damping rate in the single domain regime. In Sec. IV we present results for the current-induced damping rate in the independent local regime. We conclude in Sec. V.

II Theoretical Formalism

Fig. 1 shows a schematic view of the ferromagnetic heterostructure under investigation consisting of a 2D ferromagnet-Rashba plane bilayer attached to two semi-infinite normal (N) leads whose chemical potentials are shifted by the external bias, $V_{bias}$. The magnetization of the FM precesses around the axis specified by the unit vector, $\vec{n}_{M}$, with frequency $\omega$ and cone angle $\theta$. The FM has length, $L_{x}^{FM}$, along the transport direction. The total Hamiltonian describing the coupled conduction electron-localized spin moment system in the heterostructure in Fig. 1 can be written as,

[TABLE]

Here, $\vec{\bm{s}}_{d}(\mathbb{r})$ is the local spin moment at atomic position $\mathbb{r}$, the trace is over the different configurations of the local spin moments, $\{s_{d}\}$, ${\psi}_{\{s_{d}\}\mathbb{r}\sigma}=|\{s_{d}\}\rangle\otimes\psi^{e}_{\mathbb{r}\sigma}$ is the quasi-particle wave-function associated with the conduction electron ($\psi^{e}$) entangled to the FM states ($|\{s_{d}\}\rangle$), $J_{sd}$ is the $s-d$ exchange interaction, $\bm{1}_{s}$ is the identity matrix in spin configuration space, and $\hat{\sigma}_{x,y,z}$ are the Pauli matrices. We use the convention that, except for $\mathbb{r}$, bold symbols represent operators in the magnetic configuration space and symbols with hat represent operators in the single particle Hilbert space of the conduction electrons. The magnetic Hamiltonian $\bm{H}_{M}$ is given by

[TABLE]

where, the first term is the Zeeman energy due to the external magnetic field, the second term is the magnetic coupling between the local moments and the third term is the energy associated with the intrinsic magnetic field acting on the local moment, $\vec{\bm{s}}_{d}(\mathbb{r})$, induced by the local spin of the conduction electrons, $\vec{\bm{s}}_{c}(\mathbb{r})$.

The Rashba model of a two-dimensional electron gas with spin orbit coupling interacting with a system of localized magnetic moments has been extensively employed Haney2013 ; Park2013 ; Kim2012 to describe the effect of enhanced spin-orbit coupling solely at the interface on the current-induced torques in ultrathin ferromagnetic (FM)/heavy metal (HM) bilayers. The effects of (i) the ferromagnet inducing a moment in the HM and (ii) the HM with strong spin-orbit coupling inducing a large spin-orbit effect in the ferromagnet (Rashba spin-orbit coupling) lead to a thin layer where the magnetism and the spin-orbit coupling coexist. Haney2013

The single-electron tight-binding Hamiltonian Data_book for the conduction electrons of the 2D Rashba plane, ${H}^{\sigma\sigma^{\prime}}_{\mathbb{r}\mathbb{r}^{\prime}}$ which is finite along the transport direction $x$ and infinite along the $y$ direction is of the form,

[TABLE]

Here, $x,x^{\prime}$ denote atomic coordinates along the transport direction, $a$ is the in-plane lattice constant, and $t_{so}$ is the Rashba SOI strength. The values of the local effective exchange interaction, $J_{sd}$= 1eV, and of the nearest-neighbor hopping matrix element, $t$=1 eV, represent a realistic choice for simulating the exchange interaction of 3$d$ ferromagnetic transition metals and their alloys (Fe, Co). Zhang2004 ; Sanvito2005 ; Eastman1980 The Fermi energy, $E_{F}$=3.1 eV, is about 1 eV below the upper band edge at 4 eV consistent with the ab initio calculations of the (111) Pt surfaceKokalj1999 . Furthermore, we have used $t_{so}$=0.5 eV which yields a Rashba parameter, $\alpha_{R}$ = $t_{so}a\approx$1.4 eVÅ  ($a$=2.77 Å is the in-plane lattice constant of the (111) Pt surface) consistent with the experimental value of about 1-1.5 eVÅ Miron2011 and the ab initio value of 1 eVÅ Park2013 . However, because other experimental measurements for Pt/Co/Pt stacks reportHaazen2013 a Rashba parameter which is an order of magnitude smaller, in Fig.3 we show the damping rate for different values of the Rashba SOI.. For the results in Sec. IV, we assume a real space tight binding for propagation along $y$-axis.

The single particle propagator of the coupled electron-spin system is determined from the equation of motion of the retarded Green function,

[TABLE]

where, $\eta$ is the broadening of the conduction electron states due to inelastic scattering from defects and/or phonons, and for simplicity we ignore writing the identity matrices $\hat{1}$ and $\mathbb{1}$ in the expression. The density matrix of the entire system consisting of the noninteracting electrons (fermionic quasi-particles) and the local magnetic spins is determined (see Appendix A for details of the derivation for a single FM domain) from the expression,

[TABLE]

It is important to emphasize that Eq. (5) is the central result of this formalism which demonstrates that the effect of the local magnetic degrees of freedom is to shift the chemical potential of the Fermi-Dirac distribution function by the eigenvalues, $\varepsilon_{m}$, of $\bm{H}_{M}|m\rangle=\varepsilon_{m}|m\rangle$, i.e., the excitation energies of the FM from its ground state. Here, $|m\rangle$ are the eigenstates of the Heisenberg model describing the FM. The density matrix can then be used to calculate the local spin density operator of the conduction electrons, $[\vec{\bm{s}}_{c}(\mathbb{r})]^{mm^{\prime}}=\sum_{ss^{\prime}}\rho_{ss^{\prime},\mathbb{r}\mathbb{r}}^{mm^{\prime}}{\vec{\sigma}}_{ss^{\prime}}/2$, which along with Eqs. (2),  (4), and  (5) form a closed set of equations that can be solved self consistently. Since, the objective of this work is the damping/anti-damping (transitional) behavior of the FM in the presence of bias voltage, we only present results for the first iteration.

Eq. (5) shows that the underlying mechanism of the damping phenomenon is the flow of conduction electrons from states of higher chemical potential to those of lower one where the FM state relaxes to its ground state by transferring energy to the conduction electrons. Therefore, the FM dynamical properties in this formalism is completely governed by its coupling to the conduction electrons, where conservation of energy and angular momentum dictates the excitations as well as the fluctuations of the FM sate through the Fermi distribution function of the electrons coupled to the reservoirs. This is different from the conventional Boltzmann distribution function which is commonly used to investigate the thermal and quantum fluctuations of the magnetization.

Due to the fact that the number of magnetic configurations (i.e. size of the FM Hilbert space) grows exponentially with the dimension of the system it becomes prohibitively expensive to consider all possible eigenstates of the $\bm{H}_{M}$ operator. Thus, in the following sections we consider two opposite limiting cases of magnetic configurations. In the first case we assume a single magnetic moment for the whole FM which is valid for small FMs with strong exchange coupling between local moments and small MAE. In this case the dynamics is mainly governed by the FM rotational modes and local spin flips can be ignored. In the second case we ignore the correlation between different local moments and employ a mean field approximation such that at each step we focus on an individual atom by considering the local moment under consideration as a quantum mechanical object while the rest of the moments are treated classically. We should mention that a more accurate modeling of the system should contain both single domain rotation of the FM as well as the local spin flipping but also the effect of nonlocal correlations between the local moments and conduction electrons, which are ignored in this work.

III Single Domain Rotational Switching

In the regime where the energy required for the excitation of a single local spin moment ($\approx meV$) is much larger than the MAE ($\approx\mu eV$) the low-energy excited states correspond to rotation of the total angular momentum of the FM acting as a single domain and the effects of local spin flips described as the second term in Eq 2, can be ignored. In this regime all of the local moments behave collectively and the local moment operators can be replaced by the average spin operator, $\vec{\bm{s}}_{d}(\mathbb{r})=\sum_{\mathbb{r}^{\prime}}\vec{\bm{s}}_{d}(\mathbb{r}^{\prime})/N_{d}=s_{d}\vec{\bm{S}}/S$, where $N_{d}$ is the number of local moments and $\vec{\bm{S}}$ is the total angular momentum with amplitude $S$. The magnetic energy operator is given by $\bm{H}_{M}=-\vec{B}\cdot\bm{S}$, where, $\vec{B}=g\mu_{B}\vec{B}^{ext}+J_{sd}\vec{s}_{c}$. Here, for simplicity we assume $\vec{s}_{c}$ to be scalar and independent of the FM state. The eigenstates of $\bm{H}_{M}$ operator are then simply the eigenstates, $|S,m\rangle$, of the total angular momentum $\bm{S}_{z}$, with eigenvalues $m\omega=-S\omega,\ldots,+S\omega$, where $\omega=B_{z}$ is the Larmor frequency. Thus, the wave function of the coupled electron-spin configuration system, shown schematically in Fig. 2 is of the form, ${\psi}_{ms^{\prime}\mathbb{r}}(t)=|S,m\rangle\otimes\psi_{s^{\prime}\mathbb{r}}(t)$. One can see that the magnetic degrees of freedom corresponding to the different eigenstates of the $\bm{S}_{z}$ operator, enters as an additional auxiliary dimension for the electronic system where the variation of the magnetic energy, $\langle S,m|\bm{H}_{M}|S,m\rangle=m\omega$, shifts the chemical potentials of the electrons along this dimension. The gradient of the chemical potential along the auxiliary direction, is the Larmor frequency ($\mu eV\approx GHz$) which appears as an effective “electric field”in that direction.

Substituting Eq (5) in Eq (14)(b) and averaging over one precession period we find that the average rate of angular momentum loss/gain, which we refer to as the effective “damping rate”per magnetic moment, can be written as

[TABLE]

where,

[TABLE]

is the current along the auxiliary $m$-direction in Fig. 2 from the $m\leftrightarrow m+1$ ($\pm$ sign) state of the total $\bm{S}_{z}$ of the FM. Here, $Tr_{el}$, is the trace over the conduction electron degrees of freedom, and ${\bm{S}}_{m}^{\pm}=\sqrt{S(S+1)-m(m\pm 1)}$ are the ladder operators. It is important to note that within this formalism the damping rate is simply the net current across the $m^{th}$-layer along the auxiliary direction associated with the transition rate of the FM from state $m$ to its nearest-neighbor states ($m\pm 1$).

Fig. 3 shows the damping rate as a function of the precession cone angle, $\theta$ = cos${}^{-1}(\frac{m}{S})$, for different values of bias and for an in-plane effective magnetic field (a) along and (b) normal to the transport direction, and (c) an out-of-plane magnetic field. For cases (a) and (c) the damping rate is negative and relatively independent of bias for low bias values. A negative damping rate implies that the FM relaxes towards the magnetic field by losing its angular momentum, similar to the Gilbert damping rate term in the classical LLG equation, where its average value over the azimuthal precession angle, $\varphi=\omega t$, is of the form, $\mathcal{T}=-\alpha s_{d}\int\frac{d\varphi}{2\pi}\vec{m}\times(\vec{m}\times\vec{B})\cdot\vec{n}_{M}$, which is nonzero (zero) when the unit vector $\vec{n}_{M}$ is along (perpendicular to) the effective magnetic field. The dependence of the damping rate on the bias voltage when the effective magnetic field ${\vec{B}}$ is inplane and normal to the transport direction can be understood by the spin-flip reflection mechanism accompanied by Rashba spin-momentum locking described in Ref. MahfouziPRB2016 . One can see that a large enough bias can result in a sign reversal of the damping rate and hence a magnetization reversal of the FM. It’s worth mentioning that due to the zero-point quantum fluctuations of the magnetization, at $\theta=0,\pi$ (i.e. $m=\pm S$) we have $\mathcal{T}\neq$0 which is inversely proportional to the size of the magnetic moment, $S$.

In Fig. 4(a) we present the effective damping rate versus bias for different values of the Rashba SOC. The results show a linear response regime with respect to the bias voltage where both the zero-bias damping rate and the slope, $d\mathcal{T}/dV$ increases with the Rashba SOC. This is consistent with Kambersky’s mechanism of Gilbert damping due to the SOC of itinerant electrons, Kambersky2007 and the SOT mechanism MahfouziPRB2016 . Fig. 4(b) shows that in the absence of bias voltage the damping rate is proportional to $t_{so}^{2}$ and the effect of the spin current pumped into the left and right reservoirs is negligible. This result of the $t_{so}^{2}$ dependence of the zero-bias damping rate is in agreement with recent calculations of Costa and MunizCosta2015 and EdwardsEdwards2016 which took into account the collective excitations. In the presence of an external bias, $\mathcal{T}$ varies linearly with the SOC, suggesting that to the lowest order it can be fitted to

[TABLE]

where $c_{1}$ and $c_{2}$ are fitting parameters.

The bias-induced efficiency of the anti-damping SOT, $\Theta\equiv\hbar\omega(\mathcal{T}(V_{bias})-\mathcal{T}(0))/eV_{bias}\mathcal{T}(0)$, describes how efficient is the energy conversion between the magnetization dynamics and the conduction electrons. Accordingly, for a given bias-induced efficiency, $\Theta$, one needs to apply an external bias equal to $\hbar\omega/e\Theta$ to overcome the zero-bias damping of the FM. Fig. 5 displays the anti-damping efficiency versus the position of the Fermi energy of the FM from the bottom (-4$t$=-4 eV) to the top (4$t$=4eV) of the conduction electron band for the two-dimensional square lattice. The result is independent of the bias voltage and the Larmor frequency in the linear response regime (i.e. $V_{bias},\omega\ll t$). We find that the efficiency peaks when the Fermi level is in the vicinity of the bottom or top of the energy band where the transport is driven by electron- or hole-like carriers and the Gilbert damping is minimum. The sign reversal of the antidamping SOT is due to the electron- or hole-like driven transport similar to the Hall effect. Kittel_book

Classical Regime of the Zero Bias Damping rate — In the following we show that in the case of classical magnetic moments ($S\rightarrow\infty$) and the adiabatic regime ($\omega\rightarrow 0$), the formalism developed in this paper leads to the conventional expressions for the damping rate. In this limit the system becomes locally periodic and one can carry out a Fourier transformation from $m\equiv S_{z}$ space to azimuthal angle of the magnetization orientation, $\varphi$, space. Conservation of the angular momentum suggests that the majority- (minority-) spin electrons can propagate only along the ascending (descending) $m$-direction, where the hopping between two nearest-neighbor $m$-layers is accompanied by a spin-flip. As shown in Fig.  2 the existence of spin-flip hopping requires the presence of intralayer SOC-induced noncollinear spin terms which rotate the spin direction of the conduction electrons as they propagate in each $m$-layer. This is necessary for the persistent flow of electrons along the $\varphi$ auxiliary direction and therefore damping of the magnetization dynamics. Using the Drude expression of the longitudinal conductivity along the $\varphi$-direction for the damping rate, we find that, within the relaxation time approximation, $\eta/\omega\rightarrow\infty$, where the relaxation time of the excited conduction electrons is much shorter than the time scale of the FM dynamics, $\mathcal{T}$ is given by

[TABLE]

Here, $v_{n\vec{k}}^{\varphi}=\partial\varepsilon_{n\vec{k}}(\varphi)/\partial\varphi$ is the group velocity along the $\varphi$-direction in Fig. 2, and $\varepsilon_{n,\vec{k}}=\varepsilon^{0}(|\vec{k}|)\pm|\vec{h}(\vec{k})|$ for the 2D-Rashba plane, where $\varepsilon^{0}(|\vec{k}|)$ is the spin independent dispersion of the conduction electrons and $\vec{h}=at_{so}\hat{e}_{z}\times\vec{k}+\frac{1}{2}J_{sd}\vec{m}$, is the spin texture of the electrons due to the SOC and the $s-d$ exchange interaction. For small precession cone angle, $\theta$, the Gilbert damping constant can be determined from $\alpha=-\mathcal{T}/s_{d}\omega\sin^{2}(\theta)$, where the zero-temperature $\mathcal{T}$ is evaluated by Eq. (9). We find that

[TABLE]

where $D^{+(-)}(E)$ is the density of states of the majority (minority) band, $\gamma$ is the angle between the precession axis and the normal to the Rashba plane, and the Fermi wave-vectors ($k^{\pm}_{F}$) are obtained from, $\varepsilon_{0}(k^{\pm}_{F})=E_{F}\mp J_{sd}/2$. Eq. (10) shows that the Gilbert damping increases as the precession axis changes from in-plane ($\gamma=\pi/2$) to out of plane ($\gamma=0$), mahfouzi_SPIN2016 which can also be seen in Fig. 3.

It is important to emphasize that in contrast to Eq. (9) the results shown in Fig. 4 correspond to the ballistic regime with $\eta=0$ in the central region and the relaxation of the excited electrons occurs solely inside the metallic reservoirs. To clarify how the damping rate changes from the ballistic to the diffusive regime we present in Fig. 6 the damping rate versus the broadening, $\eta$, of states in the presence (solid line) and absence (dashed line) of bias voltage. We find that in both ballistic ($\eta/\omega\approx 0$) and diffusive ($\eta/\omega\gg 1$) regimes the damping rate is independent of the size of the FM domain, $S$. On the other hand, in the intermediate regime the FM dynamics become strongly dependent on the effective domain size where the minimum of the damping rate varies linearly with $S$. This can be understood by the fact that the effective chemical potential difference between the first, $m=-S$ and last, $m=S$ layers in Fig.3 is proportional to $S$ and for a coherent electron transport the conductance is independent of the length of the system along the transport direction. Therefore, in this case the FM motion is driven by a coherent dynamics.

IV Demagnetization Mechanism of Switching

In Sec. III we considered the case of a single FM domain where its low-energy excitations, involving the precession of the total angular momentum, can be described by the eigenstates $|m\rangle$ of $\bm{S}_{z}$ and local spin flip processes were neglected. However, for ultrathin FM films or FM nanoclusters, where the MAE per atom ($\approx meV$) is comparable to the exchange energy between the local moments (Curie temperature), the low-energy excitations involve both magnetization rotation and local moments spin-flips due to conduction electron scattering which can in turn change also $S$. In this case the switching is accompanied by the excitation of local collective modes that effectively lowers the amplitude of the magnetic ordering parameter. For simplicity we employ the mean field approximation for the 2D FM nanocluster where the spin under consideration at position $\mathbb{r}$ is treated quantum mechanically interacting with all remaining spins through an effective magnetic field, $\vec{B}$. The spatial matrix elements of the local spin operator are

[TABLE]

where, $\vec{\bm{\tau}}$s are the Pauli matrices. The magnetic energy can be expressed as, $\bm{H}_{M}(\mathbb{r})=-\vec{B}(\mathbb{r})\cdot\vec{\bm{\tau}}/2$, where, the effective local magnetic field is given by,

[TABLE]

The equation of motion for the single particle propagator of the electronic wavefunction entangled with the local spin moment under consideration can then be obtained from,

[TABLE]

The density matrix is determined from Eq. (5) which can in turn be used to calculate the spin density of the conduction electrons, $\vec{s}_{c}(\mathbb{r})=Tr(\hat{\vec{\sigma}}\hat{\bm{\rho}}_{\mathbb{r}\mathbb{r}})/2$, and the direction and amplitude of the local magnetic moments, $\vec{s}_{d}(\mathbb{r})=Tr(\vec{\bm{\tau}}\hat{\bm{\rho}}_{\mathbb{r}\mathbb{r}})/2$.

Fig. 7 shows the spatial dependence of the spin-$\frac{1}{2}$ local moment switching rate for a FM/Rashba bilayer (Fig. 1) for two bias values ($V_{bias}=\pm 0.4V$) and for an in-plane effective magnetic field (a) along and (b) normal to the transport direction, and (c) an out-of-plane magnetic field. The size of the FM island is $25a\times 25a$, where $a$ is the lattice constant. Negative local moment switching rate (blue) denotes that, once excited, the local moment relaxes to its ground state pointing along the direction of the effective magnetic field; however positive local damping rate (red) denotes that the local moments remain in the excited state during the bias pulse duration. Therefore, the damping rate of the local moments under bias voltage can be either enhanced or reduced and even change sign depending on the sign of the bias voltage and the direction of the magnetization. We find that the bias-induced change of the damping rate is highest when the FM magnetization is in-plane and normal to the transport directions similar to the single domain case. Furthermore, the voltage-induced damping rate is peaked close to either the left or right edge of the FM (where the reservoirs are attached) depending on the sign of the bias. Note that there is also a finite voltage-induced damping rate when the magnetization is in-plane and and along the transport direction ($x$) or out-of-the-plane ($z$).

Fig. 8 shows the bias dependence of the average (over all sites) damping rate for in- (a and b) and out-of-plane (c) directions of the effective magnetic field (direction of the equilibrium magnetization) and for two values of $|B|$. This quantity describes the damping rate of the amplitude of the magnetic order parameter. For an in-plane magnetization and normal to the transport direction (Fig. 8) the bias behavior of the damping rate is linear and finite in contrast to the single domain [Fig. 3(a)] where the damping rate was found to have a negligible response under bias. On the other hand, the bias behavior of the current induced damping rate shows similar behavior to the single domain case when the equilibrium magnetization direction is in-plane and normal to the transport direction (Fig. 8(b)). For an out-of-plane effective magnetic field [Fig. 8(c)] the damping torque has an even dependence on the voltage bias.

In order to quantify the efficiency of the voltage induced excitations of the local moments, we calculate the relative change of the average of the damping rate in the presence of a bias voltage and present the result versus the Fermi energy for different orientations of the magnetization in Fig 9. We find that the efficiency is maximum for an in-plane equilibrium magnetization normal to the transport direction and it exhibits an electron-hole asymmetry. The bias-induced antidamping efficiency due to spin-flip can reach a peak around 20% which is much higher than the peak efficiency of about 2% in the single domain precession mechanism in Fig. 5 for the same system parameters.

Future work will be aimed in determining the switching phase diagram MahfouziPRB2016 by calculating the local antidamping and field-like torques self consistently for different FM configurations.

V Concluding remarks

In conclusion, we have developed a formalism to investigate the current-induced damping rate of nanoscale FM/SOC 2D Rashba plane bilayer in the quantum regime within the framework of the Kyldysh Green function method. We considered two different regimes of FM dynamics, namely, the single domain FM and independent local moments regimes. In the first regime we assume the rotation of the FM as the only degree of freedom, while the second regime takes into account only the local spin-flip mechanism and ignores the rotation of the FM. When the magnetization (precession axis) is in-plane and normal to the transport direction, similar to the conventional SOT for classical FMs, we show that the bias voltage can change the damping rate of the FM and for large enough voltage it can lead to a sign reversal. In the case of independent spin-$1/2$ local moments we show that the bias-induced damping rate of the local quantum moments can lead to demagnetization of the FM and has strong spatial dependence. Finally, in both regimes we have calculated the bias-induced damping efficiency as a function of the position of the Fermi energy of the 2D Rashba plane.

Appendix A Derivation of Electronic Density Matrix

Using the Heisenberg equation of motion for the angular momentum operator, $\vec{\bm{S}}(t)$, and the commutation relations for the angular momentum, we obtain the following Landau-Lifshitz equations of motion,

[TABLE]

where, ${\bm{S}}^{\pm}={\bm{S}}_{x}\pm{\bm{S}}_{y}$ (${\sigma}^{\pm}={\sigma}_{x}\pm{\sigma}_{y}$), is the angular momentum (spin) ladder operators, $\vec{\bm{s}}^{mm^{\prime}}_{c}(\mathbb{r})=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}\vec{\sigma}_{\sigma\sigma^{\prime}}\rho^{mm^{\prime}}_{\sigma\sigma^{\prime},\mathbb{r}\mathbb{r}}$ is the local spin density of the conduction electrons which is an operator in magnetic configuration space. Here, $\rho$ is the density matrix of the system, and the subscripts, $\mathbb{r},m,\sigma$ refer to the atomic cite index, magnetic state and spin of the conduction electrons, respectively. In the following we assume a precessing solution for Eq (14)(a) with a fixed cone angle and Larmor frequency $\omega=h^{z}$. Extending the Hilbert space of the electrons to include the angular momentum degree of freedom we define ${\psi}_{ms^{\prime}i}(t)=|S,m\rangle\otimes\psi_{s^{\prime}i}(t)$. The equation of motion for the Green function (GF) is then given

[TABLE]

where, $n=(-S,-S+1,...,S)$ and the gauge transformation $\psi_{n\sigma i}(t)\rightarrow e^{in\omega t}\psi_{n\sigma i}(t)$ has been employed to remove the time dependence. The density matrix of the system is of the form

[TABLE]

where, $f_{p\hat{\mu}}(E)=f(E-p\omega-\hat{\mu})$ is the equilibrium Fermi distribution function of the electrons. Due to the fact that $p\omega$ are the eigenvalues of $\bm{H}_{M}=-g\mu_{B}\vec{B}\cdot\bm{S}$, one can generalize this expression by transforming into a basis where the magnetic energy is not diagonal which in turn leads to Eq (5) for the density matrix of the conduction electron-local moment entagled system.

Appendix B Recursive Relation for GFs

Since in this work we are interested in diagonal blocks of the GFs and in general for FMs at low temperature we have $S\gg 1$, we need a recursive algorithm to be able to solve the system numerically. The surface Keldysh GFs corresponding to ascending $\hat{g}^{u,r/<}$, and descending $\hat{g}^{d,r/<}$, recursion scheme read,

[TABLE]

where, $\hat{\Sigma}^{r}_{n}(E,k)=\sum_{\alpha}\hat{\Sigma}^{r}_{\alpha}(E-\omega_{n},k)$ corresponds to the self energy of the leads, $\alpha=L,R$ refers to the left and right leads in the two terminal device in Fig. 3 and $S^{\pm}_{m}=\sqrt{S(S+1)-m(m\pm 1)}$. Using the surface GFs we can calculate the GFs as follows,

[TABLE]

where the ascending and descending self energies are given by,

[TABLE]

The average rate of angular momentum loss/gain can be obtained from the real part of the loss of angular momentum in one period of precession,

[TABLE]

which can be interpreted as the current flowing across the layer $n$.

[TABLE]

Acknowledgements.

The work at CSUN is supported by NSF-Partnership in Research and Education in Materials (PREM) Grant DMR-1205734, NSF Grant No. ERC-Translational Applications of Nanoscale Multiferroic Systems (TANMS)-1160504, and US Army of Defense Grant No. W911NF-16-1-0487.