The hard sphere view of the outer core

Abstract

The hard sphere model for liquids attempts to capture the physical behavior of a real liquid in a simple conceptual model: a fluid of fixed size spheres that only interact repulsively when they come into contact. Is the model good enough to use for modeling internal planetary structure? To answer this question, I survey variants of hard sphere liquid theory by applying them to the Earth’s outer core to determine which of them explains wavespeeds in the outer core best. The variants explored here are the Carnahan-Starling hard sphere model, the Mansoori-Canfield extension to hard sphere mixtures, the transition metal hard sphere liquid, and the Lennard-Jones hard sphere liquid with attractive forces. With an empirical addition of a temperature dependence to the liquid’s hard sphere diameter, all of the variants explored can replicate wavespeeds in most of the radius range of the outer core. The hard sphere model for liquid transition metals explains the wavespeed best because it yields a mean liquid atomic weight of 48.8 g mol ⁻¹ at 10 wt% light element abundance in the core which is in good cosmochemical agreement with core light element models. Other variants also fit core wavespeeds but require implausibly low liquid mean atomic weight implying excessive incorporation of hydrogen or helium in the core. Applied to the detailed wavespeed structure of the Earth’s outermost outer core, the model suggests that the mean atomic weight could be reduced by up to 1.74% or the temperature could be increased by up to 400 K relative to an adiabatic profile, or there could be 8% fewer valence electrons in the liquid.

Background

According to PREM (Dziewonski and Anderson 1981), the outer core appears to be essentially in a state of adiabatic self-compression, a state representing a convectively well-mixed and chemically homogeneously material. Birch (1952) established the link between finite strain theories of material behavior and the variation of wavespeed along an adiabat by deriving the dependence of bulk modulus on pressure, finding it related to the gravitational acceleration (g) and the bulk sound speed derivative with depth (d V _Φ /d z).

In subsequent decades, physicists developed theories to describe fluids at high densities. They recognized that the dominant factors governing the behavior of high-density fluids were the repulsive interactions between the fluid constitutents (atoms or molecular species in the liquid) and the volume that those constitutents occupied in the liquid - the essential ingredients of the van der Waals equation of state for gases. A simple model captures these features: a fluid of spheres of fixed size that only interact repulsively when they come into contact (Carnahan and Starling 1969).

The Earth’s core is not far from this idealization. It is a liquid metal comprised of iron alloyed with perhaps 10% by weight of a combination of light elements (Birch 1952, 1964). Liquid metals are characterized by a bulk fluid of positively charged ionic cores separated from their valence electrons that serve to neutralize the cores’ charges (Hansen and McDonald 2013). Except for the valence electrons, this picture corresponds fairly well to a hard sphere fluid. For example, Stevenson (1980) used hard sphere theory to derive some fairly general properties of the core’s density and wavespeed dependence on depth, and Helffrich (2014) used it to estimate diffusion coefficients in the core.

One motivation for using the hard sphere model is to estimate seismic wavespeeds in planets whose interior structure is only roughly known. The basic information from planetary discovery missions such as MESSENGER (Smith et al. 2012) typically consists of a mean density and rotational moment of inertia. From this, a radial density model may be obtained (Smith et al. 2012). The radial density model then allows a hard sphere model to be applied to the liquid portions of the planet to predict seismic wavespeeds. In this way, preliminary seismic velocity models may be made to guide future planetary exploration missions that involve landers whose sensors include seismometers to determine the planet’s detailed internal structure (Dehant et al. 2012).

Another motivation for this work is to have a means to study the wavespeed profile in the Earth’s outer core directly. Because we know the aggregate elastic properties of the core quite well through decades of observational effort, it provides a good discriminator between successful and unsuccessful hard sphere theories. It is not obvious to the writer that a theory this simple should work at all. If one is suitable, it can serve as a reference for the core’s constitution from which changes from an adiabatic, well-mixed state may be assessed. Wavespeeds lower than PREM observed at the top and bottom of the outer core (Helffrich and Kaneshima 2010; Ohtaki et al. 2012; Zou et al. 2008) are not expected from traditional views of core behavior (where lower wavespeeds are associated with higher densities and vice versa), though here I focus only the top of the core in order to prove the viability of the approach.

Methods

Data

I use the PREM model (Dziewonski and Anderson 1981) for values of density (ρ), bulk sound speed (V _Φ ), and gravity (g) through the outer core. As a reference model for the outermost outer core, Tanaka (2007) found that the PREM travel time predictions for SmKS (the family of arrivals that travel through the outer core and reflect m−1 times from the underside of the core-mantle boundary (CMB)) are the best among more recent radial velocity models, thus recommending its use. Gravity is not directly parameterized by PREM but may be fit with a polynomial form to radius when integrated from the density profile. The explicit expression for gravity (m s ⁻²) in the outer core is

$$ \begin{aligned} {}g(r) =& \,2.114073\times 10^{-3} +3.818921\times 10^{-3} r-2.131285\\ &\times 10^{-7}r^{2}~, \end{aligned} $$

((1))

where r is the radius in km. A similar fit to ρ(r) to derive the pressure leads to a polynomial expression for pressure (in GPa) of

$$ \begin{aligned} {} P(r)~ =~& 3.853676\times 10^{2} -3.083633\times 10^{-2} r-1.184717\\ &\times 10^{-5}r^{2}~, \end{aligned} $$

((2))

again with radius r in km.

Wavespeeds in the outer core are closely approximated by adiabatic self-compression (Dziewonski and Anderson 1981). While PREM does not specify a core temperature, with a CMB reference temperature (T ₀) it can be projected down the adiabat using PREM data. Using the relation $T(h)=T_{0}\exp ({\gamma \int _{0}^{h}} g(z)/\mathrm {V_{\Phi }}^{2} (z)dz)$ (Stacey and Davis 2008), the pressure, density, and temperature at any depth h below the CMB may be calculated. Alfè et al. (2002) showed that the Grüneisen parameter γ is virtually constant throughout the core and has the value 1.52 that I use to calculate normal core temperatures. However, various shock wave studies and empirical models suggest a range of 1.3≤γ≤1.8 (Stacey and Davis 2008) so this is also investigated. I use a CMB temperature of 4,300 K throughout, except where the explicit thermal wavespeed variation is explored.

Methods used

There are various flavors of the hard sphere model proposed by different researchers. I focus on two properties of these models to apply them to the core. Each theoretical variant has an expression for the compressibility factor Z=P/(ρ R T), essentially providing the equation of state for the liquid. P (pressure), T (absolute temperature), and R (the universal gas constant) have their customary meanings, but the density ρ is in moles per unit volume. With a CMB temperature, all of these can be calculated from PREM. Z is related through the theory to the hard sphere diameter σ and thence to the hard sphere packing fraction η=π n σ ³/6. n is the number density, the number of particles per unit volume, related through the molecular weight of the liquid M to the PREM density.

The second property is the liquid wavespeed. This is fixed by the hard sphere packing fraction η and other free parameters of the particular theory variant. For each variant I search over the free parameters to determine the best-fitting set that reproduces PREM wavespeeds throughout the core. There are usually one or two free parameters to be searched over, so the best combination is determined computationally by a relatively undemanding grid search or line minimization. The following sections introduce the particular hard sphere theory variants.

Carnahan-Starling hard sphere liquid

The hard sphere model relates the compressibility factor Z to the hard sphere packing fraction η through a very simple formula. The best formulation of the hard sphere equation of state is Carnahan and Starling’s (1969) one, wherein

$$ Z~ =~ (1+\eta +\eta^{2} -\eta^{3})/(1-\eta)^{3}~. $$

((3))

Thus if Z is known, η may be solved for (provided Z≥1; clearly, 0≤η≤1) and σ obtained.

The hard sphere packing fraction also governs the wavespeeds in the hard sphere liquid. Rosenfeld (1999) derived the sound speed c in the hard sphere liquid in terms of Z (given by (3)), η and the molecular weight M:

$$ c^{2}~ =~\frac{RT}{M}\left[Z+\eta{\frac{dZ}{d\eta}} +\frac{2}{3}Z^{2}\right]~~~. $$

((4))

The only free parameter in this theory is the mean atomic weight, M.

The expression appears to indicate that c monotonically increases with T (and decreases with M), but this is not true because Z depends on T as well. Plotting c(T)/c(T _m ) (T _m is the melting temperature, defined for a hard sphere liquid when η=0.494) in Fig. 1, one sees the surprising result that c drops with decreasing η (increasing T) but then rises again as η decreases (T continues to increase). Thus the dependence on T is not straightforward and captures the behavior of both a liquid (wavespeed decreases with T at high η) and a gas (wavespeed increases with T at low η). For later use, it is helpful to define the quantity p(η)=(1+η+η ²−η ³)/(1−η)³ and its derivative with respect to η, p ^′(η)=2(2+2η−η ²)/(1−η)⁴.

Fig. 1

Hard sphere liquid properties. Variation of wavespeed in hard sphere liquid c with hard sphere packing fraction η. The speed is relative to the value at the melting temperature (T _m ), where there is a phase change in the liquid to solid form (at η=0.494). While Equation (4) suggests that the dependence should be monotonic on temperature (lower η implies lower density thus higher T), the co-variation of Z on T leads to the more complex trend depicted. After Rosenfeld (1999).

Hard sphere liquid mixtures

The core is impure and can be viewed as a liquid mixture of different alloying elements whose sizes and masses differ from iron. Mansoori et al. (1971) extended the hard sphere liquid model to include mixtures by adding further terms to the expression for Z:

$$ Z~ =~ \left[\left(1+\eta +\eta^{2} \right)-3\eta (y_{1} +y_{2}\eta)-y_{3}\eta^{3} \right]/(1-\eta)^{3}~. $$

((5))

The y _1,2,3 terms depend on the mixture properties. If x _i is the mole fraction of component i in the m-component mixture with $\sum _{i=1}^{m}x_{i} =1$, then

$$\begin{array}{*{20}l} \eta_{i} &=\frac{\pi}{6} n{\sigma^{3}_{i}}x_{i}~~~,~~~\eta~ =~\sum_{i=1}^{m}\eta_{i}~~~, \end{array} $$

((6a))

$$\begin{array}{*{20}l} y_{1} &=\sum_{j>i=1}^{m}\Delta_{ij} (\sigma_{i} +\sigma_{j})(\sigma_{i}\sigma_{j})^{-1/2}~~~, \end{array} $$

((6b))

$$\begin{array}{*{20}l} y_{2} &=\sum_{j>i=1}^{m}\Delta_{ij}\sum_{k=1}^{m}\left[\frac{\eta_{k}}{\eta}\right]\frac{(\sigma_{i}\sigma_{j})^{1/2}}{\sigma_{k}}~~~, \end{array} $$

((6c))

$$\begin{array}{*{20}l} y_{3} &=\left[\sum_{i=1}^{m}\left[\frac{\eta_{i}}{\eta}\right]^{2/3}x^{1/3}_{i}\right]^{3}~~~, \end{array} $$

((6d))

$$\begin{array}{*{20}l} \Delta_{ij} &= \left[(\eta_{i}\eta_{j})^{1/2} /\eta \right]\left[(\sigma_{i} -\sigma_{j})^{2} /(\sigma_{i}\sigma_{j})\right](x_{i}x_{j})^{1/2}~~~. \end{array} $$

((6e))

For a single component liquid, Δ _ii =0 implying that y ₁=y ₂=0, and y ₃=1, reducing Equation (5) to (3). In a two-component liquid, the theory’s free parameters are therefore the size ratio (relative to iron) and mass of the particle. Assuming a weight percentage of a light element in the core, these factors determine the mole fractions x _i and the mean atomic weight M in the liquid.

Transition metal hard sphere liquids

Yokoyama (2001a) showed that the hard sphere model could be successfully applied to calculate the sound velocities in liquid metals at 1 atm close to their melting temperatures. This was achieved by modifying Rosenfeld’s (1999) expression for the wavespeed dependence on packing fraction to account explicitly for the dependence of hard sphere diameter on temperature. In addition, Yokoyama (2001a) recognized that the electronic structure of the transition metals added a suite of longer-ranged but weak repulsive forces to the inter-particle interactions in the liquid. These two factors lead to modifications to the expression for wavespeeds. Including the hard sphere size dependence on T changes Equation (4) to

$$ \begin{aligned} {}c^{2}\! &=\!\frac{RT}{M}\left[{\vphantom{\frac{1}{1}}}p(\eta)+\eta p'(\eta)+\frac{2}{3}\left[3\eta p'(\eta)(\partial\log\sigma /\partial\log T)_{V}\right.\right.\\ &\qquad\qquad \left.\left.+p(\eta)\right]^{2}\vphantom{\left[3\eta p'(\eta)(\partial\log\sigma /\partial\log T)_{V}\right.}\right]~~~. \end{aligned} $$

((7))

If (∂ logσ/∂ logT) _V =0, Equation (4) is recovered. The explicit temperature dependence of σ given by Yokoyama (2001a) is $(\partial \log \sigma /\partial \log T)_{V} =-k_{1}\sqrt {T/T_{m}} /(1-k_{2} \sqrt {T/T_{m}})$ with constants k ₁=0.0444 and k ₂=0.112. However, for the hard sphere liquid, all properties should be relatable to η, and the explicit T dependence is only valid at a particular pressure, in this case 1 atm. Noting that (∂ logσ/∂ logT) _V =(T/σ)(d σ/d η) _V (d η/d Z) _V (d Z/d T) _V , an expression only involving η may be obtained: (∂ logσ/∂ logT) _V =−p(η)/(3η p ^′(η)). Substituting this into (7) eliminates the squared term and lowers wavespeeds. Rosenfeld (1999) shows this term is (∂ P/∂ T) _V . Using the Maxwell relation for V(P,T), (∂ P/∂ T) _V =α K _T which, if zero, implies that the thermal expansivity α is zero - unlikely for the core or indeed any substance in a planetary interior (Helffrich and Connolly 2009). Thus the η-only dependence through Z is unphysical and implies that a viable wavespeed model must contain factors leading to a nonzero thermal expansivity.

The electronic contributions to the wavespeed are an electron-gas energy (arising from the charge-compensating valence electrons in the liquid), the electrostatic repulsive energy of the charged ionic cores, and a s-d orbital hybridization term (B _H ) due to the irregular filling of the d orbitals among the transition metals. See Shimoji (1977) and Yokoyama (2001b) for a discussion of these terms. B _H is essentially a static correction to Z for a liquid for a particular metallic species (for the core, iron). Thus

$$ {\small{\begin{aligned} {}B_{H} =&\frac{a^{6}_{m}}{6}\epsilon_{\text{Ryd}} \left[0.031z+(0.916z^{4/3} \,+\,1.8z^{2})/a_{m} \,-\,4.42z^{5/3} /{a^{2}_{m}} \right]\\ ~& \quad-3k_{B}T_{m} \left(1+\eta_{m} +{\eta^{2}_{m}} -{\eta^{3}_{m}} \right)/(1-\eta_{m})^{3}. \end{aligned}}} $$

((8))

k _B is Boltzmann’s constant and ε _Ryd is the Rydberg energy, 2.1798741×10⁻¹⁸ J. a _m is the hard sphere radius at melting conditions nondimensionalized by dividing it by the Bohr length (5.2917725×10⁻¹¹ m), and η _m is the hard sphere packing fraction at melting (here taken to be 0.463 in contrast to Rosenfeld’s value). The particular metal is characterized by a valence electron count z, melting temperature T _m and melting density ρ _m (and number density n _m ) defining the hard sphere radius through $\eta _{m} =4/3\pi {a^{3}_{m}}n_{m}$. For iron, the valence in the liquid metal state is z=1.33 (Yokoyama 2001a).

Incorporating the electron gas (u _eg) and the ion repulsion (u _ion) contributions to the compressibility factor allows the hard sphere packing fraction η and Wigner-Seitz radius a to be determined. u _eg, u _ion, and B _H all serve to increase the effective Z, given by

$${} {\fontsize{7.8pt}{9.6pt}\selectfont{\begin{aligned} {}Z=&-\left[{\vphantom{\frac{1}{2}}}\epsilon_{\text{Ryd}} \left[0.031z+\left(0.916z^{4/3} +1.8z^{2}\right)/a-4.42z^{5/3} /a^{2}\right]-6B_{H} /a^{6}\right] \\ ~ &\,\,/(3k_{B} T) ~+~ \left(1+\eta +\eta^{2} -\eta^{3} \right)/(1-\eta)^{3}~~. \end{aligned}}} $$

((9))

With η (and a, again in multiples of the Bohr radius) obtained, the liquid metal wavespeed, incorporating (7), is

$$ {\small{\begin{aligned} {}c^{2}=&\frac{RT}{M}\left[{\vphantom{\frac{0}{0}}}\left.(\epsilon_{\text{Ryd}} [u_{\text{eg}} +u_{\text{ion}} ]+B_{H} /a^{6})/(k_{B} T)+p(\eta)+\eta p'(\eta)\right.\right. \\ ~&~\left.+\frac{2}{3}\left [3\eta p'(\eta)(\partial\log\sigma /\partial\log T)_{V} +p(\eta)\right ]^{2}\right]~~~, \end{aligned}}} $$

((10))

with u _eg=−0.031z/3−(4/9)(0.916z ^4/3+1.8z ²)/a and u _ion=(22.1/9)z ^5/3/a ². The electronic terms serve to raise wavespeeds relative to liquid without charge interactions both by increasing η for a given Z and by adding charge-related terms (u _eg, u _ion, and B _H ) to (7). Within this theoretical framework for a particular metal, the free parameters governing liquid speed are valence electron count z and mean atomic weight M.

The Lennard-Jones hard sphere liquid with attractive forces

A characteristic property of any liquid is its tendency to cohere and to maintain a fixed volume in the absence of confinement. A simple way to parameterize this is by defining an intermolecular potential based on separation r, u(r), whose gradient entails both repulsive and attractive forces. The Lennard-Jones potential is a common one used for this parameterization, given by

$$ u(r)~ =~ 4\epsilon \left[(\sigma /r)^{12} -(\sigma /r)^{6} \right] $$

((11))

(Fig. 2). ε is an interaction energy and σ is a length scale called the collision diameter. There are various ways to incorporate an interaction potential into a liquid theory, but for simplicity I choose Kolafa and Nezbeda’s (1994) parameterization based on the Barker and Henderson (1976) recipe for determining hard sphere size in the liquid. Again, the hard sphere diameter is chosen by matching Z, given in this theory by

$$ {}Z~ =~ p(\eta)+\rho^{*} (1-2\gamma\rho^{*})\exp (-\gamma{\rho^{*}}^{2})\Delta B_{2} (T)+f(T,\rho^{*})~~. $$

((12))

Fig. 2

Lennard-Jones potential. Lennard-Jones potential with ε=5. The potential scale governs the depth of the well at the minimum r/σ=2^1/6, beyond which the force between liquid constitutents is attractive. The function tends to zero at infinite separation and to infinitely strong repulsion at zero separation.

ρ ^∗=n σ ³ is the reduced number density, and η=π ρ ^∗/6. Kolafa and Nezbeda (1994) give polynomial expansions in powers of T ^1/2 and logT for Δ B ₂(T) and additionally powers of ρ ^∗ for f(T,ρ ^∗). γ=1.9290728 is an empirical fitting constant. The liquid wavespeed is given by Equation (7), with d Z/d η (=p ^′(η)) and (∂ logσ/∂ logT) _V evaluated numerically. Thus the free parameters in the method are the Lennard-Jones energy scale ε and M, the mean atomic weight of the liquid.

Results

Carnahan-Starling hard sphere liquid

The only free parameter in this theory is the molecular weight of the liquid. The seismic wavespeed is calculated with Equation (4). Minimizing the misfit to PREM wavespeeds by varying M leads to a best-fit M=6.1 g mol ⁻¹. The comparison with PREM (Fig. 3) shows the fit not only to be poor but leads to an unreasonably low molecular weight for the core. The model is evidently too simple to be a realistic description for the Earth. Equation (7) is an alternative expression for the hard sphere liquid wavespeed when there is an explicit temperature dependence to σ. If in the spirit of the hard sphere model this dependence is parameterized in terms of η, then a satisfactory fit may be obtained. However, parameterizations in η lead to d c/d T>0 which is unphysical for the larger η values found in the core (Fig. 1). To avoid this I parameterize the correction in powers of $\sqrt {T}$ as did Protopapas et al. (1973). The explicit form for the correction is in rational polynomial form, parameterized with numerator and denominator coefficients c _i and d _i :

$$ {}(\partial\log\sigma /\partial\log T)_{V}~ =~ (c_{0} +c_{1}\sqrt{T} +c_{2} T)/(1+d_{1}\sqrt{T} +d_{2} T) $$

((13))

Fig. 3

Hard sphere liquid performance. Optimized hard sphere liquid wavespeeds for the Carnahan-Starling model compared with PREM. The best-fit to PREM, using Equation (4) for the wavespeed, is for a liquid atomic weight of M=6.1g mol ⁻¹. If (∂ logσ/∂ logT) _V is used with Equation (7) for the wavespeed, a much better fit results virtually overlying PREM. Equation (13) parameterizes (∂ logσ/∂ logT) _V with (c ₀,c ₁,c ₂)=(−1.170076×10⁻¹,3.297547×10⁻³,2.274412×10⁻⁵) and (d ₁,d ₂)=(−0.0294468508,0.0002102121). The lines paralleling PREM show the wavespeeds for ±200 K change in T _cmb. Dotted lines show result of varying γ through the range shown. The result is similar to varying T _cmb and a corresponding change to (∂ logσ/∂ logT) _V .

(Figure 3). This correction is entirely an expedient tactic whose coefficient estimates come from fitting velocity residuals when M=15 g mol ⁻¹ to (13). It shows how good a hard sphere model can be if it can properly predict this dependence, however.

The applicability of the Carnahan-Starling (CS) liquid model to the core is severely limited. Fig. 4 depicts those limits, which derive from the liquid’s molecular weight. Only low molecular weight liquids (M<20g mol ⁻¹) are admissible throughout the core.

Fig. 4

Limits to hard sphere model. Limits to Carnahan-Starling (CS) hard sphere model. This diagram shows the upper limit (solid lines) for η in the core for the liquid wavespeed to be described by a real-valued (∂ logσ/∂ logT) _V labeled at right with the mean atomic weight for the liquid (M). The dashed lines show values for η for the CS hard sphere liquid obtained for compressibility factor Z calculated for PREM with CMB temperature T _cmb = 4,300 K with the labeled M at the left. At the same M, a crossing of the lines indicates an infeasible combination for the core. For any M≥20g mol ⁻¹, no CS liquid is feasible; if T _cmb = 5,300 K, the feasible M rises slightly to about 25 g mol ⁻¹. Dashed line at η=0.494 indicates hard sphere freezing and at η=0.704 indicates the hard sphere packing fraction for an hexagonally closest packed solid.

Hard sphere liquid mixtures

The free parameters in this model are the mass ratio and the radius ratio with respect to iron (at a fixed 10 wt% light element). A light element mass m at a fixed weight percent in the core fixes the mole fraction and thus M for the bulk liquid. A grid search over 1≤m≤56 and ratio 0.2≤r≤1 (corresponding to a range of element radii from He through H to Fe (Cordero et al. 2008)) yields a best-fit to PREM with r = 0.25 and m = 1.45 g mol ⁻¹. Fig. 5 shows the resulting fit, which is poor at the base of the core when η approaches freezing conditions for the hard sphere liquid. This model suggests that the core’s light element is dominantly hydrogen, surprisingly affirmed by the radius ratio for H/Fe: 0.31/1.32 = 0.235 (Cordero et al. 2008).

Fig. 5

Hard sphere mixture liquid performance. Best-fitting hard sphere mixture at 10 wt% light element compared to PREM. The second hard sphere component has radius ratio r=0.25 and molecular weight m=1.45g mol ⁻¹. The standard theory (with (∂ logσ/∂ logT) _V calculated numerically - ‘theor. mix’) fits PREM poorly. An empirical correction (‘corr. mix’) added to (∂ logσ/∂ logT) _V with (c ₀,c ₁,c ₂) = (−1.059965, 4.244515, −4.246083) and (d ₁,d ₂) = (−8.158507, 12.285397) improves the fit in the shallow core but diverges in the deeper core at packing fractions approaching hard sphere freezing. Dotted lines show result of varying γ through the range shown. The result is similar to varying T _cmb and a corresponding change to (∂ logσ/∂ logT) _V .

Transition metal hard sphere liquids

The free parameters in this model are the number of valence electrons z and the liquid molar mass M. The grid search is over 1≤z≤2.5 and 1≤M≤56. The best-fit combination is z=1.649 and M = 48.8 g mol ⁻¹. z is slightly higher than the 1 atm value at melting (z=1.33). Fig. 6 shows that the model liquid properties closely follow PREM through straight application of the theory. The average deviation from PREM throughout the core is 0.27%; this is only slightly larger than the approximately 0.25% wavespeed deviations found at the top of the outer core (Helffrich and Kaneshima 2010; Kaneshima and Helffrich 2013). It is possible to improve the fit by allowing a radial dependence for z and M, but I defer this to the discussion of this particular model’s adjustment to decrease speeds at the top of the core. This theory prescribes only a weak temperature dependence for the wavespeeds.

Fig. 6

Hard sphere metallic liquid performance. Plot of best-fitting Yokoyama hard sphere metallic liquid compared to PREM. The center line is for the T _cmb reported and ±200 K perturbations. For this model, the temperature dependence of the hard sphere liquid wavespeeds is low. The variation with γ is within the ±200 K bounds and is not shown for clarity.

The Lennard-Jones hard sphere liquid with attractive forces

The free parameters in this model are the energy scale ε and the liquid molar mass M. A grid search over the typical 1≤M≤56 and 1 K≤ε/k _B ≤5,000 K allows for two solutions. The first is a low- ε solution with ε/k _B = 20 K and M = 6.1 g mol ⁻¹. ε is so small that this is effectively no attraction at all at core temperatures, and is similar to the unmodified CS hard sphere model. The other solution has strong attractive forces with ε/k _B = 2,460 K and M = 6.5 g mol ⁻¹. The comparison to PREM is shown in Fig. 7. The mean deviation from PREM is 0.84%, but the radial trend is much flatter than PREM. This is a poor fit with the unmodified theory, but an empirical temperature-dependent correction can improve it in the shallower parts of the outer core.

Fig. 7

Lennard-Jones model performance. Plot of best-fitting Lennard-Jones hard sphere liquid compared to PREM. The dotted line (‘theory’) is the fit to the unmodified theory with ε/k _B = 2640 K and M = 6.5 g mol ⁻¹. The solid lines show an augmented model with a temperature dependent correction to (∂ logσ/∂ logT) _V with (c ₀,c ₁,c ₂) = (−7.928909×10⁻²,1.887171×10⁻³,−1.108079×10⁻⁵) and (d ₁,d ₂)=(−3.06458231×10⁻²,2.259932×10⁻⁴) in addition to the numerical differentiation of σ. The center line is for the reported T _cmb and ±200 K perturbations. PREM is effectively overlaid by the center solid line. Dotted lines show result of varying γ through the range shown. The result is similar to varying T _cmb and a corresponding change to (∂ logσ/∂ logT) _V .

Discussion

With approximate fits of each model to the core’s structure, we can ask what perturbations to it will approximate the observed wavespeed reductions from PREM. Fig. 8 shows the perturbations to individual models. The models, for comparison, are KHOCQ (Helffrich and Kaneshima 2010) and KHOMC (Kaneshima and Helffrich 2013). They deviate from PREM in the topmost approximately 350 km of the outer core by approximately 0.25%.

Fig. 8

Model comparisons to observed wavespeeds at top of outer core. Perturbations to variants of the hard sphere model compared to observed perturbations to PREM (KHOCQ and KHOMC) at the top of the outer core. (a) Hard sphere model; (b) hard sphere mixture; (c) transition metal hard sphere liquid; (d) Lennard-Jones hard sphere liquid attractive forces. Seismic models are KHOCQ and KHOMC (Helffrich and Kaneshima 2010; Kaneshima and Helffrich 2013) (solid lines). Dashed lines show various perturbations to the models that approximate seismic models: temperature deviations from adiabatic (dotted lines in (a-d)), mean atomic weight (dashed lines in (a-d)), radius ratio (dashed lines in (b)), valence electron count (dashed lines in (c)), ε/k (dashed lines in (d)).

The hard sphere model

The simplest model to evaluate is the hard sphere model (Fig. 1a). The only adjustable parameter in the theory is the mass M of the particle. Additionally, the temperature may be varied to allow conditions to deviate from adiabatic. The variation of each of these is shown in Fig. 8 by two virtually coincident lines deviating from the fitted model (see Fig. 3) in the topmost 375 km of the core. One case shows an increasingly subadiabatic (hotter) temperatures reaching about 17 K higher than the adiabat. This is a surprisingly small value but suggests that the heat flux would not be significantly reduced down the adiabat: by only approximately 4%. Alternatively, only a 2% decrease in density of the core material is needed to reduce the wavespeeds there. The mass and temperature perturbations are both approximately linear, so any linear combination of mass and temperature within these limits could also serve.

However, the reduction is from M = 15 g mol ⁻¹, an already extremely light core. For this M to represent a 10 wt% light element mix, the light element mass would be a combination of H+He and the concentration would be around 79 mol%. This is quite a large amount of hydrogen and helium given the expected core composition based on chondritic element abundances (Allègre et al. 1995), leaving no possibility of others whose presence is more likely on cosmochemical grounds.

The hard sphere mixture

The adjustable parameters in this model are radius and mass ratio relative to iron, and temperature. A hard sphere mixture is subject to an even smaller amount of subadiabatic temperature increase than a single-component hard sphere fluid, only about 10 K. The wavespeed decrease could also arise by a 3% increase in radius ratio relative to the reference ratio, 0.25, or an increase to 0.2575. Based on covalent radii for light elements, this represents a species change to Na or Mg among period 3 elements or Cr and Mn in period 4 elements: all others are significantly larger or smaller (Cordero et al. 2008). The mass must decrease, however, by about 1%, restricting the elements to Cr or Mn. The variations are again nearly linear, so linear combinations of any three of the factors could work to reduce wavespeeds in the core. However, the low m also suggests that the core is dominantly an H+He mix, unlikely for cosmochemical reasons as before.

The transition metal hard sphere liquid

This is the most intriguing of the models because it is a very good description of outer core wavespeeds all by itself (Fig. 6). The adjustable parameters are the M and z, along with T. The best-fitting values unfortunately yield the largest differences to PREM at both the top and base of the core, where the structural variations are most interesting. The basic model (Equation (10)) is therefore perturbed to minimize the difference to PREM by adding a radius dependence to the valence electron number δ z to fit PREM, given by

$$ \delta z(h)~ =~\left\{ \begin{array}{ll} &\frac{c_{0} +c_{1} h+c_{2}h^{2} +c_{3}h^{3}}{1+d_{1} h+d_{2}h^{2}}\quad,~ h>400 \\ &\qquad\frac{a_{0} +a_{1} h}{1+b_{1} h} \qquad\,\,\,\,,~ h{\leq}400\\ \end{array} \right. $$

((14))

with (c ₀,c ₁,c ₂,c ₃) = (−9.402105×10⁻², 7.577652×10⁻⁴, −1.62307×10⁻⁶, 1.039725×10⁻⁹), (d ₁,d ₂) = (−3.455461×10⁻³, 5.173742×10⁻⁶) and (a ₀,a ₁) = (−0.0958702639,0.0005057239) and b ₁ = 0.004435071 (h is in km in these expressions). This is admittedly overfitting PREM, but it leads to a simple, linear perturbation to the layer density in the outermost core shown in Fig. 9. The M decrease along with an adiabat is about 1.74%. At the same time, the valence electronic charge z decreases by 8% from 1.686 to 1.553. An unchanged z and M requires much higher temperatures to slow the layer due to the weak theoretical temperature dependence (Fig. 9): the topmost core is virtually 400 K hotter than the projection of a 4,300 K adiabat from 400 km deep in the core. This means that there is only a 5 K temperature drop throughout the topmost 400 km of the core. Hence there would be a low core heat flux, only 1.2% of the adiabatic one. This could contribute to solving the paradox of the young inner core age compared to the oldest recorded magnetic field (Gomi et al. 2013; Labrosse et al. 2001).

Fig. 9

Summary of perturbations to conditions at top of outer core. Perturbations to conditions at the top of the core to approximately fit the observed velocity perturbations to PREM for various hard sphere model variants: straight hard sphere (HS), hard sphere mixture (H S _mix ), a hard sphere liquid with electronic structure (H S _tm ), and a hard sphere liquid with attractive forces (H S _LJ ). Changes to temperature (Δ T), hard sphere diameter (Δ σ), density (Δ M), valence electron charge (Δ z), or interaction potential (Δ ε) are shown. Δ T is the temperature difference relative to the adiabatic temperature at the depth below the CMB. The other perturbations to σ, M, z, or ε are to the reference value yielding PREM (Figs. 3, 5, 6 and 7). Most hard sphere variants require only percent-level changes to intrinsic parameters or few-K changes to temperature to reduce wavespeeds at the top of the core, except for the transition metal hard sphere liquid.

The model, moreover, provides a sensible atomic mass for the core liquid, 48.8 g mol ⁻¹. If this represents a 10 wt% light element mix with iron, it maps into about a 20 mol% mixture of an element whose atomic mass is about 25 g mol ⁻¹. This puts it squarely in the range between C, N, O (on the light side - 12 to 16 g mol ⁻¹) and Si and S (on the heavy side - 28 to 32 g mol ⁻¹). Thus it is cosmochemically plausible for a quasi-chondritic Earth model. Allègre et al. (1995) suggest that there are 7.35 wt% Si, 2.3 wt% S, and 4.1 wt% O in the core, which, normalized to an average atomic weight, is 23 g mol ⁻¹, close to the value prescribed by this hard sphere model variant.

The Lennard-Jones hard sphere liquid

The adjustable parameters in this model are the energy scale ε and the atomic mass M, along with temperature. There could be a mass decrease of about 3% relative to the deeper core near the top. Alternatively, there could be a slight decrease, by 1%, to the energy scale. This would mean that the alloying element slightly reduces the intermolecular attraction between species in the liquid, but it is difficult to attribute it to any particular element’s intrinsic chemical properties such as its position in the periodic table. Temperature, on the other hand, is only 10 K higher than the adiabatic gradient, leading to slightly subadiabatic conditions and a 2.5% lower heat flux than what is conducted down the adiabat. On the other hand, the mean atomic mass is 6.5 g mol ⁻¹, implies that the core’s light element is pure hydrogen, with a mole fraction of 83%. This is again in conflict with cosmochemical expectations for the core’s composition because it excludes any other element from being in the core (Allègre et al. 1995).

Summary and analysis

Table 1 summarizes the performance of the models analyzed here. Though all hard sphere model variants may be made to fit the wavespeeds in the outer core, only one fits and provides a plausible view of the core from which sensible deviations may be explored: the transition metal hard sphere liquid.

Table 1 Summary of hard sphere model performance reproducing PREM

The main reason why most variants fail is because hard sphere models use as a reference state the ideal gas and add further Virial-like terms in powers of density that capture more complex interactions between liquid particles. Dense liquids like the core are quite far from the ideal gas reference and therefore require many terms before their behavior can be completely expressed. The Carnahan-Starling model actually captures Virial coeffients up to order 9 (Hansen and McDonald 2013), but this is clearly insufficient (Fig. 3). It serves as the basic reference for all of the hard sphere variants, so all are limited by the need for further higher order terms.

The transition metal liquid model recognizes additional repulsive forces that change the reference model to something other than an ideal gas. The particles have longer-ranged repulsive interactions that are essentially electrostatic. These act to increase wavespeeds relative to other models at a given Z (Equations (7), (9), and (10)). Hence this model leads to higher liquid mean atomic weights to counterbalance the additional force contributions and to match the wavespeeds in the core.

Conclusions

The hard sphere model for liquids may successfully be applied to the core. The standard versions - a hard sphere liquid with only repulsive forces, a liquid mixture of hard spheres, and a hard sphere liquid with attractive forces - require the addition of an explicit temperature dependence to the hard sphere diameter, but they can fit the uppermost outer core well. Among these, the mixed-sphere liquid performs worst near the base of the core where packing fractions are highest. All, however, suffer from the defect that the mean atomic weight of the core liquid is too low for any cosmochemically plausible composition because they require so much hydrogen and helium match core wavespeeds, leaving little room in the core for the heavier light elements that have sensible cosmochemical constraints (Allègre et al. 1995).

The best hard sphere model for the core is a transition metal liquid whose valence electrons form an electron sea bathing and charge-compensating the spherical ionic cores. This model yields a surprisingly good approximation to PREM wavespeeds throughout the core with a fixed number of valence electrons and atomic mass (z=1.649 and M = 48.8 g mol ⁻¹). To mix 10 wt% of a light element with iron to yield this atomic mass, that element’s mass must be about 25 g mol ⁻¹. This compares favorably to one estimate of the core’s average light element mass, 23 g mol ⁻¹ (Allègre et al. 1995). Thus one model satisfies PREM wavespeed constraints and cosmochemical constraints on the core.

Using the hard sphere electronic liquid model to examine the wavespeed deviations from PREM at the top of the outer core, it appears that a approximately 1.7% decrease in M at the top of the core (along with an 8% decrease in valence electrons z) can explain the wavespeed reduction there. With no change in M or z, a linear temperature increase over the adiabat of approximately 400 K can also explain the wavespeed reduction. Therefore, a thermally and chemically stratified layer at the top of the core is feasible within the framework of the hard sphere model within these mass anomaly and temperature limits. Its success further suggests that the model may, in future, play an integrated role in thermal histories of the core and mantle (Hirose et al. 2013; Labrosse et al. 2001).