Atomistic theory of thermally activated magnetization processes in Nd2Fe14B permanent magnet

ABSTRACT To study the temperature dependence of magnetic properties of permanent magnets, methods of treating the thermal fluctuation causing the thermal activation phenomena must be established. To study finite-temperature properties quantitatively, we need atomistic energy information to calculate the canonical distribution. In the present review, we report our recent studies on the thermal properties of the Nd2Fe14B magnet and the methods of studying them. We first propose an atomistic Hamiltonian and show various thermodynamic properties, for example, the temperature dependences of the magnetization showing a spin reorientation transition, the magnetic anisotropy energy, the domain wall profiles, the anisotropy of the exchange stiffness constant, and the spectrum of ferromagnetic resonance. The effects of the dipole–dipole interaction (DDI) in large grains are also presented. In addition to these equilibrium properties, the temperature dependence of the coercivity of a single grain was studied using the stochastic Landau-Lifshitz-Gilbert equation and also by the analysis of the free energy landscape, which was obtained by Monte Carlo simulation. The upper limit of coercivity at room temperature was found to be about 3 T at room temperature. The coercivity of a polycrystalline magnet, that is, an ensemble of interactinve grains, is expected to be reduced further by the effects of the grain boundary phase, which is also studied. Surface nucleation is a key ingredient in the domain wall depinning process. Finally, we study the effect of DDI among grains and also discuss the distribution of properties of grains from the viewpoint of first-order reversal curve.


