Calibration of material parameters based on 180∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^\circ $$\end{document} and 90∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^\circ $$\end{document} ferroelectric domain wall properties in Ginzburg–Landau–Devonshire phase field models

Ferroelectric phase field models based on the Ginzburg–Landau–Devonshire theory are characterized by a large number of material parameters with problematic physical interpretation. In this study, we systematically address the relationship between these parameters and the main properties of ferroelectric domain walls. A variational approach is used to derive closed form solutions for the polarization fields at the phase transition regions as well as for the propagation velocities of the domain walls. Introducing a modified set of material parameters, which appropriately scales different contributions to the free energy, we are able to accurately calibrate these parameters based on domain wall thickness and energy of both 180∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^\circ $$\end{document} and 90∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^\circ $$\end{document} domain walls. Moreover, the mobility parameter appearing in the Ginzburg–Landau evolution equation can be accurately calibrated based on the propagation velocity of the domain walls.

carry bound charges resulting from jumps in the polarization component normal to the domain wall, see Ref. [4].
The phenomenon that certain directions and magnitudes of polarization are more favorable than others can be modeled by introducing a polarization-dependent potential energy, the so-called Landau energy. Such type of energy, originally applied to the theory of superconductivity, has been used by Devonshire [5][6][7] to model the behavior of ferroelectric phase transitions phenomenologically. Extending the Landau energy by a term that depends on the gradient of the polarization yields the Ginzburg-Landau-Devonshire theory, which allows the modeling of stable domain walls.
In order to obtain more than one energetically stable state of polarization, the polynomial expansion of the Landau energy needs to be of at least fourth order. For such an energy expansion, Ishibashi [8,9] derives analytical solutions for the polarization fields at 90 • domain walls based on a variational approach. However, these solutions are not expressed in a closed form dependent on the material parameters. Since fourth-order polynomial expansions of the Landau energy fail to model the first-order transitions between the ferroelectric and paraelectric phase, the Landau energy is generally described by a polynomial of sixth order, see, for example, Refs. [10][11][12]. Analytical solutions for the polarization fields at 180 • and 90 • domain walls for a sixth-order energy expansion are given by Cao and Cross [13] and Refs. [14][15][16]. More recently, some authors suggested to use eighth-order polynomial expansions to capture more complex material behavior [17][18][19][20], which leads to increasing difficulties in finding analytical solutions in a closed form.
In the multidimensional setup, the derivation of analytical solutions becomes a cumbersome or even impossible task, leading to an increasing interest in the numerical solution of ferroelectric phase field models. These models couple the Ginzburg-Landau-Devonshire theory with the theory of linear piezoelectricity. The resulting set of field equations can then be solved on a computational domain by choosing suitable numerical methods, such as the finite difference method [21,22], the finite element method [23,24] or the Fourier spectral method [25][26][27][28]. Although the phase field model can be formulated in three dimensions, see, for example, Ref. [25], a two-dimensional model is used in this article for the sake of comprehensibility. All presented derivations could be, however, applied to three dimensions in an analogous manner.
Coming from the phenomenological Ginzburg-Landau-Devonshire theory, the large number of material parameters in ferroelectric phase field models and their influence on the properties of domain walls cannot be interpreted easily. This is the reason why these parameters are calibrated, for example, through atomic level simulations [29] or uncertainty analysis [30]. Schrade et al. [31][32][33][34] propose a calibration of some model parameters based on 180 • domain wall properties but do not address 90 • domain walls in detail.
In this paper, we aim at addressing the relation between the material parameters in ferroelectric phase field models and the properties of both 180 • and 90 • domain walls. In contrast to previous studies, we adopt a variational approach, which has the advantage of being extendable to higher-order polynomial expansions of the Landau energy, and introduce a modified set of material parameters, which appropriately scales different components of the free energy. Two variational formulations of the phase field model are presented in Sect. 2, namely a "static" and a "dynamic" variational formulation. Through the former, analytical solutions are derived as closed form expressions for the polarization fields at 180 • and 90 • domain walls in Sect. 3 and 4 , respectively. In Sect. 5, based on the dynamic variational formulation and in the presence of an applied electric field, we aim at improving the analytical solutions for the velocity of the domain walls available in the literature [31]. The new analytical solutions are finally used for the calibration of the mobility parameter appearing in the phase field model.
In the following, in contrast to scalars, vectors and second-order tensors are denoted by bold symbols, see, e.g., the second-order identity tensor I, with components I i j = δ i j , where δ i j is the Kronecker delta. The tensors of order higher than two are denoted by blackboard bold symbols, e.g., the fourth-order elastic stiffness tensor and the third-order piezoelectric tensor are denoted by c and e, respectively. Taking symmetry considerations into account, c and e can be alternatively written in Voigt notation as matrices c and e. Two consecutive tensors denote a tensor product with summation over all indices of the right tensor. Scalar products are denoted by a dot and imply the summation over all indices of the tensors. In some formulas, the comma convention is used to abbreviate partial derivatives Variations are denoted by δ, e.g., δu is the variation of the field u. Partial variations are labeled by a corresponding subscript, e.g., the partial variation of the functional F with respect to the field u is denoted by δ u F.

