Role of atomic-scale thermal fluctuations in the coercivity

The microscopic mechanism of coercivity at finite temperature is a crucial issue for permanent magnets. Here we present the temperature dependence of the coercivity of an atomistic spin model for the highest-performance magnet Nd2Fe14B. For quantitative analysis of the magnetization reversal with thermal fluctuations, we focus on the free energy landscape as a function of the magnetization. The free energy is calculated by the replica-exchange Wang–Landau method. This approach allows us to address a slow nucleation problem, i.e., thermal activation effects, in the magnetization reversal. We concretely observed that the thermal fluctuations lead to a downward convexity in the coercivity concerning the temperature. Additionally, through analyzing the microscopic process of the thermal activation (nucleation), we discover the activation volume is insensitive to a magnetic field around the coercivity. The insensitivity explains the linear reduction of the free energy barrier by the magnetic field in the nucleation process.


INTRODUCTION
Demands on high-performance permanent magnets have become an imperative subject of the save-energy technology for sustainable environments 1,2 . To improve the performance of magnets, there is often a great desire to understand the temperature dependence of coercivity. Experimental and theoretical researches have been extensively carried out to investigate the microscopic understanding of the coercivity mechanism from atomic-scale magnetic structure [3][4][5][6][7][8][9][10] . Recently, a scanning hard X-ray microprobe was employed to visualize the nucleation and domain wall motion of magnetization reversal processes from the grain boundary interfaces in the highest performance magnet Nd 2 Fe 14 B 8 . Until now, a study on the microscopic approach to temperature effects on coercivity, especially a slow nucleation process, has not been reported. The nucleation leads to the observation time dependence of the coercivity. This is phenomenologically known as coercivity reduction by thermal activation effect [11][12][13][14][15] .
Large reduction of coercivity from the theoretically expected value is an essential current problem which has been called "Discrepancy to theory" by Kronmüller et al. 16 . In this regard, extensive works have been done to elucidate the reason for the discrepancy 17-21 , but not yet clear from the atomistic viewpoint. The present paper studies the reduction due to thermal fluctuations. The magnetization reversal is a transition from a metastable state to a stable state by overcoming energy barriers under a reverse magnetic field. In this overcoming process, a reverse nucleus is formed. Under thermal fluctuations, the nucleation occurs stochastically, that is, thermal-activated relaxation. Thus the coercivity, i.e., the threshold field, depends on the observation time. In the permanent magnet applications, the metastable state is required to be stable for the duration of the order of a second. To handle such a slow relaxation process, a method that uses the energy landscape has been developed, also known as the minimum energy path (MEP) 22 . For the magnetization reversal, MEP has been studied in the continuum approximation models with parameters with respect to a given temperature 15,18,19,23 . However, the aforementioned method has not considered thermal fluctuations. In the present work, we propose an approach to evaluate the temperature dependence of the coercivity from a microscopic viewpoint by using an atomistic model. In contrast to the MEP, the present technique directly handles thermal fluctuations and also the atomistic information in the magnetization reversal process. In particular, we focus on the Neodymium (Nd) magnet Nd 2 Fe 14 B 24 .
The coercivity is a non-equilibrium concept concerning the collapse of a metastable state, but not expressed by an equilibrium expectation value. Thus, we need a dynamical model. The relaxation time is estimated from the free-energy barrier F B by using the form (like Arrhenius law): where τ 0 is the reciprocal of the attempt frequency, which is typically set to 10 −11 s 12 . Then, for the observation time of 1 s, setting τ = 1 s, the barrier height is given by F B = 25.3 k B T. The free energy is given by the Nd 2 Fe 14 B atomistic spin Hamiltonian, which was recently developed and successfully reproduced static properties of the Nd 2 Fe 14 B magnet (shown in "Methods") 25 .
To calculate the free energy as a function of the z-component of the magnetization F(M z ) for a huge system (up to 24.6 nm × 24.6 nm × 25.6 nm, 1,130,626 spins), the size required to handle a nucleation process, we have developed an approach for the coercivity calculation based on the replica-exchange Wang-Landau (REWL) method 26,27 . Remarkably, we can perform an efficient parallelization scheme for the original Wang-Landau Monte Carlo (WLMC) method 28 (see "Methods").

