An anisotropic equation of state for solid solutions, with application to plagioclase

SUMMARY

This paper presents a framework for building anisotropic equations of state for solid solutions. The framework satisfies the connections between elastic and thermodynamic properties required by Maxwell’s relations. It builds on a recent anisotropic equation of state for pure phases under small deviatoric stresses, adding a dependence on a vector $\mathit{n}$⁠, whose components $n_i$ contain the molar amounts of independent end-members in the solid solution. These end-members may have distinct chemical compositions, site species occupancies or electronic spin states. The high albite-anorthite (C$\overline{\text{1}}$⁠) plagioclase solid solution is used to illustrate the formulation.

1 INTRODUCTION

Solid solutions can be thought of as mixtures of different end-members, each of which have fixed structure, composition and distribution of species on sites. In traditional hydrostatic thermodynamic models, the Gibbs (or Helmholtz) energy is expressed as a function of pressure (or volume), temperature and end-member proportions (e.g. Helffrich & Wood 1989; Holland & Powell 2003; Stixrude & Lithgow-Bertelloni 2005; Myhill & Connolly 2021). Derivatives of the energy with respect to these variables yield physical properties including the volume (or pressure), entropy, thermal expansivity, Reuss bulk modulus and isobaric and hydrostatic–isochoric heat capacities.

Partial derivatives of the Gibbs (or Helmholtz) energy with respect to pressure (or volume) cannot be used to obtain anisotropic physical properties such as the thermal expansivity tensor and elastic stiffness tensor. Calculating these properties requires that partial derivatives of the energies be taken as functions of stress or strain. In Myhill (2022), I presented an anisotropic equation of state for pure, isochemical substances. The equation of state was designed so that it could be used in conjunction with any traditional hydrostatic equation of state, with the new parameters defining only the anisotropic properties. In this paper, I extend this equation of state to solid solutions. Symbols used in this paper are given in Table 1.

2 CELL TENSORS, DEFORMATION AND REFERENCE FRAMES

2.1 The standard state cell tensor

At a given reference state (usually 1 bar, 298.15 K), the unit cell of each end-member in a solid solution can be defined using vector lengths and angles (a, b, c, $\alpha$⁠, $\beta$ and $\gamma$⁠). These are reference frame-invariant - that is, none of the values depend on the orientation of the crystal relative to the frame of reference. However, a frame of reference must be defined when considering deformation of the unit cell in an equation of state, so it is convenient to also define the unit cell relative to that frame of reference. In the feldspar example used in this study, we adopt the convention of Brown et al. (2016), where the y-axis of the Cartesian reference frame is aligned parallel to the crystallographic b-axis in the standard state, the x-axis is perpendicular to the b- and c-axes, and the z-axis is chosen to satisfy a right-handed coordinate system. Using this convention, we can define a reference state molar ‘cell tensor’ ${\mathit{M}}_0$ [m mol⁻¹]:

where $N_{\text{A}}$ is Avogadro’s constant and Z is the number of unit cells per formula unit. Each column of this tensor represents a vector from the origin to an adjacent vertex of a parallelepiped. This parallelepiped has the shape of the unit cell and a volume equal to the molar volume of the material $V_0$⁠.

At standard temperature and pressure, triclinic albite (NaAlSi${}_3$O${}_8$⁠, C$\bar{\text{1}}$⁠) has unit cell parameters $a=8.1366(2)\cdot {10}^{-10}$ m, $b=12.7857(2)\cdot {10}^{-10}$ m, $c=7.1582(3)\cdot {10}^{-10}$ m, $\alpha =94.253(2{)}^{\circ }$⁠, $\beta =116.605(2{)}^{\circ }$⁠, $\gamma =87.756(2{)}^{\circ }$ and $Z=4$ (Brown et al. 2016). These values lead to the following reference state cell tensor:

yielding a molar volume $V_0=det({\mathit{M}}_0)=99.965$ cm³ mol⁻¹.

2.2 The deformation gradient tensor and its derivatives