Spontaneous polarization and strain
The idea of ferroelectric phase field models is to describe the varying material properties of the domain structure by introducing a continuous field, the so-called phase field. Ferroelectric domains are then modeled by regions with a constant phase field, and domain walls are represented by continuous changes in the phase field. Following Schrade et al [31], the phase field is chosen to be the spontaneous polarization P caused by the non-centrosymmetric ionic structure of ferroelectric materials. For the sake of simplicity, P is shortly denoted as polarization in the following. At a fixed temperature, the magnitude of polarization in a ferroelectric domain free from external influences is denoted by P 0 . In some of the following formulas, the polarization is scaled or normalized such that Besides the polarization, the non-centrosymmetric atomic structure results in non-cubic unit cells in the ferroelectric crystals. This effect is described by introducing the so-called spontaneous strain. In an unloaded ferroelectric domain, the magnitude of the normal spontaneous strain in the direction of polarization is denoted by ε 0 . The normal spontaneous strain perpendicular to the direction of polarization is assumed to be −νε 0 , where ν is a material parameter. For a material polarized in x 2 -direction, this results in the following spontaneous strain in the two-dimensional case According to Ref. [1], the spontaneous strain can be assumed to be volume preserving, resulting in ν = 0.5 in the three-dimensional case. Assuming volume preserving spontaneous strain in the two-dimensional plane strain case, ν = 1 is chosen in all simulations. Based on the experience of the authors, however, the influence of the choice of ν on the calibration process and the resulting calibration errors is marginal. For other states of polarization, the spontaneous strain is calculated by The term written in brackets in the equation above describes a tensor rotation with respect to the direction of polarization of the spontaneous strain state given in Eq. (3) with ν = 1 and the factor in front of the brackets scales the spontaneous strain such that a larger magnitude of polarization leads to larger strains and vice versa. More information regarding the spontaneous strain can be found in Sect. B of the Supplemental Material and Ref. [1].

