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

###### Abstract

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 magnetization and suggest a more general approach.

## Introduction

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.

In this work we consider the exemplar problem of identification of the switching field distribution (SFD) in magnetic particulate and granular systems. The SFD carries information about the intrinsic conditions for magnetisation reversal of individual magnetic particles and, for example, is a crucial characteristic determining the quality of high density magnetic media for current and future hard disk technologies, such as based on the bit patterned media or granular materials considered for HAMR [5, 6]. This problem is especially challenging due to the strength and complexity of the interparticle interactions. The recently developed identification schemes to extract the SFD based on the inverse problem solving approach include applications to assemblies of magnetic nanoparticles[7], and method for granular materials relevant in magnetic recording [8, 9, 10, 11] [12, 13, 14], which was motivated by the earlier methodology [15]. Simpler approaches to identify the SFD with variable degree of consistency were based on the differentiation of hysteresis loop ‘de-sheared’ to remove the contribution from magneto-static interactions [16], methods based on Preisach models [17], and methods based on analysing the transformed first order reversal curves (FORC), i.e. magnetisation curves generated by reversing the external magnetic field starting from a point on a major hysteresis loop branch (Fig. 1). The FORC methods are equivalent to the classical Preisach modelling if the measurement data display microscopic memory of states of magnetic particles after the external field excursion (wiping-out property) and a minor hysteresis loop congruency [18, 19].

The FORC methods are used broadly as a tool for qualitative, and in some cases quantitative, description of general magnetic characteristics of magnetic systems, such as of distributions of magnetic properties, mixed magnetic phases[20], clustering and long-range ferromagnetic state, magnetic characterisation of geological mixtures and minerals and the differences in magnetization reversal mechanisms[21, 22, 23, 24]. The attractiveness of the FORC method is in its simplicity and its straightforward application to a wide range of systems displaying hysteresis. The accuracy of determining the SFD quantitatively in various classes of systems is presently under intensive critical discussion [25, 26], and quantifying its range of validity, and understanding the microscopic reasons for its breakdown, is of broad interest with implications beyond the exemplar magnetic hysteresis considered here.

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 0.1). We use a kinetic Monte-Carlo (kMC) model (Methods 0.3) 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 () of 1.17, uniaxial anisotropy () with mean value of erg/cm and 3 degree dispersion of the anisotropy easy axis around the perpendicular direction to the grain plane. The grain height () is 10 nm, the mean grain size () is 8.5 nm and the saturation magnetization () is emu/cm. The calculations assume an external field rate of Oe/s at room temperature (300K). In all cases studied the intrinsic SFD can be easily calculated in the model by switching off all interactions (Methods 0.2.1) and histogramming the switching fields of individual particles along a hysteresis loop.

### 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. Fig. 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, , when the field direction is reversed and increased from back to positive saturation. This process generates a FORC attached to the major hysteresis loop at the reversal point (blue line in Fig. 1(c)). The magnetisation point at an applied field along this FORC, denoted as , is internal to the major hysteresis loop. As illustrated in Fig. 1(c), at any value of in the hysteresis region, there is an entire family of such internal magnetisation points distinguished by the reversal field of their corresponding FORCs. The FORC data are then analysed by computing the numerical second-order derivative of the functional dependence with respect to the applied field and :

(1) |

where is the saturation magnetisation of the material. It is next conventional to transform by introducing new variables and such that and , which leads to the FORC distribution represented as:

(2) |

from which the SFD can be obtained by a straightforward integration over the variable :[27]

(3) |

The interpretation of these equations is as follows. The distribution in Eq. (1) is defined in terms of the differentiation of the magnetisation attained through general applied fields , 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 . 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 and . Then, Eq. (1) can be interpreted as measuring the fraction of magnetic grains with the switching thresholds adding up to the cumulative magnetisation at the field after the field excursion from . The transformed variables and 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 and . 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 and due to the absence of interactions the 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 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 , where is the Dirac delta-function and the statistical distribution of coercive fields of RHLs of particles, which according to Eq. (3) gives the SFD directly as .

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 in Eq. (1) (Methods 0.5). 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 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 0.2.1 and 0.3). The intrinsic magnetic properties of particles in the model were set to represent a typical magnetic recording medium (Methods 0.4), including a 3 misalignment of the particle anisotropy easy axes around the applied field direction, and the driving field rate set to Oe/s of a typical experimental MOKE setup, which determined the extent of thermal activation. The inter-particle interactions were turned off.

