Magnetic-field driven domain wall evolution in rhombohedral magnetostrictive single crystals: a phase-field simulation

INTRODUCTION

Giant magnetostrictive materials, i.e., quasi-binary (Tb-Dy)Fe₂ alloys, have attracted considerable attention due to their potential applications as energy conversion components in magneto-mechanical transducers and actuators^[1-6]. During the past decades, efforts have been made to develop an anisotropy compensation system that minimizes magnetic anisotropy in order to take advantage of the giant magnetostriction at low fields^[7-10]. Newnham^[11] firstly credits the spin reorientation boundary of (Tb-Dy)Fe₂ as a magnetic equivalent to the morphotropic phase boundary (MPB) of a quasi-binary ferroelectric solid solution, i.e., the ferromagnetic MPB, which has been confirmed by Yang et al.^[12] and Bergstrom et al.^[13]. While the ferromagnetic MPB has demonstrated utility as a design approach for achieving highly sensitive magnetoelastic responses, the predominant focus in current research has been on materials exhibiting a rhombohedral structure, exemplified by Tb_0.3Dy_0.7Fe₂ (Terfenol-D)^[13,14]. This is because of the highly anisotropic magnetostriction λ₁₁₁ in the pre-transitional rhombohedral side of ferromagnetic MPB, in which giant magnetostriction could be induced by movement of non-180° domain walls or rotation of magnetic moments^[1,14-16].

As is well known, the magnetization of the rhombohedral (Tb-Dy)Fe₂ domains is expected to distribute equally along one of the eight <111> easy axes, as shown in Figure 1. The magnetostriction “jump” effect is primarily attributed to the motions of 109° and 71° domain walls under applied magnetic fields^[17-19], although the underlying internal mechanism remains unclear. Several phenomenological approximations, such as domain wall motion and magnetization rotation models, have been proposed to elucidate the magneto-mechanical behaviors of Terfenol-D^[8,20-23]. However, the prevailing mechanisms largely overlook the presence of internal magnetic (stress) fields and assume that each domain evolves independently, which means no magnetic and elastic long-range interactions among domains are considered. Furthermore, the inherent stress resulting from elastic incompatibility in domain walls is of considerable significance, particularly in the case of giant magnetostrictive materials with large intrinsic strain. Since Khachaturyan’s microelasticity theory^[24] was introduced into the phase-field method, it has emerged as a robust technique for simulating the mesoscale microstructural evolution of ferroelectric systems^[25-28]. Zhang et al. have proposed a phase-field model that integrates the Khachaturyan microelasticity theory with micromagnetic simulation, enabling the prediction of domain structure stability and temporal evolution^[29]. Subsequently, comparable micromagnetic microelastic phase-field models have been developed to explore domain evolution and magnetoelastic response in proximity to the ferromagnetic MPB^{[15,16,30,31]}.

Figure 1. Schematic orientations for eight <111> easy axes in the cubic (Tb-Dy)Fe₂ unit cell.

In this work, the micromagnetic microelastic phase-field method was employed to study the formation and evolution of domains in (Tb-Dy)Fe₂ single crystals within the rhombohedral phase region of the ferromagnetic MPB. Particular emphasis was placed on the two distinct categories of domain walls, namely 71° and 109° domain walls. In the case of the 71° domain wall, the application of an external magnetic field below the coercive field leads to the observation of a domain wall broadening effect, while exceeding the coercive field results in domain switching. Additionally, the applied external magnetic field in the [100] direction leads to 109° domain walls sweeping, which results in heterogeneous magnetization switching.

MATERIALS AND METHODS

The phase-field model employs the local magnetization vector M (m₁, m₂, m₃) as the primary order parameter, with its spatial distribution representing the magnetic domain microstructure. Therefore, the domain structure can be derived from the time evolution of the local magnetization configuration as governed by the Landau-Lifshitz-Gilbert (LLG) equation, denoted as

where M_s represents the saturation magnetization, and α and γ₀ denote the damping constant and gyromagnetic ratio, respectively.

The effective magnetic field, denoted as H_eff, is calculated by H_eff = -(μ₀M_s)^-1(δE_tot/δm) with the vacuum permeability μ₀ and the total free energy E_tot. The total free energy, expressed as^[1]

is composed of the magnetocrystalline anisotropy (E_ani), exchange (E_exc), magnetostatic (E_ms), elastic (E_el), and external energy (E_ext).

The magnetocrystalline anisotropy energy can be defined as^[20,32]

where K₁ and K₂ represent the magnetocrystalline anisotropy coefficients, and V indicates the total volume of the system.

The exchange energy is determined solely by the spatial variation of the magnetization direction and can be written as^[29]

where spatial differentiation is denoted by a comma and A is the exchange stiffness constant.

The magnetostatic energy of a system can be denoted as^[29,33]

where H_d represents the stray field resulting from the long-range interaction between magnetic moments within the system. In the context of a periodic boundary condition, the simulation system is conceptualized as a repetitive building block within 3-D space, and the stray field can be expressed as

where N_D is the demagnetizing factor dependent solely on the sample’s shape. Additionally, $$\overline{\mathbf{M}}$$ signifies the average magnetization of the simulation system, and ϕ represents the magnetic scalar potential solved utilizing the Fourier spectral method under periodic boundary conditions.

The elastic energy resulting from local deformation can be given as^[20,24,34]

where c_ijkl represents the elastic stiffness tensor, e_ij denotes the elastic strain, ε_ij signifies the total strain, and $$\varepsilon_{i j}^{0}$$ indicates the stress-free strain. In cubic magnetostrictive materials, the stress-free strain refers to the spontaneous lattice deformation associated with local magnetization and is denoted by

Where λ₁₀₀ and λ₁₁₁ are the magnetostrictive constants. Khachaturyan’s elastic theory^[24] posits that the total strain ε_ij can be expressed as the sum of homogeneous strain $$\overline{\varepsilon}$$ and heterogeneous strain $$\varepsilon_{i j}^{het}$$, i.e., $$\varepsilon_{i j}=\bar{\varepsilon}+\varepsilon_{i j}^{h e t}$$. The heterogeneous strain adheres to an integral relationship, denoted as $$\iiint_{V} \varepsilon_{i j}^{h e t} \mathrm{~d} V=0$$. The equilibrium heterogeneous strain satisfies the mechanical equilibrium condition outlined by the Euler equation for elastic displacement, given as

where the stress component σ_ij is calculated using

Given the assumption of elastic equilibrium at each evolutionary step, the strain and stress values are determined by solving the mechanical equilibrium equation $$\nabla \cdot \sigma=0$$ using a Fourier spectral method derived from Khachaturyan’s elasticity theory^[24]. When a system is exposed to a homogeneous applied stress $$\sigma_{i j}^{a}$$, the total potential energy is denoted as $$E_{p}=E_{e l}-\int_{V} \sigma_{i j}^{a} \bar{\varepsilon}_{i j} \mathrm{~d} V$$. Furthermore, the competition between stress-induced anisotropy and magnetic crystal anisotropy has the potential to influence the magnetization behavior of the material^[35].

The external energy resulting from the influence of an externally applied magnetic field, denoted as H_ex, can be established as^[36]

The LLG equation is employed to elucidate the progression of domain microstructure, and its numerical solution is achieved through the Gauss-Seidel projection method^[37]. The simulation focuses on the representative rhombohedral ferromagnetic single system, utilizing material parameters derived from a combination of experimental data and theoretical computations conducted previously^[1,38-41], as outlined comprehensively in Table 1. Figure 1 provides a schematic illustration of the unit cell of (Tb-Dy)Fe₂, which is cubic. The process of domain formation was simulated using the phase-field method, with a dimension of 512Δx × 512Δx × 1Δx. (The 3D phase-field simulation of the domain formation process can be seen in Supplementary Figure 1 in Supplementary Materials) The simulation grid Δx is 2 nm, smaller than the exchange length $$l_{e x}=\sqrt{\frac{2 \mathrm{~A}}{\mu_{0} M_{s}{ }^{2}}} \approx 4.7 \mathrm{~nm}$$^[42]. A periodic boundary condition that accommodates non-zero strain is imposed along the three coordinate axes.

RESULTS AND DISCUSSION

