A continuum constraint-free phase field model is proposed to simulate the magnetic domain evolution in ferromagnetic materials. The model takes the polar and azimuthal angles (ϑ1,ϑ2), instead of the magnetization unit vector m(m1,m2,m3), as the order parameters. In this way, the constraint on the magnetization magnitude can be exactly satisfied automatically, and no special numerical treatment on the phase field evolution is needed. The phase field model is developed from a thermodynamic framework which involves a configurational force system for ϑ1 and ϑ2. A combination of the configurational force balance and the second law of thermodynamics leads to thermodynamically consistent constitutive relations and a generalized evolution equation for the order parameters (ϑ1,ϑ2). Beneficial from the constraint-free model, the three-dimensional finite-element implementation is straightforward, and the degrees of freedom are reduced by one. The model is shown to be capable of reproducing the damping-dependent switching dynamics, and the formation and evolution of domains and vortices in ferromagnetic materials under the external magnetic or mechanical loading. Particularly, the calculated out-of-plane component of magnetization in a vortex is verified by the corresponding experimental results, as well as the motion of the vortex under a magnetic field.
Owing to their ferromagnetic property and magnetic-mechanical coupling, ferromagnetic materials are widely used in industrial applications, such as magnetic data storage, sensors and actuators, transducers and microelectromechanical systems [1–3]. The viable applications and reasonable design of the devices based on ferromagnetic materials are highly dependent on the fundamental understanding of the microstructures of these materials. At the microscale, ferromagnetic materials are composed of discretized magnetized regions in which the magnetization is uniform. These regions are generally called magnetic domains [4,5]. Accordingly, at the macroscale, the macroscopic magnetic-mechanical properties are determined by these magnetic domain evolutions which can be driven by the external magnetic field and/or mechanical loading [4–6]. Hence, developing a model for ferromagnetic materials that can describe or predict the structure and evolution of the magnetic domains is critically important for understanding, designing and engineering the macroscopic properties of the ferromagnetic devices.
Presently, the micromagnetic model and the phase field method are two widely used approaches for modelling ferromagnetic materials. Based on the well known domain wall calculation by Landau and Lifshitz in 1935 , in the 1960s Brown laid the foundation of the micromagnetics theory . The micromagnetic model uses the celebrated Landau–Lifshitz–Gilbert (LLG) equation  to describe the temporal evolution of the magnetization and of the domain structure. In the last two decades, the micromagnetic model has achieved vital success in modelling ferromagnetic materials [10–13]. However, in this model, it is difficult to consider the mechanical effects, especially the inhomogeneous stress resulting from the elastic incompatibility of the magnetostrictive strain [14,15]. By contrast, the conventional phase field model is based on the continuum thermodynamics and the kinematics of the materials. This type of model usually uses the magnetization as the order parameters [14–19]. It can easily take into account the magnetic-mechanical coupling and predict the detailed domain structure and evolution under external magnetic field or mechanical loading without a priori assumption on domain morphologies.
As in most of the work about micromagnetic models and phase field models for ferromagnetic materials, we consider here merely the evolution at a constant temperature which is far below the Curie temperature and the phase transition temperature. Thereby, the magnetization magnitude remains constant and only its direction changes during the domain evolution. This constraint can be given as 1.1 where M is the magnetization vector introduced as a continuum field variable, M is the saturation magnetization magnitude and m is the magnetization unit vector. M is constant when the temperature dependence is ignored. This is intrinsically different from the ferroelectric phase field models with polarization as the order parameter, in which no constraint on the polarization magnitude is enforced [20,21]. In the ferromagnetic materials, the magnetocrystalline anisotropy energy term is non-convex, whose function is similar to that of the Landau energy polynomial in ferroelectrics [20,21]. But it works only when the magnitude constraint of ∥M∥=M is applied. In other words, only when the constraint ∥m∥=1 is considered, this non-convex term can form a multi-well landscape which characterizes the easy axes of the magnetization. This constraint makes the determination of the phase field evolution path challenging and introduces complications in the phase field modelling and numerical implementation [22–25].
In conventional ferromagnetic phase field models, the LLG equation is usually taken as the evolution equation. There the magnitude constraint of ∥M∥=M does not explicitly appear in the formulation but is enforced in the numerical implementation [14,17,19]. This results in a rather simplified presentation. But the numerical implementation of the LLG equation with the constraint is particularly challenging. In the literature, there exist some approximation techniques, such as projection , renormalization [14,22], special test function  and Lagrange multiplier . There are also exact techniques in use of rotation in the Lie group [25,27]. Recently, Landis introduced a constraint energy term in the form of As(∥M∥−M)2, in order to handle the constraint in an energetic formulation . However, in this way, the constraint is only fulfilled at the stationary solution, and the evolution dynamics is not physically sound. Furthermore, it is not so easy to choose a reasonable value for the constraint energy coefficient As. Wang & Zhang  have used this constraint energy proposed by Landis and the time-dependent Ginzburg–Landau (TDGL) equation which is widely used in ferroelectrics to establish their phase field model. It should be noted that if the LLG equation instead of the TDGL equation is used, the contribution of the constraint energy to the LLG equation vanishes; because the equality m×m=0 always holds in the LLG equation. Therefore, the method of constraint energy is not so effective.
In this work, we develop a novel constraint-free phase field model for mechanically coupled magnetic domain evolution in ferromagnetic materials. The focus here is to construct a model which delivers the correct evolution dynamics and fulfils the constraint automatically. The idea is to take the polar and azimuthal angles (ϑ1,ϑ2), instead of the magnetization unit vector m(m1,m2,m3), as the order parameters. It is apparent that the Cartesian components of m can be expressed by (ϑ1,ϑ2) in a constraint-free formulation, i.e. , and . By using the configurational force theory and the second law of thermodynamics, a set of thermodynamically consistent constitutive relations and a generalized evolution equation for ϑ1 and ϑ2 are derived, as shown in §2. In §3, a three-dimensional nonlinear finite-element formulation of this constraint-free model is presented. Thereby no additional numerical technique is required for the magnitude constraint of ∥m∥=1. Furthermore, the node degrees of freedom are decreased by one. In §4, several numerical examples are presented. Benchmark test on precession and precessional switching is checked in the use of one-element simulation, and it demonstrates that the model reproduces the correct magnetization dynamics. Modelling on domain formation and evolution, such as vortices and 180° domain switching, shows that three-dimensional simulations are required in order to recapitulate the spatial switching path. The shape and motion of the computed vortices agree well with the related experimental results. In addition, the ferroelastic switching is also simulated. In comparison with other phase field simulations, the proposed three-dimensional constraint-free model has the merit that it can readily reveal the damping-dependent and the out-of-plane of magnetization in a thin film.
2. Continuum constraint-free phase field model
(a) Field equations
The magnetoelastic coupling is one of the most important properties for ferromagnetic materials. Thus both the mechanics and the magnetostatics should be considered in the modelling. For a ferromagnetic body with a boundary , the quasi-static mechanical equilibrium is described by 2.1 where σij is the Cauchy stress and fi is the body force. Hereafter, the Latin indices (i,j,k,l,p,q) run over the range of 1–3. The comma in a subscript denotes spatial differentiation, for example, σij,j=∂σij/∂xj in which xj is the jth Cartesian coordinate direction. The Einstein summation convention is used for the repeated indices. Two types of boundary conditions can be introduced 2.2 in which is the displacement prescribed on the boundary part , nj is the outward surface unit vector and is the surface traction on the boundary part . By assuming linear kinematics, the strain εij can be expressed as the symmetrical gradient of the displacement field ui 2.3 The Maxwell equation which governs the magnetic part has the form 2.4 where Bi is the magnetic induction. As the permittivity of the ferromagnetic materials is usually much higher than that of the free space, the stray field effect is ignored here. Thus the magnetic boundary conditions can take the form of 2.5 where is the prescribed value on the boundary part , ϕ is the scalar magnetic potential and is the given potential on the boundary part . The magnetic field Hi is defined by the negative gradient of ϕ 2.6
(b) Balance law for magnetization
The magnetic-mechanical couple in ferromagnetic materials results from the existence and rearrangement of ferromagnetic domains . The free energy of the ferromagnetic materials is also dependent on the magnetization configuration. In the continuum theory, the magnetization configuration can be described by the distribution of the magnetization vector. To obtain the distribution, one needs to first consider the balance law of magnetization. Though the magnetization dynamics has been well established , in this subsection we employ the configurational force theory [28–30] to represent the balance law, in order to support the derivation of our phase field model in the following subsections.
As it has been mentioned, we consider only the isothermal process below the Currie temperature and the phase transition temperature. Thus the magnitude of the magnetization vector remains constant, and the configuration can be depicted by the unit magnetization vector mi. The configurational force system for the configurational quantity mi includes the configurational stress tensor Σij whose power density expended on the surface is , and the internal and external configurational force vector gi and whose power density expended in the volume is and , respectively. Owing to the rotation nature of the magnetization , the angular momentum balance in the continuum aspect can be given as 2.7 where γ0 is the gyromagnetic ratio with a positive constant value of 1.76×1011/(Ts), and ϵijk is the permutation tensor. Converting the surface integration in (2.7) into volume integration and considering its validity in any volume, we can obtain 2.8
The magnetization can only change its direction, so the configurational force, which physically is the driving force on the change of magnetization, should lie perpendicular to the direction of the magnetization. In other words, the configurational force along the magnetization direction must be zero. Thus, in any volume, we have 2.9 From (2.8) and (2.9), one has 2.10 with .
As pointed out in the Introduction, although the use of m(m1,m2,m3) as the order parameters is straightforward, the numerical implementation of the constraint on the magnetization magnitude can be complicated [22–25]. Herein, we use the polar and azimuthal angles (ϑ1,ϑ2) as order parameters. Thus the components of the unit magnetization vector in the corresponding Cartesian coordinates are given as 2.11 and 2.12 In this paper, the Greek indices μ and γ run over the range of 1–2. The constraint ∥m∥=1 is satisfied automatically. The number of order parameters is also reduced from 3 to 2.
For the temperature-independent process, the Helmholtz free energy of the ferromagnetic system with magnetic-mechanical couple can be taken as . According to the second law of thermodynamics, under the isothermal condition, the external power expended on the control volume should not be less than the change rate in the Helmholtz free energy, i.e. 2.13 In this inequality, the internal configurational force gi is omitted as it has no external power. The difference between the right-hand and left-hand terms is the dissipation. After the Legendre transformation, the magnetic enthalpy can be derived as . Thus, the thermodynamic inequality (2.13) can be rewritten as 2.14 or by the Gauss law 2.15 with 2.16 Here, Πμj and ζexμ are the configurational stress tensor and the external configurational force, with respect to the order parameters ϑμ, respectively.
Left multiplication of (2.10) by matrix Aμi leads to 2.17 in which ζμ=Aμigi is the internal configurational force with respect to ϑμ, and is the antisymmetric matrix 2.18
Application of (2.17) and the form of to (2.15) leads to 2.19 Considering the field equations (2.1) and (2.4) and noticing that inequality (2.19) must be hold for all admissible processes, based on the standard arguments of rational thermomechanics, we can obtain these constitutive relations 2.20 Given these relations, the residual dissipation inequality in (2.19) in the local form can be derived as 2.21 where 2.22
The residual dissipation inequality (2.21) is satisfied by setting 2.23 where βμγ is the components of a positive semi-definite matrix which can be derived through a non-concave dissipation potential 2.24 Combining (2.17), (2.22) and (2.23), we can obtain a generalized form of the evolution equation for the order parameters ϑμ as 2.25
It should be noted that when the Rayleigh dissipation functional is adopted, (2.25) can be reduced to the standard LLG equation. In fact, inserting the Rayleigh dissipation potential  2.26 into (2.24), one has 2.27 Insertion of (2.27) into (2.25) leads to 2.28 where 2.29 and α=γ0ν is the damping coefficient which generally appears in the LLG equation. As it is well known, the LLG equation is given in terms of m in the following form: 2.30 where Heff is the effective field obtained by variational derivation of the total magnetic enthalpy with respect to m, μ0 is the vacuum permeability constant and × denotes the cross product of two vectors. By using the relation , it is not difficult to prove that (2.30) is equivalent to (2.28). For more details, one is referred to the appendix A.
(d) Magnetic enthalpy
In the micromagnetic model, the free energy determined by the magnetization configuration contains the elastic contribution , the magnetocrystalline anisotropy contribution , the exchange contribution and the magnetostatic contribution . For a cubic crystal, the following magnetic enthalpy is used 2.31 with 2.32 where Cijkl is the component of the material elastic tensor, Ae is the exchange constant, and K1 and K2 are the magnetocrystalline anisotropy constants. Spontaneous magnetization-induced strain , for cubic crystals, can be described as 2.33 where λ100 and λ111 is the magnetostriction strain along the direction 〈100〉 and 〈111〉 of a single crystal, respectively, when it is magnetized at saturation along this direction.
With the above-specified magnetic enthalpy, the constitutive relations in (2.20) can be given as 2.34 The evolution equation (2.28) can be derived as 2.35 where 2.36 The boundary condition and initial condition for ϑ can be set as 2.37 The combination of (2.31)–(2.36) constitutes the constraint-free formulation of the continuum phase field model for ferromagnetic materials. No additional constraint on the magnitude of M is required, because ∥M∥=M has been intrinsically implemented in the model.
3. Finite-element implementation
(a) Dimensionless form
A direct implementation of the constrain-free model in the physical dimension leads to equation systems with large condition numbers in the order of 1017. This can be due to the varying orders of the parameters in the magnetic and the mechanical problem. For example, the elastic constant is in the order of 1011, whereas the vacuum permeability is in the order of 10−6. Hence, a dimensionless form of the constrain-free model is implemented here. The model is normalized via the following ansatz 3.1
The magnetic enthalpy, the field equations, the constitute relations and the evolution equation take the corresponding dimensionless form as follows: 3.2 3.3 3.4 3.5 3.6 3.7 The dimensionless boundary and initial conditions are 3.8 Here, is the boundary value defined for ϑμ,j*, as it is shown in (2.37).
In the next section for the finite element formulation, only the dimensionless form of the phase field model will be employed. The superscript stars are hence omitted for notational simplicity.
(b) Finite-element formulation
The presented phase field model in the constraint-free form is implemented in the three-dimensional finite-element method. The discretization is achieved by eight-node linear elements. For more details, readers are referred to standard textbooks on finite-element methods . The displacement u(u1,u2,u3), the scalar magnetic potential ϕ, and the polar and azimuthal angles ϑ(ϑ1,ϑ2) are taken as independent variables. Thus, each node has six degrees of freedom , where the superscript I indicates the element node and the underbar denotes a matrix. Note that if m is taken as the order parameter, it requires not only additional constrain of ∥m∥=1, but also one more degree of freedom. In this sense, the constraint-free formulation favours an efficient numerical implementation.
The body force fi in (3.3) and the external configurational force in (3.7) are neglected. Actually, can be given by setting the boundary condition of ϕ. With these simplicities, corresponding to the strong forms in (3.3), (3.4) and (3.7), the following three weak forms could be formulated 3.9 where , ηϕ and are the test functions for ui, ϕ, and ϑμ, respectively.
By introducing the shape functions for the independent variables, the rate of the polar and azimuthal angles, and the test functions, the discretized equations are obtained as 3.10
The superscript I denotes the node number. , and are the shape functions for the displacements, the scalar magnetic potential and the polar and azimuthal angles, respectively. In this paper, the underbar denotes Voigt notation of the corresponding quantities. Accordingly, the strain in (2.3), the magnetic field in (2.6) and can be given as 3.11 and the constitutive relations in (3.6) as 3.12 In the three-dimensional case, the matrices in (3.11) can be given as 3.13
Note that the surface terms in (3.14) can be integrated by applying the boundary conditions in the finite-element software. Thereby, they are here ignored. The determination of the boundary conditions on the order parameter is not a trivial task. In the following, we will tacitly assume homogeneous Neumann type boundary conditions, i.e. . Note that this Neumann boundary automatically leads to the boundary mi,jnj=Aμi(ϑμ,jnj)=0 which is widely used in the conventional phase field model.
With respect to the time dependence of the residual, we use the implicit backward Euler time integration . The quantities at the previous time-step tn are denoted by and . The equation for the current time-step tn+1 3.15 should be solved to obtain . For solving these nonlinear equations, the Newton iteration scheme is performed at each time step. The corresponding iteration matrix for one element is 3.16
From the residuals, the stiffness matrix and the damping matrix can be calculated by derivation with respect to , ϕJ, and , , , respectively. Specifically, the non-zero stiffness matrix can be derived as 3.17 3.18 3.19 and the damping matrix as 3.20 All other components of the damping matrix are zero. Note that the stiffness matrix components and the damping matrix components are unsymmetric. This unsymmetric characteristic is intrinsically attributed to the precession nature of magnetization dynamics. The elements are all integrated using a standard eight point Gauss integration scheme. It should be noted that this procedure is in the element level. For the whole simulated object, the above nonlinear equations and iteration matrix should be assembled into the global ones. The model is implemented as a user element in the software FEAP .
4. Simulation results
As a typical magnetostrictive material, Fe81.3Ga18.7 is simulated by using the constraint-free phase field model and its finite-element implementation. The material parameters are taken from the literature [14,33] and are listed in table 1. For the simulation of precession in §4a, a single element test was used. For the study of domains, a thin film nanostructure is simulated, in order to obtain the characteristic features and simultaneously save computation time. The sample mesh is 10×20×1. Given that the exchange coefficient of Fe81.3Ga18.7 has a similar magnitude as Fe (), the original size of the simulated sample is around 30×60×3 nm3. The finite-element mesh size is chosen in a manner that a smooth variation of the magnetization over the whole sample should be resolved by the discretization. This means that the obtained domain wall should span over at least two or more elements. Note that continuum models have been often used to study ferromagnetic domains in nanostructures [10,13,15,17,23]. As it will be shown by the comparison with experimental observations in §4b, the choice of the sample size 30×60×3 nm3 and the application of our continuum phase field model are legitimate. For all the examples, the boundary condition ϑμ,jnj=0 is used. This leads to the boundary condition mi,jnj=Aμi(ϑμ,jnj)=0, which is commonly assumed in the literature (e.g. [8,13,15,23]).
(a) Precession and precessional switching
Precession and precessional switching are basic physical phenomena in magnetization dynamics . If there is no damping (i.e. no dissipation), the instantaneous change in m should be perpendicular to both the external field and the direction of m. In other words, the effective magnetic field (Heff) is unable to rotate m towards its direction, and m rotates itself around Heff with a constant angle, as it is shown by the thin line in figure 1. This phenomenon is the so-called precession. If there is damping, m can be gradually rotated towards the direction of Heff and forms a spiral-type path, as it is illustrated by the thick green line in figure 1. This is called as precessional switching. The spiral-type path can vary with the damping coefficient.
As a benchmark example, the constraint-free phase field model is used to reproduce precession and precessional switching. To simulate this behaviour of a monodomain, a single finite element should be sufficient. As shown in figure 2a, the initial magnetization is in the x1-direction, while the external magnetic field Hex* is applied in the x3-direction. A field of Hex*=1.0 is exerted by applying a magnetic potential difference to the two surfaces perpendicular to the x3-direction. The simulated rotation path of m is presented in figure 2b, for different damping coefficient α. In the case of α=0, namely no damping, m rotates around Hex*, and the rotation trajectory is an in-plane circle perpendicular to the x3-direction. It indicates that the precession behaviour is reproduced. In the cases of α>0, damping dependent precessional switching is observed. When α is very small, for example α=0.01, a large number of precessional loops occur till m is oriented along the direction of Hex*. When α is moderately small, for example α=0.1, much less number of loops is required. When α is relatively large, for example α=1.0, 10 and 100, less than one precessional loop is needed, and m directly rotates towards Hex*.The rotation trajectory lies nearly in one plane. These simulation results of precession and precessional switching indicate that magnetization switching dynamics has been soundly considered in the constraint-free phase field model.
(b) Domain and magnetic vortex
The proposed model is further applied to simulate the domain structure in a thin nanostructure which is free from mechanical and magnetic load, as shown in figure 3. The mechanical boundary condition for all the sample surfaces is traction-free, i.e. σ*⋅n=0. The magnetic boundary is set as B*⋅n=0 and the damping coefficient α=1.0. We start with an initial random distribution shown in figure 3a and then perform finite-element simulation till the equilibrium state is obtained. As it is shown in figure 3b, two vortices are formed in the equilibrium configuration. The vortex, a curling configuration with a core, consists of four domains separated by 90° domain walls. Most interestingly, in the vortex core, the magnetization rotates gradually out of the plane. This phenomenon makes our simulation results essentially different from the vortices given in the literature [23,25]. It should be noted that this simulated vortex with curling and out-of-plane structure is attributed to the constraint on magnetization magnitude. In the vortex core, due to the strictly satisfied constraint ∥m∥=1, the magnetization cannot vanish (must be unit) and turns perpendicular to the surface, thus avoiding the singularity. This is different in the case of ferroelectrics, where the polarization has no constraint in its magnitude and it can remain in-plane by reducing its magnitude to very small values . In fact, with ∥m∥=1, magnetization in the vortex remaining in plane would significantly increase both the exchange and the magnetocrystalline anisotropy energy. As the out-of-plane m2-direction is also the easy axis, the curling and twisting-out structure largely decreases the exchange energy and the magnetocrystalline anisotropy energy and thus is energetically favourable. Note that we have also considered a larger sample with size of 30×60×9 nm3 and found that the domain configuration has no significant difference.
This simulated vortex phenomenon is consistent with the experimental observations in nanoscale ferromagnetic thin films . Figure 3c shows the comparison of the out-of-plane component m2 along the line CC′ shown in figure 3b. A good agreement can be seen between the measured and the calculated results. Furthermore, the in-plane component m1 along the line DD′ is also considered. As demonstrated in figure 3d, the experimental and the simulated results agree well with each other. In particular, the measured and the calculated profiles show a vortex width of around 9 nm, also in accordance with the experimental measurements .
(c) Evolution of the vortex configuration under magnetic loading
Simulation is performed to show domain poling under a magnetic field. As an example, the equilibrium state in figure 3a is taken as the initial configuration. The external magnetic field Hex*=1.0 is applied antiparallel to the x3-direction. The damping coefficient is set as 1.0. As it is shown in figure 4a, the two vortices move oppositely towards the sample corners. The motion trajectory of the vortices is approximately along the in-plane diagonal of the sample. After these vortices vanish at the corners, a uniform configuration is formed gradually. The whole evolution procedure can also be demonstrated through the contour plot of the component m3, as it is shown in figure 4b. During this process, the regions with magnetization parallel to the external magnetic field expand gradually, whereas the regions with magnetization antiparallel and perpendicular to the external magnetic field shrink. It should be noted that the moving direction of the vortex has the component perpendicular to the external magnetic field, which is consistent with the experimental observation .
(d) 180° domain wall evolution under magnetic loading
As another example, the evolution of a 180° domain wall under an external magnetic field is investigated. As reported in §4a, the magnetization switching is dependent on the damping coefficient α. It has been even evidenced that the damping-induced dissipation is the driving force for magnetic domain wall motion . If there is no damping, the motion of the domain wall will not occur. In this subsection, we show that the magnetization switching of the 180° domain configuration can also be damping dependent. Figure 5 shows the distribution of out-of-plane magnetization m2. When the damping is small, e.g. α=0.01, 0.1, the magnetization near the 180° wall rotates out of the plane, as it is shown in figure 5a,b. Whereas, it can be seen that when the damping is large, e.g. α=1.0, 10.0, the rotation of magnetization is almost constrained in the plane, as shown in figure 5c,d. Intrinsically, this behaviour is attributed to the magnetization dynamics. By a LLG-type evolution equation which has only a damping term and no precession term, Miehe et al. showed in-plane rotation of magnetization . Wang et al. used the TDGL equation (not the LLG equation) that intrinsically cannot involve the damping effect and found that most magnetization vectors rotate in plane .
For a comparison, the evolution of the 180° wall in the cases of α=0.1 and 10.0 is specifically given in figure 6. In the case of α=0.1, the magnetization near the wall first rotates out of plane, destroying the sharp 180° wall configuration. Then the magnetization away from the wall also rotates out of plane, and two vortices initiate at the sample boundary. Subsequently, the vortex in the region whose magnetization is along the external magnetic field rapidly disappears. The other vortex forms and moves from one edge towards another edge under the drive of the external magnetic field. Finally, this vortex moves to the up-left corner and disappears there. The magnetization in the whole sample is then aligned along the the external magnetic field. Owing to the out-of-plane rotation of the magnetization and the formation of vortices, this evolution process is relatively complicated. During this process, the 180° wall was totally destroyed, and new magnetization configurations were constructed. By contrast, the evolution process in the case of large damping α=10.0 is much simpler, as one can see in figure 6b. The magnetization merely rotates in the plane and the 180° wall becomes diffusive. The magnetization in the left region, initially antiparallel to the external magnetic field, always rotates counterclockwise to be aligned along the external magnetic field. While the magnetization in the right region, initially parallel to the external magnetic field, rotates firstly clockwise and then counterclockwise. In summary, the poling process of a 180° domain configuration is strongly dependent on the damping coefficient α. It should be noted that a thicker sample with a mesh of 10×20×3 is also calculated and it is found that the evolution process is almost the same.
(e) Domain evolution under mechanical loading
Besides the ferromagnetic switching, the constraint-free phase field model can also simulate the ferroelastic switching. It essentially makes this phase field model different from the micromagnetic models. In the simulation, the equilibrium vortex structure is taken as the initial configuration, as shown in figure 7a. The mechanical load is applied by assigning displacement. The surface at is fixed, and the displacement at the opposite surface is set as for the case of tension and for the case of compression. The damping coefficient is 1.0. It is shown in figure 7b that the tensile loading causes the circular vortices to be elliptical, with the major axis along the tensile direction. As for the domains, the tensile loading increases the fraction of domains with the magnetization along the x3-direction at the expense of domains with the magnetization along the x1-direction, as shown in figure 7e. On the contrary, under the compressive loading, elliptical vortices with major axis perpendicular to the compressive direction appear (figure 7c). The fraction of domains with magnetization along the x1-direction remarkably increases (figure 7f). These results are expected. A tensile stress increases the regions of domains with magnetization parallel to the x3-direction, whereas a compressive stress decreases them. By examining the out-of-plane magnetization component m2 along the line through the vortex, we can quantitatively obtain the changes in vortex, as shown in figure 7g. With respect to the x1-direction, the vortex becomes narrower under tensile loading while wider under compressive loading.
In conclusion, we have proposed a continuum constraint-free phase field model for ferromagnetic materials. Unlike conventional phase field models which take m(m1,m2,m3) as the order parameters, the constraint-free phase field model uses the polar and azimuthal angles (ϑ1,ϑ2). As a result, the constraint on magnetization is satisfied automatically within the model itself, and no additional numerical strategy for the phase field evolution is needed. The model was developed in a thermodynamic framework. Based on a configurational force system for (ϑ1,ϑ2) and the second law of thermodynamics, a set of thermodynamically consistent constitutive relations and a generalized evolution equation for the order parameters (ϑ1,ϑ2) have been obtained. As it has been shown in §3 that the model leads to a straightforward three-dimensional finite-element implementation, which requires no numerical treatment on the constraint and has one less degree of freedom.
The presented numerical simulations evidence that the model can readily predict the fundamental phenomenon in ferromagnetic domain switching. Reproduction of precession and damping-dependent precessional switching behaviour shows that the model gives a physically sound switching dynamics. Results on domain structure in a free-standing thin film demonstrates that the magnetization rotates around the vortex core and gradually twists out of the plane. The magnetization distribution around the vortex and the vortex width agree well with experimental observations. Under an external magnetic field, the vortex motion has a component perpendicular to the magnetic field, which is also consistent with experimental observations. Ferromagnetic switching of 180° domain configuration is damping-dependent. Magnetization rotates out of plane and a vortex tends to form in the case of small damping, whereas magnetization is constrained into in-plane rotation in the case of large damping. Examples have been given for ferroelastic switching, i.e. domain evolution under mechanical loading. It is found that a tensile loading increases the fraction of domains along the tensile direction, whereas a compressive loading increases the fraction of domains perpendicular to the compressive direction. Accordingly, the vortex will change from a circular shape to an elliptical shape.
It can be concluded from the simulations that a three-dimensional model and implementation is indispensable, and the switching mechanism of magnetization is spatial. Furthermore, the ferromagnetic domain evolution is strongly damping dependent. These two features have not been well considered in the conventional phase field models in the literature. In our future work, the mesoscopic response of ferromagnetic materials will be investigated by using this model.
The support from the LOEWE research cluster RESPONSE (Hessen, Germany), the China Scholarship Council and the Innovation Foundation of BUAA for PhD Graduates (YWF-14-YJSY-052) is acknowledged.
With and , the LLG equation (2.30) can be written as A 1 Left multiplying (A.1) with the matrix A 2 can give A 3 Incorporating the external magnetic field , the effective field can be derived as A 4 Thus, A 5 Combining (A.3) and (A.5) yields A 6 which is the same as the evolution equation in (2.28).
- Received July 4, 2014.
- Accepted August 8, 2014.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.