Chemical potentials in non-hydrostatically stressed anisotropic phases

SUMMARY

Chemical potentials can be defined as the partial derivatives of the Helmholtz energy with respect to moles of chemical components under conditions of zero domain strain and fixed temperature. Under hydrostatic conditions, chemical potentials are dependent only on state properties. Under non-hydrostatic conditions, they also depend on a ‘chemical expansivity tensor’—a second-order tensor with unit trace that characterizes how the elastic network of the phase is compressed to accommodate new material within the local domain element. The five degrees of freedom of this tensor generate a class of chemical potentials. An important group within this class are the ‘uniaxial chemical potentials’, which quantify the Helmholtz energy change when new material is incorporated via compression along a single axis. Chemical and mechanical equilibrium is achieved when all uniaxial chemical potentials remain constant along their respective axes.The derived expressions apply to both crystalline and amorphous materials. Their utility is demonstrated through solutions to classic phase-equilibrium problems.

1 INTRODUCTION

Heterogeneous stress fields play a crucial role in metamorphic processes such as dissolution, precipitation and the development of microscale fabrics such as cleavage, schistosity and gneissic augen (Tajčmanová et al. 2014; Hobbs & Ord 2016; Wheeler 2018). A critical aspect of these fabric-forming processes is the transport of one or more chemical components. Consequently, modelling metamorphic processes requires a thorough understanding of the interplay between stress and chemical diffusion.

The mathematical treatment of chemical equilibrium under non-hydrostatic conditions has a long history. In his pioneering work, Gibbs (1878) investigated the equilibrium condition in a system comprising two reacting phases; an inviscid liquid and a one-component solid. This limited scope allowed Gibbs to derive equations for chemical equilibrium within the framework of classical hydrostatic thermodynamics and elastic theory. However, his treatment did not extend to scenarios involving solid–solid contacts, solid solutions, or chemical diffusion. During the subsequent decades, numerous authors have attempted to generalize Gibbs’ theory to systems with arbitrary numbers of chemical components and interactions between solid and liquid phases (e.g., Kamb 1961; Bowen & Wiese 1969; Larché & Cahn 1973, 1985; Robin 1974; McLellan 1980; Truesdell 1984; Mullins & Sekerka 1985; Rusanov 2005; Frolov & Mishin 2010a, b, 2012a, b). These studies share the foundational premise that the total internal energy of an isolated system is minimized at constant entropy and volume.

A central question emerging from this body of work is whether the state property known as the chemical potential can be meaningfully generalized to non-hydrostatic states. Under hydrostatic conditions, the chemical potential for each chemical component k is defined as

with all symbols detailed in Table 1. Chemical equilibrium is achieved when these chemical potentials are uniform throughout the domain. Spatial gradients in chemical potential represent thermodynamic forces that can drive chemical diffusion (e.g. De Groot & Mazur 1984). A non-hydrostatic generalization that retains the properties of its hydrostatic counterpart would be valuable for assessing equilibrium and modelling disequilibrium processes in systems with heterogeneous stresses.

Gibbs himself recognized that, under non-hydrostatic conditions, no single scalar chemical potential can uniquely characterize chemical equilibrium in a three-dimensional domain. In subsequent years, different perspectives have emerged: some researchers have argued that no useful definition of chemical potential exists for non-hydrostatic states (Kamb 1961; Wheeler 2018); others maintain that only relative chemical potentials can be defined (Larché & Cahn 1973, 1985); still others suggest that chemical potentials can always be defined, albeit without necessarily being uniform at equilibrium (Herring 1951; Mullins & Sekerka 1985); and some even propose that chemical potentials may assume tensorial forms (Rusanov 2005). In the studies that do introduce a definition of chemical potentials (e.g. Larché & Cahn 1973; Mullins & Sekerka 1985), a ‘crystallinity’ constraint is imposed and a distinction is drawn between mobile and immobile species. These features are not present in the classical hydrostatic formulation.

In recent years, the issue of chemical equilibrium under heterogeneous stress conditions has received renewed attention, driven by efforts to interpret disequilibrium textures (such as chemical zoning in minerals) and to develop more realistic models of geological processes. The use of different definitions of chemical potentials and thermodynamic driving forces has fuelled debates over whether stress gradients can generate chemical zoning (Powell et al. 2018; Tajčmanová et al. 2021; Hess et al. 2022; Hess & Ague 2023) and whether deviatoric stresses significantly influence chemical equilibria and reactions (Wheeler 2014; Fletcher 2015; Hobbs & Ord 2016; Wheeler 2020; Hess & Ague 2024). These debates highlight the need for a fresh examination of chemical potentials under non-hydrostatic conditions and an evaluation of their utility in solving thermodynamic problems. These are the goals of the current paper.

