Accurate thermoelastic tensor and acoustic velocities of NaCl

Most information about Earth’s interior is obtained by analyzing the internal propagation of seismic waves.^1–4 However, interpretation of seismic velocity models requires knowledge of equations of state (EoS) and thermoelastic properties of rock-forming minerals.^5,6 As such, accurate estimates of thermoelastic tensors of minerals are fundamental in geophysics. They are also important for the gas/oil industry. Several gas fields have been discovered with the use of seismic surveys.⁷ In the context of gas/oil exploration, rock-salt NaCl has been attracting great interest since the discovery of large reserves beneath thick layers of salt on continental shelves, the so called pre-salt fields. In this letter, we report accurate thermoelastic properties of rock-salt NaCl (B1 structure) at relevant pressures (up to 0.5 GPa) and temperatures (up to 400 K) and beyond, which should be useful for surveying pre-salt fields. First, we perform ab initio free energy calculations based on the quasi-harmonic approximation (QHA)^8,9 using state-of-the-art density-functional theory (DFT) based codes (Quantum ESPRESSO¹⁰ and VASP¹¹). Then, we use a recently proposed method to reconcile ab initio results and experimental data for thermal EoS.¹² Here, this method is further extended to reconcile ab initio results and experimental data for elastic coefficients. Experimental techniques to measure elastic properties at high pressures and temperatures (PT) are very well established. Ultrasonic interferometry^13–15 and Brillouin scattering^16–19 have provided fundamental information for interpretation of wave-velocities in Earth’s interior. These properties can be accurately determined at low pressures and temperatures. However, experimental data become scarcer at combined high pressures and temperatures, including for NaCl. In contrast, ab initio calculations of thermal elastic coefficients can reach much higher pressures and temperatures.^5,6,20–22 The limiting factors in these calculations are: a) systematic deviations originating in the choice of exchange-correlation functional in DFT and b) limitations associated with the QHA - absence of intrinsic anharmonic effects. It is therefore desirable to combine experimental data and ab initio results in a sensible manner to produce optimally predicted thermal elastic properties.

Ab initio calculations were perfomed within the QHA, with three different functionals: LDA,^23,24 GGA-PBE,²⁵ and HSE06.²⁶ We used the projector augmented wave (PAW) method²⁷ to describe the electronic wave function with a Brillouin zone sampling on a displaced 8x8x8 k-points mesh. Dynamical matrices at each sampled volume were calculated by the finite displacement method, using a 4x4x4 supercell for LDA, and GGA-PBE, and 2x2x2 for HSE06. Force constants were interpolated on a 12x12x12 mesh to produce phonon density of states (VDoS) to be used in conjunction with QHA free energies calculations:

where q is the phonon wave vector, m is the normal mode, and E is the static internal energy at volume V. The second and third terms of the r.h.s are the zero point and the thermal contributions to the free energy, respectively. Rock-salt type NaCl has cubic structure with three independent elastic coefficients. At a given volume, negative and positive strains of 1% were applied and static elastic coefficients were calculated using the relation:

Isothermal elastic coefficients can be separated into static and temperature dependent (vibrational) parts:

The vibrational contribution to $C i j k l v i b$ was calculated using the semi-analytical approach by Wu and Wentzcovitch.²² This method avoids the calculation of VDoS for strained configurations, by assuming an isotropic distribution of phonon frequencies. This method provided accurate high temperature/pressure elasticity of several minerals.^22,28–30 With the isothermal elastic coefficients, the adiabatic ones can be obtained as:

where C_V is the heat capacity at constant volume and S the entropy at P and T.²¹ It has been well documented^5,6 that the popular exchange correlation functionals LDA and PBE-GGA in general give static volumes smaller and larger than the experiments, respectively. These systematic deviations impact directly on the elastic properties. For instance, deviations in the bulk modulus correlate inversely with deviations in equilibrium volume.³¹ LDA with vibrational effects accounted by the QHA offers most accurate EoS for weakly correlated and strongly bonded materials. Several corrections have been designed to reconcile experiments and calculations and to accurately predict high temperature EoS.^12,31–33 Some of these corrections to the theoretical EoS have adjustable parameters that minimize the difference between corrected ab initio results and experimental data.^12,32 Other corrective schemes consistently rescale EoS variables to bring calculated and measured EoS parameters into close agreement.^31,33 Here we extend one of these corrective schemes¹² to the elastic coefficients. Differences between ab initio results and experimental data are pressure dependent, therefore, the 300 K corrected volume is written as:

where ΔV(P, 300) is a numerical or analytical function modeling the difference between ab initio results and experimental data. Assuming that the error in thermal expansion is small implies that ΔV(P, 300) ≅ ΔV(P, 0), then:

Therefore, corrected static volumes are obtained by subtracting the extra volume difference due to zero point motion:

which can be written as:

f(P) is then used to correct the static energy that enters in QHA calculations. This corrective procedure offers highly accurate hybrid EoS¹² in the range of validity of the QHA and here we extend it to elastic coefficients. These calculations start by computing static coefficients (first term on the r.h.s of Eq.(3)) using Eq. (2). Isothermal elastic coefficients can also be written as:

where G(P, T) is the Gibbs free energy:

Thus, the static elastic coefficients are:

Conversely, the static compliance tensor is:

Since static volumes are calculated via the Gibbs free energy:

the corrected volume must be the derivative of a “corrected” Gibbs free energy:

as in Eq. (8). Integrating Eq. (14) to obtain $G st cor$ and using Eq. (12), we arrive at the corrected compliance tensor:

from which the elastic coefficients can be calculated:

Similarly to corrected thermal EoS, as long as the thermal expansivity is well described by the uncorrected ab initio calculation and the temperature is within the range of validity of the QHA, there is no need to correct the thermal contribution to the elastic coefficient.

Calculated and measured^34,35 300 K EoS of NaCl are compared in Figure 1. LDA and GGA-PBE volumes are smaller and larger, respectively, than experimental measurements. HSE06 improves the EoS, but still produces a large deviation from experiments. Figure 2 shows the thermal expansion coefficients. GGA-PBE and HSE06 overestimate the thermal expansion, but LDA describes it well up to 600 K.^5,6,21 This is essentially the LDA limit of QHA for this material, where anharmonic effects start causing deviations of corrected results from experimental data.¹² Corrective procedures for the EoS depend on a good description of the thermal expansivity. Thus, here we perform corrections only on LDA results. Calculated adiabatic elastic coefficients are shown in Voigt notation (C₁₁ = C₁₁₁₁, C₁₂ = C₁₁₂₂ and C₄₄ = C₂₃₂₃) in Figure 3 compared to experimental data from Kinoshita et. al.,³⁶ Voronov et. al.,³⁷ Whitfield et. al.³⁸ and Yamamoto et. al.³⁹ As LDA underestimates the experimental volume, C₁₁ is overestimated, similarly to pressure.³¹ The opposite occurs with GGA-PBE and HSE06. HSE06 agrees best with experiments, but still has large errors. C₄₄ is less volume-sensitive and agrees well with experimental data.

We used the corrected EoS^12,40 with this corrective scheme to produce the elastic coefficients shown in Figure 4. The pressure dependence of C₁₁, C₁₂ and C₄₄ at 300 K is in excellent agreement with single crystal data by Kinoshita et. al.³⁶ and Voronov et. al.³⁷ Temperature dependent results at 0 GPa are also improved for C₁₁, but it appears that even at 300 K the data by Yamamoto et. al.³⁹ differs from those by Kinoshita et. al.³⁶ and Voronov et. al.³⁷ As shown in Figure 2, above about ≈500 K, anharmonic effects cause deviations from experiments and further anharmonic corrections^12,32 or fully anharmonic calculations^41,42 are needed. With these elastic coefficients, bulk and shear moduli were calculated using Voigt-Reuss-Hill (VRH) averages^43–45 and are shown in Figure 5(a). Excellent agreement with data by Kinoshita et. al.³⁶ and Voronov et. al.³⁷ is expected given the good agreement between elastic coefficients. The corrected bulk modulus is not significantly changed, compared to the uncorrected one. Experimental bulk modulus obtained from the EoS by Decker⁴⁶ and from the hugoniot by Fritz et. al.³⁵ deviate somewhat from those calculated using corrected elastic coefficients. Some differences are also observed with respect to those obtained from measured velocities in polycrystalline samples.⁴⁷ Similar agreement and deviations are observed between calculated and measured acoustic velocities.

In summary, the best agreement between corrected and measured elastic coefficients is obtained by comparing results with single crystal elasticity data, such as those by Kinoshita et. al.,³⁶ Voronov et. al.³⁷ and Whifield et. al.³⁸ The C₁₁ mean squared error is significantly reduced from 144 to 5, comparing to Voronov’s data.³⁷ Data obtained on polycrystaline NaCl⁴⁷ do not agree as well, nevertheless, they are within the Voigt-Reuss-Hill bounds. Such aggregates should be more isotropic than single crystals. The anisotropy of single crystal NaCl increases with pressure, since the deviation from the isotropic condition, C₁₁ − C₁₂ = 2C₄₄ increases with pressure. This increased anisotropy can also be seen in Fig. 5(a), where Reuss, G_R, and Voigt, G_V, bounds⁴³ are shown (K_R = K_V). For an isotropic system, G_R = G_V, a condition held only at 0 GPa. With increasing pressure, G_R and G_V depart from each other, a clear signature of the pressure enhanced anisotropy. The difference between corrected results and data on aggregates with increasing pressure, as shown in Fig. 5(b), is therefore likely to be caused by the single crystal anisotropy. Further work is needed to address this issue.

We have extended a scheme to combine/reconcile experimental data and ab initio results for thermal EoS¹² to thermoelastic coefficients. Using this approach, we were able to predict quite accurately the single crystal elastic tensor and, within uncertainties, acoustic velocities in poly-crystalline rock-salt NaCl at relevant thermodynamics conditions. The proposed corrective scheme is general, applies to other crystalline systems, and can extend data with predictive accuracy to thermodynamics conditions difficult to be reached by experiments, as long as anharmonicity does not compromise the accuracy of QHA-based predictions.