Figure 2 illustrates the evolution of domain formation within a representative volume element of the rhombohedral ferromagnetic single crystal, guided by energy minimization towards equilibrium. The initial state is characterized by a random distribution of magnetization, devoid of any predetermined assumptions. The local energy minimum, in conjunction with the presence of inhomogeneous internal stress, serves as the driving force for the initiation nucleation and subsequent growth of various domains. (The distribution of stress during the domain formation process can be seen in Supplementary Figure 2 in Supplementary Materials). At the early stage, metastable rhombohedral and tetragonal variants function as intermediates, providing low-energy kinetic pathways to final engineered domain configurations. Similar bridging domain mechanism was also identified in previous works^[16,31]. Subsequently, the tetragonal variants gradually disappear as the studied component is in the rhombohedral phase side. In the rhombohedral variants intermediate stages, as shown in Figure 2D and E, the domains of the rhombohedral phase form twins of either {100} or {110} twin planes, where the twin boundaries are 109° and 71° ferroelectric domain walls, respectively. As the evolution of rhombohedral domains progresses, there is a gradual disappearance of 71° domain walls, ultimately resulting in the survival of two 109° domain walls in the final equilibrium state. It is worth noting that Figure 2H just depicts one of the ideal states of global energy minimization. For the (Tb-Dy)Fe₂ single crystals, 71° and 109° domain walls usually coexist to form complex multi-domain patterns, which collectively affect the magnetoelastic response. To gain deeper insights into the features of both domain walls and their implications for the properties of rhombohedral ferromagnetic single crystals, the detailed evolution behavior of these domain walls under magnetic fields oriented along the $$[\overline{1}\overline{1}0]$$ and [100] directions have been investigated.

Figure 2. Domain formation process of the rhombohedral (Tb-Dy)Fe₂ single crystal. The colors correspond to the intensity of magnetization along the z direction, while the black arrows indicate the orientations of magnetization.

We first study the evolution of the 71° domain wall below the coercive field. As shown in Figure 3, the initial configuration of the simulation is a domain structure containing four 71° domain walls, which consist of [111] and $$[11\overline{1}]$$ rhombohedral twins. The periodic triangular wave magnetic field with a peak value of 30 kA/m is applied, paralleling the $$[\overline{1}\overline{1}0]$$ direction. In order to illustrate the detailed evolution of the domain wall, we magnify any of the 71° domain walls in the display. The average magnetization along [110] and [001] directions, denoted as M_[110] and M_[001], respectively, are also plotted in Figure 3. When the applied magnetic field increases along $$[\overline{1}\overline{1}0]$$ direction, which is opposite to the magnetization, the 71° domain wall becomes narrower, as shown from t = 0 to t = T/4. When the applied magnetic field increases along [110] direction, which is in the same direction of the magnetization, the 71° domain wall becomes wider, as shown from t = T/2 to t = 3T/4. The domain walls reach their maximum thickness (t = T/4) and minimum thickness (t = 3/4 T) as the field approaches the positive and negative peaks, respectively. The magnetic field dependence of M_[110] also demonstrates the narrowing and broadening of the 71° domain walls. Specifically, there is no domain switching throughout the whole process. The value of M_[001] stays constant, indicating that the center of 71° domain walls remains motionless.

Figure 3. The broadening effect of the 71° domain walls occurs when the applied magnetic field is below the coercive field. The domain wall is magnified to illustrate the fluctuations in its thickness.

Magnetization switching is observed upon surpassing the coercive field, as depicted in Figure 4. When the applied magnetic field in the $$[\overline{1}\overline{1}0]$$ direction increases and surpasses the coercive field, the M_[110] value transitions to a negative state, accompanied by the [111] and $$[11\overline{1}]$$ domains switching to the $$[\overline{1}\overline{1}1]$$ and $$[\overline{1}\overline{1}\overline{1}]$$ domains. At t = T/2, as the applied field diminishes to zero, the $$[\overline{1}\overline{1}1]$$ and $$[\overline{1}\overline{1}\overline{1}]$$ domains are stable. When the applied magnetic field in the reverse direction ([110] direction) surpasses the coercive field, M_[110] undergoes a positive jump, accompanied by the $$[\overline{1}\overline{1}1]$$ and $$[\overline{1}\overline{1}\overline{1}]$$ domains switching to the [111] and $$[11\overline{1}]$$ domains. As the reverse applied field decreases to zero at t = T, the [111] and $$[11\overline{1}]$$ domains are stable. In particular, as shown in the dashed line boxes in Figure 4, the significant fluctuation observed in the M_[001] curve indicates the complex domain switching (the domain switching can be seen in Supplementary Figures 3 and 4 in Supplementary Materials), which is different from the homogeneous polarization switching observed in ferroelectric domains^[43]. The magnetic domain switching near t = T/4 and t = 3T/4 additionally induces a sudden alteration in the magnetostriction of the single crystal (the magnetostriction can be seen in Supplementary Figure 5 in Supplementary Materials). Throughout the aforementioned domain switching, the position of the 71° domain walls remains constant.

