A continuum theory for mineral solid solutions undergoing chemo-mechanical processes

Recent studies on metamorphic petrology as well as microstructural observations suggest the influence of mechanical effects upon chemically active metamorphic minerals. Thus, the understanding of such a coupling is crucial to describe the dynamics of geomaterials. In this effort, we derive a thermodynamically-consistent framework to characterize the evolution of chemically active minerals. We model the metamorphic mineral assemblages as a solid-species solution where the species mass transport and chemical reaction drive the stress generation process. The theoretical foundations of the framework rely on modern continuum mechanics, thermodynamics far from equilibrium, and the phase-field model. We treat the mineral solid solution as a continuum body, and following the Larch\'e and Cahn network model, we define displacement and strain fields. Consequently, we obtain a set of coupled chemo-mechanical equations. We use the aforementioned framework to study single minerals as solid solutions during metamorphism. Furthermore, we emphasise the use of the phase-field framework as a promising tool to model complex multi-physics processes in geoscience. Without loss of generality, we use common physical and chemical parameters found in the geoscience literature to portrait a comprehensive view of the underlying physics. Thereby, we carry out 2D and 3D numerical simulations using material parameters for metamorphic minerals to showcase and verify the chemo-mechanical interactions of mineral solid solutions that undergo spinodal decomposition, chemical reactions, and deformation.


Introduction
When considering a deformable medium, chemical reactions may affect the solid's strength and its mechanical properties. Analogously, high mechanical strength may restrict either the volumetric shrinkage or swelling associated with the local volume changes caused by the chemical processes. Therefore, the chemical processes, associated with mass transport and chemical reactions, induce volume changes that lead to stresses around the reaction site.
Finding innovative ways of approaching the modelling of solids is an essential open research topic in science and engineering. For instance, areas such as material science and geoscience are continually searching for new models that allow to improve the properties of materials or to understand the formation of mineral assemblages, which directly relate to solids undergoing chemical processes. In particular, metamorphic petrologists report the reciprocal chemo-mechanical responses of minerals during metamorphism [1][2][3][4][5][6][7]. A variety of study cases of this coupling ranges from grain-scale pressure variations in high-temperature metamorphic rocks to the proper definition of pressure in order to define P/T conditions of mineral assemblages during a metamorphic cycle. Without loss of generality, we use the aforementioned framework to study the tempo-spatial variations of stress-assisted volume changes.
The description of solidity and its properties is crucial to describe the physical and chemical responses of solids accurately. Gibbs' comprehensive study set the foundations for the thermodynamical properties of solids [8]. However, Gibbs' solid model does not quantify the internal adjustment caused by compositional changes since the solid-state diffusion concept did not exist in his time. Herein, we seek to model multicomponent elastic solids that allow for changes in composition while remaining in the solid-state, and particularly, the impact of compositional changes on stress generation [9][10][11][12]. Larché and Cahn introduced the equilibrium conditions for deformable bodies, which change composition as a result of chemical processes [13][14][15]. For instance, dissolution and precipitation at solid-fluid interfaces change the chemical composition of the solid, which in turn induce stresses associated with volume changes. Larche-Cahn's approach models the solid as a network, which allows us to define the stress-strain relations. A solid network can be, for example, the unit cell of the crystalline structure of a mineral, which arranges the atoms in a systematic and repeating pattern. Thus, the network model of Larché and Cahn adequately describes a multicomponent solid.
The outline of this work is as follows. Section §2, we present a coupled system of partial differential equations to model mass transport and deformation processes. Finally, section §3 covers the dimensionless forms of the coupled system of chemo-mechanical equations for the multicomponent framework in conjunction with the dimensionless numbers of the coupling. Additionally, appendix §A covers a thermodynamicallyconsistent treatment to the chemo-mechanical responses of the mineral solid solution.