Deformation of a cell tensor from ${\mathit{M}}_0$ to a state $\mathit{M}$ is achieved by applying a deformation gradient tensor $\mathit{F}$ to the reference state cell tensor ${\mathit{M}}_{\mathbf{0}}$

The deformation gradient tensor $\mathit{F}$ and its time derivative $\dot{\mathit{F}}$ describe the transformation of local particle positions from their initial state ${\mathit{x}}_0$ to a final state $\mathit{x}$⁠:

The velocity gradient tensor $\dot{\mathit{L}}$ represents the change in velocity of particles relative to their current positions:

The velocity gradient tensor can be asymmetric even if the deformation gradient tensor is always symmetric (Section 2.3). Infinitesimal strain rate $\dot{\mathit{\epsilon }}$ and spin $\dot{\mathit{\omega }}$ tensors are defined as the symmetric and antisymmetric parts of $\dot{\mathit{L}}$⁠:

The infinitesimal strain rate can be decomposed into temperature and stress-related terms:

where $\mathit{\alpha }$⁠, ${\mathit{\beta }}_T$ and ${\mathbb{S}}_T$ are the thermal expansivity, isothermal compressibility and elastic compliance tensors. Stress $\mathit{\sigma }$ is positive under tension, pressure P is positive under compression and $\mathit{\tau }$ is the deviatoric stress:

Using eqs (6), (7) and (9), the thermal expansivity and isothermal compressibility can be written

where ${\beta }_{\text{RT}}$ is the isothermal Reuss compressibility. Orthotropic materials deformed under hydrostatic conditions in a rotation-free coordinate frame will have eigenvectors of $\mathit{F}$ which are constant with respect to pressure and temperature. In such cases, $\dot{\mathit{F}}$ and ${\mathit{F}}^{-1}$ are commutative, the velocity gradient tensor can be written (Haber 2018) ${\dot{L}}_{ij}={\dot{({\ln }_{\text{M}}\mathit{F})}}_{ik}$⁠, and the expressions for $\mathit{\alpha }$ and ${\mathit{\beta }}_T$ simplify considerably (Myhill 2022):

2.3 Rotation in non-orthotropic systems

Under hydrostatic conditions, conservation of angular momentum implies that any infinitesimal deformation will be rotation-free, and therefore that $\dot{\mathit{L}}$ will always be symmetric. Unfortunately, in non-orthotropic systems (monoclinic and triclinic), a symmetric $\dot{\mathit{L}}$ does not guarantee a symmetric deformation gradient tensor $\mathit{F}$⁠. Asymmetry, and the finite rotation implied by that asymmetry, arises when the eigenvectors of strain change during deformation. Fig. 1 illustrates this by means of a two-step pure-shear deformation of an initially square object. The incremental strain at each step is shown by an arrow and associated matrix, and the finite deformation gradients after each step are illustrated as a shape and matrix inside that shape. In Fig. 1(a), the box is first shortened along the y-axis, and then shortened along an axis inclined from the vertical. Note that even though both strain increments are rotation-free (i.e. symmetric), the second step produces an asymmetric deformation gradient tensor.

In the anisotropic equation of state presented by Myhill (2022), the deformation gradient tensor $\mathit{F}$ is calculated as a function of volume and temperature, and is always symmetric. This implies that increments in strain (⁠$\dot{\mathit{L}}$⁠) may be asymmetric. Fig. 1(b) illustrates the same deformation as that of Fig. 1(a), but with rotations now incorporated into the incremental strain steps, rather than the finite deformation gradient tensor. Thermal expansivity, isothermal compressibility and isothermal compliance tensors must be calculated taking these rotations into account (eqs 12 and 14 and Section 3.2.1).

3 FORMULATION

3.1 Volumetric equations of state

The anisotropic equation of state presented in this paper can be built on top of any equation of state for which molar volume can be found as a function of pressure P, temperature T and independent end-member proportions $\mathit{p}$⁠:

This can be done for almost any solution model, even those where volumes are variables in the model (e.g. Stixrude & Lithgow-Bertelloni, 2011; Myhill, 2018). In most solution models in the geological literature, the calculation of the volume is split into end-member and excess contributions:

where $p_i$ is the proportion of independent end-member i in the solution, $V_i$ is the volume of that end-member and $W_{ij}^V$ are volumetric interaction terms between end-members. End-member volumes may be determined directly from P and T, or the end-member equations of state may be formulated as a function of volume, in which case the correct volume at any given pressure may be found by iteration (e.g. Stixrude & Lithgow-Bertelloni, 2011; Myhill, 2018).

3.2 The anisotropic equation of state

3.2.1 End-members

In Myhill (2022), I showed that self-consistent anisotropic properties for pure phases could be modelled by using a fourth order symmetric anisotropic tensor $\mathbb{\Psi }(V,P_{\text{th}}(V,T))$⁠, with parameters that satisfied the condition:

In that paper, the anisotropic tensor was used to define the deformation gradient tensor and the isothermal compliance tensor:

The constraint given by eq. (20) ensures that the anisotropic equation of state remains consistent with the volumetric equation of state (Section 3.1), as a result of the mathematical identity (Petersen & Pedersen 2012):

and identity:

Myhill (2022) focused mainly on orthotropic materials (the equation of state was demonstrated by orthorhombic forsterite). For these materials, eqs (21) and (22) are consistent with each other because:

For non-orthotropic materials, this is no longer the case; rotation (eq. 13 and Section 2.3) means that:

For monoclinic and triclinic systems, eq. (22) can be replaced with eq. (26) and the following equations:

The indices in the last two equations signify different index values (e.g. $i=1$⁠, $j=3$⁠, $k=2$⁠); no summation is implied. The factor 2 in eq. (30) arises from the multiplication of individual compliances required in the conversion to the Voigt form of the compliance tensor. The modified terms occupy the positions coloured blue (eq. 29) and red (eq. 30) in the following compliance tensor:

3.2.2 Solid solutions

At any given composition, solid solutions must obey the same self-consistency rules as the end-members. One way to ensure this is to define the standard state molar cell tensor ${\mathit{M}}_0$ and deformation gradient tensor $\mathit{F}$ at any given composition as follows:

where $V_{\text{mol}}$ is the molar volume, n is the total number of moles of substance and $\mathit{p}$ is the vector of molar fractions of the independent species. The ${}_m$ subscript refers to the mth end-member in the solution. $\mathit{R}$ is a rotation matrix needed to rotate the deformed cell tensor so that it it is returned to the unit cell reference frame (Section 2.1). The following equality must be satisfied:

The reference molar volume and relative volume change at any composition can then be found from ${\mathit{M}}_0$ using the identities given by eqs (23) and (24):

which demonstrates that the decomposition is consistent with the volume of the scalar equation of state. The end-member deformation gradient tensors are already defined (Section 3.2.1) as

where the fourth rank tensor $\mathbb{\Psi }$ for each end-member m can be contracted into Voigt form and has the symmetry of the crystal structure (Fig. 2).

If we let $({\ln }_{\text{M}}\mathit{F}{)}_{\text{xs}}$ take a similar form

then we can express the deformation gradient at any composition as

The cell tensor at any molar volume and temperature can be found by combining eqs (32), (33) and (40). The thermal expansivity $\mathit{\alpha }$ and isothermal compressibility ${\mathit{\beta }}_{\text{T}}$ tensors at fixed composition can be found using eqs (12), (14) and (40), and the isothermal compliance tensor ${\mathbb{S}}_{\text{T}}$ can be calculated using eqs (27)–(30) and (41). Other anisotropic properties may be obtained from the expressions in Myhill (2022).

This formulation ensures that Maxwell’s anisotropic relations, which stem from the symmetry of mixed partial derivatives, are rigorously satisfied. For example:

4 APPLICATION TO THE PLAGIOCLASE FELDSPARS

4.1 Introduction