Introduction
The neodymium (Nd) polycrystalline magnet consisting of Nd 2 Fe 14 B [1][2][3][4][5][6][7][8][9][10][11][12][13][14] is an important highperformance permanent magnet. Because of its high coercivity, it is widely used for electric motors, electronic devices, and so forth [15]. Coercivity at finite temperatures is a key factor affecting the performance of permanent magnets. Coercivity essentially depends on the structure of grains and grain boundaries, and the nucleation of reversed magnetization and the depinning mechanisms of magnetic domain walls in the structure play an important role in coercivity [15][16][17]. Trials toward achieving higher coercivities at higher temperatures have been actively performed [18,19], but the quantitative properties of coercivity at finite temperatures have not been well understood [16].
In previous works, temperature effects have been taken into account by using the temperaturerenormalized parameters, for example the exchange stiffness constant AðtÞ and the magnetic anisotropy energy KðTÞ, which are obtained experimentally or by mean-field analyses. Coercivity at finite temperatures is, however, a phenomenon involving the breakdown of a metastable magnetic state. Magnetization reversal occurs with thermal agitation and thus it is a stochastic process [20]. To study such effects quantitatively, we need to take into account the effect of entropy. For this purpose, we must treat the temperature precisely. In thermal equilibrium state, the probability of a state i, P eq ðiÞ, is given by the canonical ensemble with the Boltzmann factor: where H is the Hamiltonian representing the system. If we use a coarse-grained Hamiltonian of the system, such as the continuum model, the definition of the degrees of freedom of the state is ambiguous and the above-mentioned probability is difficult to apply. Namely, if we study finite-temperature properties by micromagnetic simulation with continuous spin having temperature-dependent parameters, such as AðTÞ and KðTÞ, the additional noise causes the double counting of the effect of thermal agitation. To appropriately take into account the temperature effect, we must use standard statistical mechanical methods with the canonical ensemble. To obtain the Boltzmann factor, one needs an atomistic Hamiltonian, which will be explained in Section 2. Using the atomistic Hamiltonian, the thermal effect is taken into account in the Monte Carlo (MC) method by the detailed balance condition and in the SLLG method by the fluctuation dissipation relation, which guarantees realization of the thermal equilibrium state at given temperatures as explained in Refs. [21,22].
Thus, it is necessary to use an atomistic model, and its Hamiltonian must be kept unchanged with the temperature. One may consider cases in which the lattice of the system changes with the temperature (expansion or shrinkage), which causes changes in parameters in the Hamiltonian. In such cases, we must introduce some additional compromised treatments. However, in the present work, we concentrate on the cases where we can ignore such an effect.
In the present paper, we review our recent works on finite-temperature properties. As mentioned above, an atomistic Hamiltonian for a material is necessary to study its thermal properties. Thus, we constructed an explicit atomistic Hamiltonian. Nd 2 Fe 14 B which is the main phase in the neodymium magnet, has a complex structure. The unit cell contains nine different sites with 68 atoms as depicted in Figure 1(a). Although, in general, the determination of the atomistic Hamiltonian is difficult even with the most sophisticated first-principles calculations, we constructed an atomistic Hamiltonian using the latest knowledge of microscopic parameters to reproduce the known thermodynamic properties as explained in Section 2.
Using the Hamiltonian, we first studied various thermal properties of the Nd 2 Fe 14 B magnet, such as magnetization. As an important characteristic of the magnet, we also studied the temperature dependence of anisotropy energy, as well as the temperature dependence of the threshold field of the magnetization reversal (coercivity) from the temperature dependence of the free energy as a function of the direction of magnetization by constrained MC simulation [23]. The dependence of the exchange stiffness constant on the direction reflecting the anisotropy of the crystal was also studied [24]. In addition to the above-mentioned macroscopic quantities, it has been found that the model can produce the microscopic magnetic structures of the magnet, for example, domain wall profiles, and also their temperature dependence [25]. These results well reproduce the experimental data, and thus we confirmed the validity of the model. Similar works with an atomistic model for the finite-temperature properties of magnetization and the domain wall have been recently reported by Gong and coauthors [26][27][28], who obtained thermodynamic quantities and also an effective continuous model. Using the model, they obtained the dependence of stiffness constants on the direction and the coercivity of a system with a grain boundary phase as mentioned above [24].
Furthermore, the atomistic model can produce the spectrum of ferromagnetic resonance(FMR) by stochastic Landau-Lifshitz-Gilbert (SLLG) simulation [22,29]. The effects of the dipole-dipole interaction (DDI) at room temperature were also studied [30].
Hard-magnetic compounds consist of many grains, and the reversal of magnetization under a reverse magnetic field close to the coercivity occurs as a cascade of reversals of magnetizations of grains. The process has been categorized into two processes. That is, the nucleation of reversal magnetization and the propagation of reversal magnetization across the grain boundary (referred to as the depinning of the domain wall). The overall process of the reversal is very complex, and cannot be studied by atomistic simulation.
Thus, we first study the reversal phenomena in atomistic models in a nanoscale system. These reversal phenomena are governed by the emergence of nuclei of nanometer-order magnetic domains (activation volume [31]) according to the classical Arrhenius-type analysis of the time dependence of magnetization. The single-grain coercivity thus obtained gives a theoretical upper limit for the coercivity at a given temperature, although it should be reduced further by the effect of interactions between grains in bulk magnets. It is widely accepted that coercivity is not an intrinsic property of hard magnets and that nucleation occurs at a defect. The reversal of magnetization propagates to hard-magnet grains, for which the pinning of the domain wall at the boundary phase is important. In this depinning process, the propagation of reversal is also governed by the surface nucleation of each grain of a hard magnet, for which the information from nanosize reversal is important.
Note that in contrast to the above-mentioned thermodynamic quantities, which we can calculate using (1), we do not have any explicit theoretical formula to calculate coercivity. Coercivity is the threshold of the magnetic field of the metastable magnetic state in the magnetic field in the opposite direction. At T ¼ 0 without thermal fluctuation, this threshold is uniquely defined. However, at finite temperatures, relaxation occurs stochastically through nucleation. There, the relaxation time is widely distributed. In this manner, coercivity is a highly non-equilibrium property, which prevents us from calculating it theoretically. Coercivity is phenomenologically defined as the magnetic field at which the relaxation time of magnetization is 1s [32][33][34][35][36][37]. We have approached this problem by the following methods.
(1) We studied the dynamics of relaxation by the SLLG method, which incorporates thermal fluctuations [22]. In this approach, the following difficulties exist, which are general problems in microscopic molecular simulations. First, the relaxation depends on the damping constant α of the LLG equation, which is difficult to know precisely. Moreover, the maximum relaxation time that can be obtained by atomistic calculation is limited. Indeed, the time scale of spin precession in a field of 1 T is of the 10 À 12 s order, and the maximum time of simulation is up to several nanoseconds. These difficulties have prevented us from estimating the relaxation time quantitatively. However, we have been able to overcome these difficulties and obtained a quantitative estimation of the threshold field for a relaxation time of 1s [38].
(2) Alternatively, we approached relaxation phenomena from the viewpoint of the free energy barrier at finite temperatures. There have been some works focusing on the free energy barrier using the minimum energy path method [39], which enabled the energy along a path of evolution of magnetization from the metastable state to the stable state to be obtained, where the energy function contains temperaturedependent parameters. For this method, an explicit form of the free energy as a function of the configuration is necessary. Thus far, energy functionals of magnetization with temperature-dependent parameters have been used for this purpose, which may not appropriately express thermal fluctuations. On the other hand, we have developed a method of studying the metastable situation quantitatively at finite temperatures [40], for which we obtained the free energy as a function of magnetization from the atomistic Hamiltonian at a given temperature by using of the Wang-Landau method [41]. Using the explicit form of FðMÞ, we estimated coercivity with and without the activation effect, and also clarified the concept of the activation volume introduced in the literature [32][33][34]42].
With these approaches, we have estimated coercivity at room temperature around 3 T for a single grain with the shape of a 10,20 nm cube. The estimated threshold (, 3 T) gives the theoretical upper limit, which is very low compared with the zero-temperature coercivity and also smaller than the theoretical estimate of 2KðTÞ=MðTÞ for the uniform rotation of magnetization without thermal activation. Thus, the thermal fluctuation has been found to be one of the important ingredients contributing to the Kronmüller's discrepancy [5], although the existence of magnetic inhomogeneity with reduced magnetocrystalline anisotropy has been assumed for the discrepancy.
We also studied the size dependence of coercivity. At finite temperatures, the thermal fluctuation reduces the threshold field (a type of super-paramagnetism). We have confirmed that such reduction saturates when the size reaches 20 nm, and thus the above estimate is essentially valid up to a grain size of a few hundred nanometers. However, for larger grain sizes, the reduction due to DDI increases. Such dependence is also studied in Section 5. When the grain size further increases, DDI induces the formation of a multidomain magnetic structure in a grain [43]. In such a system, the mechanism of nucleation may be different from that of a nanosize grain. We systematically studied the dependence of magnetic patterns in large flat systems on parameters (anisotropy energy, DDI) and obtained a phase diagram of magnetic patterns [44]. To study the metastability of the uniformly magnetized state in such a large system whose stable state is a multidomain magnetic structure, we compared the phase diagrams obtained by thermal-quench process corresponding to the thermal demagnetization process and by field-quench process corresponding to the remanence process, and found metastability in some regions in the phase diagram. We also found that the nucleation breaking the uniformly magnetized state occurs in the middle of the surface owing to DDI effects, in contrast to the case without DDI, in which it occurs from the corners.
For the coercivity of realistic polycrystalline magnets, namely, an ensemble of grains, the interaction among the grains plays an important role. We have studied the effect of a grain boundary phase consisting of a soft magnet and also the threshold field of domain wall pinning in a sandwich (hard-soft-hard magnet) configuration [45][46][47][48][49][50][51]. We also report effects of modifications of surface anisotropy on the coercivity [52]. The effects of misalignment are also important when studying coercivity. Fujisaki et al. [53] studied the alignment dependence of coercivity using the LLG equation and explored the possibility of the Kondorsky dependence 1= cos θ on the field direction. This problem was studied in detail by Bance et al. [54]. The effects of the attachment of a soft-magnet boundary phase have been studied by micromagnetic simulation with temperature-dependent parameters [24,27].
Moreover, the ensemble effects of grains are also important. An ensemble of independent grains (a hysteron) with a distribution of their coercive fields is modeled by the so-called Preisach model [55]. Each grain has a hysteresis curve depending on the value of the field at which the field is swept back (first-order reversal curve (FORC)) [56]. The distribution of the coercive fields is related to the FORC by a mathematical formula, and the distribution is called the FORC diagram. The interaction among grains causes changes in the diagram, and the FORC diagram is used to classify the nature of magnets [57], and it is even used to classify earthenwares in archaeology, because they contain magnetic particles [58]. To study this effect, extensions of Preisach model have been explored [59,60].
The rest of the paper consists of the following: In Section 2, the model used in the present paper is explained. In Section 3, the thermodynamic properties and methods used are presented. In Section 4, approaches to examining the coercivity of nanoparticles by the SLLG and free energy methods are given. In Section 5, the coercivity of large grains where DDI is relevant is studied. In Section 6, the effects of the grain boundary and surface properties are studied. In Section 7, the coercivity of an ensemble of grains is discussed. In Section 8, the summary and perspective are given. In Appendix A, we briefly explain the exchange couplings and magnetic moments that we estimated for the parameters of the Hamiltonian.

