A model of elastic softening and displacive phase transitions in anisotropic phases, with application to stishovite and post-stishovite

SUMMARY

This paper introduces a comprehensive framework for modelling both instantaneous and time-dependent elastic softening in anisotropic materials at high pressure and temperature. This framework employs Landau Theory, minimizing the Helmholtz energy by varying isochemical parameters (⁠|$\boldsymbol {q}$|⁠) that capture structural changes, atomic ordering and/or electronic spin states. This allows for internally consistent predictions of volume, unit cell parameters, the elastic tensor, and other thermodynamic properties, while allowing large symmetry-breaking strains. The formulation is validated using the stishovite-to-post-stishovite transition. It is demonstrated that, near this transition, both stishovite and post-stishovite exhibit auxetic behaviour in several directions, with post-stishovite also displaying negative linear compressibility along the long axis of its unit cell (either the a- or b-axis). The new formulation is implemented in the open-source BurnMan software package.

1 INTRODUCTION

1.1 Isochemical variability in crystalline materials

The thermodynamic and elastic properties of crystalline phases are functions of the number of moles of independent chemical species |$\boldsymbol {n}$| that comprise them. In many solutions, phases can be adequately modelled by considering only species that are compositionally independent of each other, such that there are the same numbers of independent chemical components|$\boldsymbol {x}$| and independent chemical species|$\boldsymbol {n}$|⁠. For example, (Mg,Fe)|$_3$|Al|$_2$|Si|$_{3}$|O|$_{12}$| garnet (a two-component system, |$n_x=2$|⁠) can be reasonably modelled as a mixture of two species (⁠|$n_n=2$|⁠): pyrope [Mg]|$_{3}$|Al|$_{2}$|Si|$_{3}$|O|$_{12}$| and almandine [Fe]|$_{3}$|Al|$_{2}$|Si|$_{3}$|O|$_{12}$|⁠. However, there exist many nominally pure phases and solid solutions whose properties can only be accurately modelled by explicitly considering one or more isochemical degrees of freedom |$\boldsymbol {q}$|⁠. Three important types of isochemical degrees of freedom with examples involving a single degree of freedom are:

1. Ordering: For example, Mg–Si exchange on two distinct octahedral sites in majorite (Heinemann et al. 1997):

   $$\begin{eqnarray} q = (p_{{\text{Mg}^{\text{M1}}}} + p_{{\text{Si}^{\text{M2}}}}) - (p_{{\text{Si}^{\text{M1}}}} + p_{{\text{Mg}^{\text{M2}}}}) \end{eqnarray}$$

   where, for example, |$p_{{\text{Mg}^{\text{M1}}}}$| is the proportion of Mg on the M1 site.

2. Electronic spin state: For example, variable proportions of high spin and low spin |${\text{Fe}^{2+}}$| in ferropericlase (Wu et al. 2013) or troilite (Urakawa et al. 2004):

   $$\begin{eqnarray} q = p_{{\text{Fe}^{\text{M}}_{\text{HS}}}} - p_{{\text{Fe}^{\text{M}}_{\text{LS}}}} \end{eqnarray}$$
3. Structural flexibility: For example, tilt angle |$\Theta$| of tetrahedral SiO|$_{4}$|-units in quartz (Carpenter et al. 1998; Wells et al. 2002), or octahedral SiO|$_{6}$|-units in CaCl|$_{2}$|-structured post-stishovite (Zhang et al. 2023):

   $$\begin{eqnarray} q = \Theta \end{eqnarray}$$

Materials with isochemical degrees of freedom |$\boldsymbol {q}$| are often associated with anomalous thermodynamic and elastic properties.

1.2 Anomalous properties associated with isochemical variability

In “normal” materials, most thermodynamic properties are expected to increase or decrease monotonically with temperature and pressure. For example, thermal expansivity, heat capacity and bulk compressibility are all expected to increase with temperature and decrease with pressure. However, not all materials behave in this way. In one of the earliest observations of anomalous volumetric change, the unit cell volume of quartz increased rapidly with increasing temperature up to 570 |$^{\circ }$|C, and then actually decreased above 570 |$^{\circ }$|C (Le Chatelier 1890). This peculiar behaviour is accompanied by reduced elastic moduli (Carpenter et al. 1998; Lakshtanov et al. 2007) and a peak in heat capacity (Grønvold & Stølen 1992). The change in behaviour at 570 |$^{\circ }$|C is now known to mark a displacive, weakly first-order trigonal to hexagonal transition (e.g. Bachheimer & Dolino 1975; Dolino 1990; Antao 2016). At low temperatures, the trigonal structure contains SiO|$_{4}$| tetrahedra that are tilted relative to each other. As the structure is heated, the relative tilt angle decreases, and disappears at the transition, causing the observed anomalies below 570 |$^{\circ }$|C. It turns out that the anomalous properties above 570 |$^{\circ }$|C are also related to tilting of SiO|$_{4}$| tetrahedra, but in this case the tilting is dynamic (thermally driven), in contrast to the static tilting seen in the trigonal phase (Welche et al. 1998; Wells et al. 2002).

Since the discovery of the anomalous properties of quartz, many other phases have been found to exhibit anomalous variations in heat capacities, thermal expansivities and elastic moduli (Helgeson 1978; Thompson & Perkins 1981; Wadhawan 1982; Salje 1985; Dove 1997; Carpenter 2006). These include some of the most abundant phases in the Earth; feldspar (Mookherjee et al. 2016; Lacivita et al. 2020), garnet (Heinemann et al. 1997), orthopyroxene (Jackson et al. 2004), amphibole (Cámara et al. 2003), cristobalite (Yeganeh-Haeri et al. 1992) stishovite (Kingma et al. 1995), ferropericlase (Wu et al. 2013) and CaSiO|$_{3}$| perovskite (Stixrude et al. 2007) all exhibit anomalous properties. Many of the anomalies are associated with displacive phase transitions. In all cases, the anomalous properties arise from the presence of one or more isochemical degrees of freedom in the crystalline structure.

1.3 Aim of this paper

Anomalous thermodynamic behaviours of crystalline materials have long been effectively described by models based on Landau Theory (e.g. Landau 1935). These models, which rely on the concept of energy minimization by varying internal parameters |$\boldsymbol {q}$|⁠, are successful due to their robust thermodynamic and group-theoretical foundations. However, despite their success, current applications to study elastic behaviour have limitations. Some lack self-consistency, and others rely on simplifying assumptions, such as small symmetry-breaking strains (Tröster et al. 2002, 2014, 2017). Furthermore, most models are confined to either isobaric (Lüthi & Rehwald 1981), or isothermal conditions (Carpenter et al. 2000; Tröster et al. 2014, 2017). Self-consistent high pressure, high temperature models are restricted to modelling scalar properties such as the bulk modulus, rather than the full elastic tensor (e.g. Angel et al. 2017).

The aim of the current paper is to develop a complete anisotropic theory that extends to high-pressure, high-temperature conditions and is not limited to small symmetry-breaking strains. Models based on this theory can take pressure, temperature and composition as input, and output any desired thermodynamic and elastic properties, such as the thermal expansivity or isentropic stiffness tensors. The theory presented in this paper uses a new anisotropic equation of state for solid solutions (Myhill 2024) that is implemented in the BurnMan open source software project (Cottaar et al. 2014; Myhill et al. 2023, 2024). Landau theory is incorporated into the equation of state via non-ideal mixing between isochemical endmember twins, a method that is directly analogous to the order–disorder models found in many thermodynamic data sets (Holland & Powell 2011; Stixrude & Lithgow-Bertelloni 2024). Symbols used in this paper are provided in Table 1.

2 HISTORICAL MODEL DEVELOPMENT

