First order reversal curves and intrinsic parameter determination for magnetic materials; limitations of hysteron-based approaches in correlated systems

The generic problem of extracting information on intrinsic particle properties from the whole class of interacting magnetic fine particle systems is a long standing and difficult inverse problem. As an example, the Switching Field Distribution (SFD) is an important quantity in the characterization of magnetic systems, and its determination in many technological applications, such as recording media, is especially challenging. Techniques such as the first order reversal curve (FORC) methods, were developed to extract the SFD from macroscopic measurements. However, all methods rely on separating the contributions to the measurements of the intrinsic SFD and the extrinsic effects of magnetostatic and exchange interactions. We investigate the underlying physics of the FORC method by applying it to the output predictions of a kinetic Monte-Carlo model with known input parameters. We show that the FORC method is valid only in cases of weak spatial correlation of the magnetisation and suggest a more general approach.

Identification of the intrinsic properties of magnetic nanostructures is central to the development of applications in a wide range of topics: information storage, biomedicine, permanent magnet development, and many more. The parameter identification techniques are at the heart of large scale material characterisation to quantify the properties of nanoscopic constituents of materials. For example, the optimisation of magnetic granular materials for the current and future hard disk drive technologies, such as heat assisted magnetic recording (HAMR), or the synthesis of magnetic nanoparticles for molecular sensing and detection, imaging, and cancer therapy in biomedicine relies on the possibility of efficient and accurate identification of the physical properties of billions of magnetic nanoparticles, which requires analysis in a high dimensional parameter space and the employment of a statistical approach.
In such cases, direct measurements targeting individual particles become inefficient and infeasible. Instead, an indirect approach based on relating theoretical models to macroscopic experimental data and identifying the model parameters from the optimal fit becomes the most viable approach. This inverse problem solving methodology relies on the availability of a realistic model capable of i) reliably representing the physics of elementary constituents of a physical system, ii) accurately reproducing the macroscopic measurement data (forward problem) and iii) understanding the uniqueness properties of inverse solutions of the model, i. e. whether the identified parameter set is the only set allowing the model to accurately reproduce the measurement data [1][2][3][4] . Unfortunately, inverse problems are often ill-posed and entire manifolds of parameters allow the models to reproduce the measurement data, which effectively translates to a significant error of the parameter identification. Such errors can only be reduced by providing new information from independent measurements. As a result, the development of identification techniques for materials characterisation remains a challenging problem. Figure 1. The ideal single particle hysteresis loop has a rectangular shape corresponding to the hysteron model (a), where the only change in magnetisation is due to switching events. The hysteresis loop will have a more complex shape (b) in the presence of thermal effects and reversal components as included in our model. Example of hysteresis loops for a HDD media system (c) and the corresponding first order reversal curves in H a H b plane (50 FORC curves are used with a field step of 100 Oe, although 5 curves only are shown for clarity). (d) or H c H u plane (e) for a non-interacting system of 10000 elongated grains (1.17 aspect ratio and D = 8.5 nm) simulated at 300 K and field rate of 4.10 4 Oe/s. The system parameters are: M s = 700 emu/cm 3 , K = 7.10 6 erg/cm 3 with 3 degree dispersion of the easy axis.
In the simplest case of an assembly of bistable magnetic particles, the elementary hysteresis loop of a particle is rectangular and can be represented by a hysteron (Fig. 1a). In the absence of inter-particle interactions, when any magnetic correlations are irrelevant, the macroscopic hysteresis loop is simply a superposition of projections of magnetic moments of particles onto the field direction, ordered according to the switching events (hysteron thresholds) of individual particles. Then the SFD can be determined by de-constructing the hysteresis loop into distribution of hysterons uniquely linked with particle properties. In this case, FORC method is an inherently accurate technique for its identification. On the other hand, the presence of thermal relaxation and significant inter-particle interactions gives rise to more complex magnetisation reversal mechanisms (Fig. 1b). The emergent magnetic correlations fundamentally transform a macroscopic hysteresis loop and mask the direct information about the intrinsic switching fields of individual particles. The accuracy of the FORC method becomes parameter region-dependent and requires validation against the systematic inverse problem solving framework.
The purpose of the present article is to study the validity range of the FORC method against the large-scale computational data generated from a fully featured model of hard disk drive (HDD) media, which incorporates details of the statistical nature of inter-granular interactions, intrinsic properties of individual grains, and thermal activation. Of particular importance is the role of magnetic correlations and consequent departures from simple hysteron-based model predictions. Gilbert et al. 24 have shown that the introduction of nearest-neighbour correlations strongly modifies the FORC diagram. Here we use a fully-featured kinetic Monte Carlo model with short-and long-ranged interactions to create a complete picture of interaction effects, most importantly the balance between exchange and magnetostatic interactions. We proceed with models of increasing complexity to demonstrate firstly the effects of thermal activation for a non-interacting system. We then proceed to the study of interaction effects using the kinetic Monte-Carlo (kMC) model and a simplified model of correlation effects. By evaluating the inter-granular magnetic correlation function, we demonstrate the direct relationship between the emergence of magnetic correlations and the failure of the FORC methodology to determine the SFD, and establish the criteria for the validity of the FORC method as a quantitative approach for accurate identification of the SFD in HDD magnetic media.

