Phase-Field Modeling and Peridynamics for Defect Dynamics, and an Augmented Phase-Field Model with Viscous Stresses

This work begins by applying peridynamics and phase-field modeling to predict 1-d interface motion with inertia in an elastic solid with a non-monotone stress-strain response. In classical nonlinear elasticity, it is known that subsonic interfaces require a kinetic law, in addition to momentum balance, to obtain unique solutions; in contrast, for supersonic interfaces, momentum balance alone is sufficient to provide unique solutions. This work finds that peridynamics agrees with this classical result, in that different choices of regularization parameters provide different kinetics for subsonic motion but the same kinetics for supersonic motion. In contrast, conventional phase-field models coupled to elastodynamics are unable to model, even qualitatively, the supersonic motion of interfaces. This work identifies the shortcomings in the physics of standard phase-field models to be: (1) the absence of higher-order stress to balance unphysical stress singularities, and (2) the ability of the model to access unphysical regions of the energy landscape. Based on these observations, this work proposes an augmented phase-field model to introduce the missing physics. The augmented model adds: (1) a viscous stress to the momentum balance, in addition to the dissipative phase-field evolution, to regularize singularities; and (2) an augmented driving force that models the physical mechanism that keeps the system out of unphysical regions of the energy landscape. When coupled to elastodynamics, the augmented model correctly describes both subsonic and supersonic interface motion. The augmented model has essentially the same computational expense as conventional phase-field models and requires only minor modifications of numerical methods, and is therefore proposed as a replacement to the conventional phase-field models.


Introduction
Peridynamics [Sil00] and phase-field modeling [AA12, AGDL15,Che02] are currently the leading approaches to model the evolution of microstructure and defects. An important open question is whether there are qualitative differences in the predictions of these models that cannot be resolved simply by calibration. That is, given sufficient calibration of model parameters, can both models provide similar predictions for phenomena of interest? Or, are there settings in which these models -irrespective of the sophistication of the calibration -necessarily provide qualitatively different predictions? If the predictions are different, which -if either -could reasonably be considered to be correct ? We examine this question in the context of 1-d interface motion with inertia in a material with a nonmonotone stress-strain response. Classical elasticity has shown that: (1) a kinetic law is required, in addition to momentum balance, to obtain unique solutions for subsonic motion; (2) in contrast, momentum balance alone is sufficient to provide unique solutions for supersonic motion [AK06,Tru93]. We find, in brief, that peridynamics agrees with classical elasticity while standard phase-field models do not. Following [AK91b] for strain-gradient models, we find that different choices of regularizing parameters in a given peridynamic model gives rise to different kinetics for subsonic motion but the same kinetics for supersonic motion. In contrast, we show that standard phase-field models are unable to model, even qualitatively, the supersonic motion of interfaces; supersonic motion is shown to necessarily require unbounded stresses.
Given this clear qualitative difference between peridynamics and phase-field models, the next question is which could be considered more reliable? While peridynamics and phase-field models can model complex phenomena that are beyond the reach of classical elasticity, it is reasonable to require that these more complex models recover classical elasticity when it is applicable. We therefore propose an augmentation of phase-field models that agrees with the predictions of classical elasticity.
We highlight that a major area of application of both peridynamics and phase-field modeling is to model microstructure evolution, defect motion, dynamic fracture, and so on. While these phenomena are far more complex than 1-d interface motion, we focus on the latter for several reasons. First, there is a clear benchmark solution in 1-d interface motion, unlike more complex phenomena. Second, if predictions do not agree even in simple settings, they are unlikely to agree in more complex settings. Third, the reason for the disagreement in predictions is easier to understand in a simple setting.
Standard Phase-Field Models Are Unsuitable for Problems with Inertia. We next turn to why current phase-field models are largely unsuitable for phenomena in which inertial effects -rather than energy minimization alone -play a significant role.
Consider a classical elasticity strain energy density W ( ) that is a nonconvex function of the strain . The nonconvexity implies a non-monotone stress-strain response; consequently, there are multiple strain values for a given stress ( Figure 2). This can lead to the formation of microstructure in which regions of constant strain are separated by singular sharp interfaces.
The sharp interfaces are challenging for numerical calculations, and therefore are typically regularized. In phase-field models, this regularization is done by introducing a phase-field parameter φ to keep track of the phase or energy well, and introducing gradients of φ into the energy to penalize sharp interfaces. That is, the energy density W ( ) is replaced by • W ( , φ) = w(φ) + W con ( − 0 (φ)), and regularized by adding |∇φ| 2 . We notice that while W ( ) is nonconvex in , W con ( , φ) is convex in and the nonconvexity is introduced through the nonconvex function w(φ).
The correspondence between classical elasticity and phase-field models can be seen in Figure 1. We consider homogeneous deformations to enable us to focus on the energy density. The left panel is a plot of To model the general setting without inertia, it is typical to minimize the total energy with respect to the strain field and use steepest-descent dynamics for the evolution of φ [CK14, YD10, LPM + 15, LMTS + 18, ZKL16, AA12, BH16,VTK17]. When inertia is present, it is typical to extremize the actionor equivalently, to add the inertial term to the momentum balance -and retain the steepest-descent dynamics for φ [AD17, BRLM17, BVS + 12, AGDL15, PZM + 20, GLH + 19].
Energy minimization plays a central role in relating W and • W : it ensures that the material stays near the minimizing curve φ min ( ) in the energy landscape of Figure 1. However, in problems that include inertia, energy minimization is not relevant, and the kinetic energy has an important contribution. The material can then explore parts of the energy landscape far from the minimizing curve; however, the energy landscape away from the minimizing curve has no physical connection to the original classical elasticity energy. We see in our numerical calculations that the material does explore nonphysical regions of the energy landscape -particularly at large interface velocities (Section 4.A.1). This is a central reason for phase-field models to fail in correctly modeling supersonic interfaces . Therefore, an important element of our augmented phase-field model is a physically-motivated driving force that keeps the system away from these nonphysical regions.  The Proposed Augmented Phase-Field Model. As noted above, one problem with existing phase-field models is that they allow the material to explore unphysical parts of the energy landscape when inertia is important. A second problem is that there exists a singularity for all supersonic interfaces, i.e., the strain and stress necessarily go to infinity at some point for an interface that moves supersonically. The first problem shows up in numerical solutions of initial-value problems. Specifically, as the interface velocity approaches the sonic velocity, we find that the φ interface has a different velocity and spatial location than the interface. This implies that the spatial region between the two interfaces is in the top-left or bottom-right quadrants of the energy landscape in Figure 1 (left) -the unphysical regions. The energy • W in the unphysical regions should be infinitely high to respect the original energetic formulation; this would keep the system from exploring these regions, but also cause severe practical difficulties. However, we notice that if it were possible to set • W to infinity, this would be reflected in an additional driving force contribution to the dynamical equation for φ. We therefore augment the dynamical equation for φ by the missing driving force, and the consequence is that it nudges the system downward in the top-left quadrant, and upward in the bottom-right quadrant.
The second problem is the appearance of an unphysical singularity that forces the stress and strain to go to infinity at a point for a supersonic interface in the standard phase-field formulation (Section 5.B). We also find that a regularization of the momentum equation -in addition to the usual regularization of phase-field models -resolves this singularity. While it possible to use various regularizing stresses, we choose to use a viscous stress because this is the simplest, has a clear physical interpretation, and is readily compatible with standard numerical methods, e.g. FEM with C 0 continuity [BVS + 12, AGDL15, KMB18, GLH + 19]. We note that all real materials have some level of dissipation, and even in materials in which dissipation is generally small, it can be very important in problems of shocks that are near sonic or supersonic [Daf05,AK06]. Therefore, it is not surprising that it plays an important role in phase-field models that aim to be valid when inertial effects are significant.
A Note on Strain Gradient Models. Strain gradient models are an important class of regularized models of elasticity [AK91b,Ros95,Tru93,Tur97]. They use energetic terms of the form |∇ | 2 in the energy to penalize singularly sharp interfaces. In contrast to phase-field models, they do not introduce any extra fields but have the displacement field as the sole primary field. However, the strain gradients impose additional restrictions on the continuity of the displacement field that can be challenging for standard FEM. We notice a heuristic connection between mixed FEM for problems with higher derivatives and the replacement of strain gradient models by phase-field models. In both cases, we introduce auxiliary variables that nominally relax the smoothness requirements, and then constrain the auxiliary variables to the primary variables. Phase-field models can be considered analogous to further replacing the constraint by a penalty, which can be justified by energy minimization.
In Section 2.C, we discuss the findings in the literature on using a strain gradient model to study the problem of interface motion. In summary, the strain gradient model provides predictions that agree well with classical elasticity: model parameter-dependent kinetics for subsonic motion, and parameter-independent kinetics for supersonic motion. This raises the question of why one cannot simply use strain gradient models rather than either of peridynamics and phase-field. While strain gradient models would work well for the particular problem studied here, a number of reasons make strain gradient models unsuitable for broader application. First, while strain gradient models are useful for regularizing problems that can be described by classical elasticity, it is unclear how to use them for fracture where the displacement itself is discontinuous . Second, the dependence of the nucleation and kinetics of interfaces on model parameters is extremely opaque and practically impossible to rationally specify, in contrast to the phase-field models discussed in this paper [AD15a,AD15b,AZ05]. Third, the higher derivatives that appear in the model require nonstandard or restrictive numerical methods, compared to phase-field models that can be solved using standard FEM.
Organization. Section 2 formulates the interface motion problem, and summarizes relevant results from the literature on the solution to the problem in the settings of classical elasticity and strain gradient elasticity. Sections 3 and 4 present, respectively, the peridynamic and phase-field model solutions to the interface motion problem. Section 5 discusses the reasons that existing phase-field models are unable to model supersonic interface motion. Section 6 presents and characterizes the augmented phase-field model that is able to correctly model the subsonic and supersonic behavior of interfaces. Section 7 provides further discussion.