Energy contributions
Ferroelectric phase field models couple the Ginzburg-Landau-Devonshire theory with the theory of linear piezoelectricity, resulting in two primary variables in addition to the phase field, which are the displacement u and the electric potential φ. Their gradients define the infinitesimal strain tensor ε and the electric field E Dependent on the primary variables and their gradients, the free energy density f free is formulated. It can be written as the sum of five contributions which are the mechanical, piezoelectric, electrical, Landau and gradient energy density, respectively. The mechanical, piezoelectric and electrical energy density describe the energy stored in the mechanical and electrical fields and are given by Schrade et al. [31] as where the spontaneous strain is incorporated as eigenstrain in the mechanical energy density. The tensors c, e and are the fourth-order elastic stiffness tensor, the third-order piezoelectric tensor and the secondorder absolute permittivity tensor, respectively. In general, the material behavior described by these tensors is anisotropic with respect to the direction of polarization, suggesting a dependence of the tensors on the polarization vector p as, for example, successfully implemented in Ref. [32]. However, in order to reduce the complexity of the numerical simulations and motivated by the circumstance that the effect of anisotropy is large on the piezoelectric tensor but comparatively small on the elastic stiffness tensor and the absolute permittivity tensor, c and are here assumed to be isotropic, i.e., independent on the polarization or the crystal orientation. Due to its symmetry properties, the fourth-order elastic stiffness tensor can be written as a matrix using Voigt notation The absolute permittivity tensor is given by In contrast to the elastic stiffness tensor and the absolute permittivity tensor, the dependence of the piezoelectric tensor on both the magnitude and direction of polarization is not negligible. If the magnitude of polarization equals P 0 and the direction of polarization is the x 2 -direction, the piezoelectric tensor is expressed in Voigt notation as the following matrix e( p 1 = 0, p 2 = 1) = 0 0 e 15 e 31 e 33 0 .
Rotation and scaling of the piezoelectric tensor given above yields For more information on the piezoelectric tensor, the reader is referred to Sect. C of the Supplemental Material and Ref. [1].
As the Landau and gradient energy density depends on the crystal orientation, a rotation tensor is defined as which maps the phase field and its gradient from the coordinate system of interest x to the coordinate systemx that is aligned with the crystal lattice, ψ being the angle between the two systems. Applying the coordinate transformation results inx The Landau and gradient energy density is taken from Schrade et al. [31], who introduced two parameters related to the thickness and energy of 180 • domain walls. For the purpose of relating the energy parameters to the properties of both types of domain walls, two additional parameters are introduced in this paper. This results in a total number of four parameters denoted by κ land 180 , κ land 90 , κ The Landau energy density, which models the potential energy density associated with the state of polarization, is described by a sixth-order polynomial dependent on the polarization. Due to symmetry requirements, i.e., f land (p) = f land (−p), only polynomial terms with even exponents are considered. Furthermore, following Schrade et al. [31], the polynomial termsp 4 1p 2 2 andp 2 1p 4 2 are not considered in order to keep the calibration process feasible, resulting in The parameters α 1 , α 11 , α 111 and α 12 are chosen in such a way that the polynomial has four minima corresponding to the four energetically stable directions of polarization in the two-dimensional setting, i.e., Further, choosing yields zero Landau energy density at the energetically stable states, see Ref. [31]. For 180 • Ising-type domain walls, the termp 2 1p 2 2 in the Landau energy vanishes. Therefore, the parameter κ land 90 can be used to calibrate the 90 • domain wall properties without influencing the 180 • domain wall properties.
Finally, the gradient energy density penalizes gradients in the phase field and is taken as dependent on the partial derivatives of the polarization in which the partial derivatives are taken with respect to the coordinate systemx. In contrast to Schrade et al. [31], who assumed κ grad 180 = κ grad 90 , the partial derivativesp 1,2 andp 2,1 are penalized differently thanp 1,1 andp 2,2 . A similar treatment can be found in Ref. [15]. During the calibration process presented in this paper, 180 • domain walls are assumed to be charge neutral, as this is the energetically favorable configuration, see Ref. [4]. Neutral 180 • domain walls are aligned with the crystal lattice, resulting inp 1,1 =p 2,2 = 0. The parameter κ grad 90 can hence be used to calibrate the 90 • domain wall properties without influencing the 180 • domain wall properties.
Integrating the free energy density over the domain of interest yields the free energy The body force density b and free charge density ρ f acting on the computational domain as well as the surface traction t and free surface charge σ f acting on the Neumann boundaries, see Sect. 2.3, are contributing to the external energy Polarization switching in ferroelectric materials is a thermodynamically irreversible process, resulting in dissipation of energy. Dependent on the time derivative of the polarizationṖ, the dissipation rate is given bẏ As the time-dependent evolution of the phase field is assumed to be isotropic, the mobility parameter β is scalar.

Initial and boundary conditions
Initial and boundary conditions are imposed on the primary variables and their gradients. In the initial state, the primary variables are prescribed at every point of the domain. In order to formulate mechanical boundary conditions, the boundary of the computational domain ∂ , with outer normal unit vector n, is split into the Dirichlet boundary ∂ u with prescribed displacement u * and the Neumann boundary ∂ t with prescribed Similarly, ∂ is split into two parts in order to impose the electrical boundary conditions. The electric potential φ * is prescribed on the Dirichlet boundary ∂ φ , and the free surface charge σ * f is prescribed on the Neumann boundary ∂ σ where ∂ φ ∪ ∂ σ = ∂ and ∂ φ ∩ ∂ σ = ∅. For the phase field, homogeneous Neumann boundary conditions are assumed

