Calculation of Gilbert damping in ferromagnetic films

The Gilbert damping constant in the phenomenological Landau-Lifshitz-Gilbert equation which describes the dynamics of magnetization, is calculated for Fe, Co and Ni bulk ferromagnets, Co films and Co/Pd bilayers within a nine-band tight-binding model with spin-orbit coupling included. The calculational efficiency is remarkably improved by introducing finite temperature into the electronic occupation factors and subsequent summation over the Matsubara frequencies. The calculated dependence of Gilbert damping constant on scattering rate for bulk Fe, Co and Ni is in good agreement with the results of previous ab initio calculations. Calculations are reported for ferromagnetic Co metallic films and Co/Pd bilayers. The dependence of the Gilbert damping constant on Co film thickness, for various scattering rates, is studied and compared with recent experiments.


Introduction
The switching of magnetization in nanoscale magnetic elements is of importance for application in magnetoelectronic devices.The switching process, done with external field or current induced spin-transfer [1] is substantially affected by magnetic damping phenomena present in such systems.The following three mechanisms are usually held responsible for the relaxation of the magnetisation direction in metallic ferromagnets: Gilbert damping [2,3], spin pumping [4,5] and two-magnon scattering [6,7].The first of them originates from the spin-orbit (SO) interaction which couples the spins of electrons to the lattice, and thus leads to a torque exerted on the magnetisation as its direction varies with time.This torque is represented by the Gilbert damping term in the phenomenological Landau-Lifshitz-Gilbert (LLG) equation and it plays an important role in spin dynamics of ferromagnetic metallic systems.In particular, it affects domain wall (DW) velocity in currentcarrying ferromagnetic structures/nanostripes where fast propagation of domain walls is of great importance for applications in spintronic devices, such as racetrack memories [8].Magnetic layered structures continue to attract much attention in the context of such applications as they exhibit low critical current for depinning domain walls [9].Theory predicts an enhanced Gilbert damping constant for magnetic layered systems [10,11] and this is observed experimentally for a Pt/Co/Pt trilayer structure [12].In Ref. [13], an increase in Gilbert damping in [Co/Pt] N multilayer films with increasing stacking number N is found using the time-resolved magneto-optical Kerr effect.Magnetic anisotropy, also related to the SO coupling, is another important factor for DW structures and it can be particularly large in multilayers, like Co/Pt or Co/Ni (see, e.g., [14,15]).
In this communication, we calculate the Gilbert damping constant within Kamberský's spin-orbital torque correlation theory using a realistic nine-band tight binding (TB) a e-mail: ebarati@ichf.edu.plmodel [16,17] for a range of electron scattering rate.The calculations have been performed for bulk Fe, Co and Ni, for ferromagnetic Co films and for Co/Pd bilayers.The dependence of the Gilbert damping constant on Co film thickness and the effect of a Pd overlayer is studied and compared with recent experiments.