Formulation and Classical Results
The entire paper works in 1-d; our domain is a 1-d bar denoted Ω. Where closed-form calculations are possible, we will consider Ω to correspond to the entire real line. Where numerical calculations are required, we will use a finite bar, and take care to only consider results that are not influenced by the boundaries.
The displacement of the material point at the spatial position x at time t is denoted u(x, t); the stress by σ(x, t); the strain by = ∂ x u(x, t); and the phase by φ(x, t).
The material response is formulated to have two phases, denoted phase 1 for the low-strain phase and phase 2 for the high-strain phase. These phases can coexist in certain situations, and in those situations each phase occupies distinct regions of space separated by interfaces between them. Throughout the paper, we focus on the motion of individual interfaces. Various quantities can be discontinuous across the interfaces, but, in line with the fundamental assumptions of continuum mechanics, we require that the displacement is always continuous in space and time. We denote the location of the interface by s(t) and the velocity bẏ s(t). The jump g(x = s + , t) − g(x = s − , t) across the interface of a quantity g(x, t) is denoted g .
Our convention is to have phase 1 on the left and phase 2 on the right. Therefore,ṡ > 0 corresponds to a transformation of phase 2 to phase 1, andṡ < 0 corresponds to a transformation of phase 1 to phase 2.
For use in further sections, we introduce H l (·), a regularized Heaviside / step function that transitions from 0 to 1 over a scale l as its argument transitions from negative to positive. We emphasize that l is not a lengthscale, but is used to scale the phase φ. For computations, we use the choice H l (x) = 1 2 (1 + tanh (x/l)).
We mention recent progress by P. Rosakis and coworkers [RHA20].

2.A. Material Response
The stress-strain responseσ( ) and strain energy density W ( ) = 0σ (˜ ) d˜ are plotted as a function of strain in Figure 2, and have the following expressions: (2.1) The quantities E 1 and E 2 are the elastic moduli of phases 1 and 2; 1m and 2m are the limits of existence of phases 1 and 2; and are chosen to ensure continuity ofσ( ).
We assume for simplicity that the density ρ is constant, and that E 1 > E 2 . We define the sonic velocities c 1 = E 1 /ρ and c 2 = E 2 /ρ, corresponding to small-amplitude linearized waves in the phases; notice that c 1 > c 2 . Subsonic interfaces have −c 2 <ṡ < c 2 , and supersonic interfaces have c 2 <ṡ < c 1 . Interfaces withṡ > c 1 orṡ < −c 2 are not permitted by momentum balance and thermodynamics.