Plagioclase feldspars are a dominant mineral group in the crust. Most plagioclase feldspars are well-approximated by a mix of two different chemical components, an albitic component with composition NaAlSi${}_3$O${}_8$ and an anorthitic component with composition NaAlSi${}_3$O${}_8$⁠. They are triclinic, which makes them an ideal example solid solution to demonstrate the anisotropic equation of state introduced in this paper.

Plagioclase feldspars change space group depending on their composition, pressure and temperature of equilibration. End-member anorthite adopts the P$\overline{\text{1}}$ space group, anorthite-rich (⁠$p_{\text{an}}>0.5$⁠) plagioclases adopt the ordered I$\overline{\text{1}}$ space group, and albite-rich plagioclases adopt the disordered C$\overline{\text{1}}$ space group. The I$\overline{\text{1}}$ field shrinks relative to the C$\overline{\text{1}}$ field with increasing temperature but expands with increasing pressure, a consequence of the lower configurational entropy and molar volume of the I$\overline{\text{1}}$ phase.

Structural and elastic data for plagioclase crystals across the albite-anorthite binary have been collected by Angel et al. (1990) and Brown et al. (2006, 2016), and these data are used here to create an anisotropic model for the disordered C$\overline{\text{1}}$ (high) plagioclases. Extension to the I$\overline{\text{1}}$ space group is left for a future study, given the need to develop an anisotropic extension of Landau Theory consistent with the equation of state (Carpenter 1988; Dubacq 2022).

4.2 End-member and solution volumes

End-member properties for albite and anorthite were taken from the dataset of (Stixrude & Lithgow-Bertelloni 2022). These were then adjusted to best fit the data for C$\overline{\text{1}}$ plagioclases. As elastic properties across the solid solution have only been collected at room temperature and near-room pressure, only the $V_0$ and $K_{\text{T}0}$ parameters were modified.

The excess properties across the solid solution were modelled as a symmetric (regular) binary solution. A constant volume excess term cannot accurately reproduce both the variation of V and $K_{\text{TR}}$ across the solution, and so the volume excess was modelled by creating two new intermediate end-members, one (aban1) with $V_0$⁠, $K_0$ and $K_0^{\prime }$ equal to the average of the albite values, and a second (aban2) where $V_0$ and $K_0$ were allowed to vary. The volume was then modelled as:

The fitted end-member and mixing properties are provided in Table 2. Because aban1 is derived from the properties of the ab and an end-members, there are only three free variables for each end-member (six in total). The resulting standard state volumes and Reuss isothermal bulk moduli are plotted in Fig. 3.

4.3 Cell parameters

Anisotropic model parameters in this study are all treated as ideal; in other words, every element in the tensor ${\mathbb{W}}_{ijklmn}^{\mathbb{\Psi }}$ in eq. (41) is equal to zero. All the anisotropic properties across the binary are calculated using the end-member ${\mathit{M}}_0$ and $\mathbb{\Psi }$ functions.

The standard state molar cell tensor ${\mathit{M}}_0$ for each end-member is calculated from the molar cell parameters in Table 3 (see also Section 2.1). Because $V_0=det({\mathit{M}}_0)$⁠, there are only five free variables for each end-member (10 variables in total). The resulting cell properties across the binary are shown in Fig. 4.

4.4 Elastic properties

As the data used in the inversion in this study is all collected at room pressure and temperature, the simplest formulation for $\mathbb{\Psi }$ can be adopted, whereby:

The elements of the tensor $A_{mijkl}$ for each end-member m have the symmetries of an elastic compliance tensor, and so can be written in Voigt form (21 parameters for each end-member for triclinic symmetry). The only other requirement of this tensor is that $A_{mijkl}{\delta }_{ij}{\delta }_{kl}=1_m$ (eq. 20), leaving 20 free variables for the two end-members (40 variables in total). Values are provided in Tables 4 and 5.

4.5 Data inversion

Inversion of the data to obtain model parameters was performed in several parts:

• An approximation to the scalar $V(P,T)$ equation of state was obtained by fitting the parameters in Section 4.2 to the volume and Reuss isothermal bulk modulus data.