Model
We adopt the following atomistic Hamiltonian for the Nd magnet [23]: where s i denotes the classical spin at the ith site and js i j depends on the type of site. In a unit cell of Nd 2 Fe 14 B (Figure 1(a)), there are two types of Nd site, six types of Fe site, and one type of B site. J ij is the exchange interaction between the ith and jth sites, D i is the magnetic anisotropy constant for Fe atoms, the third term is the crystal electric field (CEF) for the magnetic anisotropy energy of Nd atoms, and h is the external magnetic field. The CEF for a rare-earth atom (ion) is given in the following form: Here, Ô m l is the Stevens operator in the classical manner, for example, Ô 0 2 ¼ 3J 2 z À J 2 . Θ l and A m l are the Stevens factor and the coefficient of the spherical harmonics of the CEF, respectively. r is the radiation diameter of an electron and hr l i is the average over the radial wave function estimated in Ref [61]. In our works, we adopt the terms of l ¼ 2; 4; 6 with only the diagonal operators (m ¼ 0), which give the dominant contribution to the tilt of the spins in the ground state.
For Fe and B atoms, s i denotes the magnetic moment at the ith site, but for Nd atoms, s i is the moment of valence (5d and 6s) electrons. For the Zeeman energy, the magnetic moment of site i is given by S i . For the Fe and B sites, S i ¼ s i . For Nd, however, not only s i but also the magnetic moment J i of the 4 f electron contributes. Here, J i consists of the orbital moment L and the spin moment s 4f , which are antiferromagnetically coupled by the spin-orbit interaction (SOI), and is given by jJ i j ¼ g T Jμ B , in which g T ¼ 8=11 is the Landé g-factor and J ¼ 9=2. s 4f and s i are ferromagnetically coupled by the Hund rule as schematically depicted in Figure 2 [62]. The total moment for each Nd atom is S i ¼ s i þ J i , which is used for the interaction with the magnetic field, while the exchange interaction is given in the form J ij s i � s j as in (2), which is antiferromagnetic between Fe and Nd. We adopted the values of Stevens coefficients A m l estimated from experimental data [63]. For the anisotropy of Fe, we considered only the single-ion anisotropy D i and adopted the values calculated by Miura et al. [64]. B makes little contribution to the magnetic Hamiltonian.
Nd 2 Fe 14 B is an itinerant magnet. Here, we express the magnetic interaction in the form of the Heisenberg model. Thus, the interaction is widely distributed owing to the itinerancy of electrons. We adopted the exchange energies between spins obtained by firstprinciples calculation with the Korringa-Kohn-Rostoker (KKR) Green's function method [65] (Akai-KKR). The interactions between spins are widely distributed as shown in Figure 1(b). In the calculation, we cut the interactions separated by longer than r cut . In most of the calculations, we used the interaction up to 3.52 Å. All data for the magnetic moments s i and exchange couplings J ij s i s j estimated by the Akai-KKR method are given in Supplementary Material 1 with a brief explanation in Appendix A.

Magnetizations
First, we study the temperature dependence of magnetization. It is known that the Nd magnet exhibits a spin-reorientation (SR) transition at T r ,135 K, below which the magnetization is tilted from the c-axis [7][8][9]63]. For the tilt in the ground state, the anisotropy energy must have the minimum point at a non-zero angle from the c-axis. Thus, we check the CEF for Nd. Substituting J z ¼ J cos θ into (3), the anisotropy energy at 0 K for Nd atoms is expressed with diagonal terms in the following form: where K 1 ¼ À 3f 2 B 0 2 À 40f 4 B 0 4 À 168f 6 B 0 6 , and so on, with positive constants f l . With the coefficients A 2 0 ¼ 295:0 K a À 2 0 , A 4 0 ¼ À 12:3 K a À 4 0 , and A 6 0 ¼ À 1:84 K a À 6 0 (a 0 is the Bohr radius) estimated by Yamada et al. [63], the first single-ion anisotropy satisfies K 1 < 0 at T ¼ 0.
Note that although A 2 0 > 0, and thus B 0 2 < 0 (∵ the Stevens factor Θ 2 < 0), the contributions of the other terms (B 0 4 and B 0 6 terms) make K 1 < 0. This potential energy causes a tilted magnetization in the ground state (at zero temperature) as shown in Figure 3, in which θ ' 0:2π gives the minimum.
In Figure 4(a), the temperature dependences of the z and xy components of magnetization are depicted [23]. The SR transition is clearly observed at the transition temperature T r ' 145K, which is close to experimentally estimated values ( ' 135K). T r ¼ 145K does not depend on the choice of the range of interaction, r cut , as explained at the end of Section 2. This fact indicates that at low temperatures around T r , the competition between the anisotropy and short-range strong interactions is relevant.
As to the ferromagnetic phase transition, the critical temperature is given as T c , 750-870 K, depending on r cut . These values are slight overestimates compared with the experimental values [4,7] of T c , 585 K. This difference is attributed to the exchange constants used in the calculation, and we may need to rescale the exchange constants. Nevertheless, the overall properties are semi-quantitatively reproduced and there are no serious problems in studying the effects of the thermal fluctuation. Therefore, we accept the present model. When we compare the results with  experimental ones, the temperature rescaled by the critical temperature should be used.
Owing to the atomistic model, we can observe atom-specific properties individually, for example, the temperature dependence of the magnetizations of Fe and Nd atoms as depicted in Figure 4(b). As shown in the figure, the magnetization of Nd decreases much faster than that of Fe as the temperature increases. We attribute this difference to the interactions. Namely, the exchange coupling between Nd and Fe is small, while Fe spins are strongly coupled with each other (the exchange constants between Nd atoms are negligible).

Angle dependence of free energy
To obtain the temperature dependence of anisotropy energies K i ðTÞ (4), we need the angle dependence of the free energy FðθÞ, for which we adopted the constrained Monte Carlo (C-MC) method. In this method, the total magnetization M ¼ ð M x ; M y ; M z Þ is fixed in a direction. For a fixed angle θ (the angle between the c-axis and the magnetization), we calculated the magnetization torque T defined by and using it, the excess free energy ΔFðθÞ ¼ FðθÞ À Fð0Þ is given by where e i and n are the unit vectors of s i and M, respectively [66]. In Figure 5(a), the obtained angle dependence of the torque and the free energy are given. By analyzing the angle dependence of the free energy using the form in (4), the temperature dependences of the coefficients K 1 ðTÞ; K 2 ðTÞ, and K 4 ðTÞ are obtained ( Figure 5(b)), which agree with those in previous works [67,68].
Using the angle dependence of the free energy on the magnetic field H ext , we obtain free energies as a function H ext at various temperatures as depicted in Figure 5(c), which gives an idea for the temperature dependence of the energy barrier for the uniform (coherent) rotation of magnetization [23,69]. However, a more comprehensive study on the temperature dependence of coercivity will be given in Section 4.2.

Domain wall
The profile of a magnetic domain wall is one of the important characteristics of magnets. We studied the profile of the domain wall of the Nd magnet and showed how well the present atomistic model reproduces the microscopic ordering property [25,70]. Because of the anisotropy of the crystal structure, the profile depends on the direction. Along the a-axis (perpendicular to the easy axis), the domain wall is of the Bloch type, while it is of the Néel type along the c-axis, as depicted in Figure 6(a). In the case of the Bloch type propagating along the x-axis, the z component of magnetization shows a step-like structure and the y component shows a bell-type structure, and respectively, and m x ðxÞ ¼ 0. Here, δ 0 is the wall parameter and the width of the domain wall is given by � ¼ πδ 0 . These dependences also apply to the Néel type. The profiles in both cases at T ¼ 300K are depicted in Figure 6(b). The profiles are well fitted by (7) and (8). The width of the domain wall is consistent with the experimentally estimated value [71]. We found that the width of the Bloch type (along the a-axis) is larger than that of the Néel type (along the c-axis). This difference originates from the differences in the stiffness constants with the direction, which we study in more detail in Section 3.3. We obtained the width for various temperatures and found that it increases with the temperature, indicating that the renormalization of the stiffness constant is faster than that of the anisotropy energy.

