Anisotropy of exchange stiffness based on atomic-scale magnetic properties in rare-earth permanent magnet Nd$_2$Fe$_{14}$B

We examine the anisotropic properties of the exchange stiffness constant, $\mathcal{A}$, for rare-earth permanent magnet, Nd$_2$Fe$_{14}$B, by connecting analyses with two different scales of length, i.e., Monte Carlo (MC) method with an atomistic spin model and Landau-Lifshitz-Gilbert (LLG) equation with a continuous magnetic model. The atomistic MC simulations are performed on the spin model of Nd$_2$Fe$_{14}$B constructed from ab-initio calculations, and the LLG micromagnetics simulations are performed with the parameters obtained by the MC simulations. We clarify that the amplitude and the thermal property of $\mathcal{A}$ depend on the orientation in the crystal, which are attributed to the layered structure of Nd atoms and weak exchange couplings between Nd and Fe atoms. We also confirm that the anisotropy of $\mathcal{A}$ significantly affects the threshold field for the magnetization reversal (coercivity) given by the depinning process.


I. INTRODUCTION
In rare-earth permanent magnets such as Nd 2 Fe 14 B, SmCo 5 , Sm 2 Fe 17 N 3 , and so on, the rare-earth elements and the transition metals combine in atomic scale to produce strong magnetic anisotropy and strong ferromagnetic order. The strong magnetic anisotropy originated from the 4f -electrons of the rare-earth atom is maintained at room temperature (RT) by an interaction with the strong ferromagnetic order of the transition metals [1,2]. However, in Nd 2 Fe 14 B which is well known as the highest performance permanent magnet at RT [3][4][5], improvement of coercivity at high temperature is required for industrial application [6] because it has a relatively low Curie temperature (∼ 585 K) as compared with the other rare-earth magnets.
Coercivity mechanism of Nd 2 Fe 14 B and the other rareearth magnets has not been fully elucidated yet. To analyze intrinsic magnetic properties of the rare-earth magnets (e.g., magnetization, magnetic anisotropy, exchange stiffness, Curie temperature, etc.), we should consider the inhomogeneous magnetic structures and thermal fluctuations in an atomic scale (∼Å). The rare-earth magnets in the practical use are polycrystalline materials composed of the main rare-earth magnet phase (∼ µm), and magnetic or non-magnetic grain boundary phases (∼ nm). Therefore, the coercivity depends on not only the intrinsic magnetic parameters, such as the anisotropy energy of the main phase but also on the magnetic properties of the grain boundary and microstructure of the main and grain phases [7,8]. For this reason, study of the coercivity requires analyses from the atomic scale (∼Å) to the macroscopic scale (∼ µm). This fact has inhibited to elucidate the coercivity in a systematic manner.
Theoretically analyses for the coercivity in most have been carried out with a continuum model that uses the magnetic anisotropy constants, K u , and the exchange stiffness constant, A as macroscopic parameters (see Eq. (2)). In the junction systems consisting of hard and soft magnetic phases [9][10][11][12], it has been indicated that both K u and A have a large effect on incoherent magnetization reversals (i.e., nucleation and depinning) which reduce the coercivity from the value of the uniform reversal. In order to elucidate the coercivity mechanism and improve the coercivity at high temperatures, it is important to clarify the temperature dependences of K u and A. Therefore, many researchers in both experimental and theoretical have studied the thermal properties of K u [13][14][15][16][17]. However, regarding A for Nd 2 Fe 14 B, experimentally, the values are observed only at several temperatures [18][19][20][21], and also there is no theoretical estimation of the temperature dependence as far as we know.
The value of A is given as a macroscopic properties of exchange stiffness of continuous magnets at each temperature. This quantity is related to the domain wall (DW) width [22] and the critical (magnetization reversal) nucleus size [7]. For some other magnetic materials, YCo 5 and L1 0 -type magnets (CoPt, FePd, FePt), Belashchenko indicated using ab-initio calculation that A depends on the orientation in the crystals [23]. And Fukazawa et al. also pointed out the anisotropic A for Sm(Fe, Co) 12 compounds [24]. Recently, Nishino et al. examined the temperature dependence of DW width of Nd 2 Fe 14 B using an atomistic spin model constructed from ab-initio calculations, and indicated that it also has an orientation dependence [25]. Although, in Ref. [23], the effect of anisotropic A on coercivity is also discussed in a phenomenological way, microscopic understanding of the temperature and orientation dependences of A is essential for the coercivity mechanism in Nd 2 Fe 14 B magnet.
In the present paper, by a comparison of the results obtained by a continuum model and those obtained by the atomistic spin model constructed from ab-initio calculations (same spin model as Ref. [17,25,26]), we study the temperature and the orientation dependence of A for Nd 2 Fe 14 B magnet. The comparison of the two different scale models is done by making use of magnetic DW energy. This scheme has been used successful for FePt [27][28][29] and Co [30] magnetic materials. We find that A has the anisotropic property in Nd 2 Fe 14 B. Moreover, the reason for the anisotropy is attributed to the weak exchange couplings of Nd and Fe atoms. Additionally, we performed micromagnetics simulations on the configuration in which a soft magnetic phase is attached to (001) plane or (100) plane of a grain of Nd 2 Fe 14 B. The simulations predict that the anisotropy of A reduces the coercivity for the former configuration because A of zdirection is larger than that of x-direction, while it enhances the coercivity for the latter configuration. The series of the calculation schemes in the present study corresponds to the multiscale analysis [31,32] that connects different scales from ab-initio calculation to coercivity as macroscopic physics.

