Temperature in micromagnetism: cell size and scaling effects of the stochastic Landau–Lifshitz equation

The movement of the macroscopic magnetic moment in ferromagnetic systems can be described by the Landau–Lifshitz (LL) or Landau–Lifshitz-Gilbert (LLG) equation. These equations are strictly valid only at absolute zero temperature. To include temperature effects a stochastic version of the LL or LLG equation for a spin density of one per unit cell can be used instead. To apply the stochastic LL to micromagnetic simulations, where the spin density per unit cell is generally higher, a conversion regarding simulation cell size and temperature has to be established. Based on energetic considerations, a conversion for ferromagnetic bulk and thin film systems is proposed. The conversion is tested in micromagnetic simulations which are performed with the Object Oriented Micromagnetic Framework (OOMMF). The Curie temperatures of bulk Nickel, Cobalt and Iron systems as well as Nickel thin-film systems with thicknesses between 6.3 mono layer (ML) and 31 ML are determined from micromagnetic simulations. The results show a good agreement with experimentally determined Curie temperatures of bulk and thin film systems when temperature scaling is performed according to the presented model.


Introduction
The classical Landau-Lifshitz equation (LL) describes the precession of a macroscopic magnetic moment (M) in a ferromagnet around the effective magnetic field (H eff ) [1,2]. () The precession of the magnetic moment (M) with time (t) around an effective magnetic field (H eff ) is governed by the gyromagnetic ratio (γ) and the damping parameter (α). The effective magnetic field includes contributions from external fields, exchange interaction, anisotropies, etc. The phenomenological damping parameter enables the system to transfer energy and angular momentum from the magnetic movement to other degrees of freedom. The LL is mathematically equivalent to the Landau-Lifshitz-Gilbert (LLG) [3] equation and only differs by the relation of the two parameters, α and γ [2]. The LL and LLG leave the length of the magnetisation vector unchanged and do not include temperature effects [4]. Therefore they are strictly valid only at the absolute zero temperature. To include effects of elevated temperatures on magnetisation dynamics, which are of importance for various research topics and effects, such as temperature gradients, spin-torque effects, laser induced change of magnetisation dynamics, all-optical switching and heat assisted magnetic recording, the LL or LLG have to be modified [2,[5][6][7]. Different approaches exist to include these temperature effects in micromagnetic simulations for elevated temperatures. For constant temperatures, theoretically or empirically determined values of the material parameters can be adjusted according to their temperature dependence and then included as input parameters into the simulation [8]. A second possibility is to use the Landau-Lifhshitz Bloch equation which was derived from the stochastic form of the Landau-Lifhshitz (sLL) equation by a mean-field approximation [2,4,9]. One of the most fundamental approaches is, to include the effects of finite temperature on a microscopic level, as it is done for the sLL. This can be achieved by transforming the LL or LLG into a stochastic differential equation of Langevin (sLL) type [10]. Here, the interaction of an isolated spin with a thermal bath is modelled by addition of a temperature dependent, stochastic Langevin field to the effective field H eff [2,4,5]. This temperature dependent field points for each time step and elementary cell into a random direction. The spatial varying fluctuation leads to a similarly varying precession of the magnetic moment around the effective field direction. By averaging over a huge number of spins, this random perturbation results in a decrease of the total magnetisation, compared to the unperturbed case at zero temperature. The average of many of these interacting spins leads to a temperature dependence of the macroscopic magnetic moment in accordance with spin wave theory [2]. Due to the high number of spins in realistic ferromagnetic structures, this approach based on one spin per unit cell results in high computational costs during micromagnetic simulation, and is therefore only applicable to rather small systems [2]. To apply the sLL to micromagnetic simulations, where the spin density per unit cell is generally higher, a conversion has to be established. To express the problem with the words of the author of the Langevin extension (thetaevolve package) for OOMMF [11]: 'For accurate numerical results (e.g. when trying to determine the Curie-temperature of a system) the change of saturation magnetisation with temperature has to be respected. It is not yet clear exactly how this is done best, but it is clear that it must be done in some way. Otherwise good results can only be expected for a density of just a single spin per cell' [12].
In the following we will derive such a conversion for bulk and two dimensional thin film systems from energetic considerations at the Curie temperature. The conversion will be applied in temperature dependent simulations with OOMMF. Nickel, Cobalt and Iron bulk and Nickel thin film systems with varying cell sizes are simulated. It will be tested if, and under which conditions, the given conversion has the power to estimate the Curie temperatures of the different systems, which is not possible in standard micromagnetism.

Temperature scaling
To determine the scaling between the effective physical temperature (T eff ) and the input parameter used as simulation temperature (T sim ) in dependence of the lattice constant (a eff ) and the length of a elementary simulation cell (a sim ) we follow arguments as given in chapter VII of [13]: above the Curie temperature (T C ) the ferromagnetic behaviour of the respective material vanishes, because the energy originating from thermal excitation overcomes the exchange interaction, which favours magnetic ordering. With other words: for a cubic crystal at the Curie temperature

( )
With k B as Boltzmann's constant, A as the exchange stiffness, M as the magnetisation and a as the characteristic length of the system. In a real world system, the characteristic length is identified with the lattice constant. In micromagnetic simulations the characteristic length is given by the edge of the elementary simulation cell. The magnetisation and exchange stiffness are independent of the characteristic length of the system, due to their definition as magnetic moment density and energy per length, respectively. Thus, following equation (2) the relation between temperature and length of two systems (simulation (sim) and effective/physical (eff) parameter) can be expressed as Hence, the temperature T sim as used in the simulation as input parameter can be determined from the physical temperature T eff and the respective lattice constant a eff and simulation cell length a sim as follows: To determine the range where scaling can be applied, one has to consider the temperature effects on the exchange length of the system, which depends on various simulation parameters [14][15][16] as will be discussed in the following.

Thermal exchange length
To obtain meaningful results from micromagnetic simulations, the deviation of the direction of two magnetisation vectors of neighbouring simulation cells has to be small. This can be achieved by decreasing the size of the simulation cells. At T=0 K for systems where the exchange interaction dominates over other anisotropies, an upper limit for the reasonable size of a simulational cell is given by the exchange length λ ex [17]: Here, μ 0 is the magnetic permeability of the free space, and M S the saturation magnetisation. The exchange length gives the distance over which the above mentioned condition of small deviations of the directions of two neighbouring magnetic moments is fulfilled. The exchange lengths of all simulated materials are in the range of nanometers (see table 1). When temperatures are above 0 K, thermal excitation's and the resulting disorder have to be considered as well. Therefore, the length of a simulation cell a sim has to be smaller than the length over which thermal fluctuations decay in space [14,15]: Whereby the characteristic length λ thex is called thermal exchange length. The interplay of simulation cell size with thermal effects was investigated by different authors [14][15][16]. For example, Tsiantos et al [14] defined λ thex in analogy to equation (5) as and the simulation time step Δt. By replacing T sim via equation (4) with T eff the dependence of Tsiantos thermal exchange length on the simulational cell length a sim , T eff , and Δt can be expressed as For given material parameters the condition of equation (6) can be strictly fulfilled only by a certain set of temperatures, and simulational parameters. Especially when it is considered that the upper limit for the simulation parameter Δt has to be chosen, such that for a given set of fixed parameters the absolute value of the equilibrium magnetisation converges, when Δt is decreased further. To asses if the given set of applied material parameter for bulk materials (table 1) allows a simulational cell size to be chosen, such that it fulfils equation (6), the thermal exchange length for different Δt and a sim (figure 1 left) and T eff (figure 1 right) was calculated. It can be seen, that especially at elevated temperatures, it cannot be avoided that the thermal exchange length becomes smaller than the simulational cell size. Nevertheless, according to Martinez et al [15] λ thex as defined by Tsiantos et al [14], does not necessarily represent the true scale over which thermally introduced fluctuations within a ferromagnet decay in space [15]. With other words, solutions may converge even if equation (6) is not strictly fulfilled. Table 1. Overview over the simulation parameters for bulk materials. Shown are from left to right, the material and its crystal structure, the length of the elementary cell, the magnetisation, anisotropy constant, exchange constant, the simulated, scaled and experimental Curie temperatures, the exchange length and the thermal exchange length at T C exp . The sources of the experimental data are given next to the respective value. All other simulation parameters are given throughout the text. Note. a For the hcp lattice, a eff is given by the shorter lattice vector a 0 =0.250 nm compared to b 0 =0.406 nm [20].
Therefore, in the following sections we will address the question if the proposed temperature scaling (equation (4)) can successfully reproduce the Curie temperatures 3 of different materials and whether it is necessary to strictly fulfil the condition given in equation (6).

Micromagnetic simulations
To test the presented model various Ni bulk and thin films were simulated with the version 1.2b2 of the Object Oriented Micromagnetic Framework (OOMMF) [11] which implements micromagnetic simulations in magnetic solids described by the finite-difference method. The temperature effects were simulated with the thetaevolve extension = + + ) was determined in all simulations from the average over the last hundred values before reaching the end of the simulations at 50ps. As initial condition the magnetisation was maximum along the x axis, which was oriented in plane. Under the influence of the temperature the magnetisation in each simulation cell gets disturbed by the fluctuating field and the total magnetisation, averaged over all cells, decreases. As an example the time evolution of a nickel film with 6.3 ML is shown for four simulated temperatures in figure 2. Within the first picosecond's, the magnetisation decreases quickly. Afterwards the approach towards the equilibrium magnetisation slows down and the influence of fluctuations becomes visible. As expected, fluctuations are stronger for higher temperatures.

