Mechano-chemical spinodal decomposition: A phenomenological theory of phase transformations in multi-component, crystalline solids
MMechano-chemical spinodal decomposition: A phenomenological theory of phasetransformations in multi-component, crystalline solids
Shiva Rudraraju , Anton Van der Ven & Krishna Garikipati
We present a phenomenological treatment of diffusion-driven martensitic phase transformationsin multi-component crystalline solids that arise from non-convex free energies in mechanical andchemical variables. The treatment describes diffusional phase transformations that are accompa-nied by symmetry breaking structural changes of the crystal unit cell and reveals the importanceof a mechano-chemical spinodal , defined as the region in strain-composition space where the freeenergy density function is non-convex. The approach is relevant to phase transformations whereinthe 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 decom-position are variationally derived from a free energy density function that accounts for interfacialenergy via gradients of the rapidly varying strain and composition fields. A robust computationalframework for treating the coupled, higher order diffusion and nonlinear strain gradient elasticityproblems is presented. Because the local strains in an inhomogeneous, transforming microstruc-ture can be finite, the elasticity problem must account for geometric nonlinearity. An evaluation ofavailable experimental phase diagrams and first-principles free energies suggests mechano-chemicalspinodal decomposition should occur in metal hydrides such as ZrH − c . The rich physics thatensues is explored in several numerical examples in two and three dimensions and the relevance ofthe mechanism is discussed in the context of important electrode materials for Li-ion batteries andhigh temperature ceramics. PACS numbers:Keywords:
I. INTRODUCTION AND FORMULATION
Spinodal decomposition is a continuous phase trans-formation mechanism occuring throughout a solid thatis far enough from equilibrium for its free energy den-sity to lose convexity with respect to an internal degreeof freedom. The latter could include the local composi-tion as in classical spinodal decomposition described byCahn and Hilliard [1], or a suitable non-conserved orderparameter as in Allen and Cahn’s theory for spinodalordering [2]. A key requirement for continuous transfor-mations is that order parameters can be formulated touniquely describe continuous paths connecting the var-ious phases of the transformation. These phases thencorrespond to local minima on a single, continuous freeenergy density surface in that order parameter space. Forclassical 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 densitysurface for all phases participating in a transformationimplies, by geometric necessity, the presence of domainsin order-parameter space where the free energy density isnon-convex. Reaching those domains through supersat-uration (by externally varying temperature or composi-tion) makes the solid susceptible to a generalized spinodaldecomposition.Many important multi-component solids undergophase transformations that couple diffusional redistribu-tion of their components with a structural change of thecrystallographic unit cell. One prominent example is thedecomposition that occurs when cubic yttria-stabilized zirconia Zr − c Y c O − c/ is quenched into a two-phaseequilibrium region between tetragonal Zr − c Y c O − c/ having low Y composition and cubic Zr − c Y c O − c/ hav-ing a high Y composition. Another occurs in Li batteryelectrodes made of spinel Li c Mn O . Discharging to lowvoltages causes the compound to transform from cubicLiMn O to tetragonal Li Mn O through a two-phasediffusional phase transformation mechanism. As withsimple diffusional phase transformations, these coupleddiffusional/martensitic phase transformations can eitheroccur through a nucleation and growth mechanism, or, ifcertain symmetry requirements are met, through a con-tinuous mechanism due to an onset of an instability withrespect to composition and/or a structural order param-eter.Here we present a treatment of coupled diffu-sional/martensitic phase transformations triggered byinstabilities with respect to both strain and composi-tion. These phase transformations are characterized bya mechano-chemical spinodal that is defined as a non-convex region of the free energy density function in thestrain-composition space. The possibility of a mechano-chemical spinodal decomposition is motivated by recentfirst-principles studies of martensitic transformations intransition metal hydrides where a high temperature cu-bic phase is predicted to display negative curvatures withrespect to strain, thus making strain a natural order pa-rameter to distinguish the cubic parent phase from itslow temperature tetragonal daughter phases [33]. In ad-dition, the coupling with composition degrees of freedom(e.g. through the introduction of hydrogen vacancies inthe metal hydrides) allows for the possibility that the free a r X i v : . [ phy s i c s . c h e m - ph ] S e p energy also exhibits a negative curvature with respect tothe composition. Mechano-chemical spinodal decomposi-tion, therefore, is a phenomenon that is likely present inmany multi-component materials but has to date beenoverlooked as a mechanism by which a high symmetryphase can decompose martensitically and through dif-fusional redistribution upon quenching into a two phaseregime. Structural phase transformations in solids drivenby instability with respect to an internal shuffle of theatoms within the unit cell have been treated rigorouslyin the literature with coupled Cahn-Hilliard and Allen-Cahn approaches [5, 6, 7, 8]; but mechano-chemical spin-odal decomposition is fundamentally different and neces-sitates a coupled treatment of both the strain and com-position instabilities.Our treatment is based on a generalized, Landau-type free energy density function that couples strainand composition instabilities. The governing equationsof mechano-chemical spinodal decomposition, obtainedby variational principles, generalize the classical equa-tions of the Cahn-Hilliard formulation [1], and of nonlin-ear gradient elasticity [9, 11] by coupling these systems.The ability to solve this complex, nonlinear, strain- andcomposition-gradient-driven, mechano-chemical systemfor sufficiently general initial and boundary value prob-lems in two and three dimensions also has been lackingheretofore. We introduce the computational frameworkto obtain such solutions in general, three-dimensionalsolids. This newfound capability allows us to then rein-force our discussion of the phenomenology with dynamicspredicted by the numerical solutions. A. The mechano-chemical spinodal in twodimensions
For accessibility of the arguments, we first considerthe two-dimensional analogue of the cubic to tetrago-nal transformation: the square to rectangle transforma-tion. The high-symmetry square lattice will serve as thereference crystal relative to which strains are measured.Lower symmetry lattices that can be derived from thesquare lattice by homogeneous strain include the rect-angle, the diamond and lattices where there are no con-straints on the cell lengths and their angles.We use the composition c , varying between 0 and1, as our order parameter for the chemistry of our bi-nary two-dimensional solid. Symmetry breaking struc-tural changes are naturally described by the strain rel-ative to the high symmetry square lattice. The elasticfree energy density is also a function of strain. In bothcontexts, 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 = (Supporting Informa-tion). In two dimensions, the relevant strain componentsare E , E and E = E . However, it is more conve- FIG. 1: (a) Square to rectangle transition (2D) inreparametrized strain space. (b) Cubic to tetragonal transi-tion (3D) in reparametrized strain space (Equation[1]). Lat-tice vectors are labelled by the corresponding coordinatedirections X , X , X and the corresponding lower symmetryphases are labelled 1, 2, 3. The deformations shown in thesub-figures are area/volume preserving, respectively. nient to use linear combinations of these components thattransform according to the irreducible representations ofthe point group of the high symmetry square lattice. Intwo dimensions, these include e = ( E + E ) / √ e = ( E − E ) / √ e = √ E . Here, e and e reduce to the dilatation and shear strain, respectivelyin the infinitesimal strain limit. The strain measure e uniquely maps the square lattice into the two rectangu-lar variants [Figure 1a]: positive and negative e generatethe rectangles elongated along the global X and X di-rections, respectively. The equivalence of the rectangularvariants under the point group symmetry of the squarelattice ( e = 0) restricts the free energy density to evenfunctions of e .If the crystalline solid has multiple chemical species,its free energy density dependence on e , e and e canchange with composition, c . Figure 2a illustrates a freeenergy density surface, F ( c, e ), for a binary solid that,at higher temperature, forms a solid solution havingsquare symmetry. In this case F is convex, a condi-tion made precise by specifying positive eigenvalues ofits Hessian matrix over the { c, e } space. Additionally, F ( c,
0) is a minimum with respect to e for fixed c , mak-ing the square phase stable for all c at this tempera-ture. However, at a lower temperature, F may lose con-vexity, inducing the notion of a mechano-chemical spin-odal region . We define it as the domain in { c, e } spaceover which the Hessian matrix admits negative eigen-values, as illustrated in Figure 2b. Here, we focus onthe conditions ∂ F /∂c < ∂ F /∂e <
0. Thesquare phase ( e = 0) remains stable at high compo-sition ( c ∼
1) with positive eigenvalues of the Hessian(Figure 2b). But it is mechanically unstable at low com-position, with ∂ F /∂e < c, e ) ∼ (0 , ∂ F /∂e > c, e ) = (0 , ± s e ). While not shownin Figure 2b we assume that F is convex with respect to e and e . Figure 2c illustrates a schematic temperature-composition phase diagram consistent with the free en-ergy densities of Figure 2a-b. The square/cubic phase β ,forms at high composition or at high temperature; therectangle/tetragonal phase α , forms at low compositionand low temperature. A large two-phase region separatesthem.Zuzek et al. [32] have obtained phase diagrams ex-perimentally that are topologically equivalent to Figure2c for the ZrH − c system. Recent first principles cal-culations [33] have demonstrated the existence of a me-chanical instability that exists in this system at low c vianon-convexity with respect to strain. Furthermore, thetwo-phase coexistence at low temperature also impliesnon-convexity with respect to composition as we havedemonstrated in supporting information. Figures 2a andb therefore represent a two-dimensional analogue of thefree energy for such systems. FIG. 2: (a) Free energy density for the 2D formulation athigh temperature, and (b) at low temperature, showing themechano-chemical spinodal along with two energy minimizingpaths (brown lines) and their corresponding minimum energystrained structures. (c) Temperature-composition phase dia-gram.
Consider our model binary solid, annealed at high tem-perature T h , to form a solid solution in the square phase, β . Its state is at point A in Figure 2a, with e = 0. Itis then quenched into the two-phase region (Figure 2c)with free energy density at point B in Figure 2b. Fora quench at sufficiently high rate, the dimensions of thesquare lattice, controlled by strains e , e and e , andthe composition remain momentarily unchanged. How-ever, since the state at point B satisfies ∂ F /∂e < ∂ F /∂c <
0, there exist thermodynamic drivingforces for segregation by strain and composition withinthe mechano-chemical spinodal.Diffusion being substantially slower than elastic relax-ation, the solid immediately deforms to either positive or negative e , driven towards a local minimum at con-stant c . These deformations due to the mechanical insta-bility will happen like many martensitic transformationswhere a mix of symmetrically equivalent rectangular vari-ants coexist to minimize macroscopic strain energy. Forfinite but moderate strain, the transformation could pro-ceed coherently, even if the two symmetrically equivalentrectangular variants coexist. We neglect non-essentialcomplexities of this process and assume that, instantlyupon quenching, finitely sized neighborhoods of the soliddeform homogeneously into one of these rectangular vari-ants at the original composition.The solid also becomes susceptible to uphill diffusionbecause ∂ F /∂c < e , sincethe valleys traversing the local minima, ∂ F /∂e = 0,between the square lattice at c = 1 and the rectangularlattices at c = 0 span intervals of negative and positive e . Mechano-chemical spinodal decomposition sets in.Composition modulations are amplified: high c regionsstrive to be more square (point C) while low c regionsstrive to be more rectangular (points D or D ). How-ever, since coherency is maintained, some neighborhoodsin the solid will be frustrated from attaining strains thatensure minima, ∂ F /∂e = 0, for local values of c . Co-herency strain-induced free energy penalties arise to alterthe driving forces for purely chemical spinodal decompo-sition. Movie S5 in the supporting information shows theevolution of the state ( c, e ) of the material points on thefree energy manifold F . B. Mathematical formulation: three dimensions
Armed with the insight conveyed by the two-dimensional study, we next lay out the three-dimensionaltreatment, in which setting the considered phase trans-formations proceed via lattice deformation and diffusion.The arguments are made more concrete by consideringthe general mathematical form of the free energy den-sity.
1. Strain order parameters
The strain measures e - e are first redefined using thefull Green-Lagrange strain tensor components in threedimensions: e = 1 √ E + E + E ) , e = 1 √ E − E ) ,e = 1 √ E + E − E ) , e = √ E = √ E ,e = √ E = √ E , e = √ E = √ E (1)In the limit of infinitesimal strains, e describes the di-latation, while e , e and e reduce to shears. The pointgroup operations of the cubic crystal leave e invariantsince it is the trace of E , while collecting each of the sub-sets { e , e } and { e , e , e } into a symmetry-invariantsubspace whose elements transform into each other. Themeasures e and e are especially suited as order parame-ters to describe cubic to tetragonal distortions. All threetetragonal variants that emerge from the cubic referencecrystal can be uniquely represented by these two mea-sures. See Figure 1b where tetragonal distortions of thecubic crystal along the X , X and X axes have beenlabelled as 1, 2 and 3, respectively. Deviations from thedashed lines in the e - e space of Figure 1b correspondto orthorhombic distortions of the cubic reference crys-tal. This role of e and e as structural order parametersto denote the degree of tetragonality, and to distinguishbetween the three tetragonal variants, complements theirfundamental purpose as arguments of the elastic free en-ergy density. FIG. 3: Mechano-chemical spinodal for the 3D formulationdepicted by contours of the free energy manifold along the e - e and e - c planes. The three energy minimizing paths (brownlines) and their corresponding energy minimizing strainedstructures and are also shown.
2. The mechano-chemical spinodal in three dimensions
For brevity we write the strains as e = { e , . . . , e } .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 compositionand strain dependence cannot be additively separated.Figure 3, for example, illustrates contour plots of F ( c, e )on the e - e and e - c planes for a binary solid havinga temperature-composition phase diagram similar tothat of Figure 2c. To generate the tetragonal α phaseas illustrated in that diagram, the homogeneous freeenergy density as a function of strain must qualitativelyfollow the contours of the e - e plane at c = 0 in Figure3. Here, the tetragonal variants are local minimizers of F in e - e space, with equal minima. In turn, to obtain the cubic β phase, F ( c, e ) must follow the contours ofthe e - e plane at c = 1. Thus, F changes smoothlyfrom convex with respect to e and e on the c = 1 planeto non-convex at c = 0 to form the three variants of thetetragonal α phase. Importantly, the planes at c = 0and c = 1 themselves must be minimum energy surfacesto represent the tetragonal α and cubic β phases,respectively. Movie S8 in the supporting informationshows the evolution of the state ( c, e , e ) of the materialpoints on the free energy manifold F . Gradient regularization of the free energy density:
Fol-lowing van der Waals [10], and Cahn & Hilliard [1] intheir treatment of non-uniform composition fields, we canextend the total free energy density beyond F by writ-ing it as a Taylor series, retaining terms that depend onthe composition gradient, ∇ c . We extend this gradientdependence to the strain measures ∇ e , as did Barsch &Krumhansl [9] following Toupin [11]. (Also see Karathaand co-workers [12, 13], for treatments using infinitesi-mal and finite strain, respectively.) These gradient de-pendences appear in a non-uniform free energy density G ( c, e , ∇ c, ∇ e ). Frame invariance of F and G is guar-anteed since the members of e are linear combinations ofthe tensor components of E . They also must be invari-ant under point group operations of the cubic referencecrystal.
3. Free energy functional
The crystal occupies a reference (undeformed) config-uration Ω with boundary Γ. The total free energy, Π, isthe integral of the free energy density F + G over thesolid with boundary contributions included. Thus, Π isa functional of the composition c and the displacementvector field u , from which the strains are derived (seeSupporting Information):Π[ c, u ] = Z Ω ( F + G ) d V − X i =1 Z Γ Ti u i T i d S. (2)where traction vector component T i is specified on theboundary subset Γ T i ⊂ Γ. Following the above authorswe include only quadratic terms in the gradients, but ina generalization we also allow cross terms between ∇ c and ∇ e in the non-uniform contribution G , which cantherefore be written as G ( c, e , ∇ c, ∇ e ) = 12 ∇ c · κ ( c, e ) ∇ c + X α,β ∇ e α · γ αβ ( c, e ) ∇ e β + X α ∇ c · θ α ( c, e ) ∇ e α . (3)Here κ is a symmetric tensor of composition gradient en-ergy coefficients, each γ αβ ( α, β = 1 , . . . , ) is a tensorof strain gradient energy coefficients, and each θ α is atensor of the mixed, composition-strain gradient energycoefficients. Note that, in general, these coefficients willbe functions of the local composition and strain. Thepoint group symmetry of the cubic reference crystal im-poses constraints on the tensor components of κ , γ αβ and θ α as well as on the form of F .While the gradient energies bestow greater accuracyupon the free energy description of solids with non-uniform composition and strain fields, they are essentialat a more fundamental level if the homogeneous free en-ergy density is non-convex. At compositions that render F non-convex, the absence of a gradient energy term willallow spinodal decomposition characterized by composi-tion fluctuations of arbitrary fineness, thus leading tonon-unique microstructures—a fundamentally unphysi-cal result [14]. With κ = , the composition gradient en-ergy penalizes the interfaces wherein composition variesrapidly between high and low limits. This ensures phys-ically realistic results, manifesting in unique microstruc-tures with a mathematically well-posed formulation.An essentially analogous situation exists with respectto the negative curvatures of F in the e - e plane atlow c , which drive the cubic lattice to distort into thetetragonal variants corresponding to the three free en-ergy wells. Consider a solid with a homogeneous freeenergy density as in Figure 3 and a strain state lying be-tween the valleys in the e - e plane at c = 0. Absent thestrain gradient energy, mechano-chemical spinodal de-composition would allow tetragonal variants of arbitraryfineness—an unphysical result, reflecting further math-ematical ill-posedness. Retention of the strain gradientenergy, ( γ αβ = ) penalizes interfaces of sharply vary-ing strain between tetragonal variants to ensure phys-ically realistic results and unique microstructures froma mathematically well-posed formulation. This is well-understood in the literature that studies the formationof martensitic microstructures from non-convex free en-ergy density functions[16–18].
4. Governing equations of non-equilibrium chemistry
The free energy for non-homogeneous composition andstrain fields, Eq. [2], must be a minimum at equilibrium.The state of a solid out of equilibrium will evolve to re-duce the free energy Π[ c, u ]. In formulating a kineticequation for the redistribution of atomic species throughdiffusion, we are guided by variational extremization ofthe free energy to identify the chemical potential, µ . De-tails of this calculation appear as supporting information.The result follows: µ = ∂ F ∂c − ∇ · ( κ ∇ c ) + ∇ c · ∂ κ ∂c ∇ c + X α,β ∇ e α · ∂ γ αβ ∂c ∇ e β + X α (cid:18) ∇ c · ∂ θ α ∂c ∇ e α − ∇ · ( θ α ∇ e α ) (cid:19) . (4)For solids where c tracks the composition of an intersti-tial element within a chemically inert host, such as Liin Li c Mn O , µ in Eq. [4] corresponds to the chemicalpotential of the interstitial element. If c tracks the com-position of a substitutional species, such as in alloys oron sublattices of complex compounds (e.g. the cationsublattice of yttria stabilized zirconia), µ is equal to thechemical potential difference between the substitutionalspecies.The common phenomenological relation for the fluxis J = − L ( c, e ) ∇ µ , where L is the Onsager transporttensor (see de Groot & Mazur [15]). For an interstitialspecies, 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 equationyields the strong form of the governing partial differen-tial equation (PDE) for time-dependent mass transport.It is of fourth order in space due to the composition gra-dient dependence of µ in Equation [4]. See SupportingInformation for strong and weak forms of this PDE.
5. Governing equations of mechanical equilibrium: Straingradient elasticity
Mechanical equilibrium is assumed since elastic wavepropagation typically is a much faster process than dif-fusional relaxation in crystalline solids. Equilibrium isimposed by extremizing the free energy functional withrespect to the displacement field. Standard variationaltechniques lead to the weak and strong forms of straingradient elasticity. The treatment is technical, for whichreason we restrict ourselves to the constitutive relationsthat are counterparts to the chemical potential equation[4] for chemistry. Coordinate notation is used for trans-parency of the tensor algebra, and summation is impliedover repeated spatial index, I . Details appear in Sup-porting Information. The final form of the equations iscomplementary to Toupin’s [11], since our derivation isrelative to the reference crystal, Ω.With the deformation gradient F being related to theGreen-Lagrange strain as E KL = ( F iK F iL − δ KL ), thefirst Piola-Kirchhoff stress tensor and the higher-orderstress tensor respectively, are given by P iJ = X α ∂ ( F + G ) ∂e α ∂e α ∂F iJ + X α ∂ G ∂e α,I ∂e α,I ∂F iJ (5) B iJK = X α ∂ G ∂e α,I ∂e α,I ∂F iJ,K (6)The higher-order stress B , which is absent in classical,non-gradient elasticity [21] (and in earlier treatments ofmechano-chemistry [4, 22]) makes the strong form of gra-dient elasticity a fourth order, nonlinear PDE in space(Supporting Information). The first three-dimensionalsolutions to general boundary value problems of Toupin’sstrain gradient elasticity theory at finite strains were re-cently presented by the authors [23]. II. RESULTSA. Two dimensional examples
We first consider a two-dimensional solid to bet-ter visualize the microstructures that can emerge frommechano-chemical spinodal decomposition. Plane strainelasticity is assumed, for which E , E , E = 0, giv-ing e , e = 0, and e = e / √
2, reducing Equa-tions [1–6] to two dimensions. The discussion on two-dimensional mechano-chemical spinodal decomposition(Figure 2) holds: as a particular parameterization of F ( c, e , e , e ) we consider a regular solution model asa function of c at zero strain. At c = 0, F ( c, e , e , e )is double-welled in e corresponding to the two rectan-gular variants. The gradient energy is G ( ∇ c, ∇ e ) witha constant, isotropic κ = , and constant γ = 0, allother gradient coefficients being zero. While the fullestpossible complexity of the coupling is not revealed bythese simplifications, the aim here is to present the essen-tial physics that is universal to mechano-chemical spin-odal decomposition, postponing details of the more com-plex couplings and tensorial forms to future communi-cations, where they will be derived for specific materialssystems. See Supporting Information for specific formsof F ( c, e , e , e ) and G ( ∇ c, ∇ e ).Figure 4 shows the evolution of microstructure overa 0 . × .
01 domain whose reference (initial) state hasthe square crystal structure and c = 1. The displace-ment component u = 10 − on the right boundary ( X =0 . u = ). Anoutward flux is imposed on the top and bottom bound-aries causing a decrease in composition, c , starting fromthe boundaries. As the composition falls, the homo-geneous free energy density F loses convexity and thestate of the material ( c, e ) enters the mechano-chemicalspinodal along e = 0 (Figure 2b). The continuing out-ward flux first drives material near the top and bottomboundaries fully into the regime where the rectangularcrystal structure is stable. As explained for Figure 2b, FIG. 4: Evolution of 2D microstructure during outflux fromthe top and bottom surfaces of a solid under plane strain.Contours show strain e . Note the legend and correspond-ing square/rectangular variant crystal structures. The defor-mation and accompanying twinned microstructure are seenclearly in the distorted mesh. the negative curvature ∂ F /∂e < e values). The coexistence of theparent, square structure with the daughter, rectangularvariants also can give rise to cross-hatched microstruc-tures that parallel the tweed microstructures describedin the work of Karatha and co-workers [12, 13]. MovieS9 in the supporting information shows the formation ofsuch a microstructure.A laminar, micro-twinned microstructure develops asthe two rectangular variants form, distinguished by thesign of e (see legend). Note that e remains continu-ous because of the penalization of strain gradients ∇ e ,although discontinuities can develop in ∇ e , itself. Be-cause our numerical framework uses basis functions thatare continuous up to their first derivatives (see Methodsand Supporting Information), ∇ e is indeed discontinu-ous in the computations. The lamination accommodatesthe strain difference between the two rectangular variantsto minimize the free energy. If strains were infinitesimal( | e | , | e | , | e | & e would correspond to shear alongdirections that are rotated π/ − . ≤ e ≤ . F .In some cases, the long-range nature of elasticity forceslike rectangular variants to align even when separatedby an untransformed square phase. This is first seenin the finger-like extensions of strain contours from therectangular variants into the square phase in Figure 4,followed by their alignment and eventual incorporationinto laminae of the same variant (top right panel and itsevolution shown as an inset). In this instance, althoughthe micro-twins end at an invariant habit plane that iscommon with the untransformed square phase, the latticeparameters of the rectangular variants differ from thoseof the square phase, inducing elastic strain energy in thelatter. The alignment lowers this strain energy, and alsois seen in Figure 5a. Note, however, that this is not auniversal feature. Movie S2 shows instances where unlikevariants come close to alignment. FIG. 5: Microstructure controlled by elastic gradient lengthscale parameter ( l e ). Shown are the final contours of e forthe simulation of (a) outflux from top and bottom surfaces,and (b) quenching. The fineness of laminae depends on thestrain gradient elasticity length scale l e ∼ r γ / qP α,β =1 , , ( ∂ F /∂e α ∂e β ) , as explored inFigure 5. For Figure 5a, the initial and boundaryconditions are the same as in the example of Figure4. Decreasing l e weakens the penalty on the straingradient ∇ e across neighboring, unlike rectangularvariants and allows more twin boundaries. Notably,self-similarity is not maintained between microstructuresfor different l e , even for the same initial and boundaryvalue problem. We understand this to be the influenceof elastic strain accommodation: To minimize the totalfree energy when the strain gradient penalty changes,the physics optimizes twin boundaries via laminations ofdifferent sizes as well as different patterns. Importantlyhowever, crystal symmetry admits non-vanishing strain gradient energy coefficients beyond γ = 0 used in thesesimulations (Supporting Information). Furthermore, thecomposition dependence of F could be more complexthan the simple regular solution model used here. Giventhe already strong effect of l e alone, we conjecture thatvarying these forms will have a significant influence onthe resulting coherent twin microstructure. The properform, while guided by crystal symmetry argumentsmust ultimately be determined experimentally or byfirst-principles statistical mechanical methods.The example in Figure 5b further explores the influ-ence of l e . In this suite of computations, the unstrainedmaterial with an initially square microstructure, convexfree energy density and composition having random fluc-tuations about c = 0 .
45 was quenched to a low tempera-ture and into the mechano-chemical spinodal. The localcomposition and strain evolve under the thermodynamicdriving forces detailed in the context of Figure 2. Wedraw attention, once again, to the changing identity ofrectangular variants, their shapes and sizes depending on l e . Also note the progressively finer lamination of rectan-gular domains with decreasing l e . Such studies suggesthow dynamic mechano-chemical spinodal decompositioncan lead to an atlas of microstructures, which in turn willdetermine material properties. Movies S1-S4 in the sup-porting information show the evolution of some of thesemicrostructures. B. A three dimensional study
This final example displays the full three-dimensionalcomplexity of microstructures resulting from themechano-chemical spinodal. We persist with the abovesimplifications aimed at presenting the essential physicsof mechano-chemical spinodal decomposition, and post-pone a full exploration of the coupling and anisotropiesin coefficients to future work on specific materials sys-tems. The free energy density function used for the three-dimensional study appears as supporting information.Figure 6 shows the equilibrium microstructure that re-sults in a solid that is initially in the cubic phase, c = 1,subject to an outflux on all surfaces. The domain is a unitcube with displacement components u , u = 0 .
01 on theboundary X = 1, with zero displacement, u = , onthe boundary X = 0. The cubic to tetragonal trans-formation takes place as c →
0. Under these condi-tions all three tetragonal variants form, with an intri-cately interleaved microstructure for strain accommoda-tion. Note the finer microstructures and changing pat-tern for smaller l e . The inset shows the surface strain-ing around a corner, delineated by the distorted meshlines. Observe the three, oriented, tetragonal variantsformed by twinning deformation from the initial cubicstructure. See Movie S7 in the supporting informationfor a detailed view of the three-dimensional structure ofthese individual tetragonal variants. Movie S10 in sup-porting information shows other three-dimensional mi-crostructures whose cross-sectional planes bear closer re-semblance to the plate-like structures predicted by two-dimensional computations. To the best of our knowledge,such computations of a cubic to tetragonal transforma-tion with twinned variants whose microstructure is con-trolled by nonlinear, strain gradient interface energies,have not been previously presented. FIG. 6: Microstructure observed in 3D for different values ofthe elastic gradient length scale parameter ( l e ). The threetetragonal variants appear in blue (variant 1), yellow (variant2) and red (variant 3) for c = 0. The transformation strainsare easily discerned in the distorted mesh. III. DISCUSSION
Several kinetic mechanisms have been proposed to de-scribe the decomposition of a solid upon entering a two-phase region [42, 43]. A commonly observed mechanismoccurs through localized nucleation followed by diffu-sional growth that is mediated by the migration of in-terfaces separating the daughter phases from the parentphase. Nucleation and growth mechanisms are commonwhen the phases participating in the transformation havelittle or no crystallographic relation to each other. Thetransformation then proceeds reconstructively, where theinterphase interfaces are disordered or at best only semi-coherent. Numerical treatments of this mechanism relyon sharp-interface methods such as a level set approach[44].Other kinetic mechanisms of decomposition involvingsome degree of coherency are also possible when the par-ticipating phases share sufficient crystallographic com-monality that order parameters can be defined describingcontinuous pathways connecting the parent and daugh-ter phases. The structural changes of the crystal thataccompany these transformations can fall in one of twocategories. One subclass of structural transformations isdriven by an internal shuffle, where the atomic arrange-ment within the unit cell of the crystal undergoes a sym-metry breaking change. The unit cell vectors of the crytalmay also undergo a change and lower the symmetry of the lattice, but this effect is secondary and in response to theatomic shuffle within the unit cell. The order parametersfor the structural change therefore describe the atomicshuffles within the unit cell and are non conserved. Thesecond subclass of structural transformations are drivenby a symmetry reducing change of the lattice vectors ofthe crystal, with any atomic rearrangement of the ba-sis of the crystallographic unit cell occurring in responseto the symmetry breaking changes of the unit cell vec-tors. The natural order parameters describing these tran-sitions are strains. The free energy in the first subclasswill exhibit negative curvatures with respect to the non-conserved shuffle order parameters, while the free energyin the second subclass will have negative curvatures withrespect to the relevant strain order parameters, as wasrecently predicted for the cubic to tetragonal transitionof ZrH [33].Decomposition reactions that require both diffusionalredistribution and a structural change due to an internalshuffle have been treated successfully and rigorously withphase field approaches that combine Cahn-Hilliard withAllen-Cahn. The Cahn-Hilliard description accounts forthe composition degrees of freedom while the Allen-Cahncomponent describes the evolution of the non-conservedshuffle order parameters required to distinguish the vari-ous phases of the transformation [5–8]. These treatmentshave included strain energy as a secondary effect servingonly as a positive contribution to the overall free energydue to coherency strains. The approach is rigorous if thefree energy density remains convex with respect to strain,with instabilities in the free energy appearing as a func-tion of concentration and the non-conserved shuffle orderparameters.Here we have presented a mathematical formulationand an accompanying computational framework for de-composition reactions that combine diffusional redistri-bution with a structural change driven by a symmetrybreaking strain of the crystallographic unit cell as op-posed to an internal shuffle within the unit cell. Strainsplay a primary role, explicitly serving as order parame-ters to distinguish variants of a daughter phase that hasa symmetry subgroup-group relationship with its parentphase due to a structural change of the crystallographicunit cell. Crucial to the description is that the drivingforce for the formation of the daughter from a super-saturated parent phase emerges from an instability withrespect to composition. The treatment therefore com-bines Cahn-Hilliard for composition with a descriptionof martensitic transformations at fixed concentration in-troduced by Barsch and Krumhansl [9]. The existenceof simultaneous instabilities with respect to strain andcomposition has not been considered as a mechanism todescribe decomposition reactions upon quenching into atwo-phase region.The possibility of such a mechano-chemical spinodalmechanism is motivated by recent first-principles studiesof the cubic to tetragonal phase transformations of ZrH ,where the free energy of the high temperature cubic formwas predicted to become unstable with respect to strainupon cooling below the cubic to tetragonal second ordertransition temperature [33]. The ZrH − c hydride can ac-commodate large concentrations of hydrogen vacancies, c , and has a phase diagram that is topologically identicalto that depicted in Figure 2c, with a two-phase regionseparating a hydrogen-rich tetragonal form of ZrH − c α from a cubic form of ZrH − c β (with c α < c β ). See Zuzeket al. [32].[45] To be consistent with the predicted freeenergies for stoichiometric ZrH (i.e. c =0) and the ex-perimental T versus c phase diagram with the form ofFigure 2c, the free energy of this hydride as a functionof composition and strain (i.e. e and e ) should be sim-ilar to those depicted in Figures 2a and b, and Figure3. Decomposition upon quenching cubic ZrH − c intothe two-phase region can therefore proceed through amechano-chemical spinodal decomposition mechanism.Not only does the phenomenological description in-troduced in this work describe a new mechanism ofdecomposition that is qualitatively distinct from previ-ous combined Cahn-Hilliard and Allen-Cahn approaches,its numerical solution also proves to be substantiallymore complex due to contributions from strain gradi-ent terms. Indeed, even numerical solutions to general,three-dimensional boundary value problems of gradientelasticity at finite strain were not available until pre-sented recently by the authors [23]. Furthermore, asother authors have pointed out before, the use of strainmetrics as order parameters makes a reliance on infinites-imal strains untenable due to the rigid rotations that ac-company the finite strains characterizing most structuraltransformations [8, 41]. The use of finite strain met-rics introduces geometric non-linearity into the problem,which while having been treated in past Cahn-Hilliardand Allen-Cahn approaches [8], adds additional numeri-cal challenges when also considering strain gradient con-tributions. These challenges have been overcome in thiswork as demonstrated by our three-dimensional numeri-cal examples.Mechano-chemical spinodal decomposition as de-scribed here can be expected in solids forming hightemperature phases that exhibit dynamical instabili-ties at low temperature. An accumulating body offirst-principles calculations of Born-Oppenheimer sur-faces have shown that many high temperature phasesare dynamically unstable at low temperature [26, 27], be-coming stable at high temperature through large anhar-monic vibrational excitations [28–31, 33], usually througha second order phase transition. While such instabilitiesare frequently dominated by phonon modes describingan internal shuffle, a subset of chemistries become dy-namically unstable at low temperature with respect tophonon modes that break the symmetry of the crystalunit cell [30, 33], an instability that can be described phe-nomenologically with strain order parameters. If, uponalloying such compounds, the high symmetry phase be-comes stable by passing through a two-phase region, itsfree energy surface will resemble that of Figures 2b and 3, and thereby make the solid susceptible to the mechano-chemical spinodal decomposition mechanism describedhere.While first-principles and experimental evidence sug-gests mechano-chemical spinodal decomposition shouldoccur in ZrH − c , we expect it to occur in a wide rangeof other chemistries as well. One possible example, asdescribed in the introduction, is the decomposition ofcubic yttria stabilized zirconia [34, 35] upon quenchingfrom the high temperature cubic phase into a two-phaseregion separating a low-Y tetragonal phase from a Y-richcubic phase. Past treatments of this transformation [5–7]relied on non-conserved order parameters to distinguishthe different tetragonal variants from each other and fromthe cubic parent phase. The elastic part of the free en-ergy density function was parameterized using linearizedelasticity and infinitesimal strains, as other authors havealso done [39, 40]. The chemical part of the free energydensity was assumed to have a negative curvature as afunction of the non-conserved order parameter at Zr-richcompositions, but made up of convex (and quadratic)potentials with respect to strain. We reiterate that sucha treatment rests on the implicit assumption that inter-nal shuffles drive the cubic to tetragonal transformation,while the free energy as a function of strain remains con-vex for all relevant values of the non-conserved order pa-rameters.The possibility also exists that a coarse-grained freeenergy density of cubic ZrO exhibits negative curva-tures with respect to strain below the cubic to tetragonaltransition temperature making the formalism developedhere an accurate description of decomposition reactionsof yttria stabilized zirconia. While the precise mecha-nism and nature of the cubic to tetragonal transition ofpure ZrO remains to be resolved [29], first-principlescalculations predict that the cubic form of ZrO is dy-namically unstable with respect to transformation to thetetragonal variant [36]. A rigorous statistical mechan-ics treatment [28, 33] is required to determine whetherthe cubic to tetragonal transition of pure ZrO is accom-panied by a change in the sign of the curvature of thefree energy density with respect to e and e . If thisproves to be the case, yttria stabilized zirconia shouldalso be susceptible to mechano-chemical spinodal decom-position upon quenching, consistent with the coherentspinodal microstructures between tetragonal and cubicphases observed in single crystal regions of quenched cu-bic Zr − c Y x O − c/ [35].A mechano-chemical spinodal may also play a rolein a variety of important electrode materials for Li-ion batteries, and intercalation compounds consideredfor two-dimensional nano-electronics. These include cu-bic LiMn O transforming to tetragonal Li Mn O [37].While most mechano-chemical spinodal transitions willoccur in three dimensions, many of the qualitative fea-tures of these transitions are more conveniently illus-trated in two dimensions. Our two-dimensional studiesshould also prove relevant to understanding mechano-0chemical phase transformations in two-dimensional lay-ered materials for nano electronics. Materials such asTaS , are susceptible to Peierls instabilities upon varia-tion of the composition of adsorbed or intercalated guestspecies that donate to or extract electrons from the sheet-like host [38]. Our phenomenological treatment by intro-duction of the concept of a mechano-chemical spinodal,coupled with gradient stabilization of the ensuing non-convex free energy in strain-composition space, thus of-fers a framework with potential for extension to a widerange of phase transformation phenomena.With regard to the fundamental thermodynamic un-derpinnings of analogous processes, the phenomenologi-cal model introduced here can also be used to describetemperature driven martensitic transformations. Thecomposition, c , would then be analogous to the internalenergy density, U , while the chemical potential, µ , wouldbe analogous to the temperature, T . Due to the presenceof temperature gradients throughout the solid, though,the starting point must be a formulation of the entropydensity, S , (as opposed to the Helmholtz free energy) asa functional of spatially varying U and displacements, u .In analogy with Equation (2), the total entropy can bewritten as a volume integral over a homogeneous entropydensity that depends on the local internal energy densityand strain ¯ S ( U , e ), as well as a non-uniform entropydensity contribution expressed in terms of gradients ininternal energy density and strain ( ∇ U , ∇ e ). Varia-tional maximization of the entropy will yield mechani-cal equilibrium equations as well as an expression for thetemperature T (strictly speaking, for its reciprocal, 1 /T ),similar to Equation (4) for the chemical potential. Dueto the presence of gradient terms, ∇ U and ∇ e , the tem-perature will not only be a function of the local internalenergy density and strain, but also gradients in these fieldvariables. As with the diffusion problem treated here,heat flow can be described with an Onsager expressionrelating the heat flux to a gradient in temperature. Thepossibility exists that the homogeneous entropy ¯ S ex-hibits instabilities with respect to internal energy density U , allowing for spinodal decomposition with respect tothe redistribution of U in a manner that is similar to thewell understood problem of spinodal decomposition withrespect to composition. The treatment introduced herecan therefore describe (upon replacing c with U and µ with T ) temperature driven martensitic transformationsfor solids that exhibit instabilities with respect to bothinternal energy density and strain. Past treatments ofthis phenomenon also relied on a Barsch and Krumhanslapproach, solving the mechanical problem together withFourier’s law of heat conduction [39, 40]. However, thesestudies were restricted to two dimensions [39–41] withsome neglecting geometric non-linearity by relying on in-finitesimal strains [39, 40]. At a more fundamental level,they did not consider the possibility of spinodal instabil-ities with respect to the redistribution of internal energydensity. IV. METHODS
The governing fourth-order PDEs of mass transportand gradient elasticity are solved numerically in weakform, wherein second-order spatial gradients appear onthe trial solutions and their variations. Numerical so-lutions require basis functions that are continuous upto their first spatial gradients, at least. Our numeri-cal framework (see Rudraraju et al. [23]) uses isogeo-metric analysis [24] and the PetIGA code framework [25]with spline basis functions, which can be constructed forarbitrary degree of continuity (Supporting Information).This framework has proven pivotal to the current work.The code used for the numerical examples in this papercan be downloaded at the University of Michigan Com-putational Physics Group open source codes webpage:
V. FUNDING
The mathematical formulation for this work was car-ried out under an NSF CDI Type I grant: CHE1027729“Meta-Codes for Computational Kinetics” and NSFDMR 1105672. The numerical formulation and computa-tions have been carried out as part of research supportedby the U.S. Department of Energy, Office of Basic EnergySciences, Division of Materials Sciences and Engineeringunder Award [1] Cahn, J. W., and Hilliard, J. E. Free energy of a nonuni-form system. i interfacial energy.
The Journal of Chem-ical Physics , 28:258–267 (1958).[2] Allen, S. M., and Cahn, J. W. A microscopic the-ory for antiphase boundary motion and its applicationto antiphase boundary coarsening.
Acta Metallurgica ,27:1085–1091 (1979).[3] Bhattacharya, K., Conti, S., Zanzotto, G., and Zimmer,J. Crystal symmetry and the reversibility of martensitictransformations.
Nature , 428:55–59 (2004). [4] Voorhees, P. W., and Johnson, W. C. The thermodynam-ics of elastically stressed crystals.
Solid State Physics:Advances in Research and Applications , 59:1–201 (2004).[5] Wang, Y., Wang, H., Chen, L. Q., and Khatchaturyan,A. G. Shape evolution of a coherent tetragonal precip-itate in partially stabilized cubic zirconia: A computersimulation.
Journal of the American Ceramic Society ,76: 3029-3033 (1993).[6] Fan, D., and Chen, L. Q. Computer simulation of twinformation during the displacive c → t phase transforma- tion in the Zirconia-Yttria system. Journal of the Amer-ican Ceramic Society , 78:769–773 (1995).[7] Fan, D., and Chen, L. Q. Possibility of spinodal decompo-sition in
ZrO − Y O alloys: A theoretical investigation. Journal of the American Ceramic Society , 78:1680–1686(1995).[8] Hildebrand, F. E. and Miehe, C. A phase field modelfor the formation and evolution of martensitic laminatemicrostructure at finite strains.
Philosophical Magazine ,92:4250-4290 (2012).[9] Barsch, G. R., and Krumhansl, J. A. Twin boundaries inferroelastic media without interface dislocations.
PhysicalReview Letters , 53:1069–1072 (1984).[10] van der Waals, J. D. Thermodynamische theorie der cap-illariteit in onderstelling van continue dichtheidsveran-dering.
Verhanderlingen der Koninklijke NederlandseAkademie , 1:3–56 (1893).[11] Toupin, R. A. Elastic materials with couple-stresses.
Archive for Rational Mechanics and Analysis , 11:385–414(1962).[12] Kartha, S., Castan, T., Krumhansl, J. A., and Sethna, J.P. Spin-glass nature of tweed precursors in martensitictransformations.
Physical Review Letters , 67:3630–3633(1991).[13] Kartha, S., Krumhansl, J. A., Sethna, J. P., and Wick-ham, L. K. Disorder-driven pretransitional tweed pat-tern in martensitic transformations.
Physical Review B ,52:803–827 (1995).[14] Hilliard, J. E.
Phase transformations . H. I. Aronson, Ed.ASM, 470–560 (1970).[15] de Groot, S. R., and Mazur, P.
Non-equilibrium Ther-modynamics . Dover (1984).[16] Ball, J. M., and James, R. D. Fine phase mixtures asminimizers of energy.
Archive for Rational Mechanicsand Analysis , 100(1):13–52 (1987).[17] M¨uller, S.
Variational models for microstructure andphase transitions , volume 1713 of
Lecture Notes in Math-ematics . Springer Berlin Heidelberg (1999).[18] Ball, J. M., and Crooks, E. C. M. Local minimizers andplanar interfaces in a phase-transition model with inter-facial energy.
Calculus of Variations and Partial Differ-ential Equations , 40(3-4):501–538 (2011).[19] Van der Ven, A., Bhattacharya. J, and Belak, A. A.Understanding li diffusion in li-intercalation compounds.
Accounts of Chemical Research , 46:1216–1225 (2013).[20] Van der Ven, A., Yu, H. C., Ceder, G., and Thornton,K. Vacancy mediated substitutional diffusion in binarycrystalline solids.
Progress in Materials Science , 55:61–105 (2010).[21] Truesdell, C., and Noll, W.
The nonlinear field theoriesof mechanics . Springer Verlag, Berlin, Heidelberg (1965).[22] Larche, F., and Cahn, J. W. Non-linear theory of ther-mochemical equilibrium of solids under stress.
Acta Met-allurgica , 26:53–60 (1978).[23] Rudraraju, S., Van der Ven, A., and Garikipati,K. Three-dimensional isogeometric solutions to generalboundary value problems of Toupin’s gradient elasticitytheory at finite strains.
Computer Methods in AppliedMechanics and Engineering , 278:705–728 (2014).[24] Cottrell, J., Hughes, T. J .R., and Bazilevs, Y. Isogeo-metric Analysis: Toward Integration of CAD and FEA.
Wiley (2009).[25] Collier, N., Dalc´ın, L., and Calo, L. M. Petiga:High-performance isogeometric analysis arXiv:1305.4452 (2013).[26] Craievich, P. J., Weinert, M., Sanchez, J. M., and Wat-son, R. E. Local stability of nonequilibrium phases.
Phys-ical Review Letters , 72:3076–3079 (1994).[27] Grimvall, G., Magyair-Koepe, B., Ozolins, V., and Pers-son, K. A. Lattice instabilities in metallic elements.
Re-views of Modern Physics , 84:945–986 (2012).[28] Zhong, W., Vanderbilt, D., and Rabe, K. A. Phase tran-sitions in
BaT iO from first-principles. Physical ReviewLetters , 73:1861–1864 (1994).[29] Fabris, S., Paxton, A. T., and Finnis, M. W. Free en-ergy and molecular dynamics calculations for the cubic-tetragonal phase transition in zirconia.
Physical ReviewB , 63:094101 (2001).[30] Bhattacharya, J. and Van der Ven, A. Mechanical in-stabilities and structural phase transitions: The cubic totetragonal transformation.
Acta Materiala , 56:4226–4232(2008).[31] Souvatzis, P., Eriksson, O., Katsnelson, M. I., and Rudin,S. P. Entropy driven stabilization of energetically un-stable crystal structures explained from first principlestheory.
Physical Review Letters , 100:095901 (2008).[32] Zuzek, E., Abriata, J. P. and San-Martin, A. and Manch-ester, F. D. The H-Zr (hydrogen-zirconium) system
TheBulletin of Alloy Phase Diagrams , 11:385 (1990).[33] Thomas, J. C., and Van der Ven, A. Finite-temperatureproperties of strongly anharmonic and mechanically un-stable crystal phases from first principles.
Physical Re-view B , 88:214111 (2013).[34] Chien, F. R., Ubic, F. J., Prakash, V., and Heuer, A. H.Stress-induced martensitic transformation and ferroelas-tic deformation adjacent microhardness indents in tetrag-onal zirconia single crystals.
Acta Materialia , 46:2151–2171 (1998).[35] Krogstad, J. A., et al. Phase stability of t’-Zirconia-basedthermal barrier coatings: Mechanistic insights.
Journalof the American Ceramic Society , 94:S168–S177 (2011).[36] G. Jomard., et al. First-principles calculations to describezirconia pseudopolymorphs.
Physical Review B , 59:4044–4052 (1999).[37] Thackeray, M. M. Manganese oxides for lithium batter-ies.
Progress in Solid State Chemistry , 25:1–71 (1997).[38] Rossnagel, K. Suppression and emergence of charge-density waves at the surfaces of layered 1 T − T iSe and1 T − T aS by in situ Rb deposition. New Journal ofPhysics , 12:125018 (2010).[39] Bouville, M., and Ahluwalia, R. Effect of lattice-mismatch-induced strains on coupled diffusive and dis-placive phase transformations.
Physical Review B ,75:054110 (2007).[40] Maraldi, M., Wells, G., and Molari, L. Phase field modelfor coupled displacive and diffusive microstructural pro-cesses under thermal loading.
Journal of the Mechanicsand Physics of Solids , 59:1596–1612 (2011).[41] Alphonse, F., Le Bouar. Y., Gaubert, A., and Salman.U. Phase field methods: Microstructures, mechanicalproperties and complexity.
Comptes Rendus Physique ,11: 245–256 (2010).[42] Balluffi, R. W., Allen, A. M. and Carter, W. C. Kineticsof materials John Wiley & Sons, Hoboken, New Jersey(2005).[43] Porter, D. A. and Easterling, K. E. Phase transforma-tions in metals and alloys. Nelson Thorne, Ltd., UK(2001). [44] Sethian, J. A. Level set methods and fast marching meth-ods: Evolving interfaces in computational Geometry,fluid mechanics, computer vision and materials science.Cambridge University Press, Cambridge, UK (1999). [45] In their phase diagrams, Zuzek et al have an invertedcomposition axis relative to our notation. In their work,the tetragonal phase is labelled ε and the cubic phase is δ . upplementary Material for Mechano-chemical spinodal decomposition: Aphenomenological theory of phase transformation in multi-component, crystallinesolids Shiva Rudraraju , Anton Van der Ven, Krishna Garikipati
PACS numbers:Keywords:
I. THE EXISTENCE OF A FREE ENERGYDENSITY FUNCTION THAT IS NON-CONVEXWITH RESPECT TO COMPOSITION
Consider the experimentally derived equilibrium phasediagram for transition metal hydrides [1] shown schemat-ically in Figure 2c of the main text. The single cubicphase, stable for all composition and strains at high tem-perature, only remains stable at high composition whenquenched to low temperature. At low composition, thetetragonal phase is stable at low temperature. Also re-call that first principle calculations have shown that thetetragonal phase is unstable with respect to strain anddecomposes into three energetically equivalent variants[2]. The homogeneous free energy density therefore issymmetric with respect to zero strain. The two-phasecoexistence region implies that the free energy densitysurface, F in strain-composition, ( e , c ) space admits acommon tangent hyperplane. For clarity, let us considerthe zero strain hyperplane in this space. Its projectionappears in Figure S1. For a point in the solid with com-position c ∗ , the common tangent construction leads tocompositions c α and c β , respectively of the tetragonaland cubic phases, with ∂ F ∂c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) c α = ∂ F ∂c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) c β . (1)Furthermore, the stability of the cubic and tetragonalphases with respect to composition imply that ∂ F ∂c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) c α > , ∂ F ∂c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) c β > . (2)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 ∂ F ∂c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) c + α > ∂ F ∂c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) c α , and ∂ F ∂c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) c − β < ∂ F ∂c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) c β . (3)But then, from Equations (1) and (3), it follows that, fora smooth F that is at least twice differentiable in theinterval ( c α , c β ), there must exist at least one value c γ with c + α ≤ c γ ≤ c − β satisfying ∂ F ∂c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) c γ < . (4) Therefore there is at least one intermediate value of com-position in the two-phase coexistence region at which thefree energy density is non-convex with respect to compo-sition. This result is simply a consequence of the meanvalue theorem of calculus.FIG. S1: Figure S1 : Projection of the free energy onto the zero strain hyperplane.
II. FINITE STRAIN KINEMATICS
In the interest of clarity and consistency with the maintext, the Euclidean basis of our Cartesian coordinate sys-tem is denoted by X i , i = 1 , . . . X i · X j = δ ij . Thisconstitutes a mild abuse of notation, because referencepositions have also been denoted by X , X , X . Thehigh symmetry, cubic reference configuration of the crys-tal is denoted by Ω, on which material points are labelledby their position vectors X = X X + X X + X X .The deformation map from Ω to the current configurationof the crystal is ϕ ( X , t ) = X + u , where the displace-ment field u was introduced in the main text. The de-formation gradient tensor is F = ∂ ϕ /∂ X = + ∂ u /∂ X .Using lower case indices to denote the components ofvectors and tensors on the current configuration, and up-per case indices for those on the reference configuration,we have the coordinate expression F iJ = ∂ϕ i /∂X J = δ iJ + ∂u i /∂X J .The Green-Lagrange strain tensor is E = ( F T F − ),where is the second-order isotropic tensor. We re-peat this relation in coordinate notation from the maintext: E IJ = ( F iI F iJ − δ IJ ). The polar decomposi-tion theorem states that the deformation gradient canbe decomposed multiplicatively as F = RU , where R is an (orthogonal) rotation tensor belonging to the space a r X i v : . [ phy s i c s . c h e m - ph ] S e p SO(3). It therefore satisfies R T R = . The other ten-sor in this decomposition, U , is symmetric and positive-definite. From the polar decomposition it follows that E as defined above is invariant to rotation of the currentconfiguration. The requirement of material frame invari-ance dictates that the dependence of the free energy den-sity on elasticity must be as a function of E . Such a freeenergy density function is invariant to rotations on thecurrent configuration. III. THE CHEMICAL POTENTIAL:VARIATIONAL DERIVATION
Since the crystal is, in general, far from chemicalequilibrium we follow the standard treatment of non-equilibrium thermodynamics to model the chemistry. Inorder to obtain a relation for the chemical potential weconsider variations on the composition field of the form c ε := c + ε ξ , where ξ is a function such that, on theDirichlet boundary, Γ c ⊂ Γ, we have c = ¯ c and ξ = 0. δδc Π[ c, u ] = d Π[ c + εξ, u ] dε (cid:12)(cid:12)(cid:12)(cid:12) ε =0 (5)Carrying out the differentiation with respect to ε andintegration by parts yields δ Π δc = Z Ω ξ ∂ F ∂c − ∇ · ( κ s ∇ c ) + ∇ c · ∂ κ ∂c ∇ c + X α,β ∇ e α · ∂ γ αβ ∂c ∇ e β + X α (cid:18) ∇ c · ∂ θ α ∂c ∇ e α − ∇ · ( θ α ∇ e α ) (cid:19) ! d V + Z Γ ξ κ s ∇ c + X α θ α ∇ e α ! · n d S, (6)where n is the outward normal vector to Γ. At chemicalequilibrium we have δ Π /δc = 0, yielding two conditions: Z Ω ξ " ∂ F ∂c − ∇ · ( κ s ∇ c ) + ∇ c · ∂ κ ∂c ∇ c + X α,β ∇ e α · ∂ γ αβ ∂c ∇ e β + X α (cid:18) ∇ c · (cid:18) ∂ θ α ∂c ∇ e α (cid:19) − ∇ · ( θ α ∇ e α ) (cid:19) d V = 0(7) Z Γ ξ " κ s ∇ c + X α θ α ∇ e α · n = 0 (8) Standard variational arguments lead to the conclusionthat, for chemical equilibrium, the quantities within thebrackets must vanish on the corresponding domains Ωand Γ. The expression contained within brackets in theintegrand of Equation [7] is the chemical potential µ . Inthe non-equilibrium setting that is of interest here, wecontinue to be guided by the above treatment and usethe same definition for the chemical potential [3]: µ = ∂ F ∂c − ∇ · ( κ s ∇ c ) + ∇ c · ∂ κ ∂c ∇ c + X α,β ∇ e α · ∂ γ αβ ∂c ∇ e β + X α (cid:18) ∇ c · ∂ θ α ∂c ∇ e α − ∇ · ( θ α ∇ e α ) (cid:19) . (9)Variants of Equation [8] appear in a proper variationalderivation of chemical equilibrium—a detail that is some-times overlooked in phase field treatments. We will im-pose chemical equilibrium of the boundary Γ, by requir-ing the term in brackets to vanish.With the standard, phenomenological relation J = − L ( c, e ) ∇ µ for the flux, the strong form of the governingPDE of mass balance follows: ∂c∂t + ∇ · J = 0 in Ω − J · n = J on Γ( κ s ∇ c + X α θ α ∇ e α ) · n = 0 on Γ c ( X ) = c ( X ) at t = 0 (10)Note that we have assumed that there is no Dirichletboundary; i.e. Γ c = ∅ , in stating the final strong formof this PDE. The conjugate boundary condition on themass flux holds on the entire boundary Γ, in Equation[10 ]. Additionally, we have the boundary condition aris-ing from requiring chemical equilibrium on Γ, in Equation[10 ]. Using standard variational arguments, the weakform of mass balance can be obtained by starting from[10]. Finally, we have the initial condition on composi-tion, Equation [10 ]. The substitution of Equation [9] in J = − L ∇ µ , followed by this phenomenological flux re-lation’s substitution in (10 ) shows that this is a fourthorder PDE in space. A. On mass transport posed in the referenceconfiguration
We note that the mathematical formulation of diffusionis commonly stated in the current (deformed) configura-tion, which is the true state of the solid. However, forcrystalline solids (with negligible extended defects) it ismore convenient to define Onsager transport coefficientsand evaluate composition gradients on the reference (un-deformed) configuration, here denoted by Ω. Diffusion inthe crystalline solids results from atomic hops betweenwell-defined crystal sites, which can always be mappedonto a reference crystal structure. Increasingly, mobilitycoefficients in multi-component solids are calculated fromfirst principles through the evaluation of Kubo-Green ex-pressions in kinetic Monte Carlo simulations on the refer-ence crystal in configuration Ω. Since the mobility tensoris reported on Ω, it is most convenient to adopt the La-grangian description for diffusion over this configurationin single crystalline solids.Our numerical implementation is based on the weakform, whose derivation we begin by first reintroducing ξ , now in the guise of a weighting function lying in thefunction space Q = { ξ ∈ H (Ω) } . We seek functions c ∈ P = { c ∈ H (Ω) } such that Z Ω ξ ∂c∂t d V + Z Ω ∇ ξ · L ∇ (cid:18) ∂ F ∂c + ∂ G ∂c (cid:19) d V + Z Ω ∇ · ( L ∇ ξ ) ∇ · ( κ ∇ c )d V − Z Γ ξJ d S − Z Γ ( L ∇ ξ ) · n ∇ · ( κ ∇ c )d S = 0(11)Equation [11] is obtained by multiplying ξ into [10 ], in-tegrating by parts twice, and using the Neumann bound-ary condition [10 ]. Note that the higher-order Dirich-let boundary condition (10 ) has not been built into thefunction spaces as is commonly done in weak forms. In-stead, it is imposed weakly using the method of Nitsche[4–6] and extending the weak form to Z Ω ξ ∂c∂t d V + Z Ω ∇ ξ · L ∇ (cid:18) ∂ F ∂c + ∂ G ∂c (cid:19) d V + Z Ω ∇ · ( L ∇ ξ ) ∇ · ( κ ∇ c )d V − Z Γ ξJ d S − Z Γ ( L ∇ ξ ) · n ∇ · ( κ ∇ c )d S − Z Γ ∇ · ( κ ∇ ξ ) ( L ∇ c ) · n d S + Z Ω τ ∇ ξ · n ∇ c · n d S = 0 (12)Equation [12] is based on the coupling tensor θ α = for all α , and for κ , γ αβ independent of c , correspondingto the numerical examples presented in the main text.This weak form maintains consistency. Note that thesecond-to-last term that has been added over the formin Equation [11] is symmetric with the term precedingit. This formulation shows better numerical performancethan the also-consistent formulation in which this addi-tional term is preceded by a negative sign. Finally, τ is a stabilization parameter that is inversely proportional tothe element size. IV. MECHANICAL EQUILIBRIUM
As explained in the main text, since elastic wave prop-agation is orders of magnitude faster than the kineticsof mass transport, it can be assumed that mechanicalequilibrium is attained instantaneously for the prevail-ing composition field. To obtain the relevant equilibriumequations, we introduce variations on the displacementfield: u ε = u + ε w , such that on the Dirichlet boundaryΓ u i , the fields satisfy u i = ¯ u i and w i = 0. We constructthe first variation of Π[ c, u ] with respect to u . δδ u Π[ c, u ] = ddε Z Ω ( F ( c, e ε ) + G ( c, e ε , ∇ c, ∇ e ε )) d V − X i =1 Z Γ Ti u i ε T i d S !(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ε =0 (13)where e ε are obtained by first replacing u by u ε in thekinematic relations for E discussed above, which arethen substituted in the relations for the reparameterizedstrains e . We recall the constitutive relations for the firstPiola-Kirchhoff stress and the higher-order stress, respec-tively: P iJ = X α ∂ ( F + G ) ∂e α ∂e α ∂F iJ + X α ∂ G ∂e α,I ∂e α,I ∂F iJ (14) B iJK = X α ∂ G ∂e α,I ∂e α,I ∂F iJ,K (15)Using these constitutive relations, Equation [13] can bemanipulated to yield the Euler-Lagrange equation, whichis also the weak form of mechanical equilibrium for astrain gradient elastic material, on Ω in terms of P and B . We retain coordinate notation in the interest of trans-parency: Z Ω ( P iJ w i,J + B iJK w i,JK ) d V − X i =1 Z Γ Ti w i T i d S = 0(16)We introduce the normal and surface gradient operators, ∇ n and ∇ s , defined by ∇ n ψ = ∇ ψ · n , and ∇ s ψ = ∇ ψ − ( ∇ n ψ ) n . (17)Starting from Equation [16] and applying integration byparts several times we have the strong form of mechani-cal equilibrium for a strain gradient elastic material. SeeRudraraju et al. [7] for detailed steps. The results of thisvariational treatment are analogous to those of Toupin[8], except that it is carried out over the reference con-figuration of the crystal, rather than its current configu-ration. P iJ,J − B iJK,JK = 0 in Ω u i = ¯ u i on Γ u i P iJ n J − ∇ n B iJK n K n J − ∇ sJ ( B iJK ) n K − B iJK ∇ sJ ∇ sK + ( b LL n J n K − b JK ) B iJK = T i on Γ T i B iJK n J n K = 0 on Γ J n Γ J n K B iJK K = 0 on Υ L i (18)Here, Γ = Γ u i ∪ Γ T i for i = 1 , ,
3, is the decomposi-tion of the boundary surface into a subset with Dirichletboundary condition on displacement component u i andNeumann boundary condition on traction component T i ,respectively. Finally, Υ is a smooth boundary edge onwhich a jump in higher-order stress traction may arise.The fourth-order nature of the governing partial differ-ential equation emerges on evaluating [14–15] and sub-stituting in [18 ].The Dirichlet boundary condition in [18 ] has the sameform as for conventional elasticity. However, its dualNeumann boundary condition, [18 ] is notably more com-plex than its conventional counterpart, which would haveonly the first term on the left hand-side. Equation [18 ]is the higher-order Neumann boundary condition on thehigher-order stress traction. The physical interpretationof this boundary condition in its homogeneous form isthat the boundary has no mechanism to impose a mo-ment on bonds at the atomic scale. Equation [18 ] spec-ifies that there is no discontinuity of higher-order stresstraction across edges. The weak form of mechanical equi-librium has already been encountered in [16]. For com-pleteness we specify functional spaces: Find u lying in S = { u ∈ H (Ω ) | u i = ¯ u i on Γ u i } such that givenfunctions T i on Γ T i for i = 1 , , P , B de-fined by Equations [14–15], Equation [16] is satisfied forall w lying in V = { w ∈ H (Ω) | w i = 0 on Γ u i } . V. ISOGEOMETRIC ANALYSIS FORHIGH-ORDER PDES
Isogeomeric Analysis (IGA) is a mesh-based numericalmethod with NURBS (Non-Uniform Rational B-Splines)basis functions. The NURBS basis functions have manydesirable properties. They are partitions of unity withcompact support, satisfy affine covariance and providecertain advantages over Lagrange polynomial basis func-tions, which are the mainstay of Galerkin finite elementmethods. These advantages include the imposition ofhigher-order, C n -continuity, positivity, convex hull prop-erties, and being variation diminishing. For a detaileddiscussion of the NURBS basis and IGA, interestedreaders 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 B-spline functions that are defined as follows: Consider pos-itive integers p and n , and a non-decreasing sequence ofvalues χ = [ ξ , ξ , ...., ξ n + p +1 ], where p is the polynomialorder, n is the number of basis functions, the ξ i are co-ordinates 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 startingwith the zeroth order basis functions B i, ( ξ ) = (cid:26) ξ i ≤ ξ < ξ i +1 , p ≥ B i,p ( ξ ) = ξ − ξ i ξ i + p − ξ i B i,p − ( ξ ) + ξ i + p +1 − ξξ i + p +1 − ξ i +1 B i +1 ,p − ( ξ )(20)The knot vector divides the parametric space into inter-vals referred to as knot spans (equivalent to elements inFEM). A B-spline basis function is C ∞ -continuous in-side knot spans and C p − -continuous at the knots. If aninterior knot value repeats, it is referred to as a multi-ple knot. At a knot of multiplicity k , the continuity is C p − k . The boundary value problems in the main textconsider only simple geometries, for which B-spline basisfunctions have sufficed. However, we outline the exten-sion to NURBS basis functions for the sake of complete-ness, noting that the numerical formulation as presentedis valid for any single-patch NURBS geometry. Usinga quadratic B-spline basis (Figure S2), a C -continuousone dimensional NURBS basis can be constructed. N i ( ξ ) = B i, ( ξ ) w i P n b i =1 B i, ( ξ ) w i (21)where w i are the weights associated with each of theB-spline functions. In higher-dimensions, NURBS ba-sis functions are constructed as a tensor product of theone dimensional basis functions: N ij ( ξ, η ) = B i, ( ξ ) B j, ( η ) w ijn b P i =1 n b P j =1 B i, ( ξ ) B j, ( η ) w ij (2D) N ijk ( ξ, η, ζ ) = B i, ( ξ ) B j, ( η ) B k, ( ζ ) w ijkn b P i =1 n b P j =1 n b P k =1 B i, ( ξ ) B j, ( η ) B k, ( ζ ) w ijk (3D)(22) VI. NUMERICAL EXAMPLES
Recall that in two dimensions, the free energy den-sity function depends only on c, e , e and e . In thisfirst communication we wish to focus only on broad phe-nomenology, for which reason we choose the simplest ver-sions of the free energy density function. Apart fromFIG. S2: Figure S2 : Quadratic (p=2) B-spline basisconstructed from the knot vector χ = [0, 0 , 0 , 1/6,1/3, 1/2, 2/3, 5/6, 1, 1, 1].the dependence on c , which is essential for the mechano-chemical spinodal, the coefficients are constant. We haveintroduced the coefficients d c and d e , which control thedepths of the free energy wells and s e , which determinesthe strain minimum in terms of e . Figure 2 of the paperdepicts this construction. We also reduce κ and γ αβ toisotropic forms: κ = κ , γ αβIJ = λ e l δ αβ δ IJ , and set thecomposition gradient-strain gradient coupling coefficients θ α = . While the fullest possible complexity of the cou-pling is not revealed by these simplifications, the aimhere is to present the essential physics that is universalto mechano-chemical spinodal decomposition, postpon-ing details of the more complex couplings and tensorialforms to future communications, where they will be de-rived for specific materials systems. Accordingly, in twodimensions, we have: F ( c, e , e , e ) =16 d c c − d c c + 16 d c c + 2 d e s e ( e + e )+ d e s e e + (2 c −
1) 2 d e s e e (23) G ( ∇ c, ∇ e ) = 12 ∇ c · κ ∇ c + 12 ∇ e · λ e l ∇ e (24)where d c = 2 . d e = 0 . s e = 0 . κ = 10 − and λ e = 10 − . The problem domain is a square of dimen-sions 0 . × . F ( c, e ) =16 d c c − d c c + 16 d c c + 3 d e s e ( e + e + e + e )+ 3 d e s e ( e + e ) + (1 − c ) d e s e e ( e − e )+ (2 c −
1) 3 d e s e ( e + e ) (25) G ( c, e , ∇ c, ∇ e ) = 12 ∇ c · κ ∇ c + 12 ∇ e · λ e l ∇ e + 12 ∇ e · λ e l ∇ e (26)where d c = 2 . d e = 100 . s e = 0 . κ = 10 − and λ e = 1 .
0. The problem domain is a unit cube. The straingradient energy coefficient tensors have been chosen to beisotropic, γ αβIJ = λ e l δ αβ δ IJ , and the three-dimensionalformulation of the main text has been simplified here bysetting the composition gradient-strain gradient couplingcoefficients θ αIJ = 0. VII. LIST OF VIDEOS
FIG. S3:
Movie S1 [MovieS1.mp4]: Evolution ofmicrostructure of the two-dimensional solid withoutward flux corresponding to the plot for l = 0 . Movie S2 [MovieS2.mp4]: Evolution ofmicrostructure of the two-dimensional solid withoutward flux corresponding to the plot for l = 0 . Movie S3 [MovieS3.mp4]: Evolution ofmicrostructure of the two-dimensional solid due toquenching from an initial composition of c = 0 .
46. Thiscorresponds to the plot for l = 0 . Movie S4 [MovieS4.mp4]: Evolution ofmicrostructure of the two-dimensional solid due toquenching from an initial composition of c = 0 .
46. Thiscorresponds to the plot for l = 0 .
01 considered inFigure 5b of the main text.FIG. S7:
Movie S5 [MovieS5.mp4]: Evolution ofmaterial points on the homogeneous free energy densitysurface, F ( c, e , e , e ) for the two-dimensionalformulation, restricted to the { c, e } sub-space. Themicrostructure that forms is similar to that in Movie S1.FIG. S8: Movie S6 [MovieS6.mp4]: Zoomed-in viewwith distorted mesh lines showing the evolution of aninitially square crystal lattice (green) into rectangularvariants (red/blue) as the chemistry-driven cubic torectangular transformation occurs, followed by thesplitting of the variants into a twinned structure.FIG. S9:
Movie S7 [MovieS7.mp4]: Isovolumes of thethree interleaved tetragonal variants and thethree-dimensional microstructure corresponding to theplot for l = 1 . Movie S8 [MovieS8.mp4]: Evolution ofmaterial points on the homogeneous free energy densitysurface, F ( c, e ) for the three-dimensional formulation,restricted to the { c, e , e } sub-space. Themicrostructure that forms is similar to that in Movie S7. FIG. S11: Movie S9 [MovieS9.mp4]: Evolution of across-hatched two-dimensional microstructure.FIG. S12:
Movie S10 [MovieS10.mp4]: Evolution ofplate-like structures in a three-dimensionalmicrostructure. [1] Zuzek, E., Abriata, J. P. and San-Martin, A. and Manch-ester, F. D. The H-Zr (hydrogen-zirconium) system
TheBulletin of Alloy Phase Diagrams , 11:385 (1990).[2] Thomas, J. C., and Van der Ven, A. Finite-temperatureproperties of strongly anharmonic and mechanically un-stable crystal phases from first principles.