2.B. Results from Classical Elasticity
The behavior of interfaces in the material described in Section 2.A has been studied, among various other topics, in the body of work of Abeyaratne and Knowles [AK06]. We briefly summarize the relevant details of their results here.
First, we consider the quasistatic setting. The balance of momentum leads to a PDE where the fields are smooth and a jump condition at the interface: This implies that the stress σ(x, t) is constant in the bar, i.e. σ(x, t) = σ 0 (t). Considering load control, i.e. σ 0 (t) specified, it is clear from Figure 2 that when σ 0 (t) ∈ [2, 5], there are an infinity of solutions that satisfy equilibrium. Specifically, any displacement field that everywhere has derivative ∂ x u ∈ {σ 0 /E 1 , σ 0 /E 2 } will satisfy momentum balance. In the context of displacement control, the situation is similar when the average strain in the bar, u(L, t) − u(0, t) L , is in the range [ 2 5 , 5]. Each discontinuity in ∂ x u corresponds to an interface across which the strain jumps. Even in the simplified setting of a single interface, the nonuniqueness persists; for instance, given σ 0 (t), the solution has the form : where s(t) can be arbitrary. The nonuniqueness in the quasistatic setting is typically resolved using energy minimization as a selection mechanism [Eri75]. For instance, for solutions of the form in (2.4), we find s by minimizing the potential energy over s. Next, we consider the dynamic setting. The balance of momentum reads: The jump condition (2.6) has an insightful graphical interpretation ( Figure 3). Writing it asṡ 2 = 1 ρ σ ∂ x u , we see that the velocity of the interface is related to the slope of the chord on the stress-strain curve connecting the stress-and strain-states ( − , σ − ) and ( + , σ + ) on either side of the interface. Further, we recall that the sonic velocity in each phase is similarly related to the slope of the stress-strain curve of the corresponding phase, i.e. c 2 1 = E 1 ρ and c 2 2 = E 2 ρ . Indeed, sonic waves are governed by precisely the jump condition (2.6), but with both end states on the same branch. The graphical interpretation tells us that any interface that connects the two branches is represented by a chord whose slope must be less than the slope of the phase 1, implyingṡ < c 1 always. However, we notice thatṡ can be smaller or larger than c 2 , and we therefore consider two regimes. First, the subsonic regime, wherein 1 < M 2 :=ṡ c 2 ≤ 1 and the interface is subsonic with respect to phase 2; and, second, the supersonic regime, wherein M 2 > 1 and the interface is supersonic with respect to phase 2. The subsonic regime −1 < M 2 ≤ 1 inherits the non-uniqueness that was evident in the quasistatic setting. In particular, initial-boundary-value problems with subsonic interfaces do not have unique solutions; further, since the problem is dynamic, energy minimization is not an appropriate selection principle. The non-uniqueness can be related to a lack of information of the kinetics of the interface. That is, we require a kinetic relation that relatesṡ to ( − , σ − ) and ( + , σ + ) to obtain unique solutions .
The kinetic relationṡ =v(f ) relates the velocity of the interface to the driving force f acting on the We could also have phase 1 on the right and phase 2 on the left. The general problem requires also nucleation criteria, but we focus on the behavior of a single already-nucleated interface throughout this paper. The slope of a chord is proportional to the square of the interface velocity, and the slope of a stress-strain branch is proportional to the square of the sonic speed for that branch. The subsonic chords have slope such that |ṡ| < c 2 , and the supersonic chord has c 2 <ṡ < c 1 . Notice that it is not possible to construct a chord that hasṡ > c 1 .
interface, and was introduced by [AK90, HL85,Tru82]. The driving force is given by the expression: It contains information about the state of the material on both sides of the interface, and is precisely the work conjugate ofṡ. The supersonic regime M 2 > 1, however, has unique solutions without a kinetic relation [AK91a, Tru93,TV10]. That is, (2.5) and (2.6) have a unique solution when the interface is supersonic with respect to phase 2. No additional kinetic relation is required, and using such a kinetic relation will generally over-constrain the problem such that there are no admissible solutions.
We mention that the chord construction -i.e., momentum balance -does not rule out interfaces with M 2 < −1, but these would have negative dissipation and hence are ruled out by thermodynamics.

2.C. Results from Strain Gradient (Viscosity-Capillarity) Models
Strain gradient models, also called viscosity-capillarity models, regularize the sharp interfaces of classical elasticity by adding a strain gradient term to account for the surface energy, and a viscous term to account for the dissipation associated with defect motion [FM06, AK91a, AK06, Ros95, Tru93, Tur97]. We use σ(x, t) =σ(∂ x u) + ρν∂ xt u − ρλ∂ xxx u to find the field equation of momentum balance: Here, ν is the coefficient of dissipation and λ is the coefficient of surface energy. The solution is sufficiently smooth due to the higher derivatives, and therefore the jump condition, (2.6), is not required. This model typically has unique solutions given appropriate initial and boundary conditions [AK91a].
An important finding in [AK91a,Tru93] is that this regularization preserves the key distinction between the subsonic and supersonic regimes. That is, the kinetics of subsonic interfaces depends sensitively on the choice of ν and λ, whereas the kinetics of supersonic interfaces is relatively insensitive to this choice. One can therefore think of ν and λ as inducing a kinetic relation when it is required for uniqueness, and providing merely a minor regularizing effect when the kinetic relation is not required for uniqueness.