• Approximate end-member cell parameters were found (Section 4.3), fixing the scalar parameters found in the previous step.

• Approximate end-member elastic parameters were found (Section 4.4) by fitting to the elastic data. Because the relationship between the model parameters and the elements of the elastic stiffness tensor are highly non-linear, it was found to be much more efficient to initially fit the ratio of isentropic compliances to isentropic Reuss compressibilities (Fig. 5), which bear a near 1:1 relationship with the tensors A (eqs 27–30 and 44).

• Finally, all 56 parameters were simultaneously inverted using all the available data and uncertainties, including the isentropic elastic tensors (Fig. 6).

4.6 Observed data versus model predictions

Overall, the fit between the observed data and model predictions is extremely good. This is perhaps not surprising, given the large number of fitting parameters (56) versus the number of data points (116). However, the good fit does suggest that treating the mixing of the anisotropic tensor as ideal (as done here) is suitable even when the solutions are volumetrically non-ideal. In addition, the model does draw out some nice contrasts between the C$\overline{\text{1}}$ and I$\overline{\text{1}}$ properties; in particular, the kink in volume and cell parameters at the phase boundary (⁠$p_{\text{an}}\sim 0.5$⁠) and abrupt drop in isothermal bulk modulus (Figs 3 and 4). Also of interest are the near linear trends in $S_{\text{N}ij}/{\beta }_{\text{SR}}$ (Fig. 5). Linearity in these trends (at constant V and T) is a prediction of the ideal anisotropic model, and so it is comforting to see that the simple model formalism does a good job at representing the data.

The second order isothermal compressibility tensors reported by (Angel 2004) and Brown et al. (2016) were not used in the creation of the anisotropic model. This is because the isothermal compressibility tensor is very closely related to the isothermal elastic tensor (eq. 25), which in turn is closely related to the isothermal elastic tensor (at 0 K the two are identical). A comparison between the isothermal compressibilities reported by Brown (2018) and the model predictions are presented in Fig. 7. Note that while the values of the compressibilities are reasonable, the trends of ${\beta }_{\text{T}2}$⁠, ${\beta }_{\text{T}3}$ and ${\beta }_{\text{T}6}$ are in poor agreement with the data, suggesting minor conflict between the high pressure unit cell data and the elastic tensor.

5 CONCLUSIONS

There has been a concerted effort in the last decade or so to provide elastic moduli for a range of materials at elevated pressure. Data at high pressure and temperature is still sparse, but as it grows, so too will the need for models that can reproduce this data in a reproducible format. The equation of state proposed in this paper has the benefit of being compact, self-consistent, and can be applied to pre-existing V(P,T) equations of state.

Not included in this paper is any treatment of isochemical variation. Isochemical variation can include order-disorder of chemical species on sites (e.g. Al and Si in plagioclase; Carpenter 1988), or structural flexibility (e.g. tetrahedral tilting in plagioclase; Mookherjee et al. 2016; Lacivita et al. 2020), or variation in proportions of spin states (e.g. iron in ferropericlase; Wu et al. 2013). Changes in isochemical state driven by changes in pressure or temperature can occur rapidly on the timescales of observations or natural phenomena such as seismic waves, and result in anomalous thermodynamic behaviour such as elastic softening. Plagioclase is one phase exhibiting such anomalous properties (Carpenter 1988; Mookherjee et al. 2016; Lacivita et al. 2020). A treatment of isochemical variations in anisotropic solid solutions will be the subject of a follow-up study.

ACKNOWLEDGMENTS

I would like to thank Bruce Hobbs and Nicolas Riel for their careful reviews of this paper, and for their encouragement while working on this and previous papers. This work was supported by NERC Large Grant MC-squared (Award No. NE/T012633/1) and STFC (Grant No. ST/R001332/1). Any mistakes or oversights are my own.

DATA AVAILABILITY

The anisotropic equation of state described in this paper is provided as a contribution to the BurnMan open source software project: https://github.com/geodynamics/burnman (Cottaar et al. 2014; Myhill et al. 2023).

REFERENCES