2.1 Types of variables

Throughout this paper, the following terminology is applied to different types of variables:

1. Natural variables: The independent variables or parameters that define thermodynamic potentials such as the hydrostatic Gibbs energy |$\mathcal {G}(P, T, \boldsymbol {x})$| or the single-crystal anisotropic Helmholtz energy |$\mathcal {F}(\boldsymbol {F}, T, \boldsymbol {x})$| (Nye et al. 1985; Holzapfel 2000; Myhill 2024). Spontaneous processes within a system will tend to decrease the thermodynamic potential toward a minimum when the natural variables of that potential are held constant.

2. Internal variables: Variables that describe the internal state of the system but are not controlled or directly manipulated from outside the system. Spontaneous changes in these variables allow the thermodynamic potential to approach a minimum. In this paper, internal variables are indicated by the symbol |$\boldsymbol {q}$|⁠.

3. Conjugate variables: The derivatives of the thermodynamic potentials with respect to their natural variables. These are related to properties of the system. For example, Cauchy stress is related to the first derivative of the Helmholtz energy density with respect to strain:

   $$\begin{eqnarray} \left(\frac{\partial \mathcal {F}^{\rho }}{\partial \varepsilon _{ij}}\right)_{T, \boldsymbol {x}} &= \sigma _{ij}, \end{eqnarray}$$

   (1)

   $$\begin{eqnarray} d \varepsilon _{ij} &= \frac{1}{2} (d F_{ik} F_{jk}^{-1} + d F_{jk} F_{ik}^{-1} ). \end{eqnarray}$$

   (2)

2.2 Landau theory

In the 1930s, Landau presented a quantitative theory predicting the shape and magnitude of anomalous peaks in heat capacity and other phenomena involving continuous transitions (Landau 1935, 1937a, b). English translations of these works can be found in Landau (2008) and Ter Haar (2013). Landau (1935) theorized that such anomalies were a natural consequence of the system minimizing energy by varying internal variables |$\boldsymbol {q}$|⁠. To demonstrate his theory, he constructed a thermodynamic potential |$\Phi$| as a polynomial in q relative to a reference state:

where |$\Phi _0$|⁠, a, b, c, |$\ldots$| are all potentially functions of the natural variables of |$\Phi$|⁠. Landau implicitly used the Gibbs energy (⁠|$\mathcal {G}(P,T,\boldsymbol {x})$|⁠; e.g. Landau 1935) and Helmholtz energy (⁠|$\mathcal {F}(V,T,\boldsymbol {x})$|⁠; e.g. Landau & Ginzburg 1950) as convenient thermodynamic potentials for different problems. Assuming that q varied rapidly to minimize the potential with changes in temperature, the heat capacity could be expressed as:

where |$q^{*}$| is the value of q that minimizes |$\mathcal {G}$| at the pressure and temperature of interest. Landau later showed that symmetry-breaking displacive transitions could be modelled using an energy expansions that were symmetric about |$q=0$|⁠, with |$\Phi _0$| representing the energy of the high-symmetry phase (Landau 1937a):

If |$\partial ^2 \Phi / \partial q \partial q> 0$| at |$q=0$|⁠, then the high-symmetry state is more stable than states where the symmetry is slightly broken. The high-symmetry state either represents the globally stable state, or at least a metastable state. Conversely, if |$\partial ^2 \Phi / \partial q \partial q< 0$|⁠, then the high-symmetry state is unstable, and two local minima in the energy will appear at positive and negative q. These minima correspond to the q values of stable ‘twins’ that are identical in every respect except that they are mirror images about the symmetry element that has been lost. If a change in conditions causes |$q=0$| to lose or gain its status as the global minimum, a symmetry-breaking transition can take place.

Landau (1937a) motivated his use of a single internal variable q to describe symmetry breaking (rather than a large vector of variables) as follows. Let the spatial probability density function of different atoms be equal to |$\boldsymbol {\rho } = \boldsymbol {\rho }(x, y, z)$|⁠. Furthermore, let the difference between the equilibrium and high-symmetry structures be equal to |$\delta \boldsymbol {\rho }$|⁠. The perturbation |$\delta \boldsymbol {\rho }$| can be split into a number of normalized irreducible functions (i.e. a basis set of functions) |$\boldsymbol {I}_i$| that satisfy the symmetry of the parent group—stransition (Stokes & Hatch 1988; Powell 2010):

The excess energy |$\Phi _{\text{xs}}$| can be given as a function of |$\delta \boldsymbol {\rho }$| or, equivalently, by the values of |$c_i$|⁠:

The vector of internal variables |$\boldsymbol {c}$| potentially has a large number of components. To simplify the system, we assume that all the accessible values of |$\boldsymbol {c}$| in any given state lie along some line in |$\boldsymbol {c}$|-space, and that the distance along this line from the origin can be parametrized by a single variable q. It is then possible to write both |$\Phi$| and |$\boldsymbol {c}$| as a function of this parameter:

The excess energy |$\Phi _{\text{xs}}$| can then be expanded in powers of q, as in eq. (5). In Landau (1937a), |$q^2 = \sum c_i^2$|⁠, but such a restrictive choice is not required.

2.3 Coupled models of elastic strain around phase transitions

The models considered by Landau focused on the evolution of systems as a function of temperature. Subsequent models, informed by work on piezoelectrics and ferroelectrics such as Rochelle Salt (Mueller 1940; Ginzburg 1945; Devonshire 1949), considered systems with more than one natural variable and also systems where conjugate variables were constrained, rather than the natural variables. The resulting ‘Landau–Ginzburg–Devonshire Theory’ (Levanyuk & Sannikov 1969, 1970, 1971) has been applied to studies of elastic effects during displacive transitions since the 1970s (Dvořák 1971; Höchli 1972; Lüthi & Rehwald 1981; Kityk et al. 2000). Under unstressed conditions (⁠|$\boldsymbol {\sigma } = \boldsymbol {0}$|⁠), and with one internal variable (q), the Helmholtz energy density in the small strain approximation can be written as a set of simultaneous equations:

where b, d, f, |$\lambda _{i,m,n}$|⁠, |$\mathbb {C}_{ik}^{0}$|⁠, |$\ldots$| are all potentially functions of temperature, |$\hat{\boldsymbol {\varepsilon }}$| is strain relative to the high-symmetry form, and ‘strain-order coupling’ coefficients |$\lambda _{i,m,n}$| are only non-zero when allowed by the symmetry of the system (Carpenter & Salje 1998). In this expression, T and |$\hat{\boldsymbol {\varepsilon }}$| are natural variables of the system. By substituting eq. (10) back into eq. (9), a classic Landau-type expression is recovered:

where the modified coefficients |$b^{\prime }$|⁠, |$d^{\prime }$| and |$f^{\prime }$| are functions of b, d, f, |$\lambda _{i,m,n}$| and |$\mathbb {C}_{ik}^{0}$|⁠.

Attempts have been made to apply this theory to phase transitions at high pressure. Some studies (Liakos & Saunders 1982; Carpenter et al. 2000; Buchen 2021; Zhang et al. 2021), use eqs (9) and (10) essentially unaltered, except that the Helmholtz energies are relabelled as Gibbs energies, and |$\mathcal {F}_0$| (now |$\mathcal {G}_0$|⁠) becomes a function of pressure rather than temperature. The derivative of the energy with respect to strain (eq. 10) is still set to zero, despite the pressure not being equal to zero. This method has been successful in reproducing observed changes in elastic tensor, but the formulation is thermodynamically inconsistent and produces conflicting values of physical properties such as the volume (see Supporting Information). To achieve thermodynamic consistency in the small strain limit, Tröster et al. (2014, 2017) introduced an additional |$\sigma _{ij} \varepsilon _{ij}$| term to eq. (9), and set |$\partial \mathcal {F}^{\rho } / \partial \hat{\boldsymbol {\varepsilon }} = \boldsymbol {\sigma }$|⁠. The use of this formulation still requires self-consistent expressions for the Helmholtz energy, cell tensor |$\boldsymbol {M}$| and isothermal stiffness tensor |$\mathbb {C}_T$| of the high-symmetry phase as a function of pressure and, ideally, temperature, which is prohibitively complex using classic expansions in finite strain (Tröster et al. 2017).