A. Atomistic Spin Model
To include the information of electronic states at an atomistic scale, we treat the (spin) magnetic moment of each atom as a classical spin and then construct a classical Heisenberg model. The atomistic Hamiltonian has the form: whereJ ij is the Heisenberg exchange coupling constants including the spin amplitudes (S i S j ) between the ith and jth sites, and e i is the normalized spin moment at the ith site. The coefficient, D i , in the second term denotes the strength of the magnetic anisotropy of Fe sites. The third term is the magnetic anisotropy of Nd sites, which is formulated by the crystal field theory for 4f -electrons [33,34], B m l l,i is the crystal electric field coefficient andÔ m l l,i is the Stevens operator. In the present study, B m l l,i takes a fixed value, whereasÔ m l l,i depends on the state of a total angular momentum of 4f -electrons, J i . Here, we fix J i parallel to the normalized spin moment on the ith site, i.e. J i = J i e i (J i = 9/2 for Nd atom). For simplicity, we consider only the diagonal terms m l = 0.
For these input parameters of Nd 2 Fe 14 B magnet, we adopt the same values in the previous studies [17,25,26].
Exchange coupling constants,J ij , were calculated with Liechtenstein's formula [35] on the Korringa-Kohn-Rostoker (KKR) Green's-function code, machikaneyama (akaikkr) [36]. In the present study, to reduce computational cost, we use only short-range exchange couplings within the range of r cut = 3.52Å, in which primary Fe-Fe and Nd-Fe interactions are included. Anisotropy terms, D i and B m l l,i , were determined from the previous first-principles calculation [37] and the experimental result [34], respectively. Consequently, the atomistic spin model uses many input parameters in a unit cell which includes 68 atoms (see Fig. 1 (a)). The previous study [17] confirmed that the model and parameters are highly reliable for the magnetic properties of the Nd 2 Fe 14 B magnet.

B. Continuum Model
Under the continuum approximation, the micromagnetic energy of the exchange couplings and the magnetic anisotropies at temperature T is expressed as follows [38]: where m denotes a normalized magnetization vector, A l (T ) is the exchange stiffness constant of each direction (x, y, z) in the crystal, and E K is the magnetic anisotropy energy which is usually expressed as: where K u (T ) (u = 1, 2, 4) are the magnetic anisotropy constants and θ is the angle of magnetization measured from z-axis. The continuum model uses the temperaturedependent parameters to express the thermal effects instead of thermal fluctuations in the atomistic model. Micromagnetic simulations [38,39] have been carried out based on the continuum approximation and Landau-Lifshitz-Gilbert (LLG) equation which describes time evolution of magnetization relaxation process (see Sec.III C). In the practical simulation, the real space in Eq. (2) is discretized with a grid whose width should be smaller than a DW width and a magnetostatic exchange length [38]. In most micromagnetic simulations for Nd 2 Fe 14 B magnet, each grid size was set to (1-2 nm) 3 , and as the input parameters, experimentally observed values of A l (T ) and K u (T ) were used.
In the continuum model, many input parameters of Eq. (1) in atomic scale are expected to be renormalized in the few macroscopic parameters of Eq. (2) at each temperature.  [5]. Two types of spin configurations for Nd2Fe14B atomistic spin model: DW orientation is along (b) the yz-plane (Type I) and (c) the xy-plane (Type II). The arrows in figures describe the antiparallel (AP) boundary condition. These crystal structures were plotted using vesta [40].

