Title Efficient methods for macroscopic magnetization simulation described by the assembly of simplified domain structure models

This article presents two methods for the fast computation of macroscopic magnetization model called assembled domain structure model. First, an efficient method for computing the demagnetizing field is proposed. Secondly, a direct searching method of equilibrium point is developed, which greatly reduces the computation time.


Introduction
Macroscopic magnetic properties of iron-core material result from multiscale magnetization processes such as the microscopic domain-wall motion and the mesoscopic domain-structure transition. Recently, to construct a physical macroscopic magnetization model, several energy based multiscale approaches [1]- [4] have been developed. For example, Ref. [1] has successfully represented the macroscopic anhysteretic magnetic property of the grain oriented silicon steel sheet including the magnetoelastic property. In Refs. [3] and [4], an assembly of simple domain structure models (SDSMs) represented the macroscopic hysteretic behavior of magnetic sheets.
The SDSM [5] is a mesoscopic magnetization model of crystal-grain scale describing domain-wall motion and magnetization rotation. The assembly of SDSMs [3] is expected to constitute a physical macroscopic magnetization model based on the local energy minimization. However, the assembly of large number of SDSMs requires long computation time because of the large computational cost for obtaining the demagnetizing field and the long transient process to an equilibrium point.
This article proposes an efficient method for the computation of demagnetizing field and develops a direct searching method of equilibrium point.

Simplified domain structure model
An SDSM with two domains [5], as shown in Fig. 1(a) is used to describe behavior of a mesoscopic magnetic particle, where the magnetization is assumed uniform in each domain i (i = 1, 2). The normalized magnetization vector in domain i is given by m i = (sinθ i cosφ i , sinθ i sinφ i , cosθ i ).
The total magnetic energy, e, is assumed to be given by the summation of Zeeman energy, the crystalline anisotropic energy, the domain-wall energy, and the magnetostatic energy as is summarized in Appendix A.1.
The magnetization is determined by finding a local energy minimum that satisfies e/X = 0 where X = (θ1, φ1, θ2, φ2, λ) and λ is the volume ratio of domain 1. In Ref. [3], a local minimum is obtained by finding an equilibrium point of artificial state equation given as (1) where β is a dissipation coefficient. A local energy minimum is obtained by the numerical integration of (1) until reaching the steady state where dX/dt = dY/dt = 0. If there are several equilibrium points of the state equation, one of them is obtained depending on the initial condition of (1) that reflects the history of past magnetization.