Kinematics of Motion
We introduce a system of partial differential equations to capture the evolution of a multicomponent elastic solid undergoing spinodal decomposition under reversible chemical reactions. In our framework, the deformations induced across the solid boundaries and compositional changes drive the stress generation process. Henceforth, we refer to this mechanism as stress-assisted volume changes. We treat the solid as a continuum body that occupies an open subset B of the Euclidean space E. A time-dependent deformation field χ : B x T → B t ⊂ E describes the motion from a configuration B onto another configuration B t . We refer to B as the reference configuration and to X as the particles in B. We set the reference configuration such that B represents an undeformed state of the solid. The deformation field characterizes the kinematics of the motion of the particles in the body, and after deformation, it assigns to each material particle X at a given t ∈ T a spatial particle x in the current configuration B t . Then, we express the deformation field as and abusing notation B t = χ t (B).
The deformation field is invertible; namely, there exists an inverse deformation field χ −1 :

Constitutive Equations
Following A, we derive a set of coupled chemo-mechanical equations. By using the Coleman-Noll procedure [16], we find a set of constitutive equations as a pair for each kinematic process. We then rewrite (73) following 74 as Thereby, we find arbitrary values forḞ,φ α R , ∇φ α R , and ∇µ α R at a given time and position such that 5 always hold.
Following the notation and the definition for the Larché-Cahn derivatives for both scalar and gradient fields as proposed by [17], the relative chemical potential µ α Rσ results from the Larcé-Cahn derivative as a consequence of incorporating the mass constraint given by (52). According to Larché-Cahn [15], the relative chemical potential expresses the chemical potential of α-th species measured relative to the chemical potential of σ-th species. This definition entails that, for saturated systems, the mass constrain given by (52) must always hold. Analogously, the relative microforce ξ α Rσ emerges from the constraint imposed in the concentration gradients by (53). As a consequence, we rewrite (5) in the Larché-Cahn sense the following terms: π α := π α σ , µ α := µ α Rσ and ξ α R := ξ α Rσ as well as the material mass fluxes  α R :=  α Rσ as all these quantities are expressed relative to the σ-th reference species. Thus, the energy imbalance is The latter implies that the following relations must hold to keep consistency with the dissipation imbalance We use a logarithmic multi-well potential together with a multi-gradient-type potential for the chemical energy, that is, This expression corresponds to the extension of the Cahn-Hilliard equation towards multicomponent systems [18,19]. The Ginzburg-Landau free energy governs the dynamics of the phase separation process undergoing spinodal decomposition. In (8), N v is the number of molecules per unit volume, k B is the Boltzmann constant, and Ω αβ represents the interaction energy between the mass fraction of the α-th and β-th species. The interaction between the species α over β 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 α = β [18][19][20]. Furthermore, Γ αβ = σ αβ αβ [force] (no summation implied by the repeated indexes) represents the magnitude of the interfacial energy between the α-th and β-th species. The parameters σ αβ and l αβ are the interfacial tension [force/length] and the interfacial thickness 1 for each pair of species (between the α-th and β-th species) [length], respectively. In [18], the authors define the force Γ αβ as N v Ω αβ ( αβ ) 2 .
Following [11], we assume the elastic solid behaves as a compressible neo-Hookean material whose elastic energy is given byψ where G and β are material parameters that relate the shear modulus and the weak compressibility of the material. β is a function of the Poisson ratio ν such that β = 2ν/1 − 2ν. In line with treatments of thermoelasticity, we assume a multiplicative decomposition of the deformation gradient [11], that is, This expression suggests that as the local species concentrations change relative to their initial distribution, the solid must undergo elastic deformation. Moreover, the swelling material parameter ω α is related to the crystalline structure of the solid and its mechanical properties. The evolution of the conserved field ϕ α R obeys a non-Fickian diffusion driven by the chemical potential differences between the species. We combine (7b) and (7c) using the balance of microforces (63) and the constitutive relation for the free energy (74) to express the relative chemical potential of the α-th species as and therefore, where We define p := − 1 3 tr[T R F ] as the mechanical pressure and emphasize that this pressure modifies the mass transport rate. Therefore, for deformable bodies undergoing mass transport, this physical quantity alters the driving force of the chemical process. Furthermore, we define p α ϕ := ω α σ J −1 ϕ p as a mechanical pressure scaled by the local variation of α-th species. Thus, we rewrite (12) as For convenience, we split the chemical potential µ α Rσ such that µ α Rσ = µ α ϕ + µ α s + p α ϕ . Thereby, The constitutive relation for the first Piola-Kirchhoff stress tensor is We also consider the off-diagonal terms in the Onsager reciprocal relations and thus, we describe the species fluxes as where M αβ are the Onsager mobility coefficients. 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 αβ ϕ α R (δ αβ − ϕ β R )I (no summation implied by the repeated indexes) where δ αβ and M αβ are the Kronecker delta of dimension n and the mobility coefficients [19], respectively.

