Introduction

In spite of being a mature field of research, studying magnetism and spin-dependent phenomena in solids still remains one of the most exciting area in modern condensed matter physics. In fact, enormous progress in technological development over the last few decades is mainly held by the achievements in spintronics and related fields1,2,3,4,5,6,7,8,9,10,11. However the theoretical description of magnetization dynamics is at best accomplished on the level of Landau-Lifshitz-Gilbert (LLG) equation that characterizes torques on the magnetization. In essence, this equation describes the precession of the magnetization, m(r, t), about the effective magnetic field, Heff(r, t), created by the localized moments in magnetic materials, and its relaxation to equilibrium. The latter, known as the Gilbert damping torque12, was originally captured in the form αm × ∂tm, where the parameter α determines the relaxation strength, and it was recently shown to originate from a systematic non-relativistic expansion of the Dirac equation13. Thus, a proper microscopic determination of the damping parameter α (or, the damping tensor in a broad sense) is pivotal to correctly simulate dynamics of magnetic structures for the use in magnetic storage devices14.

From an experimental viewpoint, the Gilbert damping parameter can be extracted from ferromagnetic resonance linewidth measurements15,16,17 or established via time-resolved magneto-optical Kerr effect18,19. In addition, it was clearly demonstrated that in bilayer systems made of a nonmagnetic metal (NM) and a ferromagnet material (FM) the Gilbert damping is drastically enhanced as compared to bulk FMs20,21,22,23,24. A strong magnetocrystalline anisotropy, present in CoNi, CoPd, or CoPt, hints unambiguously for spin-orbit origin of the intrinsic damping. A first theoretical attempt to explain the Gilbert damping enhancement was made in terms of sd exchange model in ref.25. Within this simple model, magnetic moments associated with FM layer transfer angular momentum via interface and finally dissipate. Linear response theory has been further developed within free electrons model26,27, while the approach based on scattering matrix analysis has been presented in refs28,29. In the latter scenario spin pumping from FM to NM results in either backscattering of magnetic moments to the FM layer or their further relaxation in the NM. Furthermore, the alternative method to the evaluation of the damping torque, especially in regard of first-principles calculations, employs torque-correlation technique within the breathing Fermi surface model30. While a direct estimation of spin-relaxation torque from microscopic theory31, or from spin-wave spectrum, obtained on the basis of transverse magnetic field susceptibility32,33, are also possible. It is worth mentioning that the results of first-principles calculations within torque-correlation model34,35,36,37,38 and linear response formalism39,40 reveal good agreement with experimental data for itinerant FMs such as Fe, Co, or Ni and binary alloys.

Last but not least, an intensified interest towards microscopic foundations of the Gilbert parameter α is mainly attributed to the role the damping torque is known to play in magnetization reversal41. In particular, according to the breathing Fermi surface model the damping stems from variations of single-particle energies and consequently a change of the Fermi surface shape depending on spin orientation. Without granting any deep insight into the microscopic picture, this model suggests that the damping rate depends linearly on the electron-hole pairs lifetime which are created near the Fermi surface by magnetization precession. In this paper we propose an alternative derivation of the Gilbert damping tensor within a mean-field approach according to which we consider itinerant subsystem in the presence of nonequilibrium classical field m(r, t). Subject to the function m(r, t) is sufficiently smooth and slow on the scales determined by conduction electrons mean free path and scattering rate, the induced nonlocal spin polarization can be approached within a linear response, thus providing the damping parameter due to the itinerant subsystem. In the following, we provide the derivation of a Kubo-Středa formula for the components of the Gilbert damping tensor and illustrate our approach for a two-dimensional Rashba ferromagnet, that can be modeled by the interface between NM and FM layers. We argue that our theory can be further applied to identify properly the tensorial structure of the Gilbert damping for more complicated model systems and real materials.

Microscopic Framework

Consider a heterostructure made of NM with strong spin-orbit interaction covered by FM layer as shown in Fig. 1, e.g., CoPt. In general FMs belong to the class of strongly correlated systems with partially filled d or f orbitals which are responsible for the formation of localized magnetic moments. The latter can be described in terms of a vector field m(r, t) referred to as magnetization, that in comparison to electronic time and length scales slowly varies and interacts with an itinerant subsystem. At the interface (see Fig. 1) the conduction electrons of NM interact with the localized magnetic moments of FM via a certain type of exchange coupling, sd exchange interaction, so that the Hamiltonian can be written as