2.C.1. Numerical Computation of Kinetics Using Traveling Waves
As mentioned above, the addition of strain gradient and dissipation terms effectively induces a kinetic relation. While the kinetic relations induced by (2.8) have been computed in closed-form in [AK91a], we nonetheless compute these numerically here as a means to both describe as well as verify our numerical scheme, that we will apply later on to other models studied in this paper.
We begin by assuming a traveling-wave form for the solution: u(x, t) = U (x −ṡt). Using this in (2.8), we find: wherex := x −ṡt is the traveling coordinate, and primes represent differentiation with respect tox.
We can immediately integrate (2.9) once to get: where we have nondimensionalized; introduced the strain E(x) := U (x); and C is the undetermined constant of integration.
The displacement U is not unique due to rigid-body translations, and we have no boundary conditions to fix this; therefore, we will solve (2.10) directly for the strain for various given values of M 2 .
We consider possible boundary conditions for (2.10). To do this, consider the limits E(x → ±∞), denoting these by E ±∞ . Using that the derivatives of E vanish far from the interface [AK91a], we find the equations: Subtracting the equations above, we find This equation is precisely the strain-gradient version of (2.6), and does not provide any new information in the regularized context. However, we notice that if we use E +∞ as given data, along with M 2 given, we can use (2.12) to find E −∞ , and vice-versa. That is, if we pick either of the far-field strains as data along with given M 2 , we can solve for the other far-field strain, and consequently solve also for the driving force: Thus, if we specify M 2 and either of E ±∞ , we can find the driving force and consequently the kinetic relation corresponding to M 2 without using (2.9) at all; in particular, the dissipative and surface energy contributions play no role in determining the kinetic relation. We conclude that specifying either of E ±∞ overconstrains the problem at a given M 2 .
Therefore, we do not specify either of E ±∞ and treat both as unknown quantities. Following [AD15a, AD15b, DB06], we solve (2.10) by treating the entire function E(x) as well as C as unknown and solve for them using a least-squares approach. Defining the residue: we solve by minimizing R = +∞ −∞ r(x) 2 dx over C and functions E.
The computational discretization uses a finite segment of the bar of length L in the translating framex. We use finite differences and divide the domain into N equal intervals, each of length ∆x = L/N , with nodesx i , and define E i := E(x i ). We approximate the residue by the quadrature: (2.15) and use left-or right-rather than centered-differences at the edge of the domain. We minimize R in (2.15) over the finite set of variables {E i , i = 2 . . . N } and C to find the traveling wave profile, using a standard monolithic solver. To test that we are not stuck at a local minimum, we check that the value of R is close to 0 after the minimization is complete.
Using this numerical procedure, for various choices of M 2 , we compute the strain profile and use this to infer the driving force through (2.13); this procedure is repeated for several choices of ω where ω = 2 √ λ ν . The corresponding kinetic curves and some representative strain profiles are shown in Figures 4 and 5. These match well with the closed-form expressions obtained in [AK91a] and provides us with confidence in the numerical scheme.

Interface Motion in Peridynamics
Following the model postulated by [Sil00], the 1-d peridynamic equation of motion is: where f(δu, δx) is the bond force between two volume elements with separation in the reference δx :=x − x and relative displacement δu It is useful to have a dissipative mechanism to account for the dissipation in interface motion. Instead of adding terms containing the strain rate as in continuum mechanics -which would nullify the goal of peridynamics of eliminating spatial derivatives -we follow [DB06] and add a dissipative contribution to f: where l 0 is the nonlocal length scale and ν b is the dimensionless coefficient of viscous bond-level damping. The argument forσ (the classical stress-response function from (2.1)) is the bond strain δu/δx, rather than the classical strain . This choice ensures that the stress-strain response for homogeneous deformations, computed sufficiently far from the boundaries, is identical to that chosen in strain-gradient and classical elasticity [WA05, BD18, DB06, TR14]. We noticed in strain gradient models that changing the parameters ν and λ induced different kinetic relations for subsonic interfaces. In peridynamics, we will analogously change the parameters ν b and l 0 to induce different kinetic relations. Regardless of the values chosen for ν b and l 0 , the form of the expression for the bond force in (3.2) gives us the stress-strain responseσ in the setting of homogeneous deformations. That is, changing ν b and l 0 leaves the homogeneous stress-strain response unchanged.

3.A. Numerical Computation of Kinetics Using Traveling Waves
Similar to the strain gradient approach, we seek a solution in the form of a traveling wave Substituting (3.3) into (3.1), we have: We highlight that we use a finite domain [0, L] for the numerical calculations, and set L l 0 and ensure that the interface is far the boundaries.
The residue is defined as: Here, we diverge from the method used for the strain gradient model due to the different nature of the boundary conditions in peridynamics. Specifically, boundary conditions in peridynamics are applied over a finite layer [Sil00]. In the context of a traveling wave where we do not have boundary conditions to apply, the procedure to deal with the boundaries follows [DB06,Day17]. First, we decompose our domain [0, L], denoted Ω, into boundary regions on the left and right, denoted Ω − and Ω + respectively, and an interior region, denoted I. The boundary regions have size l 0 such that there are negligible interactions between points in I and points beyond the boundary regions. We now define our solution as the minimization of We notice that R is a functional of U (x) over Ω and not I, despite the domain of integration in the definition of R; due to the nonlocal interactions, R involves a double integration.
We now discretize Ω into N equal intervals, each of length ∆x = L/N , with nodesx i , and define U i = U (x i ). Then, we approximate the integrations using the quadrature: and use standard central differences for the derivatives. We use ∆x = l 0 10 as numerical experiments show that the results are essentially converged with this level of discretization. We set the average of U over the domain to be 0 to fix the rigid translation.
We highlight that we minimize the error only over the interior I, but the variables over which we minimize are all the nodal values. We refer to [AD15a,DB06,Day17] for more discussion of this approach.
We use this numerical procedure to compute the traveling wave profile for various values of M 2 . Given the displacement profile, we can compute the driving force using (2.13), thus giving a kinetic relation that relates the interface velocity to the driving force. This procedure is repeated for several choices of ν b and l 0 . The corresponding kinetic curves and some representative displacement derivative profiles are shown in Figures 6 and 7.