Dimensionless Form of the Chemo-Mechanical Equations
Following A, the coupled system of chemo-mechanical equations readṡ subject to periodic boundary conditions.
Moreover, we introduce a free energy density To make dimensionless the energy densities, and the governing and constitutive equations, we define the following dimensionless variables where u 0 , D 0 , and l 0 account for a reference deformation state, the diffusion coefficient, and the interface thickness of a reference species, respectively. We propose the following sets of scalar and vector dimensionless numbers for the multicomponent chemo-mechanical system, that is, By inserting the dimensionless quantities in (8) and (9), we find the following dimensionless forms of the chemical energy and the mechanical energy 5 where F e = J −1/3 ϕ (I + l∇u). Likewise, the dimensionless forms of the governing and constitutive equations read subjected to the boundary conditions (19).

Numerical Simulations
In this section, we present 2D and 3D simulations of the temporal evolution of single materials modelled as a solid solution composed of three phases A 1 , A 2 , and A 3 to investigate their coupled chemo-mechanical interactions. In particular, the 2D simulation shows how interfacial interactions, together with a reversible chemical reaction between the phases, engender volumetric stresses as a result of local volume changes. The 3D simulation, on the other hand, studies a ripening mechanism. The interfacial interactions between the phases drive the phase separation process and allow for the Ostwald ripening and Gibbs-Thomson effects. We show the dimensionless temporal evolution of the dimensionless phases concentrations as well as the dimensionless displacements in each coordinate direction. By doing so, we seek to understand the interleaving between the physical and the chemical responses of the mineral solid solutions.

