Another view on Gilbert damping in two-dimensional ferromagnets.

A keen interest towards technological implications of spin-orbit driven magnetization dynamics requests a proper theoretical description, especially in the context of a microscopic framework, to be developed. Indeed, magnetization dynamics is so far approached within Landau-Lifshitz-Gilbert equation which characterizes torques on magnetization on purely phenomenological grounds. Particularly, spin-orbit coupling does not respect spin conservation, leading thus to angular momentum transfer to lattice and damping as a result. This mechanism is accounted by the Gilbert damping torque which describes relaxation of the magnetization to equilibrium. In this study we work out a microscopic Kubo-Středa formula for the components of the Gilbert damping tensor and apply the elaborated formalism to a two-dimensional Rashba ferromagnet in the weak disorder limit. We show that an exact analytical expression corresponding to the Gilbert damping parameter manifests linear dependence on the scattering rate and retains the constant value up to room temperature when no vibrational degrees of freedom are present in the system. We argue that the methodology developed in this paper can be safely applied to bilayers made of non- and ferromagnetic metals, e.g., CoPt.

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 fields[1-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, (, t), about the effective magnetic field, eff(, t), created by the localized moments in magnetic materials, and its relaxation to equilibrium. The latter, known as the Gilbert damping torque[12], was originally captured in the form α × ∂, where the parameter α determines the relaxation strength, and it was recently shown to originate from a systematic non-relativistic expansion of the Dirac equation[13]. 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 devices[14].

From an experimental viewpoint, the Gilbert damping parameter can be extracted from ferromagnetic resonance linewidth measurements[15-17] or established via time-resolved magneto-optical Kerr effect[18,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 FMs[20-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 model[26,27], while the approach based on scattering matrix analysis has been presented in refs[28,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 model[30]. While a direct estimation of spin-relaxation torque from microscopic theory[31], or from spin-wave spectrum, obtained on the basis of transverse magnetic field susceptibility[32,33], are also possible. It is worth mentioning that the results of first-principles calculations within torque-correlation model[34-38] and linear response formalism[39,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 reversal[41]. 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 (, t). Subject to the function (, 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 (, 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 aswhere first two terms correspond to the Hamiltonian of conduction electrons, on condition that the two-dimensional momentum  = (p, p) = p(cos φ, sin φ) specifies electronic states, m is the free electron mass, α stands for spin-orbit coupling strength, while σ = (σ, σ, σ) is the vector of Pauli matrices. The third term in (1) is responsible for sd exchange interaction with the exchange field (, t) = Δ(, 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()U(′)〉 = (mτ)−1δ( − ′) with the scattering rate τ (we set ħ = 1 throughout the calculations and recover it for the final results).

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.

Subject to the norm of the vector |(, t)| = 1 remains fixed, the magnetization, in broad terms, evolves according to (see, e.g., ref.[42]),where corresponds to so-called spin torques. The first term in describes precession around the effective magnetic field eff 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, (, t). It is worth mentioning that in Eq. (2) the parameter γ is the gyromagnetic ratio, while χ = (gμ/ħ)2μ0/d is related to the electron g–factor (g = 2), the thickness of a nonmagnetic layer d, with μ and μ0 standing for Bohr magneton and vacuum permeability respectively. Knowing the lesser Green’s function, G< (t; t), one can easily evaluate nonequilibrium spin density of conduction electrons induced by slow variation of magnetization orientation,where summation over repeated indexes is assumed (μ, ν = x, y, z). The lesser Green’s function of conduction electrons can be represented as G< = (G − G + G)/2, where G, G, G 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 (, t) ≈  + (t − t0)∂ with  = (, t0) is supposed to be hold with high accuracy,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 . 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,where the integral in time domain is taken along a Keldysh contour, while g(t1, t2) = g(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 , which is needed to evaluate (3), one can write down,where μ = x, y, z, while g, g, and g< are the bare retarded, advanced, and lesser Green’s functions respectively. To derive the expression (6) we made use of Fourier transformation and integration by parts.

To finally close up the derivation we employ the fluctuation-dissipation theorem according to which g< (ε) = [g(ε) − g(ε)]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 , wherewhich involves the integration over the whole Fermi sea, andwhich 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 ∂M. We can immediately establish a quantitative agreement between the result given by Eq. (8) and the previous studies within a Kubo formalism[40,43-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 gap[47]. 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 calculations[49] and current-induced torques[50].

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  =  aligned along the z axis. In the limit of weak disorder the Green’s function of a bare system can be expressed aswhere ε = p2/(2 m) is the electron kinetic energy. We put the self-energy Σ due to scattering off scalar impurities into Eq. (9), which is determined from Σ = −i(δ − ησ) (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,where ρ = mα2. Thus, thanks to the factor of delta function δ(ε − μ) = −∂f(ε)/∂ε, to estimate at zero temperature one should put ε = μ in Eq. (10). As a result, we obtain,

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, 〈gσg〉dis, and the product of two disorder averaged Green’s functions, 〈g〉disσ〈g〉dis, in Eq. (8). Thus, we further proceed with identifying the vertex part by collecting the terms linear in δ exclusively,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 and with the aid of Eq. (10) we finally derive at ε = μ,

We defined Q = −mτμρ/[2π(Δ2 + μρ)] and Q = −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 (m and m respectively, so that , ). The expressions (11) and (13) immediately suggest that the Gilbert damping at the interface is a scalar, α,where the renormalized gyromagnetic ratio and the damping parameter are,

In the latter case we make use of the fact that 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 obtainswhile taking account of vertex correction gives rise to a different result,

To provide a quantitative estimate of how large the Středa contribution in the weak disorder limit is, on condition that µ > |Δ|, we work out . Using ∂g(ε)/∂ε = −[g(ε)]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 in the clean limit. It has been demonstrated that including scattering rates δ and η does not qualitatively change the results, leading to some smearing only[52].

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 calculations[54] 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, μ). Thus, the parameter of sd exchange splitting, appropriate for our model, is  ~ 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

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.

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.