Formalism
The magnetization dynamics of a metallic nanowire system can be described by an extended phenomenological LLG equation, where m = M/M s is the magnetization vector M normalized to its saturation value M s , γ is the gyroscopic ratio and A denotes the exchange stiffness constant.The first term in Eq. ( 1) describes Larmor precession of magnetization around the effective magnetic field H eff , applied externally and/or arising due to magnetic anisotropy.It is followed by the relaxation term proportional the dimensionless Gilbert damping constant α.The remaining terms in Eq. ( 1) are effective only for non-homogenous magnetization m(x) (varying with the position x along the nanowire) as in spin waves or DWs.The fourth and the fifth terms, introduced by Zhang and Li [18], describe adiabatic and non-adiabatic spin transfer torques, respectively.They both arise due to the flow of the electric current J, with finite spin polarization P in ferromagnets, and are defined by the parameter v 0 = gµ B PJ/(2eM s ) [16] where g is the Landé factor, e electric charge, and µ B the Bohr magneton.It can be shown [19] that uniform DW motion can be driven by a current with velocity v = (β/α)v 0 .Both α and β are therefore important and both depend on SO coupling.In this EPJ Web of Conferences paper, however, we concentrate on α.Other terms can be included in a further extension of the LLG equation [20].
In the present work, the Gilbert damping constant α is calculated with the modified Kamberský formula [2] which relates the damping of the precessing magnetization to transitions between electronic states due to the torque induced by spin-orbit coupling.Kamberský's expression for the damping constant can be derived as α = ω 1 /ω 0 where ω = ω 0 + iω 1 is the complex frequency of a long wavelength spin wave damped by the SO coupling (see e.g.Ref. [20]).The Gilbert damping constant derived with such an approach for ferromagnetic films and Co/nonmagnet bilayers [21] is given by the formula where the integration over k = (k x , k y ) is done over the two-dimensional Brillouin zone (BZ), a 2d is the (001) fcc surface square lattice constant, µ s denotes the atomic magnetic moment in units of µ B and N Co is the number of atomic layers in the ferromagnet.The matrix elements include spin-orbit torque [S − , H SO ], where S − = S x − iS y and H SO , are the spin operator ( S x,y = (1/2)σ Pauli x,y ) and the SO interaction, respectively.The magnetization is in the z direction, which is assumed to be perpendicular to the film plane in the present calculations, due to either intrinsic magnetic anisotropy and/or an external magnetic field.The factor F nn is found by integrating the product of two Lorentzian functions depending on eigenenergies n , n , with band indexes n, n The presence of the negative derivative η( ) of the Fermi-Dirac distribution function confines the integration over to the vicinity of F .Thus, damping is affected solely by states close to the Fermi level F .The actual range of contributing states also depends strongly on the value Γ of the full width at half maximum (FWHM) of the Lorentzians which represents the electron scattering due to defects and electron-phonon interaction.The numerical calculation of the above integral over is done by employing the residue theorem and the Matsubara frequencies, which are the poles of the Fermi-Dirac distribution function and its derivative.
The expression for α, Eq. ( 2), accounts for two kinds of electronic transitions: intraband (n = n ) and interband (n n ), i.e., within a single energy band and between different energy bands.The former contributes a conductivitylike term (decreasing with Γ), whilst the latter gives a resistivity like term (increasing with Γ) [2,22,23].The Gilbert damping constant is calculated for metallic layered systems in the slab geometry, namely Co films and Co/Pd bilayers.The calculations are performed using a realistic nine-band tight-binding model for a range of the electron scattering rate Γ which corresponds to the inverse lifetime /τ used in [22,24].The TB Hamiltonian for layered systems [17] is constructed with two-center hopping integrals (between the first-and the second-nearest neighbours) obtained by fitting to the ab initio bulk energy bands for bulk metals [25].The SO coupling term H SO is included in the Hamiltonian which is found to be vital for the intraband terms (n = n ) to give significant contribution to the Gilbert damping [21].For the assumed perpendicular orientation of magnetisation, the numerical integral in the k-space can be done over 1/8 BZ, with proper weights at the boundary points.
The temperature dependence of the damping constant given by Eq. ( 2) emerges basically due to two factors: the derivative η of the Fermi function and the scattering rate Γ, which increases as the temperature rises because of electronphonon scattering.However, the scattering rate is treated as an independent parameter, like in previous works [2,22].Thus, the explicit dependence on temperature is only through η and it is found to be very weak.This is convenient since the convergence of the numerical integration over BZ is much faster for finite temperature T .
It is found that using T = 300 K the convergence of the Gilbert damping constant for bulk Fe, Co and Ni is reached with as few as (60) 3 k points in the full three-dimensional BZ.This number is much smaller than (100) 3 k (or more) points used in the previously reported calculations [22,23].

Bulk metals
The Gilbert damping constant α calculated for bulk bcc Fe, fcc Co and fcc Ni within the TB model is shown in Fig. 1.It is found to have a very similar dependence on the scattering rate Γ to that obtained for the three bulk ferromagnets in ab initio calculations by Gilmore et al [22].In particular, the characteristic minimum of α occurs for Γ close to 0.1 eV due to the dominance of intra-and interband contributions to α at small and large Γ, respectively.