The paper is organized as follows: Section 2 briefly explains how the non-hydrostatic formulation of the Helmholtz energy at fixed composition reduces to its hydrostatic counterpart. Chemical potentials are introduced in Section 3, first in their conventional hydrostatic form (Section 3.1) and then extended to non-hydrostatic stress states (Section 3.2). In Section 4, these chemical potentials are used to study chemical equilibration under non-hydrostatic conditions, including a classic Gibbs thought experiment and the melting-point depression of a one-component solid. Finally, Section 5 compares the new derivation with existing approaches and discusses several arguments from the literature.

2 NON-HYDROSTATIC STRESSES AT FIXED COMPOSITION

Throughout this paper, we define thermodynamic properties as derivatives of the Helmholtz energy, |$\mathcal {F}$| (⁠|$\mathcal {F} = \mathcal {E} - T S$|⁠), rather than the internal energy as given in eq. (1). The primary advantage of this approach is that, while thermodynamic equilibrium is attained by minimization of the internal energy at constant volume and entropy, it is also attained by minimization of the Helmholtz energy at constant volume and temperature (Callen 1985). By employing the Helmholtz energy, we can avoid calculations involving thermal equilibration and concentrate on the interplay between mechanical and chemical effects. Importantly, both the minimization of the internal energy and that of the Helmholtz energy comply with the Second Law of Thermodynamics (often expressed as the maximization of entropy at fixed energy and volume; see Appendix  A).

Before introducing chemical variability, we first derive an expression for the time derivative of the Helmholtz energy density (⁠|$\dot{\mathcal {F}}^{\rho }$|⁠) for a closed, one-component, purely elastic system under hydrostatic conditions from the more general non-hydrostatic expression:

The first term in eq. (2) is the rate of mechanical work. The tensor |$\dot{\boldsymbol {\varepsilon }}^{\text{ntwk}}$| is the network strain rate, that is, the strain rate of the elastic network of bonds connecting the atoms, ions and vacancies that make up the material. This strain rate is a function of |$\boldsymbol {x}$|⁠, the current network coordinates, and |$\dot{\boldsymbol {x}}$|⁠, the rate of change of those coordinates with respect to time. The network strain rate can be decomposed into isotropic (⁠|$\dot{{\varepsilon }}^{\text{ntwk,iso}}_{ij} = \frac{1}{3} \delta _{ij} \sum _{l,m} \delta _{lm} \dot{\varepsilon }^{\text{ntwk}}_{lm}$|⁠) and deviatoric (⁠|$\dot{\varepsilon }^{\text{ntwk,dev}}_{ij} = \dot{\varepsilon }^{\text{ntwk}}_{ij} - \dot{{\varepsilon }}^{\text{ntwk,iso}}_{ij}$|⁠) parts. Using the small strain relation

we can decompose the work term into volumetric and deviatoric components:

In these expressions, |$V^{\text{ntwk}}_{\text{mol}}$| is shorthand to describe the volume of a mole of units of the elastic network. In a crystalline material, one could equate |$V^{\text{ntwk}}_{\text{mol}}$| with a mole of formula units, a mole of unit cells or a mole of lattice points. In amorphous materials such as glasses, |$V^{\text{ntwk}}_{\text{mol}}$| could correspond to a mole of network nodes rather than lattice points.

Under strictly hydrostatic conditions (⁠|$\boldsymbol {\tau } = 0$|⁠), eq. (5) reduces to the classic total differential for the Helmholtz energy at fixed composition:

To derive this expression, we have also used the relationship between total volume, molar network volume and number of moles of network,

and the constraint that |$\dot{(\ln n^{\text{ntwk}})} = 0$|⁠.

3 CHEMICAL POTENTIALS

3.1 Hydrostatic conditions

Eq. (6) can be extended to open systems with multiple chemical components by adding a dependence on a vector of variables |$\boldsymbol {n}^{\rho }$|⁠, where |$n^{\rho }_k$| represents the number density of moles of chemical component k in the material (Nye et al. 1985). The relationship between |$n^{\text{ntwk}}$| and |$n_k$| can be written:

where |$s_k$| represents the change in |$n^{\text{ntwk}}$| when 1 mole of component k is added to the network. For substitutional solid solutions where each component k has the same network structure as the whole (e.g. |$\mathrm{[Mg]_{2}SiO_{4}}$| and |$\mathrm{[Fe]_{2}SiO_{4}}$| endmembers in |$\mathrm{[Mg,Fe]_{2}SiO_{4}}$| olivine), |$s_k = 1$| for each and every component k. If the component is defined as an exchange component (e.g. |$\mathrm{[Mg_{-1}Fe_{1}]}$|⁠), |$s_k = 0$|⁠.

