Mechano-chemical spinodal decomposition: A phenomenological theory of phase transformations in multi-component, crystalline solids

We present a phenomenological treatment of diffusion-driven martensitic phase transformations in multi-component crystalline solids that arise from non-convex free energies in mechanical and chemical variables. The treatment describes diffusional phase transformations that are accompanied by symmetry breaking structural changes of the crystal unit cell and reveals the importance of a mechano-chemical spinodal, defined as the region in strain-composition space where the free energy density function is non-convex. The approach is relevant to phase transformations wherein the structural order parameters can be expressed as linear combinations of strains relative to a high-symmetry reference crystal. The governing equations describing mechano-chemical spinodal decomposition are variationally derived from a free energy density function that accounts for interfacial energy via gradients of the rapidly varying strain and composition fields. A robust computational framework for treating the coupled, higher order diffusion and nonlinear strain gradient elasticity problems is presented. Because the local strains in an inhomogeneous, transforming microstructure can be finite, the elasticity problem must account for geometric nonlinearity. An evaluation of available experimental phase diagrams and first-principles free energies suggests mechano-chemical spinodal decomposition should occur in metal hydrides such as ZrH$_{2-2c}$. The rich physics that ensues is explored in several numerical examples in two and three dimensions and the relevance of the mechanism is discussed in the context of important electrode materials for Li-ion batteries and high temperature ceramics.


I. INTRODUCTION AND FORMULATION
Spinodal decomposition is a continuous phase transformation mechanism occuring throughout a solid that is far enough from equilibrium for its free energy density to lose convexity with respect to an internal degree of freedom. The latter could include the local composition as in classical spinodal decomposition described by Cahn and Hilliard [1], or a suitable non-conserved order parameter as in Allen and Cahn's theory for spinodal ordering [2]. A key requirement for continuous transformations is that order parameters can be formulated to uniquely describe continuous paths connecting the various phases of the transformation. These phases then correspond to local minima on a single, continuous free energy density surface in that order parameter space. For classical spinodal decomposition inside a miscibility gap, all phases have the same crystal structure and symmetry, and the order parameter is simply the local composition. The existence of a single, continuous free energy density surface for all phases participating in a transformation implies, by geometric necessity, the presence of domains in order-parameter space where the free energy density is non-convex. Reaching those domains through supersaturation (by externally varying temperature or composition) makes the solid susceptible to a generalized spinodal decomposition.
Many important multi-component solids undergo phase transformations that couple diffusional redistribution of their components with a structural change of the crystallographic unit cell. One prominent example is the decomposition that occurs when cubic yttria-stabilized zirconia Zr 1−c Y c O 2−c/2 is quenched into a two-phase equilibrium region between tetragonal Zr 1−c Y c O 2−c/2 having low Y composition and cubic Zr 1−c Y c O 2−c/2 having a high Y composition. Another occurs in Li battery electrodes made of spinel Li c Mn 2 O 4 . Discharging to low voltages causes the compound to transform from cubic LiMn 2 O 4 to tetragonal Li 2 Mn 2 O 4 through a two-phase diffusional phase transformation mechanism. As with simple diffusional phase transformations, these coupled diffusional/martensitic phase transformations can either occur through a nucleation and growth mechanism, or, if certain symmetry requirements are met, through a continuous mechanism due to an onset of an instability with respect to composition and/or a structural order parameter.
Here we present a treatment of coupled diffusional/martensitic phase transformations triggered by instabilities with respect to both strain and composition. These phase transformations are characterized by a mechano-chemical spinodal that is defined as a nonconvex region of the free energy density function in the strain-composition space. The possibility of a mechanochemical spinodal decomposition is motivated by recent first-principles studies of martensitic transformations in transition metal hydrides where a high temperature cubic phase is predicted to display negative curvatures with respect to strain, thus making strain a natural order parameter to distinguish the cubic parent phase from its low temperature tetragonal daughter phases [33]. In addition, the coupling with composition degrees of freedom (e.g. through the introduction of hydrogen vacancies in the metal hydrides) allows for the possibility that the free energy also exhibits a negative curvature with respect to the composition. Mechano-chemical spinodal decomposition, therefore, is a phenomenon that is likely present in many multi-component materials but has to date been overlooked as a mechanism by which a high symmetry phase can decompose martensitically and through diffusional redistribution upon quenching into a two phase regime. Structural phase transformations in solids driven by instability with respect to an internal shuffle of the atoms within the unit cell have been treated rigorously in the literature with coupled Cahn-Hilliard and Allen-Cahn approaches [5,6,7,8]; but mechano-chemical spinodal decomposition is fundamentally different and necessitates a coupled treatment of both the strain and composition instabilities.
Our treatment is based on a generalized, Landautype free energy density function that couples strain and composition instabilities. The governing equations of mechano-chemical spinodal decomposition, obtained by variational principles, generalize the classical equations of the Cahn-Hilliard formulation [1], and of nonlinear gradient elasticity [9,11] by coupling these systems. The ability to solve this complex, nonlinear, strain-and composition-gradient-driven, mechano-chemical system for sufficiently general initial and boundary value problems in two and three dimensions also has been lacking heretofore. We introduce the computational framework to obtain such solutions in general, three-dimensional solids. This newfound capability allows us to then reinforce our discussion of the phenomenology with dynamics predicted by the numerical solutions.