Anisotropy of exchange stiffness constant
Because of the anisotropy of the crystal, the exchange stiffness depends on the direction as depicted in Figure 7(a), which causes the difference in the domain wall shape as mentioned in Section 3.2 [25]. We explicitly studied the dependence of the stiffness constants A x ðTÞ along the a-axis and A z ðTÞ along the c-axis [24].
To obtain the stiffness constants at a given temperature, we used domain wall energy E dw ðTÞ and the anisotropy energy E K ðT; θÞ [72,73]. The former is expressed in the continuum model as We regard E dw ðTÞ as equal to the domain-wall free energy in the atomistic spin model, F dw ðTÞ, which is given by where E dw ðTÞ is the internal energy of domain-wall formation, which is defined as the difference between internal energies for systems with a parallel periodic condition and an antiparallel periodic boundary condition. The anisotropy energy E K ðT; θÞ was calculated by the C-MC as in Section 3.1.2. The thus obtained temperature dependence of the stiffness constants along the a-axis (A x ðTÞ) and c-axis (A z ðTÞ) is depicted in Figure 7(b). Gong et al. have obtained a similar direction dependence of the stiffness constants and extended the study to the interface exchange coupling strength between the Nd magnet and the grain boundary phase [28].
Recently, the anisotropy of the stiffness constants has been evaluated by spin-wave measurement in a single crystal sample [74], where the anisotropy is weaker than the above result at high temperatures. This difference may be due to how we define the exchange stiffness constant and/or the existence of long-range exchange interactions and so on. This difference must be clarified in the future.
The method used in our approaches is general and can be applied for other systems. For example, similar calculations have also been carried out for the magnets SmCo 5 [75] and SmFe 12 [76]. A x ðTÞ,A z ðTÞ for the former and A x ðTÞ < A z ðTÞ for the latter, which suggest significantly different anisotropic properties from Nd 2 Fe 14 B.

FMR
The dynamics of magnetization is simulated by the LLG equation. To take into account the thermal fluctuation, we used the SLLG equation 21,22 in which we add a white Gaussian random field f� i g to the LLG equation: where α i is the damping factor at the ith site, γ is the gyromagnetic constant, and is the effective field due to the exchange interaction and the anisotropy terms. We adopted the commonly accepted value of α i ¼ 0:1 for the damping constant [17]. The random field � i satisfies The strength of the noise, D k , is related to the temperature of the system by the fluctuation-dissipation relation The physical noise has a finite auto-correlation time, and the white-Gaussian noise whose auto-correlation time is zero is an extreme case for short correlation noise. Indeed, the white-Gaussian noise contains very high-frequency components which should be quantum mechanically suppressed at a finite temperature. Such situation has been studied in the literature [77]. However, it is also known that in the classical limit � h ! 0, the present treatment is consistent, and the present SLLG realizes the dynamics of probability distribution of the Fokker-Planck equation which leads the distribution to the thermal equilibrium distribution at a given temperature [21,22]. Thus, as a classical model for such process, we adopt the present scheme of dynamical model. Using the Kubo formula [29,78], we obtained the temperature dependence of the FMR resonance frequency at zero external field from the time-correlation function C α ðτÞ (α ¼ x; y; or z) of magnetization M ¼ ðM x ; M y ; M z Þ in equilibrium. The FMR spectrum Iðf Þ is given by the time-correlation function C α ðτÞ as where The time-correlation function was obtained by the SLLG method.
In Figure 8, we depict the temperature dependence of the resonance frequency of the Nd-magnet model. Here, we find that the peak of the spectrum shows a nonmonotonic dependence on the temperature and that the almost zero FMR frequency in the lowtemperature phase reflects SR transition [29]. Now we study why the resonance frequency is zero, f R ,0, in the low-temperature phase. This property is related to the tilted configuration. Here, we consider the dynamics of the total magnetization M of a minimal model for SR, that has exchange interactions and the first and second anisotropy energies (18), because the spins are tightly connected by the interaction and they move together. FMR is the response of the total magnetization to the external uniform field. The effective field for the precession motion applied to the total magnetization is given as where H Aniso is the contribution from the anisotropy term of the minimal model. The precession frequency is given by For D 2 ¼ 0, the resonance frequency is given by . Thus, f R is proportional to M z , which is the conventional temperature dependence in the usual FMR.
On the other hand, for D 2 �0, the situation is different. Considering the relations we note the relation because the following relation holds at the tilted configuration: Thus, we have an important consequence: In the present calculation, anisotropy in the plane given by Stevens operators with m�0, is not taken into account. If non-zero m terms exist, the effective field deviates from the z-axis and the resonance frequency is not zero. However, non-zero terms are small, and the large reduction in resonance frequency around the SR transition temperature obtained here ( Figure 8) is one of the characteristics of the present material.

Coercivity of nanoparticles
Coercivity is a threshold field for metastable magnetization, and the situation is schematically depicted in Figure 9. Coercivity at zero temperature is given by the  field at which the barrier disappears. However, at finite temperatures, the thermal fluctuation causes a jump over the barrier as shown in Figure 9(b). In the absence of a barrier (Figure 9(a)), magnetization relaxes smoothly in a deterministic manner. On the other hand, in the case of Figure 9(b), relaxation is triggered by a large thermal fluctuation and occurs stochastically as a type of Poisson process, where the distribution of the relaxation time is large. Magnets consist of grains, and magnetization reversal is a sequence of reversals of grains. Thus, as a fundamental process, we first studied the process in a single grain. At zero temperature, the reversal occurs as the Stoner-Wohlfarth process [79]. However, at finite temperatures, the effect of thermal agitation plays an important role [20]. To study this effect quantitatively, we adopted two complementary methods, that is, a direct SLLG simulation of the dynamics of magnetization [38] and an analysis of free energy as a function of magnetization FðMÞ obtained by a MC simulation [40].