RESULTS AND DISCUSSION
Temperature dependence of the coercivity The obtained F(M z ) for zero magnetic field (H z = 0) at T ¼ 0:46 T cal C is depicted in Fig. 1a by the red curve (T cal C is the Curie temperature of the model). Notably, this is the quantitatively correct (not schematic) form of F(M z ) for an atomistic spin model 1 representing the Nd 2 Fe 14 B magnet. Shapes of the free energy with the reverse magnetic fields H z are exactly given by F(M z ) + μ 0 H z M z . This is clearly illustrated in Fig. 1a. Close to the spinodal point, the point at which the metastability disappears, magnetization reversal occurs when the free energy barrier F B (see Fig. 1a) becomes compatible with the temperature. The H z dependence of F B is given in Fig. 1b, where we define H 0 as the field at which F B = 0. Similarly, H c is defined as the field at which F B = 25.3 k B T. We call the former "spinodal field" and the latter "coercivity with thermal activation". As mentioned above, H c corresponds to the coercive field for the observation time τ = 1.0 s. In order to see how the magnetization reversal initiates, in Fig. 1c, we depict the spatial distributions of the magnetization reversal probability at the points shown by arrows in Fig. 1a, which were obtained by the REWL method. As expected, we found out that the reversal begins at the corner. There are two reasons for this reversal: One is the decrease of a local magnetic anisotropy due to surface thermal fluctuations. Secondly, the reduction of an exchange energy by decouplings at the corner facilitates the nucleation. Here, all the corners are equivalent, and symmetric configurations are found. However, in an individual process, one of the corners is selected. Figure 1d illustrates a plot of MC snapshots at the corresponding fields.
The temperature dependencies of several formulae for the coercivities are depicted in Fig. 2, and are compared with the ideal coercivity H k . Here, we determined H k as H 0 of the same system size with the periodic boundary condition. In the calculated temperature range, we confirmed that H k takes almost same value as the magnetic anisotropy field H a ¼ 2K 1 =M s , where K 1 is the magnetic anisotropy constant 20,25 and M s is the saturated magnetization at the given temperature. The coercivity H r of permanent magnets is reduced from its limit value H k by various external influences. These reductions have been expressed in the following forms 13,16 : where H t = H 0 − H c and α = H 0 /H k . Here, H t and α represent the thermal activation field and phenomenological factor, respectively. The factor α often means the reduction due to the microstructure such as defects, and the surfaces. On the other hand, α0 ¼ H c =H k is an expression with an alternate definition for H k including H t . The second line of Eq. 2 is known as "Kronmüller equation". The terms of the demagnetization field ÀM s N d express the reduction due to the dipole-dipole (DD) interaction.  The green and purple squared lines denote the experimental measurements in sintered 29 and the hot-deformed with grain boundary diffusion of Nd-Cu alloy 14 magnets, respectively. Inset Our approach can handle the temperature dependence of α and H t explicitly. Unlike the continuum approximation, temperature is naturally introduced in our model by the Boltzmann weights in the MC simulation. This is because we used the atomistic model. The continuum approximation works with temperaturedependent input parameters in which thermal fluctuations are not included in the motion. We have also plotted the results of coercivity H r with respect to the demagnetization factors N d = 0.5 and 1.0 with temperature-dependent magnetization M s . Since the demagnetization field ÀN d M s approximately introduces the effects of the DD interaction as the uniform field, so it cannot consider the size and shape dependences of N d 30 . However, the simulation results qualitatively validate those of the experiments because an increase in temperature leads to a downward convexity in coercivity even without ÀN d M s . The effects of the DD interaction will be discussed later. In the figure, we compare the result with experimental observations. The green squares are data for a sintered sample and the purple ones for a diffused hotdeformed sample. In the former, the effect of polycrystalline causes a further reduction 31 , while in the latter grains are separated and compatible with the present calculation. In the former case, other effects would subject to the multi-grain approach in the future. Thus, the present result gives the upper limit of the coercivity at finite temperatures.
For some sintered polycrystalline magnets (corresponding to the green line in Fig. 2), Kronmüller and Durst phenomenologically estimated the decay factor as α 0 exp ¼ 0:89 À 0:93 from measured magnetic properties around room temperature 16 . In the inset of Fig. 2, we show the temperature dependence of α and α 0 . The difference of α from 1 reflects the decrease of the surface magnetic anisotropy by thermal fluctuations, and α 0 also includes the thermal activation effects. The decay factors at room temperature (here, T = 0.51 T C ) are determined as α = 0.93 and α 0 ¼ 0:73. In the case of the experiment, exchange interactions at the grain boundary interfaces suppress the interface thermal fluctuations. The suppression is one of the origins of the difference between α 0 exp and α 0 .
Nucleation mechanism Next, we consider the mechanism of the thermal activation (nucleation) process for which the concept of "activation volume" has been introduced [11][12][13][14][15] . The symbol ΔM z represents the difference of the magnetization between those at the local minimum M z = M b of the free energy and those at the local maximum M z = M t (see Fig. 1a). It is confirmed that ΔM z gives the activation volume V: Note that the unit of M s is (μ B /volume). Activation volume has been defined by By differentiating F B : it is obvious that definition Eqs. 3 and 4 are equivalent. Here, we note that F(M t(b) , 0) depends on H z since M t(b) is a function of H z . In Fig. 3a, the activation volumes obtained by the both definitions Eqs. 3 and 4 were plotted. This confirms the equivalence. From Eq. 4, the field dependence of F B is related to the nucleation process. Thus, we focus on the exponent n which is widely used for the dependence in the phenomenological form: The exponent n is 2 for coherent magnetization reversal. However, according to the authors 12,14 , the experimental value for n was given to be approximately 1.0. We evaluated the value of n in the atomistic magnetization reversal process around the coercivity. It should be noted that the value of n varies depending on the range of μ 0 H z (shown in Fig. 3a), where we fit the data. However, for a wide range of μ 0 H z , it was observed that the value of n was approximately 1.3. We, therefore, propose that the small value of n is ascribed to the peculiar shape of the free energy. Here we define δF(δM z ) for the free energy near the local minimum, i.e., δF The shapes of δF(δM z ) for the fields μ 0 H z = 0, 2.8, and 3.61 T are depicted in Fig. 3b, where a cusp is pointed by the arrow in each figure.
The cusp represents the point at which the type of the magnetization reversal changes 32 : for small value of δM z (see in Fig. 3b), the reversal process is the nucleation type, where we expect while for large value of δM z , δF(δM z ) is mainly given by the change due to a domain wall (DW) motion, where we expect This magnetization reversal mechanism is schematically pictured in Fig. 4. This difference can be understood in the following picture. In the case of the Heisenberg type with uniaxial anisotropy (corresponding to the Nd 2 Fe 14 B spin model), the critical nucleation size (ΔM 1=3 z R c ) is about the domain wall width δ dw . And thus, before the cusp point (R < R c , where R is the width of the magnetization reversing region, see Fig. 4b) the nucleation region is growing where we have Eq. 7. Around the cusp point, a fully reversed region appears. And after this point, the excess free energy δF comes from the domain wall, i.e., the surface of the fully reversed region. In this case we have Eq. 8. For comparison, we show the case of Ising model where δ dw is much smaller than R c at low temperatures, and we have the dependence Eq. 8 for almost all R except very small R. This comparison is schematically pictured in Fig. 4d.
To clearly see the change of exponent, in Fig. 3b, we show the differentiate log ðδFÞ by log ðδM z Þ with the green curve which gives roughly the power k when we assume the dependence: δF ¼ CδM k z þ D (C and D are constant). We find that k is close to 2/ 3 for large δM z and jump at the cusp. Around δM z = 0, thermal fluctuations enhance the behavior of k = 2, i.e., δFðδM z Þ / δM 2 z ; which makes it difficult to identify the behavior δF ∝ δM z (i.e., k = 1) for the nucleation process. This fact causes that we have an intermediate value (n ≃ 1.3) numerically. However, k clearly changes from 2 to 1 with δM z , which supports the present scenario.
In any case, δF has regions of upward convex and downward convex between which the cusp exists. This structure causes ΔM z to be less dependent on H z for a range of H z 33 , which can be seen from comparing the blue points in Fig. 4a, c. In this range, V is a constant and it results in the value of n = 1 by Eq. 4. Note that, the reason for the cusp explains the drastic change in the magnetization distribution between (i) and (iii) in Fig. 1c, d.
If we assume n = 1 for any H z , then the result is the widely used phenomenological equation for the thermal activation effects 13 : where V c is the value of V at H z = H c . In Fig. 5, we compare the temperature dependence of H t and H 0 t . We found a qualitatively similar dependence. The difference between H 0 t and H t becomes significant in the high-temperature range, which is attributed to the fact that the activation volume and n are not exactly constant. Thus, the difference due to the changes in the magnetization reversal type and surface effects, as shown in Figs. 3 and 4.
It is necessary to mention the exact effects of the DD interaction, which is approximately introduced as the demagnetization field in the present study. The DD interaction brings the non-uniformity of local demagnetization fields near corners and edges in grains, unlike the demagnetization field. The non-uniformity should change the nucleation process from the calculation results. However, the strength of the local field increases logarithmically with system size 30,34,35 , and the nucleation behavior is protected by the cusp structure. Thus, the nucleation mechanism we proposed is expected not to qualitatively change up to a certain large grain size (L x ≃ 1 − 10 μm). Above this size, since magnetic domain structures are formed 36 , the DD interaction should be incorporated into the spin Hamiltonian.
In summary, our atomistic approach by the free energy landscape simulations revealed essential characteristics of the thermal activation effects for permanent magnets, i.e., the downward convexity in the coercivity concerning the temperature, the microscopic definition of the activation (nucleation) volume, and also the exponent n = 1 for the well-used energy barrier formula (6). The present study should be the first step to study the characterization of coercivity at finite temperatures from microscopic information, and the method will be extended to the cases complex of grains with the DD interaction in the future.

