FirstPrinciples Lattice Dynamics Method for Strongly Anharmonic Crystals
Abstract
We review our recent development of a firstprinciples lattice dynamics method that can treat anharmonic effects nonperturbatively. The method is based on the selfconsistent phonon theory and temperaturedependent phonon frequencies can be calculated efficiently by incorporating recent numerical techniques to estimate anharmonic force constants. The validity of our approach is demonstrated through applications to cubic strontium titanate, where overall good agreements with experimental data are obtained for phonon frequencies and lattice thermal conductivity. We also show the feasibility of highly accurate calculations based on a hybrid exchangecorrelation functional within the present framework. Our method provides a new way for studying lattice dynamics in severely anharmonic materials where the standard harmonic approximation and the perturbative approach break down.
I Introduction
Lattice vibration has been one of the most important research subject in the field of condensed matter physics and materials science because of its close connection with various thermodynamical, dynamical, and transport properties of solids including phase stability at finite temperature, lattice thermal conductivity (LTC), and superconducting critical temperature of phononmediated superconductors.Ziman (1960) With the advent of the increased computational power and developments of computational methods based on densityfunctional theory (DFT),Hohenberg and Kohn (1964); Kohn and Sham (1965) firstprinciples calculation of phonons and related properties is widely performed in order to interpret experimental results and to design new materials possessing desired physical properties.
While firstprinciples methods are convenient and reasonably accurate for various applications, their validity is always limited by approximations and assumptions made therein. In the case of phonon calculations, the harmonic approximation (HA) is usually adopted, in which only the secondderivative of the BornOppenheimer (BO) energy surface is considered assuming that atomic displacements are small enough compared with interatomic distances.Wallace (1972) The HA is in many cases valid and useful to obtain phonon dispersion curves and discuss phase stability based on the vibrational free energy. Baroni et al. (2001) However, it obviously fails to describe many important properties related to the lattice anharmonicity such as thermal expansion, LTC, and temperature and volumedependence of phonon frequencies, for which we must go beyond the HA.
When the cubic and higherorder anharmonic terms of the BO energy surface are sufficiently small compared with the harmonic one, the anharmonic effects can be treated by the manybody perturbation theory (PT). Maradudin and Fein (1962) The PT has successfully been employed to explain phonon linewidths in semiconductorsNarasimhan and Vanderbilt (1991); Debernardi et al. (1995) and to understand/predict LTC values of a wide variety of materials. Broido et al. (2007); Lindsay et al. (2010); Cepellotti et al. (2015); Jain and McGaughey (2015); Zeraati et al. (2016); Esfarjani et al. (2011); Shiomi et al. (2011); Tian et al. (2012); Li and Mingo (2014); Tadano et al. (2015); Carrete et al. (2014); Togo et al. (2015); Dekura et al. (2013); Seko et al. (2015a) In these calculations, the lowestorder contribution associated with the thirdorder anharmonic terms is considered to describe the intrinsic phononphonon scattering events that is of primary importance as a phonon scattering process in semiconductors and insulators. Despite its great success, the PT is often inadequate and sometimes breaks down completely in important applications. For example, in quantum crystals such as solid heliumHorner (1967) and the superconducting sulfur hydride under pressure,Drozdov et al. (2015); Errea et al. (2015) the zeropoint motion is so significant that the anharmonic renormalization of phonon frequencies must be considered beyond the PT. Even worse, the PT is not applicable to hightemperature phases of dielectric materials due to the existence of unstable phonon modes within the HA. This seriously limits the predictive power of the firstprinciples lattice dynamics methods for emergent thermoelectric materialsDelaire et al. (2011); Zhao et al. (2014) and hybrid perovskite solar cells.Kojima et al. (2009); Frost and Walsh (2016)
To overcome the aforementioned limitations, several DFTbased methods have recently been developed for treating anharmonic effects in solids nonperturbatively, which rely on the vibrational selfconsistentfield theory,Bowman (1978); Monserrat et al. (2013) the selfconsistent phonon (SCP) theory, Hooton (1958); Souvatzis et al. (2008); Errea et al. (2014); Tadano and Tsuneyuki (2015); van Roekeghem et al. (2016a) or ab initio molecular dynamics (AIMD) methods. Hellman et al. (2011); Sun et al. (2014) They have been employed to study anharmonic renormalization of phonon frequencies as well as its effects on the phase stability,Hellman et al. (2011); Engel et al. (2015) the superconducting critical temperature,Errea et al. (2011, 2013, 2015); Sano et al. (2016); Errea (2016) and the LTC of hightemperature phases.Tadano and Tsuneyuki (2015); van Roekeghem et al. (2016b) For example, the selfconsistent ab initio lattice dynamics (SCAILD)Souvatzis et al. (2008) and the stochastic selfconsistent harmonic approximation (SSCHA)Errea et al. (2014) are SCPbased approaches which can incorporate the effect of lattice anharmonicity at the meanfield level. In these methods, anharmonic phonon frequencies, or equivalently, effective harmonic force constants are obtained selfconsistently by repeatedly calculating atomic forces in supercells with suitably chosen atomic configurations. These stochastic algorithms are quite useful since a calculation of fourthorder anharmonic force constants is unnecessary. More recently, another efficient implementation of the SCP theory was developed,Tadano and Tsuneyuki (2015) which employs anharmonic force constants calculated using the compressive sensing lattice dynamics method.Zhou et al. (2014) The temperaturedependent effective potential (TDEP) methodHellman et al. (2011); Hellman and Abrikosov (2013) is an AIMDbased approach, in which effective harmonic and cubic force constants are extracted from displacementforce data sets calculated for atomic configurations sampled by AIMD. Although the AIMDbased approaches are not valid in the lowtemperature region due to their inability to account for the zeropoint motion, they should be accurate in the hightemperature region because anharmonic effects are fully included. With these constant efforts, accurate firstprinciples modeling of lattice vibration is becoming practical even for severely anharmonic materials.
In this paper, we review the recent developments of firstprinciples lattice dynamics methods for strongly anharmonic materials, particularly focusing on one of the SCPbased approaches developed by the authors.Tadano and Tsuneyuki (2015) Using microscopic force constants estimated by supercell methods as inputs,Esfarjani and Stokes (2008); Tadano et al. (2014); Zhou et al. (2014) our method can calculate temperaturedependent phonon frequencies and LTC efficiently. We demonstrate the validity of our method by applying it to cubic SrTiO, which is realized only in the high temperature region above 105 K.Yamada and Shirane (1969) In Sec. II, we present theoretical formulation and some numerical technologies that underlie our SCP approach. In Sec. III, we show the validity of the approach by carefully comparing our numerical results with available experimental data. In addition, the possibility of highly accurate calculation based on a hybrid exchangecorrelation functional is demonstrated. Finally, we make concluding remarks in Sec. IV.
Ii Methodology
ii.1 Taylor expansion potential
Let us start by expressing the potential energy of the interacting nuclear system as a Taylor expansion with respect to the atomic displacement:
(1)  
(2) 
Here, is the displacement of atom in the th unit cell along the direction and is the thorder derivative of the potential energy by displacements, which is usually termed interatomic force constants (IFCs). In Eq. (2), the IFCs are determined at the ground state and therefore have no temperature dependence.
In the harmonic approximation (HA), all of the anharmonic terms are neglected. Then, the Hamiltonian of the system can be represented as with the harmonic phonon frequency and associated displacement operator , where and are annihilation and creation operators of phonon with the crystal momentum and the branch index , respectively. The phonon frequency can be obtained from the dynamical matrix
(3) 
where is the mass of atom and is the translational vector of the th unit cell. By diagonalizing the matrix , one obtains the squared harmonic frequency and the polarization vector as . Since the harmonic IFCs are temperature independent, the intrinsic temperature dependence of phonon frequencies and polarization vectors is absent within the HA. The atomic displacement can be represented in terms of as follows:
(4) 
Here, is the number of points in the first Brillouin zone (BZ) and is a component of . By substituting Eq. (4) into Eq. (2), we obtain
(5) 
Here and in the following, we use for the shorthand notation of , satisfying and . The term is 1 if is a vector of the reciprocal lattice and 0 otherwise. The coefficient is defined as
(6)  
(7) 
Equation (5) is the thorder potential energy represented in terms of the harmonic displacement operator .
ii.2 Anharmonic perturbation theory
To describe the intrinsic phonon scattering processes and temperature dependence of phonon frequencies, we need to consider the anharmonic terms. When the anharmonic terms are sufficiently small compared with the harmonic one, we can treat them as perturbation to the noninteracting Hamiltonian as
(8) 
Here, we omitted fifth and higherorder terms since their contributions are much smaller than the cubic and quartic terms. Let be the onephonon Green’s function and be that of the noninteracting system . Then, the following Dyson equation holds:
(9) 
with being the anharmonic selfenergy which can be estimated within a systematic diagrammatic approximation. If the selfenergy correction is small and the condition is wellsatisfied, the phonon quasiparticle picture is still valid and the phonon frequency shift and the linewidth are given as and .
Figure 1 shows the first and secondorder selfenergy diagrams associated with the cubic and quartic anharmonic terms. The tadpole (a) and the bubble (b) are the secondorder diagrams resulting from the cubic terms and their explicit formulae are given as follows:Maradudin and Fein (1962)
(10)  
(11) 
Here, is the BoseEinstein distribution function and with being a positive infinitesimal. In addition, the summation in Eq. (11) is restricted to the pairs satisfying the momentum conservation .
The tadpole diagram is real and therefore give rise to the frequency shift only. In Eq. (10), we neglect the processes involving zonecenter acoustic modes as the intermediate state because of the singularity of the matrix element for acoustic modes. Therefore, the tadpole diagram only accounts for the frequency shift due to the relaxation of internal coordinates driven by the cubic anharmonicity. The neglected contributions involving acoustic modes can be treated separately by the quasiharmonic approximation.Paulatto et al. (2015) The bubble diagram causes the dominant contribution to the phonon linewidth and the quasiparticle lifetime can be estimated as .
The loop diagram (c) is the firstorder contribution from the quartic terms, whose expression is
(12) 
Since is real, the loop diagram solely causes the frequency shift. The expression for the higherorder diagrams (d) and (e) associated with the quartic terms are not shown in this paper, but they can be derived straightforwardly using the Feynman rule.Tripathi and Pathak (1974); Monga and Pathak (1978); Procacci et al. (1992)
Firstprinciples calculations of the bubble and tadpole diagrams have been performed for explaining the temperature dependence of Raman shift and linewidth of group IV and IIIIV semiconductors,Narasimhan and Vanderbilt (1991); Debernardi et al. (1995). Recently, there are many DFTbased estimation of the bubble diagram mostly for thermal conductivity prediction described in Sec. II.4. On the other hand, the calculation of the loop diagram has only been reported for relatively simple systems until very recently.Lang et al. (1999); Lazzeri et al. (2003); Bonini et al. (2007) This is mainly because the calculation of relies on the quartic IFCs that are technically more difficult to obtain based on DFT than the harmonic and cubic IFCs.Rousseau and Bergara (2010); Errea (2016) In Sec. II.5, we will introduce some recent technical developments for estimating the quartic terms.
ii.3 Selfconsistent phonon theory
When the anharmonic terms are comparable with the harmonic term, they need to be treated nonperturbatively. The selfconsistent phonon (SCP) theory originally developed by Hooton,Hooton (1958) followed by more elegant formulations based on the manybody theory,Koehler (1966); Horner (1967); Werthamer (1970) is one of the nonperturbative approach that can treat the anharmonic renormalization of phonon frequencies. To derive the SCP equation, we rewrite the Hamiltonian [Eq. (8)] as
(13) 
Here, is the effective harmonic Hamiltonian, which can be written as with the renormalized phonon frequency and the associated displacement operator . Then, the free energy of the original system can be written as the cumulant expansion
(14) 
Here, and with the partition function , and indicates the thorder cumulant of the operator . The firstorder SCP theory only considers the first cumulant in Eq. (14), in which the following GibbsBogoliubov inequality is satisfied:
(15) 
Since the variational principle holds in the firstorder SCP theory, the solution can be found by minimizing the righthand side of the above equation with respect to the adjustable parameters such as anharmonic frequencies , polarization vectors, and internal coordinates. For brevity of the derivation, let us suppose that the polarization vectors and internal coordinates are not altered by anharmonic effects. Then, from the condition with , the following SCP equation can be readily derived:
(16)  
(17) 
By solving Eqs. (16) and (17) selfconsistently, anharmonic phonon frequencies can be obtained. The equation is valid even if some of the harmonic phonon frequencies are unstable (). These unstable modes are renormalized by the second term to give a stable phonon mode (). Since the formulation assumes the existence of a welldefined Hamiltonian , imaginary frequency () is not allowed as a solution of the equation.
The same result has also been obtained based on the Green’s function approach. Hermes and Hirata (2013, 2014); Tadano and Tsuneyuki (2015) Diagrammatically, the SCP theory corresponds to include an infinite set of diagrams that can be generated from the loop diagram as shown in Fig.2(a). Therefore, the SCP theory automatically includes the contribution from the higherorder diagram shown by Fig.1(e). Substituting for in the second term of Eq. (16) gives the correct perturbation limit, where we have .
In the above derivation of Eqs. (16) and (17), we made an assumption that the polarization vectors do not change by anharmonic effects. Such an assumption is reasonable for simple systems containing a few atoms in the primitive cell,Errea et al. (2011) but not valid for more complex structures because the polarization mixing (PM) can be significant between phonon modes belonging to the same irreducible representation. To correctly account for the PM, we have made an extension of the SCP equation to include the offdiagonal elements of the loop selfenergy and developed an efficient algorithm for solving the equation.Tadano and Tsuneyuki (2015) We have applied the extended method to cubic SrTiO and demonstrated the importance of the PM for the soft modes at highsymmetry points.Tadano and Tsuneyuki (2015) The same implementation was also employed to study anharmonic effects in the superconducting hydrogen sulfides.Sano et al. (2016) Moreover, it is possible to extend the theory to include the tadpole diagram.Hermes and Hirata (2013)
The deterministic implementation based on the reciprocal formalism [Eqs. (16) and (17)] is efficient and applicable to relatively complex structures. Also, we can easily check the finite size effect by increasing the mesh in Eq. (17) based on the Fourier interpolation, which can be significant near the structural phase transition.Cowley (1996); Tadano and Tsuneyuki (2015) However, as already mentioned in Sec. II.2, firstprinciples calculation of the quartic IFCs necessary for has been a technical challenge. To avoid the cumbersome estimation of the quartic IFCs, several stochastic implementations have been developedSouvatzis et al. (2008); Errea et al. (2014); Errea (2016); van Roekeghem et al. (2016a) and applied to study anharmonic effects in simple metals,Souvatzis et al. (2008); Lian et al. (2015) light element superconductors,Errea et al. (2013, 2015) transition metal dichalcogenide,Leroux et al. (2015) and simple perovskite materials.van Roekeghem et al. (2016b) In these approaches, effective harmonic IFCs renormalized by anharmonic effects are estimated by randomly displacing atoms in a supercell. The stochastic methods should be more accurate than the present deterministic approach [Eqs. (16), (17)] because they don’t require the expansion of Eq. (1) nor the truncation of higherorder terms done in Eq. (8). However, since a large number of DFT calculations is generally needed for obtaining convergence, the stochastic approaches are more computationally intensive than the deterministic one.
ii.4 Lattice thermal conductivity
The lattice thermal conductivity (LTC) has been intensively studied as it plays an important role for optimizing the thermoelectric figureofmerit . Most of the theoretical predictions of LTC rely on either the Boltzmann transport equation (BTE) or the molecular dynamics (MD) method. According to the formulation of BTE for phonons developed by Pierls,Peierls (1955) LTC is given as
(18) 
where is the unit cell volume and is the group velocity. The vector is a solution of the following linear equation:
(19) 
Here, and are intrinsic transition rates of phonon involving two and three phonons, respectively. For example, includes the effect of phononisotope scatterings and is dominated by threephonon scattering processes associated with the cubic lattice anharmonicity. To obtain LTC of bulk solids, Eqs. (18) and (19) must be solved on a fine mesh in the first BZ. Broido et al. first solved the linearized BTE for bulk silicon and germanium based on DFT and obtained values that agree well with experimental data.Broido et al. (2007) Owing to the increased computational power and developments of efficient implementations, it is now feasible to solve Eq. (19) for relatively complex materials using iterative algorithms,Broido et al. (2007); Fugallo et al. (2013); Chernatynskiy and Phillpot (2010) or even by using the direct method for simple systems.Chaput (2013); Togo et al. (2015)
To make the computation more feasible, the relaxationtime approximation (RTA) is often employed in the literature, whereby LTC is simplified as
(20) 
In RTA, the current vertex corrections are neglected and the transport relaxation time is approximated by the quasiparticle lifetime .Sun and Allen (2010) Therefore, RTA incorrectly treat the normal process of the threephonon scattering as resistive, which results in . The underestimation by RTA is predominant in highLTC materials such as diamond Ward et al. (2009) as well as in a variety of twodimensional materials.Lindsay et al. (2010); Cepellotti et al. (2015); Jain and McGaughey (2015); Zeraati et al. (2016) Except for these cases, RTA usually yields results similar to the full solution to BTE. Hence, RTA has successfully been employed to analyze and predict LTC of various kinds of solids.Esfarjani et al. (2011); Shiomi et al. (2011); Tian et al. (2012); Li and Mingo (2014); Tadano et al. (2015); Carrete et al. (2014); Togo et al. (2015); Dekura et al. (2013); Seko et al. (2015a)
While the predictive accuracy of the conventional BTE [Eq. (18)] and BTERTA [Eq. (20)] is reasonably high, it has a limitation for severely anharmonic systems because the method treats the anharmonic effect perturbatively. Most importantly, the conventional approach neglects the temperature dependence of phonon frequencies and eigenvectors. This treatment prohibits us to estimate LTC of high temperature phases because of the imaginary modes within the HA. To overcome this limitation, we made an extension as follows:Tadano and Tsuneyuki (2015)
(21) 
where is the SCP frequency, is the group velocity of the renormalized phonon, and . The lifetime of the renormalized phonon associated with the threephonon scattering processes can be calculated from the diagram shown in Fig. 2(b). This is similar to the original bubble diagram [Fig. 1 (b) and Eq. (11)], but the harmonic quantities are replaced with the SCP solution . We have successfully applied the new method to predict LTC of cubic SrTiO including its unusual temperature dependence, which will be discussed in detail in Sec. III. Other previous studies employed a similar approach for studying phonon transport properties in BiTe,Hellman and Broido (2014) PbTe,Romero et al. (2015) and SrTiO,Feng et al. (2015) in which the effective harmonic and cubic IFCs were obtained by the temperaturedependent effective potential method.Hellman et al. (2011); Hellman and Abrikosov (2013) These studies also indicate the importance of the renormalization of IFCs in LTC of the studied materials.
Incorporating the fourphonon scattering process into BTE [Eq. (19)] is still a big computational challenge even within RTA. In order to achieve this, the higherorder diagram shown in Fig. 1(d) must be calculated with a dense mesh, which has only been reported for the group IV semiconductors based on empirical potentials. Feng and Ruan (2016) More work is thus needed to develop robust understanding of the fourphonon scattering effects on LTC of various kinds of solids including lowdimensional systems.
MD is another powerful tool for calculating LTC of solids. Since the anharmonicity is fully included, the MDbased methods, such as the nonequilibrium MDEvans (1982); MüllerPlathe (1997) and the GreenKubo methods,Kubo et al. (1957) are in principle more accurate than BTE at high temperatures if accurate AIMD can be performed. Moreover, they are suitable for systems with static and dynamical disorders, nanostructured devices, and aperiodic structures. To achieve a reliable prediction of LTC with MD, it is necessary to employ a fairly large supercell to sample longwavelength phonons that have dominant contributions to LTC. Schelling et al. (2002); He et al. (2012); Tadano et al. (2014) In addition, the simulation length must be long enough to reduce statistical errors. These requirements make the MD methods very costly especially for highLTC materials due to large meanfreepaths and long relaxation times of phonons. Therefore, straightforward applications of AIMD have only been attempted for a few materials having low LTC.Stackhouse et al. (2010); *Stackhouse:2015ew; Yue et al. (2016) To make the computation feasible, classical force fields are commonly employed, but they can be a source of inaccuracy. To address this problem, many efforts have been made to improve the accuracy of classical force fields or to develop a new force field by using DFT results as training data.Lindsay and Broido (2010); Hata et al. (2016); Behler and Parrinello (2007); Bartók et al. (2010)
ii.5 Firstprinciples calculation of force constants
To conduct the SCP calculations described in this section, harmonic and quartic IFCs are necessary as inputs. Cubic terms are also required for performing the BTE calculations. For that purpose, two conceptually different methods are commonly employed: densityfunctional perturbation theory (DFPT) and the supercell method.
DFPT is one of the standard ab initio method for calculating harmonic IFCs.Giannozzi et al. (1991); Gonze (1997); Baroni et al. (2001) In DFPT, harmonic dynamical matrix and linear electronphonon coupling coefficients for an arbitrary point can be obtained. Recently, an efficient implementation based on DFPT has also been developed for cubic IFCs.Paulatto et al. (2013) While DFPT is very accurate and efficient since it doesn’t involve expensive supercell calculations, its implementation is relatively complicated and its extension to the quartic and higherorder IFCs is still challenging.
The supercell method is another approach to calculate IFCs, which is more expensive than DFPT but easier to implement. In this approach, harmonic IFCs are estimated from the firstorder numerical derivative of the atomic forces as
(22) 
or inversely asParlinski et al. (1997)
(23) 
which are valid when the displacement is small, say 0.01 Å. is the atomic force acting on the th atom in the th cell along the direction when the displacement of is made with all other atoms being kept at the equilibrium positions. Equation (23) is defined for every single coefficient , many of which are linearly dependent with each other because of the space group symmetry, the permutation symmetry, and the lattice periodicity.Esfarjani and Stokes (2008); Tadano et al. (2014) For the convenience of the following discussion, let us introduce the column vectors and , which consist of atomic displacements and forces in the supercell containing atoms, respectively. We also denote the th displacement pattern as and the associated force vector as . Then, denoting by a column vector comprising of unique harmonic coefficients, a set of the linear equation [Eq. (23)] can be symbolically written as
(24) 
Here, is a matrix defined as , is a matrix, , and is the number of displacement patterns. Please note that Eq. (24) is valid if the displacements are small enough to neglect anharmonic effects in . If we consider sufficient displacement patterns and make the matrix full rank, the harmonic IFCs can be obtained, in the leastsquares sense, asParlinski et al. (1997); Chaput et al. (2011)
(25) 
where is the pseudoinverse of the matrix , which can be readily obtained by the singular value decomposition. A similar equation can be formed for the cubic terms asChaput et al. (2011)
(26) 
in which is the cubic contribution to the atomic forces and with defined as . Again, Eq. (26) assumes the quartic and higherorder contributions to atomic forces are negligible, which is valid when the displacement is small. To make the matrix full rank, we need to consider many displacement patterns where at least two atoms must be displaced simultaneously. Therefore, their calculation is more expensive than the harmonic terms. To reduce the computational burden, it is common to introduce a cutoff distance and/or employ a smaller supercell for cubic IFCs. These treatments are well justified if the anharmonic interaction is shortranged, which is usually the case in covalent compounds. By contrast, longrange cubic interaction may be significant in ionic compounds,Li et al. (2015) for which a careful choice of the cutoff radius must be made to avoid erroneous numerical results.
Following the above procedure, it is straightforward to form a linear equation for quartic and higherorder IFCs, but they are difficult to solve because of a large number of DFT calculations necessary for making the matrix full rank. Moreover, an accurate DFT calculation of is difficult because they are very small in a small displacement limit. While this numerical issue may be avoided by using a larger displacement, finding an appropriate value of displacements is formidable as it depends on atomic environments. Performing AIMD simulation is a reasonable solution to this problem as it can generate physically relevant atomic configurations automatically. From displacementforce data set obtained from the AIMD trajectory, we may estimate IFCs by solving the ordinary leastsquares (OLS) problemTadano et al. (2014)
(27) 
Here, is the column vector defined as and with being . Given a sufficient number of AIMD snapshots, the unknown IFCs can be estimated as . Using this approach, we successfully extracted anharmonic IFCs of Si and MgSi from AIMD trajectories.Tadano et al. (2014) Recently, Zhou et al. proposed a more sophisticated approach called the compressed sensing lattice dynamics (CSLD) method,Zhou et al. (2014) where they estimated IFCs by the least absolute shrinkage and selection operator (LASSO)Tibshirani (1996):
(28) 
Owing to the regularization term, physically irrelevant IFCs are driven to be exactly zero and important terms are selected and calculated automatically. The method is robust against random and systematic noiseCandès et al. (2006) and the penalty term helps to avoid the overfitting inherent in OLS. The coefficient is a hyperparameter that controls the tradeoff between sparsity and accuracy of the model; all IFCs will be zero in the large limit and corresponds to OLS. An appropriate value of can be estimated, for instance, by crossvalidation as will be demonstrated in Sec. III.
In Ref. Zhou et al., 2014, Zhou et al. extracted harmonic and anharmonic IFCs of Si, NaCl, and CuSbS using Eq. (28) combined with AIMD simulation. Based on the IFCs estimated by LASSO, they calculated LTC of these materials and obtained good agreements with experimental data. Using the same approach, we successfully extracted anharmonic IFCs of SrTiO,Tadano and Tsuneyuki (2015) superconducting hydrogen sulfides,Sano et al. (2016) and thermoelectric clathrates and applied them for SCP and BTE calculations. Also, compressed sensing has successfully been employed for developing interatomic potentialsSeko et al. (2014); *Seko:2015gq and cluster expansion modelNelson et al. (2013) based on DFT.
Iii Applications to Strontium Titanate
In Ref. Tadano and Tsuneyuki, 2015, we calculated anharmonic phonon properties of cubic SrTiO (cSTO), which is stable above 105 K. All of the DFT calculations were performed using the VASP codeKresse and Furthmüller (1996); Kresse and Joubert (1999), where the generalized gradient approximation for solidsPerdew et al. (2008) (PBEsol) was employed for the exchange correlation functional. This choice was made because the PBEsol functional is known to predict the lattice constants for a wide class of materials, which is crucial for the accurate calculation of phonon properties.Wahl et al. (2008); Skelton et al. (2014) Using quartic IFCs estimated by LASSO, we performed SCP simulations and obtained overall good agreements with experimentally observed phonon spectra at room temperature. However, the quantitative accuracy based on the PBEsol functional was not very satisfactory because the calculated frequencies of the ferroelectric (antiferrodistortive) soft modes were overestimated (underestimated).
In this paper, we report more accurate computational result within the HeydScuseriaErnzerhof hybrid functional (HSE06, hereafter HSE)Heyd et al. (2003); *Heyd:2006dc. Before performing phonon calculations, we optimized the lattice constant of cSTO using the HSE functional and obtained 3.900 Å, in good agreements with the experimental value of 3.905 Å (Ref. Okazaki and Kawaminami, 1973, 293 K) and the previous HSE result of 3.904 Å.Wahl et al. (2008) The computational conditions employed in this study are basically the same as the previous PBEsol study,Tadano and Tsuneyuki (2015) except that BZ integration was performed with the 888 MonkhorstPack point grid. All of the phonon calculations reported here are performed with a 222 cubic supercell (40 atoms) using the ALAMODE package,Tadano (); Tadano et al. (2014) an open source software developed for modeling lattice thermal conductivity and anharmonicity.
iii.1 Anharmonic force constants by LASSO
To develop an accurate model based on Eq. (1), we considered anharmonic terms up to the sixth order. Here, all possible cubic terms occurring in the 222 supercell were considered. For quartic IFCs, we only considered onsite, twobody, and threebody IFCs inside a cutoff length of 6 Å and neglected less important fourbody terms. Also, only onsite and twobody terms inside the same cutoff radius were considered for the fifth and sixthorder IFCs. We then determined a set of linearly independent parameters by making full use of the available symmetry operations and the constraints for the translational invariance.Esfarjani and Stokes (2008); Tadano et al. (2014) The numbers of independent parameters were found to be 698, 2215, 43, and 125 for cubic, quartic, quintic, and sixtic terms, respectively, which in total form a parameter vector of length 3018. Calculating all elements of would be unfeasible if the conventional finite displacement method were employed as it would need a few thousand DFT calculations. To determined these IFCs by LASSO, we prepared 40 displacementforce data set using AIMD, whose detailed procedure is described in Ref. Tadano and Tsuneyuki, 2015. After that, the data set was split into 4 smaller subsets and the hyperparameter in Eq. (28) was selected by the fourfold crossvalidation.Hastie et al. (2009) The solution path was obtained by using the coordinate descent method,Hastie et al. (2015) for which each column vector of the matrix was standardized beforehand. Upon solving Eq. (28), harmonic IFCs were fixed to the values obtained by OLS.
Figure 3 shows the result of the crossvalidation obtained for PBEsol. Since very similar results were obtained for HSE, they are not shown here. The data shown in the top panel is the relative error of the atomic forces, which is defined as the square root of . In the large region, the difference between the training error and the validation error is marginal. With decreasing , the difference becomes more predominant; the training error keep decreasing, whereas the validation curve has its minimum around indicated by the dotted vertical line in the figure. This value of was selected as an optimal choice because it is expected to give best prediction accuracy against independent data sets. The bottom panel shows the number of nonzero coefficients. With the optimal value of , we obtained 2024 nonzero coefficients which is about 67% of the total elements. The accuracy of the extracted IFCs was carefully verified by applying them to the independent atomic configurations sampled by AIMD.
iii.2 SCP solution with HSE
Using the calculated harmonic and quartic IFCs, we conducted SCP calculations based on Eqs. (16) and (17). To correctly describe the PM, we considered offdiagonal elements of the loop diagram as done in Ref. Tadano and Tsuneyuki, 2015. The convergence of the anharmonic frequencies at the commensurate 222 point grid was carefully checked with respect to the intermediate grid in Eq. (17) and found sufficiently converged results with 888 points. In addition, the nonanalytic correction to the harmonic dynamical matrix was included to reproduce the LOTO splitting around the point. To this end, we employed the Born effective charges and the dielectric constant reported in Ref. Tadano and Tsuneyuki, 2015 for PBEsol. Since calculating these quantities based on HSE was very expensive, we approximately employed the Born effective charges obtained by PBEsol and the experimental dielectric constant of 5.18Zhong et al. (1994) in the following HSE calculations.
Figure 4 shows the phonon dispersion curves of cSTO calculated within the HA (dotted lines) and the SCP theory (thin solid lines). For the SCP results, we show several curves corresponding to different temperatures ranging from 200 K (0 K) to 1000 K for PBEsol (HSE) in steps of 100 K. For ease of the comparison with inelastic neutron scattering (INS) data measured at 300 K,Stirling (1972); Cowley et al. (1969) the theoretical curves at 300 K are highlighted as thick solid lines. In the harmonic phonon dispersion, unstable phonon modes () occur at and points, which correspond to the ferroelectric (FE) and antiferrodistortive (AFD) modes, respectively. The harmonic phonon frequency of the FE mode is 101 cm with HSE, whose absolute value is larger than the PBEsol result of 58 cm, indicating the enhanced instability of the FE mode by the Fock exchange term. For the AFD mode, the trend was opposite and the frequency changed from 76 cm (PBEsol) to 25 cm (HSE) as summarized in Table 1. These tendencies are consistent with the previous DFT study of Wahl et al.Wahl et al. (2008)
As can be seen in Fig. 4, the quartic anharmonicity generally increases the vibrational frequencies in cSTO, which is particularly predominant in the lowenergy optical modes. The PBEsol result obtained here is almost identical with our previous result in Ref. Tadano and Tsuneyuki, 2015. The SCP frequency of the FE mode is 133 cm at 300 K with PBEsol, which overestimates the INS data of 87.15.6 cm (Ref. Yamada and Shirane, 1969, 293 K) and the IR data of 89 cm (Ref. Servoin et al., 1980, 300 K). By contrast, the frequency of the AFD mode at R point is underestimated with PBEsol. We have found that the underestimation at R point can be partially cured by employing HSE as shown in the right panel of Fig. 4. However, was improved only slightly by HSE and the SCP frequency of 128 cm still overestimates the experimental values. In what follows, we discuss of origin of the disagreement by going beyond the SCP theory.
PBEsol  HSE  Expt.  