Dynamics of magnetization
By using the SLLG equation (12), we studied the dynamics of magnetization. In Figure 10(a), we show snapshots of the magnetization reversal from the down-spin state for α ¼ 0:1 under a reversed field, h ¼ 4:0 T, which is in the stochastic region. There, we find that nucleation occurs from a corner. Then, the reversed region expands first in the ab plane by a Bloch-type domain wall and then grows in the direction of the c-axis by a Néel-type domain wall. We observed that this tendency is independent of α. This process is attributed to the fact that the effective exchange interactions along the a-and b-axes are stronger than that along the c-axis [24,25].
In Figure 10(b), we depict the time dependence of the magnetization m z ¼ P N site i¼1 S z i =N in the relaxation processes for α ¼ 0:1 at h ¼ 8T (Figure 10(b) (left)) and at h ¼ 4:1T (Figure 10(b) (right)). There are two typical types of relaxation, that is, deterministic and stochastic. The former type occurs for a large field, and the relaxation is characterized by a multi-nucleation Avrami process. Examples of relaxation (12 samples) of magnetization are shown in Figure 10(b) (left), where the distribution of relaxation times is small. On the other hand, the latter type occurs for a small field, and the relaxation is characterized by a single nucleation. In this case, the distribution of relaxation times is very wide as shown in Figure 10(b) (right). There, a few samples do not relax, and thus it is impossible to estimate the average relaxation time until they relax. At the border between the two cases, the relaxation time increases very rapidly, and this region of the field can be regarded as the practical end of metastability, which is called the dynamical spinodal point [80].
In the latter case, the large distribution of relaxation times prevents us from estimating the average relaxation time. To overcome this difficulty, we introduced a statistical method to evaluate the relaxation time. We derived the statistical relation between the reversal probability p and the relaxation time τ. If an event (relaxation) occurs with the probability p in a unit time, the probability that the event occurs for the first time in the period ½t; t þ Δt� is pe À pt Δt: The mean relaxation time hτi is given by The probability PðtÞ that the event occurs in the period ½0; t� is PðtÞ ¼ 1 À e À pt : If we perform N simulations, the number of surving (unchanged) samples is N sv ðtÞ ¼ N À N done ðtÞ ¼ Ne À pt : Then, p (and τ) can be estimated from the slope of lnðN sv ðtÞ=NÞ versus À pt. We plot lnðN sv ðtÞ=NÞ as a function of t in Figure 11. From the slope, we can estimate hτi. In this method, we can recognize the time range of linear dependence where the dynamics is governed by a single nucleation process, which is difficult when taking the naive average of the samples. In Figure 12(a), we give the field dependence of the relaxation time with different α values. The relaxation time increases rapidly below h ' 4:2T. For large relaxation times, we expect the single exponential decay of the Arrhenius type. Thus, when we extrapolate the increasing relaxation times, we fit them including a correction term in the form of a double exponential fitting, where C ¼ B=A and d ¼ b À a. In Figure 12 This value is close to the estimation h c ' 3:3 obtained by the MC method [40], which will be reviewed in the next subsection. The dashed line in Figure 12(b) is the estimation by the MC method. We find that, although the simulation time of the SLLG method is limited and much shorter than 1s, owing to the fact that the relaxation is governed by a single nucleation process, the extrapolation of the relaxation time using τðhÞ as a fitting function is effective for estimating coercivity.

Free energy as a function of magnetization
The free energy for a given temperature T and a field H is given by FðT; HÞ ¼ À k B T ln ZðT; HÞ; ZðT; HÞ  where H 0 is the Hamiltonian without the field, S z i is the spin at the ith site, and P N i S z i is the total magnetization. If we fix the magnetization to M, then

FðT; H; MÞ ¼ À k B T ln ZðT; H; MÞ; ZðT; H; MÞ
The probability that the system has the magnetization M is given by In Figure 13(left), FðT; H; MÞ is depicted for several fields. The free energy barrier F B is defined as depicted in Figure 13(right). The field dependences of the free energy barrier F B ðHÞ for several sizes are given in Figure 14. Here we find that the size dependence of FðT; H; MÞ saturates for L > 20nm, and they are large enough to study the cases of larger sizes. Here, we define H 0 as the field where F B becomes zero, at which the magnetization becomes unstable. The relaxation time is zero at this field. At finite temperatures, the metastable state may relax with thermal fluctuation. Using the Arrhenius model, we estimate the rate of relaxation per 1s and the relaxation time τ as where N contact is the number of contacts with the thermal bath in 1s, which is usually given as 10 11 . Thus, for the relaxation time of 1s, The field that gives this value is the coercivity for the relaxation time of 1s, which we call the thermally activated coercivity H c . The thus obtained H 0 and H c at different temperatures are plotted in Figure 15 by blue and red circles, respectively. In the figure, the coercivity H k , which is obtained as H 0 in the system with the periodic boundary condition, is plotted by a dashed line. In the calculated temperature range, we confirmed that H k takes almost the same value as the magnetic anisotropy field H a ¼ 2K 1 =M s , where K 1 is the magnetic anisotropy constant [23,24] and M s is the saturated magnetization at the given temperature. In the figure, experimentally observed coercivities for a sintered magnet [81] and a hot-deformed magnet with the grain boundary diffusion of Nd-Cu alloy [42] are also plotted. Within a nanosize grain, DDI is not relevant. However, when we compare the results with the experimental data of a magnet, we must  consider the effect of DDI from other grains as a demagnetization field. The demagnetization field À N d M s approximately represents the effects of DDI as an external uniform field. Here, we do not explicitly consider the size and shape dependences [82] of the demagnetization field, but they are expected in the range of N d μ B ¼ 0.5-1.0. If we take into account the contribution of this demagnetization effect, H c gives a good agreement with that of the hot-deformed magnet in which the grains are rather isolated.
The temperature dependence of the coercivity H coercivity has been expressed in the following forms, known as the Kronmüller equation (5,33) which indicates how much the coercivity is reduced from H k À M s N d , which is naively expected. In the first form, H t ¼ H 0 À H c is the thermal activation field and the parameter α is given by H 0 =H k , while in the second form, all the thermal effects are included in the parameter α 0 ¼ H c =H k . Note that α is not the damping factor in the LLG equation. Our approach can handle the temperature dependences of α and H t explicitly. In the inset of Figure 15, the temperature dependences of these parameters are depicted. For some sintered polycrystalline magnets (corresponding to the green line in Figure 15), Kronmüller and Durst phenomenologically estimated the decay factor as α 0 exp ¼ 0.89-0.93 from magnetic properties measured around room temperature [5], which is close to our estimation.

Activation volume
Next, we consider the mechanism of thermal activation (nucleation) for which the concept of 'activation volume' has been introduced [32][33][34]42]. ΔM represents the difference in magnetization between the local minimum M ¼ M b and the local maximum M ¼ M t of the free energy (see Figure 13(right)), i.e.
The activation volume has been defined by it is obvious that (35) gives the relation Here, we note that FðM tðbÞ ; 0Þ depends on H since M tðbÞ is a function of H, but because dF=dM ¼ 0 at M ¼ M tðbÞ , the contribution from these terms is zero. If we assume a linear dependence of the barrier on H, that is, taking n ¼ 1 in the phenomenological formula then we have the widely used phenomenological equation for thermal activation effects [33]: where V c is V at H ¼ H c . In Figure 16, we compare the temperature dependences 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 hightemperature range, which is attributed to the fact that the activation volume and n are not exactly constant.