Atomistic spin model
The atomistic classical spin Hamiltonian of Nd 2 Fe 14 B magnet was recently proposed for the study of thermodynamic properties at finite temperatures by using Monte Carlo (MC) method 20,25,35,37 and stochastic Landau-Lifshitz-Gilbert equation 38 . It is given by where J ij is the exchange coupling constants between the i-th and j-th spins up to a cut-off range (in the present work 3.52 Å), and s i is the normalized spin moment at the i-th site. The second and third terms are the magnetic anisotropy of Fe and Nd sites, respectively. The coefficients J ij , D i , D l i and magnetic moments are determined from ab-initio calculations and experimental data 25 . Note that Curie temperature of the model with given parameters is T cal C ¼ 870 K, which is higher than the experimental observation (T exp C ' 585 K). This difference must be adjusted by finetuning of the parameters, but the qualitative properties are well reproduced 20,25 .

Wang-Landau method for magnetization
The partition function of a spin system is generally described as follows: where β = 1/k B T is the inverse temperature, E is the internal energy of a spin configuration R, M z is the z-component of a total magnetization and g  Fig. 4 Nucleation picture in permanent magnets. Schematic pictures of a, c energy landscape and b, d spin configuration in the magnetization reversal processes at low-temperature limit. (E, M z ) is the joint density of states. By using g(E, M z ), the magnetization dependence of the free energy F(β, M z ) and a partial partition function Z Mz can be defined by following formula: The Wang-Landau Monte Carlo (WLMC) method 28  As mentioned above, adjusting the weight for two-dimensional or multidimensional space requires massive calculation cost. To reduce the cost, for the MC sampling in the energy space, we use Boltzmann weight, i.e., wðE; M z Þ ¼ wðM z Þ expðÀβEÞ 39,40 . Then, the histogram of M z space can be written as follows (using Eq. 14): Therefore, by obtaining a flat histogram in one dimensional space of M z to adjust w(M z ), the free energy can be calculated from Eq. 13: where C is constant, w f is the adjusted w when obtaining the flat histogram. Since the continuous spin system cannot determine C by calculation at each temperature, so in this study, we proceed with the calculation of the coercivity using the relative values of F. Although the above scheme requires the simulation for each temperature due to the use of Boltzmann weight, it significantly reduces the calculation cost for free energy as a function of M z .

1/t algorithm
To obtain a flat histogram, the WLMC method adjusts the weight as wðM 0 z Þ ! e Àη wðM 0 z Þ every update attempt (even if no update) of the spin state with a random walk, where η(>0) is a modification factor, M 0 z is the total magnetization of the state after the attempt. From the formulations of the previous section, update attempts of the spin state were performed according to the following transition probability: The REWL method for M z with 1/t algorithm.

Replica exchange parallelization approach
For further large-scale simulation, we use the replica-exchange Wang-Landau (REWL) method 26,27 , which is an efficient parallelization approach. The simple idea for parallelization in the WLMC method can be realized to divide the magnetization range by a small range and to allocate the piece into each processor. Each processor simulates the WLMC method individually. However, the small range may break ergodicity and inhibit the relaxation of the spin state. The REWL method is to recover the ergodicity by exchanging the spin configurations among these processors.
For the replica-exchange (RE), it is necessary to overlap the areas allocated to each processor, like Fig. 6. Previous research 27 proposed that the efficiently overlap ratio is 62.5-75% (we set 70%). In this study, every N s spin update attempt, we try to exchange the spin configuration between neighbor processors according to the following exchange probability: w n ðYÞw m ðXÞ w n ðXÞw m ðYÞ where n(, m) is the index of the adjacent processors, X(, Y) is the spin configuration from which E and M z can be calculated, and w n (, G n ) is the values of w(, G) on the processor n. The exchange probability only depends on M z and not on E. If the spin configurations are not in the overlapping range, P R = 0. As mentioned in Eq. 16, F (and G) can only be calculated as relative values. Thus, the values of F for the processors require to be corrected and connected so that the average value of the overlapping range is equal. After this connection, the average of F for all the processors is calculated as a result. Note that, just to improve connectivity among the processors, we Y. Toga et al.
omit the edge ranges (5%, shown in Fig. 6) of each processor in the above correction and averaging (not so important).
To show the effectiveness of the RE, in Fig. 7a, b, we plot the flatness δh ¼ maxðhÞ À minðhÞ ½ =meanðhÞ as a function of t=N bin ð¼ 1=ηÞ in the simulations with and without RE. The accuracy of the REWL method is determined by the smallness of δH andη. Thus, the figure clearly indicates that RE is necessary to obtain a converged result, especially at large system size. The results in this study are set to satisfy the following convergence conditions: δH < 10 −2 andη < 10 À8 ðeVÞ. Figure 7c and d show the parallelization efficiency for the two system sizes in strong scaling with under the two convergence conditions. This result indicates that the parallelization is achieved with the ideal efficiency.
The pseudocode of the implemented methods in the present study is described in Algorithm 1. Through the methods, for large system size (up to 1,130,626 spins), we made it possible to calculate the free energy as a function of M z .

DATA AVAILABILITY
The data reported in this paper is available from the corresponding author upon reasonable request.  Fig. 6 Diagram of parallelization for the REWL method. The calculating range of the magnetization M z is partitioned by many processors, here P 1−4 . The double-headed arrows mean the replicaexchange and the shaded ranges denote the area to be omitted when connecting G among the processors (see text).