Results
We apply the FORC method to large-scale computational data generated from a fully featured model of HDD media which incorporates details of the statistical nature of intergranular interactions, intrinsic properties of individual grains, and thermal activation (Methods 4). We use a kinetic Monte-Carlo (kMC) model (Methods 4) to computationally reproduce the magnetisation behaviour of the HDD media, and a FORC method as the technique for identification of the underlying SFD. Here we consider realistic media with elongated grains with an aspect ratio (h/d) of 1.17, uniaxial anisotropy (K) with mean value of 7.10 6 erg/cm 3 and 3 degree dispersion of the anisotropy easy axis around the perpendicular direction to the grain plane. The grain height (h) is 10 nm, the mean grain size (grain diameter) (d) is 8.5 nm and the saturation magnetization (Ms) is 700 emu/cm 3 27 . The calculations assume an external field rate of 4.10 4 Oe/s at room temperature (300 K). In all cases studied the intrinsic SFD can be easily calculated in the model by switching off all interactions (Methods 4) and histogramming the switching fields of individual particles along a hysteresis loop. For all the results presented here, between 50 and 60 FORC curves are used with a field step of maximum 100 Oe, although 5 curves only are shown for clarity.
From reversal curve to FORC diagram and to SFD. The FORC method is used as a quantitative tool to investigate the SFD and interaction field distribution in granular materials. It is typically applied to the measurements of macroscopic hysteresis loops. The application of the method contains two main steps. The first step requires measurement of the first order reversal curves (FORC) and their transformation to the so-called FORC diagram ( Fig. 1). In the second step, the FORC diagram is processed such that the undesirable contribution of the inter-particle interaction is removed, which then allows accessing information about the intrinsic SFD.
FORC data, FORC diagram, and the SFD. Figure 1(c) illustrates the measurement protocol used to generate the FORC data. The starting point is the saturation of the sample by applying a large positive applied field. The field is then decreased towards the reversal field, H b , when the field direction is reversed and increased from H b back to positive saturation. This process generates a FORC attached to the major hysteresis loop at the reversal point H b (blue line in Fig. 1(c)). The magnetisation point at an applied field H a > H b along this FORC, denoted as M(H a , H b ), is internal to the major hysteresis loop. As illustrated in Fig. 1 The interpretation of these equations is as follows. The distribution ρ ab in Eq. (1) is defined in terms of the differentiation of the magnetisation M(H a , H b ) attained through general applied fields H a , H b along the hysteresis loop ( Fig. 1(c)), and it is not immediately obvious how it relates to microscopic material properties such as the distribution of intrinsic switching field thresholds of magnetic grains ρ SFD . The key to establishing this link is the notion of a magnetic particle having an elementary rectangular hysteresis loop (RHL) as shown in Fig. 1(a), with the up and down switching thresholds corresponding to the fields H a and H b . Then, Eq. (1) can be interpreted as measuring the fraction of magnetic grains with the switching thresholds H a > H b adding up to the cumulative magnetisation M(H a , H b ) at the field H a after the field excursion from H b . The transformed variables then represent the coercive and the bias fields of such RHLs ( Fig. 1(a)), and the FORC distribution ρ defined in Eq. (2) is the joint probability distribution of H c and H u . Consequently, the SFD defined by Eq. (3) is the distribution of the coercive fields of particles, i.e. their intrinsic switching thresholds.
In the ideal system of isolated magnetic particles represented by RHLs, such as the non-thermal system of non-interacting Stoner-Wohlfarth particles with the anisotropy axes aligned along the field direction, the RHLs have symmetric up and down switching thresholds ± H c and due to the absence of interactions the H u = 0 for all particles. The macroscopic hysteresis loop is a superposition of magnetic states of all particles and, due to the rectangular shape of RHLs, any magnetisation change along the hysteresis loop can occur only at applied fields corresponding to the particle switching thresholds. The differentiation in Eq. (1) filters the contribution from the 'flat' parts of RHLs and as a residual the distribution ρ ab carries an accurate representation of the switching thresholds of particles. In this case, the transformed FORC distribution in Eq. (2) can be shown to be Historically, an elementary RHL of a particle has been referred to as a hysteron in Preisach modelling 18,26 , which has served as a basis for developing the FORC method 16 . The essence of Preisach models is to represent the macroscopic hysteresis loops of materials as a superposition of RHLs with the RHL threshold distribution, termed as a Preisach distribution, defined identically as the ρ ab in Eq. (1) (Methods 4). The uniqueness of identification of the Preisach distribution has been shown to be guaranteed if the macroscopic magnetisation data satisfy the wiping-out and congruency properties 18,19 . Consequently, if the wiping-out and congruency properties are satisfied, the FORC distribution ρ ab is a valid and unique Preisach distribution. Unfortunately, the straightforward interpretation of Eqs (1-3) as given above does not apply in realistic cases when the particles are represented by non-ideal RHLs, the inter-particle interactions are relevant, or in the presence of thermal fluctuations. Moreover, general systems with hysteresis do not always display the wiping-out and congruency properties, and the accuracy and uniqueness of the identification of SFD from the FORC distributions needs to be established with respect to the relevant physical picture and by independent measurement methodologies. Such cases are analysed in detail below.
Effects on imperfect RHLs on FORC diagram. To access the effects of deviations of elementary hysteresis loop of particles from the RHLs on the accuracy of determining the SFD, we applied the kMC model to study the hysteresis loop behaviour of a reduced system of isolated magnetic particles represented as Stoner-Wohlfarth particles (Methods 4 and 4). The intrinsic magnetic properties of particles in the model were set to represent a typical magnetic recording medium (Methods 4), including a 3° misalignment of the particle anisotropy easy axes around the applied field direction, and the driving field rate set to 10 4 Oe/s of a typical experimental MOKE setup, which determined the extent of thermal activation. The inter-particle interactions were turned off. Figure 1(b) shows an example of the computed hysteresis loop of an ensemble of isolated particles, which clearly deviates from RHL behaviour ( Fig. 1(a)). The rounding features are typical of a loop with strong component from thermal activation. The computed macroscopic hysteresis loop of a system of 10000 non-interacting particles with representative FORCs is shown in Fig. 1(c). The FORC diagrams ρ ab and ρ, obtained from this loop by applying Eqs (1) and (2), are shown in Fig. 1(d,e). Note, that given the nature of their transformation, the FORC distributions ρ ab or ρ are related by the 45° rotation of the (H a , H b ) coordinate plane with scaling factor of 1/sqrt(2). Figure 1(e) shows that the FORC distribution ρ is no-longer a straight line ρ ρ δ as in the case of a system of ideal non-interacting particles with RHLs discussed above, and instead has a significant H u component even if the inter-particle interactions are absent. This is due to the particle hysteresis loop rounding seen in Fig. 1(b), when the change of magnetisation along the macroscopic loop no longer occurs only at the switching thresholds of particles, as in the ideal RHL case, but in addition includes a smooth nonlinear component from the rounding effect. The magnetisation data transformation through Eq. (1) then convolutes the residual of the differentiation of this smooth component with the actual distribution of the switching thresholds of particles, which in the FORC diagram becomes manifested as a 'fictitious' H u field distribution ( Fig. 1(e)). This poses a difficulty in the interpretation of the FORC diagram, which appears to suggest the presence of interactions in the system of non-interacting particles. Nevertheless, we find that evaluating the underlying SFD in Eq. . The slightly non-symmetric peak seen in Fig. 1(e), is often observed experimentally in systems with thermal activation.
Interactions: Mean-field correction of the FORC diagram. The effects of inter-particle interactions on the FORC diagram are illustrated in Fig. 2, which shows that the FORC diagram becomes considerably modified by interactions with respect to the non-interacting case in Fig. 1(e). Figure 2(a) shows the FORC diagram of the model based on the mean-field approximation of the full granular model (Methods 4). The mean-field interactions between grains have random strength following from a Gaussian distribution, obtained by histogramming the interaction fields of a full granular system in saturation used as a reference, weighted by the overall magnetisation M of the granular system (Methods 4). The observed FORC diagram has the shape of a rotated 'V' or 'L' as expected 24 .
To apply the FORC method to identify the underlying SFD distribution, it is first necessary to extract the mean-field interaction and recover the non-interacting particle FORC diagram. This can be achieved by introducing a correction factor α, variation of which allows to symmetrise the FORC diagram equivalently to subtracting the average interaction field acting on the system 29,30 . Specifically, varying the mean-field correction factor α transforms the field axes H a and H b in the raw FORC diagram (Fig. 2(a)) to new axes , until obtaining the optimal value of α ≡ α o when the new FORC diagram ρ becomes symmetric around the H u axis and any possible negative regions of ρ < 0 that may result from an over-or under-estimated mean-field correction become eliminated 31 . This procedure is equivalent to the hysteresis loop 'de-shearing' procedure typically applied to extract the effects of demagnetising fields from experimental hysteresis loops. In the ideal case, the optimal value of the correction factor, α o , corresponds to the mean-field interaction strength H inter of the mean field granular model (Methods 4). After applying the mean-field correction, the resulting FORC diagram shown in Fig. 2(b) resembles that of the non-interacting case Fig. 1(e), which allows to calculate the SFD by applying Eq. (3). The values found are consistent with the non-interacting cases within the statistical error corresponding to uncertainty of 5%.
Interactions: magnetic clusters. The mean-field interaction is expected to be an oversimplification as it does not account for the inter-granular magnetic correlations typically present in real systems. The presence of such correlations leads to the emergence of magnetic clusters which influences the accuracy of the FORC method. To begin investigating the magnetic clustering effect, we first consider the mean-field model discussed above reduced to an ensemble of disconnected regions of N g grains. In this 'toy model' the regions act as non-interacting clusters of N g grains interacting via equivalent mean-field-like interactions dependent on the average magnetisation within each cluster (Methods 4). In this model, a switching grain affects only the magnetisation of its own cluster, while the magnetisation of all other clusters in the ensemble remains unaffected, and the hysteresis loop is a superposition of magnetisation jumps from individual clusters. Thus, when the cluster size N g is large, approaching the system size, the behaviour recovers that of a full mean-field system discussed above. On the other hand, as N g decreases the behaviour moves away from being mean-field-like and the macroscopic loop results from a FORC diagram for interacting case with only magnetostatic interaction calculated using mean-field approach (a) and grain-grain interaction (c). The corresponding FORC diagram after mean interaction field is removed are given in (b,d) for the mean-field and grain-grain interaction models respectively. The radial magnetisation correlation function for the grain-grain interaction model is given as the inset in (c).
combined contribution of an increasing number of elementary hysteresis loops of individual clusters available in the system. These elementary loops of individual clusters have shape deviating from the RHLs, which is expected to reduce the accuracy of the FORC method. Figure 3 shows analysis of five ensembles of uniform clusters of variable N g = 4, 5, 10, 100, 500. Applying the cluster model (Methods 4) combined with the kinetic Monte-Carlo solver (Methods 4) we first computed the macroscopic hysteresis loops with FORCs for every ensemble. Then we transformed the FORC data to the underlying FORC diagram by applying Eqs (1) and (2), applied the mean-field correction α 0 to remove the interactions as discussed above -which is a standard procedure used in the practical FORC method, and computed the SFD from the corrected FORC diagram using Eq. (3). As expected, the results of extracting the SFD are accurate when N g approaches the full system size, while the accuracy of the FORC method reduces with the decreasing cluster size N g . When N g is small, the clusters contain only small numbers of grains relative to the full system size, and there are many clusters contributing to the overall macroscopic hysteresis loop. There are two main sources of error expected to contribute to the loss of accuracy of the FORC method: (1) the mean-field correction is no-longer accurate as in the mean-field model of a full granular system, and (2) distorted RHL shape of elementary hysteresis loops of individual grains, due to correlated behaviour inside each cluster. Examples of the obtained raw FORC diagrams prior applying the mean-field correction and the orignal FORC data are shown in the insets i-iv in Fig. 3. The mean-field like nature of the interaction in the cluster model (Methods 4) results in an equivalent effective shift of the switching thresholds of the grains in each cluster, which results in the observed segmentation of the 'V'-like shape FORC diagram into distinct regions along each branch. The number of these regions per branch corresponds to the number of grains per cluster, the 'V' shape of the arrangement of the segments reflects the interaction induced symmetry breaking of the up and down intrinsic switching thresholds of grains where the separation between the segments corresponds roughly to the magnitude of the mean interaction field H inter in the model. The interpretation of this FORC diagram is consistent with the recent work, where analogous segmentation effects have been studied in terms of a different model with the nearest neighbour grain interactions 24 . Increasing N g in clusters results in the increased density of segments in the FORC diagram until gradually reproducing the FORC diagram of the mean field model in Fig. 2(a). We note that in a real system the correlations vary with applied magnetic field and clusters with arbitrary size and shape are formed. This will blur the segmentation so the FORC diagrams look more conventional, as also found by Gilbert et al. 24 . However, the The cluster size in inverse proportional with the degree of correlation in the system. The larger the cluster size the closer is the model to a mean-field-like model, which is completely uncorrelated system and FORC method can be applied successfully. By increasing the correlation, the FORC method is underestimating the σ SFD . The result from the HDD model is also included (blue). Example of FORC diagram and original FORC data for different cluster size are illustrated in the insets i-iv.
Scientific RepoRts | 7:45218 | DOI: 10.1038/srep45218 calculated SFD still shows a significant error. Next we give calculations for a range of realistic parameter values which show large systematic errors in the calculated SFD.