Coercivity of large grains
Because of the uniform magnetization, DDI has a relevant effect on breaking down the uniform ordering. We have studied the effect of DDI for a thin film of the present material. Because DDI is a long-range interaction, the computational cost of simulation increases with the square of the number of spins (N). Various methods have been proposed to overcome this difficulty [83], but they are not very efficient for systems with large unit cells, such as the present material as above. For example, in the method using fast Fourier transform (FFT) a system containing 68 atoms in a unit cell requires 68 modes in Fourier space. The socalled stochastic cutoff (SCO) method has also been proposed as an alternative [84,85]. The SCO method introduces a selection of bonds by the socalled switching procedure. The selection procedure is performed stochastically, maintaining the detailed balance condition. Thus, the stationary state of the simulation is guaranteed to be the same as the equilibrium state of the original model. Because the bond update process rarely adopts long-distance weak bonds, the overall computational time is markedly decreased. As an example, for a three-dimensional system with DDI, one MC step can be computed in a time of OðβN ln NÞ, where N denotes the number of spins in the system. In this method, the complex unit cell also introduces bothersome procedures. In this simulation, we developed a modified SCO method by using of the walker's algorithm [30]. In Figure 17, we demonstrate that the computational time decreases from OðN 2 ) to O (N ln N), and that the modified SCO (MSCO) method gives a smaller coefficient than the conventional SCO method.

Effect of DDI on coercivity of nanoparticles
Thus far, we have studied coercivity without DDI.
When the system size increases, DDI becomes an important factor determining coercivity. In this subsection, we study the effect of DDI on the coercivity of nanoparticles using the method of free energy, while the effect of DDI for larger grains in which multiple magnetic domains appear will be studied in Section 5.2. As long as the nucleation starting from corners dominates, the threshold field should be independent of the size. At finite temperatures, however, super-paramagnetism reduces the threshold field, and thus the threshold increases with the size. We confirmed that the threshold field saturates when the grain size becomes larger than 20 nm, as depicted by open blue circles in Figure 18(a) where the effect of superparamagnetism stops.
The threshold fields with and without DDI are depicted in Figure 18(a) as coercivity. Here, we clearly find that DDI reduces coercivity. In Figure 18  obtained in the literatures [86][87][88][89] are plotted by black circles. We find that coercivity tends to decrease as the logarithm of the linear dimension of the grain. The data obtained in Figure 18(a) are also plotted in the figure.
There we find that the maximum point of coercivity exists as a function of grain size at finite temperature due to super-paramagnetism.

Coercivity in systems with multiple magnetic domains
For larger grains, the ferromagnetic configuration is no longer stable and a multidomain magnetic structure appears. Even in such large systems, a metastable uniform ferromagnetic state may exist. We studied this problem by using a simplified model: To investigate the parameter dependence of characteristic magnetic configurations, we surveyed the magnetic profiles under open boundary conditions with different anisotropy constants K, DDIs D, and thicknesses of systems L z . We present magnetic profiles in a phase diagram in the space ðK=J; D=JÞ for various thicknesses [44]. Five typical magnetic phases are found: out-of-plane ferromagnet, inplane ferromagnet, vortex, multidomain, and canted multidomain. We depict examples for L z ¼ 10 in Figure 19. We evaluated coercivity by comparing the magnetic profiles obtained by the following two approaches: (A) thermal-quench process in which the simulation starts from a random spin configuration at a high temperature, and then a MC update is performed at a given temperature to find a stationary state, which corresponds to the thermal demagnetization process, and (B) field-quench process in which we switch off the magnetic field and observe the evolution of the system from a saturated ferromagnetic state obtained at a high field, which corresponds to the remanence process. Figure 18. (a) Size dependence of coercivity with and without DDI. Blue and red circles denote the thermally activated coercivity with and without DDI, respectively. (b) Size dependence of coercivity (black circles) taken from literatures [86][87][88][89] adding to the data in (a). Above the dotted line, the system tends to have a magnetic multidomain structure.  Figure 19 indicates that the ferromagnetic state remains metastable in some parameter region where states obtained by the thermal-quench process are multidomain. This mechanism may give coercivity in rather large grains. The metastable ferromagnetic state collapses when D reaches a threshold. In Figure 20, we depict a configuration immediately after the collapse, where the magnetization reversal begins inside the plane. This nucleation is in significant contrast to the case of a nanoscale system, where the nucleation begins from corners [38,40].

Effect of grain boundary
The Nd 2 Fe 14 B magnet consists of hard-magnet (Nd 2 Fe 14 B) grains, each of which is covered by a grain boundary material [24,51,54,90]. Thus, it is important to study how boundary phases affect the coercivity of the grains studied in the previous section.
When we investigate nanocomposite magnets or magnets with soft-magnet defects, we consider α-iron as the boundary material. On the other hand, for sintered or hot-deformed magnets, the boundary phase consists of a Nd-rich material that shows a weak ferromagnetic property. The components of the soft phase have been studied by the concentration depth profile method [18], in which the width of the grain boundary is a few nanometers, and in most cases, it consists of a ferromagnetic soft material [15][16][17][18][19]. In the following, we study the former case in Subsection 6.1, in which the soft magnet first reverses its magnetization and reduces the coercivity of the hard magnet, which is called the spring effect. Then, we study nucleation and depinning phenomena in a sandwich structure of hard-soft-hard parts in Subsection 6.2, mainly focusing on the latter case. As a related problem, we study the effects of the surface in Subsection 6.3.

Effect of soft-magnet grain boundary on coercivity
Here, we study the spring effect due to a soft material. Because the exchange stiffness constant depends on Figure 20. Nucleation pattern just after the collapse of out-of-plane ferromagnetic state. The spin direction ðS x ; S y ; S z Þ is coded in the manner that ðS x ; S y Þ is given by color, e.g. (0.1) is red, and S z is coded by brightness of the color, i.e. the radius in the color code denotes S z from À 1 (center) to þ 1 (edge). (From reference [44]: modified.). the direction, the spring effect also depends on the direction. Thus, we studied magnetization reversal in systems with a soft phase in both directions [24]. We carried out micromagnetic simulations for two-phase models with the parameters, AðTÞ and KðTÞ for a finite temperature (T ¼ 400K), which were previously obtained in Section 3.3. The models are composed of soft and the hard mage exchange interactions among theses, respectively depicted in Figure 21(a), where the shaded parts denote the soft phases. Models A and B are the same if we do not take into account the anisotropy of A and DDI.
In Figure 21(b) (left), the dependence of coercivity on the thickness s l of the soft phase is depicted. Dashed lines denote the analytical results of the depinningtype coercivity [45,46]. In Figure 21(b) (right), the field dependences of the magnetization are plotted, in which changes in magnetization due to the nucleation in the soft layer and the depinning of the domain wall (i.e. the reversal of the hard magnet) are clearly seen.
Westmoreland, et al. [90] studied the effect of the soft magnet α-Fe by MC simulation and the SLLG equation using an atomistic Hamiltonian with thermal fluctuation. They used a simplified Hamiltonian that gives critical temperatures correctly, although the spin-reorientation transition is not realized. Using the model, they systematically studied the properties of core/shell nanocomposites with improved performance at a temperature suitable for motor applications (around 450 K).