A. The mechano-chemical spinodal in two dimensions
For accessibility of the arguments, we first consider the two-dimensional analogue of the cubic to tetragonal transformation: the square to rectangle transformation. The high-symmetry square lattice will serve as the reference crystal relative to which strains are measured. Lower symmetry lattices that can be derived from the square lattice by homogeneous strain include the rectangle, the diamond and lattices where there are no constraints on the cell lengths and their angles.
We use the composition c, varying between 0 and 1, as our order parameter for the chemistry of our binary two-dimensional solid. Symmetry breaking structural changes are naturally described by the strain relative to the high symmetry square lattice. The elastic free energy density is also a function of strain. In both contexts, strain will in general be of finite magnitude. Following Barsch & Krumhansl [9], we use the Green-Lagrange strain tensor, E, relative to the square lattice. Rotations are exactly neutralized in this strain measure; for any rigid body motion, E = 0 (Supporting Information). In two dimensions, the relevant strain components are E 11 , E 22 and E 12 = E 21 . However, it is more conve- (a) Square to rectangle transition (2D) in reparametrized strain space. (b) Cubic to tetragonal transition (3D) in reparametrized strain space (Equation [1]). Lattice vectors are labelled by the corresponding coordinate directions X 1 , X 2 , X 3 and the corresponding lower symmetry phases are labelled 1, 2, 3. The deformations shown in the sub-figures are area/volume preserving, respectively. nient to use linear combinations of these components that transform according to the irreducible representations of the point group of the high symmetry square lattice. In two dimensions, these include e 1 = ( Here, e 1 and e 6 reduce to the dilatation and shear strain, respectively in the infinitesimal strain limit. The strain measure e 2 uniquely maps the square lattice into the two rectangular variants [ Figure 1a]: positive and negative e 2 generate the rectangles elongated along the global X 1 and X 2 directions, respectively. The equivalence of the rectangular variants under the point group symmetry of the square lattice (e 2 = 0) restricts the free energy density to even functions of e 2 .
If the crystalline solid has multiple chemical species, its free energy density dependence on e 1 , e 2 and e 6 can change with composition, c. Figure 2a illustrates a free energy density surface, F (c, e 2 ), for a binary solid that, at higher temperature, forms a solid solution having square symmetry. In this case F is convex, a condition made precise by specifying positive eigenvalues of its Hessian matrix over the {c, e 2 } space. Additionally, F (c, 0) is a minimum with respect to e 2 for fixed c, making the square phase stable for all c at this temperature. However, at a lower temperature, F may lose convexity, inducing the notion of a mechano-chemical spinodal region. We define it as the domain in {c, e 2 } space over which the Hessian matrix admits negative eigenvalues, as illustrated in Figure 2b. Here, we focus on the conditions ∂ 2 F /∂c 2 < 0 and ∂ 2 F /∂e 2 2 < 0. The square phase (e 2 = 0) remains stable at high composition (c ∼ 1) with positive eigenvalues of the Hessian (Figure 2b). But it is mechanically unstable at low composition, with ∂ 2 F /∂e 2 2 < 0 for (c, e 2 ) ∼ (0, 0), and has two symmetrically equivalent, stable, rectangular phases with ∂ 2 F /∂e 2 2 > 0 at (c, e 2 ) = (0, ±s e ). While not shown in Figure 2b we assume that F is convex with respect to e 1 and e 6 . Figure 2c illustrates a schematic temperaturecomposition phase diagram consistent with the free energy densities of Figure 2a-b. The square/cubic phase β, forms at high composition or at high temperature; the rectangle/tetragonal phase α, forms at low composition and low temperature. A large two-phase region separates them.
Zuzek et al. [32] have obtained phase diagrams experimentally that are topologically equivalent to Figure  2c for the ZrH 2−2c system. Recent first principles calculations [33] have demonstrated the existence of a mechanical instability that exists in this system at low c via non-convexity with respect to strain. Furthermore, the two-phase coexistence at low temperature also implies non-convexity with respect to composition as we have demonstrated in supporting information. Figures 2a and b therefore represent a two-dimensional analogue of the free energy for such systems. Consider our model binary solid, annealed at high temperature T h , to form a solid solution in the square phase, β. Its state is at point A in Figure 2a, with e 2 = 0. It is then quenched into the two-phase region (Figure 2c) with free energy density at point B in Figure 2b. For a quench at sufficiently high rate, the dimensions of the square lattice, controlled by strains e 1 , e 2 and e 6 , and the composition remain momentarily unchanged. However, since the state at point B satisfies ∂ 2 F /∂e 2 2 < 0 and ∂ 2 F /∂c 2 < 0, there exist thermodynamic driving forces for segregation by strain and composition within the mechano-chemical spinodal.
Diffusion being substantially slower than elastic relaxation, the solid immediately deforms to either positive or negative e 2 , driven towards a local minimum at constant c. These deformations due to the mechanical instability will happen like many martensitic transformations where a mix of symmetrically equivalent rectangular variants coexist to minimize macroscopic strain energy. For finite but moderate strain, the transformation could proceed coherently, even if the two symmetrically equivalent rectangular variants coexist. We neglect non-essential complexities of this process and assume that, instantly upon quenching, finitely sized neighborhoods of the solid deform homogeneously into one of these rectangular variants at the original composition.
The solid also becomes susceptible to uphill diffusion because ∂ 2 F /∂c 2 < 0 implies a negative diffusion coefficient. However, it does not occur at constant e 2 , since the valleys traversing the local minima, ∂F /∂e 2 = 0, between the square lattice at c = 1 and the rectangular lattices at c = 0 span intervals of negative and positive e 2 . Mechano-chemical spinodal decomposition sets in. Composition modulations are amplified: high c regions strive to be more square (point C) while low c regions strive to be more rectangular (points D or D ). However, since coherency is maintained, some neighborhoods in the solid will be frustrated from attaining strains that ensure minima, ∂F /∂e 2 = 0, for local values of c. Coherency strain-induced free energy penalties arise to alter the driving forces for purely chemical spinodal decomposition. Movie S5 in the supporting information shows the evolution of the state (c, e 2 ) of the material points on the free energy manifold F .

B. Mathematical formulation: three dimensions
Armed with the insight conveyed by the twodimensional study, we next lay out the three-dimensional treatment, in which setting the considered phase transformations proceed via lattice deformation and diffusion. The arguments are made more concrete by considering the general mathematical form of the free energy density.

Strain order parameters
The strain measures e 1 -e 6 are first redefined using the full Green-Lagrange strain tensor components in three dimensions: In the limit of infinitesimal strains, e 1 describes the dilatation, while e 4 , e 5 and e 6 reduce to shears. The point group operations of the cubic crystal leave e 1 invariant since it is the trace of E, while collecting each of the subsets {e 2 , e 3 } and {e 4 , e 5 , e 6 } into a symmetry-invariant subspace whose elements transform into each other. The measures e 2 and e 3 are especially suited as order parameters to describe cubic to tetragonal distortions. All three tetragonal variants that emerge from the cubic reference crystal can be uniquely represented by these two measures. See Figure 1b where tetragonal distortions of the cubic crystal along the X 1 , X 2 and X 3 axes have been labelled as 1, 2 and 3, respectively. Deviations from the dashed lines in the e 2 -e 3 space of Figure 1b correspond to orthorhombic distortions of the cubic reference crystal. This role of e 2 and e 3 as structural order parameters to denote the degree of tetragonality, and to distinguish between the three tetragonal variants, complements their fundamental purpose as arguments of the elastic free energy density.

The mechano-chemical spinodal in three dimensions
For brevity we write the strains as e = {e 1 , . . . , e 6 }. We introduce further phenomenology by specifying that F (c, e) is only one part of the total free energy density. It is a homogeneous contribution whose composition and strain dependence cannot be additively separated. Figure 3, for example, illustrates contour plots of F (c, e) on the e 2 -e 3 and e 3 -c planes for a binary solid having a temperature-composition phase diagram similar to that of Figure 2c. To generate the tetragonal α phase as illustrated in that diagram, the homogeneous free energy density as a function of strain must qualitatively follow the contours of the e 2 -e 3 plane at c = 0 in Figure  3. Here, the tetragonal variants are local minimizers of F in e 2 -e 3 space, with equal minima. In turn, to obtain the cubic β phase, F (c, e) must follow the contours of the e 2 -e 3 plane at c = 1. Thus, F changes smoothly from convex with respect to e 2 and e 3 on the c = 1 plane to non-convex at c = 0 to form the three variants of the tetragonal α phase. Importantly, the planes at c = 0 and c = 1 themselves must be minimum energy surfaces to represent the tetragonal α and cubic β phases, respectively. Movie S8 in the supporting information shows the evolution of the state (c, e 2 , e 3 ) of the material points on the free energy manifold F . Gradient regularization of the free energy density: Following van der Waals [10], and Cahn & Hilliard [1] in their treatment of non-uniform composition fields, we can extend the total free energy density beyond F by writing it as a Taylor series, retaining terms that depend on the composition gradient, ∇c. We extend this gradient dependence to the strain measures ∇e, as did Barsch & Krumhansl [9] following Toupin [11]. (Also see Karatha and co-workers [12,13], for treatments using infinitesimal and finite strain, respectively.) These gradient dependences appear in a non-uniform free energy density G (c, e, ∇c, ∇e). Frame invariance of F and G is guaranteed since the members of e are linear combinations of the tensor components of E. They also must be invariant under point group operations of the cubic reference crystal.

Free energy functional
The crystal occupies a reference (undeformed) configuration Ω with boundary Γ. The total free energy, Π, is the integral of the free energy density F + G over the solid with boundary contributions included. Thus, Π is a functional of the composition c and the displacement vector field u, from which the strains are derived (see Supporting Information): where traction vector component T i is specified on the boundary subset Γ Ti ⊂ Γ. Following the above authors we include only quadratic terms in the gradients, but in a generalization we also allow cross terms between ∇c and ∇e in the non-uniform contribution G , which can therefore be written as G (c, e, ∇c, ∇e) = 1 2 ∇c · κ(c, e)∇c Here κ is a symmetric tensor of composition gradient energy coefficients, each γ αβ (α, β = 1, . . . , 6) is a tensor of strain gradient energy coefficients, and each θ α is a tensor of the mixed, composition-strain gradient energy coefficients. Note that, in general, these coefficients will be functions of the local composition and strain. The point group symmetry of the cubic reference crystal imposes constraints on the tensor components of κ, γ αβ and θ α as well as on the form of F .
While the gradient energies bestow greater accuracy upon the free energy description of solids with nonuniform composition and strain fields, they are essential at a more fundamental level if the homogeneous free energy density is non-convex. At compositions that render F non-convex, the absence of a gradient energy term will allow spinodal decomposition characterized by composition fluctuations of arbitrary fineness, thus leading to non-unique microstructures-a fundamentally unphysical result [14]. With κ = 0, the composition gradient energy penalizes the interfaces wherein composition varies rapidly between high and low limits. This ensures physically realistic results, manifesting in unique microstructures with a mathematically well-posed formulation.
An essentially analogous situation exists with respect to the negative curvatures of F in the e 2 -e 3 plane at low c, which drive the cubic lattice to distort into the tetragonal variants corresponding to the three free energy wells. Consider a solid with a homogeneous free energy density as in Figure 3 and a strain state lying between the valleys in the e 2 -e 3 plane at c = 0. Absent the strain gradient energy, mechano-chemical spinodal decomposition would allow tetragonal variants of arbitrary fineness-an unphysical result, reflecting further mathematical ill-posedness. Retention of the strain gradient energy, (γ αβ = 0) penalizes interfaces of sharply varying strain between tetragonal variants to ensure physically realistic results and unique microstructures from a mathematically well-posed formulation. This is wellunderstood in the literature that studies the formation of martensitic microstructures from non-convex free energy density functions [16][17][18].

Governing equations of non-equilibrium chemistry
The free energy for non-homogeneous composition and strain fields, Eq. [2], must be a minimum at equilibrium. The state of a solid out of equilibrium will evolve to reduce the free energy Π[c, u]. In formulating a kinetic equation for the redistribution of atomic species through diffusion, we are guided by variational extremization of the free energy to identify the chemical potential, µ. Details of this calculation appear as supporting information. The result follows: For solids where c tracks the composition of an interstitial element within a chemically inert host, such as Li in Li c Mn 2 O 4 , µ in Eq. [4] corresponds to the chemical potential of the interstitial element. If c tracks the composition of a substitutional species, such as in alloys or on sublattices of complex compounds (e.g. the cation sublattice of yttria stabilized zirconia), µ is equal to the chemical potential difference between the substitutional species.
The common phenomenological relation for the flux is J = −L(c, e)∇µ, where L is the Onsager transport tensor (see de Groot & Mazur [15]). For an interstitial species, L is related to a mobility [19], while it is a kinetic, intermixing coefficient for a binary substitutional solid [20]. Inserting the flux in a mass conservation equation yields the strong form of the governing partial differential equation (PDE) for time-dependent mass transport. It is of fourth order in space due to the composition gradient dependence of µ in Equation [4]. See Supporting Information for strong and weak forms of this PDE.

Governing equations of mechanical equilibrium: Strain gradient elasticity
Mechanical equilibrium is assumed since elastic wave propagation typically is a much faster process than diffusional relaxation in crystalline solids. Equilibrium is imposed by extremizing the free energy functional with respect to the displacement field. Standard variational techniques lead to the weak and strong forms of strain gradient elasticity. The treatment is technical, for which reason we restrict ourselves to the constitutive relations that are counterparts to the chemical potential equation [4] for chemistry. Coordinate notation is used for transparency of the tensor algebra, and summation is implied over repeated spatial index, I. Details appear in Supporting Information. The final form of the equations is complementary to Toupin's [11], since our derivation is relative to the reference crystal, Ω.
With the deformation gradient F being related to the Green-Lagrange strain as E KL = 1 2 (F iK F iL − δ KL ), the first Piola-Kirchhoff stress tensor and the higher-order stress tensor respectively, are given by The higher-order stress B, which is absent in classical, non-gradient elasticity [21] (and in earlier treatments of mechano-chemistry [4,22]) makes the strong form of gradient elasticity a fourth order, nonlinear PDE in space (Supporting Information). The first three-dimensional solutions to general boundary value problems of Toupin's strain gradient elasticity theory at finite strains were recently presented by the authors [23].

A. Two dimensional examples
We first consider a two-dimensional solid to better visualize the microstructures that can emerge from mechano-chemical spinodal decomposition. Plane strain elasticity is assumed, for which E 13 , E 23 , E 33 = 0, giving e 4 , e 5 = 0, and e 3 = e 1 / √ 2, reducing Equations [1][2][3][4][5][6] to two dimensions. The discussion on twodimensional mechano-chemical spinodal decomposition ( Figure 2) holds: as a particular parameterization of F (c, e 1 , e 2 , e 6 ) we consider a regular solution model as a function of c at zero strain. At c = 0, F (c, e 1 , e 2 , e 6 ) is double-welled in e 2 corresponding to the two rectangular variants. The gradient energy is G (∇c, ∇e 2 ) with a constant, isotropic κ = 0, and constant γ 22 = 0, all other gradient coefficients being zero. While the fullest possible complexity of the coupling is not revealed by these simplifications, the aim here is to present the essential physics that is universal to mechano-chemical spinodal decomposition, postponing details of the more complex couplings and tensorial forms to future communications, where they will be derived for specific materials systems. See Supporting Information for specific forms of F (c, e 1 , e 2 , e 6 ) and G (∇c, ∇e 2 ). Figure 4 shows the evolution of microstructure over a 0.01 × 0.01 domain whose reference (initial) state has the square crystal structure and c = 1. The displacement component u 1 = 10 −5 on the right boundary (X 1 = 0.01), with the remaining boundaries fixed (u = 0). An outward flux is imposed on the top and bottom boundaries causing a decrease in composition, c, starting from the boundaries. As the composition falls, the homogeneous free energy density F loses convexity and the state of the material (c, e 2 ) enters the mechano-chemical spinodal along e 2 = 0 ( Figure 2b). The continuing outward flux first drives material near the top and bottom boundaries fully into the regime where the rectangular crystal structure is stable. As explained for Figure 2b, the negative curvature ∂ 2 F /∂e 2 2 < 0 creates thermodynamic driving forces that distort the square structure, shown in green, into rectangular variants, which form as red/blue laminae (all e 2 values). The coexistence of the parent, square structure with the daughter, rectangular variants also can give rise to cross-hatched microstructures that parallel the tweed microstructures described in the work of Karatha and co-workers [12,13]. Movie S9 in the supporting information shows the formation of such a microstructure.
A laminar, micro-twinned microstructure develops as the two rectangular variants form, distinguished by the sign of e 2 (see legend). Note that e 2 remains continuous because of the penalization of strain gradients ∇e 2 , although discontinuities can develop in ∇e 2 , itself. Because our numerical framework uses basis functions that are continuous up to their first derivatives (see Methods and Supporting Information), ∇e 2 is indeed discontinuous in the computations. The lamination accommodates the strain difference between the two rectangular variants to minimize the free energy. If strains were infinitesimal (|e 1 |, |e 2 |, |e 6 | 0) e 2 would correspond to shear along directions that are rotated π/4 radians from the crystal axes. But the strains in these microstructures can be finite (reflected in the range −0.1 ≤ e 2 ≤ 0.1 in Figure 4) and involve large rotations to accommodate the rectangular variants; this necessitates a finite-strain formulation, as shown in recent studies by Finel et al [41] comparing infinitesimal-strain and finite-strain formulations for polytwinned microstructure evolution in martensitic alloys, and by Hildebrand & Miehe [8].
The micro-twinning and coherency strains are seen in the distorted mesh of discretization (inset of top center panel). The undeformed mesh cells are squares, and hence the numerical discretization strikingly delineates the kinematics of the cubic to rectangular transformation, and highlights the rectangular twins as well as the concentrated distortions along the twin boundaries. Movie S6 in the supporting information shows such mesh distortion and formation of the rectangular twins with a different form of the function F .
In some cases, the long-range nature of elasticity forces like rectangular variants to align even when separated by an untransformed square phase. This is first seen in the finger-like extensions of strain contours from the rectangular variants into the square phase in Figure 4, followed by their alignment and eventual incorporation into laminae of the same variant (top right panel and its evolution shown as an inset). In this instance, although the micro-twins end at an invariant habit plane that is common with the untransformed square phase, the lattice parameters of the rectangular variants differ from those of the square phase, inducing elastic strain energy in the latter. The alignment lowers this strain energy, and also is seen in Figure 5a. Note, however, that this is not a universal feature. Movie S2 shows instances where unlike variants come close to alignment. The fineness of laminae depends on the strain gradient elasticity length scale l e ∼ γ 22 / α,β=1,2,6 (∂ 2 F /∂e α ∂e β ) 2 , as explored in Figure 5. For Figure 5a, the initial and boundary conditions are the same as in the example of Figure  4. Decreasing l e weakens the penalty on the strain gradient ∇e 2 across neighboring, unlike rectangular variants and allows more twin boundaries. Notably, self-similarity is not maintained between microstructures for different l e , even for the same initial and boundary value problem. We understand this to be the influence of elastic strain accommodation: To minimize the total free energy when the strain gradient penalty changes, the physics optimizes twin boundaries via laminations of different sizes as well as different patterns. Importantly however, crystal symmetry admits non-vanishing strain gradient energy coefficients beyond γ 22 = 0 used in these simulations (Supporting Information). Furthermore, the composition dependence of F could be more complex than the simple regular solution model used here. Given the already strong effect of l e alone, we conjecture that varying these forms will have a significant influence on the resulting coherent twin microstructure. The proper form, while guided by crystal symmetry arguments must ultimately be determined experimentally or by first-principles statistical mechanical methods.
The example in Figure 5b further explores the influence of l e . In this suite of computations, the unstrained material with an initially square microstructure, convex free energy density and composition having random fluctuations about c = 0.45 was quenched to a low temperature and into the mechano-chemical spinodal. The local composition and strain evolve under the thermodynamic driving forces detailed in the context of Figure 2. We draw attention, once again, to the changing identity of rectangular variants, their shapes and sizes depending on l e . Also note the progressively finer lamination of rectangular domains with decreasing l e . Such studies suggest how dynamic mechano-chemical spinodal decomposition can lead to an atlas of microstructures, which in turn will determine material properties. Movies S1-S4 in the supporting information show the evolution of some of these microstructures.

B. A three dimensional study
This final example displays the full three-dimensional complexity of microstructures resulting from the mechano-chemical spinodal. We persist with the above simplifications aimed at presenting the essential physics of mechano-chemical spinodal decomposition, and postpone a full exploration of the coupling and anisotropies in coefficients to future work on specific materials systems. The free energy density function used for the threedimensional study appears as supporting information. Figure 6 shows the equilibrium microstructure that results in a solid that is initially in the cubic phase, c = 1, subject to an outflux on all surfaces. The domain is a unit cube with displacement components u 2 , u 3 = 0.01 on the boundary X 1 = 1, with zero displacement, u = 0, on the boundary X 1 = 0. The cubic to tetragonal transformation takes place as c → 0. Under these conditions all three tetragonal variants form, with an intricately interleaved microstructure for strain accommodation. Note the finer microstructures and changing pattern for smaller l e . The inset shows the surface straining around a corner, delineated by the distorted mesh lines. Observe the three, oriented, tetragonal variants formed by twinning deformation from the initial cubic structure. See Movie S7 in the supporting information for a detailed view of the three-dimensional structure of these individual tetragonal variants. Movie S10 in supporting information shows other three-dimensional mi-crostructures whose cross-sectional planes bear closer resemblance to the plate-like structures predicted by twodimensional computations. To the best of our knowledge, such computations of a cubic to tetragonal transformation with twinned variants whose microstructure is controlled by nonlinear, strain gradient interface energies, have not been previously presented.

III. DISCUSSION
Several kinetic mechanisms have been proposed to describe the decomposition of a solid upon entering a twophase region [42,43]. A commonly observed mechanism occurs through localized nucleation followed by diffusional growth that is mediated by the migration of interfaces separating the daughter phases from the parent phase. Nucleation and growth mechanisms are common when the phases participating in the transformation have little or no crystallographic relation to each other. The transformation then proceeds reconstructively, where the interphase interfaces are disordered or at best only semicoherent. Numerical treatments of this mechanism rely on sharp-interface methods such as a level set approach [44].
Other kinetic mechanisms of decomposition involving some degree of coherency are also possible when the participating phases share sufficient crystallographic commonality that order parameters can be defined describing continuous pathways connecting the parent and daughter phases. The structural changes of the crystal that accompany these transformations can fall in one of two categories. One subclass of structural transformations is driven by an internal shuffle, where the atomic arrangement within the unit cell of the crystal undergoes a symmetry breaking change. The unit cell vectors of the crytal may also undergo a change and lower the symmetry of the lattice, but this effect is secondary and in response to the atomic shuffle within the unit cell. The order parameters for the structural change therefore describe the atomic shuffles within the unit cell and are non conserved. The second subclass of structural transformations are driven by a symmetry reducing change of the lattice vectors of the crystal, with any atomic rearrangement of the basis of the crystallographic unit cell occurring in response to the symmetry breaking changes of the unit cell vectors. The natural order parameters describing these transitions are strains. The free energy in the first subclass will exhibit negative curvatures with respect to the nonconserved shuffle order parameters, while the free energy in the second subclass will have negative curvatures with respect to the relevant strain order parameters, as was recently predicted for the cubic to tetragonal transition of ZrH 2 [33].
Decomposition reactions that require both diffusional redistribution and a structural change due to an internal shuffle have been treated successfully and rigorously with phase field approaches that combine Cahn-Hilliard with Allen-Cahn. The Cahn-Hilliard description accounts for the composition degrees of freedom while the Allen-Cahn component describes the evolution of the non-conserved shuffle order parameters required to distinguish the various phases of the transformation [5][6][7][8]. These treatments have included strain energy as a secondary effect serving only as a positive contribution to the overall free energy due to coherency strains. The approach is rigorous if the free energy density remains convex with respect to strain, with instabilities in the free energy appearing as a function of concentration and the non-conserved shuffle order parameters.
Here we have presented a mathematical formulation and an accompanying computational framework for decomposition reactions that combine diffusional redistribution with a structural change driven by a symmetry breaking strain of the crystallographic unit cell as opposed to an internal shuffle within the unit cell. Strains play a primary role, explicitly serving as order parameters to distinguish variants of a daughter phase that has a symmetry subgroup-group relationship with its parent phase due to a structural change of the crystallographic unit cell. Crucial to the description is that the driving force for the formation of the daughter from a supersaturated parent phase emerges from an instability with respect to composition. The treatment therefore combines Cahn-Hilliard for composition with a description of martensitic transformations at fixed concentration introduced by Barsch and Krumhansl [9]. The existence of simultaneous instabilities with respect to strain and composition has not been considered as a mechanism to describe decomposition reactions upon quenching into a two-phase region.
The possibility of such a mechano-chemical spinodal mechanism is motivated by recent first-principles studies of the cubic to tetragonal phase transformations of ZrH 2 , where the free energy of the high temperature cubic form was predicted to become unstable with respect to strain upon cooling below the cubic to tetragonal second order transition temperature [33]. The ZrH 2−2c hydride can accommodate large concentrations of hydrogen vacancies, c, and has a phase diagram that is topologically identical to that depicted in Figure 2c, with a two-phase region separating a hydrogen-rich tetragonal form of ZrH 2−2cα from a cubic form of ZrH 2−2c β (with c α < c β ). See Zuzek et al. [32]. [45] To be consistent with the predicted free energies for stoichiometric ZrH 2 (i.e. c=0) and the experimental T versus c phase diagram with the form of Figure 2c, the free energy of this hydride as a function of composition and strain (i.e. e 2 and e 3 ) should be similar to those depicted in Figures 2a and b, and Figure  3. Decomposition upon quenching cubic ZrH 2−2c into the two-phase region can therefore proceed through a mechano-chemical spinodal decomposition mechanism.
Not only does the phenomenological description introduced in this work describe a new mechanism of decomposition that is qualitatively distinct from previous combined Cahn-Hilliard and Allen-Cahn approaches, its numerical solution also proves to be substantially more complex due to contributions from strain gradient terms. Indeed, even numerical solutions to general, three-dimensional boundary value problems of gradient elasticity at finite strain were not available until presented recently by the authors [23]. Furthermore, as other authors have pointed out before, the use of strain metrics as order parameters makes a reliance on infinitesimal strains untenable due to the rigid rotations that accompany the finite strains characterizing most structural transformations [8,41]. The use of finite strain metrics introduces geometric non-linearity into the problem, which while having been treated in past Cahn-Hilliard and Allen-Cahn approaches [8], adds additional numerical challenges when also considering strain gradient contributions. These challenges have been overcome in this work as demonstrated by our three-dimensional numerical examples.
Mechano-chemical spinodal decomposition as described here can be expected in solids forming high temperature phases that exhibit dynamical instabilities at low temperature. An accumulating body of first-principles calculations of Born-Oppenheimer surfaces have shown that many high temperature phases are dynamically unstable at low temperature [26,27], becoming stable at high temperature through large anharmonic vibrational excitations [28][29][30][31]33], usually through a second order phase transition. While such instabilities are frequently dominated by phonon modes describing an internal shuffle, a subset of chemistries become dynamically unstable at low temperature with respect to phonon modes that break the symmetry of the crystal unit cell [30,33], an instability that can be described phenomenologically with strain order parameters. If, upon alloying such compounds, the high symmetry phase becomes stable by passing through a two-phase region, its free energy surface will resemble that of Figures 2b and 3, and thereby make the solid susceptible to the mechanochemical spinodal decomposition mechanism described here.
While first-principles and experimental evidence suggests mechano-chemical spinodal decomposition should occur in ZrH 2−2c , we expect it to occur in a wide range of other chemistries as well. One possible example, as described in the introduction, is the decomposition of cubic yttria stabilized zirconia [34,35] upon quenching from the high temperature cubic phase into a two-phase region separating a low-Y tetragonal phase from a Y-rich cubic phase. Past treatments of this transformation [5][6][7] relied on non-conserved order parameters to distinguish the different tetragonal variants from each other and from the cubic parent phase. The elastic part of the free energy density function was parameterized using linearized elasticity and infinitesimal strains, as other authors have also done [39,40]. The chemical part of the free energy density was assumed to have a negative curvature as a function of the non-conserved order parameter at Zr-rich compositions, but made up of convex (and quadratic) potentials with respect to strain. We reiterate that such a treatment rests on the implicit assumption that internal shuffles drive the cubic to tetragonal transformation, while the free energy as a function of strain remains convex for all relevant values of the non-conserved order parameters.
The possibility also exists that a coarse-grained free energy density of cubic ZrO 2 exhibits negative curvatures with respect to strain below the cubic to tetragonal transition temperature making the formalism developed here an accurate description of decomposition reactions of yttria stabilized zirconia. While the precise mechanism and nature of the cubic to tetragonal transition of pure ZrO 2 remains to be resolved [29], first-principles calculations predict that the cubic form of ZrO 2 is dynamically unstable with respect to transformation to the tetragonal variant [36]. A rigorous statistical mechanics treatment [28,33] is required to determine whether the cubic to tetragonal transition of pure ZrO 2 is accompanied by a change in the sign of the curvature of the free energy density with respect to e 2 and e 3 . If this proves to be the case, yttria stabilized zirconia should also be susceptible to mechano-chemical spinodal decomposition upon quenching, consistent with the coherent spinodal microstructures between tetragonal and cubic phases observed in single crystal regions of quenched cubic Zr 1−c Y x O 2−c/2 [35].
A mechano-chemical spinodal may also play a role in a variety of important electrode materials for Liion batteries, and intercalation compounds considered for two-dimensional nano-electronics. These include cubic LiMn 2 O 4 transforming to tetragonal Li 2 Mn 2 O 4 [37]. While most mechano-chemical spinodal transitions will occur in three dimensions, many of the qualitative features of these transitions are more conveniently illustrated in two dimensions. Our two-dimensional studies should also prove relevant to understanding mechano-chemical phase transformations in two-dimensional layered materials for nano electronics. Materials such as TaS 2 , are susceptible to Peierls instabilities upon variation of the composition of adsorbed or intercalated guest species that donate to or extract electrons from the sheetlike host [38]. Our phenomenological treatment by introduction of the concept of a mechano-chemical spinodal, coupled with gradient stabilization of the ensuing nonconvex free energy in strain-composition space, thus offers a framework with potential for extension to a wide range of phase transformation phenomena.
With regard to the fundamental thermodynamic underpinnings of analogous processes, the phenomenological model introduced here can also be used to describe temperature driven martensitic transformations. The composition, c, would then be analogous to the internal energy density, U , while the chemical potential, µ, would be analogous to the temperature, T . Due to the presence of temperature gradients throughout the solid, though, the starting point must be a formulation of the entropy density, S , (as opposed to the Helmholtz free energy) as a functional of spatially varying U and displacements, u. In analogy with Equation (2), the total entropy can be written as a volume integral over a homogeneous entropy density that depends on the local internal energy density and strainS (U , e), as well as a non-uniform entropy density contribution expressed in terms of gradients in internal energy density and strain (∇U , ∇e). Variational maximization of the entropy will yield mechanical equilibrium equations as well as an expression for the temperature T (strictly speaking, for its reciprocal, 1/T ), similar to Equation (4) for the chemical potential. Due to the presence of gradient terms, ∇U and ∇e, the temperature will not only be a function of the local internal energy density and strain, but also gradients in these field variables. As with the diffusion problem treated here, heat flow can be described with an Onsager expression relating the heat flux to a gradient in temperature. The possibility exists that the homogeneous entropyS exhibits instabilities with respect to internal energy density U , allowing for spinodal decomposition with respect to the redistribution of U in a manner that is similar to the well understood problem of spinodal decomposition with respect to composition. The treatment introduced here can therefore describe (upon replacing c with U and µ with T ) temperature driven martensitic transformations for solids that exhibit instabilities with respect to both internal energy density and strain. Past treatments of this phenomenon also relied on a Barsch and Krumhansl approach, solving the mechanical problem together with Fourier's law of heat conduction [39,40]. However, these studies were restricted to two dimensions [39][40][41] with some neglecting geometric non-linearity by relying on infinitesimal strains [39,40]. At a more fundamental level, they did not consider the possibility of spinodal instabilities with respect to the redistribution of internal energy density.