HA  SCP  SCP + B  Peak in  HA  SCP  SCP + B  Peak in  
(FE)  58 ^{a)}  133  120  120  101  128  109  104  87.15.6 ^{b)}, 89 ^{c)} 
(AFD)  76 ^{a)}  31  2  18  25  66  53  52  48.71.8 ^{d)} 

Ref. Tadano and Tsuneyuki, 2015

Ref. Yamada and Shirane, 1969 (INS, 293 K)

Ref. Servoin et al., 1980 (IR, 300 K)

Ref. Cowley et al., 1969 (INS, 300 K)
iii.3 Frequency shift by cubic anharmonicity
Here, we discuss the effect of the cubic anharmonicity on phonon frequencies. In the present SCP calculations, only the quartic anharmonicity is included as diagrammatically shown in Fig. 2(a). However, frequency shifts due to the cubic terms may not be negligible. To unveil its effect quantitatively, we calculated the bubble selfenergy shown in Fig. 2(b) and estimated the frequency shift of the FE mode from . The results of are shown in Table 1 as “SCP+B” (B stands for bubble). The correction was found to be negative and significant for both of the FE and AFD modes at 300 K. For the FE mode, it was about cm with PBEsol and cm with HSE. For the AFD mode, the correction was very large cm with PBEsol and the corrected frequency was almost zero. It is interesting to note that the frequency dependence of was marginal in the lowfrequency region ( cm) and the difference between and was minor especially if the SCP frequency was small. Therefore, it may be reasonable to make the static approximation as in order to calculate the frequency shift due to the bubble diagram, which was done in Ref. Shukla and Cowley, 1998. With the correction from the bubble diagram, the theoretical values based on HSE agrees better with the experimental values. Although other theoretical and numerical aspects must be carefully investigated, our results represent the importance of accurate theoretical and computational modeling of anharmonic effects for robust understandings of lattice dynamics in severely anharmonic materials.
To make a straightforward comparison with the experimental data, we also calculated the spectral function using the following formula:
(29) 
Since has a direct connection with the dynamical structural factor measured in INS experiments, it should be more reasonable and accurate to estimate theoretical phonon frequencies from peak positions of Eq. (29). The estimated frequencies are also shown in Table 1, which differ from the “SCP+B” results. This indicates the limited accuracy of the “SCP+B” result when the condition of is not satisfied well. Figure 5 shows the spectral function at 300 K calculated along the highsymmetry lines of BZ. For comparison, we also show the SCP dispersion curves by dashed lines. It is evident from the figure that the AFD mode and highenergy optical modes above cm are strongly damped by the cubic anharmonicity. Moreover, the frequency shift by the bubble diagram is most significant in the lowlying optical modes around the point and along the RM path.
iii.4 Anomalous thermal transport in SrTiO
Finally, we calculated LTC of cSTO using Eq. (21) with 202020 grid points. In Fig. 6, we compare our computational results with experimental data from Refs. Muta et al., 2005; Popuri et al., 2014. As can be seen in the figure, we obtained overall good agreements with experimental LTC including its anomalous temperature dependence. Unlike the case of soft mode frequencies discussed in the previous section, the improvement made by the accurate HSE functional was minor. This occurs because LTC is a quantity averaged over the first BZ, for which the prediction accuracy of the less expensive PBEsol functional should be sufficient.
The anomalous temperature dependence in LTC (, 0.6–0.7) can be attributed to the hardening of phonon frequencies caused by the quartic anharmonicity. When we employed the SCP frequencies at 200 K for BTE calculations at higher temperatures, we obtained the conventional temperature dependence of as shown by the dashed line in Fig. 6. However, such a treatment cannot be justified because the LTC values are significantly underestimated at high temperatures, which was also discussed in Refs. Romero et al., 2015; Stackhouse et al., 2015; van Roekeghem et al., 2016a. We think the effect of the frequency renormalization by the quartic terms can be significant in other ultralowLTC materials such as SnSe and clathrate, which will be a topic of future work.
It is interesting to note that in cSTO more than 70% of the total thermal conductivity is carried by highfrequency optical modes above 150 cm,Feng et al. (2015); Tadano and Tsuneyuki (2015) as can be seen in Fig. 6 inset. Here, is the thermal conductivity spectrum defined as with being the contribution from each phonon mode . This can mainly be attributed to the large group velocities of these optical modes, which are larger than those of the acoustic modes.
Iv Summary and Conclusions
We reviewed our recent development of a firstprinciples framework to study lattice dynamical properties of strongly anharmonic crystals. By combining an efficient implementation based on the selfconsistent phonon (SCP) theory and the compressed sensing of anharmonic force constants, phonon frequencies renormalized by the quartic anharmonicity can be calculated nonperturbatively at various temperatures. In addition, the intrinsic phonon linewidths and frequency shifts by the cubic anharmonicity can be calculated by considering the bubble diagram on top of the SCP latticedynamics wavefunctions, which is actually essential for predicting lattice thermal conductivity (LTC) of anharmonic crystals from first principles.
We demonstrated the high predictive accuracy of the developed computational approach by carefully comparing soft mode frequencies and LTC values of cubic SrTiO with available experimental data, for which overall good agreements were obtained. In addition, we showed that the accurate, albeit expensive, calculation based on the HSE hybrid functional was feasible within the present framework, where the remarkable improvement was achieved in the frequency of the antiferrodistortive soft mode. We expect the presented approach to open up a wide range of applications and provides more insight into anharmonic effects in severely anharmonic materials, such as ferroelectric, thermoelectric, and photovoltaic materials, as well as in lightelement superconductors where the conventional harmonic approximation and the perturbative treatments break down.
Acknowledgements
This study is partly supported by JSPS KAKENHI Grant Number 16K17724 and “Materials research by Information Integration” Initiative (MI2I) project of the Support Program for Starting Up Innovation Hub from Japan Science and Technology Agency (JST). The computation in this work has been done using the facilities of the Supercomputer Center, Institute for Solid State Physics, The University of Tokyo.
References
 Ziman (1960) J. M. Ziman, Electrons and Phonons (Oxford University Press, 1960).
 Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
 Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
 Wallace (1972) D. C. Wallace, Thermodynamics of Crystals (Dover Publications, INC., 1972).
 Baroni et al. (2001) S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
 Maradudin and Fein (1962) A. A. Maradudin and A. E. Fein, Phys. Rev. 128, 2589 (1962).
 Narasimhan and Vanderbilt (1991) S. Narasimhan and D. Vanderbilt, Phys. Rev. B 43, 4541 (1991).
 Debernardi et al. (1995) A. Debernardi, S. Baroni, and E. Molinari, Phys. Rev. Lett. 75, 1819 (1995).
 Broido et al. (2007) D. A. Broido, M. Malorny, G. Birner, N. Mingo, and D. A. Stewart, Appl. Phys. Lett. 91, 231922 (2007).
 Lindsay et al. (2010) L. Lindsay, D. A. Broido, and N. Mingo, Phys. Rev. B 82, 115427 (2010).
 Cepellotti et al. (2015) A. Cepellotti, G. Fugallo, L. Paulatto, M. Lazzeri, F. Mauri, and N. Marzari, Nature Commun. 6, 6400 (2015).
 Jain and McGaughey (2015) A. Jain and A. J. H. McGaughey, Sci. Rep. 5, 8501 (2015).
 Zeraati et al. (2016) M. Zeraati, S. M. Vaez Allaei, I. Abdolhosseini Sarsari, M. Pourfath, and D. Donadio, Phys. Rev. B 93, 085424 (2016).
 Esfarjani et al. (2011) K. Esfarjani, G. Chen, and H. T. Stokes, Phys. Rev. B 84, 085204 (2011).
 Shiomi et al. (2011) J. Shiomi, K. Esfarjani, and G. Chen, Phys. Rev. B 84, 104302 (2011).
 Tian et al. (2012) Z. Tian, J. Garg, K. Esfarjani, T. Shiga, J. Shiomi, and G. Chen, Phys. Rev. B 85, 184303 (2012).
 Li and Mingo (2014) W. Li and N. Mingo, Phys. Rev. B 90, 094302 (2014).
 Tadano et al. (2015) T. Tadano, Y. Gohda, and S. Tsuneyuki, Phys. Rev. Lett. 114, 095501 (2015).
 Carrete et al. (2014) J. Carrete, W. Li, N. Mingo, S. Wang, and S. Curtarolo, Phys. Rev. X 4, 011019 (2014).
 Togo et al. (2015) A. Togo, L. Chaput, and I. Tanaka, Phys. Rev. B 91, 094306 (2015).
 Dekura et al. (2013) H. Dekura, T. Tsuchiya, and J. Tsuchiya, Phys. Rev. Lett. 110, 025904 (2013).
 Seko et al. (2015a) A. Seko, A. Togo, H. Hayashi, K. Tsuda, L. Chaput, and I. Tanaka, Phys. Rev. Lett. 115, 205901 (2015a).
 Horner (1967) H. Horner, Z. Phys. 205, 72 (1967).
 Drozdov et al. (2015) A. P. Drozdov, M. I. Eremets, I. A. Troyan, V. Ksenofontov, and S. I. Shylin, Nature 525, 73 (2015).
 Errea et al. (2015) I. Errea, M. Calandra, C. J. Pickard, J. Nelson, R. J. Needs, Y. Li, H. Liu, Y. Zhang, Y. Ma, and F. Mauri, Phys. Rev. Lett. 114, 157004 (2015).
 Delaire et al. (2011) O. Delaire, J. Ma, K. Marty, A. F. May, M. A. McGuire, M. H. Du, D. J. Singh, A. Podlesnyak, G. Ehlers, M. D. Lumsden, and B. C. Sales, Nature Mater. 10, 614 (2011).
 Zhao et al. (2014) L.D. Zhao, S.H. Lo, Y. Zhang, H. Sun, G. Tan, C. Uher, C. Wolverton, V. P. Dravid, and M. G. Kanatzidis, Nature 508, 373 (2014).
 Kojima et al. (2009) A. Kojima, K. Teshima, Y. Shirai, and T. Miyasaka, J. Am. Chem. Soc. 131, 6050 (2009).
 Frost and Walsh (2016) J. M. Frost and A. Walsh, Acc. Chem. Res 49, 528 (2016).
 Bowman (1978) J. M. Bowman, J. Chem. Phys. 68, 608 (1978).
 Monserrat et al. (2013) B. Monserrat, N. D. Drummond, and R. J. Needs, Phys. Rev. B 87, 144302 (2013).
 Hooton (1958) D. J. Hooton, Philos. Mag. 3, 49 (1958).
 Souvatzis et al. (2008) P. Souvatzis, O. Eriksson, M. I. Katsnelson, and S. P. Rudin, Phys. Rev. Lett. 100, 095901 (2008).
 Errea et al. (2014) I. Errea, M. Calandra, and F. Mauri, Phys. Rev. B 89, 064302 (2014).
 Tadano and Tsuneyuki (2015) T. Tadano and S. Tsuneyuki, Phys. Rev. B 92, 054301 (2015).
 van Roekeghem et al. (2016a) A. van Roekeghem, J. Carrete, and N. Mingo, Phys. Rev. B 94, 020303 (2016a).
 Hellman et al. (2011) O. Hellman, I. A. Abrikosov, and S. I. Simak, Phys. Rev. B 84, 180301 (2011).
 Sun et al. (2014) T. Sun, D.B. Zhang, and R. M. Wentzcovitch, Phys. Rev. B 89, 094109 (2014).
 Engel et al. (2015) E. A. Engel, B. Monserrat, and R. J. Needs, Phys. Rev. X 5, 021033 (2015).
 Errea et al. (2011) I. Errea, B. Rousseau, and A. Bergara, Phys. Rev. Lett. 106, 165501 (2011).
 Errea et al. (2013) I. Errea, M. Calandra, and F. Mauri, Phys. Rev. Lett. 111, 177002 (2013).
 Sano et al. (2016) W. Sano, T. Koretsune, T. Tadano, R. Akashi, and R. Arita, Phys. Rev. B 93, 094525 (2016).
 Errea (2016) I. Errea, Eur. Phys. J. B 89, 237 (2016).
 van Roekeghem et al. (2016b) A. van Roekeghem, J. Carrete, C. Oses, S. Curtarolo, and N. Mingo, Phys. Rev. X 6, 041061 (2016b).
 Zhou et al. (2014) F. Zhou, W. Nielson, Y. Xia, and V. Ozoliņš, Phys. Rev. Lett. 113, 185501 (2014).
 Hellman and Abrikosov (2013) O. Hellman and I. A. Abrikosov, Phys. Rev. B 88, 144301 (2013).
 Esfarjani and Stokes (2008) K. Esfarjani and H. T. Stokes, Phys. Rev. B 77, 144112 (2008).
 Tadano et al. (2014) T. Tadano, Y. Gohda, and S. Tsuneyuki, J. Phys: Condens. Matter 26, 225402 (2014).
 Yamada and Shirane (1969) Y. Yamada and G. Shirane, J. Phys. Soc. Jpn. 26, 396 (1969).
 Paulatto et al. (2015) L. Paulatto, I. Errea, M. Calandra, and F. Mauri, Phys. Rev. B 91, 054304 (2015).
 Tripathi and Pathak (1974) R. S. Tripathi and K. N. Pathak, Il Nuovo Cimento B 21, 289 (1974).
 Monga and Pathak (1978) M. R. Monga and K. N. Pathak, Phys. Rev. B 18, 5859 (1978).
 Procacci et al. (1992) P. Procacci, G. Cardini, R. Righini, and S. Califano, Phys. Rev. B 45, 2113 (1992).
 Lang et al. (1999) G. Lang, K. Karch, M. Schmitt, P. Pavone, A. P. Mayer, R. K. Wehner, and D. Strauch, Phys. Rev. B 59, 6182 (1999).
 Lazzeri et al. (2003) M. Lazzeri, M. Calandra, and F. Mauri, Phys. Rev. B 68, 220509 (2003).
 Bonini et al. (2007) N. Bonini, M. Lazzeri, N. Marzari, and F. Mauri, Phys. Rev. Lett. 99, 176802 (2007).
 Rousseau and Bergara (2010) B. Rousseau and A. Bergara, Phys. Rev. B 82, 104504 (2010).
 Koehler (1966) T. R. Koehler, Phys. Rev. Lett. 17, 89 (1966).
 Werthamer (1970) N. R. Werthamer, Phys. Rev. B 1, 572 (1970).
 Hermes and Hirata (2013) M. R. Hermes and S. Hirata, J. Phys. Chem. A 117, 7179 (2013).
 Hermes and Hirata (2014) M. R. Hermes and S. Hirata, Int. Rev. Phys. Chem. 34, 71 (2014).
 Cowley (1996) E. R. Cowley, Physica A 232, 585 (1996).
 Lian et al. (2015) C.S. Lian, J.T. Wang, and C. Chen, Phys. Rev. B 92, 184110 (2015).
 Leroux et al. (2015) M. Leroux, I. Errea, M. Le Tacon, S.M. Souliou, G. Garbarino, L. Cario, A. Bosak, F. Mauri, M. Calandra, and P. Rodière, Phys. Rev. B 92, 140303 (2015).
 Peierls (1955) R. E. Peierls, Quantum Theory of Solids (Oxford University Press, 1955).
 Fugallo et al. (2013) G. Fugallo, M. Lazzeri, L. Paulatto, and F. Mauri, Phys. Rev. B 88, 045430 (2013).
 Chernatynskiy and Phillpot (2010) A. Chernatynskiy and S. R. Phillpot, Phys. Rev. B 82, 134301 (2010).
 Chaput (2013) L. Chaput, Phys. Rev. Lett. 110, 265506 (2013).
 Sun and Allen (2010) T. Sun and P. B. Allen, Phys. Rev. B 82, 224305 (2010).
 Ward et al. (2009) A. Ward, D. A. Broido, D. A. Stewart, and G. Deinzer, Phys. Rev. B 80, 125203 (2009).
 Hellman and Broido (2014) O. Hellman and D. A. Broido, Phys. Rev. B 90, 134309 (2014).
 Romero et al. (2015) A. H. Romero, E. K. U. Gross, M. J. Verstraete, and O. Hellman, Phys. Rev. B 91, 214310 (2015).
 Feng et al. (2015) L. Feng, T. Shiga, and J. Shiomi, Appl. Phys. Express 8, 071501 (2015).
 Feng and Ruan (2016) T. Feng and X. Ruan, Phys. Rev. B 93, 045202 (2016).
 Evans (1982) D. J. Evans, Phys. Lett. A 91, 457 (1982).
 MüllerPlathe (1997) F. MüllerPlathe, J. Chem. Phys. 106, 6082 (1997).
 Kubo et al. (1957) R. Kubo, M. Yokota, and S. Nakajima, J. Phys. Soc. Jpn. 12, 1203 (1957).
 Schelling et al. (2002) P. K. Schelling, S. R. Phillpot, and P. Keblinski, Phys. Rev. B 65, 144306 (2002).
 He et al. (2012) Y. He, I. Savić, D. Donadio, and G. Galli, Phys. Chem. Chem. Phys. 14, 16209 (2012).
 Stackhouse et al. (2010) S. Stackhouse, L. Stixrude, and B. B. Karki, Phys. Rev. Lett. 104, 208501 (2010).
 Stackhouse et al. (2015) S. Stackhouse, L. Stixrude, and B. B. Karki, Earth Planet. Sci. Lett. 427, 11 (2015).
 Yue et al. (2016) S.Y. Yue, X. Zhang, G. Qin, J. Yang, and M. Hu, Phys. Rev. B 94, 115427 (2016).
 Lindsay and Broido (2010) L. Lindsay and D. A. Broido, Phys. Rev. B 81, 205441 (2010).
 Hata et al. (2016) T. Hata, G. Giorgi, and K. Yamashita, Nano Lett. 16, 2749 (2016).
 Behler and Parrinello (2007) J. Behler and M. Parrinello, Phys. Rev. Lett. 98, 146401 (2007).
 Bartók et al. (2010) A. P. Bartók, M. C. Payne, R. Kondor, and G. Csányi, Phys. Rev. Lett. 104, 136403 (2010).
 Giannozzi et al. (1991) P. Giannozzi, S. de Gironcoli, P. Pavone, and S. Baroni, Phys. Rev. B 43, 7231 (1991).
 Gonze (1997) X. Gonze, Phys. Rev. B 55, 10337 (1997).
 Paulatto et al. (2013) L. Paulatto, F. Mauri, and M. Lazzeri, Phys. Rev. B 87, 214303 (2013).
 Parlinski et al. (1997) K. Parlinski, Z. Q. Li, and Y. Kawazoe, Phys. Rev. Lett. 78, 4063 (1997).
 Chaput et al. (2011) L. Chaput, A. Togo, I. Tanaka, and G. Hug, Phys. Rev. B 84, 094302 (2011).
 Li et al. (2015) C. W. Li, J. Hong, A. F. May, D. Bansal, S. Chi, T. Hong, G. Ehlers, and O. Delaire, Nature Phys. 11, 1063 (2015).
 Tibshirani (1996) R. Tibshirani, J. R. Stat. Soc. Series B Stat. Methodol. 58, 267 (1996).
 Candès et al. (2006) E. J. Candès, J. K. Romberg, and T. Tao, Comm. Pure Appl. Math. 59, 1207 (2006).
 Seko et al. (2014) A. Seko, A. Takahashi, and I. Tanaka, Phys. Rev. B 90, 024101 (2014).
 Seko et al. (2015b) A. Seko, A. Takahashi, and I. Tanaka, Phys. Rev. B 92, 054113 (2015b).
 Nelson et al. (2013) L. J. Nelson, G. L. W. Hart, F. Zhou, and V. Ozoliņš, Phys. Rev. B 87, 035125 (2013).
 Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
 Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
 Perdew et al. (2008) J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Phys. Rev. Lett. 100, 136406 (2008).
 Wahl et al. (2008) R. Wahl, D. Vogtenhuber, and G. Kresse, Phys. Rev. B 78, 104116 (2008).
 Skelton et al. (2014) J. M. Skelton, S. C. Parker, A. Togo, I. Tanaka, and A. Walsh, Phys. Rev. B 89, 205203 (2014).
 Heyd et al. (2003) J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 118, 8207 (2003).
 Heyd et al. (2006) J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 124, 219906 (2006).
 Okazaki and Kawaminami (1973) A. Okazaki and M. Kawaminami, Mater. Res. Bull. 8, 545 (1973).
 (106) T. Tadano, “ALAMODE,” http://sourceforge.net/projects/alamode/.
 Hastie et al. (2009) T. Hastie, R. Tibshirani, and J. Friedman, “The Elements of Statistical Learning,” Springer New York, New York, NY (2009).
 Hastie et al. (2015) T. Hastie, R. Tibshirani, and M. Wainwright, Statistical Learning with Sparsity: The Lasso and Generalizations (Chapman & Hall/CRC, 2015).
 Stirling (1972) W. G. Stirling, J. Phys. C 5, 2711 (1972).
 Cowley et al. (1969) R. A. Cowley, W. J. L. Buyers, and G. Dolling, Solid State Commun. 7, 181 (1969).
 Zhong et al. (1994) W. Zhong, R. D. KingSmith, and D. Vanderbilt, Phys. Rev. Lett. 72, 3618 (1994).
 Servoin et al. (1980) J. L. Servoin, Y. Luspin, and F. Gervais, Phys. Rev. B 22, 5501 (1980).
 Shukla and Cowley (1998) R. C. Shukla and E. R. Cowley, Phys. Rev. B 58, 2596 (1998).
 Muta et al. (2005) H. Muta, K. Kurosaki, and S. Yamanaka, J. Alloys Compd. 392, 306 (2005).
 Popuri et al. (2014) S. R. Popuri, A. J. M. Scott, R. A. Downie, M. A. Hall, E. Suard, R. Decourt, M. Pollet, and J.W. G. Bos, RSC Adv. 4, 33720 (2014).