Extended Larch\'e--Cahn framework for reactive Cahn--Hilliard multicomponent systems

At high temperature and pressure, solid diffusion and chemical reactions between rock minerals lead to phase transformations. Chemical transport during uphill diffusion causes phase separation, that is, spinodal decomposition. Thus, to describe the coarsening kinetics of the exsolution microstructure, we derive a thermodynamically consistent continuum theory for the multicomponent Cahn--Hilliard equations while accounting for multiple chemical reactions and neglecting deformations. Our approach considers multiple balances of microforces augmented by multiple constituent content balance equations within an extended Larch\'e--Cahn framework. As for the Larch\'e--Cahn framework, we incorporate into the theory the Larch\'e--Cahn derivatives with respect to the phase fields and their gradients. We also explain the implications of the resulting constrained gradients of the phase fields in the form of the gradient energy coefficients. Moreover, we derive a configurational balance that includes all the associated configurational fields in agreement with the Larch\'e--Cahn framework. We study phase separation in a three-component system whose microstructural evolution depends upon the reaction-diffusion interactions and to analyze the underlying configurational fields. This simulation portrays the interleaving between the reaction and diffusion processes and how the configurational tractions drive the motion of interfaces.


Introduction
Deep in the Earth, both high temperature and pressure allow for solid diffusion and chemical reactions between rock minerals, which in turn, lead to phase transformations and induced deformation. Significantly, the transport of chemical constituents during uphill diffusion generates phase separation processes as a result of the geothermal gradient in the crust. The phases that compose these solid solutions of minerals diffuse at different rates, and when considering changes in temperature as a result of uphill diffusion, for instance during cooling, phase separation processes such as spinodal decomposition occur. For example, ternary feldspars formed by orthoclase, anorthite, and albite show spinodal decomposition during cooling. Thus, such a process controls the coarsening kinetics of the exsolution microstructure [1][2][3]. Rocks are complex systems composed of several minerals, grain-boundaries, fractures, and pore space where the chemical and mechanical properties may vary in each direction. We describe each mineral as a component of a solid solution; this interpretation of the mixture allows us to explain the coupled reactive spinodal decomposition during exsolution. As a case study, we model a phase merging process driven by interfacial responses coupled with chemical reactions as a first attempt to understand the dynamics of reactive exsolution by spinodal decomposition [4]. We derive a thermodynamically-consistent reactive n species Cahn-Hilliard model that captures the dynamics of such interactions while following in detail the configurational forces that drive this coupled kinematical process.
The multicomponent Cahn-Hilliard model is a useful tool for studying the kinetics of multiphase systems undergoing phase separation. Most importantly, this model tracks the microstructure evolution of the resulting phases to enhance our understanding of the resulting material properties. To describe the underlying physics of this problem, we consider n phase fields representing the concentration of conserved species and use a set of coupled Cahn-Hilliard equations. This representation leads to a system of n degenerate nonlinear fourth-order parabolic partial differential equations. The degeneracy is due to a nonlinear mobility tensor that can vanish depending on the phase field values. We assume that there exist n microforce balances, as similarly proposed by Fried & Gurtin [5,6], and n mass balances accounting for all the relevant chemical reactions, as similarly proposed by Clavijo et al. [7]. We then build an extended Larché-Cahn framework to account for the interdependence between the conserved species. Given the set ϕ = {ϕ 1 , . . . , ϕ n } of species, where n ∈ N, we consider n − 1 independent speciesφ = ϕ \ {ϕ σ } while the σ-th conserved species is used as a reference and determined byφ. Thus, to compute partial derivatives with respect to ϕ α and grad ϕ α , the dependence among species must be taken into account. We then redefine the partial derivative of functions depending upon ϕ and grad ϕ where ϕ σ is constrained to derive the multicomponent Cahn-Hilliard equations. Moreover, in defining these partial derivatives, we arrive at a constrained inner product on a constrained space to appropriately define the gradient energy coefficients Γ αβ .
The outline of this article is as follows. In section §2, we introduce the balances of microforces and augment them with mass balances. In section §3, we present the configurational forces and their balances, and describe how they drive the interface evolution. In section §4, we make the equations dimensionless. Section §5 exemplifies the use of configurational tractions to explain the evolution of a three-alloy mixture. The final section enumerates our conclusions and future work. Appendix A presents the mathematical derivations of this theory.

Theoretical framework
We give a brief overview of the theoretical framework that describes the isothermal evolution of n ≥ 2 reacting and diffusing chemical constituents that occupy a fixed region B of a three-dimensional point space.
2.1. Constituent content balances. We assume that a mass density α , a diffusive flux  α , and a reactive mass supply rate s α characterize the instantaneous state of each constituent α = 1, . . . , n. Also, we require that α ,  α , and s α evolve subject to a pointwise constituent content balance in the forṁ where a superposed dot denotes partial differentiation with respect to time and div denotes the divergence on B. Stipulating that the mass supply rates and the diffusive fluxes satisfy constraints of the form n α=1 s α = 0 and we sum the constituent content balance (1) over α from 1 to n to find that the total mass density = n α=1 α (3) must satisfy˙ = 0 and, thus, is constant.
2.3. Thermocompatible constitutive relations. We introduce constitutive relations for the diffusive flux  α , the reactive mass supply rate s α , the internal microforce density π α , and the microstress ξ α for each constituent α = 1, . . . , n, which allow us to close the system of evolution equations for the phase fields ϕ α , α = 1, . . . , n. These relations must be compatible with the constraints (2) and (5) and with the first and second laws of thermodynamics, which, since we consider only isothermal processes, combine to yield an inequality of the form where ψ is the specific free-energy, and µ α is the chemical potential of constituent α (for details, see Appendix A.2). We define the chemical potential using the Coleman-Noll procedure in the next section.
2.4. Thermodynamical constraints. Throughout the derivation of the constitutive relations for the multicomponent Cahn-Hilliard system, we use the Larché-Cahn derivative (49) from Appendix A.1. Using the Coleman-Noll procedure [9], we find the sufficient conditions to ensure the inequality (8) for arbitrary fields. Thus, a set of paired constitutive equations emerges for each kinematic process. We assume the following constitutive dependency of the free energy ψ within the context of isothermal processes ψ :=ψ(ϕ, grad ϕ), which specializes the free-energy (8) as follows The free-energy imbalance (10) must hold for any arbitraryφ α , gradφ α , and grad µ α fields at a given time and place. Thus, the following relations must hold where M is the mobility tensor, which must be positive semi-definite, that is, n α=1 n β=1 p α · M αβ p β ≥ 0, holds for all p. As (10) expresses all thermodynamically consistent choices, we write the terms π α := π α σ , µ α := µ α σ , and ξ α := ξ α σ relative to the Larché-Cahn construction given their explicit dependence on the Larché-Cahn derivatives as an essential consequence of (5). As a byproduct, we also write the mass flux,  α :=  α σ (x, t; grad µ α σ ), and the surface microtraction, ξ α S := ξ α Sσ (x, t; ξ α σ ) as constructions dependent on the Larché-Cahn derivative. Finally, intrinsically in these definitions, we express all quantities relative to the σ-th species.
Guided by the original Cahn-Hilliard equation [10], we assume that the evolution of the Ginzburg-Landau free energy governs the nature of phase separation undergoing spinodal decomposition. In a multicomponent system, we determine the constitutive relations in (11) from the Ginzburg-Landau free energy expressed aŝ where N v is the total number of molecules of the species α per unit volume, k B is the Boltzmann constant, and Ω αβ represents the interaction energy between the mass fraction of the α-th and β-th species, which is reciprocal; thus, Ω αβ is symmetric. The interaction energy is positive and is related to the critical temperature for each pair of species, ϑ αβ c , (between the α-th and β-th species). Following standard convention, we adopt that Ω αβ = 0 when α = β and Ω αβ = 2k B ϑ αβ c when α = β [10][11][12]. Furthermore, Γ αβ = σ αβ αβ [force] (no sum on α and β) represents the magnitude of the interfacial energy between the α-th and β-th species. The parameters σ αβ and αβ are the interfacial tension [force/length] and the interfacial thickness 1 for each pair of species (between the α-th and β-th species) [length], respectively. Cahn & Hilliard [10] define the force Γ αβ as N v Ω αβ ( αβ ) 2 .
We express the relative chemical potential of the α-th species, in the Larché-Cahn sense, by combining the expressions (11a), (11b), the microforce balance (7), and the constitutive relation for the free energy (12), we arrive at Therefore, the combination of (13) with (12) specializes to In the following, we assume an isotropic mobility M αβ = M αβ 1 being M αβ symmetric, but we consider the off-diagonal terms in the Onsager reciprocal relations. We use the standard assumption that the mobility coefficients depend on the phase composition. In particular, we express this dependency in terms of the concentration of each species. We use the definition M αβ := M αβ 0 ϕ α (δ αβ − ϕ β ) with no summation on α and β and M αβ 0 is the mobility between the α and β species, with dimension of length 4 per unit force and time [11]. Thus, (2) implies the following relation 2.5. Chemical reaction. Let ϕ α be the concentration of a species A α , such that ϕ α := [A α ]. Following Krambeck [13], we express the c-th chemical reaction, in a set of n s chemical reactions, n s ∈ N, as where υ c,α r and υ c,α p are the α-th stoichiometric coefficient of the c-th chemical reaction of the reactants and products, respectively. The number of non-zero stoichiometric coefficients υ c,α r and υ c,α p define the number of reactants n c r and products n c p in the c-th chemical reaction. For the c-th chemical reaction, zeros populate υ c,α r for α > n c r , whereas zeros populate υ c,α p for α ≤ n c r . k c + (k c − ) denotes the c-th forward (backward) reaction rate (see, for details Appendix A.3). We focus on ideal materials, then, the c-th rates of both the forward and backward reactions read Finally, the internal rate of mass supply term for all n s chemical reactions that enters in (6) is

Configurational fields
We describe the interfacial evolution, and its thermodynamics using the configurational forces proposed by Gurtin [14], which relate the integrity of the material and the movement of its defects. The configurational forces expend the power associated with the transfer of matter, which allow us to interpret them thermodynamically. Using the configurational balance for a part P by Fried [15], we have which renders, after localization, where C is the configurational stress tensor and f (e) is the internal (external) force.
Following Appendix A.4, we substitute the constitutive relation (78) in the relation (74), allows us to express the configurational stress as We obtain explicit forms for the internal and external configurational forces by combining (11) and (21) with (22), that is By considering the Larché-Cahn derivatives, we express the configurational stress (22) as a configurational stress relative to the σ-th species as follows while is the relative internal configurational force. The external configurational force is not determined using a constitutive relation; thus, it does not depend upon the choice of the reference species.
Remark 1 (Invariance of configurational balance to reference species). Let ϕ σ be the reference species. We establish the following relations for the terms appearing in the configurational stress (24) while for the internal configurational force (25) Analyzing (26) and (27), we conclude that only one term in C σ (24) depends on the reference species. Therefore, the relative configurational stress becomes Meanwhile, we specialize representation of the relative internal configurational force (25) with (28) yielding Finally, although both the configurational stress and the internal configurational force explicitly depend on the choice ϕ σ ; nevertheless, their dependencies cancel each other's contribution to the configurational balance (21),
To make the equations dimensionless, we introduce the reference energy density ψ 0 := 2N v k B ϑ and define the set of diffusion coefficients D αβ , The reference energy density relates the species mobilities with the species diffusion as proposed in [16,17]. We also define the following dimensionless variables Conventionally, the definition of the reference time T 0 for the Cahn-Hilliard system relates the diffusion coefficient, the interfacial thickness, and domain length, that is, where L 0 0 [18,19]. We set D 0 and 0 as the reference diffusion coefficient and interface thickness of a reference species, and introduce the following dimensionless numbers for the multicomponent system   Thus, by inserting the dimensionless quantities in (32), we find the following dimensionless forms in D × (0, T ), with the initial condition (33).

Numerical simulation: merging of circular inclusions
We now simulate the interactions between three species where A 1 and A 2 represent the reactants, while A 3 the reaction products. The inclusions (represented by species 1, A 1 ) are embedded in species 2, A 2 . We express the chemical reaction as which takes place at the interface producing the third species, A 3 . We state the problem as: find ϕ satisfying (37) given (33) subject to periodic boundary conditions up to the fourth derivative of ϕ with respect to x in a square open region D = (0, 1) × (0, 1). We discretize the resulting system of partial differential equations using PetIGA [20], a high-performance isogeometric analysis framework. We solve this system of equations in their primal form using a 256 × 256 element mesh of a polynomial degree 4 and continuity 3. The initial and boundary conditions are subject to periodic boundary conditions on ∂D × (0, T ), and the three subfigures on top of Figure 1 depict this initial condition. Table 1 summarizes the dimensional parameters used to obtain the dimensionless parameters in (40) and (41). The diffusion matrix for each entry α and β reads Next, for clarity, we represent α and β as matrix-columns and -rows indices, which render the remaining dimensionless parameters as follows.
where we choose D 0 = D and 0 = 23 as the reference diffusion coefficient and interface thickness of a reference species, respectively. Here, the configurational tractions drive the interfacial motion in this multicomponent system undergoing reactions. We express the configurational traction along a level curve L α * , upon which ϕ α = ϕ α * . We then introduce the normal and tangential coordinates n α and m α on L α * , with unit vectors ν α and τ α defined such that grad ϕ α = |grad ϕ α |ν α , augmented by a sign convention which ensures that rotating τ α clockwise by π/2 yields ν α . In reckoning the relative configurational stress in a {n α , m α }-frame, we arrive at with ζ := (ψ − n α=1 µ α ϕ α ), see (78) in Appendix A.4. We can now specialize (43) with a free-energy of the formψ which renders the following relative configurational stress Thus, the configurational tractions C σ ν α are In the simulations, we compute the relative physical and chemical quantities, such as the relative chemical potential, mass fluxes, microstresses, and byproducts, by setting the reaction product species, that is, A 3 as the reference phase field. This simulation shows that the configurational fields can describe the behavior of the phase evolution. However, this initial work does not exploit this tool exhaustively nor comprehensively. Figure 1 depicts the merging process of two circular inclusions of distinct size into a single one. The figure spans from the early stages until the merged inclusion becomes stationary. From left to right, we depict phases ϕ 1 , ϕ 2 , and ϕ 3 , while from top to bottom, the evolution of the three phases for the dimensionless times t = 0, 6.00 × 10 −6 , 2.25 × 10 −5 , 7.41 × 10 −5 , and 5.96 × 10 −2 . Figures 2 and 3, respectively, present e 2 · C 3 ν 2 (x 1 = 0.5, x 2 ) and e 2 · C 3 ν 3 (x 1 = 0.5, x 2 ) on the left panel, and e 2 · C 3 ν 2 (D) and e 2 · C 3 ν 3 (D) on the right panel. That is, the left panels display the profile of the relative configurational traction along x 2 , while the right panels display the vertical component of the relative configurational traction on the whole domain. These figures show the x 2 axis in red. From top to bottom, we present these configurational fields at the dimensionless times t = 0, 5.39 × 10 −6 , 6.00 × 10 −6 , 6.57 × 10 −6 , 2.25 × 10 −5 , 7.41 × 10 −5 , and 5.96 × 10 −2 . Figure 2 shows that configurational tractions between the inclusions have opposite directions pushing against one another. As the inclusions approach each other, the configurational traction profiles become antisymmetric in the region where the merging takes place (second plot, from top to bottom, in Figure 2). In this region, the ridge and the valley propagate towards each other until the interfaces merge. Later, the configurational tractions annihilate one another (third plot, from top to bottom, in Figure 2). The third species appears as the chemical reaction takes place. Figure 3 shows how the relative configurational traction C 3 ν 3 pushes apart the boundaries of the double ring, formed by this species. This traction drives the growth of the area encircled by the double ring, which occurs at the expense of the other two species through the chemical reaction. Figure 3 (second plot, from top to bottom) shows the tractions on each ring as they push against each other, which favours merging. At later stages, a single ring-like structure remains, formed by the product species. This ring lies in between the interface formed by the reactant species. Consequently, the process reaches a semblance of a steady-state when the product species obstructs further the chemical reactions. Figure 4 presents a snapshot sequence detailing the merging process from left to right and top to bottom. We use the relative configurational traction C 3 ν 2 , their streamlines (white), L 2 0.5 (black) on top of the phase field ϕ 1 at the dimensionless times t = 0, 5.39 × 10 −6 , 6.00 × 10 −6 , 6.57 × 10 −6 , and 2.25 × 10 −5 to exemplify this evolution. In these snapshots, we show the configurational traction C 3 ν 2 with black arrows. Before the merging occurs, two node sinks arise, see Figure 4b. These node sinks pull the phase ϕ 2 initiating the merging process. Soon after the node sinks are formed, Figure 4c, phase ϕ 2 migrates and leaves a 'bridge' between the inclusions. This 'bridge' is formed by phase ϕ 1 . After, the merging of the level curve L 1 0.5 (black line) occurs, see Figure 4d.   . Node sinks triggering the merging C 3 ν 2 . From left to right, the relative configurational traction C 3 ν 2 , streamlines of this traction (white), L 2 0.5 (black) on top of the phase field ϕ 1 at the dimensionless times t = 0, 5.39 × 10 −6 , 6.00 × 10 −6 , 6.57 × 10 −6 , and 2.25 × 10 −5 .

Final remarks
In this work, we present a continuum framework to model phase separation processes such as spinodal decomposition during cooling as a result of uphill diffusion. These phases are composed of solid solutions of minerals and diffuse at different rates. In this first attempt to model solid diffusion and chemical reactions between rock minerals, we neglect deformation and heat transfer. To this end, we derive a thermodynamically consistent continuum theory for the multicomponent Cahn-Hilliard equations while accounting for multiple chemical reactions. We consider multiple balances of microforces augmented by multiple constituent content balance equations within an extended Larché-Cahn framework. Moreover, we derive a configurational balance that includes all the associated configurational fields in agreement with the Larché-Cahn framework. In a simple simulation, we depict the role of the configurational tractions during the merging process coupled with a chemical reaction. Last, in upcoming works, we plan to model the contributions of deformation in the thermodynamic pressure arising from chemical processes such as mass transport, chemical reactions, and interfacial effects.

Acknowledgments
We are indebted to Professor Eliot Fried. We had many exhaustive discussions in which he gave us valuable ideas, constructive comments, and encouragement. This publication was Constraint (5), with (4), implies that the set of concentrations ϕ must be 0 < ϕ α < 1. If we vary one concentration ϕ α while holding all others fixed violates the constraint (5). Thus, the conventional partial derivative on functions such as F, on which the constraint (5) is active, is not appropriately defined. To overcome this shortcoming, Larché and Cahn [21] defined the following operation in which we choose any two concentrations ϕ α and ϕ σ from the set of variables. Then, we introduce an infinitesimal change in ϕ α , which induces the opposite infinitesimal variation onto ϕ σ , while holding all other variables unchanged. Thus, this definition satisfies (5) by construction while we express the concentration ϕ σ as In multicomponent Cahn-Hilliard systems, we incorporate cross-diffusion gradient energy coefficients Γ αβ into the free-energy definition and obtain the following free-energy densitŷ ψ(ϕ, grad ϕ) := f (ϕ) + n α=1 n β=1 Γ αβ grad ϕ α · grad ϕ β .
Elliott & Garcke in [11] prove that multicomponent systems are well-posed when Γ αβ is positive definite, among other conditions. We show that this condition is sufficient but not necessary. To do so, we extend the ideas of Larché-Cahn and define a constrained inner product on a constrained space. We consider a set of vectors {p α } subject to the following constraint n α=1 p α = 0, and use the following inner product n α=1 n β=1 Γ αβ p α · p β .
Let each entry of Λ αβ be a single number κ. Thus, due to (52), {p α } is in the null space of Λ αβ , that is, Null(Λ αβ ) = {p α }. Similarly, if each row of Λ αβ is given by the same entry κ β , we arrive to the same conclusion. For any of these cases, we have that We impose the constraint (52) with respect to the component σ to the quadratic form (53) to obtain We reinterpret this result as an inner product in an unconstrained space of dimension n − 1 with a noninvertible mapping Γ αβ → Γ αβ σ defined as Consequently, the problem is well-posed if Γ αβ σ is positive definite. Moreover, Γ αβ can be indefinite without compromising the well-posedness of the problem. Now, let Γ αβ be a diagonal matrix such that From (54), we rewrite Γ αβ as where 1 αβ is a constant matrix populated by ones and δ αβ is the Kronecker delta, both of dimension n.
Although the matrix (58) has a null diagonal, the mapping defined by (56) is identical to the one of the diagonal matrix (57) for all vectors that satisfy the constraint (52).
A.2. Thermodynamics. Here, we establish the first and second law of thermodynamics. First, we augment the species balances (6) to consider an external mass supply s α ext as well as an internal one s α arising from chemical reactions. We treat chemical reactions in a similar fashion as Gurtin & Vargas [22]. Moreover, following Gurtin [23] and Cherfils et al. [24], see also [25,26], we separate conservation statements from constitutive equations. Thus, we introduce the external power expenditure W ext to P done by the external microforces on P and microtractions on S to describe the thermodynamics of this system as follows where n is the total number of species and ξ α S = ξ α · n is the α-th microtraction. The first law of thermodynamics states the energy balance between the interleaving of internal energy and the expenditure rate of the chemical (diffusion and reaction) power. The entropy production imbalance, or the second law of thermodynamics in the form of the Clausius-Duhem inequality, states that the rate of growth of the entropy is at least commensurate with the entropy flux and its supply. Thus, we can express these two laws as where ε and η represent the internal-energy and entropy densities, respectively, q is the heat flux, r is the heat supply, and ϑ > 0 is the absolute temperature. There is no contribution of s α to the energy balance (61). Using the external power expenditure (60), the microforce balance (59), and the constituent content balance (59), we can localize the first two laws of thermodynamics (61) to Rewriting (62) 2 , we obtain We now define the free-energy density as which allow us to rewrite the equation system in terms of ϑ and ψ. To employ this transformation, we multiply (63) by ϑ and subtract the result from (62) 1 to express the pointwise free-energy imbalance as Remark 2 (Alternative derivation-Principle of virtual power). The definition of virtual power expenditure encompasses internal and external contributions. Internally to P, the power exerted by internal microforces and the microstresses; while externally to P, the power effected by the external microforces on P and microtractions on S. This definition assumes that these contributions equilibrate each other, that is, where the definitions of the internal and external virtual powers are and where {χ α } is a set of n kinematically admissible fields. Finally, we apply the divergence theorem to (66) and use standard variational arguments to localize the balance of microforces (59) to the following. For a more general approach, see [27].
A.3. Theory of reacting materials. Theoretically, the total number m of possible independent chemical reactions, where m ≥ n s ∈ N, is not arbitrary. We seek to fit our framework in the thermochemistry theory of reacting materials (see, [28] and [29,30]). Thus, we also postulate the indestructibility of the atomic substances n α=1 t αι s α m α = 0, where n a ∈ N is the number of atomic substances making up all the components A, m α is the molecular weight of the α-th component, and t αι is a non-negative integer expressing the number of atoms of the ι-th atomic substance present in the α-th component. This postulate assumes that the atomic substance are indestructible. Moreover, usually t αι is not a square matrix and rank(t αι ) = min(n, n a ). Finally, the maximum number of possible chemical reactions is m := n − rank(t αι ).
In this setting, forward reactions and their reciprocal backward reaction are not independent. Thus, we represent them as a single, effective, chemical reaction.
A.4. Configurational stress and force. We describe the configurational stress, and the internal and external forces arising in multicomponent systems. We first establish how configurational forces expend power in a migrating control volume P . We define q as the migrating boundary velocity acting on S with n being its outward unit normal. We also refer the reader to [15,31,32]. For a migrating volume P the constituent content balance (59) in the partwise form specializes tȯ P ϕ α dv − S ϕ α q · n da = − S  α · n da + P s α dv.
We use the external virtual power (68), where γ α and ξ α S are conjugate toφ α . We set as virtual field the advective termφ α + grad ϕ α · q to follow the motion of S augmented by the fact that the configurational traction Cn is power conjugate to q on S . Since ξ α S (φ α + grad ϕ α · q) = (ξ α · n )φ α + (grad ϕ α ⊗ ξ α )n · q, we arrive at an expression of the total external configurational power The relevant part of the motion of S only involves its normal component q · n . Thus, the power expended is indifferent to the tangential component of q, yielding C + n α=1 grad ϕ α ⊗ ξ α =: ζ1 , where ζ is a scalar field. Thus, the first integral of (73) becomes S ζ q · n da.
The arguments that led to the free-energy imbalance (65), allow us to analyze isothermal processes in a migrating control volume P with a velocity q. Hence, we arrive aṫ