Introduction

Soft magnetic materials (SMMs), such as Fe-Si electrical steels and nanocrystalline soft magnetic materials (NSMMs), are key components in many electrical devices that rely on magnetic flux, including inductors, transformers, motors, and generators. These devices are indispensable for the generation and conversion of energy in devices used in daily life. The core loss, which is generated by retardation of magnetization dynamics to an external magnetic field, reduces the energy efficiency of these devices. To improve the efficiency of SMMs and related devices, it is necessary to reduce the core loss while maintaining the saturation magnetization as high as possible1,2,3. Fe-Si electrical steels are a traditional family of SMMs that have been used for more than a century4. Fe-3wt%Si exhibits a high saturation magnetization of 2 T, higher than those of conventional NSMMs, also referred to as metal amorphous nanocomposites by some authors5. However, the coercivity of nonoriented Fe-3wt%Si steel is typically 50 A/m2, and the core loss is considerably larger than that of NSMMs. Thus, despite the high-saturation magnetization, the upper limit of operational induction is usually lower than 2 T for Fe-Si steel6. In the late 1980s and early 1990s, it was demonstrated that a low coercivity on the order of 1 A/m could be realized by reducing the crystallite size to 10 nm. The saturation magnetization in these NSMMs ranged from 1.3 to 1.7 T, lower than that of Fe-Si steels7,8,9,10,11,12. Most recently, it has been enhanced to 2 T by realizing magnetically soft nanostructures with a mass fraction of Fe comparable to that of Fe-Si steels13,14,15. This has made NSMMs promising candidates for improving the efficiency of electrical devices.

High-frequency magnetic fields are applied to SMMs in general applications, such as motors and generators, in which the energy loss in the SMMs is dominated by classical eddy current loss and excess loss; the excess loss comprises losses other than the hysteresis loss and classical eddy current loss16,17,18,19,20. Since NSMMs are prepared in rapidly solidified ribbon forms with large electrical resistances due to their thin thickness, the classical eddy current loss is greatly suppressed21. In contrast, the excess losses could remain high when the coherence of the domain wall displacement is affected at high frequencies. Furthermore, it has been noted that the excess losses show a tendency to increase with increasing saturation magnetostriction, suggesting that magnetostriction may play an important role in determining the domain wall motion.

Magnetostriction is a strain on the crystal lattice induced by magnetization22. It is well known that excellent soft magnetic properties are obtained when both the saturation magnetostriction and the magnetocrystalline anisotropy (K1) are small. The effective K1 in NSMMs is reduced considerably by the exchange interaction. This exchange softening effect has been described well by the random anisotropy theory23,24,25. The coercive field of NSMMs is described by extending this theory to a two-phase model, which has a residual amorphous phase26,27,28. Recently, micromagnetic simulation has reproduced the grain-size dependency of the coercive field29. Similar to magnetocrystalline anisotropy, the mean magnetostriction can be reduced by exchange coupling between randomly oriented nanocrystals30. This is another advantage of NSMMs. However, the local magnetostriction persists in each nanocrystallite, thereby disturbing the motion of the domain walls even if the nanostructure is designed to exhibit negligible mean magnetostriction. Nevertheless, little attention has been given to the effects of local and mean magnetostriction on the excess loss in NSMMs, even though the experimental data have suggested a correlation between them12,31,32. There is no simulation that clarifies the distribution of magnetostriction in NSMMs because the random orientations of crystal axes complicate calculations for the lattice dynamics. It is therefore important to formulate a micromagnetic simulation model of NSMMs that can address the local magnetostriction explicitly.

In this work, we formulate a micromagnetic simulation model of NSMMs and simulate the excess loss due to magnetostriction under an oscillating magnetic field. This model specifically formulates the effective field of magnetostriction in NSMMs where the crystal axis of each nanocrystallite is randomly distributed. The magnetostriction in each nanocrystallite is calculated in a local frame whose principal axis is parallel to the crystal axis. Then, the effective field due to magnetostriction is calculated by summing the contributions from all nanocrystallites in the absolute frame, and the dynamics of magnetization are simulated by solving the Landau–Lifshitz–Gilbert (LLG) equation together with the equation of motion for lattice deformation. Our simulator enables us to determine the conditions under which the mean magnetostriction disappears in an NSMM thin sheet.

We have found that the strain is nonuniformly distributed in NSMMs, and adjacent nanocrystallites influence each other through the stress. The domain wall motion is delayed by both the mean and local magnetostriction with respect to the oscillating magnetic field. This retardation effect of the wall motion due to magnetostriction results in excess loss. The results of this work can enable the magnetization dynamics in NSMMs to be simulated and may help identify ideal nanostructures for reducing the core loss in nanocrystalline soft magnetic alloys.

Models and methods

Formulation of a micromagnetic simulation model for NSMMs

To simulate the effect of magnetostriction on the technical magnetization process in NSMMs, it is necessary to calculate both the magnetization dynamics and the lattice dynamics due to magnetostriction33,34,35,36,37,38. The dynamics of magnetization are obtained by solving the LLG equation described by

$$\frac{{dm_i}}{{dt}} = - \gamma \varepsilon _{ijk}m_jH_k^{{{{\mathrm{eff}}}}} + \alpha \varepsilon _{ijk}m_j\frac{{dm_k}}{{dt}},$$
(1)

where γ is the gyromagnetic ratio, α is the Gilbert damping constant, εijk is the Levi–Civita tensor, mi is the i-th component of the magnetization unit vector, and \(H_i^{{{{\mathrm{eff}}}}}\) is the effective magnetic field. Hereinafter, the subscripts for examples i, j, and k are 1, 2, and 3, which represent the components in x1, x2, and x3, respectively. The effective magnetic field comprises an external magnetic field, an exchange field, a dipole field, and an effective field due to magnetostriction. The exchange and dipole fields are given by