In (Myhill 2022, 2024), I presented an alternative approach to building self-consistent models of anisotropic material properties in solid solutions through use of a fourth order anisotropic tensor |$\mathbb {\Psi }(V^{\text{mol}}, T, \boldsymbol {p})$|⁠. In Section 4, I couple this approach to Landau theory by decomposing the endmember proportions |$\boldsymbol {p}$| into natural variables |$\boldsymbol {x}$| and internal variables |$\boldsymbol {q}$| (Myhill & Connolly 2021):

where |$\boldsymbol {A}$| and |$\boldsymbol {B}$| are constant transformation matrices. However, first, we must expand on Landau’s models of heat capacity (e.g. eq. 4) to determine how anisotropic properties behave under energy-minimizing internal variables q (Section 3).

3 MATERIAL PROPERTIES AND RELAXATION

3.1 Relationships between state, material properties and thermodynamic potentials

The beauty of thermodynamics is that observable properties of materials arise from the partial derivatives of thermodynamic potentials with respect to their natural variables. Taking the Helmholtz energy density |$\mathcal {F}^{\rho }$| relative to an initially hydrostatic state as an example, first derivatives with respect to the natural variables define the Cauchy stress |$\boldsymbol {\sigma }$|⁠, entropy density |$S^{\rho }$| and chemical potentials |$\mu _i$|⁠:

Integrating over a small homogeneous volume V yields the exact differential of the extensive Helmholtz energy:

where |$\boldsymbol {\varepsilon }$| is the small domain-strain tensor, which can be determined from changes in the ‘extensive cell tensor’ |$\boldsymbol {M}$|⁠, a second order tensor describing the shape of the unit cell but with the volume of the material. The second derivatives of the Helmholtz energy density with respect to strain and temperature when the compositional natural variables |$\boldsymbol {x}$| are fixed yield the isothermal stiffness tensor (⁠|$\mathbb {C}_{\text{T}}$|⁠), the thermal stress tensor (⁠|$\boldsymbol {\pi }$|⁠, that describes how internal stress changes with temperature if deformation is not allowed to occur) and the isometric heat capacity density (⁠|$C^{\rho }_{\boldsymbol {\varepsilon }}$|⁠):

In these expressions, we have explicitly fixed only the compositional variables |$\boldsymbol {x}$|⁠, rather than the full set of species variables |$\boldsymbol {n}$|⁠. We have not specified whether the internal variables |$\boldsymbol {q}$| are fixed or allowed to vary, and therefore these equations are not yet uniquely defined. Once the behaviour of |$\boldsymbol {q}$| has been constrained (Section 3.2), other thermodynamic properties may be determined through appropriate algebraic manipulation (Davies 1974; Nye et al. 1985; Holzapfel 2000):

Under hydrostatic conditions, the following scalar properties can be defined using the properties defined above:

See Appendix  A for derivations of eqs (20), (23) and (24).

3.2 The concept of relaxation and unrelaxed properties

As the natural variables of a system change, a simultaneous change in the internal variables |$\boldsymbol {q}$| might be able to reduce the appropriate thermodynamic potential (i.e. increase the entropy of the universe). Some changes, such as tilting of structural units or changes in electronic spin state, are extremely rapid with respect to typical rates of change of strain and temperature (Slonczewski & Thomas 1970; Kimizuka et al. 2003; Zhang et al. 2018). Others, such as order–disorder processes, are often sluggish because the activation energy barrier required for the migration of chemical species is larger than the average thermal energy of these species (Seifert & Virgo 1975; Ganguly 1982; Redfern 1998, 2000; Redfern et al. 1999). If |$\boldsymbol {q}$| are not allowed to change during changes in state, then we can define a set of ‘unrelaxed’ properties of the system:

If |$\boldsymbol {q}$| responds rapidly to minimize energy compared with the timescales over which the natural variables change, then a set of modified second derivatives can be defined, and these are referred to as ‘relaxed’ properties (Section 3.3). Between these two endmember cases, there is an intermediate situation where changes in |$\boldsymbol {q}$| take place on similar timescales to changes in state, in which case the material can be described as ‘partially relaxed’ (Section 3.4).

3.3 The properties of a rapidly relaxed elastic material

In the scenario that one or more isochemical order parameters |$\boldsymbol {q}$| change instantaneously to minimize |$\mathcal {F}$| under small changes in strain and temperature, variational calculus can be used to calculate relaxed properties. Relaxed values of |$\boldsymbol {q}$| (denoted |$\boldsymbol {q}^{*}$|⁠) are a function of the other variables of the system:

Collecting small strain and temperature together as |$\boldsymbol {z} = \lbrace \boldsymbol {\varepsilon }, T\rbrace$|⁠, the change in |$\boldsymbol {q}^{*}$| with change in |$\boldsymbol {z}$| at fixed |$\boldsymbol {x}$| can be determined using the following expression (Appendix  B):

where |$\boldsymbol {R}$| is the left inverse matrix of |$\partial ^2\mathcal {F}^{\rho } / \partial q_l \partial q_m$|⁠:

The relationship between the relaxed and unrelaxed second derivatives of the Helmholtz energy density with respect to state is then given by the following expressions:

The values of |$\boldsymbol {q}^{*}$| are state properties; they have unique values as a function of |$\boldsymbol {M}^{\text{mol}}$|⁠, T, |$\boldsymbol {x}$| and do not depend on the path taken to achieve that state. Therefore, the properties calculated in eq. (36), which correspond to relaxed properties at fixed strain and temperature, can be used in combination with the equations in Section 3.1 to determine relaxed properties at fixed entropy or stress. For example, the equation for the relaxed isentropic stiffness tensor can be written:

which is directly analogous to the unrelaxed expression (eq. 20).

3.4 Time-dependent seismic relaxation

If the timescales of relaxation are comparable with the timescales of perturbations in strain, neither unrelaxed properties (Section 3.1) or relaxed properties (Section 3.3) will be able to correctly predict effective material properties. For example, passage of a seismic wave through a material may be too rapid for |$\boldsymbol {q}$| to keep pace with the changes in strain at constant entropy. To model time-dependent seismic relaxation, we can assume that relaxation of each structural parameter |$q_i$| is driven by the gradient in the internal energy density |$\mathcal {E}$| with respect to that parameter (analogous to Redfern et al. 1999, replacing the Gibbs energy with the internal energy density):

In the limit of small strains (and therefore small changes in |$\boldsymbol {q}$|⁠) and at constant entropy density and internal variables |$\boldsymbol {x}$|⁠, the change in internal energy density with respect to |$\boldsymbol {q}$| can be approximated as linear in |$\boldsymbol {\varepsilon }$| and |$\boldsymbol {q}$|⁠:

This expression can be written as a non-homogeneous linear system of equations with the form

where

If there is only one parameter q, the behaviour of q is analogous to the deformation of a forced Kelvin–Voigt body (a viscous dashpot in parallel with a spring). If there are multiple order parameters |$\boldsymbol {q}$|⁠, the solution to the non-homogeneous system of equations is (Coddington & Levinson 1984):