Fig. 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 and , obtained from this loop by applying Eqs. (1) and (2), are shown in Figs. 1(d) and (e). Note, that given the nature of their transformation, the FORC distributions or are related by the 45 rotation of the coordinate plane. Fig. 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 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’ 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. (3) based on this FORC diagram actually yields the accurate SFD as a result of the reflection symmetry of the FORC diagram around the axis, when the component of the FORC distribution simply integrates to unity after factorisation of in Eq. (3). 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 Figs. 2, which shows that the FORC diagram becomes considerably modified by interactions with respect to the non-interacting case in Fig. 1(e). Fig. 2(a) shows the FORC diagram of the model based on the mean-field approximation of the full granular model (Methods 0.2.2). 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 of the granular system (Methods 0.2.2). 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[28]. Specifically, varying the mean-field correction factor transforms the field axes and in the raw FORC diagram (Fig. 2(a)) to new axes and , until obtaining the optimal value of when the new FORC diagram becomes symmetric around the axis and any possible negative regions of that may result from an over- or under-estimated mean-field correction become eliminated. 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, , corresponds to the mean-field interaction strength of the mean field granular model (Methods 0.2.2). 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 grains. In this ‘toy model’ the regions act as non-interacting clusters of grains interacting via equivalent mean-field-like interactions dependent on the average magnetisation within each cluster (Methods 0.2.3). 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 is large, approaching the system size, the behaviour recovers that of a full mean-field system discussed above. On the other hand, as decreases the behaviour moves away from being mean-field-like and the macroscopic loop results from a 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.

Fig. 3 shows analysis of five ensembles of uniform clusters of variable , 5, 10, 100, 500. Applying the cluster model (Methods 0.2.3) combined with the kinetic Monte-Carlo solver (Methods 0.3) 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 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 approaches the full system size, while the accuracy of the FORC method reduces with the decreasing cluster size . When 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 are shown in the insets i-iv in Fig. 3. The mean-field like nature of the interaction in the cluster model (Methods 0.2.3) 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 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 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).

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 0.1) 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. Fig. 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 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 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 0.6), 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). Fig. 4(a) shows the dependence of the maximum value of the correlation function on the strength of exchange and magnetostatic fields. Representative FORC diagrams after 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. Fig. 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 Figs. 4(a) and (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 1200Oe) considering a deviation of 10% from the expected value of . Finally we map the deviation of the SFD from FORC and combine the results with magnetization 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 different modes of failure of the FORC method 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 magnetization 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 [26, 18, 16, 25], 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 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 [29] to fit the full recoding model (Methods 0.1) 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

### 0.1 Full interacting model recording media

The system consists of Stoner-Wohlfarth grains, where the volume () and geometry of the grains is generated by using a Voronoi construction. The energy of a system of grains is:

(4) |

where the first term is the uniaxial anisotropy terms with being the uniaxial anisotropy vector and the volume of a particle , the saturation magnetisation, and the particle moment normalised to unity. The values of , , and 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 .

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 with the locally varying exchange field :[30]

(5) |

where is the mean strength of the exchange interaction field, is the fractional exchange constant between the adjacent grains and with being the length of the connecting boundary, is the area of the grain , 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 . The contribution to the magneto-static interaction field is performed by a direct integration of the magneto-static surface charge [31]. The evaluation of 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.[32] Both exchange and magnetostatic interactions in Eqs. (4) are dependent on the size and shape of grains, and on the inter-granular distance.

### 0.2 Approximations of the full model

Various levels of reduction of the full model can be introduced as follows.

#### 0.2.1 Non-interacting model approximation of recoding media

In the non-interacting model, the definition of the system energy reduces from Eq. (4) to:

(6) |

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.

#### 0.2.2 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:

(7) |

where the symbol implies averaging over all grains in the system, i.e. average magnetisation at a given field , is the random interaction field given by Gaussian distribution with mean and standard deviation . 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.

#### 0.2.3 Cluster ensemble model of recording media