IV. METHODS
The governing fourth-order PDEs of mass transport and gradient elasticity are solved numerically in weak form, wherein second-order spatial gradients appear on the trial solutions and their variations. Numerical solutions require basis functions that are continuous up to their first spatial gradients, at least. Our numerical framework (see Rudraraju et al. [23]) uses isogeometric analysis [24] and the PetIGA code framework [25] with spline basis functions, which can be constructed for arbitrary degree of continuity (Supporting Information). This framework has proven pivotal to the current work. Consider the experimentally derived equilibrium phase diagram for transition metal hydrides [1] shown schematically in Figure 2c of the main text. The single cubic phase, stable for all composition and strains at high temperature, only remains stable at high composition when quenched to low temperature. At low composition, the tetragonal phase is stable at low temperature. Also recall that first principle calculations have shown that the tetragonal phase is unstable with respect to strain and decomposes into three energetically equivalent variants [2]. The homogeneous free energy density therefore is symmetric with respect to zero strain. The two-phase coexistence region implies that the free energy density surface, F in strain-composition, (e, c) space admits a common tangent hyperplane. For clarity, let us consider the zero strain hyperplane in this space. Its projection appears in Figure S1. For a point in the solid with composition c * , the common tangent construction leads to compositions c α and c β , respectively of the tetragonal and cubic phases, with Furthermore, the stability of the cubic and tetragonal phases with respect to composition imply that Recognizing that c β > c α in Figure 2c of the main text, it follows from Equation (2) that there exists c + α > c α and c − β < c β such that But then, from Equations (1) and (3), it follows that, for a smooth F that is at least twice differentiable in the interval (c α , c β ), there must exist at least one value c γ with c + α ≤ c γ ≤ c − β satisfying Therefore there is at least one intermediate value of composition in the two-phase coexistence region at which the free energy density is non-convex with respect to composition. This result is simply a consequence of the mean value theorem of calculus.
FIG. S1: Figure S1: Projection of the free energy on to the zero strain hyperplane.