The steady state response of |$\boldsymbol {q}$| to a periodic sinusoidal variation in |$\boldsymbol {\varepsilon }$| is a sinusoidal variation with a phase lag. The change in |$\boldsymbol {q}$| drives changes in stress and temperature:

where |$\mathbb {C}_{\text{S}}$|⁠, |$C^{\rho }_{\boldsymbol {\varepsilon }}$| and |$\boldsymbol {\pi }$| are evaluated at fixed |$\boldsymbol {q}$|⁠. Analytical expressions for the second order partial derivatives of the internal energy as a function of the partial derivatives of the Helmholtz energy required by eqs (42), (43), (45) and (46) are as follows (Appendix  E):

3.5 Choosing natural variables for geological materials

So far, we have used |$\lbrace \boldsymbol {\varepsilon }, T, \boldsymbol {q}, \boldsymbol {x}\rbrace$| as an independent set of variables. However, |$\boldsymbol {\varepsilon }$| is a small strain tensor, and materials experience large strains over geologically relevant pressure and temperature ranges. One notable property of geological materials is that they tend to deform elastically only up to 1 or 2 per cent non-hydrostatic strain, and so one may consider equations of state that use |$\lbrace V, \bar{\boldsymbol {\varepsilon }}, T, \boldsymbol {q}, \boldsymbol {x}\rbrace$| as an independent variable set, where:

The variable |$\bar{\boldsymbol {\varepsilon }}$| is an infinitesimal isochoric strain away from a hydrostatic state; it has only five independent variables (there is zero volume strain) and the two variable sets have the same total number of independent variables. The partial derivatives required by eqs (33), (34) and (36) can be calculated by a change of variables. At fixed |$\boldsymbol {x}$|⁠, the following expressions can be derived (Appendix  C):

where |$\mathbb {C}_{\text{T}}$|⁠, |$\boldsymbol {\beta }_{\text{T}}$|⁠, |$\beta _{\text{TR}}$|⁠, |$\boldsymbol {\alpha }$| and |$\alpha _V$| are all unrelaxed properties.

Another change of variables is required if pressure is a more convenient natural variable than volume, as is common for geological applications. At fixed |$\boldsymbol {x}$| and |$\bar{\boldsymbol {\varepsilon }}=\boldsymbol {0}$| (i.e. hydrostatic conditions; Appendix  D):

where |$K_{\text{TR}}$| and |$\alpha _V$| are unrelaxed properties.

4 BUILDING A SELF-CONSISTENT LANDAU MODEL USING AN ANISOTROPIC EQUATION OF STATE

4.1 An anisotropic equation of state

Any complete model of elastic softening must be tied to a specific equation of state. In a recent paper (Myhill 2024), I proposed a self-consistent anisotropic equation of state for solid solutions that provided expressions for the Helmholtz energy, extensive cell tensor |$\boldsymbol {M}$| and isothermal elastic stiffness tensor |$\mathbb {C}_T$| as a function of volume V, temperature T and endmember proportions |$\boldsymbol {n}$| under strictly hydrostatic conditions. If we now consider small isochoric strains |$\bar{\boldsymbol {\varepsilon }}$| away from the hydrostatic state at constant temperature (eqs 50 and 51), we obtain the following expressions for the Helmholtz energy density and unit cell tensor:

Eqs (58) and (59) are consistent with a small strain expansion of the Helmholtz energy density about a reference volume at fixed temperature and composition:

The hydrostatic part of the Helmholtz energy density |$\mathcal {F}^{\rho \text{hyd}}(V^{\text{mol}}, T, \boldsymbol {p})$| can be calculated from any scalar equation of state. If pressure is a more convenient independent variable than volume, it can be calculated from the Gibbs energy:

In Myhill (2024), self-consistent expressions for |$\boldsymbol {M}$| and |$\mathbb {C}_{\text{T}}$| were provided by first defining a fourth order tensor |$\mathbb {\Psi }$|⁠, expanded as a Taylor series in molar proportions:

with parameters that satisfy the condition:

at all compositions and at all conditions. In orthotropic materials (see Section 2 of Myhill 2024, for a formulation for non-orthotropic materials):

The unrelaxed thermal expansivity |$\boldsymbol {\alpha }$|⁠, isothermal compressibility |$\boldsymbol {\beta }_{\text{T}}$| and isothermal compliance tensor |$\mathbb {S}_{\text{T}}$| at fixed composition can be found using:

where again, these expressions are only correct for orthotropic materials. These partial derivatives provide all the unrelaxed properties introduced in Section 3.2 and used in Section 3.3. Other anisotropic properties may be obtained from the expressions in Section 3.1.

4.2 Compositional derivatives of the Helmholtz energy

In addition to unrelaxed properties taken at constant |$\boldsymbol {q}$|⁠, the expressions for relaxed properties in Section 3 also require a set of compositional partial derivatives. Assuming the scalar equation of state is formulated as a function of P and T (Section 3.5), the derivatives required to evaluate eqs (52), (53) and (54) are:

which can then be substituted into eqs (55), (56) and (57). The first three partial derivatives are defined by the scalar equation of state; only the last partial derivative is a function of the anisotropic tensor |$\mathbb {\Psi }$|⁠. Using eq. (16)

For orthotropic materials (where the eigenvectors of the terms in eq. 72 are equal to each other):

For non-orthotropic materials, eq. (72) must be evaluated numerically, rather than analytically.

4.3 Step-by-step guide to calculating relaxed properties

The formalism of the anisotropic solution model Myhill (2024) consists of appending an anisotropic tensor (⁠|$\mathbb {\Psi }(V^{\text{mol}}, T, \boldsymbol {x}, \boldsymbol {q})$|⁠) to any scalar hydrostatic equation of state (e.g. |$\mathcal {G}^{\text{hyd}}(P, T, \boldsymbol {x}, \boldsymbol {q})$|⁠). This is particularly convenient for Landau-type models, because the scalar model can be used to calculate equilibrium |$\boldsymbol {q}^{*}$| values at any state (Section 2.2); it is unnecessary to solve a system of equations to determine the symmetry-breaking strain (Section 2.3).

Given a previously constructed anisotropic model (e.g. |$\mathcal {G}(P, T, \boldsymbol {x}, \boldsymbol {q})$| and |$\mathbb {\Psi }(V, T, \boldsymbol {x}, \boldsymbol {q})$|⁠), the relaxed properties of any orthotropic material can be obtained at any values of the natural variables P, T and |$\boldsymbol {x}$| by the following steps:

• Find the value of |$\boldsymbol {q}$| that minimizes the Gibbs energy of the scalar equation of state at fixed P, T and |$\boldsymbol {x}$|⁠.

• Calculate the unrelaxed scalar properties of the material at the state given by (i) :

  $$\begin{eqnarray} && V = \left( \frac{\partial \mathcal {G}}{\partial P} \right)_{T, \boldsymbol {x}, \boldsymbol {q}} \text{, } S = -\left( \frac{\partial \mathcal {G}}{\partial T} \right)_{P, \boldsymbol {x}, \boldsymbol {q}} \text{, }\\ && \beta _{\text{RT}} = -\left( \frac{\partial \ln V}{\partial P} \right)_{T, \boldsymbol {x}, \boldsymbol {q}} \text{, } \alpha _{V} = \left( \frac{\partial \ln V}{\partial T} \right)_{P, \boldsymbol {x}, \boldsymbol {q}} \text{, } C_P = \left( \frac{\partial S}{\partial \ln T} \right)_{P, \boldsymbol {x}, \boldsymbol {q}} \text{, } \\ && \left( \frac{\partial ^2 \mathcal {G}}{\partial q_i \partial q_j} \right)_{P, T, \boldsymbol {x}} \text{, } \left( \frac{\partial S}{\partial q_i} \right)_{P, T, \boldsymbol {x}} \text{ and } \left( \frac{\partial V}{\partial q_i} \right)_{P, T, \boldsymbol {x}} \text{.} \end{eqnarray}$$

