Multiphase field modeling of grain boundary migration mediated by emergent disconnections

Knowledge about grain boundary migration is a prerequisite for understanding and ultimately modulating the properties of polycrystalline materials. Evidence from experiments and molecular dynamics (MD) simulations suggests that the formation and motion of disconnections is a mechanism for grain boundary migration. Here, grain boundary migration is modeled using a multiphase field model based on the principle of minimum dissipation potential with nonconvex boundary energy, along with a stochastic model for thermal nucleation of disconnection pairs. In this model, disconnections arise spontaneously in the presence of an elastic driving force, and that their motion mediates boundary migration. The effect is due to the fact that the formation of the disconnections pairs results in a stress concentration, causing the elastic driving force to exceed the threshold value and driving the propagation of the disconnection along the interface. The model is applied to study the propagation/annihilation of single disconnection pairs, the relaxation of a perturbed interface, and shear coupling at various temperatures. The results are consistent with the current understanding of disconnections, and capture the effect of thermal softening.


I. INTRODUCTION
Microstructural evolution of structural materials affects properties of materials ranging from strength and ductility [1][2][3], to radiation resistance [4][5][6], and processes ranging from grain growth and recrystallization to severe plastic deformation [7,8]. If and when evolution does occur, it must do so in a controlled and predictable manner in order to guarantee reliability during the material's lifetime. It is often desirable to prevent certain behaviors such as significant coarsening (e.g. in the case of nanocrystalline materials [9]), refinement (e.g. in the case of transformer steels [10]), or general evolution (e.g. in GB engineered materials [11,12]). These behaviors can occur in response to a variety of thermal, mechanical, chemical, or radiative loadings. Mechanical coarsening, in particular, has attracted substantial attention as a basic mechanism that couples microstructure evolution and mechanical response [3,13,14]. Coarsening is often balanced by the eventual decrease in motion (which typically happens prior to the annihilation of all grains but one) is referred to as stagnation. Both mechanical coarsening and stagnation have yet to be thoroughly explained due to the complex interplay between material properties, GB properties, and elastic driving force [15,16].
Microstructure evolution is primarily mediated through the motion of grain boundaries (GBs). GB migration is a complex, non-equilibrium, dissipative process that exhibits a vast range of behaviors depending on crystallography and loading conditions, and there is still little consensus on what constitutes intrinsic GB migration properties [17]. Current knowledge of GB migration behavior has been built up predominantly through atomistic simulations [18][19][20][21][22]. The seminal work of Cahn and Mishin [18] established the existence of shear coupled boundary motion as a means for GBs to mediate permanent deformation. Since then, shear coupled boundary migration has been quantified for a broad range of materials and boundary characters [23][24][25][26][27]. Shear coupled GB motion is generally understood to be mediated through the nucleation, propagation, and annihilation of disconnections, boundary defects that carry both a Burgers vector and a step height [28][29][30][31].
Although MD is the de facto method for determining boundary migration properties, molecular dynamics is fundamentally limited in its ability to treat problems with large (mesoscale) length and time scales. This makes it infeasible to scale atomistic methods to large microstructures, and prompts the need for multiscale modeling. Microstructure modeling at the mesoscale has historically treated boundaries as either immobile (and therefore irrelevant to microstructure evolution problems), or using isotropic properties and curvature-driven flow (typically resulting in unrealistic soap-bubble-like behavior with no nontrivial stable solution) [32,33]. Curvature-driven flow models can be improved with the inclusion of strongly orientation dependent boundary energy, which endows the microstructure with a length scale and produces complex structures such as microfaceting [34,35]. However, even with these enrichments, the full complexity of boundary migration behaviors such as GB mediated plasticity or GB stagnation is not fully captured.
Recently, several alternative approaches to modeling GB migration at the mesoscale have been proposed. By explicitly accounting for disconnections and their modes, it was shown that GB migration can be captured with greater accuracy using a continuum level model [22,23,[36][37][38]. Alternatively, the Kobayashi-Warren-Carter (KWC) method [39][40][41], which has recently been coupled to crystal plasticity [42][43][44], was shown to be able to successfully capture shear coupling as well as grain rotation, although it is limited in that it relies on a simplified model for grain boundary energy. At the larger scale, recent work proposed that multiple modes of GB migration can be accurately captured using the principle of minimum dissipation potential [45], suggesting that GBs have intrinsic properties that can be compactly represented in a single dissipation potential and extended to continuum scales.
By combining the dissipation potential model with the traditional multiphase method, it was shown that disconnection-like structures can arise as a consequence of the dissipation potential and nonconvex GB energy [46]. The work presented here builds on this central idea: modeling disconnections as emergent phenomena that arise as a consequence of nonconvex energy minimization and the principle of minimal dissipation potential. Consequently, the contribution of this work is twofold. First, a number of enhancements to the original phase field disconnections model are proposed to account for thermal nucleation of disconnections, and to extend thresholding to curvature terms as well as elastic terms. Second, the enhanced model is applied to a variety of boundaries and GB migration loading conditions. It is used to investigate the behavior of single disconnections, disconnection-mediated relaxation, and thermally activated disconnection pairs under shear loading.
The remainder of the paper is structured as follows. In Section II, the model is developed starting from the principle of minimum dissipation potential, and integrating the multiphase field model, special linearized elasticity, grain boundary energy, second order regularization, and disconnection nucleation. The computational methods are described in Section III. The results are presented in Section IV for three copper boundaries subjected to three types of loading: a single disconnection pair, a sinusoidal interface relaxation, and thermally-activated shear coupling. The results of the model are summarized, along with the model limitations, in Section V.
Notation used in this work is generally standard. Tensor equations are expressed in invariant notation, with ∇ indicating the gradient, ∇· the vec-tor/tensor divergence, ∇ 2 the Hessian, and ∆ the Laplacian. Functions with arguments indicated using braces (e.g. W [u]) should be understood as functionals that depend on the argument and its spatial or temporal derivatives.

II. DIFFUSE BOUNDARY MODEL
This section presents a model for disconnectionmediated grain boundary migration that builds on the phase field disconnections model presented in [46]. The premise of the model is that disconnections are natural mediators of GB migration, arising automatically as a consequence of the concerted nonconvexity of grain boundary energy as well as elastic energy. Therefore, the model does not presume the existence of disconnections; rather, they arise spontaneously as an emergent phenomenon.
The master governing equation is taken to be the principle of minimum dissipation potential. The idea of an extremal principal for non-equilibrium thermodynamics, though an inherently heuristic construction, has proved extremely useful in the modeling of a wide variety of non-equilibrium mechanical processes including plasticity [47][48][49] and viscosity [50,51]. Recently the construction of minimum dissipation potential was applied to planar GBs in an attempt to identify the "intrinsic" GB properties, and it was determined that the so-called "dissipation energy" along with the traditional mobility were able to capture a variety of planar GB migration behaviors at the mesoscale [45]. Following [46], the minimum dissipation potential construction is used: where W is the free energy, φ * is the dual dissipation potential, and u is the elastic displacement field. In [45] the internal variable was simply the scalar interface position z; here, the interface descriptor is replaced with a vector of n order parameters η (where N is the number of grains), where η satisfy the usual properties of a multiphase field model. It should be noted that C 4 continuity on η is required in order to perform second order regularization, and it is also required that the displacement field u satisfy all natural and essential boundary conditions. The free energy W depends on the gradient of the displacement field ∇u, the order parameter η, its gradient ∇η and Hessian ∇ 2 η. It is decomposed into the following components: which are the chemical potential, the boundary energy, the corner penalization, and the mechanical strain energy, respectively. The chemical potential is taken from [52,53] and has the form with µ = 3.26. W M is minimized when exactly one order parameter is unity and the rest are zero. This drives the grain segregation and can be interpreted either as a mixing energy or simply as a Lagrange multiplier. The boundary energy is where n n is the GB normal to grain n along the boundary. The coefficient k(n n ) is generally given by the quadratic mixing rule although for the present work only two grains are used and so k is simply k 12 . The boundary term k mn corresponds to the boundary energy between grains m, n and, again following Moelans et al. [53], is given by where is the diffuse grain boundary width and σ mn is the orientation-dependent, strongly nonconvex grain boundary energy. The calculation of σ mn will be discussed further in section II A. The strong nonconvexity of σ mn is well-known to generate microfacets in the boundary [54]. Numerically this presents a challenge since the facets have no inherent length scale, resulting in numerical instability. The solution is to add an additional curvature penalization, represented here by W C : whereβ is a regularization parameter and K 23 = 1 2 (κ 2 2 + κ 2 3 ), where κ 2 2 and κ 2 3 are the second and third principle curvatures of η, calculated from ∇ 2 η. (In 2D κ 2 3 = 0.) This form for W C is advantageous for the the proposed model, as it depends only on the physical curvature of the boundary and not the diffuse curvature.
The elastic energy is given by the following quartic mixture rule U n (∇u) = 1 2 (I + ∇u − F GB n ) : C n (I + ∇u − F GB n ). (9) Cubic linear elasticity was used, with the elastic modulus tensor C n rotated to correspond to the known crystallographic orientation of the grain. Prior work [46] used a linear mixing rule, which does not sufficiently localize the stress field to produce the necessary driving force for disconnection motion. The form for the elastic strain energy U n follows the "special linear elasticity" convention [45], which uses a second order Taylor expansion around the grain boundary-induced deformation F GB . This allows the computationally advantageous small strain to be used for elasticity calculations, since the deviation from the large F GB deformations is large. The selection of grain boundary deformation tensor F GB is used to encode the shear coupling factor β. In all of the applications considered here, we consider a single value for β, and then let The value of β is calculated by crystallography and, by the prior determination (via atomistic simulations [24,55]) of the most likely deformation mode for the boundary. For any rational GB there is a countably infinite number of shear coupling factors, and it is possible for a boundary to move by one or more of those factors [17,45]. This model currently aims to capture one shear coupling mode only. However, it should be emphasized that the value of β, though determined based on prior experience with GB migration, is not fundamentally a GB property; rather, it is a property arising from the crystallography of the elastic medium. Indeed, the elastic contribution of this model could be modified to account for a broad range of GB type deformations by accounting for the multiplicity of GB-shears. We leave this to future work. The next component of the model is the dual dissipation potential, which models the rate of energy dissipation as a function of the derivative of the order parameter. Here, a second order form is used: where φ 0 is the dissipation energy and φ 1 = 1/L is the rate-dependent coefficient, defined to be in which M is the GB mobility and is the diffuse boundary width [53]. The mobility used in this model is constant, but in reality is a function of GB character, orientation, and temperature. This dependence is not explored here. The governing equations for η and u are given by solutions to (1). The equilibrium displacement field u * are given by the Euler-Lagrange equations, which are the governing equations of elasticity subject to appropriate boundary conditions. The computational methods for solving the elasticity equations are given in Section III. The Euler-Lagrange equation for η is given by where δ/δη is the variational derivative and ∂/∂η indicates the subderivative [56]. The subderivative becomes necessary due to the lack of smoothness induced by the φ 0 |η| term in the dissipation potential. Accounting for the form of φ * enables the evolution equation to be expressed in the following way: which amounts to a thresholding scheme for evolving η. Unlike [46], which split η into elastic and inelastic components, here, both are evolved together. It was determined that there was no appreciable affect on the diffuse boundary. In the case where φ 0 → 0, the familiar Allen-Cahn equations are recovered. The variational derivatives for the chemical potential and elastic energy portions of the free energy are easily computed. It has been shown that the variational derivative for the boundary energy term in 2D reduces to where θ is the orientation of the normal vector, and κ 2 is the second principal curvature. The variational derivative for the second order regularization term is wherex 1 ,x 2 ,x 3 are the coordinates in the basis corresponding to the principal curvatures of η. In the sharp interface limit in 2D, where the eigenbasis corresponds to the angle θ n , the fourth derivatives are computed to be The above derivatives can be readily extended to 3D when working in the eigenbasis of ∇ 2 η [34]; however, the 2D implementation is used for all examples considered here.   Table I.

A. Grain boundary energy
The energy of a grain boundary is a complex function of the five-dimensional space of macroscopic grain boundary character. GB energy in FCC materials exhibits many sharp minima ("cusps") that can be highly localized in both the space of orientation relationships and interface inclinations. These cusps are of preeminent importance at small scales, where they become the primary driver of interface morphology. In this work, it is necessary to compute the boundary energy as a function of interface normal, σ(n), in order to account for GB orientation dependence.
Over the past several decades, a multitude of GB energy calculation methods have been proposed ranging from the analytic expressions of Read and Shockley [58] to large-scale molecular statics (MS) simulations [24]. The lattice-matching model for GB energy is an analytic model that uses optimal transportation theory and lattice geometry to estimate boundary energy [59]. Lattice matching is particularly useful for bicrystal configurations in which it is difficult to obtain a periodic unit cell for MS simulations, which is often the case when calculating orientation-dependent GB energy for high Σ boundaries.
Because the GB energy exhibits sharp cusps, the curvature σ (θ) is generally infinite at those points. This creates numerical issues and discretization dependency due to the dependence of the η evolution equation (16)  The GB energy is mollifed using a Gaussian (standard deviation 2.5 • ) to regularize the energy. After regularization, the energy is also rescaled and offset using coefficients C 1 , C 2 , respectively, to account for the change in cusp magnitude resulting from the regularization.

GB1 GB2 GB3
Tilt axis 0 0 1 0 0 Three different boundaries in Copper are considered. All are symmetric tilt GBs; asymmetric tilt GB migration behavior is substantially more complex [60][61][62] and out of the scope of the present work. Boundary 1, 2 and 3 have sigma values Σ13, Σ37, and Σ7. Boundary 1 and 2 are 0 0 1 tilt boundaries with {5 1 0} and {7 5 0} boundary planes, respectively, and Boundary 3 is a 1 1 1 tilt boundary with an {1 2 3} boundary plane. Previous work has examined these and other boundaries using both molecular dynamics [63] and a disconnection-based model for grain boundary migration [17,23,36]. Parameters for these values were determined manually and from literature (Table II).
Using lattice-matching, the boundary energy is calculated for each boundary ( Figure 1). An angle of 0 • corresponds to an unrotated interface. All of the boundaries exhibit cusps at 0 • , as expected, although the magnitude of the cusps varies substantially. The 1 1 1 {1 2 3} STGB has a very small cusp at 0 • , but very large cusps at ±90 • . On the other hand the 0 0 1 {7 5 0} boundary exhibits a moderate cusp at 0 • but a multiplicity of cusps at a variety of other angles. The differences in these energy landscapes will contribute to the differing types of disconnections observed in the phase field model results.

B. Nucleation model
As with plasticity by dislocation motion, boundary motion relies on disconnection nucleation in addition to propagation. Disconnections can nucleate in many ways; grain boundary dislocations [64,65], triple junctions [66], or other interfaces [67] can act as disconnection sources. However for large, pristine boundaries, disconnection nucleation becomes primarily a thermally activated process. To capture this, a thermally activated mechanism for disconnection nucleation must be introduced.
Due to the large energy barrier inherent to disconnection formation, the phase field model is not capable of generating nucleation events spontaneously. Therefore, many phase field models include a zeromean noise term to provide perturbations that mimic the effect of thermal fluctuations, making it possible to access high-energy states. On the other hand, in this present work, it was determined that the energy barrier to disconnection formation was high enough so that thermal noise was insufficient to perturb the boundary to form a disconnection pair.
In general, it was found that the boundary invariably destabilized before it was able to spontaneously form any disconnections.
In order to capture disconnections, therefore, a new model for explicit pair nucleation is proposed. Introduce a probability density function p(x, t) that is the probability of disconnection formation within a region with volume v 0 over a characteristic time interval τ . Let p : where k B is Boltzmann's constant, T is the temperature, and E an energy given as a function of space and time. This energy is now connected to the phase field model through the following: E 0 is a numerical parameter corresponding to the creation of a disconnection pair. η 1 , η 2 are the order parameters (in this treatment we consider a bicrystal system only). is a very small numerical parameter (O(10 −20 ) to prevent division-by-zero errors, and n is a positive integer. The form of (20) is designed to restrict nucleation events to the boundary, which is characterized by the existence of non-zero values for η 1 , η 2 . The value for n determines the width of the region in which nucleation events are allowed.
Finally, E → ∞ if the order parameter deviates outside the [0, 1] range. This prevents numerical instabilities that can sometimes arise due to grain boundary energy orientation dependence.
As stated in its initial definition, p is the probability of nucleation with respect to a pre-defined, arbitrary spatial region and time interval. In practice, this must be converted to an effective probability for a pre-defined, possibly larger (or smaller) region of time and space. Because an adaptive mesh is used, along with temporal subcycling, it is essential to calculate probabilities in a manner that is consistent between various levels of refinement in space and time. Towards that end, the following heuristic is introduced to estimate the probability of nucleation for some arbitrary region B ⊂ Ω over an interval [a, b] ⊂ R: P is thus a functional on R 3 ×R corresponding to the probability of a nucleation event occurring in a certain domain and over a certain interval. We pause briefly to point out some of the interesting mathematical properties of P . Although P is decidedly not a measure or a distribution, it does possess some analagous properties: In other words, 0 ≤ P ≤ 0, and P possesses the particularly useful property that the probability of nucleation P over a region corresponds to the combined probability of nucleation P in the corresponding subsets. The above definitions are now used to construct an algorithm for nucleating a disconnection. Let x 0 be a grid point on a level with timestep ∆t and (in two dimensions) grid spacings ∆x 1 , ∆x 2 . At each timestep, the probability P of disconnection nucleation can then be computed, which is using a quadrature rule to compute the integral. (It should be noted that on fine levels, p varies smoothly enough such that single-point quadrature is sufficient; on coarse levels, p(x 0 , t) has been pre-averaged through the AMR process.) A uniform distribution random number is then computed, and if P exceeds that number, then a nucleation event occurs. Another random number is generated to determine the permutation of the nucleation, and then the order parameters are modified thus (assuming a grain 1 permutation): The mappings are naturally permuted depending on the type of pair that was activated.
Nucleation energy [68] E0 = 0.0457 Regularization parameter = 10 −20 Reference area-time τ v0 = 10 −3 Nucleation width κ = 0.005 Nucleation model parameters are either determined by calibration or from literature (Table III). Temperature is treated as an input, so that the effect of temperature can be determined in the examples.

III. COMPUTATION
All phase field methods require sufficient grid resolution in order to capture diffuse boundaries without mesh dependency, usually requiring four to eight grid spacings across the boundary [69]. This can induce a computational bottleneck, since most diffuse boundary models are considered accurate only in the limit as the diffuse boundary width goes to zero. Since one of the primary advantages of the proposed method is its ability to scale well beyond the spatial and temporal timescales of molecular dynamics, it is necessary to employ methods to increase performance and eliminate wasteful computations.
This work uses a block-structured adaptive mesh refinement (BSAMR) strategy to selectively refine the grid near the grain boundary. Unlike other AMR methods (such as quad/octree division), BSAMR treats each refinement level completely separately, eliminating the need to explicitly track connectivity information. Each refined region overlays a coarse region, and both are evolved independently and then periodically synchronized by averaging the fine region onto the coarse region. For explicit time integration, there are two main advantages. First, the block-structured nature of the grid allows for memory-efficient data organization, and the minimal amount of connectivity reduces the amount of parallel overhead for very efficient parallel scaling. Second, it enables the use of temporal subcycling, so that each level has its own timestep and fine levels experience several iterations for each iteration experienced on the coarse level. This is particularly useful in the present work, for which the temporal integration involves a fourth order spatial derivative. By using a subcycling ratio of 32, the mesh can be refined arbitrarily without violating the CFL condition.
The implicit component, i.e. the solution to the elasticity equations, requires special treatment on a BSAMR grid. Unlike most elastic solvers, this method solves the strong form of the elasticity equations directly, using a finite difference discretization. The strong form method, when combined with BSAMR, is advantageous for implicit solutions using the geometric multigrid method, as the structure of coarsened levels can be situated tidily within the framework of the AMR refinement layers. Additional care is needed to treat the boundary between levels ("coarse-fine boundary") during the elastic solve; here, the "reflux-free" method, proposed in [70,71], is used.
A staggered method for solving the implicit (elasticity) equations and explicit (η) equations is used. Because of the small timestep required due to the fourth-order spatial derivative, it is determined that where r is the refinement theshold and ∆V the AMR grid size.

IV. RESULTS
In this section, three examples of disconnection migration are considered. In all cases, the ground state is a 2D plane strain bicrystal occupying a domain with dimensions 8x8 in nondimensionalized units. The top and bottom grains have eigenstrains F GB 1 , F GB 2 as described in the methods section, where the F GB encodes the shear coupling factor in the model. Neumann boundary conditions are used on the left and right faces for both the order parameters η 1 , η 2 and the displacements u; Dirichlet conditions are prescribed for the top and bottom. In order to avoid numerical instability, each simulation begins with isotropic curvature-driven flow until t = 1; this is necessary to avoid instabilities resulting from the initially very high interface sharpness. At t = 1 anisotropy is enabled, and the timestep is decreased from dt = 0.05 to dt = 0.0005; then at t = 1.1, elasticity is enabled and a positive shear displacement is applied to the top face at a prescribed rate.

A. Single disconnection pair
In this section the behavior of a single disconnection pair under an externally applied loading is considered. Starting with a bicrystal, a single nucleation pair is generated at the beginning of the simulation. In each case, two possible disconnection pairs are considered: "up-down" (a disconnection pair extending into the top grain) and "down-up" For the two boundaries with positive coupling factor, an "up-down" disconnection pair repels, whereas the "down-up" disconnection pair annihilates. On the other hand, the boundary with negative coupling factor exhibits the reverse behavior. Because of the small disconnection size, the "down-up" pair additionally dissociates into smaller disconnections as they move.
(extending into the bottom grain). The former corresponds to shear coupling action that moves the boundary up; the latter, to moving the boundary down. Initial tests without the presence of external shear shows that the disconnection pair either annihilated or did not move; as expected. The disconnection pairs are then subjected to an external load. The three bicrystals are sheared with an initial value of 0.0125 that linearly increases to 0.05 over an interval of ∆t = 10.
For the "up-down" nucleation in the 0 0 1 {5 1 0} boundary, the disconnections immediately move away from each other. The rate at which they move is dependent on the elastic driving force from the shear deformation. The profile of the disconnections is constant over the course of their motion, and there is no motion of the boundary due to curvature. On the other hand, the "down-up" nucleation produces a very different result: the two nucleations move towards each other and annihilate immediately. Both of these results are consistent with the known behavior of disconnections in the 0 0 1 {5 1 0} boundary, and the resultant net positive shear coupling.
The second boundary, 0 0 1 {7 5 0}, has a negative coupling factor. When the "up-down" disconnection pair was nucleated, the driving force (combined with the general attraction) causes the pair to annihilate almost immediately. On the other hand, the "down-up" pair repels for a net effect of moving the interface downward. Unlike the other two cases, the 0 0 1 {7 5 0} downward disconnection pair dissociates into four smaller disconnections. This difference in behavior can be attributed to the different energy landscape for this particular boundary: whereas the other two boundaries have only two cusps at non-0 • locations, Boundary 2 has cusps at θ = ±45 • as well as at θ = 0 • , ±90 • . Therefore, there is less energetic cost associated with disconnection formation. This phenomenon of "disconnection dissociation" into what might be called "partial disconnections" has not, to our knowledge, been observed in molecular dynamics or experiment. Therefore we leave this as a model prediction, urging caution in the interpretation of these results until they can be confirmed by MD or experiment.
Finally, the 1 1 1 {1 2 3} boundary (Boundary 3), which has a coupling factor, exhibits similar behavior to Boundary 1. The primary difference between the two is that Boundary 3 creates much smaller disconnections. This effect can again be attributed to the difference in grain boundary energy landscape.
These results demonstrate that motion by steps, which we identify as "disconnections" is the mechanism by which the boundary moves. In all of the cases presented here, the planar boundary is immobile; only the steps are able to move. This illustrates the proposed mechanism by which disconnections propagate: steps induce stress that in turn increase the driving force above the threshold value. It should also be noted that, unlike in [46], there is no curvature-driven motion. This is a consequence of the updated flow rule that thresholds the curvature terms as well as the elastic driving force terms.
It is important to clarify that the step heights of the disconnections for the different boundaries is not directly connected to the crystallography of the bicrystal, but are governed by the model's corner energy parameter. The Burgers vector, which is equal to the (prescribed) coupling factor times the step height, is thus similarly determined by the corner energy in the model. In other words, by adjusting the corner energy, both the disconnection step height and Burgers vector will change accordingly (as seen in [34]) while maintaining a constant shear coupling factor. Consequently the corner energy coefficient is effectively a modeling parameter which could be varied between boundaries in order to more accurately capture the physical corner energy. Because model does not distinguish between facet corners and disconnection corners, additional enrichment will be required to capture the interaction between facet corners and disconnections. This, and the calibration of the model's corner parameter to physical corner energy (c.f. [73]) are left to future work.

B. Relaxation of sinusoidal perturbation
The three boundaries are subjected to an initially sinusoidally perturbed interface to determine the effect of disconnection migration on interface relaxation. Although no load is applied, the primary driving force in all three cases was the elastic mismatch. Each boundary is allowed to evolve until T = 40 ( Figure 3). The temperature for all simulations is sufficiently low so that no nucleation events occurred during the interval, beyond those spontaneously generated from the initially smooth cosine boundary.
All boundaries relax initially due to the substantial curvature at the extrema of the boundary, but the curvature-driven motion quickly stagnate. Subsequent motion of the interface is driven primarily by the motion of disconnections. The 0 0 1 {5 1 0} boundary exhibits a variety of step sizes. Generally, large steps move slowly and eventually dissociate into smaller, faster moving steps.
The second boundary 0 0 1 (7 5 0), unlike the first, retains a faceted configuration even after motion has begun to stagnate. This is unsurprising, as the persistence of angled facets is almost certainly due to the presence of sharp cusps at ±45 • as well as 0 • and ±90 • (Figure 1). This boundary is the only boundary to have a negative coupling factor, but as expected, there is no apparent effect on the evolution.
The third boundary 1 1 1 (123) differs from the first two in that it relaxes substantially faster. This is likely due in part to the greater amount of initial curvature-driven flow resulting from the higher GB energy at θ = 0 • . The boundary also exhibits small cusps, similar to those for the previous case, and consistent with the single pair nucleation results. An aspect of particular interest is the transient early-time corners at the extrema, which is markedly different from the usual self-similar curvature-driven flow. In fact, the results appear qualitatively similar to those recently obtained by Zhang et al. using a continuum model for GB migration [37], although the referenced work was considering disconnections with no shear character.
In all cases, the boundaries remain nearly perfectly flat between steps. This is particularly apparent in the trace plots, since they are exaggerated in the vertical direction and virtually no slope is visible. It is also of interest that all of the boundaries appear to stagnate well before they reach their final equilibrium position. The implication is that the relaxation is due entirely to disconnection motion, and that the driving force is purely elastic. Also, it appears that the energy barrier (dissipation energy) is sufficient to hold these boundaries in a nontrivial meta-stable state. Finally, the substantial difference between motion-by-curvature and motion by elastically-driven disconnection clearly indicates the importance of accounting for this mechanism at the mesoscale.

C. Thermally-activated shear coupling
Shear coupling has been a subject of interest in microstructure evolution for multiple decades, and is one of the motivating examples of the present work. It is of particular interest in mechanical modeling because it is a dissipative mechanism by which permanent microstructure evolution takes place (i.e. plasticity), yet it is distinct from the traditional mechanism of dislocation mediated plasticity. Some have coined the mechanism "grain boundary-mediated plasticity"; one might also think of it as "disconnection-mediated plasticity." Indeed, there are a number of analogues between crystal plasticity and grain boundary plasticity [45] that are consistent with this work.
As discussed in the introduction, it is believed that disconnections act as mediators of grain boundary migration. It has been observed in molecular dynamics simulations that both the simulation domain size and the temperature substantially affect the shear coupling behavior. This may be the result of a number of possible thermally activated mechanisms, including the activation of alternate disconnection modes (and corresponding shuffle patterns) [74]. Since disconnection nucleation is a thermally activated mechanism, it is also a mechanism by which temperature can affect migration propensity. Indeed, the earlier example (Section IV A) illustrates that without the nucleation of disconnections, the boundary is immobile (at least up until the elastic driving force exceeds the dissipation potential even without disconnection-induced stress concentration).
Thermal nucleation of disconnection pairs is certainly not the only mechanism by which they can form; there are myriad other potential disconnection sources. The sinusoid relaxation example (Sec-tion IV B) is one possible example of this. However, for a large pristine boundary, such as considered in most molecular dynamics simulations, thermal nucleation is the primary means of disconnection generation. Such boundaries are the most natural examples against which this work should be compared.
All three boundaries are subjected to a displacement-driven shear stress. The prescribed shear strain begins at zero and is linearly increased to 0.125 over a time interval of ∆t = 38 (Recall that the strain, which would normally be too large to be valid in traditional elasticity, can be acceptable when working with special linearized elasticity.) The simulation was repeated for each boundary with temperatures ranging from T = 50 to T = 800. Boundary 1, as expected, exhibits strong positive shear coupling under the applied load. As the boundary is subjected to the increasing elastic driving force, it moves upwards to relieve the stress. Examination of the boundary migration clearly shows that the mechanism for motion is the creation, motion, and annihilation of steps ( Figure 5). While under little or no stress, disconnections generally annihilated, or in some cases persists for a few timesteps. As stress increases, "up-down" disconnection pairs (corresponding to shear coupling action that moves the boundary up) began to repel sporadically, exhibiting a kind of stick-slip behavior. Eventually, with increasing stress, all nucleated disconnection pairs begin moving quickly and with increasing speed over the course of the simulation. For high temperature simulations, near the end, disconnections move and annihilate so quickly that they become indistinguishable. A detail of the boundary migration is presented here (Figure 4), and the reader is referred to the supplementary material for an animation of the boundary migration along with the accompanying stress-strain curve.
The stress-strain curves measured for Boundary 1 (Figure 6) indicate the effect of thermally activated disconnection pair nucleation on boundary migration. Boundaries that are "cold" (T = 50, 200) nucleate only a couple of disconnection pairs. As a result, the boundary is immobile until the stress reached σ = 0.8. At this point, the driving force exceeds the dissipation energy value and causes the boundary to move even without disconnections. This is manifested as "elastic-perfectly-plastic" behavior, and would correspond to grain boundary sliding. Generally, most of the boundaries below 200 all exhibit the elastic-perfectly-plastic behavior. A notable exception is the boundary at T = 200, which slightly exceeds the yield stress of the T = 50 case. Upon inspection, this is the result of a single disconnection that nucleates at ε = 0.125 and induces a small amount of curvature, which effectively slows the boundary motion at that time. As the temperature increases, the maximum attained stress reduces as well and the stress-strain curves exhibit softening similar to thermal softening in plasticity. As disconnections nucleate more frequently, it naturally becomes easier for the boundary to move, resulting in a lower overall yield stress. The reader is again referred to the supplementary material for additional visualizations of Boundary 2 motion. Boundary 2 exhibits negative shear coupling, as expected ( Figure 7). Behavior is generally similar to that of Boundary 1, with the main difference being the sign of the dominant disconnection pairs and also the size of the disconnections themselves. Because of the different energy landscape, there is a lower energetic penalty for the creation of steps, and so a greater number of steps is created resulting in a net greater curvature of the boundary. Moreover, because of the reduced size of the disconnections, each nucleation event effectively creates a greater number of smaller disconnections.
The most substantial difference observed in the migration of Boundary 2 is the substantially lower yield stress (Figure 8). Whereas the maximum yield stress in Boundary 1 ranges from 0.06 to 0.09, the yield stress in Boundary 2 ranges from 0.03 to 0.06. (Both boundaries are subjected to the same loading conditions.) Furthermore, Boundary 2 moves considerably further than Boundary 1 under the applied strain. In fact, for all of the cases with T > 200, the boundary reaches the end of the domain, ending the motion and causing the stress curves to collapse onto a single line. There are a number of factors that affect the yield stress, including the crystallographic elastic anisotropy as well as the grain boundary energy. However, in this case, the predominant factor is the shear coupling factor β which is lower for Boundary 2 (|β| = 0.33) than for Boundary 1 (|β| = 0.4). The lower coupling factor means that Boundary 2 must move a greater distance than Boundary 1 in order to relieve the same amount of stress; therefore, the driving force from the same elastic loading condition is greater on Boundary 2 than for Boundary 1. The motion of Boundary 3 was sufficiently similar  T=50  T=100  T=200  T=300  T=400  T=500  T=600  T=800 FIG. 9. Stress-strain data for Boundary 3 ({1 1 1}(1 2 3)) at various temperatures (qualitatively) to that of Boundary 1 that visualizations of its motion are not included here; however, animations of the motion of Boundary 3 are included in supplementary material. The behavior of Boundary 3 is entirely consistent with the conclusions drawn from the motion of the prior two boundaries. The positive coupling factor (β = 0.69) induced upwards motion of the boundary, as expected. For low temperatures, a nearly elasticperfectly-plastic behavior is observed (Figure 9), although it appears that disconnection nucleation has a similar effect as the T = 200 case for Boundary 1disconnections induced curvature that temporarily obstructed boundary motion.
Because of the larger shear coupling factor for Boundary 3, a substantially higher yield stress is observed, ranging from 0.6 to 1.4. This, and the fact that the boundary moved a relatively small amount under the prescribed strain, is consistent with the above explanation.

V. CONCLUSIONS
This work aims to understand and model disconnection-mediated grain boundary migration as an emergent phenomenon. Rather than building disconnections into the model as a fundamental entity, they arise as a consequence of energy nonconvexity and the principle of minimum dissipation potential. The model combines several mesoscale modeling approaches: the multiphase field model for microstructure, the minimum dissipation potential for grain boundaries, the lattice-matching model for nonconvex grain boundary energy, and the higher order regularization for faceted boundaries. The model is applied to three copper symmetric tilt grain boundaries, and three cases are considered: (i) nucleation of a single disconnection pair, (ii) relaxation of a sinusoidally perturbed boundary, and (iii) thermally activated shear coupling. In all three cases, the results are consistent with experimental and atomistic observation of boundary migration via disconnections.
A number of simplifications were made in the development of this model. For each boundary con-sidered here, only one disconnection mode was enabled. This improves the performance and ease of implementation for the model, but artificially limits the types of migration that the boundary can experience (in particular, multi-mode migration as observed in [22]). It should also be noted that the mobility of the boundary was taken to be constant across all three boundaries and also with respect to boundary orientation as well as temperature. By accounting for character and temperature-dependent mobility it may be possible to improve the model's predictions. Temperature effects were also neglected for the elastic moduli for the grain boundary energy, both of which can change substantially over the temperature range examined. More generally, the model does not account for non-stress driven disconnection migration. It also does not yet account for other mechanisms such as grain boundary sliding or grain rotation, both of which are important effects in microstructure evolution.
The authors gratefully acknowledge Lawrence Berkeley National Laboratory (LBL) subcontract #7473053, which supported the development of the computational methods used in this work. In addition, MG acknowledges support from the California Institute of Technology Summer Undergraduate Research Fellowship (SURF) for Summer 2019 and Summer 2020.