II. FINITE STRAIN KINEMATICS
In the interest of clarity and consistency with the main text, the Euclidean basis of our Cartesian coordinate system is denoted by X i , i = 1, . . . 3, X i · X j = δ ij . This constitutes a mild abuse of notation, because reference positions have also been denoted by X 1 , X 2 , X 3 . The high symmetry, cubic reference configuration of the crystal is denoted by Ω, on which material points are labelled by their position vectors X = X 1 X 1 + X 2 X 2 + X 3 X 3 . The deformation map from Ω to the current configuration of the crystal is ϕ(X, t) = X + u, where the displacement field u was introduced in the main text. The deformation gradient tensor is F = ∂ϕ/∂X = 1 + ∂u/∂X. Using lower case indices to denote the components of vectors and tensors on the current configuration, and upper case indices for those on the reference configuration, we have the coordinate expression The Green-Lagrange strain tensor is E = 1 2 (F T F −1), where 1 is the second-order isotropic tensor. We repeat this relation in coordinate notation from the main text: E IJ = 1 2 (F iI F iJ − δ IJ ). The polar decomposition theorem states that the deformation gradient can be decomposed multiplicatively as F = RU , where R is an (orthogonal) rotation tensor belonging to the space SO(3). It therefore satisfies R T R = 1. The other tensor in this decomposition, U , is symmetric and positivedefinite. From the polar decomposition it follows that E as defined above is invariant to rotation of the current configuration. The requirement of material frame invariance dictates that the dependence of the free energy density on elasticity must be as a function of E. Such a free energy density function is invariant to rotations on the current configuration.