Applicability of the FORC method to a full granular model of recording media.
To investigate the accuracy of the FORC method in determining the SFD in realistic magnetic recording materials we use the full granular kinetic Monte-Carlo model with exchange and magnetostatic interactions (Methods 4) to simulate the underlying hysteresis loops and FORCs. Such general interactions introduce magnetic correlations between grains, which lead to correlated behaviour when magnetic grains begin to switch in unison in clusters of size equal to the characteristic correlation length. This leads to magnetisation jumps (Barkhausen noise) along the hysteresis loop, analogous to the case of the cluster model discussed above. Figure 2(c) shows the corresponding FORC diagram.
The inset in the figure shows the radial correlation function, suggesting the presence of significant short range grain-grain correlations in a typical recording medium which are absent in the non-interacting and mean-field granular systems. To remove the contribution from interactions, we first subtract the mean-field correction after finding the optimal α o as discussed above, as is typically done in practical applications of the FORC method. The corrected FORC diagram shown in Fig. 2(d) deviates from the non-interacting case shown in Fig. 1(e), which is due to the fact that the mean-field interaction mis-represents the full exchange and magnetostatic interactions. Consequently, applying Eq. (3) we find that the FORC method underestimates the σ SFD by as much as 60%. Thus the presence of significant magnetic correlations results in the loss of accuracy of the FORC method. The question of main interest is to understand the relationship between the extent of correlations and the accuracy of the SFD determined by the FORC method.
To study this issue in simulations, we systematically varied the strength of exchange and magnetostatic interactions, in each case computing the underlying radial pair correlation function between the grains (Methods 4), and evaluated the reference SFD directly by histogramming the intrinsic field thresholds of grains during switching for comparison with the SFD obtained by the FORC method through Eqs (1-3). Figure 4(a) shows the dependence of the maximum value of the correlation function on the strength of exchange and magnetostatic fields. Representative FORC diagrams before applying the mean-field correction are shown in the insets (i)-(v). Magnetic correlations increase with the strength of one of the interaction types increasing relative to the other, while they remain negligible in the weakly interacting case or in the interaction compensating region corresponding to the region with similar total magnitudes of exchange and magnetostatic interactions. The contour lines quantify the correlation strength. Figure 4(b) shows the corresponding relative accuracy of the SFD determined by the FORC method, measured relative to the SFD determined directly from the kMC model. The comparison of Fig. 4(a,b) reveals close agreement between the correlation strength and the accuracy of the FORC method for determining the SFD. Errors can also be attributed to the fact that the model is thermal and RHL are not perfect, or the effects of the slight misalignment of anisotropy axes of grains combined with the interactions. However, as shown in Fig. 1(e) for the thermal effects, these factors are relatively small and the largest discrepancy is caused by the interactions and specifically by the interaction-induced spatial correlations. The accuracy of the FORC method is the highest in the weakly correlated interaction regions. Nevertheless, depending on the required accuracy of determination of the SFD, Fig. 4(b) indicates the range of parameter space in which this can be achieved. The FORC method is limited to very small field (up to 1200 Oe) considering a deviation of 10% from the expected value of σ SFD . Finally we map the deviation of the SFD from FORC and combine the results with magnetisation correlation data to draw a validity diagram for using FORC as a quantitative tool.