C. Methods for determining Ku(T ) and A(T )
In the present study, the macroscopic parameters are evaluated by comparing the energies of the DW and of the magnetic anisotropy obtained by the above two models at each temperature. We use a Monte Carlo (MC) scheme to calculate Helmholtz free energies on the atomistic spin model. In order to evaluate A, we consider the two quantities obtained from DW energy: E dw (T ) and F dw (T ) [22,28]. E dw is the DW energy in the continuum model, which is expressed by the following formula: To obtain A(T ), we need the values of E dw (T ) and E K (T, θ). We regard E dw (T ) to be equal to the DW free energy in the atomistic spin model, F dw (T ), which is ex-pressed by the following formula: The DW internal energy, E dw (T ), is defined as an energy difference between the internal energies obtained in different boundary conditions with and without DW. For this calculation, we fix the boundaries in the direction antiparallel (AP) or parallel (PA). Additionally, since the crystal structure of Nd 2 Fe 14 B has anisotropy, we consider two models with the two fixed boundaries depicted in Fig. 1(b) and (c). Thus, we have two types of DWs [25].
One is a DW along the yz-plane (type I), and the other is along the xy-plane (type II). We set the periodic boundary condition in the yz(xy)-plane for the model of type I (type II).
To calculate E K (T, θ), we adopt the constrained MC (C-MC) method for the atomistic spin model. The C-MC method samples spin configurations under the condition of a fixed angle of total magnetization, which allows us to calculate the angle dependence of the spin torque and the free energy [41]. By considering that E K is equal to the free energy of magnetic anisotropy, we calculate E K from the atomistic model. The previous study [17] showed the validity of the C-MC methods for calculating (4) can be integrated analytically and gives the well known relation: In the present paper, the MC and C-MC simulations run 100000-200000 MC steps for equilibrium and the following 100000-200000 MC steps for taking statistical averages, where a MC step corresponds to one trial for each spin to be updated. We calculated the average of magnetic anisotropy energies and DW energies from 8-20 different runs with different initial conditions and random sequences.