III. THE CHEMICAL POTENTIAL: VARIATIONAL DERIVATION
Since the crystal is, in general, far from chemical equilibrium we follow the standard treatment of nonequilibrium thermodynamics to model the chemistry. In order to obtain a relation for the chemical potential we consider variations on the composition field of the form c ε := c + ε ξ, where ξ is a function such that, on the Dirichlet boundary, Γ c ⊂ Γ, we have c =c and ξ = 0.
Carrying out the differentiation with respect to ε and integration by parts yields where n is the outward normal vector to Γ. At chemical equilibrium we have δΠ/δc = 0, yielding two conditions: Standard variational arguments lead to the conclusion that, for chemical equilibrium, the quantities within the brackets must vanish on the corresponding domains Ω and Γ. The expression contained within brackets in the integrand of Equation [7] is the chemical potential µ. In the non-equilibrium setting that is of interest here, we continue to be guided by the above treatment and use the same definition for the chemical potential [3]: Variants of Equation [8] appear in a proper variational derivation of chemical equilibrium-a detail that is sometimes overlooked in phase field treatments. We will impose chemical equilibrium of the boundary Γ, by requiring the term in brackets to vanish. With the standard, phenomenological relation J = −L(c, e)∇µ for the flux, the strong form of the governing PDE of mass balance follows: Note that we have assumed that there is no Dirichlet boundary; i.e. Γ c = ∅, in stating the final strong form of this PDE. The conjugate boundary condition on the mass flux holds on the entire boundary Γ, in Equation [10 2 ]. Additionally, we have the boundary condition arising from requiring chemical equilibrium on Γ, in Equation [10 3 ]. Using standard variational arguments, the weak form of mass balance can be obtained by starting from [10]. Finally, we have the initial condition on composition, Equation [10 4 ]. The substitution of Equation [9] in J = −L∇µ, followed by this phenomenological flux relation's substitution in (10 1 ) shows that this is a fourth order PDE in space.