Discussion
Our study presents the analysis of the accuracy and the limitations of hysteron-based analysis of the experimental data obtained by the FORC method for quantitative determination of the SFD. This is achieved by using as a benchmark a succession of computational models of increased complexity starting from the system of non-interacting particles towards the realistic full model of magnetic granular media in magnetic recording, which includes exchange and magnetostatic interactions and various relevant sources of the material disorder. In terms of the analysis, similar to Pike et al. 16 we make the distinction between the raw FORC data, the FORC diagram and the usual interpretation of the FORC data based on the Preisach model. Using model calculations we show that, while the FORC diagram in principle contains information about the SFD and the interactions, application of the RHL interpretation does not reliably deconvolve the SFD and interaction effects. This is attributed to the spatial magnetisation correlations which are an important feature of many materials, including magnetic recording media, and which are not included in the RHL approach.
We reveal that the applicability of the FORC method for the quantitative analysis of the SFD is limited to the parameter range where inter-particle spatial correlations are insignificant, i.e. when exchange and magnetostatic interactions are weak or when they compensate. The accuracy of the FORC method decreases in the presence of significant correlations resulting from the correlated switching with multiple grains reversing their magnetic state in unison at a given field threshold. These correlated grains behave as a single entity thus hiding any information about the intrinsic switching thresholds of individual particles into the correlated reversal. The information cannot be recovered by the differentiation of the first order magnetisation reversal curves in the way of the FORC method simply because Eq. (1) no longer provides access to the switching thresholds of individual RHLs of particles but instead provides access to the fields corresponding to the magnetisation jumps along first order reversal curves, which are equivalent to the switching thresholds of the correlated particle clusters as a single entity and not as individual particles in the cluster.
Moreover, general hysteresis loops do not satisfy the wiping out and congruency properties. Then the FORC diagrams cannot be interpreted as Preisach distributions and no longer guaranteeing unique SFD 16,18,25,26 , which necessitate careful analysis in establishing its physical relevance. Recovering the intrinsic switching thresholds of individual particles from the first order reversal curve data then requires further deconvolution based on the Scientific RepoRts | 7:45218 | DOI: 10.1038/srep45218 refined models capable of accounting for the detailed structure of inter-particle interactions. Our work applies such fine-scale models and based on the evaluation of microscopic correlation functions establishes quantitatively the range of validity of the FORC method for determining the SFD with relevance to magnetic recording media (Fig. 4).
Generally speaking, to identify the SFD in the material parameter range beyond the applicability of the FORC method requires inverse problem solving techniques based on the physically realistic models, which allow reproducing the relevant correlated switching of particles. Moreover, besides identifying accurate models suitable for interpreting the experimental data, such methods also require establishing uniqueness properties of the identified solutions. We have implemented a direct approach employing optimisation techniques based on the grid-search method 32 to fit the full recoding model (Methods 4) to the computed hysteresis loop data, and uniquely recovered the expected SFD in the entire parameter range. Thus, the most reliable, albeit computationally expensive approach, seems to be to essentially carry out by a direct fit to the experimental FORC data using a microscopic approach, including the detailed calculation of the interactions such as presented here for the specific example of perpendicular recording media.