Reversible Chemical Reaction of Random Distributed Phases
The reversible chemical reaction between the phases is where the stoichiometry vectors υ αβ and αβ are given by We study the stress-assisted volume changes triggered by the chemical processes. Therefore, we do not drive the deformation by distorting the analysis domain. We set the external body forces and external microforces such that b = 0 and γ α = 0, respectively. Moreover, we neglect all inertial effects. Consequently, the spatial velocity is v = 0. We assume that the three phases diffuse at the same rate. Therefore, we only consider one diffusion coefficient. In this simulation example, we set k + > k − . Thus, the forward chemical reaction occurs faster than the backward one. Our initial condition serves as the reference configuration, which we choose as an undeformed state of the body. The mass supply of each phase, captured by the reaction term s α , results solely from internal contributions as the system (26) reacts. The initial spatial distribution of the phase concentrations is random such that ϕ α takes values between ϕ α ± 0.05 where we assume ϕ α is 1/n. We calculate the concentration of A 3 following the mass constraint given by (52). As mentioned before, we apply this mass constraint when we compute the relative quantities resulting from the Larché-Cahn derivative.
By doing so, we guarantee the consistency of the process. Furthermore, there is no mass flux at the solid boundaries. Figure 1 (a) shows the spatial distribution of the initial phases concentrations in conjunction with the initial displacements. We set the parameters in the chemical energy such that we obtain a triple-well function. This function allows us to model the phase separation process. The reactive system seeks to minimize its global free energy. Thus, the reaction drives the given initial concentrations for the phases A 1 , A 2 , and A 3 towards the concentrations at the well points. Table 1 summarises the parameters used to build up the dimensionless numbers as outlined in (22). Hence,the diffusion matrix for each entry α and β as well as the dimensionless numbers are given by where we choose D 0 = D and 0 = as the reference diffusion coefficient and interface thickness of a reference phase, respectively. The final system of coupled chemo-mechanical equations iṡ In (29), we use the Larché-Cahn derivative with A 3 as the reference phase. We solve the system of partial differential equations in its primal form (30) and (31). We state the problem as follows: find {ϕ, u} ∈ C 2 (P) such that (25) given (19) subject to periodic boundary conditions up to the second derivative of ϕ, and u with respect to X in a square open region B = (0, 1) × (0, 1). We use PetIGA [21], a high-performance isogeometric analysis framework built on top of PETSc [22]. We use a mesh with 64 × 64 elements of a polynomial degree 2 and continuity 1.
We denote H 2 as the Sobolev space of square integrable functions with square integrable first and second derivatives and (·, ·) P as the L 2 inner product over an arbitrary material part P with boundary S. We multiply the Lagrangian version of the phases mass balance (49) by a test function α , which belongs to H 2 , using the definition for the material mass fluxes (17) and integrating by parts, the primal variational formulation can then be given by: Furthermore, the weak formulation of the Lagrangian version of the linear momenta balance reads where we multiply (58) by a test function w i . At early stages Figure 1 (b), t < 4.04 × 10 −6 , the solution goes through an initial spinodal decomposition. This spontaneous phase separation process occurs due to ϑ αβ c > ϑ. Otherwise, the mixture would only diffuse without unmixing. As the phases A i , i = 1, 2, 3, diffuse as a result of their separation, the solid undergoes elastic deformation due to the mass transport. Analogously, the pressure p α ϕ alters the rate at which the phases diffuse. The deformation arises solely from mass transport since the reversible chemical reaction has no significant impact. From Figure 3, we verify that the phases masses do not change substantially in the range 0 < t < 4.04 × 10 −6 . As a consequence, there is no nucleation and growth of phases. With regards to the interfacial energies, in the range 0 < t < 4.04 × 10 −6 the energies decrease gently up to a point, t 1.0 × 10 −7 , where the interfacial energies become constant. A small change in the interfacial energies means that either the phase separation has not evolved significantly or there is no substantial coarsening. Furthermore, at the early stages, the deformation is small as the displacements (see Figure 1) are not large since the phase separation has not evolved such that the phases are unmixed. As expected, the displacements u x and u y in the solid move following the mass transport.
Later on, in the range between 4.04 × 10 −6 < t < 9.03 × 10 −6 , the phase separation becomes prominent as it allows to form spatial domains rich in each component (see Figure 1 (c)). In particular, the phase A 1 remains partially unmixed as there are no rounded inclusions with large concentration (see Figure 1 (c)). On the contrary, for phases A 2 and A 3 , rounded inclusions with large concentrations appear. The deformation results from the mass transport itself as there is no substantial influence of the reversible chemical reaction (see Figure 3). Larger displacements are collocated with the larger inclusions for A 2 and A 3 . This collocation is a consequence of the mass flux towards these points as the inclusions grow. The enlargement associated with the growth of the inclusions induces deformation. F ϕ captures this behavior. The interfacial energies remain roughly constant, which implies that there is a balance between both the creating and disassemble of phases interfaces. As expected, in the time interval between 0 < t < 9.03 × 10 −6 , the tendency of the global free energy is monotonically decreasing as the system goes to a steady state of maximum entropy (see Figure 4). This free energy encompasses the contribution from both the chemical and the mechanical energies. At the early stages, as the evolution towards a steady-state goes on, the system favors phase separation.
For instance, Figure 1 (d) shows the evolution of small A 1 inclusions in the range between 9.03 × 10 −6 < t < 1.26 × 10 −5 . As suggested before, there is mass flux towards these points that allows the inclusion to grow. Consequently, there must be deformation associated with mass transport. Figure 1 (d) also shows the larger displacements in the regions where the inclusions are. On the other hand, A 3 inclusions are large and close enough to start merging. This phenomenon is associated with the minimization of the global free energy as the system reduces its interfacial energies. Nevertheless, when considering the system undergoing chemical reactions, the interfacial energies evolve according to the chemical reaction, which in this modeling example corresponds to a reversible chemical reaction. The growth of more A 3 as a result of the reversible chemical reaction creates more A 3 interface. Figure 3 shows an increase in A 3 interfacial energy as well as its mass. On the other hand, the decomposition of A 3 into A 1 and A 2 following (26) must increase A 1 and A 2 masses, and their interfacial energies. However, this behavior does not last since the rate of creation of A 3 is faster than the decomposition into A 1 and A 2 . This is due to k + >> k − .
From t = 5.64 × 10 −5 , the system shows the merging of large inclusions and the action of the reversible chemical reaction (see Figure 2 (a)). At this stage, the creation of A 3 predominates. One can verify such an assertion by checking the masses. However, there is no interfacial energy growth since the larger inclusions are merging. Therefore, the interfacial energies decrease. Moreover, the A 1 , A 2 , and A 3 inclusions pass from rounded to square-like structures. Such behavior results from the dependency of the chemical potential upon the pressure p α ϕ caused by deformation and the Gibbs-Thomson effect associated with the curvature ∆ϕ α R . Along the boundaries of the inclusions, the driving force of the mass transport changes, which may generate inclusions of asymmetric morphologies. Moreover, the creation of A 1 , A 2 , and A 3 following (26) engenders a mechanical pressure associated with nucleation and growth. Figure 2 (a) also shows the large displacements at the phase boundaries that account for the impact of the chemical reaction on the stress generation (pressure). Figure 2 (b) shows the state at t = 4.26 × 10 −4 . The system forms a chain-like structure composed of the phases A 1 and A 2 , which is surrounded by the phase A 3 . This structure emerges as a result of the merging processes and the reversible chemical reaction. The reversible chemical reaction continues to take at the boundary between phases A 1 and A 2 . Moreover, the phase A 3 decomposes into A 1 and A 2 . At this point in the evolution, the phase A 3 composes almost the whole solid due to k + >> k − . The masses and interfacial energies for phases A 1 and A 2 decrease (see Figure 3). However, for phase A 3 , the mass increases while reducing its interfacial energy (see Figure 3). The phases are unmixed, whereby their concentrations correspond to the concentrations at the well points in the triple-well function. The displacements u x and u y are in the range of the previous stages. However, they move as the phases diffuse as a result of the relation between mass transport and deformation. Figure 2 (b) depicts the larger displacements are in line with the chain-like structure. Figure 2 (c) shows the evolution at t = 5.72 × 10 −4 . The minimization of the global free energy as the system goes to the steady-state reduces the thickness of the chain-like structure. Eventually, the chain is composed of interleaved inclusions of phases A 1 and A 2 . As the inclusions of A 1 and A 2 become smaller, the phase A 3 encloses A 1 and A 2 . The interfacial energies and masses keep decreasing since the action of the reversible chemical reaction has not ceased (see Figure 3). The displacement field shows the interaction between the mechanical and chemical processes. We observe smaller displacements inside the chain-like structure and larger ones outside this region.
In the time interval between 5.72 × 10 −4 < t < 1.84 × 10 −3 , the inclusions of phases A 1 and A 2 shrink as the reversible reaction continues (see Figure 2 (d)). When the reversible chemical reaction ceases, the inclusions are rounded. The chemical reaction is active between 1 × 10 −5 < t < ×10 −2 . The larger (smaller) displacements appear around the inclusions of phase A 1 (A 2 ). These phases partially consume while the phase A 3 gains mass. Finally, rounded structures composed of the three phases diffuse in the solid to generate displacements associated with mass transport. Figure 2 (e) portrays such a behaviour. Along with the whole evolution, free energy always behaves monotonically decreasing (see Figure 4).

