## Abstract

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*(*m*_{1},*m*_{2},*m*_{3}), 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.

## 1. Introduction

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 [7], in the 1960s Brown laid the foundation of the micromagnetics theory [8]. The micromagnetic model uses the celebrated Landau–Lifshitz–Gilbert (LLG) equation [9] 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
**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 [26], renormalization [14,22], special test function [24] and Lagrange multiplier [23]. 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 *A*_{s}(∥**M**∥−*M*)^{2}, in order to handle the constraint in an energetic formulation [18]. 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 *A*_{s}. Wang & Zhang [15] 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**(*m*_{1},*m*_{2},*m*_{3}), 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. *ϑ*_{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 *σ*_{ij} is the Cauchy stress and *f*_{i} 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}/∂*x*_{j} in which *x*_{j} is the *j*th Cartesian coordinate direction. The Einstein summation convention is used for the repeated indices. Two types of boundary conditions can be introduced
*n*_{j} is the outward surface unit vector and *ε*_{ij} can be expressed as the symmetrical gradient of the displacement field *u*_{i}
*B*_{i} 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
*ϕ* is the scalar magnetic potential and *H*_{i} is defined by the negative gradient of *ϕ*

### (b) Balance law for magnetization

The magnetic-mechanical couple in ferromagnetic materials results from the existence and rearrangement of ferromagnetic domains [4]. 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 [8], 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 *m*_{i}. The configurational force system for the configurational quantity *m*_{i} includes the configurational stress tensor *Σ*_{ij} whose power density expended on the surface is *g*_{i} and *γ*_{0} is the gyromagnetic ratio with a positive constant value of 1.76×10^{11}/(*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

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

### (c) Thermodynamics

As pointed out in the Introduction, although the use of **m**(*m*_{1},*m*_{2},*m*_{3}) 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
*μ* 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 *g*_{i} 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 *Π*_{μ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
*ζ*_{μ}=*A*_{μi}*g*_{i} is the internal configurational force with respect to *ϑ*_{μ}, and

Application of (2.17) and the form of

The residual dissipation inequality (2.21) is satisfied by setting
*β*_{μγ} is the components of a positive semi-definite matrix which can be derived through a non-concave dissipation potential *ϑ*_{μ} as

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 [9]
*α*=*γ*_{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:
**H**^{eff} 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

### (d) Magnetic enthalpy

In the micromagnetic model, the free energy determined by the magnetization configuration contains the elastic contribution *C*_{ijkl} is the component of the material elastic tensor, *A*_{e} is the exchange constant, and *K*_{1} and *K*_{2} are the magnetocrystalline anisotropy constants. Spontaneous magnetization-induced strain _{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
** ϑ** can be set as

**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 10^{17}. 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 10^{11}, 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

The magnetic enthalpy, the field equations, the constitute relations and the evolution equation take the corresponding dimensionless form as follows:
*ϑ*_{μ,j*}, as it is shown in (2.37).

### Remark

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 [31]. The displacement ** u**(

*u*

_{1},

*u*

_{2},

*u*

_{3}), 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

*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 *f*_{i} in (3.3) and the external configurational force *ϕ*. With these simplicities, corresponding to the strong forms in (3.3), (3.4) and (3.7), the following three weak forms could be formulated
*η*^{ϕ} and *u*_{i}, *ϕ*, 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

The superscript *I* denotes the node number.

Inserting the (3.10)–(3.13) into the weak forms (3.9) and taking the integration over the volume

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. *m*_{i,j}*n*_{j}=*A*_{μi}(*ϑ*_{μ,j}*n*_{j})=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 *t*_{n} are denoted by *t*_{n+1}

From the residuals, the stiffness matrix *ϕ*^{J},

## 4. Simulation results

As a typical magnetostrictive material, Fe_{81.3}Ga_{18.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 §4*a*, 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 Fe_{81.3}Ga_{18.7} has a similar magnitude as Fe (^{3}. 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 §4*b*, the choice of the sample size 30×60×3 nm^{3} and the application of our continuum phase field model are legitimate. For all the examples, the boundary condition *ϑ*_{μ,j}*n*_{j}=0 is used. This leads to the boundary condition *m*_{i,j}*n*_{j}=*A*_{μi}(*ϑ*_{μ,j}*n*_{j})=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 [34]. 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 (**H**^{eff}) is unable to rotate **m** towards its direction, and **m** rotates itself around **H**^{eff} 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 **H**^{eff} 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 2*a*, the initial magnetization is in the *x*_{1}-direction, while the external magnetic field **H**^{ex*} is applied in the *x*_{3}-direction. A field of *H*^{ex*}=1.0 is exerted by applying a magnetic potential difference to the two surfaces perpendicular to the *x*_{3}-direction. The simulated rotation path of **m** is presented in figure 2*b*, for different damping coefficient *α*. In the case of *α*=0, namely no damping, **m** rotates around **H**^{ex*}, and the rotation trajectory is an in-plane circle perpendicular to the *x*_{3}-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 **H**^{ex*}. 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 **H**^{ex*}.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***⋅

**=0 and the damping coefficient**

*n**α*=1.0. We start with an initial random distribution shown in figure 3

*a*and then perform finite-element simulation till the equilibrium state is obtained. As it is shown in figure 3

*b*, 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 [20]. 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

*m*

_{2}-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 nm

^{3}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 [35]. Figure 3*c* shows the comparison of the out-of-plane component *m*_{2} along the line CC^{′} shown in figure 3*b*. A good agreement can be seen between the measured and the calculated results. Furthermore, the in-plane component *m*_{1} along the line DD^{′} is also considered. As demonstrated in figure 3*d*, 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 [35].

### (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 3*a* is taken as the initial configuration. The external magnetic field *H*^{ex*}=1.0 is applied antiparallel to the *x*_{3}-direction. The damping coefficient is set as 1.0. As it is shown in figure 4*a*, 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 *m*_{3}, as it is shown in figure 4*b*. 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 [36].

### (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 §4*a*, 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 [37]. 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 *m*_{2}. 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 5*a*,*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 5*c*,*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 [25]. 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 [15].

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 6*b*. 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 7*a*. The mechanical load is applied by assigning displacement. The surface at *b* 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 *x*_{3}-direction at the expense of domains with the magnetization along the *x*_{1}-direction, as shown in figure 7*e*. On the contrary, under the compressive loading, elliptical vortices with major axis perpendicular to the compressive direction appear (figure 7*c*). The fraction of domains with magnetization along the *x*_{1}-direction remarkably increases (figure 7*f*). These results are expected. A tensile stress increases the regions of domains with magnetization parallel to the *x*_{3}-direction, whereas a compressive stress decreases them. By examining the out-of-plane magnetization component *m*_{2} along the line through the vortex, we can quantitatively obtain the changes in vortex, as shown in figure 7*g*. With respect to the *x*_{1}-direction, the vortex becomes narrower under tensile loading while wider under compressive loading.

## 5. Conclusion

In conclusion, we have proposed a continuum constraint-free phase field model for ferromagnetic materials. Unlike conventional phase field models which take **m**(*m*_{1},*m*_{2},*m*_{3}) 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.

## Funding statements

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.

## Appendix A

With

- Received July 4, 2014.
- Accepted August 8, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.