$$h=\frac{{p}^{2}}{2m}+\alpha \,{({\boldsymbol{\sigma }}\times {\boldsymbol{p}})}_{z}+{\boldsymbol{\sigma }}\cdot {\boldsymbol{M}}({\boldsymbol{r}},t)+U({\boldsymbol{r}}),$$
(1)

where first two terms correspond to the Hamiltonian of conduction electrons, on condition that the two-dimensional momentum p = (px, py) = p(cos φ, sin φ) specifies electronic states, m is the free electron mass, α stands for spin-orbit coupling strength, while σ = (σx, σy, σz) is the vector of Pauli matrices. The third term in (1) is responsible for sd exchange interaction with the exchange field M(r, t) = Δm(r, t) aligned in the direction of magnetization and Δ denoting sd exchange coupling strength. We have also included the Gaussian disorder, the last term in Eq. (1), which represents a series of point-like defects, or scatterers, 〈U(r)U(r′)〉 = (mτ)−1δ(r − r′) with the scattering rate τ (we set ħ = 1 throughout the calculations and recover it for the final results).

Figure 1
figure 1

Schematic representation of the model system: the electrons at the interface of a bilayer, composed of a ferromagnetic (FM) and a nonmagnetic metal (NM) material, are well described by the Hamiltonian (1). We assume the magnetization of FM layer depicted by the red arrow is aligned along the z axis.

Full size image

Subject to the norm of the vector |m(r, t)| = 1 remains fixed, the magnetization, in broad terms, evolves according to (see, e.g., ref.42),

$${\partial }_{t}{\boldsymbol{m}}={\boldsymbol{f}}\times {\boldsymbol{m}}=\gamma {{\boldsymbol{H}}}_{{\rm{eff}}}\times {\boldsymbol{m}}+\chi {\boldsymbol{s}}\times {\boldsymbol{m}},$$
(2)

where f corresponds to so-called spin torques. The first term in f describes precession around the effective magnetic field Heff created by the localized moments of FM, whereas the second term in (2) is determined by nonequilibrium spin density of conduction electrons of NM at the interface, s(r, t). It is worth mentioning that in Eq. (2) the parameter γ is the gyromagnetic ratio, while χ = (gμB/ħ)2μ0/d is related to the electron g–factor (g = 2), the thickness of a nonmagnetic layer d, with μB and μ0 standing for Bohr magneton and vacuum permeability respectively. Knowing the lesser Green’s function, G< (rt; rt), one can easily evaluate nonequilibrium spin density of conduction electrons induced by slow variation of magnetization orientation,

$${s}_{\mu }({\boldsymbol{r}},t)=-\,\frac{i}{2}\,{\rm{Tr}}\,[{\sigma }_{\mu }{G}^{ < }({\boldsymbol{r}}t;{\boldsymbol{r}}t)]={Q}_{\mu \nu }{\partial }_{t}{m}_{\nu }+\ldots ,$$
(3)

where summation over repeated indexes is assumed (μ, ν = x, y, z). The lesser Green’s function of conduction electrons can be represented as G< = (GK − GR + GA)/2, where GK, GR, GA are Keldysh, retarded, and advanced Green’s functions respectively.

Kubo-Středa Formula

We further proceed with evaluating Qμν in Eq. (3) that describes the contribution to the Gilbert damping due to conduction electrons. In the Hamiltonian (1) we assume slow dynamics of the magnetization, such that approximation M(r, t) ≈ M + (t − t0)∂tM with M = M(r, t0) is supposed to be hold with high accuracy,

$$H=\frac{{p}^{2}}{2m}+\alpha {({\boldsymbol{\sigma }}\times {\boldsymbol{p}})}_{z}+{\boldsymbol{\sigma }}\cdot {\boldsymbol{M}}+U({\boldsymbol{r}})+(t-{t}_{0})\,{\boldsymbol{\sigma }}\cdot {\partial }_{t}{\boldsymbol{M}},$$
(4)

