Role of atomic-scale thermal fluctuations in the coercivity

Abstract

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 Nd₂Fe₁₄B. 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₂Fe₁₄B⁸. 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.¹⁶. In this regard, extensive works have been done to elucidate the reason for the discrepancy^{17,18,19,20,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)²². 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₂Fe₁₄B²⁴.

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):

$$\tau ={\tau }_{0}{{\rm{e}}}^{\beta {F}_{{\rm{B}}}},$$

(1)

where τ₀ is the reciprocal of the attempt frequency, which is typically set to 10⁻¹¹ s¹². Then, for the observation time of 1 s, setting τ = 1 s, the barrier height is given by F_B = 25.3 k_BT. The free energy is given by the Nd₂Fe₁₄B atomistic spin Hamiltonian, which was recently developed and successfully reproduced static properties of the Nd₂Fe₁₄B magnet (shown in “Methods”)²⁵.

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²⁸ (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}_{{\rm{C}}}^{{\rm{cal}}}$ is depicted in Fig. 1a by the red curve (${T}_{{\rm{C}}}^{{\rm{cal}}}$ 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 representing the Nd₂Fe₁₄B magnet. Shapes of the free energy with the reverse magnetic fields H_z are exactly given by F(M_z) + μ₀H_zM_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₀ as the field at which F_B = 0. Similarly, H_c is defined as the field at which F_B = 25.3 k_BT. 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.

Fig. 1: Free energy landscape simulation.

a Free energies as a function of M_z for the Nd₂Fe₁₄B isolated grains at $0.46\,{T}_{{\rm{C}}}^{{\rm{cal}}}$, whose size is (L_x, L_y, L_z) = (14.1, 14.1, 14.6) nm (212,536 spins). Red line is for H_z = 0 and other lines are for applying H_z. b Free-energy barriers as a function of μ₀H_z for four system sizes: L_x = 10.6, 14.1, 21.1, and 24.6 nm (L_y = L_x, L_z = 1.038L_x), evaluated from F(M_z, H_z) similar to Fig. 1a. c The distributions of reversed probability P_d at Fe spins sliced by (110) plane, which correspond to the points (i–iii) in Fig. 1a. d Four MC snapshots of reversed Fe spins for the three points.

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₀ 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}=2{K}_{1}/{{\mathcal{M}}}_{\mathrm{s}}$, where K₁ is the magnetic anisotropy constant^20,25 and ${{\mathcal{M}}}_{\mathrm{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:

$$\begin{array}{rcl}{H}_{r}&=&\alpha {H}_{\mathrm{k}}-{H}_{\mathrm{t}}-{{\mathcal{M}}}_{\mathrm{s}}{N}_{\mathrm{d}}\\ &=&{\alpha }^{\prime}{H}_{\mathrm{k}}-{{\mathcal{M}}}_{\mathrm{s}}{N}_{\mathrm{d}},\end{array}$$

(2)

where H_t = H₀ − H_c and α = H₀/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, ${\alpha }{\prime}={H}_{\mathrm{c}}/{H}_{\mathrm{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 $-{{\mathcal{M}}}_{\mathrm{s}}{N}_{\mathrm{d}}$ express the reduction due to the dipole–dipole (DD) interaction. 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 temperature-dependent 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 ${{\mathcal{M}}}_{\mathrm{s}}$. Since the demagnetization field $-{N}_{\mathrm{d}}{{\mathcal{M}}}_{\mathrm{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³⁰. 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}_{\mathrm{d}}{{\mathcal{M}}}_{\mathrm{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 hot-deformed sample. In the former, the effect of polycrystalline causes a further reduction³¹, 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.

Fig. 2: Temperature dependence of coercivity.

The blue line μ₀H₀ and the red line μ₀H_c were calculated from Fig. 1b for 21.1 nm × 21.1 nm × 21.9 nm (713,172 spins) isolated grain at each temperature. The colored area depicts the coercivity μ₀H_c under the demagnetization fields in the range of demagnetization factor N_d = 0.5–1.0. The green and purple squared lines denote the experimental measurements in sintered²⁹ and the hot-deformed with grain boundary diffusion of Nd-Cu alloy¹⁴ magnets, respectively. Inset shows α = H₀/H_k and ${\alpha }^{\prime}={H}_{\mathrm{c}}/{H}_{\mathrm{k}}$.

For some sintered polycrystalline magnets (corresponding to the green line in Fig. 2), Kronmüller and Durst phenomenologically estimated the decay factor as ${\alpha }_{\exp }^{\prime}=0.89-0.93$ from measured magnetic properties around room temperature¹⁶. In the inset of Fig. 2, we show the temperature dependence of α and ${\alpha }^{\prime}$. The difference of α from 1 reflects the decrease of the surface magnetic anisotropy by thermal fluctuations, and ${\alpha }^{\prime}$ also includes the thermal activation effects. The decay factors at room temperature (here, T = 0.51 T_C) are determined as α = 0.93 and ${\alpha }^{\prime}=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 ${\alpha }_{\exp }^{\prime}$ and ${\alpha }^{\prime}$.

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:

$$V=\Delta {M}_{\mathrm{z}}/{{\mathcal{M}}}_{\mathrm{s}}.$$

(3)

Note that the unit of ${{\mathcal{M}}}_{\mathrm{s}}$ is (μ_B/volume). Activation volume has been defined by

$$V=-\frac{1}{{\mu }_{0}{{\mathcal{M}}}_{\mathrm{s}}}\frac{\partial {F}_{{\rm{B}}}}{\partial {H}_{\mathrm{z}}}.$$

(4)

By differentiating F_B:

$$\begin{array}{rcl}{F}_{{\rm{B}}}({H}_{\mathrm{z}})&=&F({M}_{\mathrm{t}},{H}_{\mathrm{z}})-F({M}_{\mathrm{b}},{H}_{\mathrm{z}})\\ &=&F({M}_{\mathrm{t}},0)-F({M}_{\mathrm{b}},0)+{\mu }_{0}{H}_{\mathrm{z}}({M}_{\mathrm{t}}-{M}_{\mathrm{b}}),\end{array}$$

(5)

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.

Fig. 3: Magnetic field dependence of activation volume.

a The activation volume V as a function of μ₀H_z for two system sizes at $0.46\,{T}_{{\rm{C}}}^{{\rm{cal}}}$. Circle points are evaluated from Eq. 3 while solid lines are evaluated from Eq. 4. Squared lines are the exponent n of free energy barrier in formula 6. b Field dependence of the shape of free energy as a function of δM_z for L_x = 24.6 nm. Blue points are the top of the free energy barrier, which corresponds to M_t, and black arrows denote the cusp in the free energy. Green line shows the value on the right axis.

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:

$${F}_{{\rm{B}}}\propto {(1-{H}_{\mathrm{z}}/{H}_{0})}^{n}.$$

(6)

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 μ₀H_z (shown in Fig. 3a), where we fit the data. However, for a wide range of μ₀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(δM_z) = F(M_b − δM_z) − F(M_b), where δF(ΔM_z) = F_B. The shapes of δF(δM_z) for the fields μ₀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³²: for small value of δM_z (see in Fig. 3b), the reversal process is the nucleation type, where we expect

$$\delta F(\delta {M}_{\mathrm{z}})\propto \delta {M}_{\mathrm{z}},$$

(7)

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

$$\delta F(\delta {M}_{\mathrm{z}})\propto \delta {M}_{\mathrm{z}}^{2/3}+{\rm{const.}}$$

(8)

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₂Fe₁₄B spin model), the critical nucleation size ($\Delta {M}_{\mathrm{z}}^{1/3}\equiv {R}_{\mathrm{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.

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. a, b Heisenberg type spin model with uniaxial magnetic anisotropy. c, d Ising type. The energy landscapes under several magnetic fields are drawn (the lowest line is for H_z = 0). Dotted lines are the reverse magnetic field terms μ₀H_zM_z and blue points are the top of the energy barriers.

To clearly see the change of exponent, in Fig. 3b, we show the differentiate $\mathrm{log}\,(\delta F)$ by $\mathrm{log}\,(\delta {M}_{\mathrm{z}})$ with the green curve which gives roughly the power k when we assume the dependence: $\delta F=C\delta {M}_{\mathrm{z}}^{k}+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.,

$$\delta F(\delta {M}_{\mathrm{z}})\propto \delta {M}_{\mathrm{z}}^{2},$$

(9)

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³³, 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¹³:

$${H}_{\mathrm{t}}^{\prime}=\frac{25.3{k}_{{\rm{B}}}T}{{\mu }_{0}{{\mathcal{M}}}_{\mathrm{s}}{V}_{\mathrm{c}}},$$

(10)

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}_{\mathrm{t}}^{\prime}$. We found a qualitatively similar dependence. The difference between ${H}_{\mathrm{t}}^{\prime}$ 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.

Fig. 5: Thermal activation reduction.

The shifts of the coercivity due to the thermal activation effects, evaluated from the two ways: H_t = H₀ − H_c and Eq. 10.

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³⁶, 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.

Methods

Atomistic spin model

The atomistic classical spin Hamiltonian of Nd₂Fe₁₄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³⁸. It is given by

$${\mathcal{H}}=-2\mathop{\sum }\limits_{\mathrm{i}\,{<}\,\mathrm{j}}{J}_{\mathrm{ij}}{{\boldsymbol{s}}}_{\mathrm{i}}\cdot {{\boldsymbol{s}}}_{\mathrm{j}}-\mathop{\sum }\limits_{\mathrm{i}\in {\rm{Fe}}}{D}_{\mathrm{i}}{({s}_{\mathrm{i}}^{z})}^{2}-\mathop{\sum }\limits_{\mathrm{i}\in {\rm{Nd}}}\mathop{\sum }\limits_{\mathrm{l} = 2,4,6}{D}_{\mathrm{i}}^{\mathrm{l}}{({s}_{\mathrm{i}}^{z})}^{\mathrm{l}},$$

(11)

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}_{\mathrm{i}}^{\mathrm{l}}$ and magnetic moments are determined from ab-initio calculations and experimental data²⁵. Note that Curie temperature of the model with given parameters is ${T}_{{\rm{C}}}^{{\rm{cal}}}=870\ {\rm{K}}$, which is higher than the experimental observation (${T}_{{\rm{C}}}^{\exp }\simeq 585\ {\rm{K}}$). This difference must be adjusted by fine-tuning 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:

$$Z=\mathop{\sum }\limits_{{\boldsymbol{R}}}{{\rm{e}}}^{-\beta H({\boldsymbol{R}})}=\mathop{\sum }\limits_{E}\mathop{\sum }\limits_{{M}_{\mathrm{z}}}g(E,{M}_{\mathrm{z}}){{\rm{e}}}^{-\beta E},$$

(12)

where β = 1/k_BT 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(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}_{{M}_{\mathrm{z}}}$ can be defined by following formula:

$$F(\beta ,{M}_{\mathrm{z}})=-{\beta }^{-1}{\mathrm{ln}}\,{Z}_{{M}_{\mathrm{z}}},$$

(13)

$${Z}_{{M}_{\mathrm{z}}}=\mathop{\sum }\limits_{E}g(E,{M}_{\mathrm{z}}){{\rm{e}}}^{-\beta E}.$$

(14)

The Wang–Landau Monte Carlo (WLMC) method²⁸ allows us to calculate g(E, M_z) in the basis of a relation: h(E, M_z) = g(E, M_z)w(E, M_z), where h(E, M_z) is the histogram of the spin states obtained by MC sampling with a weight w(E, M_z). If a flat histogram (h is constant) can be obtained by arbitrarily adjusting w(E, M_z) in a MC simulation, g(E, M_z) is proportional to the reciprocal of the adjusted w(E, M_z).

As mentioned above, adjusting the weight for two-dimensional or multi-dimensional 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}_{\mathrm{z}})=w({M}_{\mathrm{z}})\exp (-\beta E)$^39,40. Then, the histogram of M_z space can be written as follows (using Eq. 14):

$$\begin{array}{rcl}h({M}_{\mathrm{z}})=\mathop{\sum }\limits_{E}h(E,{M}_{\mathrm{z}})\\ =w({M}_{\mathrm{z}}){Z}_{{M}_{\mathrm{z}}}.\end{array}$$

(15)

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:

$$\begin{array}{rcl}F(\beta ,{M}_{\mathrm{z}})=-{\beta }^{-1}{\mathrm{ln}}\,\frac{C}{{w}_{\mathrm{f}}({M}_{\mathrm{z}})}\\ =C+{\beta }^{-1}{\mathrm{ln}}\,{w}_{\mathrm{f}}({M}_{\mathrm{z}}),\end{array}$$

(16)

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}_{\mathrm{z}}^{\prime})\to {{\rm{e}}}^{-\eta }w({M}_{\mathrm{z}}^{\prime})$ every update attempt (even if no update) of the spin state with a random walk, where η(>0) is a modification factor, ${M}_{\mathrm{z}}^{\prime}$ 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:

$$P(E,{M}_{\mathrm{z}}\to {E}^{^{\prime\prime} },{M}_{\mathrm{z}}^{^{\prime\prime} }) =\min \left(1,\frac{w({E}^{^{\prime\prime} },{M}_{\mathrm{z}}^{^{\prime\prime} })}{w(E,{M}_{\mathrm{z}})}\right)$$

(17)

$$=\min \left(1,{{\rm{e}}}^{-\beta \left[{E}^{^{\prime\prime} }-E-G({M}_{\mathrm{z}}^{^{\prime\prime} })+G({M}_{\mathrm{z}})\right]}\right),$$

$$G({M}_{\mathrm{z}})={\beta }^{-1}{\mathrm{ln}}\,w({M}_{\mathrm{z}}),$$

where E″ and ${M}_{\mathrm{z}}^{^{\prime\prime} }$ are the values of a trial state. Therefore, in the actual simulation, G is adjusted instead of w as $G({M}_{\mathrm{z}}^{\prime})\to G({M}_{\mathrm{z}}^{\prime})-\tilde{\eta }$, where $\tilde{\eta }={\beta }^{-1}\eta$ (>0). This redefined modification factor $\tilde{\eta }$ has the unit of energy, which is convenient to adjust G with consistent values regardless of temperature.

Adjustment of G at every update attempt accelerates the convergence; however, it breaks a detailed balance condition. Thus, it is necessary to reduce $\tilde{\eta }$ so as not to affect the results of the MC simulation. To perform the reduction, the present study adopt the so-called 1/t algorithm⁴¹ which resolves systematic error due to histogram flatness condition in the original WLMC method. The algorithm reduces the $\tilde{\eta }$ as αN_bin/t, where α is proportionality coefficient (here, we set α = 1.0 eV), t is the number of update attempts of the spin state and N_bin is the number of M_z levels in the simulation range of interest, ${M}_{\mathrm{z}}\in [{M}_{\mathrm{z}}^{\min },{M}_{\mathrm{z}}^{\max }]$. Note that, since the WLMC method requires to treat M_z discretely for the MC sampling with Eq. 17, in this study, we adopt the grid where h(M_z) and G(M_z) are discretized with the bins of m_w = 0.3 − 0.5μ_B width, and ${N}_{{\rm{bin}}}=({M}_{\mathrm{z}}^{\max }-{M}_{\mathrm{z}}^{\min })/{m}_{w}$. As a practical matter, during small t, the changes in $\tilde{\eta }$ is so large that the simulation does not convergence. Thus, the initial value of $\tilde{\eta }$ is set to be 1.0 eV, and it is reduced until $\tilde{\eta }\le {N}_{{\rm{bin}}}/t$ (eV) as $\tilde{\eta }\to \tilde{\eta }/2$ at every random walker visit to the all the magnetization bins.

Algorithm 1.

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²⁷ 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:

$${P}_{R}(X\in n\leftrightarrow Y\in m) =\min \left(1,\frac{{w}_{\mathrm{n}}(Y){w}_{\mathrm{m}}(X)}{{w}_{\mathrm{n}}(X){w}_{\mathrm{m}}(Y)}\right)$$

(18)

$$=\min \left(1,{{\rm{e}}}^{-\beta \left[{G}_{\mathrm{n}}(X)+{G}_{\mathrm{m}}(Y)-{G}_{\mathrm{n}}(Y)-{G}_{\mathrm{m}}(X)\right]}\right),$$

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.

Fig. 6: Diagram of parallelization for the REWL method.

The calculating range of the magnetization M_z is partitioned by many processors, here P₁₋₄. The double-headed arrows mean the replica-exchange and the shaded ranges denote the area to be omitted when connecting G among the processors (see text).

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 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 $\delta h=\left[\max (h)-\min (h)\right]/{\rm{mean}}(h)$ as a function of $t/{N}_{{\rm{bin}}}(=1/\tilde{\eta })$ in the simulations with and without RE. The accuracy of the REWL method is determined by the smallness of δH and $\tilde{\eta }$. 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⁻² and $\tilde{\eta }\;<\;1{0}^{-8}\ ({\rm{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.

Fig. 7: Parallelization efficiency of replica exchange.

The flatness δH as a function of $t/{N}_{{\rm{bin}}}(=1/\tilde{\eta }\,(1/{\rm{eV}}))$ at $0.46\ {T}_{{\rm{C}}}^{{\rm{cal}}}$ of the Nd₂Fe₁₄B spin model with the isolated grain for two sizes a L_x = 10.6 nm and b L_x = 14.1 nm. Black (Blue) lines are the results without (with) RE, which are performed by 864 processors. In the simulations, the calculated ranges of the magnetization are ${M}_{\mathrm{z}}\in [-0.02{M}_{\mathrm{z}},\ 0.76{M}_{\mathrm{z}}^{\mathrm{t}}]$ for L_x = 10.6 nm and ${M}_{\mathrm{z}}\in [0.3{M}_{\mathrm{z}}^{\mathrm{t}},\ 0.76{M}_{\mathrm{z}}^{\mathrm{t}}]$ for L_x = 14.1 nm with m_w = 0.5μ_B, where ${M}_{\mathrm{z}}^{\mathrm{t}}$ is the maximum value of M_z in each spin system. c t_c and d ${t}_{\mathrm{c}}^{\prime}$ are the computational times as a function of the number of processors N_proc, which satisfy the convergence conditions δh < 10⁻² and t/N_nbin > 10⁸, respectively. Arrows in a and b denote the points corresponding to t_c and ${t}_{\mathrm{c}}^{\prime}$.

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.