III. RESULTS AND DISCUSSION
A. Magnetic anisotropy and domain wall energy As mentioned above, the evaluation of the exchange stiffness constant requires precise calculations of the magnetic anisotropy constants and the DW energy. First, we show the temperature dependence of the magnetic anisotropy constants in Fig. 2, which are calculated using the C-MC method. There, we used the system of 10×10×7 unit cells imposing the periodic boundary conditions. In uniaxial anisotropy, we study the quantities: Figure 2 shows that F A (T ) rapidly decreases with the temperature below 400 K. This temperature dependence is understood as a result of fragile thermal properties of Nd atoms. That is, while at low temperature, the anisotropy of Nd atoms is dominant, the anisotropy of Nd atoms decreases with the temperature more rapidly than that of Fe atoms, because the exchange field of a Nd atom (sum of exchange couplings that connect to a Nd atom) is much smaller (about 20-25 %) than that of a Fe atom (see also the detailed discussion in Ref. [17]). Additionally, owing to the higher order terms of the Nd anisotropy, B 0 4 , B 0 6 , the magnetization of Nd 2 Fe 14 B is tilted from the z-axis at low temperatures. The deviation from uniaxial anisotropy occurs in the region of −2K 2 < K 1 < 0 [42]. Note that Eq. (4) is defined under the condition of uniaxial anisotropy. And thus, in the present paper, we discuss the DW energy and the exchange stiffness in the region of T ≥ 200 K.
As a matter of practice, finite-size effect of numerical results about the magnetic anisotropy is significant for the evaluation of the exchange stiffness and the DW width, especially at high temperatures. To avoid this problem, we extrapolate F A (T ) and π 0 dθ E K (T, θ) [these are used in Eqs. (7) and (4)] to the thermodynamic limit (N → ∞) using linear functions in N −1/3 . Figure 3 shows the typical fitting results which indicate that the MC results are well fitted with the linear functions. The extrapolated results for F A near the Curie temperature, T C , are summarized in the inset of Fig. 2, where we fix as K u = 0 above T C (= 870 K). We use these extrapolated results to evaluate the exchange stiffness and the DW width.
Next, we focus on the DW energy. Figure 4 shows the temperature dependence of the DW internal energy, E dw , for the two DW types (see Fig. 1). We also plot E dw for three different system sizes in the directions perpendicular to the DWs, i.e. L x =14, 17, and 21 (L z =10, 12, and 15), for type I (type II). From these results, it is confirmed that these systems are large enough to calculate E dw . To analyze in detail the temperature dependence, we divided E dw into the magnetic anisotropy term (E K dw ) and the exchange term (E J dw ), and plot the contributions in Fig. 4. The DW energies of type I and type II show a qualitatively similar temperature dependence. Both energy terms naturally vanish for T ≥ T C , because DW does not appear even in the AP boundary condition. For T < T C , the anisotropy term decreases monotonically as the temperature increases whereas the exchange term increases and takes a peak value at the temperature, T h , slightly below T C .
The temperature dependences are interpreted as follows. The monotonic decrease of E K dw is merely due to the decrease of the thermally averaged magnetic moments, which also corresponds to the decrease of F A (T ) in Fig. 2. The magnetic anisotropy energy depends on the angle from z-axis of each spin, whereas the exchange coupling energy depends on the relative angle of spin pairs. Now, the DW energy is defined as the difference of the internal energies between the PA (E PA ) and the AP dw is the part coming from the exchange term, thus it is influenced by the difference between the two boundary conditions in the fluctuations of the relative angles of spin pairs. In the configuration of the AP boundary condition, DW exists. In the DW, the ferromagnetic order is weakened because the spin configuration is forcibly twisted, and the DW has a spiral non-collinear structure with the perpendicular (xy) component at low temperatures. This fact gives the difference of the energies E J dw . We expect that E J dw due to the formation of the noncollinear structure would be enhanced near the critical point where the width of DW increases.
At a temperature T h which is slightly lower than T C , the profile of perpendicular (xy) moments is destroyed by the thermal fluctuation. As a thermally averaged magnetic structure, the non-collinear structure (the spiral structure) becomes a collinear structure. At this point, E J dw takes the peak value at T h and then decreases rapidly and approaches zero towards T C . This break of the non-collinear structure was discussed in previous studies as the disappearance of an xy-magnetization inside the DW [28,43,44].
By those internal DW energies and Eq. (5), we can calculate DW free energies, F dw . In Fig. 5, we plot the temperature dependence of F dw and E dw (the same data shown in Fig. 4) for both DW types. The differences between F dw and E dw correspond to the contribution of magnetic entropy. The difference in the DW energy between type I and type II naturally indicates an anisotropy concerning to the direction in which DWs are generated. The DW prefers to be generated in the configuration of type II. These observations imply the magnetization re-FIG. 5. Temperature dependence of DW free energy, F dw and internal energy, E dw (same as Fig. 4 (a) and (b)), for each DW type. versal starts from the z-plane. Moreover, the generation of the DW in the magnets would also depend on the grain boundary phase and the dipole-dipole interaction. We will discuss these properties in more realistic magnetization reversal process in Sec. III C.