Variational formulations
Based on the energy contributions presented in Sect. 2.2, two variational formulations are introduced in the following. More details can be found in Miehe et al. [35]. In the first variational formulation, in the following referred to as static variational formulation, the energy of the system is defined as the sum of the free energy and the external energy The objective is to find the fields u, φ and P satisfying the given Dirichlet boundary conditions such that the static energy F stat is stationary, which leads to the necessary conditions However, due to the non-convexity of the Landau energy density, this energy has multiple stationary points and hence the problem is ill-posed. A well-posed variational problem is formulated by defining the dynamic energy rate asḞ The aim of the dynamic variational formulation is then to find the stationary point of the dynamic energy rate with respect to the rates of the fieldsu,φ andṖ satisfying the given initial conditions and Dirichlet boundary conditions. This leads to the necessary conditions

Discretization
In general, the stationary points of the presented energy functionals cannot be found analytically. Therefore, numerical solution strategies are required, one of which is described in this section. The numerical approach requires discretization in time and space. For the discretization in time, the time interval of interest is split into a finite number of equally sized intervals [t n , t n+1 ]. The time step is then defined by t = t n+1 − t n , and the time derivative of the polarization is discretized by using the implicit Euler schemė For the discretization in space, the finite element method is chosen. The computational domain is divided into a finite number of elements and at each element the unknown fields and their variations are approximated by a linear combination of nodal unknowns and ansatz functions. Here, the Bubnov-Galerkin ansatz is chosen, in which the unknown fields are approximated in each element by where n is the number of nodes per element, N i are the ansatz functions andũ i ,φ i andP i are the nodal unknowns. The variations of the unknown fields δu e , δφ e and δ P e are discretized in the same manner. In the simulations presented in this article, quadrilateral elements with bilinear ansatz functions are used.
As the static variational formulation is ill-posed, the finite element method is applied to the well-posed dynamic variational formulation. Discretization of the necessary conditions in Eq. (29) in time and space as shown in Eqs. (30) and (31) leads to a nonlinear algebraic system of equations at every time step, that is solved iteratively by using the Newton-Raphson method. The integrals are computed numerically by applying Gaussian quadrature with two integration points per element in each parametric direction. More information on the numerical solution strategy can be found in Refs. [23,24,[31][32][33][34]].

Material parameters
For the comparison between the analytical solutions derived in this paper and the corresponding numerical simulations, the ferroelectric material lead zirconate titanate PZT-5H is considered exemplarily. All material parameters, except the mobility parameter and the parameters to be calibrated, are taken from Schrade et al. [31], see Table 1.

Numerical Solution
The region of continuous change in polarization between two ferroelectric domains polarized in opposite directions is known as a 180 • domain wall. In this section, the formation of a 180 • domain wall is simulated numerically by applying the finite element method to the dynamic variational formulation. Although it has been shown that 180 • domain walls exhibit mixed Ising-Bloch-Néel-type character, see Ref. [3], the Bloch-and Néel-type characteristics are assumed to be small throughout this work and are not captured by the presented phase field model, i.e., both the numerical treatment presented in this section and the analytical treatment presented in the following section result in pure Ising-type domain walls.
A square computational domain with a side length of 30 nm is discretized by using 301 × 10 equally sized quadrilateral elements. Over the whole domain, the crystal lattice orientation is assumed to be aligned with the coordinate system, i.e., ψ = 0 • . In the initial state, the fields u and φ are set to zero at all nodes. Subsequently, two ferroelectric domains are generated by assuming opposite initial polarization directions in the left and right half of the computational domain. The mechanical boundary conditions are chosen to be statically determinate, i.e., u is set to zero at the node in the lower left corner and u 2 is set to zero at the node in the lower right corner. Homogeneous Dirichlet boundary conditions are chosen for the electrical potential. No body force or free charge density is considered.
The material parameters shown in Table 1 are considered, and the remaining parameters of the Landau and gradient energy density as well as the mobility parameter are chosen as The simulation is executed with t = 2 × 10 −11 s until a steady state is reached. In the converged state, the first component of the polarization P 1 is close to zero and the second component P 2 is smoothly changing between the two ferroelectric domains, see Fig. 1. For more information on u and φ, the reader is referred to Sect. D of the Supplemental Material.