The exact differential for the Helmholtz energy density under strictly hydrostatic conditions can be written as (Callen 1985, eq. 7.25)

where the chemical potential of component k under hydrostatic conditions, |$\mu ^{\text{hyd}}_k$|⁠, is identical to that in eq. (1) but equivalently defined as:

In these expressions, V is the total volume of the homogeneous domain, that is, if |$\textrm {d}(\ln V) = 0$|⁠, then a positive value of |$\textrm {d}n_k$| will pack more chemical component k into the same volume. As we shall see in Section 3.2, the distinction between domain volume and material volume must also be considered in the generalization to non-hydrostatic conditions.

A Legendre transformation is often used to convert between Helmholtz energy (eq. 9) and a hydrostatic ‘Gibbs energy’ that is minimized by spontaneous processes at fixed P and T (Callen 1985):

It can be shown by change of variables that the hydrostatic chemical potentials can be defined using several other partial derivative expressions in addition to eqs (1) and (10), for example:

In the last expression, |$\bar{\mathcal {F}}_{k}$| and |$\bar{V}_{k}$| are the partial molar Helmholtz energy and volume at fixed network strain and temperature

The material added to the network under these constraints shares the same network structure and deformation as the existing material. If the material added has a different bulk composition to the existing material, chemical diffusion takes place at a local scale to maintain local chemical homogeneity.

3.2 Non-hydrostatic conditions

3.2.1 Chemical strain and the chemical expansivity tensor

Having reviewed the meaning of various chemical terms under hydrostatic conditions, let us now add the vector of compositional variables |$\boldsymbol {n}^{\rho }$| directly to eq. (2):

If we restrict this expression to hydrostatic states, following the steps in Section 2, we obtain the expression:

The total number of moles of network n is now allowed to change, so unlike in Section 2, we cannot replace the change in molar volume (⁠|$\ln V^{\text{ntwk}}_{\text{mol}}$|⁠) with the change in domain volume (⁠|$\ln V$|⁠) via eq. (7). Therefore, although correct, this expression does not reduce to the desired hydrostatic form (eq. 9), that is, we cannot equate |$\bar{\mathcal {F}}_k$| with a chemical potential |$\mu _k$|⁠. To find an expression that does include a generalized chemical potential, we must replace the network strain rate in eq. (18) with a strain rate that satisfies the relation

where V is the domain volume. This expression is satisfied if |$\dot{\boldsymbol{\varepsilon }}$| corresponds to the domain strain rate, which includes components of strain due to network deformation and chemical addition or removal from the domain:

While |$\dot{\boldsymbol {\varepsilon }}^{\text{ntwk}}$|⁠, |$\dot{T}$| and |$\dot{\boldsymbol {n} }^{\rho }$| in eq. (18) are all unambiguously defined independent variables of |$\dot{\mathcal {F}}$|⁠, we have not yet clearly defined |$\dot{\boldsymbol {\varepsilon }}^{\text{chem}}$|⁠. To do this, we first multiply both sides of eq. (21) by |$\delta _{ij}$| and then substitute in eqs (4), (7) and (20):

For a given change in |$n^{\text{ntwk}}$|⁠, this equation does not uniquely define |$\dot{\boldsymbol {\varepsilon }}^{\text{chem}}$|⁠, but it does impose the constraint:

or, by rearranging the first equality in eq. (24) and using eqs (8) and (17):

The ‘chemical expansivity tensor’ |$\boldsymbol {\Xi }$| defined in eq. (24) describes how the addition of material expands a domain element if there is no network strain, or, equivalently, how the addition of material squeezes the elastic network in a domain of fixed shape (the second and last equalities in eq. 24). A visualization of chemical strain, network strain and domain strain is shown in Fig. 1.

3.2.2 Non-hydrostatic chemical potentials and their relationship with state properties

Having introduced the concepts of domain strain and the chemical expansivity tensor, we can now write down an expression for the rate of change of the Helmholtz energy density and chemical potentials |$\boldsymbol{\mu }^{\boldsymbol {\Xi }}$| for a given value of |$\boldsymbol {\Xi }$|⁠:

The chemical expansivity tensor |$\boldsymbol {\Xi }$| is not a state property of the system; it must be defined independently before the chemical potentials |$\mu ^{\boldsymbol {\Xi }}_k$| become well-defined. We can see how the chemical potentials relate to state properties by substituting eqs (21) and (25) into eq. (26):

By comparison with eq. (18):

Rearranging this equation yields expressions for the chemical potentials (using eq. 23):

Under hydrostatic conditions, eq. (31) reduces to eq. (15), as required.