B. Exchange stiffness constant
The exchange stiffness constants, A, for the two directions can be evaluated by using Eq. (4) and the numerical results for the magnetic anisotropy and the DW energy, whose temperature dependence are shown in Fig. 6. Here, we define the value of A calculated from the configuration of type I (type II) as the exchange stiffness constant of the x (z)-direction, A x(z) . Reflecting the anisotropy of the DW energy, the exchange stiffness constant naturally has the anisotropy depending on the direction in the crystal. For the comparison of A with experimental values, it is reasonable to normalize the temperature dependence with T C , because the spin model overestimates (T MC C = 870 K) compared to experiment (T EXP C ≃ 585 K). In addition, as pointed in the mean field approximation [45], A is roughly proportional to T C , so that we also rescale the values of A: The inset of Fig. 6 shows the rescaled data and experimental values (green bar) at RT [20,21]. Although the experimental values have some variation (6.6-12.5 pJ/m), the calculation results are well consistent with the lower experimental values at RT. With A and K u which have been obtained, the DW width, δ dw , is calculated from the following relation [22,46] : To confirm the validity of our evaluation procedure for A, in Fig. 7, we compared our results (circle) with those of the previous study (square) [25].
In the previous study, δ dw was evaluated directly from the snapshots of spin configurations using the same atomistic Hamiltonian of the present study. Note that, we set the periodic boundary condition in the yz(xy)-plane for the model of type I (type II), whereas the previous calculation was performed under the open-boundary conditions for both models. Our results of δ dw qualitatively agree with the previous results, although they tend to take a smaller value because the thermal fluctuation becomes smaller than the previous study due to the difference of boundary conditions. The comparisons with the previous experiments and the numerical study as mentioned above guarantee our results concerning to A. Now, we study its thermal properties. The temperature dependence of A is often discussed in relation to magnetization by using the following expression: where M (T ) is the amplitude of the magnetization (not shown, see Ref. [17]). Under a mean field approximation with homogeneous spin systems, the exponent n is 2 [29,45]. However, under more accurate methods, the exponent takes a value different from 2, for example, FePt: n = 1.76 [29], hcp-Co: n = 1.79 [30]. Thus we estimate n to examine the thermal properties for the Nd 2  Fig. 6, respectively. Although the fitting results have a significant variation, the important point here is that n z always takes larger values than n x regardless of the fitting range (we checked it). This relation implies the macroscopic exchange coupling in the z-direction is weaker than that in x-direction, not only for the coupling strength but also for thermal tolerance. As the reason for the anisotropy of A to crystal orientation, the crystal structure of Nd 2 Fe 14 B is naturally invoked. Nd 2 Fe 14 B has the layered structure of Felayer and NdFe-layer (B has little effect on magnetic properties) along the z-axis as shown in Fig. 1 (a). In Nd 2 Fe 14 B, exchange couplings (J ij in Eq. (1)) are mainly contributed by bonding between Fe and Fe atoms,J Fe-Fe , and between Nd and Fe atoms,J Nd-Fe (|J Nd-Nd | is negligibly small). Each bond ofJ Fe-Fe has much larger amplitudes thanJ Nd-Fe (J Fe-Fe : −4.35 -22.34 meV,J Nd-Fe : −0.16 -3.55 meV). Therefore, it is anticipated that the anisotropy of A comes from the inhomogeneous distribution of the exchange couplings and the atom positions in the crystal structure. To support this anticipation in detail, we examine the relation between the anisotropy and the strength ofJ Nd-Fe . Figure 8 shows the ratio of A z to A x for four models with different values ofJ Nd-Fe . Beside the default model (1.0J Nd-Fe , the ratio is calculated from Fig. 6), we also calculate other three cases with all the bondsJ Nd-Fe reduced by half (0.5J Nd-Fe ), increased by half (1.5J Nd-Fe ), and doubled (2.0J Nd-Fe ). It is clearly found that A gets close to isotropic (i.e., A z = A x ) asJ Nd-Fe increase in the whole temperature range.
As another feature, the ratio slowly decreases with the temperature for all the cases, which indicates that the temperature dependence of A x and A z are different. This difference corresponds to the difference between n x and n z in Eq. (8). As the temperature increases, the contribution ofJ Nd-Fe to A becomes smaller than that ofJ Fe-Fe because the spin moments of Nd atoms are more easily broken by thermal fluctuations compared with those of Fe atoms. From the above analysis, we conclude that the reason of the anisotropy of A in Nd 2 Fe 14 B comes from the weakness ofJ Nd-Fe and the layered structure of Nd atoms.