Analytical solutions
Analytical solutions for the polarization field at the 180 • domain wall can be found in a quasi-one-dimensional setup, i.e., the solution field is assumed to be only dependent on x 1 , see Fig. 1. In contrast to the finite element setting, an infinite domain is assumed. No external forces are considered and the influence of the mechanical, piezoelectric and electrical energy density on the solution is neglected. Under these assumptions, the exact solution for the polarization field is given by [13] as (2α 11 + 6α 111 ),α = α 11 + 3α 111 α 11 + 2α 111 .
As an alternative to Cao and Cross' solution, another analytical solution based on a slightly simplified ansatz is derived in the following by making use of the static variational formulation. The motivation is twofold: first, Cao and Cross' solution is restricted to sixth-order polynomial expansions of the Landau energy density, and second, using Cao and Cross' solution for the parameter calibration in the later sections of this paper leads to very complex mathematical expressions. Under the introduced assumptions, the static energy becomes The objective is to find a polarization field such that the static energy is stationary. Instead of searching in the space of all functions, a smaller ansatz space is considered, namely the space of all polarization fields that can be expressed by where λ 180 is a real-valued parameter. It follows that the gradient of the polarization field has only one nonzero component Inserting the ansatz into the static energy and using Eq. (18) yields The integrals in the expression above can be calculated analytically, which yields an energy expression dependent on the parameter λ 180 . The necessary condition for stationarity thus reads ∂F stat ∂λ 180 The resulting algebraic equation is finally solved for the unknown parameter The derived solution is not unique. Depending on the sign of the parameter λ 180 , the polarization of the domain wall changes from a positive to a negative value or vice versa. In Fig. 2, Cao and Cross' and the new analytical solution as well as the numerical solution from the previous section are plotted. It can be seen that both analytical solutions are very similar to the finite element results. Cao and Cross' solution is based on the assumption that the mechanical, piezoelectric and electrical energy densities are zero. These energy densities, however, are considered in the numerical simulation, leading to small discrepancies between Cao and Cross' solution and the numerical results. The new analytical solution deviates slightly from the other solutions because of the chosen simplified ansatz.

Calibration
Running the finite element simulation based on the phase field model yields a 180 • domain wall that is characterized by a certain output domain wall thickness out 180 and output domain wall energy F out 180 . The objective of this section is to calibrate the parameters in the model in such a way that these two output parameters can be directly regulated by two input quantities in 180 and F in 180 . As the first component of the polarization is zero at the 180 • domain wall, the rightmost contributions of the Landau and gradient energy density vanish. This means, κ land 180 and κ grad 180 are the only parameters that influence the domain wall, which can also be seen in the analytical solutions in Eqs. (34) and (40). These parameters are calculated based on the requirement that the output quantities be equal to the corresponding input ones.
The domain wall thickness is typically defined by taking the tangent to the phase field profile at the center of the domain wall and measuring the distance between the intersections of the tangent with the limits −P 0 and P 0 (Fig. 2) in 180 Further, the domain wall energy is defined by the integral of the Landau and gradient energy density over the domain wall whereα I is given in "Appendix A." If, for example, the material in Table 1 The equations above contain information on how the calibration parameters need to be chosen such that the new analytical solution or Cao and Cross' solution exhibits the desired domain wall thickness and energy. Due to the discrepancies between the analytical and numerical solution, however, the output domain wall thickness and energy of the numerical solution might deviate from the corresponding input quantities. The calibration accuracy is estimated by defining a calibration error