3.2.3 The isotropic strain chemical potential

Substituting |$\Xi _{ij} = \delta _{ij} / 3$| into eq. (31) yields a chemical potential that is both a property of state and also not dependent on deviatoric stress. This ‘isotropic strain chemical potential’ has a simple form:

which is identical to the hydrostatic expression (eq. 15). All other chemical potentials are related to |$\boldsymbol {\mu }^{\boldsymbol {I}/3}$| by the relation:

The isotropic strain chemical potential is particularly useful in the context of generalizing the Gibbs–Duhem Equation (Appendix  B):

4 CHEMICAL EQUILIBRIUM ACROSS A PHASE BOUNDARY

We are now in a position to discuss the relationship between chemical potentials and chemical equilibrium. Imagine a system closed to chemical exchange but open to heat exchange, surrounded by a much larger, isolated temperature bath. Inside the system are two homogeneous subsystems (which might be two different phases), separated by a single planar surface with normal |$\hat{\boldsymbol {n}}$|⁠. Material can cross this interface, and the interface itself can move, but material cannot be introduced through any other surfaces. Under these restrictions, the ‘chemical expansivity tensor’ |$\boldsymbol {\Xi }$| (Section 3.2.1) for both subsystems is uniquely defined:

Minimization of the Helmholtz energy of the system (and therefore maximization of the entropy of the system and its surroundings, Appendix  A) is achieved when the ‘uniaxial chemical potentials’ |$\mu ^{\hat{\boldsymbol {n}}\otimes \hat{\boldsymbol {n}}}_{k}$| corresponding to this tensor are equal in both subsystems. These potentials are uniquely defined by substituting eq. (35) into eq. (30):

In one-component solids, |$\bar{\mathcal {F}} = \mathcal {F}^{\text{sol}}_{\text{mol}}$| and |$\bar{V} = V^{\text{sol}}_{\text{mol}}$| (the molar Helmholtz energy and molar volume). Uses of eq. (36) are demonstrated in the following two examples.

4.1 Example 1: Gibbs’ thought experiment

Gibbs (1878) (Chapter IX, his pages 343–360) envisaged a one-component solid held in a state of non-hydrostatic stress, with orthogonal faces aligned with the principal stresses of the crystal and in contact with separate liquid pools. He imagined that the pools were isolated from each other by membranes touching an edge of the solid, and that movable pistons were used to independently control the pressure in each pool (Fig. 2a). In his thought experiment, Gibbs imagined dissolution and precipitation taking place at the solid–liquid interfaces until the system was in chemical equilibrium. His expressions for the chemical potential of the solid component in the liquid at equilibrium (his eqs 393–395) can be written:

where the properties on the left-hand side belong to the solid. Eq. (37) is identical to eq. (36) after equating |$\mu ^{\text{liquid}}$| with |$\mu ^{\hat{\boldsymbol {n}}\otimes \hat{\boldsymbol {n}}}$| and |$P^{\text{liquid}}$| with |$-\sum _{i,j} \hat{n}_i \sigma _{ij} \hat{n}_j$|⁠. The uniaxial chemical potentials taken normal to an interface are constant across that interface at equilibrium (Figs 2b and c).

4.2 Example 2: melting of one-component stressed and unstressed bodies

Another thought experiment was proposed by Sekerka & Cahn (2004) to investigate the magnitude of melting point reduction by inducing non-hydrostatic stresses. Let us first define a linearized high-pressure equation of state for a one-component solid that can experience small network strains (⁠|$\boldsymbol {\varepsilon }^{\text{ntwk}}$|⁠) and changes in temperature (⁠|$\delta T$|⁠) away from a hydrostatic reference state (⁠|$P_0$|⁠, |$T_0$|⁠):

Now, consider a single crystal of this solid that is in equilibrium with a liquid of the same composition (that is, where the temperature is equal to the hydrostatic melting point) along one planar face with normal |$\hat{\boldsymbol {n}}$|⁠. At the hydrostatic reference state, chemical equilibrium satisfies the following equation:

The crystal is then stressed while keeping the pressure of the liquid constant. The temperature is allowed to change to maintain equilibrium. The chemical potential of the liquid is reduced by |$S^{\text{liq}}_{\text{mol}} \delta T$|⁠:

This liquid chemical potential must be equal to the uniaxial chemical potential in the solid (eq. 36). Mechanical equilibrium must be maintained between the solid and the liquid, such that:

Equating eqs 41 and (42):

Substituting eq. (38) into the right-hand side of eq. (44):