C. Effect of anisotropic exchange on coercivity
Let us consider the effects of the anisotropy of A on the coercivity. In actual rare-earth magnets which are composed of main rare-earth magnet phase and (magnetic or non-magnetic) grain boundary phase, magnetization reversal is considered to occur by nucleation near the interface and by the DW propagation. Thus, to study magnetization reversal in such a process, we carried out micromagnetic simulations for the two-phase models composed of the soft magnetic phase and the hard magnetic phase, depicted in Fig. 9 (a) and Fig. 9 (b). The soft phase represents the grain boundary phase. The two models (a) and (b) are the same if we do not take into account the anisotropy of A and the dipole-dipole interaction.
The simulations are based on the finite-difference method and the LLG equation [38,39]: where M i is the magnetization vector of ith cell, γ is the gyromagnetic ratio constant, and α is Gilbert damping constant. Both models (a) and (b) are discretized with a cubic cell of (1.0 nm) 3 and we set |γ| = 2.21×10 5 m/A·sec (the value of free electron) and α = 1 (coercivity does not depend on α). Effective magnetic field on ith cell, H eff i , is defined as the derivative of the micromagnetic energy, E cont i (obtained from Eq.(2)), with respect to M i [47]: where, µ 0 is the magnetic permeability of the vacuum, H ext is an external magnetic field, j represents six nearest neighbor cells of the ith cell, m i = M i /|M i |, e u is the unit vector of the easy axis, d ij is the distance between the centers of ith and jth cells (i.e. d ij = 1.0 nm), A ij is the exchange stiffness costant, and K i 1 is the magnetic anisotropy constant (terms of K i 2 and K i 4 were omitted for simplicity).
In the present study, we determine the input model parameters in Eq.(10) from the MC results at 400 K. The model parameters in the hard phases are set to |M i | = 1.38 T, K i 1 = 2.63 MJ/m 3 , A ij = 12.21 pJ/m for the pairs of i z = j z , and A ij = 9.10 pJ/m for the pairs of i z = j z , where i z (j z ) is the position of the i(j)th cell in the zaxis. In the soft phases, the model parameters are set to |M i | = 1.38 T, K i 1 = 0 MJ/m 3 , and A ij = 9.10 pJ/m for all the pairs of (i, j). The difference of A ij in the direction for the hard phases reflects the anisotropy of A in the MC results. In addition, we assume A ij = 9.10 pJ/m for the bonds connecting the soft/hard interfaces.
By applying the fourth-order Runge-Kutta algorithm [48] to the LLG equation with the aboveconstructed models, we simulated the magnetization reversal dynamics and then evaluated the coercivity. In such the two-phase models, a magnetization reversal is expected to start from the soft phase. Thus, we also examine the influence of soft phase thickness, s l . In Fig. 9 (c), we plot s l dependence of the coercivity for the models (a) and (b). It is clearly seen that the coercivity of the model (a) is weaker than that of the model (b) regardless of s l . The relation of the coercivities between the two models is also consistent with that of the magnitude of DW energy in each direction (see Fig. 5). Since the models (a) and (b) are equivalent in the absence of anisotropy of A ij in the hard phase, we can conclude that the difference of the coercivity is attributed to the anisotropic A. It is also seen that as s l increases, the coercivity decreases and gradually approaches a certain value for each model, and the difference between two models increases. As we will see in the following two paragraphs, this behavior is explained as a change in magnetization reversal mechanism from nucleation type to depinning type. Here we define the nucleation type as a magnetization reversal of the whole system occurs by a nucleation which starts from nucleation at the surface of the soft phase without depinning at the interface of the soft and hard parts, while in the depinning type, the reversed magnetization in the soft phase is pinned at the interface until the magnetic field reaches the threshold of the depinning. Figure 9 (d) shows the hysteresis loops (showing the only upper part in the figure) for the model (a) with four different s l . When the soft phase is thin (s l = 1, 3 nm), the magnetization of the hard phase is reversed at the same time as the nucleation at the soft phase (i.e., the nucleation type), whereas when the soft phase is thick (s l = 5, 7 nm), coercivity is determined not from the nucleation but from the depinning (i.e., the depinning type). The change between the nucleation type and depinning type were also pointed out in the previous studies for a one-dimensional model in which a soft phase of finite width is sandwiched between hard phases [9][10][11][12] and a corner defect model [49].
In the one-dimentional model an analytical solution of the coercivity in the limit s l → ∞, which means the depinning type, was also proposed as follows [10,11]: where |M S(H) |, K S(H) 1 , and A S(H) are the model parameters of Eq. (10) in soft (hard) phase. To apply Eq. (11) to the three-dimensional models, as the value of A H , we use the value of A ij in the direction perpendicular to the soft/hard interface, and calculate H dwp . Dashed lines in Fig. 9 (c) indicate H dwp for the models, which seems to explain well the lower limit of depinning type coercivity for our three-dimensional models. Therefore, it is understood that the anisotropy of A has a large effect on the coercivity of the depinning type compared with that of the nucleation type.
Finally, we study the influence of the dipole-dipole interaction on the coercivity in the two models. Magnetic field due to the dipole-dipole interaction is incorporated in H eff i as the following form: where K(r ij ) is the demagnetization tensor [39], r ij is the distance vector between i and j cells. Here the strength of the dipole-dipole interaction is determined automatically according to the distance of the pair and the magnetization vector of the cells. We calculate H dip i in O(N c logN c ) computational time (N c is the total number of cells) by solving convolution integral using the fast Fourier transform method [50]. In the present study, we set the models (a) and (b) in Fig. 9 are cubic regardless of s l and M H = M S , and thus the models (a) and (b) have the same shape magnetic anisotropy.
Upper figures of Fig. 10 show s l dependence of coercivity for the two models (a) and (b) with dipole-dipole interaction. The red circles and the red dotted lines (H dwp ) represent the values of the coercivity which are calculated using the same input parameters (|M i |, K i 1 , and A ij ) as those in Fig. 9 (c). Because dipole-dipole interaction prefers to construct the DW along z-axis (type I), in the model (b), the coercivity decreases compared with that in Fig. 9 (c). Conversely in the model (a), the coercivity increases in the region of s l ≥ 3nm. The dipole-dipole interaction inhibits to construct the DW (nucleation) in xy-plane.
These behaviors are confirmed from the hysteresis loops in lower figures in Fig. 10. In the model (b), magnetization reversal is clearly separated into the two parts, i.e., the small jump at the lower magnetic field where only the soft phase is reversed (nucleation), and the jump at the higher magnetic field where the DW is depinned at the interface (depinning). In contrast, in the model (a) they are not clearly distinguished. That is, in the model (a), the magnetization reversal mechanism becomes to approach the nucleation type from the depinning type by the dipole-dipole interaction. The dependence of the coercivity on the arrangement of the soft phase were similarly discussed in the most recent study (not including the anisotropy of A) [51].
It is difficult to clarify the effect of anisotropic A on the coercivity under the dipole-dipole interaction by a simple comparison the models (a) and (b). Thus, we exchanged the values of A x,y and A z in the hard phase. Namely, we set the input parameters of the hard phase as A ij = 9.10 pJ/m for i z = j z , and A ij = 12.21 pJ/m for i z = j z . The values of coercivity under these conditions are plotted by the blue circles and the blue lines (H dwp ) in the upper figures of Fig. 10. The anisotropy of A has a similar effect on the coercivity as the case without dipoledipole interaction. However, in the model (a), the difference in coercivity is relatively small. The anisotropic A strongly affects the coercivity of the depinning type com-pared with the nucleation type. On the other hand, The magnetization reversal in the model (a) is the nucleation type rather than the depinning type. Therefore, we may conclude the dipole-dipole interaction in the model (a) suppresses the effect of the anisotropy of A.