Sandwich structure
To realize stronger coercivities at higher temperatures, it is necessary to study the effect of grain boundaries on the coercivity. For this purpose, a prototype hardsoft-hard magnet model, in which outer hard magnets are in contact with a middle soft magnet, has been intensively studied [45][46][47][48][49][50]91,92]. This model captures the essence of nucleation and depinning in inhomogeneous systems, and has been frequently used in analyses of the phenomena in various experimental and theoretical studies of magnetic materials [42,93] including GMR sensors [94].
Sakuma et al. investigated the threshold fields for nucleation and depinning in a hard-soft-hard magnet continuum model at zero temperature [45,48]. Solving a one-dimensional nonlinear equation for the model with the exchange stiffness constant and magnetocrystalline anisotropy energy, they presented a phase diagram of the threshold fields as a function of the ratios between the stiffness and anisotropy constants of the soft and hard magnets. However, thermal fluctuation effects, which are also essential for coercivity, were not studied.
Mohakud et al. [49] studied the temperature dependence of the corresponding phase diagram for the hard-soft-hard magnet model in the simple cubic lattice of the Heisenberg model with single-ion anisotropy by solving the SLLG equation (21,22). They showed various parameter dependences at different temperatures, and the threshold fields were found to be significantly affected by the thermal effect. Westmoreland et al. [50] studied this problem using atomistic and continuous spin models with temperature-dependent parameters for a (Nd magnet)-(α Fe)-(Nd magnet) system at finite temperatures and found that the thermal fluctuation reduced coercivity.
These properties were studied using the atomistic model for the Nd 2 Fe 14 B magnet [51]. We do not have precise information on the detailed structure of the soft region, and we adopted the same lattice structure with reduced coupling constants and anisotropy energies. At about half of the critical temperature, the grain boundary is close to being paramagnetic, e.g. T C has been found to be about 200 � C by X-ray magnetic circular dichroism (XMCD) [95]. Here, the details of the structure are not very relevant, and how coercivity of the hard magnet is reduced by the reversed magnetization in the grain boundary phase is an important issue. We depict the sandwich structure along the a-and c-axes in Figure 22. We define F and E=F, respectively, as the ratio of exchange interactions and that of anisotropy energies between the soft-and hard-magnet phases. In Figure 23(a), we plot for F ¼ 0:5 the threshold fields of nucleation, namely, for the process ( þ þ þ ) to ( þ À þ ) (circles) and also for the process ( þ À þ ) to ( À À À )(triangles), in both cases. In Figure 23(b), the threshold fields of depinning, namely, for the process ( þ À À ) to ( À À À ), are plotted.
It is found that the thermal fluctuation effects are considerably large in the Nd magnet, and at 300 K, the threshold fields of nucleation and depinning are much reduced and the E dependences are changed from those at T ¼ 0K [51]. Figure 23(a) and (b) show that the dependence of the threshold fields on the direction (Néel or Bloch type) is clear in the nucleation case, while the thresholds for depinning do not depend on the type. The threshold fields for the process for the process ðþ þ þÞ to ðþ À þÞ of the Bloch type are larger than those of the Néel type, which should be attributed to the stronger effective exchange interaction. The threshold fields for the process ðþ À þÞ to ðÀ À À Þ are almost the same and constant for E < 0:4, where the depinning of the reversed magnetization in region II gives the threshold. The threshold for the depinning process (right) is also almost the same and constant. In these cases, surface nucleation in region I or III is important for the depinning process, in which the strength of the magnetic anisotropy in region II is not essential, and the necessary field for nucleation is almost the same for the creation of Bloch and Néel domain walls.

Effect of surface properties
As studied in the previous subsection, domain wall depinning is a process in which the reversed magnetization invades the neighbor hard grain ( Figure 24). In this process, surface nucleation occurs under the influence of the contacting soft phase and the external field. The latter acts on the entire hard grain, while the former acts at the surface. Thus, the properties at the surface have significant effects on the coercivity of the magnet [96].
In this subsection, we study how the modification of the surface affects coercivity. For this study, we also used a system consisting of 12 � 12 � 9 unit cells. We studied the effects of modification of both the c-plane and the a-plane [52]. Here, we show the case of c-plane modification, where open and periodic boundary conditions are used along the c-axis and the a-and b-axes for the (001) surface, respectively.
It has been pointed out that at the surface, the anisotropy of a Nd atom is of the easy-plane type [97][98][99], in contrast to the easy-axis type in the bulk. On the other hand, surface easy-axis anisotropy may be enhanced by the substitution of strongly anisotropic atoms, for example, Dy [7,100]. Thus, we studied the three generic cases: (1) the anisotropy of modified layers is zero, (2) the anisotropy of modified layers is of the easy-plane type, and (3) the anisotropy of modified layers is of the enhanced easy-axis type. Concretely, we investigated coercivity with the following three settings of the anisotropy parameters ðÃ 0 2 ;Ã 0 4 ;Ã 0 6 Þ for the surface Nd atoms: (1) no anisotropy in Nd atoms: Ã 0 2 ¼Ã 0 4 ¼Ã 0 6 ¼ 0; (2) in-plane anisotropy in Nd atoms: Ã 0 2 ¼ À A 0 2 < 0 and Ã 0 4 1 A 0 6 ¼ 0; where Ã 2 0 is negative and the amplitude is of the same order as that in the bulk [97][98][99], and (3) doubly reinforced anisotropy in Nd atoms: The dependence on the number of modified layers (n) is also important. In Figure 25(a), the definition of Figure 23. (a) Threshold fields for nucleation from ðþ þ þÞ to ðþ À þÞ (circles) and from ðþ À þÞ to ðÀ À À Þ (triangles) for Bloch and Néel domain walls at 300 K. (b) Threshold field for depinning from ðþ À À Þ to ðÀ À À Þ for Bloch (squares) and Néel (crosses) domain walls. Here, the temperature is 300 K, and L 1 ¼ L 2 ¼ L 3 ¼ 12 (unit cells), and the ratio of exchange interactions of the soft and hard magnetic phases is 0.5 (F ¼ 0:5). Parameter E is the ratio of the anisotropy energies of the soft and hard magnetic phases multiplied by F. (From reference [51]: modified.). n is depicted. In Figure 25(b), the threshold fields for the magnetization reversal at T ¼ 0:46T C in the three cases are plotted as a function of n. At this temperature, the threshold field is much reduced from the value at T ¼ 0 [52,96], and also we find that single-layer modification ðn ¼ 1Þ has little effect due to thermal fluctuation. However, as n increases, the effect becomes relevant. We suppose that if the width of the modified layer is comparable to the size of the activation volume, i.e. about half of the domain-wall width (n,3), the modification becomes relevant. Thus, surface coating would be a useful method to increase coercivity.
In general, the properties of the surface of a grain are very important factors determining coercivity, and it has been reported that coercivity depends on the procedures on the surface of grains, for example, sputtering, heat treatment and grain boundary diffusion [100][101][102]. Thus, not only the modification of anisotropy, but also the surface roughness should be studied, which will be reported elsewhere [103].

