Simulation study of the Gilbert damping in Ni80Fe20/Nd bilayers: comparison with experiments

We present an experimental and computational investigation the Neodymium thickness dependence of the effective damping constant ( αeff ) in Ni 80 Fe 20/Neodymium (Py/Nd) bilayers. The computational results show that the magnetic damping is strongly dependent on the thickness of Nd, which is in agreement with experimental data. Self consistent solutions of the spin accumulation model and the local magnetisation were used in the simulations. It was not possible to obtain agreement with experiment under the assumption of an enhanced damping in a single Py monolayer. Instead, it was found that the enhanced damping due to spin pumping needed to be spread across two monolayers of Py. This is suggested to arise from interface mixing. Subsequently, the temperature dependence of the effective damping was investigated. It is found that, with increasing temperature, the influence of thermally-induced spin fluctuations on magnetic damping becomes stronger with increasing Nd thickness.


Introduction
Recently, spintronic devices based on magnetic dynamics, such as magnetic random-access memory, spin valve (SV), Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.magnetic tunnel junction and so on, have attracted extensive attention [1][2][3].The Gilbert damping parameter α as a characteristic of magnetisation dynamics is a major focus of research due to its effect on the speed of the magnetisation reversal.As a typical transition magnetic (TM) alloy, Permalloy has a high saturation magnetisation and permeability, low anisotropy and coercivity.Permalloy (Ni 80 Fe 20 , Py) has a small magnetic damping, which is due to the fact that the crystal field will quench the most of orbital momentum, resulting in a small spin orbit coupling.So, it is necessary to increase the magnetic damping to increase the magnetic response.There are mainly two methods to regulate magnetic damping, i.e. doping and magnetic heterostructures [4][5][6][7][8].Rare earth elements (RE) have a large orbital moment, which could increase the spin-orbit coupling and magnetic damping in RE-TM systems.In 2002, Tserkovnyak et al [9,10] proposed the spin pumping effect: in ferromagnetic/nonmagnetic (FM/NM) bilayer, the magnetisation precession acts as spin pump which transfers angular momentum from the FM layer into the NM layer which is then dissipated in the NM layer.This dissipation increases the relaxation of the magnetic moment precession in the FM layer and leads to an increase of magnetic damping.The enhanced magnetic damping is represented by the scattering matrix of the FM layer.
In this paper, we present experimental and theoretical investigations of the Nd thickness dependence of the effective damping constant (α eff ) in Ni 80 Fe 20 /Neodymium (Py/Nd) bilayers.In experiments, the magnetic damping can be easily obtained by vector network analyser ferromagnetic resonance (VNA-FMR) in comparison with Brillouin light scattering and time-resolved magneto-optical Kerr effect.VNA-FMR is a standard and widely used tool to measure the magnetic damping by analysing the relationship between frequency and linewidth of the resonance spectrum.The FMR linewidth has two contributions: one is the intrinsic damping coming from the Py film; the other arises from a spin pumping effect which pumps a spin current from the Py layer to Nd the layer.
In experiments, there are a number of reasons for enhanced magnetic damping, including the contribution of enhanced spin orbit coupling at the interface, and dispersions in material properties which results in inhomogeneous line broadening.The important mechanism of the spin pumping effect is rarely considered, whereas spin pumping, and especially its dependence on the thickness of a non-magnetic layer, could be considered an interesting method of damping enhancement and control.So, we use the spin accumulation methodology to simulate the bilayer structures based on the driftdiffusion formalism of Zhang et al [11].Based on previous works, we introduce the bulk and interface magnetic damping to investigate the thickness dependence in Py/Nd bilayer films [12][13][14].So, the appropriate values of the bulk and interface damping are essential for a theoretical calculation for comparison with the experimental data.Py itself has a very low anisotropy.So, we set the same anisotropy energy values for bulk Py layers.Then the high interface magnetic damping values in interface Py comes from the spin pumping effect [15].The simulation of our bilayer structure is modelled by using the stochastic Landau-Lifshitz-Gilbert (LLG) equation at the atomistic level and the VAMPIRE software package [16,17].Meanwhile, with increasing temperature, the other scattering mechanism, i.e. thermal-spin fluctuations leading to spin-wave excitations are naturally taken into account in the atomistic approach.In order to study the influence of thermal-spin fluctuations in the Py/Nd bilayer, the temperature variation of the effective Gilbert damping is also studied.