where first four terms in the right hand side of Eq. (4) can be grouped into the Hamiltonian of a bare system, H0, which coincides with that of Eq. (1), provided by the static magnetization configuration M. In addition, the expression (4) includes the time-dependent term V(t) explicitly, as the last term. In the following analysis we deal with this in a perturbative manner. In particular, the first order correction to the Green’s function of a bare system induced by V(t) is,

$$\delta G({t}_{1},{t}_{2})={\int }_{{C}_{K}}\,dt\,\int \,\frac{{d}^{2}p}{{\mathrm{(2}\pi )}^{2}}{g}_{{\boldsymbol{p}}}({t}_{1},t)V(t){g}_{{\boldsymbol{p}}}(t,{t}_{2}),$$
(5)

where the integral in time domain is taken along a Keldysh contour, while gp(t1, t2) = gp(t1 − t2) [the latter accounts for the fact that in equilibrium correlation functions are determined by the relative time t1 − t2] stands for the Green’s function of the bare system with the Hamiltonian H0 in momentum representation. In particular, for the lesser Green’s function at coinciding time arguments \({t}_{1}={t}_{2}\equiv {t}_{0}\), which is needed to evaluate (3), one can write down,

$$\begin{array}{rcl}\delta {G}^{ < }({t}_{0},{t}_{0}) & = & \frac{i}{2}\,{\int }_{-\infty }^{\infty }\,\frac{d\varepsilon }{2\pi }\,\int \,\frac{{d}^{2}p}{{\mathrm{(2}\pi )}^{2}}\{{g}_{{\boldsymbol{p}}}^{R}{\sigma }_{\mu }\frac{\partial {g}_{{\boldsymbol{p}}}^{ < }}{\partial \varepsilon }-\frac{\partial {g}_{{\boldsymbol{p}}}^{R}}{\partial \varepsilon }{\sigma }_{\mu }{g}_{{\boldsymbol{p}}}^{ < }\\ & & +\,{g}_{{\boldsymbol{p}}}^{ < }{\sigma }_{\mu }\frac{\partial {g}_{{\boldsymbol{p}}}^{A}}{\partial \varepsilon }-\frac{\partial {g}_{{\boldsymbol{p}}}^{ < }}{\partial \varepsilon }{\sigma }_{\mu }{g}_{{\boldsymbol{p}}}^{A}\}\,{\partial }_{t}{M}_{\mu },\end{array}$$
(6)

where μ = x, y, z, while gR, gA, and g< are the bare retarded, advanced, and lesser Green’s functions respectively. To derive the expression (6) we made use of Fourier transformation \({g}_{{\boldsymbol{p}}}=\int \,d({t}_{1}-{t}_{2}){g}_{{\boldsymbol{p}}}({t}_{1}-{t}_{2})\,{e}^{i\varepsilon ({t}_{1}-{t}_{2})}\) and integration by parts.

To finally close up the derivation we employ the fluctuation-dissipation theorem according to which g< (ε) = [gA(ε) − gR(ε)]f(ε), where f(ε) = [eβ(ε−μ) + 1]−1 stands for the Fermi-Dirac distribution with the Fermi energy μ. Thus, nonequilibrium spin density of conduction electrons (3) within linear response theory is determined by \({Q}_{\mu \nu }={Q}_{\mu \nu }^{(1)}+{Q}_{\mu \nu }^{(2)}\), where

$$\begin{array}{rcl}{Q}_{\mu \nu }^{\mathrm{(1)}} & = & \frac{1}{4}\,{\rm{Tr}}\,[{\sigma }_{\mu }\,{\int }_{-\infty }^{\infty }\,\frac{d\varepsilon }{2\pi }\,\int \,\frac{{d}^{2}p}{{\mathrm{(2}\pi )}^{2}}\{\frac{\partial {g}_{{\boldsymbol{p}}}^{R}}{\partial \varepsilon }{\sigma }_{\nu }{g}_{{\boldsymbol{p}}}^{R}-{g}_{{\boldsymbol{p}}}^{R}{\sigma }_{\nu }\frac{\partial {g}_{{\boldsymbol{p}}}^{R}}{\partial \varepsilon }\\ & & +\,{g}_{{\boldsymbol{p}}}^{A}{\sigma }_{\nu }\frac{\partial {g}_{{\boldsymbol{p}}}^{A}}{\partial \varepsilon }-\frac{\partial {g}_{{\boldsymbol{p}}}^{A}}{\partial \varepsilon }{\sigma }_{\nu }{g}_{{\boldsymbol{p}}}^{A}\}f(\varepsilon )],\end{array}$$
(7)