Ferromagnetic films
The Gilbert damping constant obtained for ferromagnetic (001) fcc Co films is larger than for bulk Co and increases with decreasing film thickness as seen in Fig. 2. The oscillatory behavior of the damping constant originates mainly from interband term and it is attributed to quantum-well states (QWS) with energies close to the Fermi level F .The occurrence of QWS in metallic films is a well-known phenomenon which leads to oscillations of interlayer exchange coupling and magnetic anisotropy with varying thicknesses of ferromagnetic films and/or non-magnetic layers [26][27][28].

Ferromagnet/nonmagnet bilayers
The Gilbert damping constant has also been calculated for Co/Pd bilayer systems with different thicknesses of Co.The expression for Gilbert damping in bilayer systems depends on the SO coupling constants of both ferromagnetic and nonmagnetic metals when the SO interaction is represented in the atomic orbital basis in the TB model.The values of the SO coupling constant for the d electrons in these metals are taken from Ref. [29].The Gilbert damping constant α calculated for the Co/Pd(2 ML) bilayer is shown for three different scattering rates Γ in Fig. 3.The dependence of α on the Co thickness follows the approximate relation α ≈ α b + α s /N Co (5) which implies that the total Gilbert damping constant in Co/Pd bilayers consists of the bulk-like contribution α b from the ferromagnet, and α s which the sum of the contributions coming from the Co/vacuum and Co/Pd interfaces present in the system.The approximate relation ( 5) is also found to hold for the Co film (Fig. 2).A linear dependence of α on 1/N Co has recently been observed experimentally for Pt/Co/Pt trilayers [12].
As seen in Fig. 3, the Gilbert damping constant is remarkably enhanced in the Co/Pd bilayers in comparison with the Co film, particularly strongly for Co thicknesses less than 5 ML.This enhancement comes from the large SO coupling in the Pd overlayer which causes the interface term α s to be an order of magnitude larger than the Co bulk term α b and thus dominate for small Co thicknesses.The hybridization at Co/Pd interface has been shown [21] to be less important for increasing magnetic damping in Co/Pd systems.Oscillations are also present in these bilayers but the details are modified by the presence of Pd.They are not visible in Fig. 3 because they are relatively small in comparison to the average value of α.

Breakdown of Gilbert damping constant in k -space
Expression (2) for the Gilbert damping constant has the form of the integral over the BZ where each contribution g(k) is a double sum over band indexes n, n and includes all coefficients appearing in Eq.
(2).Thus, it is found how different k points in the BZ contribute to the Gilbert damping constant.The results obtained for an (001) fcc Co(8 ML) film (Fig. 4) show that there are "hot" regions in the k-space distribution of the Gilbert damping constant.The regions where the contribution g(k) has largest values (peaks) are located at points where there are pairs of eigenstates with energies close to Fermi level, as is clear from Eq. ( 4).These are located mostly near the corners of two-dimensional BZ.

Conclusions
Within EPJ Web of Conferences film thickness are in evidence.These originate largely from the interband term and we attribute these oscillations to QWS.The Gilbert damping constant for Co/Pd bilayers is found to be strongly affected by the presence of the nonmagnetic overlayer due to a large interface contribution.

Fig. 4 .
Fig. 4. k points contributions to the Gilbert damping constant in slab Co(8 ML) and with the scattering rate Γ = 0.01.
a realistic tight-binding formalism and based on Kamberský's torque correlation model new results for the Gilbert damping constant in bulk bcc Fe, fcc Co, and fcc Ni ferromagnetic films, and Co/Pd bilayers are presented.The results obtained for bulk ferromagnets show good agreement with previous ab initio calculations.The results for Co films and Co/Pd bilayers show enhancement of Gilbert damping with respect to bulk Co and an increase of the damping constant, depending on scattering rate, with decreasing layer thickness.Furthermore, characteristic oscillations of the Gilbert damping constant as a function of 18003-p.3