The spin accumulation model
The concept of spin injection into a normal metal leading to a spin accumulation of electron polarisation was put forward by Aronov [18] and has been enriched and developed by many research groups [11,19,20].This spin accumulation decays on a length scale characterised by the spin difusion length (sdl).According to the formula derived by Zhang et al [11], in a magnetic heterostructure the equation of motion for the spin accumulation is, where m is the spin accumulation, M d is the unit vector for the local magnetisation of the FM layer, J sd is the exchange energy between the electron spin and the local magnetisation, h is the reduced Planck constant, τ sf is the spin-flip relaxation time of the conduction electrons.In order to obtain a stationary solution for spin accumulation m, the spin-flip relaxation time of the conduction electrons τ sf should be shorter than the timescale of the magnetisation changes.
The equation of motion for the local magnetisation is where γ is the gyromagnetic ratio, α is the magnetic damping constant, H eff is the effective field, which includes the external field, exchange, magnetostatic and crystalline anisotropy field.J sd m is the field arising from the conduction electron/magnetic spin exchange.
In general, the magnetisation direction and the spin accumulation are not collinear.The spin accumulation be separated into two components: longitudinal mode (parallel to the local magnetisation direction) and transverse mode (perpendicular to the local magnetisation direction).So, the two components spin accumulation can be expressed as and where In the general case the solution for the spin accumulation can be obtained using a transfer matrix method [19].However, according to Heide et al [21] only the longitudinal accumulation exists for a single FM layer, while the transverse accumulation is required for a 'SV' structure.That means we can simply consider the longitudinal component of spin accumulation in our structure.The current is perpendicular to the plane of the films.As shown in figure 1(a), the pumped spin accumulation whose magnetisation direction is at the positive z direction is given by the following equation, (5) where m 0 = λ sdl 4D0 j mx (0) is the pre-exponential factor for the pumped spin accumulation.j mx (0) is the spin current transmitted from the FM to the NM, given by [9,10] with g ↑↓ r the spin mixing conductance.The reflected spin accumulation at the vacuum boundary of the non-magnet is given by where A is the pre-exponential factor for the reflected spin accumulation, d NM is the thickness of NM layer.So, the spin accumulation in the steady state as follows, Then the spin current density in x direction calculated from the spin accumulation is given by At the right vacuum boundary of the NM layer, we set the current to zero, j m (d) = 0, Then we obtain the pre-exponential factor: A = m0 λ sdl e d NM λ sdl .Hence, and at the FM/NM interface (x = 0), the net spin current is given by, According to spin pumping theory, the effective magnetic damping including the intrinsic damping of FM layer and the net spin current at the FM/NM interface.So, the interface damping of the FM layer is given by, where α max is the enhanced magnetic damping and is proportional to the interface spin-mixing conductance G ↑↓ .Therefore, the interface damping increases to an asymptotic value with the thickness of NM as shown in figure 1(b).For larger thicknesses, the increase in damping took the interfacial value to supercritical values (α > 1) which led to fictitious dynamics.Therefore, we distribute the interfacial damping among the two Py layers nearest to the interface as the spin-pumping current is an interface effect.This is an interesting observation which suggests some intermixing at the Py/Nd interface.
The magnetization dynamics was simulated using an atomistic model via integration of the LLG equation.
The effective field, B eff,i = − 1 µs ∂H ∂Si , can be calculated from the Heisenberg spin Hamiltonian describing the energy contributions of the magnetic system given by: where S ij is the unit vector of spin on site (i, j), J ij is the nearest neighbour exchange integral between spin sites i and j, k u is the uniaxial anisotropy constant, γ is the absolute value of gyromagnetic ratio, α is the damping constant of the material.The terms on the right-hand-side of equation ( 15) represent the contribution of the exchange energy, the anisotropy energy and the Zeeman energy, respectively.Atomistic calculations were carried out using the VAMPIRE package [17].

The thickness dependence of effective damping
The Py (6 nm)/Nd (t nm) structure is built as a bulk facecentred-cubic crystal with lattice constant of 3.54 Å.The structure, with dimensions of 20 × 20 × (6 + t) nm 3 , is considered as a trilayer system consisting of a 2-monolayers-Py (interface-Py) in contact with the bulk-Py and varying thickness of Nd layer as shown in figure 2(a).To verify the correctness of the magnetic parameters of the trilayer structure and the method, the thickness dependence of the effective damping at 300 K is calculated first and then compared with the experimental work (as shown in supplementary S1).The magnetic parameters of the trilayer are shown in table 1.Here we introduce an anisotropy and exchange energy [22][23][24].The interface magnetic damping of the Py in trilayer system obtained from the experimental data using equation ( 13) in this simulation are shown in figure 2(b), and the interface magnetic damping of the Py are shown in supplementary S2.Nd is a light RE material with antiferromagnetic transition temperature of 19.3 K [25,26].It shows that Nd exhibits paramagnetism at room temperature.
In simulations, the external field (10 T) is applied along the z direction.The magnetization is initially oriented in the yz plane at an angle of 30 deg to the applied field (z-direction) from where it undergoes damped precession into the field direction.A typical time evolution of the x and z component of the magnetisation precession obtained from the atomistic model are shown in figure 3(a).For different thickness of Nd, the precession frequency of the Py layer and the relaxation time are extracted, as the net magnetisation of the system evolves back to equilibrium state, using the following equation, where A, f, τ are the amplitude of the magnetisation, the precession frequency and the relaxation time, respectively.The effective magnetic damping (α eff ) can be determined from the damped precession parameters as α eff = 1 2π fτ .Using this approach, the precession frequency and damping were calculated as a function of the Nd layer thickness.
By comparing the experimental data with the simulations, for zero Nd thickness, as shown in figure 3(b), the damping value is equal to the input Py value, α eff = 0.0099.Further, this value is consistent with our experimental result and other  groups [13,27].This result shows that the calculated damping value of the bulk-Py is reasonable.By varying the interface magnetic damping according to the thickness of Nd, we also obtained the effective damping values as a function of Nd thickness.Comparison with the experimental data gives good agreement as shown in figure 3