Coercivity of magnets as an ensemble of grains
Magnets consist of grains, and the coercivity of a magnet must be studied by taking into account this feature. For the cooperativity of a magnet, the distribution of properties of grains (alignment, size, and shape) and also the interactions among grains including DDI play important roles.
As studied in previous sections, each grain has coercivity, which may be modeled by the threshold fields H 1 from up to down and H 2 from down to up. Such a unit representing a modeled grain is called a hysteron. Let us consider a magnetization process starting from a saturated state (all the hysterons are up) at a large positive field and reduce the field. In the process, each hysteron is reversed at its H 1 , and the total magnetization changes with the field. If, at a field H r , the field stops decreasing, and then increases, reversed hysterons remain in the down state until the field reaches H 2 of the grains. Thus, the total magnetization in the reverse process from H r MðH r ; HÞ shows a hysteresis that depends on H r , which is called first-order reversal curve (FORC).
If hysterons are independent, the ensemble of hysterons is called the Preisach model [55]. There, the joint distribution of H 1 and H 2 , PðH 1 ; H 2 Þ, is obtained from MðH r ; HÞ by the following relation: In real magnets, grains interact, and thus we cannot use the relation. However, the distribution obtained by the above relation is called a FORC diagram, which is often used to characterize the features of magnets [56][57][58]. It is an interesting problem to study how the interaction among the grains affects the distribution. The effect of the interaction has been taken into account in a mean-field-type analysis using the so-called moving (Preisach) model [59]. We are studying this dependence with an extended Preisach model (interacting Preisach model), in which the distributions of alignment, size, and shape are taken into account, which will be reported in the future [60].

Summary and discussion
We have reviewed works on finite-temperature properties by using the atomistic model. We first constructed an atomistic Hamiltonian describing the microscopic nature of the present material (2), where the intraatomic electronic structure of the Nd atom ( Figure 2) and the exchange interactions among the atoms were concretely set up. Then, the temperature dependences of various equilibrium properties were obtained. Then, the coercivity of nanograins was studied by the stochastic LLG method ( Figure 10) and also by the free-energy landscape method using the Wang-Laudau MC algorithm ( Figure 13). The temperature dependence of coercivity was also estimated quantitatively ( Figure 15). A large reduction in coercivity was found, which gives the upper limit of coercivity at a given temperature. The magnetization reversal of individual grains has been experimentally observed [104] and the effect of thermal fluctuation on the process has become a realistic problem. The effect of DDI was studied by a newly developed SCO method, and we also studied the mechanism of coercivity in large systems with magnetic multidomain structures in Section 5.2. In this manner, we have studied coercivity from the microscopic scale to nanoscale. The final goal should be the study of coercivity at the macroscopic scale. For magnets, 'macroscopic' does not simply mean a large size. Indeed, the simple enlargement of a system induces the multidomain state as mentioned above. Macroscopic magnets consist of an ensemble of grains. To study coercivity in the ensemble, we must take into account the cooperative behavior of grains. Thus, we studied the effects of grain boundary phases consisting of soft magnets in Section 6.1, and we also studied the temperature dependence of threshold fields for nucleation and depinning in the sandwich (hard-soft-hard) structure in Section 6.2. Moreover, in Section 6.3, we studied the effect of surface modification by changing the anisotropy energies of a few layers at the surface. There we found that the modification causes significant effects on coercivity. With this information on interactions among grains, in Section 7, we discussed possible extensions of the Preisach model [59,60] to study the effects of interactions on the FORC diagram which provides the characteristics of magnets [56,57].

Outlook for coercivity at several scales
Let us summarize the works from the so-called multiscale viewpoint. As shown in the present review, there are various scales at which coercivity can be studied. The main physical origin of coercivity is the anisotropy energy and the exchange energy, which is at the electronic scale (,1Å). These quantities are studied by first-principles calculations, and reviewed for the present model in Section 2. The next scale is that of a single grain (, 20 nm). The coercivity of single small grains at finite temperatures was quantitatively estimated for the first time in the works presented in this review. It was found that the reduction due to the thermal fluctuation is rather large even in such simple small systems. In larger grains, DDI makes a uniformly magnetized state unstable and generates a magnetic multidomain (maze) structure. Indeed, magnetic domains have been observed in grains with a size of more than 100 nm, which is even smaller than the size of grains in sintered magnets (typically micrometers). However, the sintered magnets still show coercivity. We reviewed a work on coercivity in such cases. This problem must be studied more carefully in the future. Finally, we must consider coercivity at the macroscopic scale, that is, magnets that are an assembly of grains. At this scale, the cooperative motion of grains due to DDI and also the interaction with the neighboring sites play important roles in coercivity. The former effect has been studied as a demagnetization effect, but the effect may not be so simple because of the peculiar form of DDI. Namely, depending on the relative positions of spins, DDI may generate negative and positive fields. The effects of the latter have been studied in relation to domain wall depinning, and some attempts to study this phenomenon are reviewed in section 6. Studying all the effects microscopically is difficult, although several attempts have been made [53,54,105,106]. Thus, the combination of studies on the effects of boundary phases including atomistic firstprinciples approaches [107,108] and studies on ensemble effects such as the interacting Preisach model should be accelerated to understand the coercivity of magnets.
Finally, we mention the numerical calculations. To properly study the thermal effects of the present material, we used the SLLG method, the free-energylandscape MC method, and a modified SCO method on the atomistic model. At the present computational capacity, we can use the atomistic model up to a size of 50 nm, which may include a single grain or a few grains. However, for further studies, we need to introduce a coarse-grained model. The relation between the atomistic model of the present work and the continuum magnetization model used in the so-called micromagnetic simulation (LLG equation) is an important problem. We constructed a continuum magnetization model by obtaining the stiffness constant and the anisotropy energy at a given temperature as reviewed in Section 3.3 [24]. However, to directly perform a simulation with the thermal fluctuation, it is necessary to introduce a model consisting of local variables. A renormalization analysis has been attempted by introducing coarse-grained magnetization with a variable length, where magnetization distribution was obtained by the Wang-Landau method. For a unit cell, the magnetization distribution PðMÞ is given by Using this variable-length magnetization, the Hamiltonian of the system may be approximated by As we see in Figure 26, the magnetization of an isolated unit cell is about 90μ B . However, in a crystal consisting of unit cells combined by ferromagnetic interactions, it was found that the average length of magnetization is increased by the interactions among them to about 116μ B , which agrees with the experimental value. To refine the parameters used to reproduce atomistically obtained properties, the tuning of renormalized parameters is necessary, which depends on the quantity on which we are focusing; details are under investigation. For real-time dynamics, the conventional MC method using the detailed balance is not adequate and LLG-type real-time dynamics is necessary. We have introduced the SCO algorithm to accelerate the calculation of systems with DDI. However, the applicability of the SCO method, which is based on the time evolution given by the detailed balance, to the LLG equation is not known. We have developed the socalled time-quantified MC method in which real dynamics such as the precession of spins are correctly simulated [109]. The developments of these techniques are also desirable for further studies.