Interface Motion in Existing Phase-Field Models
Two existing phase-field models will be studied in this paper. While there are several differences between these models, they share 2 key features: (1) both have gradient regularizations of the phase-field parameter, and not in the momentum balance; (2) neither has a mechanism to prevent the system from accessing unphysical regions of the energy landscape. These features are common to all phase-field formulations that we are aware of.
The first model -or closely-related variations -is completely standard, and has been used in the overwhelming majority of prior works, e.g. [Che02, YD10, AA12, AD17, BH16, PMB + 20, LPM + 15, LMTS + 18, ZKL16, AHKB20, CMB20]; we will refer to it as the "standard phase-field model" for short. In standard phase-field models, the energy landscape is formulated based on equilibrium principles of energy minimization, and has the overall structure of the form shown in Figure 1. While the model gives rise to unique evolution of microstructure -i.e., nucleation and kinetics of interfaces is contained in the model -the relation between the model parameters and the nucleation and kinetics of interfaces is completely opaque. Further, nucleation and kinetics are also coupled in the sense that changing the model parameters typically changes both the nucleation and kinetic behavior simultaneously.
The second phase-field model that we will study was motivated by the goal of overcoming the shortcomings discussed just above; we will refer to it as the "dynamic phase-field model" for short. Specifically, the dynamic phase-field model is formulated to provide a transparent separation between energetics, kinetics, and nucleation. There are distinct model parameters that correspond to each of these, i.e., we can specify the equilibrium response, the kinetics of interfaces, and the nucleation of interfaces independently. This model was proposed recently in [AD15a]. For both models, we perform initial-value numerical calculations to examine the propagation of interfaces. The results of the initial-value problems, and the analysis in Section 5, show that supersonic interfaces cannot be predicted by either of these models. Therefore, we do not present our attempts to solve traveling wave problems.

4.A. Standard Phase-Field Model
The standard phase-field model is formulated primarily on the basis of appropriately constructing the energy landscape to obtain the correct equilibrium response. The evolution is obtained by assuming that φ follows a steepest descent dynamics, coupled to static or dynamic momentum balance.
We use large computational domains and take care not to consider results that have been affected by the boundaries. Therefore, we refer to these problems as initial-value problems. Strictly speaking, the computational domain is finite and one could refer to them as initial-boundary-value problems. However, they aim to mimic a problem on an unbounded domain where the boundaries play no role.
The form of the energy of a standard phase-field model coupled with piecewise linear elasticity is: where w(φ) is a nonconvex energy that favors the formation of interfaces, while the gradient term 1 2 α|∂ x φ| 2 regularizes them. α is the gradient energy coefficient, and has broadly the same physical interpretation as λ in Section 2.C.
We choose to have phase 1 indicated by φ = 0 and phase 2 by φ = 1. Therefore, we choose for w(φ) the expression: The term φ 2 (1 − φ) 2 brings in the nonconvexity of the energy landscape, with minima at φ = 0 and φ = 1. Θ is a large constant, chosen to be 10 3 for this work. The term φ∆Ψ accounts for the fact that these phases have energy minima at different heights, with the difference quantified by ∆Ψ that was introduced in (2.1). We discuss the relation between the phase-field energy and the classical elastic energy in more detail in Section A. For the elastic response, we set 0 (φ) ≡ 0, and E(φ) = (1 − H l (φ − 0.5))E 1 + H l (φ − 0.5)E 2 , which transitions smoothly from E(φ 0.5) = E 1 to E(φ 0.5) = E 2 . This mimics the transition between the branches described in (2.1): we have a linear response with modulus E 1 when φ ≈ 0 and a linear response with modulus E 2 when φ ≈ 1. We notice, however, that the elastic energy 1 2 E(φ)(∂ x u) 2 is defined for any combination of φ and ∂ x u. Consequently, it is possible that the system has a strain value corresponding to phase 1 while φ ≈ 1; see Figure 8.
From (4.1), we find the stress and driving force: The evolution equations are given by: where we have used the standard steepest descent assumption for φ, and the constant µ is the mobility.

4.A.1. Kinetics of Interfaces from Initial-Value Problems
The kinetics of interfaces are studied through numerical solutions of initial-value problems. We use standard explicit time-stepping to solve linear momentum balance along with the evolution equation for φ for the kinetic response. Our domain is a long finite bar with an interface at the center of the bar. Our initial conditions correspond to a displacement / strain field that is not at equilibrium, and the interface has a non-zero driving force across it. Therefore, the interface will move in the direction of the driving force. This is the analog of the classical Riemann problem [DB06]. Figure 9 shows representative results for interfaces that are well below the sonic speed. These results show the system behaving as we would expect: faster acoustic waves going in both directions from the initial interface, with signatures only in the strain profile; and the slower interface which has a signature in both the strain and φ profiles.
An example with a large driving force is shown in Figure 10. An important feature that we notice is that interfaces in strain and φ are not at the same location nor do they move with the same velocity. It is therefore not possible to usefully define an interface velocity. The strain interface has barely moved and is subsonic, while the φ-interface is moving faster than both c 2 as well c 1 . Most importantly, the material between the strain interface and the φ interface is in phase 1 defined through the location of the strain interface but in phase 2 defined through the location of the φ interface. It is therefore in an unphysical part of the energy landscape.

