Models of Advance Recording Systems: A Multi-timescale Micromagnetic code for granular thin film magnetic recording systems

Micromagnetic modelling provides the ability to simulate large magnetic systems accurately without the computational cost limitation imposed by atomistic modelling. Through micromagnetic modelling it is possible to simulate systems consisting of thousands of grains over a time range of nanoseconds to years, depending upon the solver used. Here we present the creation and release of an open-source multi-timescale micromagnetic code combining three key solvers: Landau-Lifshitz-Gilbert; Landau-Lifshitz-Bloch; Kinetic Monte Carlo. This code, called MARS (Models of Advanced Recording Systems), is capable of accurately simulating the magnetisation dynamics in large and structurally complex single- and multi-layered granular systems. The short timescale simulations are achieved for systems far from and close to the Curie point via the implemented Landau-Lifshitz-Gilbert and Landau-Lifshitz-Bloch solvers respectively. This enables read/write simulations for general perpendicular magnetic recording and also state of the art heat assisted magnetic recording (HAMR). The long timescale behaviour is simulated via the Kinetic Monte Carlo solver, enabling investigations into signal-to-noise ratio and data longevity. The combination of these solvers opens up the possibility of multi-timescale simulations within a single software package. For example the entire HAMR process from initial data writing and data read back to long term data storage is possible via a single simulation using MARS. The use of atomistic parameterisation for the material input of MARS enables highly accurate material descriptions which provide a bridge between atomistic simulation and real world experimentation. Thus MARS is capable of performing simulations for all aspects of recording media research and development. This ranges from material characterisation and optimisation to system design and implementation.


Introduction
Magnetic recording using hard disk drives remains the dominant technology for cloud-based information storage. As data centres consume sufficient energy to represent a significant contribution to global warming there is an imperative to improve energy efficiency and minimise the required number of data centres by means of increasing storage density. In order to achieve the properties required for magnetic information storage, the storage medium must be granular in nature. The essential required property for the storage of binary information is the presence of a large magnetic anisotropy (K): the material property which provides an energy barrier to switching of the magnetisation and thereby creates a two-state magnetic system. The increase of areal density generally proceeds by a scaling of properties, particularly a reduction of the grain size to ensure adequate signal to noise ratio (SNR) given the reduction in bit size. This necessitates an increase in the value of K to ensure thermal stability, which makes the writing of individual bits more difficult due to an increase in write field: the magnetic 'trilemma' [1]. Current technology (perpendicular magnetic recording) is already running into write-field limitations and a step change of technology is required for future products. Based on the ASRC (Advanced Storage Research Consortium) road map, there are two future technologies: 1) heat assisted magnetic recording (HAMR) and 2) Bit patterned magnetic recording (BPMR), which combined can lead to Heated-Dot magnetic recording [2]. To facilitate the development and optimisation of the present-day and future magnetic recording technologies advanced models with greater levels of complexity are required. These models must capture the dynamics at both short and long timescales over a large number of bits/grains whilst accurately describing the variation of magnetic fields, temperatures and temperature dependent parameters.
HAMR is currently the most promising new technology to provide recording densities significantly greater than those available via current standard perpendicular magnetic recording [3]. HAMR provides increased areal densities through the utilisation of high coercivity materials, potentially enabling up to 4T b/in 2 [4]. The initial difficulty with using high coercivity materials is the requirement of increased writing head field gradients. The solution to this is to temporarily reduce the coercivity of the medium, thus enabling writing with reduced field strengths. This process is achieved by applying a laser pulse to heat the material and cause a reduction in the material anisotropy. The physics of HAMR, involving heating up to or beyond the magnetic ordering (Curie) temperature T c remains challenging and involves models with a thermodynamic basis beyond the scope of those used in the typical micromagnetic approach.
Atomistic models provide this level of detail and have been used to provide temperature dependent magnetic properties and reversal mechanisms in recording media [5,6,7,8,9,10,11,12,13,14]. Although atomistic simulations can provide exceptional detail of the underlying physical processes which govern their macroscopic properties these simulations carry a significant computational cost. This cost limits atomistic simulations of recording media to a lengthscale of a few grains and a timescale of nanoseconds. This significant computational cost limits the investigation of statistical variation in particle properties and inter-particle interaction or temperature/field profiles over a track of recording bits. However, to be able to design, test and optimise any present-day and future magnetic recording technology, it is vital to capture the recording of thousands of grains from the sub-nanosecond timescale to a data storage timescale (5-10 years).
Furthermore, experimental characterisations and tests are performed at the nanosecond timescale for FMR and millisecond timescale for standard measurements such as hysteresis loops, thermal decay, First Order Reversal Curves (FORC) and thermoremanence [15,16]. Over long timescales, thermally activated transitions over the energy barriers can lead to loss of recorded information. Transition times are governed by the Arrhenius-Néel law [17], which gives a characteristic time τ −1 = f 0 exp(KV/kT ), where V is the grain volume and K is the magnetic anisotropy constant. The pre-exponential factor is dependent on the value of K and is in the region of GHz to THz. Ensuring thermal stability of the written information for 5-10 years requires large energy barriers (KV/kT > 80) and good SNR. Typically short and long timescale investigations are performed separately, mainly due to the described fundamental difference between the simulation methodology. Nonetheless, there are numerous cases where both timescales are of interest, the simplest example of this is the effect of the writing process on nearby bits (nanosecond timescale) and data longevity ( timescale of years).
Here we present the developed multi-timescale micromagnetic code, MARS (Models of Advanced Recording Systems), which has the functionality of utilising various solvers to best accommodate the required simulation time frames. MARS includes both short and long timescale solvers to make such investigations more simple and to open up the possibility to more easily access the effect of dynamic processes on the long timescale behaviour. MARS includes a stochastic-Landau-Lifshitz-Gilbert (sLLG) equation solver along with a stochastic-Landau-Lifshitz-Bloch solver, specifically the sLLB-II (which is shown to have greater accuracy at the Curie temperature [18]), to simulate short timescale dynamics of the magnetisation. A kinetic Monte-Carlo (kMC) solver is also included in order to simulate the long timescale behaviour. Fig. 1 illustrates the benefits of the multi-timescale micromagnetic approach over computationally expensive atomistic simulations. The addition of the kMC solver enables the simulation of timescales far beyond the standards set using dynamical micromagnetic solvers alone.
MARS has been designed to be used alongside atomistic simulations by utilising parameterisation obtained via atomistic simulation to describe material properties, this allows for highly accurate descriptions of simulated materials. There is no limit on the number of materials MARS can utilise in a single sim- Figure 2: Example of the granular structure obtained via the centroidal Voronoi tessellation (a). A subsection has been chosen to illustrate Lloyd's relaxation algorithm after: zero (b), one (c), two (d) and three (e) iterations. The seed points are indicated by the crosses with the centroids represented by the circles. As more iterations are performed the angular nature of the grains is reduced.
ulation. This enables multi-layered systems such as those used for exchange-coupled-composites to be simulated. Details of granular system generation and the numerical solvers are provided in section 2. MARS includes a set of commonly used simulation types for easy use, these include: HAMR writing, time evolution and data read back, thermoremanence, FORC and FMR. A selection of results obtained via these simulations are provided and discussed in section 3.