which involves the integration over the whole Fermi sea, and

$${Q}_{\mu \nu }^{\mathrm{(2)}}=\tfrac{1}{4}\,{\rm{Tr}}\,[{\sigma }_{\mu }\,{\int }_{-\infty }^{\infty }\,\frac{d\varepsilon }{2\pi }\,(\,-\tfrac{\partial f(\varepsilon )}{\partial \varepsilon })\,\int \,\tfrac{{d}^{2}p}{{\mathrm{(2}\pi )}^{2}}\{{g}_{{\boldsymbol{p}}}^{R}{\sigma }_{\nu }{g}_{{\boldsymbol{p}}}^{R}+{g}_{{\boldsymbol{p}}}^{A}{\sigma }_{\nu }{g}_{{\boldsymbol{p}}}^{A}-2{g}_{{\boldsymbol{p}}}^{R}{\sigma }_{\nu }{g}_{{\boldsymbol{p}}}^{A}\}],$$
(8)

which selects the integration in the vicinity of the Fermi level. Generally, the form of Qμν belongs to the class of Kubo-Středa formula, and, in essence, represents the response to the external stimulus in the form of ∂tMν. We can immediately establish a quantitative agreement between the result given by Eq. (8) and the previous studies within a Kubo formalism40,43,44,45,46 which allow a direct estimation within the framework of disordered alloys. Formally, the expression (7) corresponds to the so-called Středa contribution. Such a term was originally identified in ref.47 when studying quantum-mechanical conductivity. Notably, in Eq. (7) each term represents the product of either retarded or advanced Green’s functions. In this case the poles of the integrand function are positioned on the same side of imaginary plane, making disorder correction smaller in the weak disorder limit (see, e.g., ref.48). Meanwhile, having no classical analog this contribution appears to be important enough when the spectrum of the system is gapped and the Fermi energy is placed exactly in the gap47. It is worth mentioning that the contribution due to Eq. (7) has never been discussed in this context before. In the meantime, Kubo-Středa expression for the components of the Gilbert damping tensor has been addressed from the perspective of first-principles calculations49 and current-induced torques50.

Results and Discussion

Let us apply the formalism developed in the previous section to a prototypical model: we work out the Gilbert damping tensor for a Rashba ferromagnet with the magnetization m = z aligned along the z axis. In the limit of weak disorder the Green’s function of a bare system can be expressed as

$${g}_{{\boldsymbol{p}}}^{R}(\varepsilon )={(\varepsilon -{H}_{0}-{{\rm{\Sigma }}}^{R})}^{-1}=\tfrac{\varepsilon -{\varepsilon }_{{\boldsymbol{p}}}+i\delta +\alpha {({\boldsymbol{\sigma }}\times {\boldsymbol{p}})}_{z}+({\rm{\Delta }}+i\eta ){\sigma }_{z}}{{(\varepsilon -{\varepsilon }_{{\boldsymbol{p}}}+i\delta )}^{2}-{\alpha }^{2}{p}^{2}-{({\rm{\Delta }}+i\eta )}^{2}},$$
(9)

where εp = p2/(2 m) is the electron kinetic energy. We put the self-energy ΣR due to scattering off scalar impurities into Eq. (9), which is determined from ΣR = −i(δ − ησz) (see, e.g., ref.51). In particular, for |ε| > |Δ| we can establish that δ = 1/(2τ) and η = 0 in the weak disorder regime to the leading order.

Without loss of generality, in the following we restrict the discussion to the regime μ > |Δ|, which is typically satisfied with high accuracy in experiments. As previously discussed, the contribution owing to the Fermi sea, Eq. (7), can in some cases be ignored, while doing the momentum integral in Eq. (8) results in,