The temperature dependence of the effective magnetic damping
Further, we investigated the thickness dependence of the effective magnetic damping as a function of the temperature.
In order to obtain statistically accurate results, we performed 10 different integrator-random-seeds for each parameter set, and then took the average and variance to obtained the mean magnetic damping.Following [16], the thermal-spin fluctuations are introduced via a random field calculated as, where Γ(t) is Gaussian distribution, T is the temperature in kelvin.k B is the Boltzmann constant, µ s is the atomic magnetic Figure 4(a) shows the temperature dependence of the magnetic damping at different Nd thickness.From the data it can be seen that the temperature has a weak influence on the mean effective magnetic damping as also noted by Sampan-a-pai et al [22] for a single layer of CoFeB well below T c .Interestingly, similarly to [22] we also observe statistical scatter in the calculated damping, albeit at a lower level.The standard deviation increases with Nd thickness, as shown in figure 4(b), suggesting that the effect is proportional to the damping.This is likely due to statistical variations arising from the finite size of the sample as in [22] and will decrease with increasing system size.To further investigate the influence of the temperature on effective magnetic damping, the magnetic dynamic relaxation process at different temperatures are shown in figure 5 for a Py (6 nm)/Nd (16 nm) bilayer.As shown in figures 5(a) and (b), (for the same random seed), the magnetisation precession response differs at different temperatures.At low temperature (20 K), as shown in figure 5(a), the influence of the temperature on the simulation is very small and the fitting curve is consistent for the whole time series.However, at high temperature (500 K), as shown in figure 5(b), the thermal-spin fluctuation plays a role on relaxation time.While the fitted damping is very similar in both cases, at the higher temperatures the dynamics at longer times becomes driven by the response of the system to the random thermal fields.

Conclusions
In summary, we have investigated the thickness and temperature dependence of the effective magnetic damping in Py (6 nm)/Nd (t nm) bilayer by using the spin accumulation model in an atomistic simulation.The approach used was to simulate the precession to equilibrium after rotation through a small angle.The calculation used an atomistic model to calculate the spin pumping into the non-magnetic capping layer.For large Nd layer thickness, the (supercritical) damping, when associated with a single Py atomic layer, was found to give rise to discrepancies with experiment and it was found necessary to spread the elevated damping of two Py layers.This problem and its solution is an interesting prediction of the atomistic model and probably arises from interface diffusion.Therefore, the calculations were carried out on a system of Py (bulk)/Py (2-monolayer)/Nd (t) layers.Through a series of calculations, we determined the damping value of the bulk-Py and the 2-monolayer-Py, with the Nd as a nonmagnetic material.The thickness dependence of the effective magnetic damping value from the simulations is consistent with the experimental data.Further, the effect of temperature on the mean magnetic damping was also carefully investigated.With increasing temperature, the influence of thermal fluctuations on mean magnetic damping increases, which is probably a finite size effect.Our results show that the combination of atomistic and spin accumulation models is a powerful approach to the investigation of spin pumping in heterostructures, especially where interfacial mixing is likely.

Figure 1 .
Figure 1.(a) Schematic view of FM/NM bilayer.The dark blue and the orange represent the bulk FM layer and interface FM layer, respectively.The light blue represents the NM layer.Precession of the magnetisation direction m(t) of the FM pumps a spin current into the NM layer.m+ and m− are pumped and reflected spin accumulation in NM layer, respectively.(b) The curve of effective Gilbert damping (α eff = α 0 + αmax) comes from spin pumping theory.α 0 is the intrinsic damping.αmax is the maximum Gilbert damping increase.

Figure 2 .
Figure 2. (a) Schematic view of the simulated (6 nm) Py/(t nm) Nd bilayer.The structure dimensions are 20 × 20 × (6 + t) nm 3 , where the underlayer (red represent Ni atoms and blue represent the Fe atoms) and the orange inserted layer represents the bulk and interface Py layer, respectively.The grey square represents the Nd layer.(b) The variation of the interface damping and the thickness of Nd.The lines are a guide to the eye.
(b) indicating the validity of the atomistic model with spin accumulation and the input parameters.

Figure 3 .
Figure 3. (a) The time evolution of the magnetisation trace of the Py (6 nm)/Nd (1 nm) structure at 300 K after an initial magnetisation angle of 30 • in the y-z plane back to equilibrium state.The circle represents the simulated magnetisation precession while the lines are the fitting curve using equation (16).(b) The variation of the Gilbert damping and thickness of Nd.The lines are guided to the eye.

Figure 4 .
Figure 4. (a) is the mean magnetic damping of Py (6 nm)/Nd (t nm) bilayers with different temperature and Nd thickness.The lines are guided to lines.(b) The temperature dependence of the standard deviation (σ) of the distribution of the simulated magnetic damping constants as a function of the film thickness.The lines are guided to the eye.

Table 1 .
Model parameters used in the simulation.