that is applied to the multi‐particle wave function Ψ(r). In density functional theory (DFT), the Hamiltonian operator is broken down into contributions of the kinetic energy T of non‐interacting electrons, the electrostatic potential energy EC, and the exchange‐correlation energy EXC that describes the electron‐electron interaction, all of which are functionals of the electron density n(r) that changes with the location r (Hohenberg and Kohn, 1964; Kohn and Sham, 1965). The total energy can then be expressed in terms of the electron density:
While the functionals of kinetic and potential energies can be evaluated exactly, formulations for the exchange‐correlation energy rely on approximations (Karki et al., 2001a; Perdew & Ruzsinszky, 2010; Stixrude et al., 1998). These approximations allow computing the total energy E for a given arrangement of atoms, i.e., for a crystal structure. The forces acting between atoms, also called Hellmann‐Feynman forces, can be found by evaluating the changes in energy that result from small perturbations of the atomic arrangement (Baroni et al., 2001). When the perturbation of the atomic arrangement is chosen to correspond to a homogeneous strain, the resulting stresses (Nielsen and Martin, 1985) and hence elastic properties can be derived (Baroni et al., 1987a, 2001; Wentzcovitch et al., 1995). Alternatively, an external pressure can be applied to the simulation cell, and the resulting forces and stresses be minimized by relaxing the atomic positions and the shape of the simulation cell (Wentzcovitch et al., 1995, 1993). Note that, apart from being intentionally displaced or being relaxed to their equilibrium configuration, atomic nuclei remain static in these calculations. The results of such static DFT calculations can in principle be compared with those of experiments at ambient temperature as thermal contributions at 298 K are expected to be fairly small for most materials.
In most cases, the exchange‐correlation contribution is evaluated based on either the local density approximation (LDA) or on generalized gradient approximations (GGA). The LDA approximates the exchange‐correlation potential at each point r by the exchange‐correlation potential of a uniform electron gas with the density n(r) at this point. GGA also take into account gradients in electron density when evaluating the exchange‐correlation potential. Both LDA and GGA have proven successful in reproducing the elastic properties of many mantle minerals where experimental results are available for comparison. However, LDA appears to systematically overestimate bonding forces, i.e., leading to slightly higher elastic moduli and smaller volumes, while GGA tends to underestimate bonding forces with opposite effects on elastic properties and volumes (Karki et al., 2001a). Another aspect becomes important when computing properties of compounds containing transition metal or lanthanide atoms such as iron and its cations that are present in many mantle minerals. Both LDA and GGA take into account mutual repulsion between electrons in an inhomogeneous but diffuse gas of electrons. This approximation breaks down for electrons in the less diffuse d and f orbitals of transition metal and lanthanide atoms as these electrons become effectively localized by the additional repulsion they experience from electrons in similar orbitals on the same and neighboring atoms (e.g., Cox, 1987; Hubbard, 1963). This localization is the reason for the insulating or semi‐conducting behavior of some transition metal oxides, such as FeO, despite of their partially filled d orbitals.
By far the most common strategy to account for the repulsion between localized d electrons is to add an energy term EU that depends on the Hubbard parameter U and on the occupation numbers of d orbitals (Anisimov et al., 1997, 1991; Cococcioni, 2010; Cococcioni and de Gironcoli, 2005). This approach, referred to as LDA+U or GGA+U, can reduce discrepancies between predicted and experimentally observed elastic properties (Stackhouse et al., 2010) and led to substantial improvements in treating spin transitions of ferrous and ferric iron with DFT calculations (Hsu et al., 2011, 2010a, 2010b; Persson et al., 2006; Tsuchiya et al., 2006). Despite important progress in modeling the electronic properties of iron cations in oxides and silicates, elastic properties extracted from DFT computations still seem to deviate significantly from experimental observations in particular across spin transitions in major mantle minerals (Fu et al., 2018; Shukla et al., 2016; Wu et al., 2013), highlighting persistent challenges in the treatment of localized d electrons.
While most first‐principle calculations assume the atomic nuclei to be static, thermal motions of atoms at finite temperatures can be addressed by coupling DFT to molecular dynamics (MD) (Car & Parrinello, 1985) or by computing vibrational properties using density functional perturbation theory (DFPT) (Baroni et al., 2001, 1987b; Giannozzi et al., 1991). Sometimes referred to as ab initio molecular dynamics, DFT‐MD allows computing elastic properties at pressures and temperatures that span those in Earth’s mantle (Oganov et al., 2001; Stackhouse et al., 2005b). Within the limitations imposed by the finite sizes of systems that can be simulated, DFT‐MD includes anharmonic effects that go beyond the approximation of atomic vibrations as harmonic oscillations and become discernible at high temperatures (Oganov et al., 2001; Oganov and Dorogokupets, 2004, 2003). Alternatively, the frequencies of lattice vibrations can be derived from DFPT for a given volume (Baroni et al., 2001, 1987b) and then used in the quasi‐harmonic approximation (QHA) to compute temperature‐dependent elastic properties (Karki et al., 1999, 2000; Wentzcovitch et al., 2004, 2006, 2010a). The pressure‐temperature space for which the QHA remains valid for a given material can be estimated from the inflexion points (∂2α/∂T2)P = 0 in computed curves of the thermal expansivity α(P, T) as the QHA appears to overestimate thermal expansivities at higher temperatures (Carrier et al., 2007; Karki et al., 2001b; Wentzcovitch et al., 2010b). This criterion suggests that the QHA should remain valid throughout most of Earth’s mantle for some materials while others are expected to deviate from purely harmonic behavior (Wentzcovitch et al., 2010b; Wu & Wentzcovitch, 2011). The results of QHA‐DFT computations can be corrected for anharmonic contributions by adding a semi‐empirical correction term to match experimental observations (Wu and Wentzcovitch, 2009). Anharmonic effects can also be addressed in DFPT computations (Baroni et al., 2001; Oganov & Dorogokupets, 2004). For example, a recent study on MgO combined DFPT calculations with infrared spectroscopy and IXS to relate experimentally observed indications of anharmonicity, such as phonon line widths, to multi‐phonon interactions (Giura et al., 2019). Such efforts demonstrate the possibility to assess anharmonicity in first‐principle calculations. Among other advancements, all these developments facilitate the routine application of DFT computations to study the elastic properties of structurally and chemically complex minerals at high pressures and high temperatures (Kawai and Tsuchiya, 2015; Shukla et al., 2015; Wu et al., 2013; Zhang et al., 2016).
3.5 PARAMETER UNCERTAINTIES
The results of a series of experiments or computations are typically inverted into sets of parameters that describe the variation of elastic properties as a function of finite strain and temperature. Different definitions of finite strain result in different functional forms of expressions for elastic properties. For equations of state, different definitions and assumptions give rise to a remarkable diversity in EOS formulations (Angel, 2000; Angel et al., 2014; Holzapfel, 2009; Stacey & Davis, 2004). The variation of components of the elastic stiffness tensor has been formulated using both the Lagrangian and Eulerian definitions of finite strain (Birch, 1947; Davies, 1974; Thomsen, 1972a). For interpolations between individual observations, i.e., experiments or computations, differences between formulations should be small since, in a semi‐empirical approach, finite‐strain parameters are chosen to best reproduce the observations. Extrapolations of elastic properties beyond observational constraints, however, can be extremely sensitive to the chosen formalism. Both different definitions of finite strain and expansions to different orders of finite strain can result in substantial deviations between different finite‐strain models when extrapolated beyond observational constraints. Owing to the prevalence of elasticity formalisms based on Eulerian finite strain (Davies & Dziewonski, 1975; Ita & Stixrude, 1992; Jackson, 1998; Sammis et al., 1970; Stixrude & Lithgow‐Bertelloni, 2005), this source of uncertainty is not commonly taken into account when computing seismic properties