A. On mass transport posed in the reference configuration
We note that the mathematical formulation of diffusion is commonly stated in the current (deformed) configuration, which is the true state of the solid. However, for crystalline solids (with negligible extended defects) it is more convenient to define Onsager transport coefficients and evaluate composition gradients on the reference (undeformed) configuration, here denoted by Ω. Diffusion in the crystalline solids results from atomic hops between well-defined crystal sites, which can always be mapped onto a reference crystal structure. Increasingly, mobility coefficients in multi-component solids are calculated from first principles through the evaluation of Kubo-Green expressions in kinetic Monte Carlo simulations on the reference crystal in configuration Ω. Since the mobility tensor is reported on Ω, it is most convenient to adopt the Lagrangian description for diffusion over this configuration in single crystalline solids.
Our numerical implementation is based on the weak form, whose derivation we begin by first reintroducing ξ, now in the guise of a weighting function lying in the function space Q = {ξ ∈ H 1 (Ω)}. We seek functions c ∈ P = {c ∈ H 1 (Ω)} such that Equation [11] is obtained by multiplying ξ into [10 1 ], integrating by parts twice, and using the Neumann boundary condition [10 2 ]. Note that the higher-order Dirichlet boundary condition (10 3 ) has not been built into the function spaces as is commonly done in weak forms. Instead, it is imposed weakly using the method of Nitsche [4][5][6] and extending the weak form to Equation [12] is based on the coupling tensor θ α = 0 for all α, and for κ, γ αβ independent of c, corresponding to the numerical examples presented in the main text. This weak form maintains consistency. Note that the second-to-last term that has been added over the form in Equation [11] is symmetric with the term preceding it. This formulation shows better numerical performance than the also-consistent formulation in which this additional term is preceded by a negative sign. Finally, τ is a stabilization parameter that is inversely proportional to the element size.