In the cluster model, the full model is divided into clusters of grains, with grains inside a -th cluster interacting via a mean-field like interaction, while the clusters being non-interacting. The full energy expression given in Eq. (4) reduces to a sum through individual clusters as with:

(8) |

where the symbol implies that the summations occur through the magnetic grains with the cluster , and is the average magnetisation of the cluster . The interaction field is defined identically as in the mean-field case in Section 0.2.2.

### 0.3 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[33, 34]. The effective local fields of particles are given by Eqs. (4) as . The time dependent transition for a particle moment to switch between the up (‘1’) and down (‘2’) states is , where the relaxation time constant is a reciprocal sum of the transition rates and dependent on the energy barriers seen from the ‘1’ and ‘2’ states via the standard Néel-Arrhenius law[35]: . The is the Boltzmann constant and the temperature. According Eq. (4), the depend on the intrinsic particle properties, such as and .

### 0.4 Simulation parameters of realistic recording media

Throughout this study we consider a thin film system, with elongated grains (1.17 aspect ratio), log-normal volume distribution (33%) and log-normal anisotropy distribution (5%). The uniaxial anisotropy has a 3 dispersion of easy axis around the axis perpendicular to the film. The system properties are: mean anisotropy erg/cm, saturation magnetisation emu/cm, grain height nm and the mean grain size nm. The calculations are done for an external field rate of Oe/s at room temperature 300K.

### 0.5 Rectangular hysteresis loop (RHL) model: Preisach modelling

If a granular system can be viewed as a collection of grains having rectangular hysteresis loops (RHL), with coercive field and bias field given by probability distribution, then the macroscopic hysteresis loop of the system can be obtained as a superposition of the RHLs and magnetisation represented as:

(9) |

where is the saturation magnetisation, , , and are the integration limits dependent on the applied field along the FORC attached do decreasing major hysteresis loop at the reversal field , i.e. . Applying the Leibniz integral rule to differentiate the integral we obtain:

(10) |

Given that Eq. (9) is inherently a superposition from switching events of individual grains, Eq. (10) establishes the relation between the applied fields and , and the intrinsic switching thresholds of particles, which can be labeled equivalently as (threshold of a grain flipping up along the FORC at the field ) and (threshold for a grain flipping down before generating FORC at the field ). Given that and (Fig. 1(a)), the above equation can be rewritten as:

(11) |

which agrees with the definition of the FORC distribution given in Eqs. (1) and (2). If the system displays the wiping-out and congruency properties, Eq. (9) can be shown to be a unique Preisach distribution associated with the granular system represented by magnetisation .

### 0.6 Magnetic correlation

To investigate the coupling between grains due to correlated behaviour, we computed the radial correlation function as following:

(12) |

where and , are pairs of of grains separated by a distance . The correlation data plot in Fig. 4(b) shows the correlation function .

### 0.7 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 . From this point, the FORC is obtained by increasing the field back to saturation. The magnetisation is recoded at fields along the FORC at the reversal field , . The FORC diagram is then evaluated using Eqs. (1) and (2), from which the SFD can be calculated using Eq. (3).

## References