Analytical solutions
For the quasi-one-dimensional setting, analytical solutions for the polarization field at the 90 • domain wall are given by Cao and Cross [13] as well as Ishibashi [8]. However, as Cao and Cross' solution is restricted to certain choices of the parameters in the Landau energy density and Ishibashi's solution is not expressed in a closed form, a new analytical solution is derived in the following by using the static variational formulation. The derivation is based on the same assumptions as explained in Sect. 3.2. Throughout the derivation, the first component of the polarization is assumed to be constant over the domain wall. The ansatz space is chosen to be the space of all polarization fields that can be described by The only nonzero component of the gradient of the polarization is thus The crystal lattice orientation is not aligned with the coordinate system. Therefore, a rotation of 45 • needs to be applied to the input arguments of the Landau and gradient energy density as shown in Sect. 2.2. This yields under consideration of the chosen ansatz and Eq. (18) the following static energy  In order to find the stationary points with respect to the parameter λ 90 , the integrals are evaluated analytically and the partial derivative with respect to the unknown parameter is set to zero To give an example, for the material in Table 1 The accuracy of the calibration can once again be estimated based on the definition of the calibration error in Eq. (49). In contrast to the 180 • case, the calibration error of the 90 • domain wall quantities appears to be highly dependent on the choice of the input parameters. Therefore, a parameter study is performed, in which the calibration errors are calculated for multiple combinations of the input quantities in 90 and F in 90 . The resulting errors for the domain wall thickness E( 90 ) are illustrated by the grey bars in Fig. 5.
The errors are caused by the assumptions made during the derivation of the analytical solution. In particular, the mechanical, piezoelectric and electrical energy, which are considered in the numerical simulations, were not considered during the derivation of the analytical solution, causing errors in the calibration process. Furthermore, the assumption that the first component of the polarization is constant does not hold true for any choice of input quantities. Input quantities leading to a first component of the polarization that is at every point larger or at every point smaller than assumed in the analytical solution are indicated in Fig. 5 by blue and red squares, respectively. In this case, the errors appear to be comparatively large. However, if the first component of the polarization is close to the one assumed in the analytical solution, i.e., at some points larger and at some points smaller, the input quantities are indicated by green squares, resulting in comparatively smaller errors between −3% and 3%.
The question arises if there is a way to predict analytically whether the calibration error is large or small for a given combination of input quantities. Theoretically, the calibration error should be small if the input quantities are chosen in such a way that the assumption of constant first component of the polarization holds true. The assumption is expected to be approximately fulfilled if the input quantities are chosen such that, for , the gradient of the Landau energy density vanishes at the center of the domain wall This results in a relation between the input quantities in 90 in which α II is a positive real-valued factor that depends on the parameters in the Landau energy density. For the parameters shown in Eqs. (44) and (58), it is α II = 32 15α 1 + 20α 11 + 23α 111 320α 1 + 480α 11 + 563α 111 .
The resulting relation between the input quantities is illustrated by a dashed line in Fig. 5. Close to the line, the errors are smaller than far away from the line. Using Eqs. (46) and (58) instead leads to the factorα II shown in Eq. (A.2) in "Appendix A." The corresponding relation between input quantities is plotted as a dotted line in Fig. 5, showing small differences to the dashed line.
Performing an analogous analysis for the calibration error of the domain wall energy E(F 90 ) yields similar results. However, most of the errors are about two times smaller than the calibration errors of the domain wall thickness.
The error analysis described above shows that the assumption of the constant first component of the polarization has a large influence on the calibration error. If the parameters are chosen such that the first component of the polarization is approximately constant, the calibration errors are marginal. This observation suggests that neglecting the mechanical, piezoelectric and electrical energy during the derivation of the analytical solution, while considering them in the numerical simulations, has a comparatively low impact on the calibration error.

Domain wall dynamics
If subjected to external stimuli, the domains in the bulk ferroelectric change in size and shape. This means that the domain walls, which separate the ferroelectric domains, have to move. In this section, the time-dependent behavior of a 180 • and a 90 • domain wall is studied both numerically and analytically. By making use of the dynamic variational formulation, the objective is to improve the analytical solution for the domain wall velocity and the corresponding calibration of the mobility parameter presented by Schrade et al. [31], whose derivations are based on the Euler-Lagrange equations of the phase field model.