Elastic stiffness tensors are always positive semidefinite and the molar entropy of fusion |$S^{\text{fus}}_{\text{mol}}$| is also always positive, and therefore the melting temperature for a stressed solid must always be lower than that of its hydrostatic counterpart (Sekerka & Cahn 2004). To see how much lower, consider an elastically isotropic phase strained elastically to |$\sim$|1 per cent extension and compression normal to the liquid–solid interface. This strain is about an order of magnitude smaller than the theoretical elastic limit in silicate minerals at low temperature (e.g. Fujimura et al. 2023), but much larger than the strains at which viscous flow dominates at mantle temperatures on laboratory timescales (e.g. Durham & Goetze 1977). For isotropic minerals, bulk modulus |$K_T$| and shear modulus G define the isothermal stiffness tensor:

To provide an example calculation for |$\delta T$|⁠, we use the properties of periclase at the 1 bar melting point of 3105 K: a molar volume of 13 cm|$^3$| mol⁻¹ (Dubrovinsky & Saxena 1997), entropy of fusion of 25 J K⁻¹ mol⁻¹ (Chase et al. 1974) and shear modulus of 55 GPa (Hama & Suito 1999). Substituting into eq. (47):

that is, for these parameters, 1 per cent elastic strain results in a 6 K reduction in the melting point.

If the liquid contains more than one component, the solid will progressively dissolve in the fluid as the temperature rises. In this case, |$S^{\text{fus}}_{\text{mol}}$| can be replaced with the partial molar entropy of dissolution of the solid into the liquid, but eq. (47) remains otherwise the same. The consequence, as pointed out by Gibbs (1878) and numerous other researchers after him, is that a material held in a non-hydrostatic state of stress |$\boldsymbol {\sigma }$| is always less stable than a coexisting material of the same composition and phase held at a hydrostatic pressure P defined by the normal stress across a shared interface (⁠|$P = -\sum _{i,j} \hat{n}_i\sigma _{ij}\hat{n}_j$|⁠).

A final note: In the derivation of eq. (47), it was not necessary to define the state of strain in the solid. The only requirement was that the strain does not induce a change in stress resolved along |$\hat{\boldsymbol {n}}$|⁠, that is:

5 DISCUSSION

5.1 Comparison with previous works

Progress toward defining chemical potentials under nonhydrostatic stress began with Gibbs (1878), who used his classical thermodynamics to define equilibrium between stressed one-component solids and potentially multicomponent liquids. In his exposition, chemical potentials were only defined in the liquid. In Section 4.1, I demonstrate that the equilibrium relation derived by Gibbs is equivalent to equality in uniaxial chemical potentials.

Gibbs’ work left the challenge of extending his theory to multicomponent, non-hydrostatic systems unresolved. Herring (1951) generalized the concept of non-hydrostatic chemical potentials (his eq. 1):

but he did not unambiguously define |$\boldsymbol {\varepsilon }$| or the Gibbs energy |$\mathcal {G}$| under non-hydrostatic conditions. He later proposed a more specific equation (Herring 1953, as reported by Mullins & Sekerka 1985):

This equation involves the total domain volume, rather than a molar volume, and therefore could be argued to be equivalent to eq. (27) proposed here. However, it remains unclear whether the ‘shear’ variable held constant in the partial derivative refers to domain or material shear, and the equation itself does not indicate any role for deformation in determining the value of |$\mu _i$|⁠. Kamb (1961) and Paterson (1973), apparently unaware of Herring’s papers, reviewed the developments in non-hydrostatic thermodynamics since Gibbs (1878), and argued that it was not possible to meaningfully associate chemical potentials or Gibbs free energy with non-hydrostatically stressed solids. This led to the prevailing belief that there could be no comprehensive definition of chemical potentials for stressed solids.

A major post-Gibbsian theory was developed by Larché & Cahn (1973, 1978a, b, 1985), who circumvented the issue of defining individual chemical potentials by introducing relative chemical potentials in the context of lattice-based crystalline systems, alongside a ‘network constraint’ to maintain the number of lattice sites. They maintained the view of Kamb (1961) that it was unnecessary and perhaps impossible to define individual chemical potentials. This framework has since been used by many (e.g. Voorhees & Johnson 2004; Frolov & Mishin 2010a, b, 2012a, b; Powell et al. 2018; Hess et al. 2022). Larché & Cahn (1973) use the network strain rather than domain strain in their integral for the energy (their eq. 9, the equivalent of eq. 26 in this paper), and as a result their partial molar densities are not all independent (their eq. 16). This forces them to consider changes in the relative molar densities of chemical components (e.g. |$n^{\rho }_2 - n^{\rho }_1$|⁠), and leads them to define ‘diffusion potentials’ (Larché & Cahn 1978b, their eqs 18 and 19), rather than individual chemical potentials. In a sense, then, their theory is similar to the one developed here, but enforcing |$s_k = 0$| for all chemical components. It is interesting to note that Larché & Cahn (1985) cite Herring (1951), but do not comment on his definition of individual chemical potentials.