Ripening of Spherical Inclusions
We carry out a numerical simulation of a 3D configuration of three spherical inclusions. The spherical inclusions are composed of phases A 1 and A 2 while phase A 3 serves as an interstitial phase. We study the stress-assisted volume changes triggered by the mass transport of the spherical inclusions associated with interfacial effects. We expect Ostwald ripening as a result of the differences in the size of the inclusions. We do not consider external contributions from body forces and external microforces. Consequently, we set b = 0 and γ α = 0. Regarding the kinematics of the motion, we set v = 0 as we do not take into account inertial effects during these quasi-steady processes. We do not allow for chemical reactions between the phases, and therefore, the reactions rates k + and k − are zero. The latter entails that s α = 0.0. Hence, the stresses emerge solely from the mass transport associated with the interfacial interactions between the phases. Without loss of generality, the initial condition serves as the reference configuration. We choose this reference state as an undeformed configuration of the body. The initial and boundary conditions are given by in P, subject to periodic boundary conditionson ∂P × (0, T ). (32) As in the previous examples, the range of the chemical and physical parameters are representative processes in geosciences. We use the same parameters as in Table 1 in §4.1. Nevertheless, in this example we set the reaction rates k + = 0 and k − = 0 as mentioned above. The three phases diffuse with the same diffusion coefficient. Furthermore, the dimensionless numbers correspond to (28). As mentioned before, the phase A 3 serves as an interstitial phase following the mass constraint given by (52). Figures 5 (a) show the initial condition for the phases distribution and displacements, respectively. The system of equations to solve is given byφ where we use the phase A 3 as the reference species. We solve the system of partial differential equation (33) in its primal form (30) and (31). We state the problem as follows: find {ϕ, u} ∈ C 2 (P) such that (25) given (19) subject to periodic boundary conditions up to the second derivative of ϕ, and u with respect to X in a cubic open region P = (0, 1) × (0, 1) × (0, 1). We use the PetIGA [21] to solve the 64 3 element mesh of a polynomial degree 2 and continuity 1. At early stages Figure 5 (b), t < 1.182×10 −6 , part of the phase A 1 deposits on the surface of the inclusion A 2 as a new spherical inclusion A 1 of same radius appears. As the evolution proceeds, between the time range 1.182 × 10 −6 < t < 2.79 × 10 −6 , the inclusion of phase A 2 becomes smaller as its mass goes into the solution. There is no deposition of phase A 2 at this point in the evolution. Regarding the deformation, the stresses associated with the volume changes are small since the displacements do not change substantially (see Figure 5 (c)). Nevertheless, after t > 5.476 × 10 −6 , rings composed of phases A 1 and A 2 appear around A 1 . This occurs as the solution gets supersaturated and the mass of phases A 1 and A 2 migrate to the surface of the more energetically stable structures in the system, which in this case correspond to the spherical inclusions (see Figure 5 (d)). Such mass transport induces volumetric stresses and concomitant displacements around the spherical inclusions (see Figure 5 (d)). As the system tries to minimise its free energy, the masses of phases A 1 , A 2 , and A 3 in the solution separate and merge to form new spherical structures. The phase A 2 locates around the spherical inclusions of A 1 . As mentioned before, the fact that the phases are diffusing induces volumetric stresses. Consequently, we see displacements where the phases are separating and merging (see Figure 5 (d)). Later on, t > 5.626 × 10 −3 , the spherical inclusions composed of phase A 1 and A 2 merge to form a more energetically stable distribution of elongated structures (see Figure 5 (e)). Finally, the phase A 2 wraps the phase A 1 , and A 3 acts as an interstitial phase. The structure at steady state emerges as a result of the coupled chemo-mechanical interactions of the three-component system where the source of stress generation solely results from the mass transport of the phases (see Figure 5 (e)). We develop a thermodynamically-consistent model that describes the evolution of chemically active mineral solid solutions. The theoretical foundations of the framework rely on modern continuum mechanics, thermodynamics far from equilibrium, and the phase-field model, which allow us to derive a set of coupled chemo-mechanical partial differential equations. Using 2D and 3D numerical simulations, we show the evolution of the stress-assisted volume changes triggered by mass transport and chemical reactions. We also include a constraint system using the Larché-Cahn derivative, which takes into account the mass constraint imposed by the solid crystalline structure. Finally, we conclude that the influence of mechanical effects upon chemically active mineral solid solutions must be taken into account and most importantly, they play a substantial role in the evolution of these geosystems.

