A Phase Prediction Method for Pattern Formation in Time-Dependent Ginzburg-Landau Dynamics for Kinetic Ising Model without a priori Assumptions on Domain Patterns
Ryoji Anzaki, Shin-ichi Ito, Hiromichi Nagao, Masaichiro Mizumaki, Masato Okada, Ichiro Akai
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p A Phase Prediction Method for Pattern Formation inTime-Dependent Ginzburg-Landau Dynamics for Kinetic Ising Modelwithout a priori
Assumptions on Domain Patterns
Ryoji Anzaki, Shin-ichi Ito,
1, 2
Hiromichi Nagao,
1, 2
Masaichiro Mizumaki, Masato Okada, and Ichiro Akai Earthquake Research Institute, The University of Tokyo,1-1-1, Yayoi, Bunkyo-ku, Tokyo 113-0032, Japan Graduate School of Information Science and Technology,The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, 113-8656 Tokyo, Japan Japan Synchrotron Radiation Research Institute (JASRI/SPring-8), 1-1-1 Kouto, Sayo, Hyogo 679-5198, Japan Graduate School of Frontier Sciences, The University of Tokyo, Kashiwa, Chiba 277-8561, Japan Institute of Pulsed Power Science, Kumamoto University, 39-1 Kurokami 2, Kumamoto 860-8585, Japan (Dated: September 2, 2020)We propose a phase prediction method for the pattern formation in the uniaxial two-dimensionalkinetic Ising model with the dipole-dipole interactions under the time-dependent Ginzburg-Landaudynamics. Taking the effects of the material thickness into account by assuming the uniformnessalong the magnetization axis, the model corresponds to thin magnetic materials with long-rangerepulsive interactions. We propose a new theoretical basis to understand the effects of the materialparameters on the formation of the magnetic domain patterns in terms of the equation of balance governing the balance between the linear- and nonlinear forces in the equilibrium state. Based onthis theoretical basis, we propose a new method to predict the phase in the equilibrium state reachedafter the time-evolution under the dynamics with a given set of parameters, by approximating thethird-order term using the restricted phase-space approximation [R. Anzaki, K. Fukushima, Y.Hidaka, and T. Oka,
Ann. Phys. , 107 (2015)] for the φ -models. Although the proposedmethod does not have the perfect concordance with the actual numerical results, it has no arbitraryparameters and functions to tune the prediction. In other words, it is a method with no a priori assumptions on domain patterns. I. INTRODUCTION
Magnetic materials are of great interest even beforethe beginning of the application of quantum physics tothe solid-state physics . The domain patterns are essen-tial for understanding the magnetic materials since themacroscopic properties of magnetic materials are largelyaffected by the domain patterns . In the light of recentprogress in experimental methods to observe the mag-netic domain patterns, it is now convincing that one mayobtain information on the magnetic dynamics, e.g., thematerial parameters, the external magnetic fields, andthe size of the magnetic materials from the domain pat-terns. Among such experimental methods, X-ray mag-netic circular dichroism (XMCD) and the detection ofthe Kerr effect via visible light are well-known meth-ods to detect the magnetization normal to the surfaceof the materials. In the realm of theoretical- and sim-ulation physics, researchers have already made progresstowards this aim. Jagra and Kudo et al . performednumerical simulations using similar models to reproducethe magnetic domain patterns on two-dimensional mag-netic materials. The latter proposed a relation betweenthe sweep rate of the external magnetic field and the fi-nal magnetic domain patterns in the equilibrium state.They utilized the two-dimensional kinetic Ising spin sys-tem with spins on the square lattice lying on the xy -plane, while the magnetization is restricted in the z -direction, which is normal to the xy -plane. Assumingthat the high-wavenumber components of the Green’s function of the dipole-dipole interaction play few roles,they succeeded in explaining the various domain pat-terns resulting from different sweep rates by solvingthe time-dependent Ginzburg-Landau (TDGL) equationnumerically . Iwano et al . adopted a numerically eval-uated effective two-dimensional Green’s function for thedipole-dipole interaction.In the early history of the researches of the TDGLdynamics, the probability density functions (PDFs) ofthe spin systems under the TDGL dynamics have beenstudied by Kawasaki in the 1970s. Suzuki et al . also studied the same system using the Markov chain.Their major interests were to obtain the global charac-teristics of the spin configuration, e.g., the dynamic mag-netic susceptibility and the critical exponents , usingthe analytical tools including the diagrammatic methods.In the 1980s, Grant et al . investigated the similar sys-tem in a context of the phase separation, and developeda theory using the spatial wavenumbers of the fields. Onthe other hand, Kawasaki also investigated the kinkdynamics in the one-dimensional TDGL model, whoseachievements have been inherited to the researches onthe dynamical phase-transition in the TDGL dynamicsof the XY-model . In the late 1980s, numerical sim-ulations have been performed using the TDGL equationin the real-space .In the realm of magnetic materials, the explana-tions of the magnetic domain patterns have been devel-oped for decades . The Kooy-Enz model and itsvariants assume simple domain patterns specified byfunctions with one or more parameters and minimize thetotal energy (the sum of the contributions from the do-main and the domain wall) with respect to the param-eters. The forms of the functions that determine thedomain patterns are chosen a priori so that the entireproblem simplifies into an optimization problem of realfunctions. Garel et al. analyzed the behavior of thesimilar system under finite temperature T and externalmagnetic field H thermodynamically and plotted the T - H phase diagram with three phases named uniform, bub-ble , and striped . These phases are also defined by simpleanalytic functions with a few parameters.In this Paper, we take a new strategy that does notinvolve any a priori defined functions. The effects of thematerial thickness and other parameters to the TDGLpattern formation are explained by the newly proposed equation of balance that describes the balances betweenthe linear- and nonlinear forces in the equilibrium statereached after appropriate numerical time-evolutions witha realistic initial condition. This equation enables us topredict the phase that a specific TDGL equation with agiven set of parameters forms in the equilibrium state.In the language of the magnetism, we can predict thepattern of the magnetic domain formed in thin magneticmaterials for a given set of the TDGL parameters withthe proposed method.In this Paper, we use a numerical method to constructan effective two-dimensional Green’s function by analyti-cally averaging the dipole-dipole interactions along the z -direction for each grid point on the xy -plane as proposedin Ref. , enabling us to take the effects of the thicknessinto account more precisely. II. MODEL AND METHODS
We utilize the Ising-like spin model with TDGLdynamics , also referred as the kinetic Isingmodel . We prepare an array of complex variables { φ ( r ) } with r being an element of two-dimensional dis-crete space D = { ( x, y ) : 1 ≤ x, y ≤ L } ⊂ Z for apositive integer L . Each variable φ is regarded as a mag-netic dipole restricted in the z -direction, while the vector r represents a coordinate on the xy -plane, normal to the z axis. Note that x - and y -components of spins are setto zero in this model. Introducing the saturation magne-tization ρ >
0, the TDGL equation for the spin systemabove with time parameter t is,d φ ( r )d t = W ( r | φ ] + B ( t ) , (1)where B is the explicitly time-dependent external mag-netic field (restricted in the z -direction), and W ( r | φ ] isa function of r and a functional of φ , defined as W ( r | φ ] = α [ φ ( r ) − ρ − φ ( r )] + β ∇ φ ( r ) − γF [ φ ] , (2) F [ φ ] = Z d r ′ G ( r − r ′ ) φ ( r ′ ) . (3) The terms containing α , β , and γ correspond to theanisotropy-, exchange- and the dipole-dipole interactions,respectively. The last term is represented via the Green’sfunction for the magnetic dipole-dipole interaction G ( − ).By moving into the wavenumber space by the (non-unitary) Fourier transform h f i k = L − X r ∈D f ( r )e i k · r , (4)the TDGL above becomes,d φ k d t = W k [ φ ] + B k ( t ) , (5)with W k [ φ ] , φ k and B k being the spatial Fourier transfor-mation of W ( r | φ ] , φ ( r ) and B ( t ). Performing the Fouriertransform, one obtains W k [ φ ] = α (cid:28) φ − φ ρ (cid:29) k − β | k | φ k − γL · G k φ k . (6)Here, the Fourier transformation of G is introduced viathe convolution theorem, and the prefactor L is due tothe choice of the Fourier transform Eq. (4).The effects of thickness are not apparent but intro-duced via the Fourier transform of the Green’s functionof the dipole-dipole interaction G − as already performedin . Hereafter, we assume that spin variables have thesame value along with the z -direction for each r . Thethickness (the spatial extension along the z -direction) ofthe material is assumed to take a positive value A > z -coordinate 0 ≤ z ≤ A , wedefine an effective two-dimensional Green’s function un-der the conditions specified above, as G ( r ; A ) = 1 A Z A d z Z A d z ′ G ( r ; z − z ′ ) , (7)with, G ( r ; ∆ z ) = 1( | r | + ∆ z ) / − z ( | r | + ∆ z ) / . (8)The integral Eq. (7) can be performed analytically, and G ( r ; A ) = 2 A | r | − p | r | + A ! . (9)Note that in the limit A → G ( r ) converges to theinverse-cubic law point-wise.One may consider the continuum limit , which corre-sponds to the case when the correlation length measuredin the unit of the grid spacing becomes positive infinity.In that case, the Fourier transform of the Green’s func-tion is obtained from the real-space function G ( − ) andhas the analytical form G k ( A ) = 1 πA − e − A | k | | k | . (10)The weight of the Fourier transformation is taken tobe (2 π ) − d , where d = 2 is the spatial dimension. Inthis limit, the right-hand side of the equation of motionEq. (6) becomes W k [ φ ] = α (cid:28) φ − φ ρ (cid:29) k − β | k | φ k − γ (2 π ) G k ( A ) φ k . (11)Note that this representation is formally obtained simplyby a replacement L → π . III. NUMERICAL SIMULATIONS
In the xy -plane, we use the non-unitary fast Fouriertransform (FFT) corresponding to Eq. (4) to constructthe modes φ k = h φ i k and the wavenumber representa-tion of the Green’s function G k . We adopt the periodicboundary condition for x - and y -direction, hence the en-tire topology of the simulation space is a torus. Thespacing of the grid on the xy -plane is set to unity. We in-troduce a randomness of the coefficient of the anisotropyas α → α Λ( r ) with Λ( r ) = 1 + λ ( r ) / λ ( r ) ∼ N (0 , .
3) independently and identically forall r , as in Ref. . The external magnetic field intensityis represented by the rectified linear unit function R ( − )as B ( t ) = R ( B − v B t ) with B , v B ≥ φ ( r ) randomly in a range − . ≤ φ ( r ) ≤ − φ k in thecubic term, are now circumvented by this method sim-ply by performing the algebraic operation φ ( r ) [ φ ( r )] for each r ∈ D . The time-evolutions are performed effi-ciently with the ETD2/RK4 method , one of the multi-step exponential integrator methods, with relatively largetime step δt = 0 .
2. The scalability to the system size L is quite good, with the computational time roughly pro-portional to L , up to the largest case considered here( L = 512). In Fig. 1, one can see qualitatively differentfinal magnetic domain patterns depending on differentvalues of the thickness of A .In this system we can define the two-dimensionaltranslational- and rotational symmetry in the coordinatespace and the Z symmetry of the spin. Thus the pat-terns above can be naturally classified into four phases according to these symmetries: Symmetric [Fig. 1 (a)],T-breaking [(b)], TZ-breaking [(c)], and Z-breaking [(d)]phases, with “T” standing for “translational and rota-tional” while “Z” standing for Z . Literatures e.g. use more descriptive terms referring (b) and (c), such as“labyrinth” and “sea-island”, respectively. (a) Symmetric phase ( A = 1 . - 50 - 25 0 25 50- 50- 2502550 x y (b) T-breaking phase ( A = 1 . - 50 - 25 0 25 50- 50- 2502550 x y - 0.50- 0.2500.250.50 0 100 200 300- 0.50.00.51.01.5 time MCB (c) TZ-breaking phase ( A = 2 . - 50 - 25 0 25 50- 50- 2502550 x y - 0.50- 0.2500.250.50 0 100 200 300- 0.50.00.51.01.5 time MCB (d) Z-breaking phase ( A = 3 . - 50 - 25 0 25 50- 50- 2502550 x y FIG. 1. Left Panels: The magnetic domain patterns at theend of the simulation for various A . Right Panels: The timedependencies of the external magnetic field (marked by B),the average magnetization h φ i (marked by M), and the cor-relation (cid:10) δφ (cid:11) (marked by C). Simulation size: L = 512, du-ration of the time-evolution: t max = 300, the external mag-netic field intensity: B ( t ) = R ( B − v B t ) with B = 1 . v B = 0 .
01 and R ( x ) = max( x, α = β = 1 , γ = 0 .
2. The white regions has positive magneti-zation, while the black region has negative magnetization.
IV. NORMALIZATION OF THE TDGLDYNAMICS
The physical- or dimension-full TDGL equation Eq. (1)is to be normalized by the linear temporal- and spatial co-ordinate transformations to compare with other results.One of the most convenient choices is to eliminate thedimensionfull saturation magnetization ρ . In this case,using the new time variable τ , spatial coordinate ζ , thelaplacian ∂ and magnetization ϕ , one finds the normal-ized TDGL equation for the dynamics under the externalmagnetic field swept from B > v B is,d ϕ d τ = ϕ − ϕ + ∂ ϕ − p F [ ϕ ] − R ( p − p τ ) , (12) F [ ϕ ] = Z d ζ ′ G ( ζ − ζ ′ ) ϕ ( ζ ′ ) . (13)Here the linear functional F denotes the dipole-dipoleinteraction, while the coefficients p i ( i = 1 , , ,
4) aredefined as, p = γ √ αβ , p = B √ ρα , p = v B ρα , p = r αβ A. (14)The last parameter p represents the normalized thick-ness and is used to construct the effective two-dimensional Green’s function as in Eq. (7). V. EQUATION OF BALANCE AND RPSA
If the spin configuration φ is in the equilibrium state,˙ φ k = 0 for all k , and the generic (whether it is normalizedor not) equation of motion [Eqs. (1,2)] simplifies into aset of simultaneous time-independent equations. We nowintroduce a new idea equation of balance (EOB), that isthe equation of motion in the equilibrium state with zeroexternal magnetization, as shown below. (cid:10) φ (cid:11) k = Q k h φ i k ; Q k = α − β | k | − γL G k α . (15)Let us rewrite Eq. (15) in the average magnetization h φ i and the modes h δφ i k with δφ ( r ) = φ ( r ) −h φ i . By notingthat h δφ i = 0 and h c i k = 0 for any constant c and k = 0, we immediately obtain the relations governingthe balance between the first-, second-, and the third-order moments of the field variables at the equilibriumstate. For k = 0, (cid:2) Q − (cid:10) δφ (cid:11) (cid:3) h φ i = (cid:10) δφ (cid:11) + h φ i , (16)and for k = 0, D k h δφ i k = (cid:10) δφ (cid:11) k + 3 (cid:10) δφ (cid:11) k h φ i , (17)with D k = Q k − h φ i . (18)These equations do not specify the equilibrium stateuniquely. This lack of uniqueness is obvious if one notesthat the entire dynamics led to the equilibrium state isnot included in the EOB. Thus the EOB must be under-stood as the restrictions that an equilibrium state mustsatisfy.Since the third-order moment in the EOB can hardlybe estimated, we apply the restricted phase-space ap-proximation (RPSA) to the equation above. In ourcurrent context, it is equivalent to a replacement (cid:10) δφ (cid:11) k → (cid:10) δφ (cid:11) h δφ i k ( k = 0) . (19) (20) In general, the RPSA truncates the interaction terms (cid:10) δφ (cid:11) k systematically, and known to be exact in specialmodels, e.g O ( N ) scalar model with N → ∞ .In the diagrammatic notation, the RPSA is a restric-tion of the convolution in the Ginzburg-Landau pseudofree energy [corresponding to the equation of motionEq.(6)] as shown below.In the diagram above, “perm.” indicates the permu-tations of the vertices. Noting that h δφ i = 0, it mustbe emphasized that the RPSA approximates (cid:10) δφ (cid:11) by0. Thus we obtain the RPSA-EOB as shown below. For k = 0, (cid:2) Q − (cid:10) δφ (cid:11) (cid:3) h φ i = h φ i . (21)For k = 0, h Q k − h φ i − (cid:10) δφ (cid:11) i h δφ i k = 3 (cid:10) δφ (cid:11) k h φ i . (22)We confirm that the RPSA agrees with the results of thetime-evolutions in the relatively small ( L = 128) systemwith α = 3 . , β = 2 . γ = 2 /π . The average magne-tization h φ i obtained via the RPSA equation of balanceusing the numerical results of Q and (cid:10) δφ (cid:11) matcheswith the simulation, except for the range 1 < A < . (cid:10) δφ (cid:11) has nonzero values. This is due to the factthat the RPSA neglects the third-order moment (cid:10) δφ (cid:11) of the distribution. VI. PHASE-PREDICTION METHOD BY THERPSA-EQUATION OF BALANCE
The RPSA-EOB described in the previous Section isapplicable for phase predictions of the TDGL dynamics.We simply predict the types of the patterns based onthe nonzero modes under the restrictions imposed by theRPSA-EOB [Eqs. (21,22)]. We use the continuum EOM[Eq. 11] to use the analytic form of the Green’s function[Eq. (10)]. Note that this choice causes a modification tothe definition of Q k : Q k = α − β | k | − γ (2 π ) G k ( A ) α , (23)with G k ( A ) being the continuum limit of theGreen’s function shown in Eq. (10). The method PhasePrediction is schematically shown below. This isa procedure that maps a set of parameters ( α, β, γ, A )to the output
Phase ∈ {
Symmetric, T-breaking,*Z-breaking } , with *Z-breaking means either TZ-breaking or Z-breaking. procedure PhasePrediction ( α, β, γ, A ) g ( k ) = (1 − e − Ak ) / ( πA k ) Q ( k ) = 1 − α − βk − α − γ (2 π ) g ( k ) if Q (0) < then if max( Q ) < then Phase ⇐ Symmetric else Phase ⇐ T-breaking end if else Phase ⇐ *Z-breaking end if end procedure VII. DISCUSSION
The phase diagram for the normalized TDGL dynam-ics [Eq. (12)] predicted by the method
PhasePrediction described in Sec. VI is shown in Fig. 2. The overall ten-dency matches our physical instinct well. As p becomeslarge, the demagnetization effect from the dipole-dipoleinteractions supersedes the anisotropy to yield the sym-metric phase, while for larger p , it is partly relaxed bythe thickness to have more complexed structures.We also compared the numerical results of time-evolution with the phase prediction. The results areshown in Table I. The computational cost for each sam-ple point ( p , p ) is significantly small compared to thecorresponding numerical time-evolutions.The agreement between the time-evolution and thephase prediction is good, except in the cases ( p , p ) =(0 . , . , (0 . , . , (0 . , . Q ) and Q (0) at these sam-ple points. Since the RPSA neglects the third-order mo-ments in the EOBs, the results obtained by the RPSA-EOB based method may differ from the time-evolutionfor small | Q (0) | and | max( Q ) | . This mismatch may im-prove by further developments of the approximation; inother words, it is considered that the third-order mo-ments play crucial roles in the region where the mismatchis seen.Note that the external magnetic sweep rate v B is animportant parameter in the pattern formation. It is FIG. 2. The phase diagram of the TDGL dynam-ics corresponding to Eq. (12) estimated by the method
PhasePrediction proposed in Sec. VI. Orange area (markedwith “*Z”): TZ- or Z-breaking phase, blue area (“T”): T-breaking phase, green area (“Symmetric”): symmetric phase.TABLE I. Comparison with numerical time-evolutions with
PhasePrediction . T, Z, TZ, and S denotes the phases ob-served from the time-evolution (T-, Z-, TZ-breaking, andsymmetric phases, respectively), while ( s, t ) with s, t ∈ { + , −} denotes the signs of the Q and max( Q ), respectively. Notethat PhasePrediction translates (++) to *Z (= Z, TZ), ( − +)to T, and ( −− ) to S. The numerical time-evolutions are per-formed for the system with 512 grid points using the nor-malized TDGL equation Eq.(12). Points marked with “ ⋆ ”indicate the mismatches with the numerical results. p \ p ⋆ ( −− ) S( −− ) S( −− ) S( −− ) S( −− )1.5 Z(++) T( −− ) S( −− ) S( −− ) S( −− )2.0 Z(++) TZ( − +) S( −− ) S( −− ) S( −− )2.5 Z(++) Z(++) T ⋆ ( −− ) S( −− ) S( −− )3.0 Z(++) Z(++) T ⋆ ( −− ) S( −− ) S( −− ) reported that the domain formation is largely affectedby v B . Our results here must be understood as an ap-proximated result, not only in the RPSA but also in theelimination of the effects of the magnetic sweep rate. Inthe v B → ∞ limit, our method will have the same results,while that of the time-evolutions can be quite different.Although this method does not have a perfect concor-dance with the numerical simulations, it has no a pri-ori parameters or functions in any form, but only ap-proximated in a systematic, physically reasonable way .Hence it is considered as a method without any a priori assumptions on the domain patterns. This fact meansthat one can add new features, e.g., tuning parameters,without doubting the physical meaning of this method,provided the approximation is reasonable. VIII. CONCLUSIONS
The long history of the research in the magnetism andthe mathematical structure of the TDGL dynamics showa wide variety of approaches to the pattern formation inthe magnetic materials .Although most of the existing methods use artificialfunctions that specify the magnetic domain patterns, wefocus on the equation of balance (EOB) that a magneticmaterial must satisfy in its equilibrium state. Apply-ing the restricted phase-space approximation (RPSA) to EOB enables us to predict the phase in the equilib-rium state. The prediction matches the actual numericaltime-evolution results qualitatively without any tuningparameters. Although the prediction is not perfect, ourmethod has no a priori assumptions, i.e., it does not in-volve any artificial function or experimentally justifiedparameters but only approximated systematically. Thusit is very extensive, applicable for various applications.Another aspect that must be noted is that the object ofthe new method is not limited to the magnetic systems; itis applicable for a vast class of natural/social phenomenathat seemingly have nothing in common but described by the equation of motion of type Eqs.(1,2).One of such applications is parameter estimation inmaterial- and statistical physics. Using the Bayesian in-ference methods, we can estimate the parameters of asystem with huge degrees of freedom by relatively smallobservation/numerical data, e.g., . Our method is ex-pected to serve for such parameter estimations in varioussystems as the theoretical- and numerical basis by givinginformation on the phase for each parameter using a fewcomputational costs, with physically justifiable reasons. ACKNOWLEDGEMENTS
This work was mainly supported by JST CRESTGrant Numbers JPMJCR1761 and JPMJCR1861 andpartially supported by JPMJCR1763 of Japan Scienceand Technology Agency. The key ideas in this studycame through the activities of JSPS KAKENHI GrantNumbers JP19K14671, JP17H01703, JP17H01704,JP18H03210, JP19H05662, and JP20K21785. The travelexpense needed to discuss among co-authors was par-tially supported by ERI JURP 2020-A-05, 2018-B-01,and 2019-B-04. J. H. Van Vleck, Rev. Mod. Phys. , 27 (1945). C. Kittel, Rev. Mod. Phys. , 541 (1949). M. Suzuki, N. Kawamura, M. Mizumaki, Y. Terada,T. Uruga, A. Fujiwara, H. Yamazaki, H. Yumoto,T. Koyama, Y. Senba, et al. , in
J. Phys. Conf. Ser. , Vol.430 (2013). P. N. Argyres, Phys. Rev. , 334 (1955). J. Reif, J. C. Zink, C.-M. Schneider, and J. Kirschner,Phys. Rev. Lett. , 2878 (1991). E. A. Jagla, Phys. Rev. E , 046204 (2004). K. Kudo and K. Nakamura, Phys. Rev. B , 054111(2007). K. Iwano, C. Mitsumata, and K. Ono, J. Appl. Phys. ,17D134 (2014). K. Kawasaki, Prog. Theor. Phys. , 359 (1974). K. Kawasaki, Prog. Theor. Phys. , 1064 (1974). K. Kawasaki, Prog. Theor. Phys. , 84 (1974). M. Suzuki and G. Igarashi, Prog. Theor. Phys. , 1070(1973). M. Grant, M. San Miguel, J. Vials, and J. D. Gunton,Phys. Rev. B , 3027 (1985). K. Kawasaki and T. Ohta, Physica A , 573 (1982). T. Yasui, H. Tutu, M. Yamamoto, and H. Fujisaka, Phys.Rev. E , 036123 (2002). N. Fujiwara, H. Tutu, and H. Fujisaka, Phys. Rev. E ,066132 (2004). T. M. Rogers, K. R. Elder, and R. C. Desai, Phys. Rev. B , 9638 (1988). C. Kittel, Phys. Rev. , 965 (1946). B. Kaplan and G. Gehring, J. Magn. Magn. Mater ,111 (1993). G. Bochi, H. J. Hug, D. I. Paul, B. Stiefel, A. Moser,I. Parashikov, H.-J. G¨untherodt, and R. C. O’Handley,Phys. Rev. Lett. , 1839 (1995). C. Kooy and U. Enz, Philips Research Reports (1960). A. Lisfi and J. C. Lodder, J. Phys. Condens. Matter ,12339 (2002). T. Garel and S. Doniach, Phys. Rev. B , 325 (1982). K. Kudo, M. Mino, and K. Nakamura, J. Phys. Soc. Japan. T. Yokota, J. Magn. Magn. Mater. , 532 (2017). T. Tom´e and M. J. de Oliveira, Phys. Rev. A , 4251(1990). S. Krogstad, J. Comput. Phys. , 72 (2005). C. B. Muratov, Phys. Rev. E , 066108 (2002). R. Anzaki, K. Fukushima, Y. Hidaka, and T. Oka, Ann.Phys. , 107 (2015). S. Ito, H. Nagao, T. Kurokawa, T. Kasuya, and J. Inoue,Phys. Rev. Mater.3