180 • domain wall
An electric field pointing in x 2 -direction is applied to a 180 • domain wall by changing the electrical boundary conditions of the simulation shown in Sect. 3.1. In particular, the electric potential is chosen such that it differs by φ = − 0.03 V between the top and bottom boundary, resulting in an electric field with E 2 = 10 6 V m . At the left and right boundary, homogeneous Neumann boundary conditions are chosen. In order to reduce influences from the boundary, the width of the computational domain is doubled while keeping its height unchanged. Consequently, 601 × 10 finite elements are used.
During the simulation, the domain wall moves to the right until it vanishes at the boundary, which leads to a homogeneous polarization field in the converged state, see Fig. 6. It is observed that the shape of the phase field profile at the domain wall does not change significantly during the process.
The domain wall velocity v 180 is defined as the change in time of the position of the domain wall, where this position is conventionally chosen as the x 1 -coordinate of the zero of the second component of the polarization at x 2 = 0. Fig. 6 shows the domain wall velocity as a function of the position of the domain wall. The initially assumed polarization profile at the domain wall causes a comparatively low velocity at the beginning of the simulation. However, after taking its energetically stable profile shape, the domain wall moves with an increased velocity, which remains approximately constant for some time, until it reaches the boundary towards the end of the simulation, where a rapid acceleration is observed.
An analytical solution for the velocity of the 180 • domain wall is given by Schrade et al. [31], whose derivation is based on the Euler-Lagrange equations of the phase field model evaluated at the center of the domain wall. However, this solution strategy results in errors that are larger than 35%. Therefore, a new analytical solution for the domain wall velocity is derived in the following by making use of the dynamic variational formulation, which allows a global treatment of the problem.
Since no body or surface forces are applied, the dynamic energy rate reduces to the sum of the dissipation rate and the time derivative of the free energẏ Considering a quasi-one-dimensional setup, the dissipation rate is given bẏ asṖ 1 vanishes. Under the assumption that the domain wall moves with a constant velocity v 180 while not changing its profile shape over time, the time derivative of the polarization can be expressed aṡ The dissipation rate hence becomesḞ Because of the assumption that the shape of the polarization profile at the domain wall does not change over time, all contributions of the free energy, except the coupling term between the electric field and the polarization, do not depend on the position of the domain wall in space. Hence, these contributions do not lead to a change in free energy over time, resulting in the following time derivative of the free energy in the quasi-one-dimensional setup as both P 1 and E 1 vanish. The electric field is assumed to be constant in time and space. Hence, considering Eq. (65) leads to Finally, the dynamic energy rate is written aṡ The objective is to find a domain wall velocity v 180 such that the expression above is stationary, which leads to the necessary condition ∂Ḟ dyn ∂v 180 Solving the resulting equation leads to The integral in the equation above can be calculated, for example, by making use of the new analytical solution given in Eqs. (37) and (40) or by using Cao and Cross' solution shown in Eqs. (33) and (34) ∞ −∞ P 2 2,1 dx 1 = where it can be seen thatα ≥ 1 must hold true. It is observed from Eq. (71) that the domain wall velocity is inversely proportional to the mobility parameter β. This relation can hence be used to calibrate the mobility parameter Referring to the corresponding finite element results in Fig. 6, it is seen that the numerically calculated domain wall velocity far away from the boundary differs by less than 1% from the one chosen during the calibration. Hence, the calibration shows an increased accuracy compared to that given by Schrade et al., which eliminates the need for an additional correction parameter as proposed in [31].

90 • domain wall
As follows, the behavior of a 90 • domain wall under an applied electric field is studied. The boundary conditions of the simulation shown in Sect. 4.1 are changed such that the electric field points in x 2 -direction with E 2 = 10 6 V m . Further, the width of the computational domain is doubled and 301×50 finite elements are used. Similar to the 180 • domain wall, the 90 • domain wall moves to the right until it vanishes, leading to a homogeneously polarized domain in the final state. It is observed that, in contrast to the simulation of the 180 • domain wall, the velocity of the 90 • domain wall slightly varies over the x 2 -direction, leading to a change in shape of the polarization profile at the domain wall in time. For the sake of simplicity, the behavior of the 90 • domain wall is studied at x 2 = 0 only. Fig. 7 shows the corresponding evolution of the 90 • domain wall and its velocity. It is observed that, after the jump caused by the initial condition, the 90 • domain wall slightly decelerates at the beginning of the simulation, which can be explained by the change in profile shape over time. When approaching the boundary, the domain wall accelerates rapidly as already observed for the 180 • case.

Conclusions
Analytical solutions were derived for the phase field profiles at 180 • and 90 • domain walls and for their propagation velocities. These solutions were used in order to calibrate the material parameters occurring in the Landau and gradient energy as well as the mobility parameter. Table 2 gives a comparison between the different calibration strategies for an exemplary set of input quantities. In contrast to the previous literature that focuses on the calibration of the material parameters only based on 180 • domain wall properties, in this work we calibrated the parameters based on the properties of both types of domain walls, in both cases with satisfactory accuracy. Further, the analytical solutions for the domain wall velocities derived in this paper show an increased accuracy compared to the literature, leading to an improvement of the mobility parameter calibration.