- 1. Tarantola, A. Popper, Bayes and the inverse problem. Nat. Phys. 2, 492–494 (2006). URL http://www.nature.com/doifinder/10.1038/nphys375.
- 2. Bertero, M. & Boccacci, P. Introduction to inverse problems in imaging (CRC press, 1998).
- 3. Sarvas, J. Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem. Physics in medicine and biology 32, 11 (1987).
- 4. Tarantola, A. Inverse Problem Theory and Methods for Model Parameter Estimation. Other titles in applied mathematics (Society for Industrial and Applied Mathematics, 2005). URL https://books.google.co.uk/books?id=kEboSYSU-nAC.
- 5. Hongzhi Yang et al. Method for Measurement of Switching-Field Distribution of Heat-Assisted Magnetic Recording Media. IEEE Magn. Lett. 6, 1–4 (2015). URL http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=6947939%****␣main2.tex␣Line␣375␣****http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=7079384.
- 6. Weller, D. et al. A HAMR Media Technology Roadmap to an Areal Density of 4 Tb/in. IEEE Trans. Magn. 50, 1–8 (2014). URL http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=6589966.
- 7. Ruta, S. I. Computational modelling of interaction effects in Fe 3 O 4 nanoparticle systems for comparison with experiments. Msc thesis, Dep. Physics, Univ. York, UK (2013).
- 8. Liu, Y., Dahmen, K. & Berger, a. Determination of intrinsic switching field distributions in perpendicular recording media: Numerical study of the H(M,M) method. Phys. Rev. B 77, 054422 (2008). URL http://prb.aps.org/abstract/PRB/v77/i5/e054422$\$nhttp://link.aps.org/doi/10.1103/PhysRevB.77.054422.
- 9. Hovorka, O., Liu, Y., Dahmen, K. & Berger, A. On the ability to determine intrinsic switching field distributions from hysteresis loops in the partially correlated magnetization reversal regime. Journal of Magnetism and Magnetic Materials 322, 459–468 (2010). URL http://linkinghub.elsevier.com/retrieve/pii/S0304885309010002.
- 10. Hovorka, O., Pressesky, J., Ju, G., Berger, A. & Chantrell, R. W. Distribution of switching fields in magnetic granular materials. Appl. Phys. Lett. 101 (2012).
- 11. Pisana, S. et al. Effects of grain microstructure on magnetic properties in FePtAg-C media for heat assisted magnetic recording. J. Appl. Phys. 113, 043910 (2013). URL http://scitation.aip.org/content/aip/journal/jap/113/4/10.1063/1.4788820.
- 12. Hauet, T. et al. Reversal mechanism, switching field distribution, and dipolar frustrations in Co/Pt bit pattern media based on auto-assembled anodic alumina hexagonal nanobump arrays. Phys. Rev. B 89, 174421 (2014). URL http://link.aps.org/doi/10.1103/PhysRevB.89.174421.
- 13. Tabasum, M. R. et al. Magnetic force microscopy study of the switching field distribution of low density arrays of single domain magnetic nanowires. J. Appl. Phys. 113 (2013). 1304.6441.
- 14. Wang, T. et al. Magnetic behavior in an ordered Co nanorod array. Nanotechnology 19, 455703 (2008). URL http://www.ncbi.nlm.nih.gov/pubmed/21832792.
- 15. Tagawa, I. & Nakamura, Y. Relationships between high density recording performance and particle coercivity distribution. IEEE Trans. Magn. 27, 4975–4977 (1991). URL http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=278712.
- 16. Pike, C. R., Roberts, A. P. & Verosub, K. L. Characterizing interactions in fine magnetic particle systems using first order reversal curves. J. Appl. Phys. 85, 6660 (1999). URL http://scitation.aip.org/content/aip/journal/jap/85/9/10.1063/1.370176.
- 17. Mayergoyz, I. D. Mathematical Models of Hysteresis (Springer New York, New York, NY, 1991). URL http://link.springer.com/10.1007/978-1-4612-3028-1.
- 18. Mayergoyz, I. D. Hysteresis models from the mathematical and control theory points of view. J. Appl. Phys. 57, 3803 (1985). URL http://scitation.aip.org/content/aip/journal/jap/57/8/10.1063/1.334925.
- 19. Stancu, A., Pike, C., Stoleriu, L., Postolache, P. & Cimpoesu, D. Micromagnetic and Preisach analysis of the First Order Reversal Curves (FORC) diagram. J. Appl. Phys. 93, 6620 (2003). URL http://link.aip.org/link/JAPIAU/v93/i10/p6620/s1{&}Agg=doihttp://scitation.aip.org/content/aip/journal/jap/93/10/10.1063/1.1557656.
- 20. a.P. Roberts, Pike, C. & Verosub, K. First-order reversal curve diagrams: a new tool for characterizing the magnetic properties of natural samples. J. Geophys. Res. 105, 28461–475 (2000). URL http://www.agu.org/pubs/crossref/2000/2000JB900326.shtml.
- 21. Pike, C. R., Roberts, A. P., Dekkers, M. J. & Verosub, K. L. An investigation of multi-domain hysteresis mechanisms using FORC diagrams. Phys. Earth Planet. Inter. 126, 11–25 (2001).
- 22. Muxworthy, A. R. Assessing the ability of first-order reversal curve (FORC) diagrams to unravel complex magnetic signals. J. Geophys. Res. 110, B01105 (2005). URL http://doi.wiley.com/10.1029/2004JB003195.
- 23. Roberts, A. P., Heslop, D., Zhao, X. & Pike, C. R. Understanding fine magnetic particle systems through use of first-order reversal curve diagrams. Rev. Geophys. 52, 557–602 (2014). URL http://doi.wiley.com/10.1002/2014RG000462.
- 24. Gilbert, D. A. et al. Quantitative decoding of interactions in tunable nanomagnet arrays using first order reversal curves. Sci. Rep. 4, 4204 (2014). URL http://www.nature.com/srep/2014/140226/srep04204/full/srep04204.html.
- 25. Dobrota, C.-I. & Stancu, A. What does a first-order reversal curve diagram really mean? A study case: Array of ferromagnetic nanowires. J. Appl. Phys. 113, 043928 (2013). URL http://scitation.aip.org/content/aip/journal/jap/113/4/10.1063/1.4789613.
- 26. Dobrotă, C.-I. & Stancu, A. Tracking the individual magnetic wires’ switchings in ferromagnetic nanowire arrays using the first-order reversal curves (FORC) diagram method. Phys. B Condens. Matter 457, 280–286 (2015). URL http://linkinghub.elsevier.com/retrieve/pii/S0921452614007947.
- 27. Winklhofer, M. & Zimanyi, G. T. Extracting the intrinsic switching field distribution in perpendicular media: A comparative analysis. Journal of Applied Physics 99, 08E710 (2006).
- 28. Papusoi, C., Srinivasan, K. & Acharya, R. Study of grain interactions in perpendicular magnetic recording media using first order reversal curve diagrams. J. Appl. Phys. 110, 083908 (2011). URL http://scitation.aip.org/content/aip/journal/jap/110/8/10.1063/1.3652846.
- 29. Hovorka, O., Liu, Y., Dahmen, K. A. & Berger, A. Simultaneous determination of intergranular interactions and intrinsic switching field distributions in magnetic materials. Appl. Phys. Lett. 95 (2009). URL http://scitation.aip.org/content/aip/journal/apl/95/19/10.1063/1.3263732.
- 30. Peng, Y. et al. Cluster size and exchange dispersion in perpendicular magnetic media. Journal of Applied Physics 109, 123907 (2011). URL http://scitation.aip.org/content/aip/journal/jap/109/12/10.1063/1.3596806.
- 31. Newell, A. A generalization of the demagnetizing tensor for nonuniform magnetization. Journal of Geophysical … 98, 9551–9555 (1993). URL http://onlinelibrary.wiley.com/doi/10.1029/93JB00694/full.
- 32. Liu, Y., Hovorka, O., Berger, A. & Dahmen, K. a. Effects of nonuniform exchange and magnetostatic interactions on the determination of intrinsic switching field distributions. Journal of Applied Physics 105, 123905 (2009). URL http://scitation.aip.org/content/aip/journal/jap/105/12/10.1063/1.3149696.
- 33. Chantrell, R., Walmsley, N., Gore, J. & Maylin, M. Calculations of the susceptibility of interacting superparamagnetic particles. Physical Review B 63, 024410 (2000). URL http://link.aps.org/doi/10.1103/PhysRevB.63.024410.
- 34. Ruta, S., Chantrell, R. & Hovorka, O. Unified model of hyperthermia via hysteresis heating in systems of interacting magnetic nanoparticles. Sci. Rep. 5, 9090 (2015). URL http://www.nature.com/articles/srep09090http://www.ncbi.nlm.nih.gov/pubmed/25766365. 1412.3814.
- 35. Néel, L. Théorie du traînage magnétique des ferromagnétiques en grains fins avec applications aux terres cuites. Ann. géophys (1949).

## Acknowledgements

We are grateful to Prof. A. Berger, Prof. A. Stancu, Prof. G.T. Zimanyi, M. Strungaru and A. Meo for helpful discussions and comments. This work made use of the facilities of N8 HPC provided and funded by the N8 consortium and EPSRC (Grant No. EP/K000225/1) co-ordinated by the Universities of Leeds and Manchester and the EPSRC Small items of research equipment at the University of York ENERGY (Grant No. EP/K031589/1).