• Calculate the value of |$\mathbb {\Psi }(V^{\text{mol}}, T, \boldsymbol {x}, \boldsymbol {q})$| (eq. 63) and its derivatives with respect to |$V^{\text{mol}}$|⁠, T and |$\boldsymbol {q}$|⁠.

• From the derivatives determined in (ii) and (iii), calculate the unrelaxed thermal expansivity |$\boldsymbol {\alpha }$|⁠, isothermal compliance tensor |$\mathbb {S}_{\text{T}}$| and |$(\partial \boldsymbol {\varepsilon } / \partial \boldsymbol {q})_{V, \bar{\boldsymbol {\varepsilon }}, T}$| (eqs 69, 71 and 75).

• Calculate the unrelaxed isothermal stiffness tensor |$\mathbb {C}_T$|⁠, thermal stress tensor |$\boldsymbol {\pi }$| and isometric heat capacity density |$C^{\rho }_{\boldsymbol {\varepsilon }}$| (eqs 21, 23 and 24) and use them to construct the unrelaxed block tensor |$\left(\partial ^2 \mathcal {F}^{\rho } / \partial z_i \partial z_j\right)_{\boldsymbol {q}, \boldsymbol {x}}$| (first term on the RHS of eq. 36).

• Calculate |$\left(\partial ^2 \mathcal {F}^{\rho } / \partial q_i \partial z_j\right)_{\boldsymbol {\varepsilon }, T, \boldsymbol {x}}$| and |$\left(\partial ^2 \mathcal {F}^{\rho } / \partial q_i \partial q_j\right)_{\boldsymbol {x}}$| using eqs (52), (53) and (54) via eqs (55), (56) and (57).

• Compute |$\left(\partial q^{*}_k / \partial z_j\right)_{\boldsymbol {x}}$| from eqs (33) and (34).

• Calculate the relaxed anisotropic properties by substitution into eq. (36).

For the specific example of stishovite in Section 5, Step (i) is achieved using eq. (78), Step (ii) through the derivatives defined by the Stixrude & Lithgow-Bertelloni (2005) equation of state, and Step (iii) through eqs (80) and (82).

5 EXAMPLE: THE STISHOVITE TO POST-STISHOVITE TRANSITION

5.1 Overview

Stishovite is a high-pressure mineral whose composition in natural rocks is close to SiO|$_{2}$|⁠. It is stable in mafic and felsic rocks under the P–T conditions of the mantle transition zone and most of the lower mantle (Murakami et al. 2003; Grocholski et al. 2013; Sun et al. 2019). It has abundances of around 10 to 25 per cent by weight of rocks with basaltic compositions (Hirose & Fei 2002; Ricolleau et al. 2010; Ishii et al. 2022). It has also been identified in meteorites and their shocked target rocks (Chao et al. 1962).

Stishovite experiences a weakly first-order displacive phase transition as pressure is increased (Kingma et al. 1995; Andrault et al. 1998; Hemley et al. 2000). This symmetry-breaking transition transforms stishovite from the tetragonal rutile structure (space group P4|$_2$|/mnm) to the orthorhombic CaCl|$_{2}$| structure (space group Pnnm) by rotating octahedral SiO|$_{6}$| groups around an axis parallel to the c-axis of the unit cell. Experimentally derived estimates of transition pressures range between 40–60 GPa at room temperature (Andrault et al. 1998, 2003; Hemley et al. 2000; Grocholski et al. 2013; Buchen et al. 2018; Zhang et al. 2021, 2023), and around 70–80 GPa at 2200 K (Shieh et al. 2005; Nomura et al. 2010; Fischer et al. 2018), which corresponds to a depth of 1800 km depth along a model mantle geotherm (Brown & Shankland 1981). The variation in transition pressure estimates may arise from factors including: differences in pressure standards (Andrault et al. 1998, 2003), sample and grain scale deviatoric stresses (Andrault et al. 2003; Buchen et al. 2018; Wang et al. 2023), chemical contamination (see effects of H and Al reported by Criniti et al. 2023), and hysteresis caused by the weak first-order character of the transition (Hemley et al. 2000).

In the stishovite to post-stishovite transition, either the a- or the b-axis becomes longer than the other. At the transition, the amount of energy required to change post-stishovite (⁠|$Q = \delta$|⁠) to stishovite (⁠|$Q = 0$|⁠) and to “anti”-post-stishovite (⁠|$Q = -\delta$|⁠) is equal to zero. This causes the isothermal or isentropic sum of stiffnesses |$\mathbb {C}_{11}$| + |$\mathbb {C}_{22}$| − 2|$\mathbb {C}_{12}$| to drop to zero.

Several experimental studies have provided data on stishovite and the stishovite-to-post-stishovite transition. This includes unit cell constraints at high temperature (Ito et al. 1974) and high pressure (Table 2 and Supporting Information), measurements of predictions of elastic moduli (Table 3), and also direct seismic velocity measurements on polycrystalline stishovite (e.g. Chen et al. 2024). The elastic tensor data can be used to calculate unit cell compressibilities and bulk moduli for direct comparison with the high-pressure unit cell data (Table 4).

The bulk moduli calculated directly from the variation of the volume with pressure are typically smaller than the bulk moduli calculated from the variation of the a and b axes (Table 2). These differences arise from incompatibilities between the functional forms of the volumetric and axis length models used to calculate the two bulk moduli, and suggest that the model-free uncertainty on the bulk modulus of stishovite is significantly larger than the formal uncertainties published in any given paper. Overall, the various estimates of the (Reuss) bulk modulus of stishovite imply a true value somewhere between 290 and 320 GPa, with the c-axis being somewhere between 0.4 and 0.6 times as compressible as the a-axis (Tables 2 and 4).

5.2 Model formulation

Papers modelling the elastic anisotropy of the stishovite to post-stishovite transition (Carpenter et al. 2000; Hemley et al. 2000; Buchen et al. 2018; Buchen 2021; Zhang et al. 2021, 2023) have previously used the athermal Landau formalism introduced by Carpenter & Salje (1998). Here, I use a variety of experimental data to parametrize a fully thermal model using the theory described in Sections 3 and 4.

5.2.1 Scalar model

Stishovite and post-stishovite can be modelled as a solid solution of post-stishovite (⁠|$Q=1$|⁠) and anti-post-stishovite (⁠|$Q=-1$|⁠), where both endmembers have the same properties, but with a- and b-axes exchanged. The high symmetry stishovite structure corresponds to an equimolar mix of these two endmembers. Although the displacive transition from stishovite to post-stishovite is strictly first order (Hemley et al. 2000), it can be modelled as second order using a ‘2-4’ Landau-type model (Carpenter et al. 2000), where the ‘2’ and ‘4’ indicate that the excess energy of the solution relative to |$\mathcal {G}(Q=0)$| can be modelled as the sum of a |$P-T$| dependent quadratic and constant quartic terms. Writing the excess energy |$\mathcal {G}_{\text{xs}}$| relative to the energy of the endmembers (⁠|$\mathcal {G}(P, T, Q=\pm 1)$|⁠):

The equilibrium value of Q is:

The excess Gibbs energy can be rewritten as a function of the proportions of post-stishovite (⁠|$p_{1}$|⁠) and anti-post-stishovite (⁠|$p_{-1}$|⁠):