Acknowledgments
We are indebted to Professor Eliot Fried. We had many exhaustive discussions in which he gave us valuable ideas, constructive comments, and encouragement.

A.1 Measure of Strain
In deforming bodies undergoing mass transport and chemical reactions, the particles move relative to each other as a result of external forces and compositional changes. A description of this movement measures the relative displacement of the particles. We use a Lagrangian description of the displacement field u which defines the kinematics of the motion, that is, and the deformation gradient F = ∇χ t = ∇u + I, where I denotes the second-order identity tensor. To ensure an admissible deformation, that is, a continuum body cannot penetrate itself, the Jacobian of the deformation gradient must fulfill the following constraint The velocity of a material particle X as a function of the motion is and its counterpart in the current configuration is v def = ∂χ(X, t) ∂t In (38), we describe the velocity at time t of a material particle located at x = χ t (X).
Given the definition of the deformation gradient and the spatial velocity, the right Cauchy-Green stress, the Green-Lagrange strain, the rate of strain tensors and the spatial velocity gradient are given by We apply the change of variables theorem to relate the reference and current configurations. In particular, these equations account for the deformation of an infinitesimal line, area, and volume element, that is,

A.2 Fundamental Balances
We derive a set of balance equations in the form of partial differential equations that define how the mass, linear and angular momenta, internal energy, and entropy vary in time as the solid-species system endures mechanical and chemical processes. As suggested in [9][10][11][12], three primary fields govern the coupled chemomechanical responses of the solid: the deformation χ(X, t), the species concentration ϕ α R (X, t) per unit of reference volume, and the chemical potential µ α R (X, t) per unit of reference volume where α denotes the α-th species that composes the solid.
Let P ⊂ B be an arbitrary control volume in conjunction with its boundary S = ∂P. Analogously, consider P t as a bounded control volume of B t such that P t = χ(P) with boundary S = ∂P t . According to Cauchy's theorem, the traction t on a surface da ⊂ S and whose normal n points outwards is t = T(x, t)n. This traction characterizes the force exerted by the rest of the body B t \ P t on P t through da ⊂ S [10,11]. The traction t depends linearly pointwise on the normal n through Cauchy's stress tensor T [23]. Applying the Equation (44) to the identity t R da R = t da, we find the force acting on the surface element da as a function of the surface element da R [9,10]. This identity leads to the nominal stress tensor T R , that is, the first Piola-Kirchhoff, As mentioned above, the chemo-mechanical interactions take place through an elastically deforming solid composed by a network and constituent species. Consequently, we formulate balances of mass conservation for both the solid and the constituent species. Thus, we define ϕ α R as the local concentration of the α-th species per unit of undeformed configuration together with a spatial species outflux  α . In agreement with the balance of mass conservation, the rate of mass change of the α-th species in the control volume P has to be equal to the contribution from the mass supply, typically caused by chemical reactions between the species, and the net mass flux through the boundary S, that is, where s α is the mass supply expressed in the reference configuration. The mass supply is composed of two terms, an external contribution due to external agents and internal contributions caused by chemical reactions. Thereby, Using the divergence theorem, we transform the surface integral of the species flux into a volume integral of the divergence of the species flux as followṡ 17 The Lagrangian description of Equation (49) iṡ where we use the Piola transform. Thus, the material species flux is then  α R = F −1 (J α ). Finally, the localized version of Equation (50) is given byφ The concentration of each species is linearly dependent on the other, via the following constraint, where n stands for the total number of species. The mass constraint that Equation (52) expresses must hold when the solid is solely composed of the diffusing species. Herein, we restrict our attention to the case where mass transport by vacancies is not feasible. Henceforth, a superimposed dot (˙) stands for the material time derivative, for instance,φ α R is the material time derivative of the concentration species. Given the conservation of the solid mass, we define ρ and ρ 0 as the solid density in the current and reference configuration, respectively. Then, the balance of solid mass reads In (54), we convert the volume integral in the current configuration into its counterpart in the reference configuration by employing the relation (45). Finally, we use the localization theorem that leads to the local conservation of solid mass ρ 0 = Jρ.
Neglecting all inertial effects; namely, we assume the spatial velocity v is invariant through the time, the balance of conservation of linear momentum reads The balance of linear momentum relates forces to changes in the motion of the body. Such balance involves the traction t acting on a surface element da as well as a body force b. Conventionally, the body force b accounts for forces resulting from gravitational effects. Through the divergence theorem and the definition (??), we express the surface integral in (56) as a volume integral and after some straightforward manipulations in (57), the localized Lagrangian form of the balance of linear momentum is Div The balance of conservation of angular momentum is After using the definition of the balance of linear momentum, the divergence theorem, and the localization theorem, we end up with T = T . The previous relation implies the symmetry of the Cauchy's tensor [24,25]. Finally, the localized Lagrangian form of the balance of angular momenta is Following the line of thought introduced by Gurtin and Fried [26][27][28][29], we separate balances of conservation laws from constitutive equations. As a consequence, we include a balance of microforces, that is where the vector ξ α and the scalar π α (γ α ) correspond to the α-th microstress and α-th the internal (external) microforce, respectively. In general, the microstresses and microforces are quantities associated with microscopic configurations of atoms. We express the balance of microforces in a Lagrangian form and after applying the locatization theorem, the balance for microforces reads where ξ α R = F −1 (Jξ α ).