Effects of the cell size
To asses the general power of the Langevin extension to estimate the Curie temperature without any scaling, a simulation of a bulk Nickel system with 'atomic' dimensions of the simulation cell (a sim =a eff =0.354 9 nm) and no temperature scaling (T sim =T eff ) a system size of 30×30×10 cells, a time step of Δt=0.1 fs and damping parameter α=0.05 [21,22] was performed (figure 3 red squares). All other simulation parameters of the simulated bulk systems are given in table 1. This system was compared with a similar system (figure 3 black circles) with increased cell size (a sim =1nm) and time step (Δt=1 fs). After temperature scaling for the system with a sim =1 nm according to equation (4), both simulations show the same decrease of the magnetisation in dependence of the effective temperature ( figure 3).
To be able to perform simulations for systems with increased size it is often beneficial to be able to use shorter simulation time steps and increased simulation cell sizes. Therefore, further simulations were performed with increased number of cells (50×50×10), and settings optimised for higher simulation speed (Δt=10 fs, α=0.5) for four different cell sizes (figure 3: 1 nm pink stars, 2 nm green cross, 4 nm blue triangles, 8 nm black diamonds). The T eff dependent magnetisation for cell sizes of 1 nm and 2 nm systems and Δt=10 fs show good agreement with the magnetisation dependence of the simulations with shorter time steps despite the fact that a sim >λ thex . All systems approach the region of the phase transition, where M(T) drops below 5% of M(0 K), between 630 K-670 K. This results are comparable with the Curie temperature of Ni (T 627 K C Ni = ) [19]. We note here, that simulational and experimental results are in agreement despite the fact that a T sim thex C l > ( ). For a sim =4 nm the magnetisation decreases slower, with transition temperatures above 700 K. The cell length of 8 nm, which is even above the non-thermal exchange length λ ex (table 1), leads to an even stronger overestimation of the Curie temperature above 1000 K.

Bulk ferromagnets
To test weather this parameter set (1 nm cell size, 50×50×10 cells, Δt=10 fs, α=0.5) is capable to estimate the Curie temperature of other ferromagnetic materials temperature dependent simulations were performed for Nickel, Cobalt and Iron. The results are displayed together in figure 4 with Bloch's T 3/2 law and Curie temperatures from the literature [19]. After temperature scaling according to equation (4), the simulation predicts the phase transition of Nickel around T 630 K

Thin ferromagnetic films
The behaviour of ferromagnetic systems with thickness values of a couple of mono layer (ML) is of high interest due to the possibility to tune magnetic properties such as the crystalline anisotropy or Curie temperature by changing their thickness [23,24]. For example, the type and orientation of crystalline anisotropy of thin Ni films increases by a factor of twenty in thin film systems compared to the bulk, and between 6-40 ML it is best described as uniaxial, which leads to a perpendicular easy axis [24,25]. This makes them a promising material for applications in data storage or electronics [2,26]. Due to their reduced Curie temperature [23,27] inclusion of temperature effects becomes even more important than for bulk systems. To test weather the temperature scaling is applicable to low dimensional systems, or if finite size effects have to be included, the behaviour of Ni  thin films with 6 to 36 ML (1 nm to 6 nm thickness) were simulated. Within films with a thickness below the exchange length, an uniform magnetisation can be assumed [24]. Therefore in a first simulation run only one cell in z direction was simulated, the whole system was 100×100×1 cells big, the time step Δt=10 fs, and the damping parameter α=0.5 and K=1.5×10 5 Jm −3 . During a second simulation run the cell length was set constant (1 nm) in z direction, and the number of cells was varied (nz=1 -6), with 100×100×nz cells. When the simulated Curie temperatures are compared with experimental values, it has to be noted that they strongly depend on substrate and growth conditions [23,24,27]. For example, T C of a 10ML Ni film grown on a Cu(110) or Cu(111) surface can vary approximately over 50 K [27]. Therefore we compare our results with a model developed by Zhang et al [27], for the thickness dependent (n, no of mono layers) change of the Curie temperature (T C (n)) of thin films: With T C ¥ ( ) as the bulk value of the Curie temperature, N 0 as a material dependent constant, and λ as the 'temperature shift exponent' which depends on the spin-spin interaction of the material, and the condition that n>N 0 . For Nickel, Zhang et al determined N 0 =4.7 ML and λ≈1 [27]. The simulations result and the model according to equation (10) are shown in figure 5. It is evident that the simulated thickness dependence of the Curie temperature in the single cell simulation (figure 5, black squares) is best described by linear dependence (linear fit, R 2 =0.99-not shown) instead of the model given in equation (10). In contrast, the multi cell simulation (figure 5, blue circles) with homogeneous cell length follows the behaviour of equation (10) (red curve) for Nickel and is therefore a good approximation to the model, derived from experimental data.

Summary
To describe the movement of the magnetic moment of a ferromagnet at elevated temperatures, the Landau Lifhshitz equation can be extended into a stochastic equation of the Langevin type. Temperature effects are included by addition of a fluctuating field with white noise properties. The application of such a stochastic equation in micromagnetic simulations needs a discretization of the geometry of the magnet into different simulation cells. When this fundamental simulation cells include more than one spin [2], a scaling of the temperature should be performed. Based on considerations of the energy of magnetic interactions at the transition from ferromagnetic to paramagnetic phase, a scaling relation (equation (4)) in linear dependence of the simulation cell length was proposed. Under considerations of the thermal effects on the exchange length, the application of the proposed model was investigated. The scaling model was applied in micromagnetic simulations for ferromagnetic materials of general interest, Nickel, Cobalt and Iron. It was tested for different cell sizes, temperatures and bulk as well as thin film systems. Within the simulated parameter range, the best agreements between simulated and experimentally determined bulk Curie temperatures was achieved for Nickel (within 1%), followed by Cobalt (overestimated by 8%) and Iron (overestimated by about 20%). The simulation of Nickel films in the range of 6 ML to 36 ML reproduced quantitatively the change of the Curie temperature with film thickness as predicted by theoretical models [27] when simulation cell sizes of 1 nm are considered. Generally, simulation cell sizes below 2 nm gave the best results for the determination of the Curie temperature for the simulated systems even without fulfilling the condition regarding the thermal exchange length as defined in equation (7). This is in agreement wit the results by Martinez et al [15] which concluded that in this form, λ thex '...does not represent the true scale over which thermal fluctuations decay in space.' [15] Thus, the question how to obtain a general expression for the determination of the maximal length of a simulational cell size remains open [15].
In conclusion, the presented scaling model for temperature scaling in micromagnetic simulation with the stochastic Landau-Lifshitz equation is, despite its simplicity, able to predict the temperature dependent decrease of the magnetisation up to the Curie temperature. Nevertheless, care has to be taken to choose simulation cell size and simulational time steps to achieve convergence of the results.