In this study, I use the equation of state of (Stixrude & Lithgow-Bertelloni 2005) to describe the thermodynamic properties of the |$Q=0$| and |$Q=1$| states, with parameters fit to experimental data (Table 5). For the |$Q = 0$| endmember, all equation of state parameters were optimized apart from the heat capacity-controlling Debye temperature, which was fixed to the value suggested by Stixrude & Lithgow-Bertelloni (2022). For the |$Q = 1$| endmember, the standard state Helmholtz energy, volume and Debye temperature were allowed to differ from the values given for the |$Q = 0$| endmember. These three parameters represent the first-order controls on the position and slope of the transition. The absolute value of Q is arbitrary, so the Helmholtz energy difference was scaled so that a value of |$Q = 1$| lies within the pressure range of the data. The resulting parameters produce a value of |$\Delta \mathcal {G}(P_{\text{tr}}, T_{\text{tr}}) = -593$| J mol⁻¹. The Gibbs energy as a function of Q at different pressures is shown in Fig. 1. The modelled evolution of the equilibrium value of Q is also shown. It is strongly correlated with the rotation angle of SiO|$_{6}$| octahedra about the c-axis (Zhang et al. 2023; Fig. 2).

5.2.2 Anisotropic model

The unrelaxed anisotropic properties of post-stishovite (⁠|$Q = 1$|⁠) are modelled via a fourth order tensor |$\mathbb {\Psi }$| (Section 4). After some initial curve fitting to find a functional form for c-axis compressibility using pyeq3 (Phillips & Myhill 2022), the following formulation was chosen for the elements of |$\mathbb {\Psi }$| (given in Voigt form):

which results in the following definition of the unrelaxed compliance tensor:

This is an attractive form for the compliance tensor, because it implies that ratios between compliances converge on constant values in both infinite compression and tension, making it easier to ensure elastic stability via Born’s stability criteria (Born 1940; Mouhat & Coudert 2014).

The |$\mathbb {\Psi }$| tensor of anti-post-stishovite (⁠|$Q = -1$|⁠) is identical to those of post-stishovite (⁠|$Q = 1$|⁠), but with the first and second rows and columns exchanged, and the fourth and fifth rows and columns exchanged. The |$\mathbb {\Psi }$| tensor for the solution is treated as linear in Q:

Therefore, the tensor for high-symmetry (tetragonal) stishovite is given by:

This study assumes that |$\mathbb {\Psi }$| is independent of temperature. This choice is made for pragmatic reasons; there is not yet experimental elastic data for post-stishovite at high temperature. Experimental justification for the temperature-independence of |$\mathbb {\Psi }$| includes the observation that the length of the c axis appears to be dominantly a function of volume (Fig. 3).

To facilitate data inversion, three further functions are defined (for |$i = \lbrace 1, 2, 3\rbrace$|⁠):

These parameters can be fit using high pressure unit-cell data, so these functions facilitate separating the fitting procedure into distinct volume, unit cell and elastic inversions (Section 5.3). If the diagonal elements of |$\mathbb {\Psi }$| are also fit to the experimental data, the off-diagonals can be determined by expressions:

5.3 Data processing and inversion

5.3.1 Data pre-processing, constraints and uncertainties

Before running the inversion, the reported pressures of the experiments were adjusted by converting from the pressure calibrants used in the original studies (including these references not cited elsewhere in the paper: Olinger & Jamieson 1970; Decker 1971; Mao et al. 1978; Hazen & Finger 1982; Jamieson 1982; Anderson et al. 1989; Holmes et al. 1989; Angel et al. 1997; Sata et al. 2002; Tsuchiya 2003; Dewaele et al. 2004, 2008, see experimental table in Supporting Information for more details) to a single consistent set of pressure standards (Fei et al. 2007). The one exception to this are the pressures reported from the room temperature elastic study of Zhang et al. (2021), which are internally corrected using the density and elastic moduli from that paper, using the equation:

after first converting from |$\mathbb {S}_{\text{S}}$| to |$\mathbb {S}_{\text{T}}$| (eq. 20, a very small change for stishovite at room temperature). Room-temperature molar volumes, axis lengths and axis compressibilities are shown in Fig. 4.

After pre-processing, the data was inspected for internal consistency. This led to the following selection of data:

• Volume data from Andrault et al. (2003), Nishihara et al. (2005), Wang et al. (2012) and Fischer et al. (2018). Given the very small uncertainties in the reported volumes, the weighted misfits for the high pressure data were calculated as |$(P_{\text{obs}} - P_{\text{model}})/P_{\text{unc}}$|⁠. Pressure uncertainties |$P_{\text{unc}}$| were chosen to be equal to either the reported uncertainty, or |$0.1 + 0.01 \cdot P$| [GPa], whichever was the bigger. This downweighted data points with small pressure uncertainties.

• Relative 1 bar volume data from Ito et al. (1974). Relative volumes were used rather than absolute volumes because of a large discrepancy in volume at room temperature. Weighted misfits |$(V_{\text{rel,obs}} - V_{\text{rel,model}})/V_{\text{rel,unc}}$| were obtained using the reported volume uncertainties.

• Data on a-, b- and c-axis lengths were taken from Andrault et al. (2003) and Zhang et al. (2021). Uncertainties in the axis lengths were taken to be equal to either the reported uncertainty, or 0.05 per cent of the value of the axis length, whichever was larger. This downweighted data with unusually small uncertainties. Data from other sources was not included because it either overlapped the selected studies, or because (in the case of Fischer et al. 2018), the difference between a- and b-axis lengths was in poor agreement with other studies (see discussion in Section 5.4 and Supporting Information)

• Elastic moduli were taken from Jiang et al. (2009) and Zhang et al. (2021). Only some of the data from Zhang et al. (2021) was used, because the linear compressibilities in that study are not consistent with the corresponding components of the compliance tensor (Fig. 4, two right hand panels). For this reason, we limit the fitting of the elastic data to the following data:

  • All elastic moduli of Jiang et al. (2009) with nominal 10 GPa uncertainties on each elastic modulus.

  • Elastic moduli of Zhang et al. (2021) apart from |$\mathbb {C}_{\text{S}33}$|⁠, and their reported uncertainties, up to 25 GPa. |$\mathbb {C}_{\text{S}33}$| was excluded because an initial fit to all of the data showed that |$\mathbb {C}_{\text{S}33}$| in particular was inconsistent with other data, and much more closely matched the trend from Jiang et al. (2009).

  • The size of the difference between |$\mathbb {C}_{\text{S}11}$|–|$\mathbb {C}_{\text{S}22}$|⁠, |$\mathbb {C}_{\text{S}13}$|–|$\mathbb {C}_{\text{S}23}$| and |$\mathbb {C}_{\text{S}44}$|–|$\mathbb {C}_{\text{S}55}$| from Zhang et al. (2021) in the post-stishovite field.

• The isentropic bulk modulus and shear modulus from the study of Chen et al. (2024) are compared to the model Voigt–Reuss–Hill values, calculating a weighted misfit using the reported uncertainties from the paper.

• The transition pressure at 3000 K was assigned a value of 90 |$\pm$| 1 GPa, in accordance with the results of Fischer et al. (2018). See Section 5.4 for a discussion of this constraint.

5.3.2 Inversion procedure

Weighted least squares inversion for the anisotropic properties of stishovite was split into the following stages:

• Initial fitting of the |$\mathcal {G}(P, T)$| equation of state

• Initial fitting of the unit cell properties |$\boldsymbol {M}(V, T, q)$|

• Initial fitting of the relaxed elastic moduli |$\mathbb {C}_{\text{S}}(V, T, q)$|

• Joint inversion of all the experimental data to refine the initial fit parameters.

This piecewise refinement improved convergence performance. The optimized model parameters obtained from the final inversion are given in Tables 5 and 6. They are also provided in the BurnMan open source software project (Cottaar et al. 2014; Myhill et al. 2023, 2024), along with the complete fitting and plotting scripts used in the writing of this paper.