A.3 Laws of Thermodynamics and Free Energy Inequality
We separate conservation statements from constitutive equations as suggested by . To describe the thermodynamics of this system, we introduce a power expenditure W ext = W ext (P) + W ext (P) externally to P and P done by the external microforce and force on P, and the microtraction and traction on S W ext (P) = n α=1   By neglecting all inertial effects and body forces, we use the first law of thermodynamics to characterize the balance between the rate of internal energyε and the expenditure rate of the chemo-mechanical power, caused by external forces, species transport, and chemical reactions. The first law is then, There is no contribution of s α int to the energy balance (65). The entropy imbalance, in the form of the Clausius-Duhem inequality, states that the rate of growth of the entropy η is at least as large as the entropy flux q/ϑ plus the contribution from the entropy supply q/ϑ, that is, where q, r, and ϑ stand for the spatial heat flux, heat supply, and temperature, respectively. The localized Lagrangian version of (65) and (66) reaḋ andη ≥ −Div ϑ −1 q R + ϑ −1 r, where q R = F −1 (Jq) is the material heat flux. Moreover, W ext is Rewriting (67) and (68), and multiplying (68) by ϑ, we obtaiṅ and ϑη ≥ ϑ −1 ∇ϑ · q R − Div q R + r.
The Helmholtz free energy results from applying the Legendre transform to the internal energy while replacing the entropy of the system by the temperature as an independent variable., i.e.,ψ =ε−θη−ϑη. Consequently, we obtainψ Introducing the balances of both mass conservation and microforces into (72), the free energy imbalance under isothermal conditions iṡ

A.4 The Principle of Material Frame Indifference
Throughout the derivation of the constitutive behavior of the multicomponent solid, we use the Larché-Cahn derivative for both scalar and gradient fields as expressed by [17], together with the mass constraint given by (52). We assume the following constitutive dependence of the free energy ψ ψ =ψ(ϕ R , ∇ϕ R , F) =ψ ch (ϕ R , ∇ϕ R ) +ψ el (F e (F, ϕ R )).
The objectivity principle requires the constitutive relation (74) to be invariant under a superposed rigid body motion or equivalently, independent of the observer. We can relate two different displacement fields χ and χ * as follows χ * (X, t) = Q(t)χ(X, t) + c(t), where Q(t) represents a rotation tensor and c(t) the relative translations. Therefore, the transformation of the potential (74) following (75) implies which ensures consistency with the dissipation inequality (73) and the principle of frame-indifference.