4.B. Dynamic Phase-field Formulation
The dynamic phase-field model aims to transparently separate energetics, kinetics, and nucleation. In this section, we present the model equations and only those aspects that are directly relevant; we refer to [AD15a,AD15b] for the details of the formulation and the characterization. The kinetics in this model is similar to that proposed earlier by [AZ05] though they have not focused on the energetics or nucleation.
We construct the phase-field energy as follows: This corresponds to an expansion of the energy about some strain A 0 that need not correspond to the stressfree strain. As a consequence of the piecewise-linearity of the energy, the final expressions turn out to be independent of the choice of 1 0 , 2 0 ; using (2.1) with 1 0 < 1m and 2 0 > 2m , we have: We collect expressions and simplify (4.7) to get: Figure 9. Initial-value problems for the standard phase-field model with interfaces moving at subsonic velocities. (a) interface moving to the right, into the soft phase; b) interface moving to the left, into the stiff phase. In both cases, we notice acoustic waves propagating in both directions. The interface can be identified by using that it appears in the plots of both φ and ∂ x u, whereas the acoustic waves have a signature only in ∂ x u.
The key difference is that w(φ) in (4.1) for the standard energetic phase-field model is replaced by H l (φ − 0.5)∆Ψ above. We note that E(φ) = (1 − H l (φ − 0.5))E 1 + H l (φ − 0.5)E 2 above, as in the standard phase-field model. The stress response function and driving force in this phase-field model are: and the evolution equations are given by: wherev φ n is the velocity of the interface and controls the interface kinetics; and G controls the interface nucleation. In general, bothv φ n and G can be functions of any quantity such as f, σ, and their rates; thermodynamics imposes some weak conditions on their dependence on f .
We assume the simplest linear kinetics, i.e.v φ n = κf with no dependence on other quantities, and G ≡ 0. Figure 11 plots out the local energy density and the stress. A key feature of this energy is that it is flat in (a) (b) Figure 10. Initial-value problem for the standard phase-field model close to the sonic velocity. The strain interface has barely moved and is subsonic, while the φ-interface is moving faster than both c 2 as well c 1 ! While waves above c 1 are not permissible by momentum balance, this only constrains the strain and not the evolution of φ. Further, the material between the strain interface and the φ interface is in phase 1 defined through the location of the strain interface but in phase 2 defined through the location of the φ interface. It is therefore in an unphysical part of the energy landscape. the φ-direction away from the transition at φ ≈ 0.5. This feature is critical in decoupling nucleation from the structure of the energy landscape. In particular, this energy landscape simply does not permit nucleation of new phases, regardless of the level of stress / strain. The only available mechanism for the nucleation of a new phase is through the term denoted G in (4.14). This effectively decouples the nucleation of new phases from the equilibrium energetic response, and provides a simple and transparent mechanism to specify the precise conditions for nucleation through the functional dependence of G on f, σ, ∂ x u, their rates, and so on.

Unphysical Features of Existing Phase-Field Models
We discuss here the 2 main reasons that existing phase-field models are unable to properly handle situations with inertia. The first reason is related to energetics, namely that the system explores unphysical regions of the energy landscape; further, we highlight that the energy landscape -even in the energy minimizing setting -can have unexpected behavior at large strains. The second reason is that the momentum balance, as typically formulated, leads to a strain singularity for supersonic interfaces.

5.A. Energy Landscape at Large Strains
In the 1-d piecewise-linear phase-field energy considered in this paper, the elastic energy has the form 1 2 E 1 2 and 1 2 E 2 2 + ∆Ψ in phases 1 and 2 respectively. If E 1 > E 2 , we notice that phase 2 always has lower energy when 2 is large enough, irrespective of the value (or even sign) of ∆Ψ . This is simply because the quadratic growth eventually wins. If is positive and large, then phase 2 is lower energy; however, even if is negative with large magnitude -considering, e.g., the case of antiplane deformation -phase 2 again eventually has lower energy. This unphysical behavior is shown in Figure 12. This conclusion can be readily generalized to higher dimensions and more general energies. In general, it is a consequence of defining the energy density as a function of φ and , and this higher-dimensional energy landscape has several unphysical regions.
In the relatively simple setting considered here, this pathology in the energy landscape effectively introduces an extra, unphysical, transformation that could potentially activate to give unphysical results. For instance, when we apply large driving forces to attempt to drive the interface at a high velocity, we find that the transformation occurs from the soft phase 2 at high positive strain to the soft phase 2 at high negative strain! We then observe the unphysical situation of a large strain jump across the interface, but φ does not evolve at all because the material is in the same phase.

5.B. Inability to Model Supersonic Interfaces
As observed in Section 4, neither of the existing phase-field models could model supersonic interfaces. We see below that the momentum equation, as used in those models, forbids it.
Consider a steady supersonic interface moving at a given velocityṡ. Using the traveling-wave form u(x, t) = U (x −ṡt) in the momentum equation gives: where we have integrated once to go from the second to the third step above; the constant therefore is independent ofx = x −ṡt, but can be a function ofṡ and material parameters. Since the interface is supersonic, we have that c 2 <ṡ < c 1 , implying that E 2 < ρṡ 2 < E 1 . Further, we have that E(φ) varies smoothly between E 1 and E 2 as φ transitions from phase 1 to phase 2. Furthermore, since φ(x) is smooth, there will be some pointx * in the traveling wave coordinates at which ρṡ 2 = E(φ). Atx * , we consequently have that ρṡ 2 − E(φ) = 0, and combining this with ρṡ 2 − E(φ) U = const., we have two possibilities: (1) either U → ∞ atx * ; or (2) the constant must be 0, and consequently U must be zero everywhere exceptx * . Neither of these possibilities is acceptable for the interfaces considered here, and therefore we conclude that these phase-field models cannot model supersonic interfaces.
Notice that our argument depends on the continuous variation of E(φ) with respect to its argument. If this variation was discontinuous, then we would have to track the moving surface across which E(φ) was discontinuous, and that would nullify the most important advantage of phase-field modeling, namely that we do not have to track singularities.

Augmented Dynamic Phase-Field Model for Microstructure Evolution with Inertia
To address the issues identified in Section 5, we propose an augmented phase-field model that has these 2 extra terms: 1. a local dynamical term, corresponding to G in (4.14), that moves the system -along the φ direction only, to avoid disrupting momentum balance -away from unphysical regions in the energy landscape. This term corresponds to accounting for the missing physics of a driving force that would drive the evolution away from high-/infinite-energy forbidden regions. 2. a viscous dissipative stress that regularizes the singularities identified in Section 5.B, and accounts for the missing physics of dissipative mechanisms that are always active and particularly important at defects and singularities. Figure 13(b) shows the energy landscape and the unphysical regions that we aim to avoid. To avoid introducing new terms in the momentum balance that could lead to spurious unphysical artifacts for shock and acoustic waves, we only modify the evolution of φ by introducing driving forces in (4.14) through the term denoted by G. In both the standard phase-field model and the dynamic phase-field model, this involves adding a local contribution to the evolution equation for φ, i.e., ∂ t φ = . . . + G, where G can be a function of any quantity in the system. G is 0 when we are in physical regions of the energy landscape; if we enter an unphysical region -detected by examining the values of φ and -we set G to be nonzero. That is, in short, G is a function of φ and . In more general settings, G could be a function of the rates, specific components of the stress or strain, and so on. Returning to our specific setting, if the system enter the northwest quadrant of Figure 13(b), G is activated and takes on a negative value to push the system downwards into the southwest quadrant. Similarly, if the system enter the southeast quadrant, G is activated and takes on a positive value to push the system upwards to the northeast quadrant. The specific expression that we use for G is:

6.A. Augmented Driving Force to Drive Evolution from Forbidden Regions of the Energy Landscape
, and G 0 is a magnitude that we set as large as possible while not losing numerical stability. The expression above is plotted out in Figure 13(b).
We note that if we interpret G as a nucleation mechanism in the sense of [AD15a], our choice of G would appear to violate thermodynamic requirements. However, this apparent violation of thermodynamics is itself an artifact of defining the energy without accounting for the unphysical regions. A correct definition of the energy landscape would set the energy to ∞ in the unphysical regions, but this would be difficult to use for practical computation. However, in such a landscape where the energy is infinite in the unphysical regions, our choice of G is acceptable to thermodynamics as well as bringing in the missing physics of a driving force to keep the system away from the unphysical regions.

6.A.1. Quasistatic Characterization
We characterize the behavior of the phase-field model with the strain/phase constraint -but without inertia -to illustrate its effect. Specifically, we solve numerically the equations: with load-control and a time-varying applied load. Because we are in 1-d, this corresponds to simply prescribing σ(t); we set this to be a piecewise linear function of time to model loading and unloading. We start at zero stress and zero strain with φ = 0 (phase 1) in the entire specimen; load it until the entire specimen has transformed to phase 2; and then unload back to zero stress and strain. The results are shown in Figures 14 and 15 for the standard and dynamic phase-field models respectively. For the standard phase-field model, without the strain/phase constraint G, we find that the nucleation of the forward transformation is controlled by the height of the energy barrier, while reverse transformation does not occur at all. On the other hand, with the term G included, nucleation of both forward and reverse transformations occur at precisely the values that we set in the definition of G in (6.1). For the dynamic phase-field model, with the term G included, we find similar desirable behavior. Without G, there is no nucleation in either direction, as we desire from the dynamic phase-field model, showing that nucleation is uncoupled and controlled only by the term G.

6.B. Viscous Dissipative Stresses in the Balance of Momentum
We recall the simple argument in Section 5.B that showed that supersonic interfaces could not be modeled by the existing phase-field models. In brief, momentum balance in the form ρ∂ tt u = ∂ x σ, in combination with a traveling wave ansatz, gave ṡ 2 ρ − E(φ) U = const.. Since ṡ 2 ρ − E(φ) is 0 at some location for a supersonic wave, it follows that U blows up at that point, or is zero essentially everywhere.
If we regularize this equation by using a stress response that includes higher derivatives such as ∂ xxt u or ∂ xxxx u, or indeed any number of other possibilities, we find that our simple argument no longer holds. We choose to add a linear dissipation of the form ∂ xxt u, because it is simple to use and rooted in the physics. With G, we are able to precisely and independently prescribe the critical stress for forward and reverse transformations. Without G, the system transforms from phase 1 to phase 2 at a critical stress that depends on the energy barrier, that is controlled by Θ. The reverse transformation does not occur at all, even when we go to negative strain; at large negative strains, phase 2 is again stable per the discussion in Section 5.  The equation of momentum balance will then have the form: We will see in the numerical characterization of the augmented model that this regularization, in combination with the augmented driving force, is sufficient to predict supersonic interfaces.

6.C.1. Augmented Standard Phase-Field Model
Using the model from Section 4.A, we have the following expressions for the free energy, stress, and driving force: With the viscous stress and the augmented driving force, the evolution equations are given by: Numerical solutions of initial-value problems with supersonic interfaces are shown in Figures 16, 17, and 18, for the case with only the augmented driving force, the case with only viscous stress, and the case with both mechanisms respectively. We present the results with only the augmented driving force and only the viscous stress to show that both of these mechanisms are essential.  Figure 16. Initial-value problems for the standard phase-field model with supersonic interfaces, using only the augmented driving force G. (a) ∂ x u, (b) φ. We see that φ has undesirable oscillations, but the interfaces in strain and φ move together, showing that the system is not in unphysical regions of the energy landscape.
We next perform traveling wave calculations to find the kinetic relations for various choices of model parameters. The overall numerical approach is similar to that described in Section 2.C.1. The primary difference is that we have 2 simultaneous equations to solve in (6.6). Our solution strategy is to compute the residual for each equation, square both residuals, and then add the squared residuals and integrate over the domain to obtain a single functional that we can minimize. Figure 19 shows that the kinetics is sensitive to the model parameters for subsonic interfaces, but is not for supersonic interfaces; Figure 20 shows some representative strain profiles.  Figure 17. Initial-value problems for the standard phase-field model with supersonic interfaces, using only viscous stresses. (a) ∂ x u, (b) φ. We see, as before, that the strain interface is subsonic while the φ interface is well above supersonic. Therefore, the system explores unphysical regions of the energy landscape.
(a) (b) Figure 18. Initial-value problems for the standard phase-field model with supersonic interfaces, using both the augmented driving force and viscous stresses. (a) ∂ x u, (b) φ. The evolution is precisely as we desire, in that the strain and φ interfaces move together supersonically with no undesirable oscillations. Notice that the strain interfaces are smeared out due to the additional dissipative regularization. Note that the "blip" in ∂ x u at t = 1 is a transient that has not yet stabilized into a steadily-moving interface.

6.C.2. Augmented Dynamic phase-field model
Using the model from Section 4.B, we have the following expressions for the free energy, stress, and driving force: With both the viscous stress and the augmented driving force included, the evolution equations are given by: ρ∂ tt u = ∂ x σ + νρ∂ xxt u (6.9) ∂ t φ = |∂ x φ|v φ n + G(∂ x u, φ) (6.10) and we use linear kineticsv φ n = κf as in Section 4.B. We perform numerical computations of both initialvalue problems and traveling wave problems. The results of the initial-value problems are qualitatively identical to those obtained with the augmented standard phase-field model reported in Section 6.C.1, and we do not present the details. In summary, neither viscous stresses nor the augmented driving force by themselves lead to good results -and the bad results are qualitatively similar to those reported in Section 6.C.1 -while the use of both mechanisms together provide the desired results.
The traveling wave calculations of interface kinetics also show the desired behavior. Namely, the kinetic relation for subsonic interfaces is sensitive to model parameters, while the kinetics of supersonic interfaces is not ( Figure 21); representative strain profiles are shown in Figure 22.