Figure 4. Magnetization switching occurs upon surpassing the coercive field. The domains before and after switching are depicted on the left.

Figure 5 shows magnetic field dependence of the 109° domain wall, which consists of $$[11\overline{1}]$$ and $$[\overline{1}11]$$ rhombohedral twins. At step = 0, a constant magnetic field of 20 kA/m is promptly applied along the [100] direction. Under the applied magnetic field, the favored $$[11\overline{1}]$$ domain expands, displacing the 109° domain walls towards the disfavored $$[\overline{1}11]$$ domain. The M_[100] gradually increases from the beginning to step = 7,000, corresponding to the linear movement of 109° domain walls. After that, the curves of M_[100] and M_[010] vibrate violently and remain at constant values. This phenomenon arises due to the convergence and subsequent merger of the two 109° domain walls, leading to their eventual disappearance. The $$[\overline{1}11]$$ domain undergoes a complete transition to $$[11\overline{1}]$$, facilitated by the movement of the 109° domain walls, resulting in the formation of a single magnetic domain. Heterogeneous strain occurs in the domain switching process to accommodate the local deformations. If zero-strain periodic boundary conditions are implemented, the mechanical clamping effect prevents the merging of the two 109° domain walls, thereby inhibiting the formation of a single magnetic domain. It is important to acknowledge that the slower movement of the 109° domain wall can also be induced by applying fields (2 kA/m for example) significantly lower than those examined in this study. (The 109° magnetic domain walls move fast under the magnetic field exceeding 20 kA/m, as shown in Supplementary Figure 6). This implies the absence of a critical magnetic field necessary to induce the displacement of 109° domain walls. This occurs because the $$[11\overline{1}]$$ domain is favored by the [100] applied field, whereas the $$[\overline{1}11]$$ domain is in the disfavored region with higher energy. The two domain regions compete to drive the 109° domain wall in between from the favored domain region to the disfavored one. Generally, this local heterogeneous domain switching can be used in domain engineering^[16] to tailor the magnetic domain morphology with enhanced magnetostrictive properties.

Figure 5. The moving of 109° domain walls and the occurrence of heterogeneous magnetization switching in response to a constant magnetic field of 20 kA/m. The merging and disappearance of two domain walls resulted in the formation of a single domain.

CONCLUSIONS

This study utilized micromagnetic microelastic modeling through the phase-field method to investigate the process of domain formation and the evolution of domain walls in (Tb-Dy)Fe₂ single crystals situated in the vicinity of the rhombohedral region of the ferromagnetic MPB. Particular emphasis was placed on the 71° and 109° domain walls, as their alterations in domain structure under applied magnetic field are crucial to the magnetelastic response of giant magnetostrictive materials. In the case of low applied magnetic field, a phenomenon of domain wall broadening was noted on the 71° domain walls, whereas in the case of high applied magnetic field, homogeneous magnetization switching took place without any observable movement of domain walls. The magnetization switching also helps understand the magnetostriction "jump" effect of the rhombohedral single crystal. The act of sweeping the 109° domain walls resulted in the occurrence of heterogeneous magnetization switching through the movement of domain walls. The detailed analysis of the two domain evolution mechanisms provided insightful understanding into the engineered domain structures of rhombohedral (Tb-Dy)Fe₂ single crystals. Furthermore, these findings offer valuable guidance for the future design of magnetostrictive materials through domain engineering.