5.4 Inversion results

5.4.1 Volumes, unit cells and elastic moduli

Modelled and observed volumes are shown in Fig. 5. Excess volumes due to the transition are less than 1 per cent of the total volume over the pressure range of the mantle, and the transition cannot be easily be identified in a |$V(P)$| plot. Temperature suppresses the transition to higher pressures, and also to lower volumes, an observation that is entirely dependent on the data of Fischer et al. (2018). Transition pressures increase from 50 GPa at room temperature to 59, 74 and 90 GPa at 1000, 2000 and 3000 K respectively.

Modelled and observed a-, b- and c-axis lengths are shown in Fig. 6. The optimized model is in good agreement with the experimental data in the stishovite field and in the post-stishovite field at room temperature except in the immediate vicinity of the transition where some hysteresis is apparent. The difference in length of the a and b axes is much greater in the model than in the high temperature data from Fischer et al. (2018). The small difference between the a- and b-axis lengths reported by Fischer et al. (2018) is striking; at 1200 K their highest pressure datum is around 23 GPa higher than the stishovite-post-stishovite transition, yet |$\ln (b/a) \sim 0.012$|⁠; over three times smaller than the equivalent point at room temperature (⁠|$\ln (b/a) \sim 0.04$|⁠; see Supporting Information Figure). Furthermore, |$\ln (b/a)$| reported by Fischer et al. (2018) does not seem to be strongly dependent on proximity to the transition. Fischer et al. (2018) argued that the small splitting could be due to the more hydrostatic conditions within their high temperature experiments, but this explanation is unsatisfying given that the room temperature unit cell behaviour of post-stishovite is similar across many different studies (Fig. 3). Alternative explanations for the small splitting observed by Fischer et al. (2018) include a moderate pressure overestimation in the post-stishovite field, or intragranular stresses within the diamond anvil cell (see also Wang et al. 2023). In the first case, the model transition curve may be accurate. In the second, the equilibrium transition might be at higher pressure than modelled. If the transition pressure constraint (Section 5.3.1) is removed, the best fit to the data is obtained by increasing the Debye temperature by |$\sim$|10 K relative to the value in Table 5.

Modelled and observed room temperature elastic moduli are shown in Fig. 7. The optimized model generally agrees well with the published elastic moduli (Jiang et al. 2009; Zhang et al. 2021), with both a good fit to zero pressure values and to their pressure gradients. Model values for |$\mathbb {C}_{\text{S}33}$| are in good agreement with the data of Jiang et al. (2009), but in poor agreement with Zhang et al. (2021). I note in passing that the high pressure (Zhang et al. 2021) moduli were apparently derived from Brillouin data on single platelets, which may plausibly cause systematic errors in estimates of elastic moduli.

The inversion in this paper involves fitting 29 equation-of-state parameters. Eight of these are for the scalar model that determines the equilibrium value of q (Table 5). Seven more parameters describe the standard state unit-cell of post-stishovite (⁠|$Q=1$|⁠) and its evolution with respect to volume, and the other fourteen describe the other elastic properties (Table 6). The total number of parameters in this model is comparable to the 27 parameters in the mutually independent isothermal |$V(P)$| and |$\mathbb {C}(P)$| models of (Zhang et al. 2021).

5.4.2 Seismic velocities

The optimized model was used to calculate seismic velocities, linear compressibilities and extrema in Poisson’s ratios as a function of crystal orientation 0.1 GPa below (Fig. 8) and 0.1 GPa above (Fig. 9) the stishovite to post-stishovite transition. As expected, |$V_{S2}$| velocities approach zero along the |$[1, 1, 0]$| direction from both sides of the transition. The linear compressibility is positive in all directions below the transition, but above the transition, the model predicts that linear compressibility along the longest axis (b) is negative; squeezing the structure hydrostatically results in expansion of the b-axis. This negative compressibility can also be seen in the unit cell properties (Fig. 4). It persists until about 10 GPa above the transition.

Another interesting observation from the model is that near the transition, the minimum Poisson ratio in almost every direction is negative; that is, when stretched along any direction, crystals of stishovite and post-stishovite also get longer in an orthogonal direction (Figs 8 and 9). This auxetic behaviour decreases in intensity away from the transition. It is increasingly restricted to particular orientations further away from the transition; there are almost no directions of auxeticity 40 GPa below the transition, and auxeticity more than 10 GPa above the transition is restricted to a broad, persistent band surrounding the plane perpendicular to the b-axis. A pressure sequence of figures similar to Figs 8 and 9 is provided as Supporting Information.

Finally, Fig. 10 shows the isotropic velocities |$V_P$|⁠, |$V_{\Phi }$| and |$V_S$| calculated from the Reuss, Voigt and Voigt–Reuss–Hill (VRH) isentropic bulk and shear moduli. Note the small discontinuity in the VRH compressional and bulk sound velocities at the transition that results from the discontinuous change in the isentropic Reuss bulk modulus. The prominent reduction in |$V_S$| associated with the transition has previously been associated with the observation of seismic scatterers in the lower mantle (Kaneshima & Helffrich 1999; Niu 2014; Kaneshima 2019).

6 CONCLUSIONS

The current study shows that the anisotropic solution model of Myhill (2024) can be used to model materials undergoing displacive phase transitions. Using a comparable number of parameters to existing isothermal treatments, the equation of state provides a full, internally consistent set of thermodynamic and elastic properties at high temperature and pressure. The model retains thermodynamic consistency even when symmetry-breaking strains are large (eq. 82).

The new anisotropic solution model can be built on any traditional scalar model |$\mathcal {G}(P, T, \boldsymbol {x})$|⁠, allowing elastic models to build on pre-existing data sets. In the context of displacive phase transitions, the method proposed in this paper can be used with existing Bragg & Williams (1934)-type models that have been incorporated into large thermodynamic data sets (e.g. Holland & Powell 2011; Stixrude & Lithgow-Bertelloni 2024). Given the relatively small number of parameters required to describe the stishovite model in this paper, it should be possible to extend the method used here to vectors of compositional variables |$\boldsymbol {x}$| and internal variables |$\boldsymbol {q}$|⁠.

This paper highlights the need to simultaneously consider volumetric, unit cell and elastic data in order to produce internally consistent data sets and models. Further experimental work to determine the unit cell dimensions and elastic moduli of post-stishovite at different pressures, temperatures and compositions would be most welcome.

ACKNOWLEDGMENTS

Paul Asimow introduced me to seismic relaxation many years ago at the 2014 CIDER workshop. Thanks to him, all the participants, and the National Science Foundation for providing such a fantastic environment for learning and research.

Many thanks also to three anonymous reviewers, who provided constructive criticism of earlier versions of this paper.

This work was supported by Natural Environment Research Council Large Grant MC-squared (Award No. NE/T012633/1) and Science and Technology Facilities Council (Grant No. ST/R001332/1). Any mistakes or oversights are my own.

DATA AVAILABILITY

The anisotropic equation of state described in this paper has been contributed to the BurnMan open source software project: https://github.com/geodynamics/burnman (Cottaar et al. 2014; Myhill et al. 2023, 2024).

REFERENCES

APPENDIX A: DERIVATION OF THERMODYNAMIC RELATIONS INVOLVING CONJUGATE PROPERTIES

A1 Isentropic and isothermal stiffness tensors

where eq. (A3) uses the triple product rule, eq. (A4) uses the symmetry of double derivatives,

and eq. (A6) takes the derivative with respect to strain holding the entropy constant.

A2 Thermal expansivity, isothermal compliance and thermal stress tensors

By the triple product rule:

A3 Isostress and isostrain heat capacities