IV. MECHANICAL EQUILIBRIUM
As explained in the main text, since elastic wave propagation is orders of magnitude faster than the kinetics of mass transport, it can be assumed that mechanical equilibrium is attained instantaneously for the prevailing composition field. To obtain the relevant equilibrium equations, we introduce variations on the displacement field: u ε = u + εw, such that on the Dirichlet boundary Γ ui , the fields satisfy u i =ū i and w i = 0. We construct the first variation of Π[c, u] with respect to u.
where e ε are obtained by first replacing u by u ε in the kinematic relations for E discussed above, which are then substituted in the relations for the reparameterized strains e. We recall the constitutive relations for the first Piola-Kirchhoff stress and the higher-order stress, respectively: Using these constitutive relations, Equation [13] can be manipulated to yield the Euler-Lagrange equation, which is also the weak form of mechanical equilibrium for a strain gradient elastic material, on Ω in terms of P and B. We retain coordinate notation in the interest of transparency: We introduce the normal and surface gradient operators, ∇ n and ∇ s , defined by ∇ n ψ = ∇ψ · n, and ∇ s ψ = ∇ψ − (∇ n ψ)n.
Starting from Equation [16] and applying integration by parts several times we have the strong form of mechanical equilibrium for a strain gradient elastic material. See Rudraraju et al. [7] for detailed steps. The results of this variational treatment are analogous to those of Toupin [8], except that it is carried out over the reference configuration of the crystal, rather than its current configuration.
Here, Γ = Γ ui ∪ Γ Ti for i = 1, 2, 3, is the decomposition of the boundary surface into a subset with Dirichlet boundary condition on displacement component u i and Neumann boundary condition on traction component T i , respectively. Finally, Υ is a smooth boundary edge on which a jump in higher-order stress traction may arise. The fourth-order nature of the governing partial differential equation emerges on evaluating [14][15] and substituting in [18 1 ]. The Dirichlet boundary condition in [18 2 ] has the same form as for conventional elasticity. However, its dual Neumann boundary condition, [18 3 ] is notably more complex than its conventional counterpart, which would have only the first term on the left hand-side. Equation [18 4 ] is the higher-order Neumann boundary condition on the higher-order stress traction. The physical interpretation of this boundary condition in its homogeneous form is that the boundary has no mechanism to impose a moment on bonds at the atomic scale. Equation [18 5 ] specifies that there is no discontinuity of higher-order stress traction across edges. The weak form of mechanical equilibrium has already been encountered in [16]. For completeness we specify functional spaces: Find u lying in S = {u ∈ H 2 (Ω 0 )| u i =ū i on Γ ui } such that given functions T i on Γ Ti for i = 1, 2, 3 and tensors P , B defined by Equations [14][15], Equation [16] is satisfied for all w lying in V = {w ∈ H 2 (Ω)| w i = 0 on Γ ui }.

V. ISOGEOMETRIC ANALYSIS FOR HIGH-ORDER PDES
Isogeomeric Analysis (IGA) is a mesh-based numerical method with NURBS (Non-Uniform Rational B-Splines) basis functions. The NURBS basis functions have many desirable properties. They are partitions of unity with compact support, satisfy affine covariance and provide certain advantages over Lagrange polynomial basis functions, which are the mainstay of Galerkin finite element methods. These advantages include the imposition of higher-order, C n -continuity, positivity, convex hull properties, and being variation diminishing. For a detailed discussion of the NURBS basis and IGA, interested readers are referred to Cottrell et al. [9]. However, we briefly present the construction of the basis functions.
The building blocks of the NURBS are univariate Bspline functions that are defined as follows: Consider positive integers p and n, and a non-decreasing sequence of values χ = [ξ 1 , ξ 2 , ...., ξ n+p+1 ], where p is the polynomial order, n is the number of basis functions, the ξ i are coordinates in the parametric space referred to as knots (equivalent to nodes in FEM) and χ is the knot vector. The B-spline basis functions B i,p (ξ) are defined starting with the zeroth order basis functions and using the Cox-de Boor recursive formula for p ≥ 1 [10] The knot vector divides the parametric space into intervals referred to as knot spans (equivalent to elements in FEM). A B-spline basis function is C ∞ -continuous inside knot spans and C p−1 -continuous at the knots. If an interior knot value repeats, it is referred to as a multiple knot. At a knot of multiplicity k, the continuity is C p−k . The boundary value problems in the main text consider only simple geometries, for which B-spline basis functions have sufficed. However, we outline the extension to NURBS basis functions for the sake of completeness, noting that the numerical formulation as presented is valid for any single-patch NURBS geometry. Using a quadratic B-spline basis ( Figure S2), a C 1 -continuous one dimensional NURBS basis can be constructed.
where w i are the weights associated with each of the B-spline functions. In higher-dimensions, NURBS basis functions are constructed as a tensor product of the one dimensional basis functions:

VI. NUMERICAL EXAMPLES
Recall that in two dimensions, the free energy density function depends only on c, e 1 , e 2 and e 6 . In this first communication we wish to focus only on broad phenomenology, for which reason we choose the simplest versions of the free energy density function. Apart from the dependence on c, which is essential for the mechanochemical spinodal, the coefficients are constant. We have introduced the coefficients d c and d e , which control the depths of the free energy wells and s e , which determines the strain minimum in terms of e 2 . Figure 2 of the paper depicts this construction. We also reduce κ and γ αβ to isotropic forms: κ = κ1, γ αβ IJ = λ e l 2 e δ αβ δ IJ , and set the composition gradient-strain gradient coupling coefficients θ α = 0. While the fullest possible complexity of the coupling is not revealed by these simplifications, the aim here is to present the essential physics that is universal to mechano-chemical spinodal decomposition, postponing details of the more complex couplings and tensorial forms to future communications, where they will be derived for specific materials systems. Accordingly, in two dimensions, we have: G (∇c, ∇e 2 ) = 1 2 ∇c · κ∇c + 1 2 ∇e 2 · λ e l 2 e ∇e 2 (24) where d c = 2.0, d e = 0.1, s e = 0.1, κ = 10 −6 and λ e = 10 −7 . The problem domain is a square of dimensions 0.01 × 0.01.
In keeping with the above-explained simplifications aimed at presenting the essential physics of mechanochemical spinodal decomposition, the free energy density G (c, e, ∇c, ∇e) = 1 2 ∇c · κ∇c + 1 2 ∇e 2 · λ e l 2 e ∇e 2 + 1 2 ∇e 3 · λ e l 2 e ∇e 3 (26) where d c = 2.0, d e = 100.0, s e = 0.1, κ = 10 −2 and λ e = 1.0. The problem domain is a unit cube. The strain gradient energy coefficient tensors have been chosen to be isotropic, γ αβ IJ = λ e l 2 e δ αβ δ IJ , and the three-dimensional formulation of the main text has been simplified here by setting the composition gradient-strain gradient coupling coefficients θ α IJ = 0. surface, F (c, e 1 , e 2 , e 6 ) for the two-dimensional formulation, restricted to the {c, e 2 } sub-space. The microstructure that forms is similar to that in Movie S1.