Granular model
The micromagnetic approach consists of treating magnetic grains as macrospins with associated magnetisation, m. To ensure accurate modelling of these systems it is crucial to generate accurate granular structures. The typical method for generating granular structures is via Voronoi tessellation [19,20,21,22]. The classical Voronoi algorithm starts with the creation of seed points throughout the system, cell walls are then created such that they lie halfway between two seed points. While the general process is the same there exist various methods to determine the initial locations of the seed points [23,24].
A major drawback with the classical Voronoi construction is evident when a distribution of grain sizes is required. When there is a local increase in seed density the construction can generate unrealistically angular cells. This occurs due to the only constraint on cell construction being the requirement that the cell wall must be equidistant between seed points.
To overcome this drawback, one can utilise centroidal Voronoi tessellation, this modified Voronoi process has been shown to be most effective when combined with Lloyd's algorithm [25,26]. The process involves iterative relaxation of the constructed granular system by replacing the initial seed points with the centroids (typically known as the centre of mass) of the generated cells until convergence is achieved. Fig. 2 shows the system generated via MARS using centroidal Voronoi tessellation followed by the changes produced by applying Lloyd's algorithm over three iterations.
A key limitation to all seed-based construction techniques is that they produce Gaussian cell size distributions [27,28]. In re-ality the grain size distribution has been shown to be described best by lognormal or Gamma distributions [29,4,30]. To enable the construction of systems following these distributions a method for performing a Voronoi construction using 'hard discs' instead of seeds is available, called the Laguerre-Voronoi method. The main challenge to the Laguerre-Voronoi method is the initial packing of the randomly sized hard discs. There are numerous methods for packing these hard discs, for MARS a custom "Drop and Roll" method has been implemented. This method provides a high level of contact between neighbouring discs resulting in a greater packing fraction while being more computationally efficient [31]. Using this algorithm MARS is capable of generating packing fractions of at least 80% for randomly distributed disc sizes. Once the system has been packed the cells are generated using the robust open source VORO++ package developed by Rycroft [32]. Fig. 3 shows the difference between the structures generated via the centroidal Voronoi tessellation and the Laguerre-Voronoi tessellation, the corresponding grain size distributions are shown in Fig. 4. The implementation of the Laguerre-Voronoi tessellation method within MARS enables the generation of realistic systems with a high level of control over the grain size distributions. For micromagnetic simulations periodic boundary conditions are used for the Voronoi tessellation in order to remove edge effects, Fig. 5 shows a system created via a Laguerre-Voronoi tessellation with the periodic boundaries indicated by the dashed lines.

System energy and effective fields
For a magnetic film composed of N individual grains, which we can consider as macrospins, the energy of the system can be  and Laguerre-Voronoi tessellation (bottom). The input distribution was lognormal with µ = 1.65 and σ = 0.55. The seed based Voronoi is unable to provide the desired distribution, instead providing a Gaussian. The Laguerre-Voronoi was able to produce the desired distribution type with only a small change of parameters, which is to be expected due to the random nature of the packing process.
written as: where V i is the volume of grain i with uniaxial anisotropy energy density K i , easy axis directionê i , bulk saturation magnetisation M i s and normalised magnetisation vector m i = M i /M i s . The anisotropy energy density, saturation magnetisation and fractional exchange constant J i j are all temperature dependent. The first term is the Zeeman energy which describes the interaction of the grains with an external applied field H app . The second term is the anisotropy contribution to the energy which describes the preferred alignment direction for the grains magnetisation. The third term describes the exchange interaction between nearest neighbour grains which can be expressed in terms of the local exchange field H i exch The fourth term is the long-range magnetostatic interaction, within the dipole approximation, which couples grains at sites r i and r j at a distance r i j = r j − r i across the whole system. The effective field acting on each grain is obtained from the energy of the system (Eq. 1) and is given by: The individual terms are described in the following. H i exch describes the coupling between different grains, belonging either to the same layer or to different layers as in the case of exchange-coupled composite (ECC) media. In the case x y Figure 5: Example of the periodic system generated via MARS using a Laguerre-Voronoi tessellation. The grains containing points are those generated, with the dashed lines indicating the periodic repetitions. of a granular medium the exchange results from the intergranular medium. Although this is engineered to ensure exchange decoupling, this is not necessarily complete: in fact in the case of media for perpendicular recording the exchange, which balances the effects of the magnetostatic field, is a part of the material design. Under the reasonable assumption that the intergranular exchange is proportional to the contact area between the grains, Peng et. al. have shown that the exchange H i exch is given by [33]: J i j is the fractional exchange constant between the adjacent grains and L i j is the contact length between grains i and j, A i is the area of grain i, denotes the average value and H sat is the exchange field strength at saturation, which is generally derived from experiment. Sokalski et.al. [34] investigated experimentally the exchange coupling between thin layers of CoCrPt separated by an oxide, finding an exchange strength which decayed experimentally with oxide layer thickness. Ellis et.al. [11] found a similar relation using an atomistic model based on the presence of ferromagnetic impurities in the oxide layer. This study also showed the presence of higher order (biquadratic) exchange and importantly demonstrated that the intergranular exchange decayed to zero rapidly with increasing temperature, suggesting that intergranular exchange does not play a major role in the HAMR recording process. It is important to note, given the likely origin of intergranular exchange in the presence of impurity magnetic spins in the intergranular layer, that J i j could vary significantly. According to Peng et.al. [33] this can lead to exchange weak links which act as pinning sites and reduce the sizes of clusters arising from magnetostatic interactions.
H i dmg is calculated using the dipole approximation: where W i j is the demagnetisation tensor of the system: V is the volume of the grain, r i jα is the displacement between grains i and j, with the subscript α = x, y, z denoting the component of the displacement. As W i j is dependent only on the position and size of the grains this matrix can be determined prior to simulation internally or via a separate external code. Improved methods to determine the W matrix are available but these produce additional computational cost. One such method is surface charge integration as discussed in [35]. MARS is capable of accepting the W matrix as an input enabling fast implementation of alternative methods for magnetostatic determination. The temperature dependence of H i ani is described using the following expression: whereê i is the unit vector aligned along the easy axis, K i is the anisotropy and M i s is the zero temperature saturation magnetisation of grain i. Here we exploit the fact that we can express the temperature dependence of K via the dependence on the magnetisation m described via Callen-Callen scaling [36], which allows K(T ) to be expressed as: K 0 is the anisotropy energy density at 0 K and η is determined via experiment or atomistic parameterisation. Typically the exponent η = 3 for uniaxial anisotropy and η = 2 for 2-site anisotropy appropriate for FePt [37].

Dynamical and kinetic Monte-Carlo Solvers
MARS utilises three separate solvers to describe the magnetisation, the Landau-Lifshitz-Gilbert, Landau-Lifshitz-Bloch and kinetic Monte Carlo. The LLG and LLB solvers require very short timesteps, in the region of picoseconds and femtoseconds respectively, to function and provide dynamic information about the magnetisation. The kMC is a probabilistic solver which sacrifices the dynamic information in order to enable much larger timesteps. This enables the simulation of long timescale phenomena, for example the long-term decay of written information arising from thermal activation. Each of the solvers presented in this work have been rigorously tested to ensure correct implementation and high accuracy. The details of these tests are provided in Supplementary Notes 1-3, with Supplementary Figures 1-7 showing comparisons between produced and expected test results.
The LLG equation of motion for each grain i, including stochastic effects, is given by: The first term describes the quantum mechanical precessional motion of the magnetisation around the effective field H i eff , while the second represents the phenomenological relaxation of the magnetisation towards H i eff [38]. The Gilbert damping α couples the spin system with the environment, considered to act as the thermal bath, and determines how fast the system relaxes towards equilibrium. H i th is the thermal field, this stochastic field accounts for the thermally driven behaviour of the macrospin and is described by a non-correlated white noise Gaussian function.
where: i, j label the magnetisation on the respective sites; α, β = x, y, z ; k B = 1.381 · 10 −16 ergK −1 is the Boltzmann constant; T is the temperature; δ µγ is the Kronecker delta and δ(t − t ) is the delta function. In this formulation the noise is considered to be spatially and temporally uncorrelated, i.e., white noise. This approach works at relatively low temperatures where one can consider the grain to be fully magnetically saturated and to exhibit coherent rotation with all atomic spins remaining parallel. Under these circumstances the equation of motion need only model transverse dynamic processes and the LLG equation is valid. However, as temperature increases and approaches the Curie point this is no longer true and the Landau-Lifshitz-Bloch (LLB) equation must be used instead. The LLB introduces a longitudinal relaxation of the macrospin which accounts for the loss of magnetisation and divergence of the longitudinal susceptibility as the temperature approaches the Curie point. The LLB equation was first derived by Garanin [39,40]. The LLB equation is a single (macrospin) representation of the dynamical behaviour of a single grain and differs from the LLG equation in its inclusion of longitudinal relaxation of the magnetisation. Although the LLB equation parameters to be outlined below were originally derived from mean-field theory, these can be obtained from atomistic calculations. As shown by Chubykalo et. al. [41], the LLB equation gives excellent agreement with atomistic model calculations, essentially validating its use in calculations of HAMR, which involve heating beyond the Curie temperature. The implementation of the stochastic LLB solver used by MARS follows the work of Evans et al. [18] (sLLB-II) and for each grain i reads: where m i is the reduced magnetisation M i /M(T = 0), m i is the length of m i and γ e is the electron gyromagnetic ratio. The first and third terms are the precessional and damping terms for the transverse component of the magnetisation, as in Eq. 8, while the second and fourth terms are introduced to account for the longitudinal relaxation of the magnetisation with temperature. The stochastic LLG and LLB solvers both utilise the Heun integration scheme. The benefits of the Heun scheme are two-fold.
First it provides second order accuracy in ∆t for the deterministic part, thus rendering it more numerically stable than Euler type schemes. Second, it yields the required Stratonovich solution of stochastic differential equations. The damping of the magnetic moment is split into longitudinal α and transverse α ⊥ components given by: Where λ is the thermal bath coupling, a temperature independent phenomenological parameter, that is the same as that used in atomistic spin dynamics. The transverse damping is related to the Gilbert damping by the expression: ζ ⊥ and ζ ad are the diffusion coefficients that account for the thermal fluctuations. The thermal noise terms are described by Gaussian functions with zero average and a variance proportional to the strength of the fluctuations: As temperatures approach and exceed the Curie point, Eq. 6 produces a fictitious longitudinal component of the anisotropy. This leads to a reduction in the longitudinal relaxation of the magnetisation as a function of temperature. To overcome this issue the anisotropy field can also be described as a function of the transverse susceptibility χ ⊥ [40]: where m i x and m i y are the components of the reduced magnetisation vector andx,ŷ are the unit vector along these directions, respectively. Unlike Eq. 6 this form of the anisotropy assumes that the easy axis lies along the z-axis however it is valid for all temperature ranges and is therefore the most suitable description for LLB applications. For soft materials the determination of χ ⊥ is extremely challenging and thus both forms of the anisotropy are available for use with the LLB solver to enable the simulation of both hard and soft materials.
The LLB equation includes an additional field term, H i intragrain , within the effective field. This term accounts for the exchange between the atoms within grain i, controls the length of the magnetisation and is given by Here m i is length of the reduced magnetisation m i of grain i, and m e (T ) is the equilibrium magnetisation. The term H intragrain encapsulates the new physics introduced by the LLB equation. It incorporates the longitudinal fluctuations of the magnetisation while maintaining a mean value m e (T ). It is important to note that the fluctuations diverge asχ diverges close to T c . This is responsible for the onset of the linear reversal model close to T c . Clearly specification of the temperature dependence of the LLB parameters is of paramount importance: a factor complicated by the effects of finite size on the magnetic properties.

Atomistic parameterisation
The granular model requires characterisation of the temperature dependence of the magnetisation, anisotropy and susceptibilities. These quantities are obtained via fitting of atomistic data, obtained using the VAMPIRE package [42]. A key benefit of atomistic parameterisation is the improved accuracy of the modelled material's behaviours as well as the ability to simulate granular systems which include a segregant between the grains as is typically the case in recording media. There are two available methods for fitting the magnetisation. The first is fit- Both methods are capable of producing the characteristic behaviour of the temperature dependent magnetisation. A comparison of these two methods is given in Fig. 6. For bulk systems a strong criticality is expected and the critical exponent fit reproduces the sharp transition to zero magnetisation at the Curie point. However, as grain sizes decrease finite size effects become significant which cause a reduction in the criticality of the transition. The result of finite size effects is a small but nonzero magnetisation above the Curie point. The polynomial fit is capable of reproducing this behaviour and provides greater agreement with atomistic data for small grains (i.e. 5 nm) than the critical exponent fit.
The susceptibility χ is a measure of the strength of the fluctuations of the magnetisation. The components of the susceptibility, according to the spin fluctuation model, can be obtained by the fluctuations of the same magnetisation components as follows [43]:χ Whereχ α = χ α /M s V is the reduced susceptibility and is in units of field −1 . N is the number of spins in the system with magnetic moment µ s . Here m α is the ensemble average of the reduced magnetisation component α = x, y, z and longitudinal. x, y, z are the spatial Cartesian components of the magnetisation, while longitudinal describes the length of the magnetisation.χ describes the strength of the fluctuations of the magnetisation component along the easy-axis direction, which for our system is z.χ ⊥ refers to the fluctuations orthogonal to the easy axis and thus on the x-y plane. Forχ andχ ⊥ , we use a similar approach to Ellis [43] and we fit the inverse of the susceptibility 1/χ ,⊥ : 1 Where C i and D i are the fitting parameters and T c is the Curie point, obtained by determining the temperature at which the susceptibilities intersect. Fig. 7 shows the susceptibilities and fits obtained from atomistically parameterised FePt, the Curie point of this system is 685.14 K.
Low anisotropy systems and systems of reduced dimensions cannot retain the alignment of the magnetisation along the easy axis up to T c . In such casesχ is a mix of the spatial components and becomes difficult to determine. A workaround is to avoid calculatingχ directly and to obtainχ from the longitudinal susceptibilityχ l , following the discussion presented in [44]. Unfortunately a similar method cannot be used forχ ⊥ making it difficult to determine the anisotropy for soft systems when the anisotropy field is given by Eq. 14.
Alternatively, if the anisotropy field is described as in Eq. 6, the reduced anisotropy is given by k(T ) = K(T )/K 0 = m(T ) γ , as discussed by Callen-Callen [36]. MARS implements both a standard Callen-Callen fitting and an extended version. The extended version utilises three temperature regions each with their own fit parameters such that there are no discontinuities. This extended fitting method enables greater accuracy in the reproduction of the anisotropy as a function of temperature. This  Figure 7: Fit obtained for parallel and perpendicular susceptibility using the inverse method similar to that of Ellis [43]. The Curie point of this system is the temperature where the susceptibilities first intersect, which for this data is 685.14 K approach should provide more useful results in the case of soft materials, where extractingχ ⊥ can prove difficult. Fig. 8 is a comparison of the fits obtained using the standard and extended Callen-Callen fitting methods. Once all these parameters are determined, the granular model is fully parameterised regarding the material properties.

Kinetic Monte-Carlo solver
Finally we turn to the solver for long-timescale simulations: the kinetic Monte Carlo (kMC) solver. In the kMC approach as given in [45], the switching probability is dependent on the measuring time t m as follows, where the relaxation time τ is given by the Arrhenius-Néel law [17] τ where f 0 is the attempt frequency, usually assumed around 10 9 s −1 for these systems, and ∆E 12,21 are the energy barriers for switching between states. Typically, large energy barriers of > 60k B T are required in order to ensure long-term thermal stability of written bits. τ −1 is given by τ −1 = τ −1 12 + τ −1 21 . To model the physical effect of the easy axis dispersion, easy axes are chosen randomly within a Gaussian dispersion of angle about the normal. The total energy barrier including the effect of anisotropy dispersion can be written as where H T the total effective field, g(ψ) = [cos 2/3 ψ + sin 2/3 ψ] −3/2 and κ(ψ) = 0.86 + 1.14g(ψ) are the numerical approximations given by Pfeiffer [46]. The kMC is capable of simulations on the timescale of years and is valid for systems where the energy barrier is much larger than the thermal energy. In order to function, the kMC requires calculation of the magnetisation states corresponding to the energy minima along with the energy barrier separating the two states. Stationary states are found by solution of the Stoner-Wohlfarth model of coherent rotation [47] for which the free energy is given by The final step [45] is to ensure that, after switching, the populations of the energy minima obey Boltzmann statistics. This approach leads to the condition that, if the reversal transition is allowed, the moment is then assigned to either energy minimum with a probability, with i = 1, 2 labelling the minima, thereby ensuring that the population of the two states obeys the Boltzmann distribution in thermal equilibrium [45]. This is important in order to include the 'backswitching' mechanism which leads to dc or 'remanence' noise. Taking account of both distributions the transition probability is determined using: Where P 2 is the probability of the magnetic moment jumping to the second minima and ∆E is the energy barrier separating the two minima. To determine if a switching event occurs a random number between zero and unity is generated and compared to P 2 . If it is less than P 2 the magnetic moment orientation is assigned corresponding to the second minimum otherwise it is assigned to the first minimum. t m is the measurement time. During the measurement time the external properties such as magnetic field and temperature are assumed constant, such that Eq. 24 can be applied.

Curie temperature dispersion
The HAMR process involves heating through T c which, as a result, becomes an important material parameter. More particularly, simulations by Li and Jhu [48,49] have shown that the dispersion of T c is a serious limitation for the ultimate storage density achievable for HAMR. Here we consider an irreducible contribution to the dispersion of T c which arises directly from the diameter dispersion. It is well known that finite size effects lead to a reduction of T c , demonstrated experimentally for FePt by Rong et.al. [50]. A theoretical investigation based on an atomistic model by Hovorka et. al. [51] showed that the variation M(T ) was well described by the finite size scaling law where D is the grain diameter, λ is the so-called phenomenological shift exponent and d 0 is the microscopic length scale close to the dimension of a unit cell of the lattice structure. The exponent λ is related to the correlation length universal critical exponent ν and it is expected that λ = ν −1 . However, small grains can exhibit departure from universality so we prefer the form of Eq. 25 as a functional form to represent the diameter dependence of T c . Clearly a dispersion of diameter maps onto the dispersion of T c . Assuming a lognormal distribution of D, with logarithmic mean D m and variance σ 2 D it has been shown [51] that the dispersion of T c is given by the distribution function with ∆T c = T ∞ c − T c . Eq.(26) is a lognormal distribution function with logarithmic mean T m = λ(ln(d 0 (T ∞ c ) 1/λ ) − D m ) and variance σ 2 T = λ 2 σ 2 D . Through Eq. 25, with d 0 , λ and T ∞ c determined either from experiment or atomistic model calculations, a T c value can be assigned to an individual grain and Eq. 26 used to calculate the standard deviation of the T c dispersion.

Simulations
This section presents a demonstration of some of the bundled simulation types available in MARS at release. These simulations serve to show the capability of MARS and provide some of the most common simulation types. The first example simulation is the writing of a pseudorandom bit sequence (PRBS) pattern for HAMR. The second example is the determination of switching probabilities for characterising Curie point dispersion via comparison with experimental data. Finally an example of FMR simulations for use in determining system damping is provided.

Writing and reading processes for heat assisted magnetic
recording Development of improved heat assisted magnetic recording requires investigations into the thermal reversal of the grains along with the effect of bit spacing on data stability and writing performance. MARS contains three separate HAMR focused simulations. The first applies a temporally dependent laser and field profile onto a single bit, allowing investigations of thermal reversal properties of materials. The second simulation, consists of a bit within a surrounding system of grains, utilising a temporal and spatially dependent laser and field profile, allowing for investigation into the influence of the written bit on any surrounding grains. The third simulation models realistic data writing, via a square-wave or user specified binary sequence. This writing can occur for single or multiple tracks. Using this third simulation, comprehensive investigations into the entire HAMR writing process can be achieved.
Read back of a system can also be simulated in MARS. The system is first discretised into 1 nm pixels, this enables more precise measurement of the magnetisation and produces a more realistic read back signal. A read head is then placed at the edge of the track. This head is scanned along the track and the magnetisation detected within the head is averaged and recorded as a function of down track position. In order to investigate data decay over time the simulation utilises the kMC solver to simulate the system for years and performs read back at specified intervals. Fig. 9 shows the output obtained via the realistic HAMR writing and reading simulation available in MARS. For this simulation two systems with energy barriers of KV/k B T = 79 and 62 were used. Each system consisted of 1,300 grains with a lognormal grain size distribution with an average diameter of 8 nm. The material parameters used were M s = 1,051.65 emu/cm 3 , T c = 693.5 K, λ = 0.1 and K u = 9.2 · 10 7 erg/cm 3 and K u = 7.1 · 10 7 erg/cm 3 respectively. A 31-bit PRBS was written to the systems. The sequence used was 1111100011011101010000100101100. Both systems were then evolved over time for ten years and read back to generate a second read back signal. This simulation was performed 100 times with different random seeds, the read back signal was then obtained for each simulation in order to obtain the spatial, i.e. of the medium, signal-to-noise ratio (SNR) as a function of time. There are two contributions to the spatial SNR: transition and remanence [22]. These SNRs describe the noise present at and away from the bit transitions respectively. The latter is caused by grains with KV/k B T such that switch back, whereas the former is a measure of how good a transition between bits is written. Transition noise is the noise expected to be dominant in high density media due to the reduce dimensions the bit size. The calculation of SNR is currently not included in the MARS software package, even though there is the plan to include it in future releases. The extraction of SNR is based on the ensemble wave form analysis developed by Seagate [52,53], a method that allows to separate and extract the different noise components: transition and remanence. The approach proceeds as follows: a track is written multiple times, 10 times in the present case, with each track read back once only since reader noise is not included. The signals are first synchronised via cross-correlation and afterwards the average "noise-free" signal is obtained by averaging over the signal of the 10 written tracks. Then, the total spatial noise is calculated as the variance of the average "noise-free" signal. To extract the remanence and transition noise appropriate windowing functions are applied to the total spatial noise. Eventually, the SNRs are calculated as the 10 log 10 of the ratio between the total signal power and the respective noise power. Fig. 10 shows the SNR components for the high energy barrier system (KV/k B T = 79), and as expected the transition noise is the bottleneck.

HAMR switching probability
The heating in HAMR systems occurs in a narrow temperature region near T c , thus the performance of HAMR is highly dependent on the T c dispersion of the recording medium [19,48,49]. To determine the T c dispersion experimentally the procedure requires measurement of the thermo-remanence magnetisation as a function of the applied laser pulse peak power [16]. The application of this laser pulse occurs over the long timescale thus LLB simulations are extremely computationally expensive, to the point of impracticability. MARS overcomes this problem by utilising the multi-timescale nature of the code. The solver used by the code is determined automatically based on the laser application time: if this time is of the order of microseconds or greater then the kMC is used otherwise the LLB is used. The heating and cooling phases require dynamic information to be modelled and hence are always performed via the LLB solver. The ability for MARS to automatically select the most suitable solver allows for the simulations to be run in batch jobs without specification of the solver or editing of the source code. Fig. 11(c,d) are the results of a thermoremanence simulations used to investigate the switching probability of single FePt grains as a function of peak temperature for HAMR-like heat pulses for various pulse lengths and under different applied fields. The considered pulse length falls within the few nanoseconds regime, thus the LLB solver has been employed in these simulations. These results are compared with the switching probabilities obtained by performing atomistic simulations for the same system and setup, presented in panels (a,b). The agreement between the atomistic parameterised LLB and atomistic simulations, which have also been used for the parameterisation, proves the ability of MARS in describing such processes, as discussed in more detail in Ref. [13]. Fig. 12 presents an example of the average magnetisation dynamics for 100 FePt grains under the application of an external −1 T field along the z-axis, and a 1.8 ns heat pulse with peak temperature T peak = 600 , 700 K. The dashed brown line shows the Gaussian profile used to model the time dependence of the heat pulse: where t pulse corresponds to 1/6 of the total pulse length and T a is the ambient temperature, the temperature at which the grains are in absence of the heat pulse, 300 K in this case. For low T peak and relatively fast pulses the magnetisation reversal of the grains cannot be achieved, however, despite this a partial reversal can be observed when the temperature approaches T peak . However, as the grains cool down the magnetisation is restored along the initial direction. On the other hand, when T peak approaches T c of the grains, all grains are reversed.

Ferromagnetic resonance
Ferromagnetic resonance (FMR) is an important technique used for the measurement of magnetic properties for bulk and thin-film magnetic media [54]. FMR simulations enable investigation into the dynamic properties of a material such as the damping [55] as well as the static properties such as saturation magnetisation and uniaxial anisotropy [56]. MARS provides a simple system to perform frequency and field swept FMR simulations over a range of temperatures via its LLG and LLB solvers. Fig. 13 shows the power spectrum obtained via FFT vs. the applied field frequency obtained via FMR simulations of a single system at various temperatures. The data have been fitted to a Lorentzian function: where A is the amplitude, w is the full width at half maximum, α is the damping parameter and f 0 is the resonant frequency. The resonance frequency can be determined via the Kittel formula [57].
where B is the in-plane FMR field amplitude. From this fitting the system's damping parameter and resonant frequency are extracted. Fig. 14 shows the damping and resonant frequency extracted as a function of system temperature, there is very good agreement between the extracted results and the analytical values.

Conclusion
We have developed an open-source multi-timescale micromagnetic code (MARS) for simulating granular thin films. The primary focus of MARS is to enable detailed modelling and  Short timescale simulations are possible via the LLG and LLB dynamic solvers. The LLB enables the simulation of systems up to and exceeding their Curie points, which is crucial for recording system such as HAMR. The long timescale simulations are performed via the kMC stochastic solver. The combination of these solvers allows for complex multi-timescale simulations to be developed and performed. MARS has also been developed for use in material characterisation to aid research into the development and optimisation of recording media materials. The implementation of atomistic parameterisation enables highly accurate material descriptions within MARS and produces very good agreement between re-sults obtained using MARS and those obtained atomistically. Thus MARS is highly useful for bridging the gap between atomistic simulation and real world experimentation. Furthermore to ensure an accurate description of granular media the numerous methods of Voronoi construction have been investigated and compared resulting in the Laguerre-Voronoi method being implemented in order to ensure realistic grain shapes as well as lognormal grain size distributions. Detailed descriptions of the models incorporated in MARS have been provided along with the various methods used to model temperature dependent material parameters.
Finally we have provided example results obtained via MARS for numerous published or on-going studies to show the versatility and capabilities of MARS. These studies cover the entire range of HAMR development starting from materials characterisation and optimisation through to HAMR writing and finally data storage and read back.

Acknowledgements
We are grateful to Dr. Mara Strungaru for her aid in the development and testing of the FMR simulation feature within MARS.
The authors gratefully acknowledge the funding from Transforming Systems through Partnership Programme of the Royal Academy of Engineering under Grant No. TSP1285 and Seagate Technology (Thailand).
J.C and P.C. gratefully acknowledge the financial support from Thailand Science Research and Innovation (TSRI) Numerous simulations for this work were undertaken on the Viking Cluster, which is a high performance compute facility provided by the University of York. We are grateful for computational support from the University of York High Performance Computing service, Viking and the Research Computing team.

.1 Verification of the Heun scheme
The dynamic solvers utilise the Heun integration scheme due to its ability to provide the required Stratonovich solution of stochastic equations. However, first this method must be shown to accurately solve the deterministic equations of motion. For this test analytical solutions were derived for specific cases for both the LLG and LLB solvers. For this test a single grain with an initial magnetisation along the x-axis under the influence of an applied field directed along the z-axis [1] is used, the time evolution for the LLG solver is then given by: For the LLB solver the time evolution is: The time evolution for the LLG and LLB along with the obtained errors are shown in supplementary figure 1. The maximum error obtained for the LLG is around 10 −2 , while for the LLB it is around 10 −5 . This is to be expected as the timestep utilised in the LLB is 1 fs and 1 ps in the LLG. The form of these errors is characteristic of a correctly implemented Heun integration scheme.