Discussion
Dynamic interfaces in a non-monotone stress-strain material are predicted by classical elastodynamics to have two regimes: subsonic where the evolution is nonunique, and supersonic where the evolution is unique . We show that peridynamics preserves this key feature of interfaces ( Figure 6), but existing phase-field model do not; in fact, supersonic interfaces are not admitted by existing phase-field models. We propose an augmentation of phase-field models, using a viscous stress in the momentum balance and a local dynamical term that keeps the system out of unphysical regions of the energy landscape. We demonstrate that the augmented phase-field models recover the feature that subsonic interfaces are sensitive to model parameters while supersonic interfaces are not (Figures 19 and 21).
A contribution of this paper is in providing a critical qualitative test that distinguishes between the predictions of peridynamics and phase-field models. Prior work has shown that both peridynamics and phasefield models trivially recover the homogeneous deformation limit, and rigorous work for both peridynamics An interesting analogy to this appears in the work of Eshelby on liquid crystals, where he found that the configurational force on a disclination is identical to the real force [Esh80,CF02,SH20]. The coexistence of regimes with uniqueness and nonuniqueness within a given problem have been explored also in the context of soft materials, e.g. [Coh19,Kno02]. and phase-field models show that they recover the energy-minimizing Griffith theory of brittle fracture without inertia [Lip14,AT90]. Recent numerical works compare against experiment [MBB20, DBWW20]; while valuable and complementary, these leave open the question if the failure to reproduce an experiment is a calibration issue or a difference in the fundamental structure of the models. This paper shows that there are fundamental differences in the dynamic setting that cannot be bridged by calibration of model parameters.
A key shortcoming of existing phase-field models that is discussed in this paper is that the energy landscape is expanded by introducing the phase-field φ, and this expanded energy landscape has numerous unphysical regions. For equilibrium problems, the system is governed by energy minimization, and hence it does not explore these unphysical regions. However, in dynamic problems with inertia, the potential energy is balanced against the kinetic energy. Therefore, energy-minimization formulations of the energy are found to be inadequate to avoid the unphysical regions, and we require additional physics to prevent the system from exploring the unphysical regions. In short, if new equations and variables are introduced, we require additional physics to ensure that they behave appropriately in all regimes of application. Our proposed augmentation of phase-field models includes this additional physics in the form of an additional driving force.
We have examined in detail two types of phase-field models, but we expect similar results from other phase-field models, e.g. [KMB18], because the key features that lead to these findings are similar in those other models. In turn, we expect that the augmentation proposed in this paper will prove useful in augmenting also those other models. It will similarly be interesting to examine if other regularized models for interfaces, e.g. [Cla17,Cla19,Day17], show similar results. Related to this, an important next step is to test the augmented model in realistic higher-dimensional problems, such as dynamic fracture. For instance, a key test is to perform numerical calculations of dynamic fracture with an augmented phase-field model and compare it quantitatively to the predictions from an existing phase-field model. In addition to the augmentation directly affecting the crack growth dynamics, it will also affect the kinetics of nonlinear waves that govern the transport of elastic energy around growing defects [Mar06, Gao97, CM15, KAILP17, FRMV17].

Software Availability
A version of the code developed for this work is available at https://github.com/janelchua/Phase-field-and-Peridynamics.git

A. Correspondence between Strain Energy Densities of Classical Elasticity and Phase-Field Models
The Maxwell stress is an important physical quantity that characterizes phase transformations modeled by nonconvex energies. It is the value of the mechanical stress at which the driving force on the interface (2.7) is 0 in the quasistatic setting. Because the driving force for interface motion is 0, the (infinitesimal) motion of the interface in either direction does not change the energy of the system. Therefore, equivalently, it is the value of the mechanical stress at which there is no energetic preference for either phase [AK06]. It is essential that the energetics of all the phase-field models proposed here give rise to the same Maxwell stress as predicted by classical elasticity. We show here that our construction of these models satisfies this requirement. First, we compute the Maxwell stress for the classical elasticity model in (2.1). Let the driving force f vanish at the Maxwell stress σ M . From (2.7), we get the condition: Since we have assumed quasistatics, we have that σ M = E 1 − = E 2 + , where ± are the uniform strains on either side of the interface. Further, using (2.1), we have that W ( + ) = 1 2 E 2 +2 + ∆Ψ and W ( − ) = 1 2 E 1 −2 . Substituting these in (A1), we find the expression for the Maxwell stress: We now consider the energetics of the standard phase-field model discussed in Section 4. The local part of the energy density is w(φ) + 1 2 E(φ) 2 . Decompose w(φ) =w(φ) + φ∆Ψ . The nonconvex energyw(φ) is given a specific form in that section, but for our discussion we only need that it contains two minima (φ = 0, 1) and that these minima are at the same height. We notice that the energy difference between the the two phases at zero stress is ∆Ψ .
Next we consider the Maxwell stress predicted by the standard phase-field model from (4.1), (4.2). Using the interpretation of the Maxwell stress as the stress at which the potential energy difference between the phases is zero gives: where we have used σ M = E 1 + = E 2 − , and ± are the uniform far-field strains.
We highlight an important approximation in our calculation above. Under stress, the energy minima are not precisely at φ = 0, 1, and hence the difference in energy between the phases is not precisely ∆Ψ . This calculation is exact only in the limit that Θ is large.
We next examine the energy of the dynamic phase-field model from (4.10). The local part of the energy is H l (φ − 0.5)∆Ψ + 1 2 E(φ) 2 . The difference between the energies at zero strain is ∆Ψ . Therefore, the value of the Maxwell stress in this model is also identical to the values in the standard energetic phase-field and classical elasticity. Further, we notice that this result is exact because the difference between the energies is independent of stress, i.e., H l (φ − 0.5) goes to 0 and 1 regardless of the stress.