where eq. (A12) uses the symmetry of double derivatives,

Eq. (A14) uses the triple product rule, and eq. (A16) takes the derivative with respect to temperature holding the stress constant.

APPENDIX B: CALCULATING RELAXED THERMODYNAMIC PROPERTIES

Relaxed thermodynamic properties are determined by minimizing an appropriate thermodynamic potential, allowing some or all isochemical structural variables to vary freely. Similar mathematical derivations can be used to quantify the effective isotropic thermodynamics of multiphase assemblages by varying the pressure and temperature (Stixrude & Lithgow-Bertelloni 2022), or consider the anisotropic properties of a single phase, as we do here.

Let us seek expressions for the relaxed second order derivatives of the Helmholtz energy density. We first define a set of functions |$P_l(\boldsymbol {M},T,\boldsymbol {q})$| that describe the partial derivatives of the Helmholtz energy density with respect to |$\boldsymbol {q}$|⁠:

where |$\boldsymbol {M}$| is the extensive metric tensor, with elements having units of m, and |$\boldsymbol {q}$| are the unitless isochemical structural vectors. The minimization of |$\mathcal {F}$| is achieved when

The relaxed properties |$\mathbb {C}^{*}_{\text{T}}$|⁠, |$\boldsymbol {\pi }^{*}$| and |$c^{*}_{\varepsilon }$| depend on the way in which |$\boldsymbol {q}^{*}$| varies due to small changes in |$\boldsymbol {\varepsilon }$| or T (here collected as |$\boldsymbol {z} = \lbrace \boldsymbol {\varepsilon }, T\rbrace$|⁠). This dependence can be determined by applying the chain rule:

and solving for |$\partial q^{*}_m / z_j$|⁠:

where R is the left inverse matrix:

The relaxed physical properties evaluated at constant T are then determined by repeated application of the chain rule:

where the second term in the first expression vanished by the requirement that |$\partial \mathcal {F}^{\rho }/\partial q_k = 0_k$| at equilibrium (eq. B2). Differentiating again yields

This expression is similar to the expression under fixed entropy constraints presented by Slonczewski & Thomas (1970), their eq. (26). The elements of the unrelaxed block matrix (first term on the RHS of the last equation) can be evaluated directly from the anisotropic equation of state. Expressions for |$\partial ^2 \mathcal {F}^{\rho } / \partial \boldsymbol {q} \partial \boldsymbol {q}$| (eq. B8) and |$\partial ^2 \mathcal {F}^{\rho } / \partial \boldsymbol {q} \partial \boldsymbol {z}$| (eq. B13) can be derived by change of variables (Appendix  C).

APPENDIX C: CHANGE OF VARIABLES AND DERIVATIVES OF THE HELMHOLTZ ENERGY

Here, we derive the partial derivatives in Section 3.5 by change of variables. The variables required are:

where |$\boldsymbol {M_{\text{ref}}}$| is a constant reference cell tensor about which small strains can be performed. The equation of state developed by Myhill (2024) gives the partial derivatives of |$\mathcal {F}(f^{\prime },T^{\prime },\boldsymbol {n}^{\prime },\bar{\boldsymbol {\varepsilon }}^{\prime })$|⁠, where |$f = \ln V^{\text{mol}}$|⁠. The ‘primes’ are used to specify the variables of the equation of state. Making explicit the dependences on the variables in C1, we have:

The partial differential operators for the set of variables in C1 are thus:

In order to write the second derivatives in their most compact form, we first need to derive some of the partial derivatives in the three expressions above. The small strain identity is

Other derivatives are calculated considering the total derivative of the small strain tensor |$\boldsymbol {\varepsilon }$|⁠:

1. Taking the partial derivative with respect to |$\varepsilon _{kl}$| at constant T and |$\boldsymbol {q}$|⁠:

   $$\begin{eqnarray} \frac{\partial \bar{\varepsilon }^{\prime }_{ij}}{\partial \varepsilon _{kl}} = \frac{\partial \varepsilon _{ij}}{\partial \varepsilon _{kl}} - \frac{\partial \varepsilon _{ij}}{\partial f^{\prime }} \frac{\partial f^{\prime }}{\partial \varepsilon _{kl}} = \delta _{ik}\delta _{jl} - \frac{\beta _{\text{T}ij}}{\beta _{\text{TR}}}\delta _{kl} \end{eqnarray}$$

   (C7)
2. with respect to |$q_{i}$| at constant |$\boldsymbol {\varepsilon }$| (⁠|$=0$|⁠) and T:

   $$\begin{eqnarray} \frac{\partial \bar{\varepsilon }^{\prime }_{mn}}{\partial q_i} = -\frac{\partial \varepsilon _{mn}}{\partial q^{\prime }_i} \end{eqnarray}$$

   (C8)
3. with respect to T at constant |$\boldsymbol {\varepsilon }$| (⁠|$=0$|⁠) and |$\boldsymbol {q}$|⁠:

   $$\begin{eqnarray} \frac{\partial \bar{\varepsilon }^{\prime }_{mn}}{\partial T} = -\frac{\partial \varepsilon _{mn}}{\partial T^{\prime }} = -\left(\frac{\alpha _{ij}}{\alpha _V} - \frac{\beta _{\text{T}ij}}{\beta _{\text{TR}}}\right) \alpha _V. \end{eqnarray}$$

   (C9)

We also have, for a hydrostatic reference state:

Using these expressions, we can compactly write the second derivatives of interest:

APPENDIX D: CHANGE OF VARIABLES: SCALAR GIBBS DERIVATIVES TO HELMHOLTZ DERIVATIVES

This appendix provides partial differential operators and partial derivatives required to calculate the second order partial differentials of the Helmholtz energy |$\mathcal {F}$| under hydrostatic conditions. As done for the Helmholtz energy under non-hydrostatic conditions in Appendix  C, we define a set of variables for the hydrostatic Helmholtz energy (⁠|$V^{\prime },T^{\prime },\boldsymbol {n}^{\prime }$|⁠), and a modified set for the hydrostatic Gibbs energy (⁠|$P^{\mathcal {G}},T^{\mathcal {G}},\boldsymbol {n}^{\mathcal {G}}$|⁠). Any function can then be written:

The partial differential operators for the Helmholtz variables are:

In order to express the second derivatives in their most compact form, we will use the following thermodynamic identities:

A final relation is obtained by taking the total derivative of V and differentiating with respect to |$\boldsymbol {n}^{\prime }$|⁠:

Finally, the required second derivatives can be derived:

APPENDIX E: CHANGE OF VARIABLES: ANISOTROPIC HELMHOLTZ DERIVATIVES TO INTERNAL ENERGY DERIVATIVES

This appendix provides partial differential operators and partial derivatives required to calculate the second order partial differentials of the internal energy |$\mathcal {E}$| as functions of partial derivatives of the Helmholtz energy |$\mathcal {F}$|⁠. We define a ‘natural’ set of variables for the Helmholtz energy in the small strain case (⁠|$\boldsymbol {\varepsilon },T,\boldsymbol {n}$|⁠), and a modified set for the internal energy (⁠|$\boldsymbol {\varepsilon }^{\mathcal {E}},S^{\mathcal {E}},\boldsymbol {n}^{\mathcal {E}}$|⁠). Any function can be written

The partial differential operators for the modified set of variables are:

In order to express the second derivatives in their most compact form, we will use the following thermodynamic identity:

Further relations are obtained by differentiating S with respect to |$\boldsymbol {\varepsilon }^{\mathcal {E}}$|⁠, |$S^{\mathcal {E}}$| and |$\boldsymbol {x}^{\mathcal {E}}$|⁠:

Finally, the required second derivatives can be derived: