Equations of State for the Deep Earth: Some Fundamental Considerations

: None of the 40 + equations that have been proposed to describe material properties at the pressures of the Earth’s core and mantle have escaped serious criticism. In this paper, some basic algebraic and thermodynamic constraints are reviewed, with the conclusion that the next step should be a re-examination of the relationship between the dependence of the bulk modulus, K , on pressure, P , that is K (cid:48) ≡ d K / dP , and the normalized (dimensionless) pressure, P / K . A linear relationship between 1/ K (cid:48) and P / K terminating at the inﬁnite pressure asymptote, at which these quantities become equal, has been used for analysing properties at extreme pressure, but may be inadequate for calculations requiring precise derivatives of an equation of state. A new analysis indicates that d ( 1/ K (cid:48) ) /d ( P / K ) increases with compression (or P / K ), but there are, at present, no reliable equations representing this. Relationships between higher derivatives of K and the thermodynamic Grüneisen parameter o ﬀ er the prospect of a resolution of the problem and hence a new generation of fundamentally-based equations of state. Although an earlier conclusion that a completely general ‘universal’ equation is not possible, in principle, is conﬁrmed in this study, the fundamental relationships present strong constraints for the forms of other proposed equations.


Introduction
Of the physical properties of geological materials that are central to studies of the lower mantle and core, a basic one is the bulk modulus, K = ρdP/dρ = dP/dlnρ, representing the variation of pressure, P, with density, ρ. Seismological observations of K(P) provide data that are required for the equation of state, which is the starting point for any theoretical description. This paper considers the fundamental developments that have been introduced since the pioneering contributions by Birch [1] and Anderson [2]. The regular increase of K with P, dK/dP ≡ K , over any depth range of a uniform phase and composition, is a primary equation of state parameter. Here, the fundamental theoretical constraints on its behaviour are explored, in particular, what we can learn from the extrapolation of K' and higher derivatives, K ≡ d 2 K/dP 2 , K and to their infinite pressure asymptotes, K ∞ and so on. A basic constraint [3] is (K P/K) ∞ = (dlnK/dlnP) ∞ = 1.
This is an algebraic identity that all proposed equations must satisfy. Most of them do so incidentally and not by design, but there is an advantage in making Equation (1) the starting point for the development of an equation of state. The significance of this to an understanding of the deep regions of the Earth is that in the core and deeper parts of the mantle, K is nearer in value to K ∞ than it is to K 0 , the zero pressure value, which is commonly used as an equation of state parameter. This means that K ∞ imposes a useful constraint on any equation in which it appears. As previously pointed The use of this equation has been limited by the fact that f is a parameter that depends on statistical details of atomic thermal vibration and is assigned different values in alternative theories. K T is the isothermal value of K. For most of the present discussion, it is not distinguished from the adiabatic modulus obtained from seismology. The relative difference diminishes with compression at all pressures and is not noticeable when applied to Equation (2). C V is the constant volume-specific heat. This equation is the generalized free volume formula, so named because it was first derived from free volume theory with f = 2 by Vashchenko and Zubarev [7]. A comparison of the equation with values of γ for common materials under ambient conditions shows that f < 2, and there is no requirement that it should be constant, independent of compression, or have the same value for different materials. A derivation of the V-Z formula (Equation (2) with f = 2) [8] with graphs of its pressure dependence for various equations of state shows that the formula assumes atomic thermal vibration that is random and uncorrelated with neighboring atoms. This is not valid. Lacking a reliable theory for f, Stacey and Hodgkinson [6] adopted an empirical value, f = 1.67, selected to make values of γ for the lower mantle of Equation (2) coincide with the values obtained by an alternative formula, the 'acoustic gamma', which assumes that all elastic moduli are adequately represented by the seismic P-and S-wave speeds. That calculation indicated a negligible pressure dependence. Equation (2) has an advantage for calculations of the kind considered here because, if f can have any value and any compression dependence, then the equation is completely general and, through f, incorporates the effect of the rigidity modulus. This generality extends to the infinite pressure extrapolation, in which f is canceled out by Equations (1) and (2), reducing it to Although f remains finite in this situation, it does not appear in Equation (3), which is, therefore, free of the uncertainty in f, but retains the generality of Equation (2) and is an identity that can be used to bring γ into the quest for fundamental equation of state constraints.
The two relevant derivatives of γ are In the infinite pressure extrapolation, q vanishes, as explained in Section 5 of [6], but λ ∞ is finite positive and is a useful parameter. Using Equation (2) with the assumption that f is constant, Shanker et al. [9] derived the equation Later derivations [6,10,11] removed the restriction on f and, like Equation (3), Equation (6) is an identity. It can be rewritten as the first of an infinite series of equations for dimensionless derivatives expressed in terms of K ∞ and λ ∞ , all of which are identities: and so on, with alternating signs. The quantities on the left-hand sides of these equations are all formally indeterminate, in the sense of being ratios of vanishing quantities. This makes it difficult to use the equations individually, but ratios of them, e.g., , can provide useful tests of validity for proposed equations, as in Section 5 of Stacey and Hodgkinson [6]. The fact that the infinite pressure asymptotes of derivatives of K can all be represented in terms of two constants invited an investigation of the possibility of writing a universal equation as a Taylor expansion of K as a function of P/K, anchored at the infinite pressure limit. All the coefficients in the series can be written in terms of K ∞ and λ ∞ , but are found to be indeterminate for all terms except the first, with denominators involving (1 − K P/K) ∞ , which becomes zero. In view of this result, Stacey and Hodgkinson [6] inferred that no universal form of high-pressure equation of state is possible, contrary to some expectations [12,13], and a semi-independent check of this is presented in Section 2.
The difficulty with higher derivatives of K does not affect the first derivative of K with P/K. This is conveniently written as because 1/K and P/K become equal at P → ∞ , as in Equation (1). The gradient in Equation (10) must be positive because λ ∞ < K ∞ , a condition that is basic to a constraint on equations of state, discussed in Section 4, additional to the K ∞ > 5/3 constraint imposed by a thermodynamic argument [3,4]. It needs to be recognised that, although K ∞ and λ ∞ are referred to here as constants, they are not universal constants, but are properties unique to each material. Numerical values or expressions for them are listed in Table 1 of Stacey and Hodgkinson [6] for a selection of equations that survived the culls prompted by basic tests, such as K ∞ > 5/3. In juggling Equation (10), it is useful to note another identity:

A Further Check of the Indeterminacy of Derivative Properties
If we attempt to avoid the problem of indeterminacy arising from denominators involving (7) to (9), we can differentiate derivatives with respect to P, yielding the result An indeterminacy problem also arises in this case because all the logarithms become infinite (positive or negative), but it is nevertheless possible to manipulate the infinite pressure limits by multiplying or dividing by complete expressions and applying cancellations of common factors within any expression. Therefore, the ratio of the first two expressions is in unity, with the common term, lnP, having been canceled, leaving Since, in the infinite pressure limit, both (KK ) and K 2 K vanish while their ratio remains finite, by the Bernoulli-l'Hospital rule [14], This reduces Equation (13) to 1 = 1, a symptom of a circular argument. It means that the algebraic identities involved in the manipulations cause the indeterminacy of a Taylor expansion using these derivatives, reinforcing the conclusion that there can be no universal form for equations of state.

Inner Core Problems
Of the deep regions of the Earth, the one most in need of theoretical input to the understanding of physical properties is the inner core. This is because it is the least accessible observationally. Direct P-wave propagation through it is well-observed and this gives a measure of the combination of elastic moduli (K + 4µ/3), but requires additional information to separate K from the rigidity modulus, µ. In principle, this would be obtained from the speed of shear waves, but they are not easily observed in the inner core. Conversion from and to compressional waves in the fluid outer core is extremely weak. Models of the Earth, such as Preliminary Earth Model (PREM) [15], rely largely on free oscillation records, which do not effectively resolve the fine structure. The modest variations in K and µ over the limited pressure range in the inner core contribute to the difficulty in discerning their behavior, with the result that models of the inner core structure are seen to be inconsistent with a plausible equation of state.
In the PREM tabulation, the problem is made most immediately obvious by a sharp decrease in K at the inner core boundary, from an outer core value of 3.75 to 2.34 in the inner core. This could not be caused by the liquid-solid transition, even allowing for a modest compositional difference between the inner and outer cores. At core pressures, both the solid and liquid have close-packed atomic structures and the atoms interact in the same way in both so that, for the same composition, the density difference would be 2% and the differences in K and K would both be~1.6% ( [3], Section 13). The anomalously low K in PREM is compensated for by an anomalously high value of dµ/dP, to give an inner core gradient of the observed P-wave modulus, d(K + 4µ/3)/dP, that is not anomalous. This may be only circumstantial evidence of a flaw in the PREM inner core model, but it can be checked by considering the compensating changes in the K and µ profiles, without invoking a comparison with the outer core. For a fixed phase of uniform material, the ratio µ/K is described by the µ-K-P equation ( [3], Section 8): where B is a negative constant. In the PREM inner core tabulation, µ/K increases strongly with P/K, as in the very short line in Figure 18.7 of [16], giving B ≈ + 0.68, compared with the expected value,~−1.2, derived from second-order elasticity theory ( [16], Section 18.8). A positive value of B is incompatible with any plausible equation of state for a material in a single stable phase. The fact that iron and the core are normal, in the sense of following Equation (15) with negative B, is evidenced by the very low value of µ/K in the inner core,~0.12, compared with the value for iron alloys under ambient conditions,~0.48. There is no reason for supposing that the normal behavior prevails down to core pressures, but is reversed in the inner core. As Stacey and Davis [16] argued, the anomalous behavior of µ/K and the implausibly low value of K in the PREM inner core, are simultaneously avoided by a simple empirical adjustment, leaving d(K + 4µ/3)/dP unchanged (Appendix F, Table F2 in [16]), because they are two aspects of the same problem However, this adjustment can be questioned because it does not involve any re-examination of the primary (seismological) data. Several authors have addressed the problem by re-assessing the data, but it has not yet been satisfactorily resolved. Another earth model, ak135 [17], derived with a greater emphasis on body wave travel times, permits an independent check. Its inner core value of d(µ/K)/d(P/K) is 0.68, coinciding with the PREM value, and the value of K is at least as low as that of PREM, with the further complication that K is also low at the bottom of the outer core. Deuss et al. [18] concluded that the inner core gradient of shear wave speed could not be satisfactorily resolved and attributed the average value, 3.6 k·ms −1 , to the whole depth range. If this is assumed, then d(µ/K)/d(P/K) ≈−1.25, coinciding with the theoretical expectation. This may be taken as evidence that a constant shear wave speed is a better approximation than mode variations, but other factors must be considered, especially the implausible gradient of K. Another difficulty is the density. Masters and Gubbins [19] developed a powerful method of choosing free mode periods that selectively emphasized specific depths and concluded that it indicated an inner core boundary (ICB) density increment of about 820 kg·m −3 , 200 kg·m −3 more than in the standard models. Using normal mode data, de Wit et al. [20] also favored a high ICB density increment, with a range of uncertainty extending to the Masters and Gubbins value. Such a revision reduces, but does not eliminate, the K and µ/K problems, which still need to be resolved.
Heterogeneity and anisotropy of the inner core introduce both difficulties and uncertainties to equation of state studies, but also offer some insight on its origin and development. The crystal structure of the inner core alloy has generally been considered to be hexagonal close-packed (hcp), but doubts have arisen, with rival suggestions for alternatives and even multiple alternatives [21]. One suggestion was for a subtle modification of the hcp structure (dhcp). This would not cause a dramatic revision of core physics, but another approach, proposed by Belonoshko et al. [22] and subsequently strongly advocated by the same group, was for a return, at the inner core pressure and temperature, to a body-centered cubic (bcc) phase, familiar for iron under ambient conditions. This would necessitate some rethinking of core properties and Stacey and Davis [16] argued that a case against it was that the core equation of state did not satisfactorily extrapolate to the familiar properties of laboratory bcc iron. This bears directly on the anisotropy, which is not easily explained as a property of bcc iron, but has a more natural explanation with the anisotropic hcp structure. Yoshida et al. [23] argued that it is caused by deformation because the inner core grows by the freezing of added material on to an equatorial band and the consequent excess ellipticity relaxes by progressive deformation, causing crystalline alignment. With this interpretation, there is an interesting implication in the irregularity of the anisotropy: polar wander of the inner core. Following Curie's principle of symmetry, which disallows any effect with lower symmetry than the combination of its causes, we see that the symmetry constraints on inner core formation are the superposition of axial rotation on the sphericity. This is a close analogy to the control on the geomagnetic dynamo which, when averaged over the randomizing effect of secular variation, has axial symmetry. Similarly, the driver for equatorial growth of the inner core and consequent deformation has axial symmetry, but its rotational axis moves about in response to the varying electromagnetic torques of outer core dynamo action. Although the average alignment may be axial, the anisotropy of any part depends on its position and orientation at the time of formation or deformation.
The outer core presents a somewhat simpler problem. Internal fluid motion at speeds of tens of km/s ensures homogeneity and a temperature gradient that departs insignificantly from the adiabat. This gives some confidence to the application of equation of state principles, leading to the conclusion that K ∞ = 3.05 ± 0.10, which, by Equation (3), means γ ∞ = 1.36 ± 0.05. It is clear that geophysical data cannot come anywhere near to accommodating the estimate γ ∞ = 0.5, which has been favored by some shock wave theories (e.g., [24]).