Unsatisfied by the lack of a unique chemical potential in non-hydrostatic thermodynamics, McLellan (1980) and Mullins & Sekerka (1985) attempted to generalize the concept of individual chemical potentials for stressed systems. McLellan (1980) introduced a Legendre transformation from which he could define a generalized Gibbs function. This transformation involved an ‘extensive strain tensor’ |$\Gamma _{ij} = V_0 \varepsilon _{ij}$|⁠, where each component had units of m|$^3$| rather than m. The resulting exact differential of the Gibbs function is (his eq. 20.1.17):

Unfortunately, |$\boldsymbol {\Gamma }$| is not a state property of the system unless it is referenced to some standard state, leading McLellan to admit that his Gibbs function was ‘not truly extensive’. Taking a different approach, Mullins & Sekerka (1985) proposed that individual chemical potentials could be defined for crystalline materials by starting from an exact differential for the internal energy per unit cell (their eqs 6 and 30):

where |$\boldsymbol {q} = \boldsymbol {M}^{\text{unit}T} \boldsymbol {M}^{\text{unit}}$|⁠. A puzzling aspect of this formulation is their claim that their vector |$\boldsymbol {n}^{\text{unit}}$| has n independent chemical variables for n components, despite being defined per unit cell. They note that “the extra independent variable is related to the concept of crystallinity and has no counterpart in a fluid”. No such special treatment was required in the derivations in Section 3.2.

The ‘network’ or ‘crystallinity’ constraints of Larché & Cahn (1985) and Mullins & Sekerka (1985), are mirrored in this study through the separation of domain strain into network and chemical components (Section 3.2). This separation relies on the ability to track the translation and deformation of network structure even during changes in chemical composition, but it does not strictly require a crystalline lattice and does not preclude a component of viscous flow. Even archetypal fragile liquids, such as water, can sustain a fixed ‘network’ and elastic stresses if timescales are sufficiently short (Schulz et al. 2020). Many liquids and glasses are inherently too fragile to maintain a fixed structure during chemical diffusion (Faupel et al. 2003), but some highly polymerized silicate glasses can maintain a fixed ‘network’ during diffusion by forming robust diffusion pathways (Greaves 1985; Greaves et al. 1991). Chemical diffusion in such glasses can even generate compressive stresses (Karlsson et al. 2010), which would not be possible in a fragile lattice.

Particularly interesting in light of the current work are two rarely cited papers by Bayly (1983, 1985). Bayly postulated a direction-dependent chemical potential (Bayly 1985, his eq. 1), defined as:

(symbols modified to mirror those used in this work). This expression is identical to eq. (36) if |$\mathcal {E}_k$|⁠, |$S_k$| and |$V_k$| are taken to be the partial derivatives of internal energy, entropy and volume with respect to mole number |$n_k$| at constant network strain and temperature. Bayly’s postulate was apparently without formal derivation, arising only as a logical extension of the thought experiment of Gibbs (1878). Section 3 rigorously formalises these variables, providing a consistent framework for defining chemical potentials that reduces to the correct expressions for both hydrostatic multicomponent systems and nonhydrostatic one-component systems.

Finally, Rusanov (2005) introduced a tensor-valued chemical potential for each chemical component, along with ‘change-of-mole-number tensors’ as their conjugate variables to account for the anisotropic addition and removal of material within a domain element (see Section 2.2 of his work). These tensor-valued chemical potentials are similar in spirit to the ‘chemical expansivity tensor’ presented in this study. However, the change of mole number tensors cannot function as independent thermodynamic variables. Instead, the net effect of anisotropic material addition manifests itself in the overall network strain rate, while variations in material addition from component to component manifest solely as changes in concentration—any evidence of how each component was added to the domain is lost by local chemical diffusion.

5.2 Disequilibria in heterogeneously stressed systems

In recent years, a substantial scholarly debate has arisen between Wheeler (2014, 2018, 2020) and several critics, including Fletcher (2015), Hobbs & Ord (2016) and Powell et al. (2018). Central to this debate is the extent to which variations in chemical potential in stressed solid systems influence geological processes. Wheeler (2014, 2018, 2020) argues that deviatoric stresses can have a large influence on mineral stability. His critics argue that deviatoric stresses have a minor influence on mineral stability. Both sides draw extensively on the work of Gibbs (1878) and Larché & Cahn (1985) to support their case.

The differences in views are more minor than they appear. Wheeler, in his initial review of this paper, observes that “if we apply stress to a rock, it is true that stresses parallel to a particular interface may change, while normal stresses remain fixed. The effect on chemical potential is then small.” This perspective aligns with that of Fletcher (2015), Hobbs & Ord (2016) and Powell et al. (2018), who, citing the works of Sekerka & Cahn (2004) and Frolov & Mishin (2010a), argue that the influence of deviatoric stress on mineral stability at any given interface is minimal.

The core argument of Wheeler (2014, 2018, 2020) is that spatially separated variations in chemical potentials in liquids within heterogeneously stressed media can be large, and that these differences drive reactions. This is undeniably true, and, as Wheeler notes, there is a long history of supporting observations (e.g. Rutter 1983). Much more controversial is the idea that one can usefully define a chemical affinity for a reaction between spatially separated grains (Wheeler 2014). As highlighted in different ways by Fletcher (2015), Hobbs & Ord (2016) and Powell et al. (2018), this idea violates the principle of locality. As an example, consider a case in which an interconnected fluid interacts with many grains, located at considerable distances from one another. Although the fluid serves as a continuous medium that physically connects these grains, the communication of chemical signals—namely, changes in chemical potential—is limited by chemical diffusion, which is an inherently local phenomenon. As Powell et al. (2018) explains, “the scenario considered by Wheeler should be envisaged as a set of smaller, spatially separate systems involving different equilibria, with chemical potential and pressure gradients between the systems”.

A related point of contention is the claim that chemical equilibrium cannot exist in a stressed solid (Wheeler 2018, 2020). In idealized systems with specific boundary conditions, chemical equilibrium can certainly be achieved in a heterogeneously stressed medium, as exemplified by the thought experiment of Gibbs (1878) (Section 4.1). However, such restrictive boundary conditions are rare in natural systems. Viscous deformation and chemical diffusion in natural systems act to dissipate stress and chemical potential gradients. The question then is whether chemical diffusion can be fast relative to structural relaxation by viscous flow. The answer to this question is dependent on the system and phases involved; in a fragile, inviscid liquid, viscous relaxation will be much faster than chemical diffusion. In contrast, in a robust crystalline structure like garnet, viscous relaxation is very slow relative to diffusive exchange of similar components such as Mg and Fe. In such systems, quasi-chemical equilibrium may be reached despite the maintenance of significant deviatoric stresses. Whatever the system, chemical potentials and their spatial gradients offer a natural framework to model and track the progression towards equilibrium.

5.3 A generalized law of chemical diffusion?

Under hydrostatic conditions, spatial gradients in chemical potentials result in a thermodynamic force:

Assuming a linear relationship between |$\boldsymbol {F}_k$| and the resulting chemical fluxes |$\boldsymbol {J}_j$| results in a phenomenological law of chemical diffusion:

where |$J_{ij}$| [mol m⁻² s⁻¹] is the flux of component i in direction j and |$\boldsymbol {D}$| is a diffusivity [mol|$^2$| J⁻¹ m⁻¹ s⁻¹]. For full details of these equations and the other constraints that must be considered in building a model of chemical diffusion, see (De Groot & Mazur 1984, their Section XI.2).

Although beyond the scope of this study, it is tempting to speculate on what a generalized law of chemical diffusion under non-hydrostatic conditions might look like. The uniaxial chemical potentials (eq. 36) can be written

and this equation bears a striking similarity to the equation for the normal stress (⁠|$\sigma ^{\hat{\boldsymbol {n}}\otimes \hat{\boldsymbol {n}}} = \sum _{l,m} \hat{n}_l \left(\sigma _{lm} \right) \hat{n}_m$|⁠). By analogy, then, we might consider uniaxial chemical potentials as the normal components of ‘chemical tractions’ acting on a surface. In the examples in Section 4, only this normal component was considered, but in a continuous medium, we might also consider a component of chemical potential acting parallel to the surface. The resultant ‘chemical traction’ would then have the form:

Integrating this traction over the surface of a domain element and using the divergence theorem, we might expect the resultant thermodynamic force to have the form:

This expression satisfies the two examples in the previous section, because the principal stresses and chemical forces in both examples are constrained to be parallel and perpendicular to |$\hat{\boldsymbol {n}}$|⁠. The forces |$\boldsymbol {F}_k$| are invariant under rotation, and reduce to eq. (56) under hydrostatic conditions. Further work is needed to confirm that this formulation is correct, and to build a complete model to simulate chemical diffusion under nonhydrostatic conditions.

6 CONCLUSIONS

In this paper, I have introduced a general definition of chemical potentials in multicomponent, non-hydrostatically stressed phases (Section 3.2):