Assembled domain structure model
The macroscopic magnetization model is constituted by assembling the SDSMs [3], as shown in Fig. 1 (b), which is called the assembled domain structure model (ADSM). Each SDSM composing the ADSM is called a cell. The Zeeman energy, the anisotropic energy, and the domainwall energy of cells are independently summed up to obtain the components of total energy e in the ADSM.
In the same way as in the micromagnetic simulation (see Appendix A.2), the magnetostatic energy is given as ( 2 ) where i is the cell index and h st is the normalized demagnetizing field; h st is given as where s(i−i') is the normalized demagnetizing coefficient matrix between the cells i and i' (see Appendix A.2) and m(i) is the normalized magnetization of the cell i.
The state variable vector X consists of X(i) (i = 1, …) in each cell. A local energy minimum point is obtained by solving (1).

Computation using field decomposition
The convolution (3) is efficiently executed by using the fast Fourier transform (FFT) [6], which, however, often requires a large computational cost even with the use of FFT.
A simple way to reduce the convolution computation is updating the demagnetizing field only once at every p timesteps in the numerical integration of Eq. (1), where p is an integer. However, this procedure often results in the instability of numerical integration.
The demagnetizing field hst(i) in a cell can be divided into the two components h stin (i) and h stex (i) that are generated by the own cell i and by the other cells, respectively. They are given as The components of s(0) are often large. This is why the demagnetizing field should be updated at every time-step. Compared with s(0), s(i) (i  0) is relatively small, which implies that the h stex may be updated less frequently than h stin . Consequently, it is reasonable to update h stin at every time-step and h stex at every p ( 2) time-steps in the numerical integration.

Computational results
A magnetic material having dimensions l x  l y  l z is analyzed with the ADSM, where l x : l y : l z is set to 2: 1: 0.01. The material is divided into 32  16  1 cells. A cubic crystalline anisotropy is assumed with κ = 2K / (μ 0 M S 2 ) = 0.01. Figure 2 shows the relation between the execution period p for the convolution and the computation time. When p = 20, the computation time is reduced by 64 % compared with that with p = 1. When p > 50, the computation time slightly increases with p because the infrequent update of demagnetizing field deteriorates the convergence to an equilibrium point. Figure 3 shows the MH curves along the <110> axis set along the y-direction, which are obtained with p = 1 and 20. The numerical integration is executed by the forward Euler scheme. The property obtained with p = 20 coincides with that with p = 1.

Direct solution using unit cell property
If the magnetization property of unit cell is known, the macroscopic magnetization property can be synthesized from the magnetizations of unit cells. Figure 4 shows an example of the magnetization property of SDSM along the easy-axis direction. The magnetization state is classified into three types as below: S + : the single domain state with positive magnetization, S − : the single domain state with negative magnetization, WM: the state of 180 domain-wall motion.
The three magnetization states above exist within the respective intervals of the normalized applied field h as below (see Appendix A.3): When the SDSM is used as a cell in the macroscopic model, by including h stex into the applied field h as the magnetization of each cell can be determined using the unit cell property, where h eff (i) is called effective field. When h eff (i) moves outside the interval of present magnetization state, the magnetization state transition occurs in cell i. The alternating magnetization property is obtained by changing the applied field h step by step to find corresponding equilibrium point as follows.
Among all the cells where h eff (i) moves outside the present interval, the cell is chosen where h eff (i) is the most distant from the present interval. The magnetization state in the chosen cell is changed to another state so as for h eff (i) to be included in the corresponding interval as shown in Fig. 4. When h S < −h WM , the transition to the WM state does not occur as shown in Fig. 4(b). The magnetizations are also corrected in the cells having the 180 domain-wall motion state in accordance with the change of h eff (i) so as to satisfy Eq. (A.15). After the state transition and the magnetization correction, the demagnetizing field is recalculated and the procedure above is repeated until the demagnetizing field converges. The flowchart of above procedure is shown in Fig. 5.
It is basically possible for the direct solution method to be applied to the 2D or 3D analysis if the 2D or 3D cell property is known. However, the classification of magnetization states and the transition process among them are complex in the 2D/3D analysis and are not derived in a straightforward manner generally.

Computational results
The cells aligned one-dimensionally as in Fig. 6 Figure 7 shows the magnetization curves obtained by the original ADSM solving Eq. (1) using 1, 8 and 128 cells, whereas the properties shown in Fig. 8 is given by the direct solution. The simulation is started from the demagnetization state (see Appendix A.4).
The magnetization curve obtained by the direct method coincides with that obtained by the original ADSM in the case of single cell. In the case of 128 cells, however, the coercive force given by the direct solution is smaller than that given by the original ADSM. This suggests that the ADSM solving Eq. (1) sometimes fails to judge the convergence to an equilibrium point correctly and stops the time-integration of Eq. (1) incorrectly before escape from an unstable equilibrium point becomes evident. This is because the escaping process often requires very long transient time. In fact, the coercive force obtained by the original ADSM is not very robust and is affected by the computational condition such as the convergence criterion and the amplitude of h. Figure 9 shows the distributions of m x (i) and h eff (i) / H S (i = 1, …, 9) given by the original ADSM using 128 cells when h = 0 → −0.2. Even though h eff (5) is slightly smaller than H S (h eff (5) / H S < 1), cell 5 stays in the single domain state (m x (5) = 1) when h = −0.2. Figure 10 shows those distribution obtained by the direct solution method. The cells are in the positive and negative single domain states when 1 ≤ h eff (i) / H S and h / H S ≤ − H WM / H S (= −1.28), respectively. As a result, the ADSM with the direct solution reconstructs the very small coercive force compared with the anisotropy field that is often observed in soft magnetic materials. For example, the coercive force of oriented silicon steel is generally less than 1/100 of anisotropy field of Fe, which is not easy to be predicted by the micromagnetic simulation. However, the mechanism of small coercivity should be discussed further in future study.

Conclusion
First, this article presents an efficient method for the demagnetizing field computation using the decomposition into near and far fields. It is also expected that the decomposition allows the near filed to be integrated by an implicit scheme. Second, the search of an equilibrium point is greatly accelerated by the direct solution method using the magnetization property of unit cell. If the unit cell property is unknown, the magnetization state transition should be switched based on the bifurcation point detection as was discussed in Ref. [7].

A.1. Simplified domain structure model
The SDSM locally minimizes the total magnetic energy, e, normalized by the crystalline anisotropy energy [5]; e is given as where e ap is the Zeeman energy, e an is the crystalline anisotropic energy, e w is the domain-wall energy, and e st is the magnetostatic energy. The normalized crystalline anisotropic energy is given as ) , where f an represents the angular dependence and λ is the volume ratio of domain where w = 4l k /D, l k = (A/K) 1/2 is the characteristic length relevant to the exchange energy, A is the exchange stiffness constant and D is the width of the two domains; w represents the influence of energy cost to have domain walls. It is assumed that the demagnetizing field is uniform in the SDSM and that it is approximated as the multiplication of demagnetizing factors and the average magnetization. The magnetostatic energy is given as where (m x , m y , m z ) = m = λm 1 + (1λ)m 2 , s x = k x /κ, s y = k y /κ and s z = k z /κ; k x , k y , and k z are the demagnetizing factors. The effect of domain shape becomes abstracted by the approximation of magnetostatic energy. For example, the domain size affects only D and the direction of domain wall is not taken into account ( Fig. 1 does not reflect a real domain structure).
The total energy e becomes a local extremum when Eq. (1) is satisfied. Its solution gives a local minimum for e when all the eigenvalues of  2 e/X 2 are positive.

A.2. Assembly of domain structure models
To obtain the energy components of the ADSM, the Zeeman energy, the anisotropic energy and the domain-wall energy of cells are independently summed up. The normalized magnetostatic energy e st is computed as follows.
The demagnetizing field H st in the ADSM is obtained in the same way as in the micromagnetic simulation [6]; H st at cell i is given as ) ( ) ( ) ( where i and i' are cell indexes and N is the demagnetizing coefficient matrix [3]. For example, the cell index i is ordered as i = n x n y (K1) + n y (J1) + I ( I = 1, …, n x , J = 1,…, n y , K = 1,…, n z ) (A.7) where (I, J, K) and (n x , n y , n z ) are the cell indexes and the numbers of cells along the x-, y-and z-directions, respectively. Using (I, J, K), the demagnetizing coefficient matrix N(i) =