Angular dependence of the coercivity
Verification of the implementation of the uniaxial anisotropy for all solvers is performed via a simple test which makes use of the Stoner-Wohlfarth model. Here, the Stoner-Wohlfarth particle describes the behaviour of a single grain at zero Kelvin under the influence of an applied field. The angular dependence of the coercivity is very well known [2] and thus makes a very useful property for comparison and verification of a code. This test verifies the deterministic behaviour of the LLG and LLB by ensuring the easy axis profile provides a coercivity of H k = 2K Ms as is known analytically. The grain is initially magnetised along the z-axis and the applied field strength is varied from H = 1.5H k to H = −1.5H k and back, in steps of 0.005H k . This process is repeated for a range of field directions ranging from 90 • to 0 • . The projection of the magnetisation on the field direction is then plotted as a function of the applied field strength as shown in supplementary figure 2 for the LLG, LLB and kMC solvers. As expected the same profiles are obtained irrespective of which of the three solvers are used and these profiles also agree exactly with those of Stoner and Wohlfarth's solution. upplementary Figure 2: Alignment magnetisation for a single grain under an applied field for various angles from the easy axis. The 0 • and 90 • profiles are simulated with a small deviation from the labelled value. This is done as at perfect alignment, in the absence of thermal noise, there is zero torque which prevents any change in the magnetisation alignment.

Boltzmann distributions for an ensemble of non-interacting grains
To verify the stochastic behaviour of the sLLG and sLLB-II solvers a test which utilises thermal effects is required. The most simple test available is the reproduction of the Boltzmann distribution for an ensemble of non-interacting grains, with uniaxial anisotropy in the absence of an applied field. The Boltzmann distribution for magnetisation as a function of polar angle for the LLG is simple to obtain and is given by: where θ M is the angle of magnetisation with respect to the easy axis. The ensemble is allowed to evolve over time for 10 6 steps and the angle of magnetisation is recorded at each step. The results obtained by the sLLG and the analytical solution are shown in supplementary figure 3 for an initial magnetisation direction parallel to the easy axis direction. There is excellent agreement between these results, showing correct implementation of thermal effects in the sLLG. For the LLB equation the free energy of the system is defined as [3,4]: Where the first term provides uniaxial anisotropy and the second controls the magnetisation length. The expected Boltzmann distribution is of the form: Probability distributions along m z and m x for different temperatures along with the corresponding analytical solutions are plotted in supplementary figure 4, all results show strong agreement with the analytical solutions. The simulated system has a Curie point of 685.14 K, includes anisotropy with an initial random uniformly distributed magnetisation and has been evolved over time for 1 ns (10 6 steps). The results show the ability of the sLLB-II solver to describe the loss of ferromagnetic behaviour as the temperature approaches and exceeds the Curie point, indicated by the formation of a single peak once the temperature exceeds the Curie point.  The LLG conserves magnetisation length assuming unity for all temperatures, thus the model starts to break down at temperatures approaching T c . The LLB allows for variable magnetisation lengths enabling it to simulate systems at and beyond their Curie points. This test was devised to verify the treatment of the magnetisation length by the LLB. The temperature dependence of the magnetisation is determined via parameterisation obtained using atomistic simulations. Two systems are simulated, one with a hard magnetic material the other a soft material. The hard material is modelled via the perpendicular susceptibility while the soft material utilises the Callen-Callen method. The systems are heated from 0 K to above the highest Curie point of the materials. At each temperature the system is simulated until an equilibrium is reached at which point the magnetisation is recorded and the temperature increased. supplementary figure 5 shows the results obtained for the hard and soft materials via LLB simulation compared to those obtained by atomistic simulation. There is excellent agreement between the results, showing that MARS is capable of reproducing the correct behaviour of the magnetisation length.

Longitudinal relaxation
A system of non-interacting identical grains is initially aligned at 30 • from the z-axis, placing it in a nonequilibrium state. The system is then allowed to thermally relax, in the absence of an applied field, for a range of system temperatures and the magnetisation length is recorded at each step. supplementary figure 6 shows the comparison of the results obtained via MARS and atomistic simulations. The results show excellent agreement. As expected there is a small discrepancy for very short timescales where the rate of relaxation is greatest [3], however, both results show identical equilibrium lengths within only 5 ps.

Coercivity as a function of sweep rate
In the presence of thermal effects the coercivity depends on the applied field sweep rate. The coercivity as a function of the sweep rate was determined empirically by Sharrock [5] and then derived theoretically by Chantrell [6] under the assumption of constant attempt frequency and an easy axis parallel to the applied field. A simple test for the kMC solver consists of the simulation of hysteresis profiles for a range of sweep rates. Supplementary figure 7 shows the comparison between the simulated results and the theoretical prediction.
The results agree strongly with the theory verifying the implementation of the kMC for thermal systems.