Is There a Best Equation of State?
It may appear that this question is answered by the conclusion of Sections 1 and 2 that no universal equation of state is possible, in principle, but that sidesteps the question. We are still interested to know whether one form of equation is most appropriate for all regions of the Earth, and all circumstances, or whether different equations should be applied, for example, to the core and mantle, or to different phases of the same material. Even if a single equation is acceptable for several materials, the parameters and coefficients will differ, but that is another question. The requirement K 0 > K ∞ > 5/3 is a condition for the plausibility of an equation. However, it does not constrain K ∞ sufficiently to give reliable values of the derivatives of K that describe thermal properties, as well as compositional differences. Stacey [5] noted that, for both lower mantle and core equation of state fits, K ∞ = (3/5)K 0 within the uncertainties of the data. If this is used, it is a very strong constraint, but it is empirical, not based on fundamental principles, and may be specific to terrestrial materials or influenced by the equation of state used in the analysis. However, it complies with the general rule proposed by Keane [25]: The search for further tests and constraints continues. Stacey and Hodgkinson [6] showed that λ ∞ must be positive with the restriction where f ∞ is the infinite pressure limit of the parameter f in Equation (2). It probably differs very little from the lower mantle value of f, 1.67, but cannot plausibly be outside the range 1.5 < f ∞ < 1.8, and this imposes a bound on λ ∞ : By imposing this bound on Equation (10), we have a new equation of state constraint: The significance of this to the present discussion is that it is quite a tight constraint on the gradient of the approach of a 1/K vs. P/K equation to the infinite pressure asymptote and is more restrictive than the λ ∞ and K ∞ constraints. In particular, it calls into question the reciprocal K-primed (RKP) equation [4], which was designed to take advantage of the relationship between K and P/K by building on Equation (1): Table 1 gives values of the quantities that are compared in Equation (19) for a selection of equations. A tick in the final column indicates an equation that passes the test imposed by Equation (19), regardless of its performance in the λ ∞ and K ∞ tests. Crosses mark equations that clearly fail this test and those that are marginal or doubtful are assigned a~symbol. The consistent ticks for the Birch equations in Table 1 make them appear as favorites for selection as the preferred choice, but with some qualifying comments, including the observation that for the highest pressures, the reciprocal K-primed equation (Equation (20)) is more satisfactory, despite the new doubt about it. The Birch equations are generally presented as successive approximations to an infinite series, but should not be considered that way, as they behave quite differently at high pressure. It is misleading to regard them as successive improvements of the same basic equation. The value of K ∞ increases by 1/3 for each increase in order, 7/3 for the second order, 3 for the third order, and so on, becoming infinite for the complete series. As a result, the third-order Birch equation performs best for the core and the second order for the lower mantle, and the fourth and higher orders cannot be favored at all. Since the second-order equation is a special case of the Mie equation, the general Mie equation comes into contention as a rival. This conclusion is proposing nothing more than a stop gap measure as these equations have frequently been noted to have serious flaws [12,[31][32][33], making them unusable outside the ranges of fitted data (if at all). Both the third-and fourth-order Birch equations, fitted to terrestrial data, lead to negative values of P and K not far from the terrestrial pressure range [16], making their derivatives unsatisfactory for thermal calculations even within that range. There appears to be no short cut to a choice of preferred equation, or modification of an existing one, but Section 5 presents a suggestion for a possible new approach.

Concluding Comments
We may be glimpsing a step forward on the erratic path to an analytical form of a high-pressure equation of state, one that will describe thermal, as well as mechanical, properties. Although most of the equations that have been presented so far are supported by logical arguments, and even sophisticated theoretical support, they are empirical, in the sense of relying on assumptions and intuitions that make them vulnerable to questioning. This note pursues a recent development of the questioning process, using only algebraic and thermodynamic principles. One of the outcomes is the conclusion of Stacey and Hodgkinson [6], extended in Section 2, that these principles do not lead to a completely general or 'universal' form of equation. However, they present criteria that all equations need to satisfy, as summarized by Equations (1) and (3), the thermodynamic bound on K ∞ and the related relationship for λ ∞ . Although these principles have sufficed to discount the majority of the 40+ proposed equations of state, their effectiveness is restricted by the appeal to infinite pressure asymptotes. While these introduce some rigorous conclusions by simplifying basic algebraic and thermodynamic relationships, they offer only indirect constraints on properties in the accessible range. Equation (19) presents a further development that addresses this limitation by extending the fundamental approach to finite pressures.
Recognizing that all equations of state must have end points represented by Equation (1), Stacey [4] compared several simple relationships between K and P/K with terrestrial data and concluded that a reciprocal relationship was the only one that gave a reasonable fit. This was the basis of the RKP equation (Equation (20)), which became a serious candidate for the choice of 'best' equation. Equation (19) and Table 1 now show that it needs improvement and, most importantly, the nature of its shortcoming. Equation (19) gives the gradient of a graph of 1/K vs. P/K for the approach to infinite pressure, showing it to be steeper than the gradient presented by Equation (20), which assumes a uniform gradient from 1/K 0 at P = 0. This means that d(1/K )/d(P/K) increases with P/K and is not constant, as Equation (20) assumes. The graph should be a curve that is concave upwards. An independent indication of that is a graph of 1/K vs. P/K for MgO at 300 K, as calculated by Stacey and Hodgkinson [6] from a pressure-density tabulation in a report of a 'first principles' calculation by Wu et al. [34]. This graph is concave upwards, supporting the conclusion from Equation (19) that Equation (20) needs to be modified in this way.
There appears to be the possibility of a laboratory check of the curvature of a 1/K vs. P/K graph from measurements of derivatives of K at low to moderate pressures. For the RKP equation (Equation (20)), d(1/K )/d(P/K) = (−KK )/ K 2 (1 − K P/K) = 1 − K ∞ /K 0 at all pressures and at low pressure, is thus equal to −K 0 K 0 /K 0 2 . If the suggestion that K ∞ /K 0 ≈ 3/5 is accepted as general, then for a linear 1/K vs. P/K relationship (the RKP equation), A lower value of this quantity is required for a 1/K vs. P/K graph to be concave upwards, as expected by Equation (19), and a higher value would indicate a curve that is convex upwards. This looks like a useful test, but measurements of K 0 are difficult and available values appear too imprecise to provide a convincing check. This problem needs a more thorough investigation before we can be clear that the RKP equation should be superseded and, until then, it must be considered the best available equation for very high pressures.
It remains to be seen how much further the fundamental approach can go. The Taylor expansion attempt failed, but that was conducted with the aim of finding a 'unique' equation. Perhaps the derivatives in Equations (7) to (9) or (13) could still yield a useful conclusion with a less ambitious target. Now that Equation (19) has qualitatively shown what to expect, further effort on the algebraic-thermodynamic approach will be better directed. Although the 1/K vs. P/K relationship will not necessarily give the needed inspiration, it will indicate whether an idea is on track and appears to be central to further progress. The aim is to establish algebraic relationships that give the reliable derivatives required for the calculation of thermal properties, such as γ, as well as minimizing the use of adjustable constants in data fitting and identifying conflicts between observations and theory.