A phase-field model for ferroelectrics with general kinetics. Part I: Model formulation

When subjected to electro-mechanical loading, ferroelectrics see their polarization evolve through the nucleation and evolution of domains. Existing mesoscale phase-field models for ferroelectrics are typically based on a gradient-descent law for the evolution of the order parameter. While this implicitly assumes that domain walls evolve with linear kinetics, experiments instead indicate that domain wall kinetics is nonlinear. This, in turn, is an important feature for the modeling of rate-dependent effects in polarization switching. We propose a new multiple-phase-field model for ferroelectrics, which permits domain wall motion with nonlinear kinetics, with applications in other solid-solid phase transformation problems. By means of analytical traveling wave solutions, we characterize the interfacial properties (energy and width) and the interface kinetics of straight domain walls, as furnished by the general kinetics model, and compare them to those of the classical Allen--Cahn model. We show that the proposed model propagates domain walls with arbitrarily chosen nonlinear kinetic relations, which can be tuned to differ for the different types of domain walls in accordance with experimental evidence.


Introduction
Ferroelectrics are a prime example of multifunctional materials with electro-mechanical coupling. They are largely used as sensors and actuators for their piezoelectric properties (Uchino, 2018) as well as in ferroelectric random-access memory (Ishiwara et al., 2004), which uses the reversibility of spontaneous polarization states. Below a critical temperature, referred to as the Curie temperature, ferroelectric ceramics possess a polar structure with a spontaneous electric polarization, which can be reversed under the application of electro-mechanical loading. Polarization switching is accommodated by transformations between lattice variants of different orientations-e.g., the six orientations of the tetragonal variants for barium titanate (BaTiO 3 ) at room temperature-and between phases of distinct crystal symmetry types-e.g., the tetragonal and rhombohedral phases in some compositions of lead zirconate titanate (PZT).
At the mesoscale (tens of nanometers to hundreds of microns), polarization reversal occurs through the emergence and evolution of a microstructure formed by the different variants (and phases, where applicable) referred to as ferroelectric domains separated by domain walls. In defect-free ferroelectric single-crystals, the microstructure that develops during switching vanishes at the end of that process, which leads to a single polarization domain in the sample (Little, 1955). By contrast, in polycrystalline ferroelectric ceramics, grains with different orientations exert mechanical constraints, so that microstructural arrangements with multiple ferroelectric domains remain after polarization switching (Arlt and Sasko, 1980). At the macroscale, a few representative curves are commonly used to characterize polarization switching. First, bipolar measurements are obtained by cycling the applied electric field between ±e max . These yield the TDGL equations-which were re-derived from variational principles (Zhang and Bhattacharya, 2005a) or microforce balance (Su and Landis, 2007)-with electrostatic Gauss' law and the mechanical balance equations to simulate the evolution of ferroelectric domain patterns in response to an applied electric field, i.e., polarization switching (Zhang and Bhattacharya, 2005a,b;Su and Landis, 2007;Vidyasagar et al., 2017;Indergand et al., 2020).
While these later works successfully accounted for features of the quasistatic p-e-hysteresis and strain butterfly curves, the TDGL model reaches its limitations when it comes to modeling rate effects as well as the material response to pulse switching experiments. For instance, TDGL simulations of polarization switching in response to applied electric field pulses of magnitude e sw yield a switching time τ s that scales as τ s ∝ 1/e sw (Indergand, 2019) in contradiction to the experimental scaling τ s ∝ e 1/esw (Merz, 1956;Schultheiß et al., 2019), known as Merz law.
The above scaling that results from the TDGL model is expected, since the TDGL formulation postulates a linear relation between the time evolution of the order parameter (polarization p) and the (negative of the) variational derivative of the free-energy density Ψ, i.e., with the inverse mobility coefficient µ > 0. As will be discussed in Section 3, when the TDGL model is viewed as the regularization of a sharp-interface model, the evolution law (1.1), also known as the Allen-Cahn equation (Allen and Cahn, 1979), amounts to assuming a linear kinetic relation. As reviewed in Appendix A, this point is precisely in contrast with experimental evidence, which shows that domain wall motion is governed by nonlinear kinetic relations. Hence, modeling time-dependent switching with a diffuse-interface approach requires a new phase-field formulation, one that permits nonlinear domain wall kinetics. This is the objective of the present work, which focuses on the theoretical foundation. Part II of this study will present numerical applications. Domain evolution in ferroelectrics is a specific case of the general class of problems of structural phase transformations. These include, e.g., stress-induced phase transformations (such as the martensite-austenite transformation) and deformation twinning. In this more general setting, the question of how to formulate a regularized model that accounts for the nonlinear kinetics of interfaces in solids was considered before. Hou et al. (1999) used a level-set method based on the Hamilton-Jacobi equation (Osher and Sethian, 1988) to model interface propagation with nonlinear and anisotropic kinetics. Alber and Zhu (2013) proposed a phase-field model for two-phase solid-solid transformations-termed the hybrid model because it shares properties of a Hamilton-Jacobi and a parabolic equation-which permits nonlinear kinetics (see also Zhu, 2005, 2007). The model proposed in Section 4 of this article may be interpreted as a reformulation of that of Alber and Zhu (2013), which is specialized for ferroelectrics and extended to transformations between multiple phases (as opposed to only two in the works of Alber and Zhu (2005). On a different front, Agrawal and Dayal (2015a,b) recently proposed a phase-field model that admits nonlinear kinetics. Their evolution equation closely resembles the one of Hou et al. (1999), reformulated as a phasefield model instead of a level-set method. However, the model of Agrawal and Dayal (2015a) does not lend itself easily to extensions to multiple phases (see Agrawal (2016)). Finally, Tůma et al. (2018) included a mixed-type dissipation potential-combining viscous and rate-independent contributions-in the variational formulation of a phase-field model for displacive phase transformations such as martensitic transformations. This furnishes a nonlinear kinetic law in the corresponding sharp-interface model, one that prescribes a threshold on the driving force for interface motion. This model was later extended by including this feature in a multi-phase-field model for solid-solid transformations, which rests upon a micromorphic formulation (Rezaee-Hajidehi and Stupkiewicz, 2021). While these two formulations introduce in a neat fashion a specific nonlinearity in the kinetic relation (specifically, a threshold on the driving force), our goal here is to include a kinetic relation with general nonlinearity so as to accurately account for experimentally measured kinetics. Within this setting of diffuse-interface models with nonlinear kinetics for solid-solid transformations, the one we introduce in Section 4 is the first to encompass multiple-phase transformations and embed arbitrary nonlinearity. It finds a particularly relevant application in ferroelectric domain evolution, since experimental evidence confirms that the kinetic relation for domain walls is nonlinear. Further, the proposed model allows to select different kinetics for the different types of domain walls (e.g., 90 • and 180 • domain walls in tetragonal ferroelectrics), as it is indicated by experimental results.
The remainder of this article is organized as follows. In Section 2 we introduce a sharp-interface model for ferroelectrics, which constitutes the basis for discussing the characteristics of the diffuse-interface models of subsequent sections. Section 3 reviews the characteristics of classical models based on the Allen-Cahn equation with regards to the properties (interfacial energy and thickness) and kinetics of 90 • and 180 • domain walls. In Section 4, we introduce a new phase-field model for ferroelectrics with general kinetics, and we discuss its properties with focus on the kinetics of the two types of domain walls. The main features of the new phase-field model are discussed in the conclusion in Section 5. Empirical data on the nonlinear kinetics of domain walls in barium titanate coming from the physics literature are briefly reviewed in Appendix A. Numerical applications will be reported in Part II of this study.

The sharp-interface model for ferroelectrics
We begin by formulating a mesoscale continuum model for ferroelectrics. Domain walls are represented by sharp interfaces, across which some or all of the electro-mechanical fields exhibit discontinuities. The motion of a domain wall is governed by a function, termed the kinetic relation, which relates its normal velocity to the thermodynamic driving traction exerted on the interface. The sharp-interface model introduced in this section serves as a reference to discuss the behavior of the diffuse-interface models of Sections 3 and 4.
In Section 2.1, we present the mechanical balance laws, Maxwell's equations and evolution equation for a singular surface, formulated for a general continuum subject to electro-mechanical coupling. In Section 2.2, we introduce an expression for the electric enthalpy of ferroelectrics, from which the constitutive relations derive, and we discuss the functional form of the kinetic relation appropriate for domain wall motion in ferroelectrics. The resulting model is specialized to rigid ferroelectrics in Section 2.3.
Note that we use the following notation of tensor calculus. For any vector fields v(x, t) and w(x, t), second-order tensor fields R(x, t) and S(x, t), and third-and fourth-order tensors T and U, respectively, we have in a Cartesian basis where indices following a comma indicate spatial derivatives, and we use Einstein's summation convention.

General principles
We consider a deformable body B occupying the region Ω ⊂ R 3 in the presence of electrostatic fields filling R 3 . The deformation of B is described by the displacement field u(x, t) : Ω × R → R 3 . Given the small strains encountered in ferroelectric ceramics, we assume linearized kinematics and introduce the infinitesimal strain tensor ε(x, t) = (∇u) sym = 1 2 ∇u + (∇u) T . In addition, we restrict ourselves to quasistatics both in the mechanical sense (i.e., inertia is neglected) and in the sense of electro-quasistatics (i.e., magnetic induction is neglected) and assume that no free charges exist within Ω. The former is justified by the slow rates of interest, while the latter is a common assumption for undoped ferroelectrics. Let S(t) be a regular surface of discontinuity for some or all of the electro-mechanical fields, which separates Ω into Ω − (t) and Ω + (t). For x ∈ S(t) we denote by n(x, t) the unit normal to S(t) pointing into Ω + .
Classically, one derives the field equations in Ω \ S and jump conditions on S by localizing the integral form of mechanical balances and electrostatic equations on these two sets (Abeyaratne and Knowles, 1990;Jiang, 1994b;Gurtin et al., 2009;Kessler and Balke, 2006;Mueller et al., 2006). In the following, we directly use the local forms. Mechanical equilibrium reads div σ = 0 in Ω \ S(t), σ n = 0 on S(t), (2.2) 4 where σ(x, t) denotes the Cauchy stress tensor. For a generic field w(x, t) defined in Ω \ S we define the jump w(x, t) = w + (x, t) − w − (x, t), where w + (x, t) and w − (x, t) are the limiting values of w at point x on S(t), given by When it comes to the electrostatic fields, we denote by e(x, t) the electric field, by d(x, t) the electric displacement, and by p mat (x, t) the polarization 2 field. All three are defined over R 3 × R and related by definition through where 0 is the permittivity of vacuum. In the absence of free charges, Gauss' law becomes while the Maxwell-Faraday equation in the quasistatic regime yields curl e = 0 in Ω \ S(t), e × n = 0 on S(t). (2.6) In the context of ferroelectrics, surface S typically represents domain walls (separating domains with different polarizations) or interphases (separating different phases, e.g., tetragonal and rhombohedral phases). We use the index α ∈ D = {1, . . . , N } to identify the N different domains or phases. By taking e and ε as the independent state variables, the constitutive relations for σ and d derive from the electric enthalpy density W α (e, ε) associated with domain α through, respectively, σ = ∂W α ∂ε and d = − ∂W α ∂e . (2.7) We associate with the interface S a constant excess electric enthalpy density γ, more simply referred to as interfacial energy. The combination of energy balance and entropy imbalance permits deriving the rate of entropy production δ(x, t) per unit area of S at x ∈ S(t) as where V n (x, t) is the normal velocity of S (taken positive in the direction of n), and f (x, t) is the thermodynamic driving traction (Jiang, 1994b;Mueller et al., 2006) 3 given by f = n · C n + γκ. (2.9) In (2.9), κ(x, t) denotes twice the mean curvature of S and C(x, t) is the Eshelby tensor provided by 2 Note that we here use the notation pmat for the polarization (sometimes referred to as material polarization, see e.g., Schrade et al. (2013)). Indeed, we reserve the notation p for what we call the spontaneous polarization, i.e., the polarization in a ferroelectric variant at zero stress and electric field (see Section 2.2).
3 Note that in Jiang (1994a) the driving traction at the interface is expressed with a different free-energy density, which reads, using our notations, ∼ W α(pmat, ε) with associated constitutive relations σ = ∂ ∼ W α/∂ε and e = ∂ ∼ W α/∂pmat. It is related to the electric enthalpy Wα(e, ε) through Wα(e, ε) = ∼ W α(pmat, ε) + 0 2 e · e − e · d. Rewriting the driving force (4.8) in Jiang (1994a) in terms of Wα, one obtains with the help of (2.2)-(2.6) the first term of (2.9). Alternatively, this expression has been derived variationally by Mueller et al. (2006). The term γκ in (2.9), is associated with curvature-driven motion and originates from the fact that S is endowed with excess electric enthalpy. Derivations of this classical term can be found, e.g., in Gurtin (2000) for surfaces in three dimensions or in Gurtin (1993) for curves in two dimensions.
in the domain α, with I the identity tensor.
Notice from (2.9) that f and V n are work-conjugate. The motion of S(t) is hence specified by a relation between V n and f , called the kinetic relation. The latter may be interpreted as an additional constitutive relation, which embeds information about the microscopic processes underlying the motion of S and takes the form V n =V n (f ), (2.11) whereV n : R → R is a function subject to the restriction which ensures the positivity of the entropy production rate δ (required by the second law of thermodynamics). In general, several surfaces of discontinuity between the different domains or phases evolve simultaneously.
Interface energy γ and kinetic relationV n (f ) are characteristic of each interface S, i.e., they generally depend on the two phases separated by S.

Constitutive laws
We apply the model presented in Section 2.1 to the case of a ferroelectric ceramic, which exhibits, below the Curie temperature, a tetragonal crystalline structure with N = 6 variants in three dimensions (3D), as is the case, e.g., for barium titanate (BaTiO 3 ) and some compositions of lead zirconate titanate (PZT).
Bulk constitutive behavior. Each variant possesses a spontaneous polarization p α which, in the orthonormal system (ẽ 1 ,ẽ 2 ,ẽ 3 ) aligned with the principal crystallographic directions, is defined by where p 0 is the magnitude of the spontaneous polarization of the material. Due to the distortion between the cubic and tetragonal phases, each variant exhibits a spontaneous strain ε s α , which we write as (Kamlah, 2001) where ε a < 0 and ε c > 0 are, respectively, the spontaneous longitudinal strain along the a-axes (orthogonal to the direction of spontaneous polarization) and along the c-axis (along the direction of spontaneous polarization) 4 . With these features, each ferroelectric variant is modeled as a linear elastic, dielectric, and piezoelectric material with the following electric enthalpy: where α , C α , and E α are the second-order dielectric tensor, fourth-order elasticity tensor, and third-order piezoelectric tensor (Brainerd, 1949), respectively, which satisfy the material symmetries of the tetragonal phase. The constitutive relations resulting from (2.7) are The dynamics of domain walls has been investigated both experimentally and theoretically in bulk singlecrystals since the 50's and in epitaxial ferroelectric thin films since 2000. These works, reviewed in Appendix A, show that ferroelectrics exhibit two regimes of domain wall dynamics. At low electric fields, the domain wall velocity is described by an inverse exponential law, where V ∞ is a characteristic velocity and e a an activation field which varies as the inverse of temperature. This kinetics has been interpreted as resulting from thermally-activated nucleation and growth of new domains along the domain wall (Miller and Weinreich, 1960;Shin et al., 2007;Liu et al., 2016) or as a creep process of an interface moving in a pinning potential (Tybell et al., 2002;Jo et al., 2009). As the electric field increases above a threshold field e t , the wall kinetics experiences a transition to a different regime described by a power law, where θ is an exponent that depends on the dimensionality of the system: θ = 1.4 has been reported for bulk BaTiO 3 (Stadler and Zachmanidis, 1963), while θ = 0.7 was found for PZT epitaxial thin films (Jo et al., 2009). Different theoretical explanations for that second regime have been proposed. On the one hand, in line with the work of Miller and Weinreich (1960), this second regime is interpreted in terms of the nucleation of domains at the domain boundary with thicknesses of multiple crystalline unit cells (Stadler and Zachmanidis, 1963). On the other hand, in the picture of an interface moving in a pinning potential, the change in kinetics is seen as a pinning-depinning transition to a non-activated regime, i.e., the wall moves without assistance of thermal fluctuation (Jo et al., 2009). Relations (2.17) and (2.18) were obtained on plate-like samples with out-of-plane polarization and electric field. In this configuration, the driving traction reduces to f = 2p 0 e, i.e., it is proportional to the applied electric field e. Therefore, the experimental evidence discussed above indicates that the motion of domain walls is governed by the following kinetic relation: where V l , V h , f a , f t and θ are referred to as the domain wall kinetics parameters and satisfy V l exp(−f a /f t ) = V h to ensure the continuity ofV n at f t . Such a kinetic relation is plotted in Figure 1 with the parameters corresponding to 180 • walls in BaTiO 3 . As briefly reviewed in Appendix A, experiments indicate that the kinetic relation for 90 • domain walls differs from that of 180 • walls either in the coefficients of (2.19) or in the functional form ofV n (f ), depending on the ferroelectric material.

Specialization to rigid ferroelectrics
Our goal is a phase-field model that displays the nonlinear kinetics presented in Section 2.2, in contrast to the linear kinetics that the classical Allen-Cahn equation furnishes. To this end, we make two simplifying assumptions: first, the ferroelectric ceramic is assumed mechanically rigid, so that mechanical couplings are neglected, and, second, the permittivity is assumed isotropic. These assumptions allow us to work in a clean setting, where the new features of the general kinetics model are clearly revealed. As they are not necessarily valid in practice, they will be lifted in future work, when the objective is to provide quantitative predictions. In this simplified setting, the electric enthalpy density (2.15) of each phase α reduces to that of a dielectric, W α (e) = − 2 e · e − p α · e, (2.20) where the scalar permittivity is independent of the phase. Consequently, the constitutive relation (2.16) 2 reads d = e + p α . (2.21) In view of (2.5) 2 , (2.6) 2 and (2.15), the driving traction (2.9) is rewritten in the simplified form where e = (e + + e − )/2.

A phase-field model based on the Allen-Cahn equation
Given the numerical difficulties related to the implementation of sharp-interface models, the phase-field method has been a method of choice for modeling the evolution of domains in ferroelectric ceramics. Classically, the order parameter is chosen as the material polarization p mat (x, t) (see e.g., Zhang and Bhattacharya, 2005a;Su and Landis, 2007;Vidyasagar et al., 2017) or the spontaneous polarization p(x, t) (Schrade et al., 2013(Schrade et al., , 2014, and evolves according to the Allen-Cahn equation.
Anticipating that the phase-field model with general kinetics-introduced in Section 4-uses p as the order parameter, we introduce a phase-field model based on the Allen-Cahn equation written in terms of p, with features of the sharp-interface model of Section 2.3. It serves as a regularized benchmark model to later clarify the new features of the general kinetics model. In particular, we show that the Allen-Cahn equation yields a linear kinetic relation in the limit of low applied electric fields and that the regularization length has an upper bound for properly approximating the sharp-interface model.
For notational clarity, we use a hat (ˆ) to highlight all quantities (regularized energy, domain wall properties, and kinetics) associated with the Allen-Cahn-based phase-field model as opposed to their counterparts for the general kinetic model later discussed in Section 4. We formulate the model in Section 3.1 and characterize the properties (interface width, energy, and kinetics) of 180 • and 90 • domain walls in Section 3.2 before a brief summary in Section 3.3.

Model formulation
A regularized electric enthalpy. The phase-field model is based on a regularized versionΨ(e, p, ∇p) of the electric enthalpy (2.20), which induces continuous variations of p across domain walls. We here adopt the electric enthalpy density of Schrade et al. (2013) specialized to rigid ferroelectrics, which is given bŷ whereĈ s andĈ g are constants, andψ s (p) is a multi-well separation potential with minima at the six spontaneous polarization states p α . Letting p = p/p 0 be the dimensionless spontaneous polarization,ψ s is defined aŝ being a dimensionless sextic polynomial. As shown in Schrade et al. (2013), the parameters {a i } i=0...4 are uniquely determined by requiring (i) that the minima of ψ s are located at {p α } α=1...6 , (ii) that the 180 • barrier is normalized, i.e., and (iii) by specifying the relative height 0 <ĥ 90 < 1 and locationχ of the 90 • barrier through 5 ψ s (χ(ẽ 1 +ẽ 2 )) =ĥ 90 and ∂ψ s ∂p 1 (χ(ẽ 1 +ẽ 2 )) = 0. Governing equations. In the regularized phase-field model all fields are continuous over B, so the differential forms of Gauss' law (2.5) 1 and Faraday's equation (2.6) 1 apply over all Ω. Akin to (2.7), the electric displacement derives fromΨ as (3.5) These equations are supplemented by an evolution law for p, which is usually taken as the Allen-Cahn gradient-descent equation: In view of (3.1), (3.6) becomes A numerical characteristic electric field. In the sharp-interface model, the spontaneous polarization takes the discrete values p α corresponding to the different domains or phases. In the regularized model, based on the electric enthalpy (3.1), the spontaneous polarization p is a continuously varying variable that plays the role of a vector phase field. Away from the interface, the Laplacian term in (3.7) approximately vanishes and equilibrium values of p are defined by Whereas large departures of p from the values p α are normal in the transition regions of domain walls, the spontaneous polarization p shall remain close to one of the p α inside the domains. In view of (3.8), this is the case under the condition that |e| ê c , (3.9) whereê c =Ĉ s /p 0 is a numerical characteristic electric field. In practice, (3.9) is ensured by choosing a sufficiently largeĈ s . As we will see in Section 3.2, that is equivalent to selecting a sufficiently small regularization length. Note that the fulfillment of (3.9) prevents the nucleation of new domains of polarization from (3.7). Therefore, domain nucleation, which is frequently observed in experiments, must be included through an explicit additional condition.

Analytical solutions for 180 • and 90 • domain walls
In this section, we compute analytical solutions for straight 180 • and 90 • domain walls, both in equilibrium and in steady-state motion. When adopting the form of traveling waves, these solutions provide relations between the numerical parametersĈ s ,Ĉ g , and µ of the phase-field model and the properties of domain walls (width, interfacial energy, and kinetic relation).

Equilibrium profile of a neutral 180 • domain wall
We consider the one-dimensional (1D) case of a 180 • domain wall with spontaneous polarization p(x, t) = p 2 (x 1 , t)ẽ 2 varying from the equilibrium close to p 0ẽ2 at x 1 → −∞ to the one near −p 0ẽ2 at x 1 → +∞, as schematically shown in Figure 2(a). This configuration is referred to as neutral 180 • domain wall due to the continuity of the normal component of p across the wall, which implies the absence of bound charges in the wall. In view of (3.5) and noting that the polarization profile is divergence-free, we can assume a uniform electric field e = eẽ 2 while satisfying Gauss' law (2.5) 1 .
We look for a traveling wave solution of the form where ξ = x 1 − v 180 t is the co-moving coordinate and v 180 the velocity of the wall. We denote byψ 180 (p 2 ) = ψ s (p 2ẽ2 ) the section of the separation potential relevant for this situation. Inserting (3.10) into (3.7) yields an ordinary differential equation forp 2 (ξ): where (3.11) 2 are the far-field boundary conditions at ±∞, and the prime denotes the derivative of a singlevariable function with respect to its variable.
The equilibrium profile: e = 0. In the absence of an electric field, (3.11) admits a static analytical solution with v 180 = 0 (Cao and Cross, 1991): Based on this equilibrium profile, we relateĈ s andĈ g to the interface width and energy. We define the widthl of the 180 • domain wall (see Figure 2) aŝ l = 2/p eq 2,ξ (0) (3.13) and its interfacial energy asΓ is the sharp-interface solution of the 180 • domain wall. Noting that e is uniform and equal in both the phase-field and sharp-interface formulations and that the profilep eq (ξ) is symmetric, (3.14) reduces tô The first integral of (3.11) in the equilibrium case (e = 0, v 180 = 0), obtained by multiplying (3.11) byp eq 2,ξ and integrating from ξ = −∞ to some arbitrary ξ, yields 2Ĉ sψ180 p eq 2 (ξ) =Ĉ g p 2 0 p eq 2,ξ (ξ) 2 .
( 3.17) From (3.17) we extractp eq 2,ξ (ξ) = 2Ĉ sψ180 p eq 2 (ξ) /(Ĉ g p 2 0 ), which we exploit for the change of variable ξ → p 2 in (3.16). This yieldsΓ a 0 + a 1 p 2 + a 2 p 4 + a 4 p 6 dp is a numerical coefficient. Further, combining (3.17) and (3.13) lets us express the widthl aŝ Conversely, (3.18) and (3.19) allows us to determineĈ g andĈ s fromΓ andl viâ (3.20) We note here that the interfacial energyΓ shall be taken as the physical interfacial energy γ introduced in the sharp interface model, whilel is viewed as a numerical parameter. Indeed, we take the perspective that the objective of the phase-field model is to properly account for the evolution of the pattern formed by domains or phases. In particular, the accuracy of the polarization profile within the diffuse interface (as compared to the physical profile) is of minor importance; instead, we aim to accurately capture the time evolution of the trace of domain walls. An accurate value ofΓ to represent the interfacial energy is important to the extent that it is expected from the sharp interface model (see (2.9) and (2.19)) to affect the evolution of a curved interface. This latter point can be confirmed numerically. By contrast,l does not directly affect the evolution of the domain wall, but its choice results from a compromise between the following points: • Because a diffuse interface needs sufficient numerical resolution (e.g., it should be discretized with typically no less than three points across its width in FFT-based schemes (Vidyasagar et al., 2017)), l cannot be taken too small.
• On the other hand,l cannot be taken too large either, because the characteristic electric fieldê c = C s /p 0 =Γ/(ηlp 0 ) introduced in Section 3.1 must be sufficiently large for (3.9) to be satisfied with the physical electric field (notably to ensure that p remains close to one of the p α inside domains).

Equilibrium profile of a charged 180 • domain wall
The other characteristic 180 • domain wall is the so-called charged wall, in which the polarization is orthogonal to the domain wall (head-to-head or tail-to-tail). This configuration is usually not considered in the literature (see e.g. Cao and Cross, 1991;Schrade et al., 2013;Flaschel and De Lorenzis, 2020), because in the absence of free charges it is associated with high electric fields (of the order of p 0 / ). As such, it implies high energies (see (3.1)) and is hence considered unstable in general. While it is true that ferroelectric domains evolve in such a way as to minimize the existence of charged 180 • domain walls, these cannot be completely disregarded. Indeed, charged portions of domain walls exist if one considers, e.g., two domains of antiparallel polarization separated by a closed boundary (see Figure 3). This is the natural scenario of a newly nucleated domain in its parent domain.
Let us consider the 1D configuration of a charged 180 • domain wall with spontaneous polarization p(x, t) = p 1 (x 1 , t)ẽ 1 varying from the equilibrium close to p 0ẽ1 at x 1 → −∞ to the one near −p 0ẽ1 at x 1 → +∞ (see Figure 2(b)). We assume that the applied electric field is zero in theẽ 2 -direction, which implies, e(x 1 , t) = e 1 (x 1 , t)ẽ 1 . As a result, the electric displacement is along theẽ 1 -direction and, because of (2.5) 1 , it is uniform. Assuming further that it is time-independent, we write d(x 1 , t) = dẽ 1 where d is constant. This allows us to write the non-zero electric field component as (3.21) We look for a traveling wave solution of the form where ξ = x 1 − v c 180 t is the co-moving coordinate, v c 180 the velocity of the charged wall, andp 1 (ξ) is a smooth, differentiable function. Inserting (3.22) into (3.7) and resorting to (3.21) yields the ordinary differential equation forp 1 (ξ),

Velocity of 180 • domain walls
We proceed to consider the steady-state velocity of 180 • domain walls, beginning with neutral ones and later concluding the solution for charged ones. Taking the first integral of (3.11) 1 , we express the velocity v 180 of the domain wall as which involves the unknown functionp 2 (ξ). Under the assumption that (3.9) is satisfied, the traveling wave problem (3.11) with e = 0 is a perturbation of the case e = 0. Hence, to first order in e/ê c the velocity (3.26) can be approximated using the equilibrium profilep eq 2 , which yields . We point out that the above analysis, of course, presumes that stable motion of a domain wall at constant speed is feasible. While this is realistic for neutral 180 • domain walls, the high energy associated with their charged counterparts precludes their stable motion over significant times in general. 6 For an analogous derivation with mechanical coupling accounted for, see Indergand et al. (2020).

Equilibrium profile of a neutral 90 • domain wall
For studying 90 • domain walls, the basis of spontaneous polarizations {p α } α∈(1,6) defined in (2.13) and the multi-well potential (3.2) are rotated by +π/4 around the directionẽ 3 . We consider a 1D head-to-tail 90 • domain wall with spontaneous polarization p(x, t) = p 2 (x 1 , t)ẽ 2 + p 1 (x 1 , t)ẽ 1 , varying from the equilibrium close to p 0 (ẽ 1 +ẽ 2 )/ √ 2 at x 1 → −∞ to the one near p 0 (ẽ 1 −ẽ 2 )/ √ 2 at x 1 → +∞ (see Figure 4(a)). When looking for an equilibrium or a traveling wave solution of (3.7), the p 1 -component appears to vary across the domain wall. As a result, for Gauss' law (2.5) 1 to be satisfied, the electric field is non-uniform and the polarization profile does not admit an analytical solution.

(3.30)
Here and in the following, superscript * indicates quantities associated with 90 • domain walls.

Equilibrium profile of a charged 90 • domain wall
The charged 90 • domain wall corresponds to the head-to-head (or tail-to-tail) configuration (see

Velocity of 90 • domain walls
Following the same procedure as for 180 • domain walls, the first integral of (3.28) allows us to write v 90 as an expression analogous to (3.26), where one substitutesψ 90 forψ 180 . In the regime e ê c ,p 2 can be approximated by the equilibrium profile of (3.28), which yields whereĉ 90 =l/(2η * µp 2 0 ) and f 90 = √ 2ep 0 denote, respectively, the kinetic coefficient and driving traction associated with the 90 • domain wall. Under the condition (3.25), the charged 90 • domain wall shows the same linear kinetics, with coefficientĉ 90 , as its neutral counterpart.

Summary of the properties of 180 • and 90 • domain walls
Interface energy and width. We have computed exactly the interfacial energyΓ and widthl of a neutral 180 • domain wall and found that these can be arbitrarily prescribed by choosingĈ g andĈ s according to (3.20). While we expect to use the physical domain wall energy forΓ,l is viewed as a numerical parameter to be chosen in light of considerations discussed above (the more stringent condition onl being (3.25)). The charged 180 • domain wall has approximately the same equilibrium properties as the neutral one, as long as (3.25) is satisfied.
For the classical head-to-tail 90 • domain wall, we have derived an approximation of its interfacial energŷ Γ 90 and widthl 90 , whose accuracy is best forχ = 0.5. The ratioΓ 90 /Γ can be set according to its physical value by independently settingĥ 90 . Indeed, (3.32) indicates that this ratio is equal to that of the integrals over the 90 • and 180 • barriers in the multi-well potential ψ s (p). Finally, the widthl 90 is automatically determined by (3.33) and cannot be set independently. Note that, because we viewl 90 as a numerical regularization length, there is no need to prescribe it independently. However, where it is viewed as the physical width of a 90 • domain wall, the phase-field framework can be modified such thatl 90 is chosen in agreement with its physical value (Flaschel and De Lorenzis, 2020).
Kinetics. We have found that, in the regime |e| ê c , the Allen-Cahn evolution equation (3.6) implies approximately a linear kinetic relation between the velocity of domain walls and the driving traction. In addition, the 180 • and 90 • domain walls have different kinetic coefficients related by the ratioη * /η. Whereas the choice of inverse mobility µ allows us to independently set one kinetic coefficient, the second one is automatically determined by the ratio of interfacial energies (see (3.32) and (3.35)).
This linear kinetics is not representative of the complex nonlinear kinetics of domain walls, which experiments and atomic-scale modeling discussed in Section 2.2 show. Yet, a proper account of the kinetics of domain walls is important in two respects. First, different kinetics yield different microstructure evolution 7 , hence a proper account of the kinetics is needed for obtaining an accurate description of switching at the mesoscale. Second and more importantly, rate-dependent effects in the macroscopic footprints of polarization switching-such as those discussed in Section 1-are a direct consequence of the kinetics of domain walls. Having a phase-field model that permits domain wall motion with nonlinear kinetics requires us to revise the evolution equation (3.6). This is the subject of the general kinetics model introduced in Section 4.

General kinetics model
In this section, we introduce a new phase-field model for ferroelectrics, which regularizes the sharpinterface model of Section 2 while conserving nonlinear and independent kinetics for 180 • and 90 • domain walls. We first formulate the model in Section 4.1, before computing analytical traveling wave solutions for 180 • and 90 • domain walls in Section 4.2. The main properties of straight domain walls are summarized and compared to those of the Allen-Cahn model in Section 4.3. A numerical implementation of the present model will be performed in Part II, which must also account for the special behavior of curved domain walls and triple and quadruple points at which multiple phases meet.

Model formulation
A regularized electric enthalpy. In the general kinetics model, we use a multi-phase field ϕ(x, t) ∈ (0, 1) N , where N is the number of domains, and for tetragonal ferroelectrics we use N = 6. Let ϕ α denote the volume fraction of domain α, so that the spontaneous polarization is defined as (4.1) We introduce the electric enthalpy Ψ, which regularizes W (e, p), as where C s and C g are two positive constants, and ψ s is a multi-well potential with minima at (or in the neighborhood of) the six values of ϕ associated with the spontaneous polarization states: ϕ = (1, 0, . . . , 0) and its permutations. There exists several possibilities for the choice of ψ s , and the optimal one depends on the values of the physical parameters. Here, we discuss two simple forms for ψ s that admit analytical insight, even though the framework is sufficiently general to allow for other choices as well. For analytical derivations, we introduce where h 90 denotes the height of the 90 • barrier relative to the 180 • barrier, C t > 0 is a numerical parameter and T = {(α, β, γ) : α, β, γ ∈ D, α = β = γ}, and we consider n = 1 and n = 2. As we shall see below, ψ (1) s has the advantage that the effective potential for charged walls has minima that remain at the spontaneous polarization states but the drawback of presenting discontinuities in its derivatives. On the other hand, with ψ (2) s discontinuities in the derivatives are absent, but the zeros and minima of the effective potential for charged walls are shifted from the spontaneous polarization states.
The contribution τ (ϕ) to ψ (n) s (ϕ) serves to avoid the spurious appearance of a third phase in a domain wall by penalizing the simultaneous occurrence of more than two phases. 8 Despite the apparently complicated expression of ψ where G αβ = G βα : R → R is the kinetic function associated with the transformation between phases α and β, and h β→α denotes the relative variational derivative of Ψ, (4.8)

Inserting (4.2) in (4.8) yields
h β→α e, ϕ, For tetragonal ferroelectrics, we distinguish between two types of transformations: the motion of 180 • and 90 • domain walls. Hence, we take the function G αβ as G 180 for {α, β} ∈ I 180 and G 90 for {α, β} ∈ I 90 . By the second law of thermodynamics, G 180 and G 90 must be odd functions that satisfy the constraints (4.10) Indeed, as we will show in Section 4.2, these two functions specify the kinetic relations of 180 • and 90 • domain walls in the sense of (2.11) in the sharp-interface model. Condition (4.10), like (2.12), ensures the positivity of the dissipation in the phase-field model. Equation (4.7) is used with an initial condition ϕ(x, 0) that satisfies 10 6 α=1 ϕ α (x, 0) = 1 for all x ∈ Ω. (4.11) Using that G 180 and G 90 are odd functions, one can show that (4.7) guarantees 6 α=1φ α = 0, whereby (4.11) holds for all times if it holds for the initial condition.
Comparison with other phase-field models with nonlinear kinetics. The evolution equation (4.7) is an extension to multiple phases of the hybrid model introduced by Alber and Zhu (2013) in the setting of stress-driven phase transformations between two phases. In particular, in a domain wall between two phases, say α = 1 and β = 2, for which ϕ 2 (x, 0) = 1 − ϕ 1 (x, 0) and ϕ γ (x, 0) = 0 for γ = 3 . . . 6, (4.7) reduces to (4.12) (4.12) 1 is akin to Equation (1.3) of Alber and Zhu (2013) (there formulated in the context of stress-driven transformations). Note that the double-well potential proposed by Alber and Zhu (2013) is a variation of (4.5).
Similarly, (4.12) 1 is reminiscent of the evolution equation (2.7) of Agrawal and Dayal (2015a), which derives from a conservation law for the number of interfaces. However, in their work, the free-energy density was regularized differently from ours and that of Alber and Zhu (2013), without resort to a double-well potential but through a regularized Heaviside function. Further, Agrawal and Dayal (2015a) included a source term in the evolution equation for ϕ to account for the nucleation of new domains. In our case, in the absence of such a source term, (4.7) only propagates domain walls, while for nucleation new domains shall be introduced explicitly according to an appropriate criterion, which we do not discuss here.

Analytical solutions for 180 • and 90 • domain walls
In this section, we compute analytical traveling wave solutions for the propagation of straight 180 • and 90 • domain walls as obtained from the general kinetics model. In doing so, we follow in one dimension a procedure outlined, e.g., in Fried and Gurtin (1994) for Allen-Cahn-type models. In particular, we demonstrate the nonlinear kinetics of domain walls that the general kinetics model permits, and we compare our findings to the Allen-Cahn model.
By setting the origin toφ(0) = 1/2, (4.22) forms an initial value problem forφ(ξ), starting at ξ = 0 and to be solved for both increasing and decreasing ξ. Its integration yieldš for ξ < −πl/4, cos 2 π 4 − ξ l for − πl/4 ≤ ξ ≤ πl/4, 1 for ξ > πl/4, With regards to kinetics, in (4.21) 2 2p 0 e = f corresponds to the driving traction (2.22) of the straight, neutral 180 • domain wall in the sharp-interface model. Hence, G 180 directly provides the kinetic relation for 180 • domain walls (i.e., (2.11) in the sharp-interface model). Therefore, the general kinetic model allows us to choose the kinetic relation arbitrarily (in the set of monotonously increasing odd functions that satisfy (4.10)). In particular, for modeling ferroelectrics, a physics-based choice for G 180 is to adopt the functional form (2.19) derived from experiments. In addition, note that for the neutral 180 • domain wall, the kinetic relation is exactly given by (4.21) 2 for all e, whereas with the Allen-Cahn model the linear kinetic relation (3.27) is only an approximation valid under the condition (3.9).

Charged 180 • domain wall profile and kinetics
We consider the same charged domain wall as in Section 3.2.2, i.e., p(x 1 , t) = p 1 (x 1 , t)ẽ 1 , varying between p 1 = p 0ẽ1 at x 1 → −∞ and p 2 = −p 0ẽ1 at x 1 → +∞ (see Figure 5(c)). For the same reasons as those invoked in Section 4.2.1, only ϕ 1 and ϕ 2 are non-identically zero, and they satisfy ϕ 1 + ϕ 2 = 1. Using the single phase field ϕ = ϕ 2 , we rewrite the polarization as Further, with the same assumptions on the electric field and electric displacement as in Section 3.2.2, the electric field is given by (3.21). We write the evolution equations for ϕ 1 and ϕ 2 as an equation for ϕ, (4.29) By combining (4.9), (4.28) and (3.21), h 1→2 becomes Parameter ζ is equivalent toζ involved in (3.25) in the Allen-Cahn model and represents the ratio of the electric field p 0 / induced by Gauss' law in a charged 180 • domain wall to the numerical characteristic electric field e c = Γ/(ηlp 0 ). As is apparent from Figure 6, the effective potential for n = 1 is indeed a double-obstacle potential as long as ζ < 2 with minima and zeros coinciding at ϕ = 0 and ϕ = 1. In contrast, for n = 2, minima of the potential are shifted inwards and we denote by ϕ m (ζ) the smallest minimizer of ψ c 180 . We look for a traveling wave solution of (4.29) of the form while assuming, without loss of generality,φ ,ξ ≥ 0. Inserting (4.33) into (4.29), we obtain the equation for is substituted for l and the interfacial energy is replaced by (4.37) While we see from (4.37) that the regularization introduces a modification of the interfacial energy between neutral and charged domain walls, this turns out not to be an issue in practice. From the point of view of the sharp-interface model, an accurate interfacial energy is important, as it contributes to the driving force f through the effect of curvature in (2.22). In practice we have ζ < 2, as required to have a double-obstacle potential. If ζ is low (typically ζ < 0.2) the change in interfacial energy remains less than 10%. Therefore we expect an accurate account of the curvature contribution in (2.22). By contrast, when ζ is large (typically ζ = 1 − 1.5), the artificial change in interfacial energy becomes larger (up to 50% relative change for ζ = 1.5). However, when ζ is large, Γ is small (see (4.32)), and the contribution γκ to the driving force (2.22) becomes negligible -compared to that of the bulk energies-almost everywhere except in local zones of high curvature. Hence, in this case the regularization introduces significant variations in Γ, but these are of little importance. In practice, the contribution of curvature of domain walls to the evolution of domains is mostly negligible in ferroelectric switching as well as in other solid-solid phase transformations. Therefore, the fact that the general kinetic model with n = 1 allows us to take large values of ζ up to 1.5 is a significant advantage over the Allen-Cahn model and over general kinetic model with n = 2 (discussed below), which both require ζ 1 (typically ζ < 0.2).

90 • domain walls: profile and kinetics
Neutral domain wall. As in Section 3.2.4, the basis of spontaneous polarizations {p α } α∈(1,6) defined in (2.13) is rotated by +π/4 around theẽ 3 -axis. Consider, e.g., the same 1D head-to-tail domain wall as in Section 3.2.4, i.e., between p 1 = p 0 / √ 2(ẽ 1 +ẽ 2 ) at x 1 → −∞ and p 4 = p 0 / √ 2(ẽ 1 −ẽ 2 ) at x 1 → +∞ (see Figure 7(b)). This wall involves ϕ 1 and ϕ 4 satisfying ϕ 1 + ϕ 4 = 0, while the other phase fields are identically zero. For simplicity we define ϕ = ϕ 4 . A specialty of the polarization profile obtained from the general kinetics model is that the rotation of the polarization vector across the wall is predetermined by the decomposition (4.1), i.e., * whereby p rotates with a constant component alongẽ 1 . This is in contrast to the rotation of the polarization vector resulting from the Allen-Cahn model, which depends on the multi-well potential ψ s given by (3.2), in particular on the parameterχ introduced in (3.4), which determines the location of the 90 • barrier.
As a consequence, whereas in Section 3.2.4 the computed properties (wall width and interfacial energy) of the neutral 90 • domain wall are only approximations relying on the assumption that the p 1 -component remains constant (a reasonable assumption forχ = 0.5), for the general kinetics model these same properties characterize the neutral 90 • domain wall exactly.
Under the application of a uniform electric field e = eẽ 2 and with a derivation in all points analogous to that developed in Section 4.2.1, we conclude that the 90 • domain wall can be described by a traveling wave solution, where √ 2p 0 e corresponds to the driving traction (2.22) associated with the 90 • configuration under consideration. Again, the polarization profile is given by (4.26) and (4.27) with a width l 90 = l/ √ h 90 replacing l and an interfacial energy Γ 90 = √ h 90 Γ, which can be set to its physical value through a proper choice of h 90 .
Charged domain wall. The transition from the neutral to the charged 90 • domain wall (see Figure 7(c)) shows the same features as the analogous transition for the 180 • domain wall. We discuss specifically the case n = 1, for which the kinetic relation is exactly given by the function G 90 , provided that ζ < 2 √ 2h 90 is satisfied (like ζ < 2 for 180 • walls). The modified values of interfacial width and energy due to the regularization with n = 1 are Interface energy and width. We have computed exactly the interfacial energy and width of the different types of domain walls obtained with the general kinetics model. The expressions are analogous to those of the Allen-Cahn model. In particular, the choice of C s and C g allows us to set independently the interfacial energy and width (i.e., the regularization length) of 180 • domain walls, while h 90 determines the ratio of interfacial energy of 90 • and 180 • domain walls. Furthermore, the expressions for the interfacial properties of walls remain exact under an applied electric field, in contrast to the Allen-Cahn model for which one only has approximate expressions of interfacial properties when domain walls are not in equilibrium. This is notably related to the fact that, with the general kinetics model, values of the spontaneous polarization within domains are unchanged in the presence of an electric field. Lastly, a special feature of the general kinetics model is that for 90 • domain walls the rotation of the polarization occurs with one of its components kept constant, as predetermined by the decomposition (4.1). This is in contrast to what happens with the Allen-Cahn model, where the rotation of p depends on the location of the 90 • barrier, given byχ in the multi-well potential.
Kinetics. We have found that for straight domain walls the general kinetics model allows us to prescribe independently the kinetics of 180 • and 90 • domain walls and that these kinetic relations can be nonlinear. The kinetic relations (2.11) between the velocity of the interface (velocity of the traveling wave solutions) and the driving traction (computed from the values of the fields away from the diffuse interface) are defined through the functions G 180 and G 90 for 180 • and 90 • domain walls, respectively.
Interestingly, in the general kinetics model, the kinetic relation of straight walls is exactly given by G 180 and G 90 irrespective of the regularization length l, provided that ζ < 2 and ζ < 2 √ 2h 90 are satisfied (i.e, l is smaller than a critical length determined by the material parameters). This is in contrast to the Allen-Cahn model, for which the accuracy of the linear kinetic relations (3.27) and (3.35) for straight 180 • and 90 • domain walls decreases as the regularization length l increases (because, as l increases with constant material parameters, e c in (3.9) decreases, which reduces the accuracy of (3.27) and (3.35)).
As will be discussed in Part II of this study, when it comes to curved domain walls, the general kinetics model allows us to introduce the kinetic relations via G 180 and G 90 in an exact sense only in the limit of l → 0 (see also the asymptotic analysis of Alber and Zhu (2013) for the phase-field model that inspired the present general kinetics model).

Conclusion
We have formulated a phase-field model for ferroelectrics that permits domain wall evolution with nonlinear kinetics, which we term the general kinetics model. Nonlinear domain wall kinetics is a crucial feature for properly accounting for rate effects in ferroelectric switching and is nonetheless absent from existing diffuseinterface ferroelectrics models. Seen in the general class of phase-field models for structural transformations (which includes switching between ferroelectric variants, twinning, solid-solid phase transformations), the general kinetics model is, to the best of our knowledge, the first diffuse-interface model with arbitrary nonlinear interface kinetics formulated for systems with multiple (i.e., more than two) phases. It has been devised as an extension to multiple phases of the hybrid model of Alber and Zhu (2013) (adapted here to the setting of ferroelectric switching).
To shed light on the new features of the general kinetics model, as compared to the classical Allen-Cahnbased phase-field model, we have worked in the simplified setting of rigid ferroelectrics, an assumption which only accounts for the electrostatic contribution to domain evolution and neglects the mechanical one. 12 In this setting, we have compared two regularized models: the classical Allen-Cahn-based phase-field model and the newly introduced general kinetics model. By computing analytical traveling wave solutions for the propagation of straight 180 • and 90 • domain walls, both neutral and charged, we have established the following properties of the two phase-field models: 1. Interfacial width and energy: For both diffuse-interface models, the numerical parameters associated with the regularization of the electric enthalpy permit to set independently the interfacial energy of 180 • and 90 • domain walls and the width of 180 • domain walls. While the former are material parameters characterizing the walls, we see the latter as a numerical regularization length, whose choice is limited by practical constraints. With the Allen-Cahn model, interfacial width and energy weakly depend on the applied electric field and remain approximately constant in the regime (3.9) of |e| ê c , withê c =Γ/(ηlp 0 ) a numerical characteristic electric field dependent on the regularization lengthl. By contrast, with the general kinetics model, the interfacial width and energy of neutral domain walls are constant and independent of the applied electric field. As for charged domain walls (which encounter nonuniform electric fields induced by Gauss' law), both models require, for wall properties to remain approximately independent of the regularization length, that the magnitude of the electric field associated with Gauss' law remains small compared to a characteristic electric field. This translates into an upper bound onl and l, which reads asηlp 2 0 /( Γ ) 1 (see (3.25)) for the Allen-Cahn model and ηlp 2 0 /( Γ) 1 for the general kinetics model (with quantitative estimates given in Sections 4.2.2 and 4.2.3).
2. Kinetics of domain walls: For the motion of straight 180 • and 90 • domain walls, we have shown that the Allen-Cahn formulation furnishes a regularization of a sharp-interface model with a linear kinetic relation between wall velocity and associated driving traction. In addition, while the kinetic coefficient associated with one type of domain walls (e.g., 180 • domain wall) can be freely set with a proper choice of the inverse mobility µ in (3.6), that of the other type of wall (e.g., 90 • domain wall) is automatically prescribed by the ratio of interfacial energies of the two types of domain walls. Furthermore, the Allen-Cahn model only furnishes an approximately linear kinetics behavior valid in the regime where both (3.9) and (3.25) are satisfied 13 .
The general kinetics model, by contrast, provides a regularization of the sharp-interface model with arbitrary kinetic relations given by the functions G 180 and G 90 for 180 • and 90 • domain walls, respectively. More specifically, for straight walls of all types, we have shown with traveling wave solutions that the relations between wall velocities and associated driving tractions are exactly given by G 180 and G 90 , provided that the (little restrictive) conditions ζ < 2 and ζ < 2 √ 2h 90 are satisfied. Overall, they again furnish upper bounds on the regularization length l.
In the present article, we have focused on analytically tractable features of the Allen-Cahn and general kinetics models. For that reason, we have restricted our attention to the properties and motion of straight domain walls. With this focus, we made the important clarification that with the Allen-Cahn model, the target kinetics is attained only in the limit ofl → 0, while the general kinetics model furnishes exactly the chosen kinetics, irrespective of the regularization length, provided the latter satisfies specific upper bounds. We point out that, although our discussion was specific to ferroelectric domain wall motion, the same general kinetics model bears potential for other applications involving structural transformations.
In a follow-up publication, we will report a numerical implementation of the general kinetics model and address the propagation of curved walls as well as the behavior of triple and quadruple points where multiple phases meet.
Appendix A. Experimental evidence on domain wall dynamics 180 • domain walls. Experimental and theoretical works have probed the kinetics of domain walls both in bulk single-crystals and in thin films (for an extensive review of these works, see Tagantsev et al. 2010). Here, we summarize results pertaining only to bulk BaTiO 3 . The velocity of non-ferroelastic 180 • domains walls was measured as a function of the electric field in single-crystal BaTiO 3 by Miller and Savage (1958, 1959; Savage and Miller (1960) for electric fields in the range 0.1 to 2 kV/cm and by Stadler and Zachmanidis (1963) at large fields between 2 and 450 kV/cm. These experiments were performed on platelike samples with polarization and electric field oriented along the out-of-planeẽ z -direction. As a result, the electric field e = eẽ z is uniform throughout the sample and the driving traction given by (2.9) reduces to f = 2p 0 e for an interface with negligible curvature. For electric fields in the range 0.1 to 2 kV/cm, an inverse exponential dependence of the domain wall velocity on the electric field was inferred: where V l is a characteristic velocity of the order of 100 cm/s, and the activation field e a is typically about 4 kV/cm at room temperature . The exact value of these parameters depends on temperature, sample thickness, the presence of defects, and the nature of the electrodes. These observations suggest that for 180 • domain walls in the range f ≤ 1 × 10 5 J/m 3 (cf. p 0 = 0.26 C/m 2 ), the driving traction follows a kinetic relation of inverse-exponential type, given by where f a = 2p 0 e a is an activation driving traction. Above 2 kV/cm the velocity was found to follow a power law V 180 ∝ e θ with θ = 1.4 (Stadler and Zachmanidis, 1963), which indicates for f > 1 × 10 5 J/m 3 a kinetic relation of the form with V h a characteristic velocity for large driving tractions.
90 • domain walls. While experimental data for the motion of 90 • domain walls in BaTiO 3 are lacking, observations of the kinetics of ferroelastic domain walls in other materials (Rochelle salt, galodinium molybdate, potassium dihydrogen phosphate as reviewed in Tagantsev et al. (2010), Section 8.3.3) point to a different kinetic relation specific to ferroleastic domain walls. Indeed, irrespective of the material, it appears that ferroelastic walls remain static for applied electric fields smaller than a threshold field, above which their velocity evolves linearly with the field. This suggests a threshold-type linear relation: where f 0 is a threshold driving traction and k a kinetic coefficient. At the same time, ab initio calculations by Liu et al. (2016) have shown that 90 • domain walls in PbTiO 3 exhibit the same two-regime kinetics as that observed for 180 • domain walls in BaTiO 3 and described by (A.2) and (A.3) in the regimes of small and large electric fields, respectively. While we unfortunately do not have a full set of data for the kinetics of the different kinds of domain walls in one particular material, current observations clearly indicate that the kinetics of domain walls is nonlinear and dependent on the type of domain wall (i.e., non-ferroelastic vs. ferroelastic).