$$\begin{array}{*{20}{l}}H_i^{{{{\mathrm{exc}}}}}\left( {{{\boldsymbol{x}}}} \right) = \frac{{2A}}{{\mu _0M_s}}\nabla ^2m_i\left( {{{\boldsymbol{x}}}} \right), \\ H_i^{{{{\mathrm{dip}}}}}\left( {{{\boldsymbol{x}}}} \right) = M_s{\int} {d\boldsymbol{x}^\prime K_{ij}\left( {{{{\boldsymbol{x}}}} - {{{\boldsymbol{x}}}}^\prime } \right)m_j\left( {{{{\boldsymbol{x}}}}^\prime } \right),}\end{array}$$
(2)

where x is the position vector, μ0 is the vacuum permeability, Ms is the saturation magnetization, A is the exchange constant, and Kij(x−x′) is the demagnetization tensor39,40,41,42,43. The effects of crystal lattice deformation on the exchange and dipole fields were neglected because the strain of the crystal lattice is as small as ~10−5 22. Since the aim of this report is to clarify the effect of magnetostriction on the dynamic technical magnetization process, we have purposefully left out other well-known influential parameters on the magnetization process. It is known that the magnetocrystalline anisotropy and eddy current affect the domain wall motion and energy loss in SMMs. Furthermore, there is an interplay between the magnetocrystalline anisotropy and the magnetostriction. However, we have focused only on the effects of magnetostriction using micromagnetic simulation. These conditions will unequivocally uncover the effect of magnetostriction on the technical magnetization process.

The effective field due to magnetostriction is given by the derivative of the total mechanical elastic energy density ε in terms of the magnetization, i.e.,

$$H_i^{{{{\mathrm{ms}}}}} = - \frac{1}{{\mu _0M_s}}\frac{{\partial {{{\mathcal{E}}}}}}{{\partial m_i}},$$
(3)

where ε is the sum of the magnetoelastic energy density εmagel and the elastic energy density εel. Because the crystal axis of each nanocrystallite in an NSMM is randomly distributed, it is necessary to introduce the rotation matrix Rij that aligns the crystal frame (the coordinate system established by the crystal axes) with the absolute frame to calculate ε.

Figure 1 shows the step-by-step transformation from the crystal frame to the absolute frame upon sequential application of rotation matrices. The rotation matrix that aligns the crystal frame with the absolute frame is given by

$$R_{ij} = R_{ik}^\psi R_{kl}^\theta R_{lj}^\phi = \left( {\begin{array}{*{20}{c}} {\cos \psi \cos \theta \cos \phi - \sin \psi \sin \phi } & {\cos \psi \cos \theta \sin \phi + \sin \psi \cos \phi } & { - \cos \psi \sin \theta } \\ { - \sin \psi \cos \theta \cos \phi - \cos \psi \sin \phi } & { - \sin \psi \cos \theta \sin \phi +\cos \psi\cos \phi} & {\sin \psi \sin \theta } \\ {\sin \theta \cos \phi } & {\sin \theta \sin \phi } & {\cos \theta } \end{array}} \right),$$
(4)

where the definitions of the rotational angles {θ, ϕ,ψ} and the rotation matrices \(\left\{ {R_{ij}^\theta ,R_{ij}^\phi ,R_{ij}^\psi } \right\}\) are illustrated in Fig. 1. Throughout this paper, we use the summation convention with repeated indices.

Fig. 1: Relation between the crystal frame and absolute frame.
figure 1

The crystal frame is represented by three rotational angles, θ, ϕ, and ψ. The axes of the absolute frame are represented by x1, x2, and x3. The sequential application of the rotation operators \(R_{ij}^\phi\), \(R_{ij}^\theta\), and \(R_{ij}^\psi\) aligns the crystal frame with the absolute frame.

Full size image

Given that the nanocrystallites in the vast majority of NSMMs are bcc-Fe, cubic symmetry is assumed for the elastic energy density. The magnetoelastic energy density is given by

$$\begin{array}{ll}{{{\mathcal{E}}}}_{{{{\mathrm{magel}}}}} = \left( {c_{12}^\prime - c_{11}^\prime } \right)\left( {\xi _{11}^\prime e_{11}^\prime + \xi _{22}^\prime e_{22}^\prime + \xi _{33}^\prime e_{33}^\prime } \right) \\\qquad\quad- 4c_{44}^\prime \left( {\xi _{12}^\prime e_{12}^\prime + \xi _{23}^\prime e_{23}^\prime + \xi _{31}^\prime e_{31}^\prime } \right),\end{array}$$
(5)

where \(c_{11}^\prime\), \(c_{12}^\prime\) and \(c_{44}^\prime\) are elastic constants for the cubic symmetry materials in matrix notation and \(\xi _{ij}^\prime\) is the strain35. Hereinafter, the prime symbol denotes the variables in the crystal frame. The stress-free strain \(e_{ij}^\prime\) is defined as