$$\tfrac{1}{m\tau }\,\int \,\frac{{d}^{2}p}{{\mathrm{(2}\pi )}^{2}}{g}_{{\boldsymbol{p}}}^{R}(\varepsilon ){\boldsymbol{\sigma }}{g}_{{\boldsymbol{p}}}^{A}(\varepsilon )=\tfrac{{{\rm{\Delta }}}^{2}}{{{\rm{\Delta }}}^{2}+2\varepsilon \rho }{\boldsymbol{\sigma }}+\tfrac{{\rm{\Delta }}\delta }{{{\rm{\Delta }}}^{2}+2\varepsilon \rho }({\boldsymbol{\sigma }}\times {\boldsymbol{z}})+\tfrac{{{\rm{\Delta }}}^{2}-\varepsilon \rho }{{{\rm{\Delta }}}^{2}+2\varepsilon \rho }({\boldsymbol{\sigma }}\times {\boldsymbol{z}})\times {\boldsymbol{z}},$$
(10)

where ρ = mα2. Thus, thanks to the factor of delta function δ(ε − μ) = −∂f(ε)/∂ε, to estimate \({Q}_{\mu \nu }^{(2)}\) at zero temperature one should put ε = μ in Eq. (10). As a result, we obtain,

$${Q}_{\mu \nu }^{\mathrm{(2)}}=-\,\frac{1}{4\pi }\frac{m}{{{\rm{\Delta }}}^{2}+2\mu \rho }(\begin{array}{ccc}2\tau \mu \rho & {\rm{\Delta }} & 0\\ -{\rm{\Delta }} & 2\tau \mu \rho & 0\\ 0 & 0 & 2\tau {{\rm{\Delta }}}^{2}\end{array}).$$
(11)

Meanwhile, to properly account the correlation functions which appear when averaging over disorder configuration one has to evaluate the so-called vertex corrections, which from a physical viewpoint makes a distinction between disorder averaged product of two Green’s function, 〈gRσνgA〉dis, and the product of two disorder averaged Green’s functions, 〈gR〉disσν〈gA〉dis, in Eq. (8). Thus, we further proceed with identifying the vertex part by collecting the terms linear in δ exclusively,

$${{\boldsymbol{\Gamma }}}^{\sigma }=A{\boldsymbol{\sigma }}+B({\boldsymbol{\sigma }}\times {\boldsymbol{z}})+C({\boldsymbol{\sigma }}\times {\boldsymbol{z}})\times {\boldsymbol{z}},$$
(12)

provided A = 1 + Δ2/(2ερ), B = (Δ2 + 2ερ)Δδ/(Δ2 + ερ)2, and C = Δ2/(2ερ) − ερ/(Δ2 + ερ). To complete our derivation we should replace σν in Eq. (8) by \({{\rm{\Gamma }}}_{\nu }^{\sigma }\) and with the aid of Eq. (10) we finally derive at ε = μ,

$${Q}_{\mu \nu }^{\mathrm{(2)}}=(\begin{array}{ccc}{Q}_{xx} & {Q}_{xy} & 0\\ -{Q}_{xy} & {Q}_{xx} & \\ 0 & 0 & -m\tau {{\rm{\Delta }}}^{2}/\mathrm{(4}\pi \mu \rho )\end{array})\mathrm{.}$$
(13)

We defined Qxx = −mτμρ/[2π(Δ2 + μρ)] and Qxy = −mΔ(Δ2 + 2μρ)/[4π(Δ2 + μρ)2], which unambiguously reveals that account of vertex correction substantially modifies the results of the calculations. With the help of Eqs (3), (11) and (13) we can write down LLG equation. Slight deviation from collinear configurations are determined by x and y components (mx and my respectively, so that \(|{m}_{x}|\), \(| {m}_{y}| \ll 1\)). The expressions (11) and (13) immediately suggest that the Gilbert damping at the interface is a scalar, αG,

$${\partial }_{t}{\boldsymbol{m}}=\tilde{\gamma }{{\boldsymbol{H}}}_{{\rm{eff}}}\times {\boldsymbol{m}}+{\alpha }_{G}{\boldsymbol{m}}\times {\partial }_{t}{\boldsymbol{m}},$$
(14)

where the renormalized gyromagnetic ratio and the damping parameter are,

$$\tilde{\gamma }=\frac{\gamma }{1+\chi {\rm{\Delta }}{Q}_{xy}},{\alpha }_{G}=-\,\frac{\chi {\rm{\Delta }}{Q}_{xx}}{1+\chi {\rm{\Delta }}{Q}_{xy}}\approx -\,\chi {\rm{\Delta }}{Q}_{xx}\mathrm{.}$$
(15)