IV. CONCLUSION
Regarding the exchange stiffness constant, A, of Nd 2 Fe 14 B, we examined the temperature and orientation dependences using the Monte Carlo simulations with the atomistic spin model constructed from the ab-initio calculation. We also conducted the coercivity calculations based on the micromagnetics (LLG) simulations using the continuum model with the parameters obtained by the atomic scale MC results at finite temperatures. In this way, we confirmed that the lattice structure in the atomic scale affects the coercivity as macroscopic physics. We found that A(T ) depends on the orientation of the crystal with respect to not only the amplitude but also the exponent n x(z) in the scaling behavior: and n x < n z . It is quantitatively confirmed that the anisotropic properties of A come from the weak exchange couplings between Nd and Fe atoms and the layered structure of Nd atoms. Moreover, we also found that the anisotropic A(T ) affects the coercivity of the depinning mechanism.
We focused on only Nd 2 Fe 14 B magnet in the present paper. However, the essence of anisotropic A comes from the weak exchange coupling between rare-earth atoms and transition metals and the layered structure of rare-earth atoms. Thus, the features discussion for Nd 2 Fe 14 B are probably applicable to other rare-earth magnets. In fact, it was pointed out by ab-initio calculations that A have strong anisotropy for YCo 5 [23] and Sm(Fe,Co) 12 [24].
Let us consider the coercivity from the viewpoint of the exchange spring magnet [52,53] which is composed of hard and soft phase and expected to realize the highest performance magnet. Because realization of high performance requires a large coercivity and a large thickness of soft phase, the model (a) with A x < A z (Fig. 10 (a) blue line) is the most suitable conditions in our modelings. Most strong permanent magnets, Nd 2 Fe 14 B, YCo 5 and also L1 0 -type magnet (CoPt, FePd, FePt) cannot reproduce the same condition because A x > A z [23], whereas Sm(Fe,Co) 12 would do because the anisotropy as A x < A z [24]. Therefore, Sm(Fe,Co) 12 and other R(Fe,Co) 12 -type compounds (R is a rare-earth element) may have higher potential to realize strong performance exchange spring magnet rather than the other magnets.
Finally, we point out another source of the anisotropy. In a recent experiment, it has been observed that the grain boundary phase takes different crystal structures and chemical compositions depending on the orientation with the Nd 2 Fe 14 B main phase, i.e., the Nd-rich crys-talline paramagnetic phase form on the xy-plane of the main phase, whereas the Fe-rich amorphous ferromagnetic phase in the plane parallel to the z-axis [54]. In the present paper, we have studied the effect of anisotropy of A on coercivity and of the orientation of the interface with the soft phase changes by using the same interaction for the interface. However, if the chemical structure is different, the exchange interaction would be a difference due to another source of the anisotropy, which is studied with information of the structure in the future.
In the continuum model, Eq. (2), we used the values of K u and A obtained by the MC simulations which are the values of the bulk system. There, changes in atomic scale of magnetic anisotropy [55][56][57] and exchange coupling [58][59][60][61] near the interface or surface were not taken into consideration. The influences of interface and surface are important for the coercivity [6]. Thus, the accuracy multiscale analysis needs further development of connecting scheme from atomistic spin model to macrospin model would be necessary.
As another point noted is the range of exchange interaction. In the atomistic spin model, Eq. (1), of the present study, we omitted the long-range contribution of J ij for simplicity and reduction of calculation cost. However, the recent study [62] reported that RKKY-type exchange coupling significantly effects on the DW width, and pointed out the importance of the long-range contribution. We also found that, by incorporating long-range exchange couplings up to 10.6Å, the difference between type I and type II of E dw at 400 K reduces from 11.3 % (in Fig. 5) to 4.4 %. However, to make study the effect clearly, we need precise information of the interaction at the long distance, and we postpone to study this problem later.
Although there are still problems that must be concerned, we believe that the present paper will be helpful to elucidate coercivity mechanism in rare-earth permanent magnets.

ACKNOWLEDGMENTS
We acknowledge collaboration and fruitful discussions with Taro Fukazawa, Taichi Hinokihara, Shotaro Doi, Munehisa Matsumoto, Hisazumi Akai, and Satoshi Hirosawa. This work was partly supported by Elements Strategy Initiative Center for Magnetic Materials (ESICMM) under the auspices of MEXT; by MEXT as a social and scientific priority issue (Creation of New Functional Devices and High-Performance Materials to Support Next-Generation Industries; CDMSI) to be tackled by using a post-K computer. The computation was performed on Numerical Materials Simulator at NIMS; the facilities of the Supercomputer Center, the Institute for Solid State Physics, the University of Tokyo; the supercomputer of ACCMS, Kyoto University.