The tensor |$\boldsymbol {\varepsilon }$| that is held fixed in eq. (61) corresponds to the strain of the domain. The rate of change of this tensor has components due to network strain |$\dot{\boldsymbol {\varepsilon }}^{\text{ntwk}}$| and a chemical strain due to changes in molar amounts of the chemical components |$\dot{\boldsymbol {n}}$| within the domain element:

The chemical expansivity tensor |$\boldsymbol {\Xi }$| describes the change in the shape of the domain as material is added or removed at constant network strain:

This tensor must satisfy the equality |$\operatorname{Tr} (\boldsymbol {\Xi }) = 1$|⁠, but it is unconstrained by the state of the material; it is instead determined by how material is exchanged between the local domain and its surroundings. If material can only be added to or removed from a domain through an interface with normal |$\hat{\boldsymbol {n}}$|⁠, the chemical expansivity tensor |$\boldsymbol {\Xi }$| is equal to |$\hat{\boldsymbol {n}}\otimes \hat{\boldsymbol {n}}$| (Section 4).

ACKNOWLEDGEMENTS

I would like to thank John Rudge, John Wheeler and an anonymous reviewer for their careful reviews of this work. Their comments greatly improved the paper. I also thank Bruce Hobbs for his encouragement while working on this and previous papers. Any errors or omissions are my own. This work was supported by the Natural Environment Research Council (NERC Large Grant MC-squared; Award No. NE/T012633/1) and the Science and Technology Facilities Council (STFC Grant No. ST/R001332/1).

DATA AVAILABILITY

This paper is not associated with any data.

REFERENCES

APPENDIX A: HELMHOLTZ MINIMIZATION AND ENTROPY MAXIMIZATION WITH CONSTRAINTS

Spontaneous events within closed, fixed volume, constant temperature systems act to reduce the Helmholtz energy. This follows directly from the second law of thermodynamics, that is, the natural tendency of the universe to evolve toward a more disordered (and therefore more probable) state over time. The connection between entropy maximization in an isolated system and Helmholtz energy minimization in a closed, thermally equilibrated system is explained in detail by Callen (1985). In chapter 1, he introduces fundamental thermodynamic concepts, followed by the Gibbs–Duhem equality in chapter 3 (section 2), classifications of processes (possible versus impossible, quasi-static versus reversible) in chapter 4 (sections 1 and 2), Legendre transformations, internal parameters and energy minimization in chapter 5, and the Helmholtz potential in chapter 6 (section 2).

The current paper is particularly concerned with constrained minimization of the Helmholtz energy. The constraint in this case is the requirement that all deformation in the system is either elastic or accommodated by growth or dissolution of existing phases—no viscous relaxation is allowed. This constraint imposes restrictions on internal parameters of the system; that is, unconstrained variables that adjust to reduce the Helmholtz energy while keeping system shape, temperature and composition unchanged. In this paper, the internal parameters of interest represent the cell tensors of individual phases, which encode information about lattice orientation, molar volume and non-hydrostatic strain (Myhill 2022, 2024).

The result of constrained minimization is that the Helmholtz energy reached will in general be higher than the global minimum. As a result, the systems considered here are metastable relative to those in which viscous relaxation can occur. Completely neglecting intraphase viscous or plastic relaxation in the solids considered here represents an idealized assumption. However, the fundamental physics remain unchanged whatever the rheology of the coexisting phases: the process of mechanical and chemical equilibration emerges from the balance between energy minimization and kinetic constraints, with the system evolving toward the lowest Helmholtz energy state permitted by the constraints.

APPENDIX B: THE GIBBS–DUHEM EQUATION

The hydrostatic Gibbs–Duhem Equation (Duhem 1897; Gibbs 1906) describes the relationship between changes in chemical potential and changes in pressure and temperature (see the historical perspective in Dais 2021):

To generalize this expression to non-hydrostatic states, we start from the total differential for the Helmholtz energy density given by eq. (18):

Splitting the change in mole numbers into mole proportions and total mole number (⁠|$n^{\rho }_k = n^{\rho } p_k$|⁠) yields the following expression:

Holding network strain, temperature and mole proportions constant and then integrating from 0 to |$n^{\rho }_k$| moles of each material per unit volume, we obtain an expression for the total Helmholtz energy density:

and so, by the chain rule

Equating this expression with eq. (B2):

and by using the same substitution as that used in eq. (5):

integrating over a homogeneous volume

The total change in volume is equal to the sum of elastic and chemical contributions:

Substituting:

Using the isotropic strain chemical potential defined in eq. (32):

After some rearranging, we obtain the non-hydrostatic Gibbs–Duhem equation:

Eq. (B15) reduces to eq. (B1) under hydrostatic conditions as required.