$$e_{ij}^\prime = \left\{ {\begin{array}{*{20}{c}} {\frac{3}{2}\lambda _{100}^\prime \left( {m_i^{\prime 2} - \frac{1}{3}} \right)\left( {i = j} \right)} \\ {\frac{3}{2}\lambda _{111}^\prime m_i^\prime m_j^\prime \left( {i \,\ne\, j} \right),} \end{array}} \right.$$
(6)

where \(m_i^\prime = R_{ij}m_j\) is the magnetization unit vector and \(\lambda _{100}^\prime\) and \(\lambda _{111}^\prime\) are magnetostriction constants22. When the magnetostriction constants have the same value, the magnetostriction of the single crystal of the magnetic material is isotropic. The elastic energy density due to strain is defined as34

$$\begin{array}{ll}{{{\mathcal{E}}}}_{{{{\mathrm{el}}}}} = \frac{1}{2}c_{11}^\prime \left( {\xi _{11}^{\prime 2} + \xi _{22}^{\prime 2} + \xi _{33}^{\prime 2}} \right) + 2c_{44}^\prime \left( {\xi _{12}^{\prime 2} + \xi _{23}^{\prime 2} + \xi _{31}^{\prime 2}} \right) \\\qquad+\, c_{12}^\prime \left( {\xi _{22}^\prime \xi _{33}^\prime + \xi _{33}^\prime \xi _{11}^\prime + \xi _{11}^\prime \xi _{22}^\prime } \right).\end{array}$$
(7)

After a straightforward calculation, the total mechanical elastic energy ε = εmagel + εel in tensor notation is rewritten as

$${{{\mathcal{E}}}} = \frac{1}{2}c_{ijkl}\left( {\xi _{ij} - e_{ij}} \right)\left( {\xi _{kl} - e_{kl}} \right) - {{{\mathcal{E}}}}_0,$$
(8)

where cijkl is the elastic constant tensor in the absolute frame, \(\xi _{ij} = R_{ik}^t\xi _{kl}^\prime R_{lj}\), \(e_{ij} = R_{ik}^te_{kl}^\prime R_{lj}\), and \({{{\mathcal{E}}}}_0 = 0.5c_{ijkl}e_{ij}e_{kl}\). The superscript t represents the transpose of a matrix. The elastic constant tensor in the absolute frame is denoted by

$$c_{ijkl} = R_{io}^tR_{lm}^tc_{opnm}^\prime R_{pj}R_{nk},$$
(9)

where \(c_{ijkl}^\prime\) is the elastic constant tensor in the crystal frame constructed by \(c_{11}^\prime\), \(c_{12}^\prime\), and \(c_{44}^\prime\)35. It is necessary to use the relation of the stress-free strain, i.e., \(e_{11} + e_{22} + e_{33} = 0\), to derive Eq. (8), as shown in the supplemental file. We neglect the small correction term ε0 that is related to crystal anisotropy. Thus, in the crystal frame, the effective field due to magnetostriction in the rotated x1-direction is written as

$$\begin{array}{ll}H_1^{\prime {{{\mathrm{ms}}}}} = \frac{3}{{\mu _0M_s}}\lambda _{100}^\prime \left( {c_{11}^\prime - c_{12}^\prime } \right)\left( {\xi _{11}^\prime - e_{11}^\prime } \right)m_1^\prime \\\qquad\quad+ \frac{6}{{\mu _0M_s}}\lambda _{111}^\prime c_{44}^\prime \left[ {\left( {\xi _{12}^\prime - e_{12}^\prime } \right)m_2^\prime + \left( {\xi _{31}^\prime - e_{31}^\prime } \right)m_3^\prime } \right],\end{array}$$
(10)

where \(\xi _{ij}^\prime\) is the total strain in the crystal frame. We assumed that the changes in cell volume are negligible, i.e., \(\xi _{11} + \xi _{22} + \xi _{33} = 0\), to obtain Eq. (10). The magnetostrictive magnetic field components in the rotated x2- and x3-directions are obtained through permutations of indices 1, 2, and 3. After calculating the rotated effective field \(R_{ij}^tH_j^{\prime {{{\mathrm{ms}}}}}\), we obtain the magnetostrictive field in the absolute frame \(H_i^{{{{\mathrm{ms}}}}}\).

To calculate the effective field due to magnetostriction, we have to obtain the total strain \(\xi _{ij}^\prime = R_{ik}\xi _{kl}R_{lj}^t\). Accordingly, we decompose ξij into the homogeneous strain \(\bar \varepsilon _{ij}\) and heterogeneous strain ηij as44

$$\xi _{ij} = \bar \varepsilon _{ij} + \eta _{ij}.$$
(11)

Hereinafter, the overbar symbol represents the spatial mean. Note that \(\bar \varepsilon _{ij}\) is independent of the position xi. The heterogeneous strain satisfies the condition such that the integration of ηij over the whole system disappears, i.e.,

$${\int} {d{{{\boldsymbol{x}}}}\eta _{ij} = 0.}$$
(12)

Assuming that the strain dynamics are much faster than the magnetization dynamics because the frequency of the external magnetic field is smaller than the eigenfrequency of the magnetic materials, \(\bar \varepsilon _{ij}\) is obtained by using the energy stability condition34. After performing the partial differentiation of the total elastic energy \(E = {\int} {d{{{\boldsymbol{x}}}}{{{\mathcal{E}}}}}\) with respect to \(\bar \varepsilon _{ij}\), we obtain

$$\frac{{\partial E}}{{\partial \bar \varepsilon _{ij}}} = V\left[ {\bar c_{ijkl}\bar \varepsilon _{kl} + \bar \sigma _{ij}^{{{{\mathrm{het}}}}} - \bar \sigma _{ij}^{{{\mathrm{m}}}}} \right],$$
(13)

where V is the system volume, \(\bar c_{ijkl}\) is the mean elastic constant, \(\bar \sigma _{ij}^{{{{\mathrm{het}}}}}\) is the mean stress due to the heterogeneous strain, and \(\bar \sigma _{ij}^{{{\mathrm{m}}}}\) is the stress-free strain. After calculating the magnetization unit vector Rijmj, the stress-free strain is obtained by Eq. (6). Then, we calculate the stress using \(c_{11}^\prime\), \(c_{12}^\prime\), and \(c_{44}^\prime\) in the crystal frame. The stress due to the heterogeneous strain is similarly calculated. Thus, these stresses are rewritten as

$$\bar \sigma _{ij}^{{{{\mathrm{het}}}}} = \frac{1}{V}{\int} {d{{{\boldsymbol{x}}}}R_{in}^tc_{nmkl}^\prime \eta _{kl}^\prime R_{mj},\quad \bar \sigma _{ij}^{{{\mathrm{m}}}} = \frac{1}{V}{\int} {d{{{\boldsymbol{x}}}}R_{in}^tc_{nmkl}^\prime e_{kl}^\prime R_{mj},} }$$
(14)

where \(\eta _{ij}^\prime = R_{ik}\eta _{kl}R_{lj}^t\) is the heterogeneous strain in the crystal frame. Under the energy stability condition, the homogeneous strain is obtained by

$$\bar \varepsilon _{ij} = \bar s_{ijkl}\left( {\bar \sigma _{kl}^{{{\mathrm{m}}}} - \bar \sigma _{kl}^{{{{\mathrm{het}}}}}} \right),$$
(15)

where \(\bar s_{ijkl} = \bar c_{ijkl}^{ - 1}\) is the mean elastic compliance constant. The homogeneous strain is independent of the heterogeneous strain in a single crystal. In contrast, in an NSMM, the homogeneous strain is dependent on the heterogeneous strain; therefore, we have to consider the contributions from the heterogeneous strain to calculate the homogeneous strain. Introducing the displacement of the crystal lattice ui, the heterogeneous strain is expressed as44

$$\eta _{ij} = \frac{1}{2}\left[ {\frac{{\partial u_i}}{{\partial x_j}} + \frac{{\partial u_j}}{{\partial x_i}}} \right].$$
(16)

The time evolution of the displacement is determined by the equation of motion:

$$\rho \frac{{d^2u_i}}{{dt^2}} = \frac{\partial }{{\partial x_j}}\left[ {\sigma _{ij}^{{{{\mathrm{hom}}}}} + \sigma _{ij}^{{{{\mathrm{het}}}}} - \sigma _{ij}^{{{\mathrm{m}}}}} \right] - \zeta \frac{{du_i}}{{dt}},$$
(17)

where ρ is the mass density, \(\zeta\) is the damping constant, \(\sigma _{ij}^{{{{\mathrm{hom}}}}}\) is the stress due to the homogeneous strain, \(\sigma _{ij}^{{{{\mathrm{het}}}}}\) is the stress due to the heterogeneous strain, and \(\sigma _{ij}^{{{\mathrm{m}}}}\) is the stress-free strain given by

$$\begin{array}{ll}\sigma _{ij}^{{{{\mathrm{hom}}}}} = R_{in}^tc_{nmkl}^\prime \bar \varepsilon _{kl}^\prime R_{mj},\quad \sigma _{ij}^{{{{\mathrm{het}}}}} = R_{in}^tc_{nmkl}^\prime \eta _{kl}^\prime R_{mj},\\\sigma _{ij}^{{{\mathrm{m}}}} = R_{in}^tc_{nmkl}^\prime e_{kl}^\prime R_{mj}.\end{array}$$
(18)

respectively. In heterogeneous structures, Eq. (17) is solved by using the perturbation method36,37,38. However, the elastic constant tensor varies significantly in nanocrystalline materials, so we directly simulate the dynamics of the displacement. After solving Eqs. (15) and (17), we calculate the total strain ξij using Eq. (11) and thus solve for the total strain in the crystal frame \(\xi _{ij}^\prime\). Substituting \(\xi _{ij}^\prime\) into Eq. (10), we obtain the effective magnetic field due to magnetostriction \(H_i^{{{{\mathrm{ms}}}}}\).

The flow of our simulation model is summarized in Fig. 2. First, we set material parameters, such as the saturation magnetization and the magnetostriction constants, and initialize the configurations of the magnetization and the lattice displacement. After we initialize the model, we calculate the dynamics of the magnetization and lattice displacement. The calculations of the magnetization dynamics are outlined within the blue box. The total strain ξij is calculated by using Eq. (11); then, the effective fields are calculated by Eqs. (2) and (10). The time evolution of the magnetization is calculated based on the LLG Eq. (1). During each time step of the LLG equation, the displacement of the lattice is assumed to be constant. The calculations of the lattice dynamics are outlined within the red box. The strains are calculated by Eqs. (15) and (16), and then the force is calculated by Eq. (18). Finally, the time evolution of the displacement is calculated by Eq. (17). During the lattice dynamics calculation, the magnetization is assumed to be fixed. We save the elapsed time of the magnetization and the lattice dynamics. When the elapsed time exceeds output time interval, the magnetization and the displacement are saved in files. Once the elapsed time exceeds the end time, the simulation is terminated.

Fig. 2: Flow chart of the micromagnetic simulation model.
figure 2

The blue region represents the magnetization dynamics calculations. After starting the simulation, we set the material parameters such as the saturation magnetization and the magnetostriction constants and initialize the configuration of the magnetization and the lattice displacement. After calculating the total strain ξij using Eq. (11), the effective fields are calculated by Eqs. (2) and (10). The time evolution of the magnetization is calculated based on the LLG equation (Eq. (1)). The red region represents the lattice dynamics calculations. The strains are calculated by Eqs. (15) and (16). Then, the force is calculated by Eq. (18). Finally, the time evolution of the lattice is calculated by Eq. (17). During the lattice dynamics calculation, the magnetization is assumed to be fixed. After these procedures, the simulation progresses to the next time step t + Δt. We output the magnetization and displacement configurations at each output time interval, and the simulation stops when the elapsed time exceeds the end time of the simulation.

Full size image

Simulation methods

We simulate the magnetization dynamics by solving the LLG equation using the effective field due to magnetostriction. The two (three)-dimensional model of size 2048 × 2048 × 2 nm3 (256 × 256 × 256 nm3) consists of 32,947 (17,671) nanocrystallites whose mean diameter is 12.7 nm (12.1 nm) to clarify the mean and the local magnetostriction. The nanocrystallite diameter has a Gaussian distribution whose standard deviation is 0.92 nm and 0.64 nm for the two- and three-dimensional models, respectively, and the crystal axes are oriented randomly. The distribution of the nanocrystallites is obtained by a simple molecular dynamics method45. We place the mass points randomly, and repulsive forces act between the mass points. After performing the molecular dynamics simulation under periodic boundary conditions, we obtain the center of each nanocrystal. Finally, we divide the simulation model using Voronoi decomposition.

We choose the α-Fe parameters for the nanocrystallite (except for the anisotropy and the magnetostriction constants) as follows: μ0Ms = 2.15 T, A = 25 pJ/m, \(c_{11}^\prime\) = 2.41 × 1011 N/m2, \(c_{12}^\prime\) = 1.46 × 1011 N/m2, and \(c_{44}^\prime\) = 1.12 × 1011 N/m2. The magnetostriction constants vary in a range from −20 × 10−5 to 20 × 10−5, and we neglect the crystal anisotropy to focus on the effects of the magnetostriction. For the time evolution of ui, we set the mass density ρ and the damping constant \(\zeta\) to 7.87 × 103 kg/m3 and 1.5 × 1012 kg/m3 s, respectively. The damping constant dampens the oscillation of ui to obtain stable strain within the time step of the magnetization dynamics.

The simulation models are discretized into calculation cells of 2.0 × 2.0 × 2.0 nm3 to simulate the magnetization and lattice dynamics using the finite difference method. Because the calculated cell size is smaller than the diameter of the nanocrystal, we can consider the changes in the strain within the nanocrystal. We employ the Runge–Kutta–Fehlberg method (RKF45) with an adjustable time step for the time evolution of the magnetization and the displacement46. We set the time evolution accuracies to 1.0 × 10−4 and 1.0 × 10−8 for the magnetization and lattice dynamics, respectively, and the safety factor is 0.95 in both cases. These accuracies and safety factors enable us to simulate the magnetization dynamics in magnetic materials47,48. The magnetization and lattice dynamics are simulated under periodic boundary conditions. In the case of the two-dimensional simulation model, we adopt the plane strain condition \(\xi _{33} = \xi _{23} = \xi _{31} = 0\). We parallelize the simulations using a message-passing interface49,50,51,52,53, and the magnetization dynamics inside the magnetic materials are reproduced by our simulator.

Results and discussion

Mean strain and distribution of the magnetostriction of uniformly magnetized NSMMs

To clarify the inhomogeneity of the local strain, we calculate the mean and distribution of the magnetostriction for uniformly magnetized NSMMs consisting of α-Fe nanocrystals. The magnetostriction constants of α-Fe are \(\lambda _{100}^\prime\)= 20.7 × 10−6 and \(\lambda _{111}^\prime = - 21.2 \times 10^{ - 6}\). We perform the simulation using two- and three-dimensional NSMMs, which represent the thin film and bulk, respectively. The system is divided into nanocrystals, as shown in Fig. 3a, b, which represent the two- and three-dimensional models. The nanocrystallites are distributed over the NSMM using a simple molecular dynamics method45. The diameter of the nanocrystallite has a Gaussian distribution whose standard deviation and mean are 0.92 nm (0.64 nm) and 12.7 nm (12.1 nm), respectively, in the two (three)-dimensional model. The two- and three-dimensional models contain 32,947 and 17,671 nanocrystals, respectively. When the number of nanocrystallites exceeds 10,000, the mean magnetostriction converges. The directions of the crystal axes of the nanocrystallites are randomly distributed. The NSMM model consisting of nanocrystallites is discretized into cubic cells with a side length of 2 nm, as shown in Fig. 3c. Each cell has its own displacement vector u. The magnetization unit vectors in all unit cells are assumed to be (1, 0, 0). The simulations are performed following the flow chart depicted in Fig. 2.

Fig. 3: Simulation models and magnetostriction distribution in the saturated states.
figure 3

a Two-dimensional simulation model of a uniformly magnetized NSMM. The Cartesian coordinates of the absolute frame are also shown. The size is 2048 × 2048 × 2 nm3, and the magnetization unit vector is m = (1, 0, 0). The simulation model consists of 32,947 nanocrystallites whose mean diameter, D, is 12.7 nm in the x1x2 plane, as shown in the inset. b Three-dimensional simulation model with a size of 256 × 256 × 256 nm.3 The Cartesian coordinates of the absolute frame are also shown. The simulation model contains 17,671 nanocrystallites whose mean diameter, D, is 12.1 nm, as shown in the inset. The magnetization unit vector is uniformly (1,0,0) throughout the model, the same as the two-dimensional model. c Schematic explanation of the displacement vector u. We discretize the simulation model consisting of nanocrystallites into cubic cells with a side length of 2 nm, where each cell has its own displacement vector. After deformation due to magnetostriction, the cell center moves by u. d Color map displaying the distribution of the strain ξ11 in the two-dimensional model. e Three-dimensional plot corresponding to the partial region of (d), which is enclosed by the white dotted square.

Full size image

The calculated values of the mean strain \(\bar \xi _{11}\) for the two- and three-dimensional systems are −9.89 × 10−6 and −8.86 × 10−6, respectively, while the experimental value of the mean strain is −8.9 × 10−6 22. The good agreement of the mean strain between the three-dimensional result and experimental value of the bulk sample shows the validity of the simulation model. Compared with the sum of the stress-free strain, i.e., \(2/5\lambda _{100}^\prime + 3/5\lambda _{111}^\prime\) (−4.4 × 10−6), our results are much closer to the experimental values, which implies the importance of the stress acting on adjacent nanocrystals. Color maps of the strain in the two-dimensional system are shown in Fig. 3d, e. The distribution of the strain is not uniform, and most nanocrystallites have negative strain because the magnetostriction constant \(\lambda _{111}^\prime\) is negative, which results in negative mean strain. As shown in Fig. 3e, which represents the three-dimensional plot of the region marked with the white dotted square in Fig. 3d, the strain varies even inside each nanocrystallite because of the deformation due to the stress from the other nanocrystals. Our simulator can treat the local magnetostriction by calculating the lattice deformation in the NSMM.

Static magnetoelastic properties of the stripe domain

In most applications, NSMMs have a multidomain magnetic structure, where the domains are separated by domain walls54,55,56,57. Because the orientation of the magnetization is nonuniform, it is necessary to consider the effects of the magnetic domain on the magnetostriction. We consider the simplest domain, i.e., stripe domain, in the two-dimensional model of the NSMM. The real NSMMs have defects, impurities, and complex domain structures that affect domain wall motion. The basic effects of the magnetostriction on the domain wall motion are clarified by removing interplays between these factors. The stable stripe domain structure of the NSMM was obtained as follows. The simulation model was first divided into four regions in the x2-direction, and the magnetization direction in these regions was set alternately in the +x1− or −x1-direction. We then calculate the stripe domain structure of a single crystal with the same material parameters as the initial magnetization configuration of the NSMM to prevent the artifactual generation of Bloch points that disturb the domain wall motion. We obtain the stable stripe domain structure after relaxing the NSMM model for 25 ns, as shown in Fig. 4a.

Fig. 4: Static properties of the magnetostriction in NSMM with stripe domain.
figure 4

a Stripe domain structure in NSMM when the mean ξ11 vanishes (\(\lambda _{100}^\prime = 20.0 \times 10^{ - 5}\), \(\lambda _{111}^\prime = - 7.22 \times 10^{ - 5}\)). The colors represent the polar angle of the magnetization vector measured from the x1 axis in the absolute frame. b Color map of the mean strain \(\bar \xi _{11}\) in the \(\lambda _{100}^\prime \lambda _{111}^\prime\) plane. The dotted line of \(\lambda _{111}^\prime = - 0.36\lambda _{100}^\prime\) represents the condition in which \(\bar \xi _{11}\) disappears. c Color map of the mean strain \(\bar \xi _{22}\) in the \(\lambda _{100}^\prime \lambda _{111}^\prime\) plane. The dotted-dashed line of \(\lambda _{111}^\prime = - 0.49\lambda _{100}^\prime\) represents the condition in which the strain \(\bar \xi _{22}\) disappears. d \(\lambda _{100}^\prime\) dependence of \(\bar \xi _{11}\) and \(\bar \xi _{22}\) for \(\lambda _{100}^\prime = - 0.36\lambda _{100}^\prime\) and \(\lambda _{111}^\prime = \lambda _{100}^\prime\). e, f Color maps of ξ11 when \(\lambda _{111}^\prime = \lambda _{100}^\prime\) and \(\lambda _{111}^\prime = - 0.36\lambda _{100}^\prime\), respectively, in the x1x2 plane. In both panels, \(\lambda _{100}^\prime = 20.0 \times 10^{ - 5}\). g Profile of m1 at x1=0 µm as a function of x2. The blue and red colors represent the cases of \(\lambda _{100}^\prime = 20.0 \times 10^{ - 5}\), \(\lambda _{111}^\prime = - 7.22 \times 10^{ - 5}\) and \(\lambda _{100}^\prime = 20.0 \times 10^{ - 5}\), \(\lambda _{111}^\prime = 20.0 \times 10^{ - 5}\), respectively. Semitransparent regions represent the domain wall widths Wm, and dotted lines are tangential lines. h, i Color maps of \(\left\langle {\xi _{11}} \right\rangle _{{{{\mathrm{DW}}}}}\), which is averaged within a diameter of the domain wall width Wm corresponding to (e) and (f), respectively.

Full size image

Figure 4b presents a color map of the mean strain \(\bar \xi _{11}\) in the \(\lambda _{100}^\prime \lambda _{111}^\prime\) plane. \(\bar \xi _{11}\) disappears on the dotted line of \(\lambda _{111}^\prime = - 0.36\lambda _{100}^\prime\) due to the random distribution of crystal axes. Figure 4c presents the same plot for the mean strain \(\bar \xi _{22}\), demonstrating that \(\bar \xi _{22}\) similarly disappears on the dotted-dashed line of \(\lambda _{111}^\prime = - 0.49\lambda _{100}^\prime\) due to the random distribution of crystal axes. As shown in these figures, the magnetostriction constants \(\lambda _{100}^\prime\) and \(\lambda _{111}^\prime\) should have opposite signs to considerably reduce the mean strain. Figure 4d plots the mean strains as a function of \(\lambda _{100}^\prime\). The blue and red circles connected with solid lines represent the strains \(\bar \xi _{11}\) and \(\bar \xi _{22}\), respectively, satisfying \(\lambda _{111}^\prime = - 0.36\lambda _{100}^\prime\); under this condition, \(\bar \xi _{11} = 0\), and \(\bar \xi _{22}\) is negligible. The blue and red triangles connected with dotted lines represent \(\bar \xi _{11}\) and \(\bar \xi _{22}\), respectively, for isotropic magnetostriction, i.e., \(\lambda _{111}^\prime = \lambda _{100}^\prime\). Under this condition, the random distribution of crystal axes does not help suppress \(\bar \xi _{11}\). \(\bar \xi _{22}\) under this isotropic magnetostriction condition is approximately one order of magnitude smaller than \(\bar \xi _{11}\) because the strains in the domain and in the domain wall compete with each other. Moreover, the strain in the domain is tensile, while that in the domain wall is compressive. Figure 4e displays a color map of ξ11 in the x1x2 plane for isotropic magnetostriction, where \(\lambda _{111}^\prime = \lambda _{100}^\prime = 20.0 \times 10^{ - 5}\). When the strain has spatial variation, the effective field of the magnetostriction can disturb the motion of the domain wall. The region plotted in this figure is 0.5 µm × 0.5 µm at the center of the system, where the domain wall appears at ~x2 = 1 µm. Over this entire region, ξ11 is positive, and its mean is ~10 × 10−5. Figure 4f presents the same plot for \(\lambda _{100}^\prime = 20.0 \times 10^{ - 5}\) and \(\lambda _{111}^\prime = - 7.22 \times 10^{ - 5}\), where \(\bar \xi _{11} = 0\). The positive and negative values of ξ11 are scattered at random, reflecting the random distribution of the crystal axes.

According to the exchange softening effect, the effective field is averaged within the exchange length, which is approximately on the scale of the domain wall width. Figure 4g shows a profile of m1 as a function of x2. The semitransparent colors represent domain wall regions whose width Wm is defined using a tangential line58. The estimated Wm is 99.6 nm (76.1 nm) when the mean ξ11 vanishes (the magnetostriction is isotropic). These domain wall widths are comparable to the experiment where magnetic anisotropy exists59. Hence, the exchange averaging effect of the local magnetostriction can be calculated using our simulation model even though the magnetocrystalline anisotropy is neglected. The effects of the magnetostriction on the domain wall motion are averaged within a diameter of the domain wall width Wm. Figure 4h and i shows the strain \(\left\langle {\xi _{11}} \right\rangle _{{{{\mathrm{DW}}}}}\), which is averaged within a diameter of Wm when the magnetostriction is isotropic and the mean magnetostriction vanishes, respectively. While \(\left\langle {\xi _{11}} \right\rangle _{{{{\mathrm{DW}}}}}\) for isotropic magnetostriction is almost invariant (the standard deviation is \(0.232 \times 10^{ - 5}\)), \(\left\langle {\xi _{11}} \right\rangle _{{{{\mathrm{DW}}}}}\) changes on the same scale as the domain wall when the mean magnetostriction vanishes (the standard deviation is \(0.490 \times 10^{ - 5}\)). The spatial variation of ξ11 could act as the pinning center for the domain wall motion and should trigger the loss of the excess energy if the mean strain within the exchange length is nonuniform. The simulation for the static stripe domain clarifies the condition of zero mean magnetostriction and the existence of local magnetostriction.

Effect of magnetostriction on the excess energy loss in the stripe domain

The magnetic cores in the electronic devices are subjected to an oscillating magnetic field that moves the domain walls. We simulate the motion of the domain wall in a two-dimensional NSMM with a stripe domain structure under an oscillating magnetic field. The frequency of the oscillating magnetic field is assumed to be 100 MHz to simulate the domain wall motion within an acceptable computational time, and the amplitude of the field is assumed to be 397.9 A/m, at which the domain wall motion is directly observed. Figure 5a–c shows snapshots of the motion of the stripe domain for \(\lambda _{100}^\prime = 20.0 \times 10^{ - 5}\) and \(\lambda _{111}^\prime = - 7.22 \times 10^{ - 5}\), i.e., \(\bar \xi _{11} = 0\). The width of the stripe domain oscillates with the fluctuation of the domain wall due to local magnetostriction, as explicitly shown in the Supplementary animation. The motion of the domain wall is delayed by magnetostriction with respect to the external magnetic field, as shown in Fig. 5d. The blue curve represents the results with \(\lambda _{100}^\prime = 20.0 \times 10^{ - 5}\) and \(\lambda _{111}^\prime = - 7.22 \times 10^{ - 5}\), where \(\bar \xi _{11} = 0\), whereas the red curve represents the results for isotropic magnetostriction with \(\lambda _{100}^\prime = \lambda _{111}^\prime = 20.0 \times 10^{ - 5}\). The oscillating magnetic field is indicated by the dotted curve. For isotropic magnetostriction, the oscillation of the domain wall exhibits a long delay of ~10% of the oscillation period from the oscillation of the magnetic field. In addition, the oscillation amplitude for the isotropic case \(\lambda _{100}^\prime = \lambda _{111}^\prime = 20.0 \times 10^{ - 5}\) is smaller than that for the zero mean magnetostriction case \(\lambda _{100}^\prime = 20.0 \times 10^{ - 5}\) and \(\lambda _{111}^\prime = - 7.22 \times 10^{ - 5}\), suggesting that the permeability is lower for the former case. In Fig. 5d, the blue curve shows a short delay from the oscillation of the magnetic field. This indicates that a delay between the applied field and polarization occurs even when the mean magnetostriction is zero.

Fig. 5: Dynamic properties of the magnetostriction in NSMM with stripe domain.
figure 5

Snapshots of the motion of the stripe domain in the NSMM at (a) 37.5 ns, (b) 40.0 ns, and (c) 42.5 ns. d Oscillations of the magnetization and the external fields in the x1-direction. The blue curve represents the results with \(\lambda _{100}^\prime = 20.0 \times 10^{ - 5}\) and \(\lambda _{111}^\prime = - 7.22 \times 10^{ - 5}\), where \(\bar \xi _{11} = 0\), while the red curve represents the results for isotropic magnetostriction with \(\lambda _{100}^\prime = \lambda _{111}^\prime = 20.0 \times 10^{ - 5}\). The oscillating magnetic field is indicated by the dotted curve. e \(\bar \xi _{11}\) dependence of the phase differences between the external field and the domain motion. f \(\bar \xi _{11}\) dependence of the oscillation amplitudes of the domain wall motion. g \(\bar \xi _{11}\) dependence of the excess losses. In (e–g), blue and red colors represent the results for zero \(\bar \xi _{11}\) (\(\lambda _{111}^\prime = - 0.36\lambda _{100}^\prime\)) and isotropic magnetostriction (\(\lambda _{111}^\prime = \lambda _{100}^\prime\)), respectively.

Full size image

This delay forms as the result of two phenomena: intrinsic damping, expressed by the Gilbert damping constant, and pinning of the domain wall by local magnetostriction. To separate these two contributions, the intrinsic damping and domain-wall pinning, we plot the \(\lambda _{100}^\prime\) dependence of the phase difference in Fig. 5e, oscillation amplitude in Fig. 5f, and excess loss in Fig. 5g, where the blue and red colors represent the results for zero-mean magnetostriction (\(\lambda _{111}^\prime = - 0.36\lambda _{100}^\prime\)) and isotropic magnetostriction (\(\lambda _{111}^\prime = \lambda _{100}^\prime\)), respectively. As shown by the symmetric curves in Fig. 5e, the phase difference is an even function of \(\lambda _{100}^\prime\) and increases with increasing \(\lambda _{100}^\prime\). The phase difference increases by 500% at \(\lambda _{100}^\prime = 20.0 \times 10^{ - 5}\) when the magnetostriction is isotropic, while the phase difference at \(\lambda _{100}^\prime = 0\) is approximately 0.018×2π, which represents the contribution from intrinsic damping. Additionally, the results for zero-mean magnetostriction are less sensitive to \(\lambda _{100}^\prime\). In contrast, the changes in the oscillation amplitude are smaller than those in the phase difference, as shown in Fig. 5f. Even when the magnetostriction is isotropic, the amplitude is reduced by only 22% at \(\lambda _{100}^\prime = \lambda _{111}^\prime = 20.0 \times 10^{ - 5}\), and the amplitude is almost constant when the mean magnetostriction is zero. Figure 5e, f shows that magnetostriction mainly affects the retardation of the motion of the domain wall.

This retardation and the decrease in the oscillation amplitude are responsible for the excess loss. The excess loss We is obtained by

$$W_e = {\int}_0^T {{{{\boldsymbol{H}}}}_{{{{\mathrm{ext}}}}} \cdot \frac{{d{{{\boldsymbol{B}}}}}}{{dt}}dt,}$$
(19)

where \({{{\boldsymbol{B}}}} = \mu _0\left( {{{{\boldsymbol{H}}}}_{{{{\mathrm{ext}}}}} + {{{\boldsymbol{M}}}}} \right)\) is the magnetic flux density, T is the period of the oscillating magnetic field, and \({{{\boldsymbol{M}}}} = M_s{{{\boldsymbol{m}}}}\). When calculating the excess loss, we adjust the amplitude of the external magnetic field to maintain the same mean oscillation amplitude m1 = 0.17. This is consistent with the common practice in experiments where the core loss is measured at a constant maximum flux density. Furthermore, the excess losses for isotropic magnetostriction and zero-mean magnetostriction are shown in Fig. 5g as functions of \(\lambda _{100}^\prime\). As we discussed in Fig. 5g, the excess loss shows a tendency to increase with an increase in the absolute value of \(\lambda _{100}^\prime\). This increase in the excess loss qualitatively explains the correlation between the anomalous loss factor and the magnetostriction constants in the experiments31,60. This dynamic loss induced by magnetostriction can be explained by the following mechanism. Namely, when the domain wall displacement occurs, the strain within the domain adjacent to the forefront of the wall will be distorted. The cost of the magnetic energy for this process becomes high when the magnetostriction is large. At the same time, the distorted strain will be relaxed at the tale of the moving wall with dissipation of the elastic energy. Hence, the magnetostriction generates excess loss for NSMMs with isotropic magnetostriction in which the distribution of the magnetostriction is flat, as shown in Fig. 4h. Even when the mean magnetostriction is zero, the excess loss increases due to local magnetostriction, as shown in Fig. 4i. The motion of the domain wall is disturbed by local magnetostriction, and the retardation causes excess loss. The excess loss at \(\lambda _{100}^\prime = 0\) is ~51.5 J/m3, which represents the contribution from intrinsic damping (Gilbert damping).

Conclusion

We formulated a micromagnetic simulation model of NSMMs including the effects of magnetostriction and simulated the strain and motion of the stripe domain to clarify the effects of magnetostriction on excess loss. We simulated both the magnetization dynamics and the lattice dynamics to calculate the excess loss in the NSMM. The rotation of the crystal frame is necessary to calculate the effective field due to magnetostriction because of the random distribution of the crystal axes within the NSMM. The magnetization dynamics are simulated by solving the LLG equation using the effective field due to the magnetostriction obtained by simulating the lattice deformation. The two- and three-dimensional simulation models of the NSMM both contain over 10,000 nanocrystallites whose mean diameters are ~12 nm.

The mean magnetostriction of the uniformly magnetized NSMM (−8.86 × 10−6) is in good agreement with the experimental value (−8.9 × 10−6), which confirms the validity of the simulation model. The distribution of the strain is nonuniform, and the strain varies even inside a single nanocrystallite due to the stress from adjacent nanocrystals. The simulation clarifies the condition under which the mean magnetostriction of the stripe domain in the two-dimensional model is eliminated; the mean magnetostriction disappears when \(\lambda _{111}^\prime = - 0.36\lambda _{100}^\prime\). Moreover, even when the mean magnetostriction disappears, the positive and negative strain values are scattered in the NSMM. The identification of the condition to remove the mean magnetostriction enables the improvement of the performance of the NSMM.

Mean and local magnetostriction both delay the motion of the stripe domain under an oscillating external magnetic field. Magnetostriction mainly affects the retardation rather than the amplitude of the motion of the domain wall. This retardation and the decreased amplitude are responsible for the excess loss even when no eddy current is present. Furthermore, even when the mean magnetostriction disappears, the excess loss increases as the magnetostriction constant increases due to local magnetostriction. The simulation results show that considering the mean and local magnetostriction is essential for minimizing the excess loss in the NSMM. The simulation models enhance the performance of NSMMs by reducing the effects of magnetostriction.

The micromagnetic simulation clarifies the effects of magnetostriction on the excess loss of NSMMs. However, the present simulation model does not include the residual amorphous phases whose volume fraction and chemical composition may be controlled through the crystallization kinetics of amorphous precursors61,62,63. In future work, we will perform micromagnetic simulations to reveal the effects of the residual amorphous phase on the energy loss of NSMMs.