Methods
Full interacting model recording media. The system consists of N Stoner-Wohlfarth grains, where the volume (V) and geometry of the grains is generated by using a Voronoi construction. The energy of a system of N grains is:  where the first term is the uniaxial anisotropy terms with = K K k i i i being the uniaxial anisotropy vector and V i the volume of a particle i, M s the saturation magnetisation, and = m m M / i i s the particle moment normalised to unity. The values of K i , k i , and V i are drawn from random distributions relevant to modern granular magnetic recording materials, as described below. The second term represents the Zeeman term describing the interaction of grains with the applied field H ap .
The third term in Eq. (4) describes the exchange interaction between the nearest neighbour grains. The exchange interaction in granular materials for magnetic recording is dependent on the extent of the grain boundary and is of randomised character, which can be expressed as where H exch is the mean strength of the exchange interaction field, J ij is the fractional exchange constant between the adjacent grains i and j with L ij being the length of the connecting boundary, A i is the area of the grain i, and 〈⋅〉 represent averages over all pairs of grains. The last term in Eq. (4) represents the magneto-static interaction between the grains and is represented as ij . The contribution to the magneto-static interaction field H mag ij is performed by a direct integration of the magneto-static surface charge 33 . The evaluation of H mag ij by full integration over the surface charge accounts for the correction resulting from the dipolar interaction over-estimating the magneto-static interaction in the proximity of a grain 34 . Both exchange and magnetostatic interactions in Eq. (4) are dependent on the size and shape of grains, and on the inter-granular distance.
Approximations of the full model. Various levels of reduction of the full model can be introduced as follows.
Non-interacting model approximation of recoding media. In the non-interacting model, the definition of the system energy reduces from Eq. (4) to: The kinetic Monte-Carlo modelling of this system allows to study thermal relaxation aspects in ensemble of non-interacting Stoner-Wohlfarth particles, and may serve as a reference for gauging the effects of interactions in the full interacting model. Mean-field model of recoding media. In the mean-field model the interactions between magnetic grains are introduced in a uniform ways. The energy expression given in Eq. (4) can be reduced: where the symbol m k implies averaging over all grains in the system, i.e. average magnetisation at a given field H ap , H inter i is the random interaction field given by Gaussian distribution with mean H inter and standard deviation σ inter . We found that Gaussian distribution represents well the distribution of interaction fields in the full HDD model at saturating fields. This allows us to calibrate the mean-field interaction strength to be consistent with the full model at saturating fields, which is a point used as a reference.
Cluster ensemble model of recording media. In the cluster model, the full model is divided into clusters of N g grains, with grains inside a j-th cluster interacting via a mean-field like interaction, while the clusters being non-interacting. The full energy expression given in Eq. (4)  where the symbol ∑ ∈ i j implies that the summations occur through the magnetic grains with the cluster j, and m k j is the average magnetisation of the cluster j. The interaction field H inter i is defined identically as in the mean-field case in Section 4.
Modelling hysteresis with thermal activation: Kinetic Monte-Carlo approach. The thermal fluctuation and external field driven magnetisation behaviour of interacting magnetic particles as described in the model in the Methods Section 0.1 is modelled by using kinetic Monte-Carlo approach 35 . The time dependent transition for a particle moment m i to switch between the up ('1') and down ('2') states is τ = − − P t 1 exp( / ) i i , where the relaxation time constant τ i is a reciprocal sum of the transition rates τ + i and τ − i dependent on the energy barriers ∆E i 1,2 seen from the '1' and '2' states via the standard Néel-Arrhenius law 37  where j = x, y, z and m R ( ) j , + m R r ( ) j are pairs of of grains separated by a distance r. The correlation data plot in Fig. 4(b) shows the correlation function C Z (r). FORC method. The measurement protocol to produce a first order reversal curve (FORC) begins by first applying a large field to saturate the sample, then decreasing the field to a certain value H b . From this point, the FORC is obtained by increasing the field back to saturation. The magnetisation is recoded at fields H a along the FORC at the reversal field H b , H a > H b . The FORC diagram is then evaluated using Eqs (1) and (2), from which the SFD can be calculated using Eq. (3). The numerical data, M(H a , H b ) is fitted by two-variable polynomial of the smoothing factor, SF, to a discrete set of magnetisation values. The FORC diagram is obtained similar to the FORCinel software, where the smoothing is done using Locally Weighted Polynomial Regression Method 38 . The FORC is re-interpolated after the deshearing, a smoothing factor of 3 is used and the same resolution is keep between the original data and the final-corrected results 31 .