In the latter case we make use of the fact that \(m\chi \ll 1\) for the NM thickness d ~ 100 μm–100 nm. In Eq. (14) we have redefined the gyromagnetic ratio γ, but we might have renormalized the magnetization instead. From physical perspective, this implies the fraction of conduction electrons which become associated with the localized moment owing to sd exchange interaction. With no vertex correction included one obtains

$${\alpha }_{G}=\frac{m\chi }{2\pi \hslash }\frac{\tau \mu \rho {\rm{\Delta }}}{{{\rm{\Delta }}}^{2}+2\mu \rho },$$
(16)

while taking account of vertex correction gives rise to a different result,

$${\alpha }_{G}=\frac{m\chi }{2\pi \hslash }\frac{\tau \mu \rho {\rm{\Delta }}}{{{\rm{\Delta }}}^{2}+\mu \rho }\mathrm{.}$$
(17)

To provide a quantitative estimate of how large the Středa contribution in the weak disorder limit is, on condition that µ > |Δ|, we work out \({Q}_{\mu \nu }^{(1)}\). Using ∂gR/A(ε)/∂ε = −[gR/A(ε)]2 and the fact that trace is invariant under cyclic permuattaions we conclude that only off-diagonal components μ ≠ ν contribute. While the direct evaluation results in \({Q}_{xy}^{(1)}=3m{\rm{\Delta }}/[2({{\rm{\Delta }}}^{2}+2\mu \rho )]\) in the clean limit. It has been demonstrated that including scattering rates δ and η does not qualitatively change the results, leading to some smearing only52.

Interestingly, within the range of applicability of theory developed in this paper, the results of both Eqs (16) and (17) depend linearly on scattering rate, being thus in qualitative agreement with the breathing Fermi surface model. Meanwhile, the latter does not yield any connection to the microscopic parameters (see, e.g., ref.53 for more details). To provide with some quantitative estimations in our simulations we utilize the following set of parameters. Typically, experimental studies based on hyperfine field measurements equipped with DFT calculations54 reveal the sd Stoner interaction to be of the order of 0.2 eV, while the induced magnetization of s-derived states equals 0.002–0.05 (measured in the units of Bohr magneton, μB). Thus, the parameter of sd exchange splitting, appropriate for our model, is \({\rm{\Delta }}\) ~ 0.2–1 meV. In addition, according to first-principles simulations we choose the Fermi energy μ ~ 3 eV. The results of numerical integration of (8) are presented in Fig. 2 for several choices of sd exchange and scattering rates, τ. The calculations reveal almost no temperature dependence in the region up to room temperature for any choice of parameters, which is associated with the fact that the dominant contribution comes from the integration in a tiny region of the Fermi energy. Figure 2 also reveal a non-negligible dependence on the damping parameter with respect to both Δ and τ, which illustrates that a tailored search for materials with specific damping parameter needs to address both the sd exchange interaction as well as the scattering rate. From the theoretical perspective, the results shown in Fig. 2 correspond to the case of non-interacting electrons with no electron-phonon coupling included. Thus, the thermal effects are accounted only via temperature-induced broadening which does not show up for µ > |Δ|.

Figure 2
figure 2

Gilbert damping, obtained from numerical integration of Eq. (8), shows almost no temperature dependence associated with thermal redistribution of conduction electrons. Dashed lines are plotted for Δ = 1 meV for τ = 1 and τ = 10 ns, whereas solid lines stand for Δ = 0.2, 0.3, and 1 meV for τ = 100 ns.

Full size image

Conclusions

In this paper we proposed an alternative derivation of the Gilbert damping tensor within a generalized Kubo-Středa formula. We established the contribution stemming from Eq. (7) which was missing in the previous analysis within the linear response theory. In spite of being of the order of (μτ)−1 and, thus, negligible in the weak disorder limit developed in the paper, it should be properly worked out when dealing with more complicated systems, e.g., gapped materials such as iron garnets (certain half metallic Heusler compounds). For a model system, represented by a Rashba ferromagnet, we directly evaluated the Gilbert damping parameter and explored its behaviour associated with the temperature-dependent Fermi-Dirac distribution. In essence, the obtained results extend the previous studies within linear response theory and can be further utilized in first-principles calculations. We believe our results will be of interest in the rapidly growing fields of spintronics and magnonics.