First-Principles Lattice Dynamics Method for Strongly Anharmonic Crystals

We review our recent development of a first-principles lattice dynamics method that can treat anharmonic effects nonperturbatively. The method is based on the self-consistent phonon theory and temperature-dependent 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 exchange-correlation 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 subjects 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 temperatures, lattice thermal conductivity (LTC), and the superconducting critical temperature of phonon-mediated superconductors. 1 With the advent of increased computational power and the development of computational methods based on density functional theory (DFT), 2,3 the first-principles 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 first-principles methods are convenient and reasonably accurate for various applications, their validity is always limited by the approximations and assumptions made therein. In the case of phonon calculations, the harmonic approximation (HA) is usually adopted, in which only the second derivative of the Born-Oppenheimer (BO) energy surface is considered assuming that atomic displacements are sufficiently small compared with interatomic distances. 4 The HA is in many cases valid and useful for obtaining phonon dispersion curves and discussing phase stability based on the vibrational free energy. 5 However, it fails to describe many important properties related to the lattice anharmonicity such as thermal expansion, LTC, and the temperature and volume dependences of phonon frequencies, for which we must go beyond the HA.
When the cubic and higher-order anharmonic terms of the BO energy surface are sufficiently small compared with the harmonic one, the anharmonic effects can be treated by the many-body perturbation theory (PT). 6 The PT has successfully been employed to explain phonon linewidths in semiconductors 7,8 and to understand/predict LTC values of a wide variety of materials. [9][10][11][12][13][14][15][16][17][18][19][20][21][22] In these calculations, the lowestorder contribution associated with the third-order anharmonic terms is considered to describe the intrinsic phonon-phonon scattering events, which is of primary importance in 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 helium 23 and superconducting sulfur hydride under pressure, 24,25 the zero-point motion is so significant that the anharmonic renormalization of phonon frequencies must be considered beyond the PT. Furthermore, the PT is not applicable to high-temperature phases of dielectric materials due to the existence of unstable phonon modes within the HA. This seriously limits the predictive power of first-principles lattice dynamics methods for emergent thermoelectric materials 26,27 and hybrid perovskite solar cells. 28,29 To overcome the aforementioned limitations, several DFTbased methods have recently been developed for treating anharmonic effects in solids nonperturbatively, which rely on the vibrational self-consistent-field theory, 30,31 the self-consistent phonon (SCP) theory, [32][33][34][35][36] or ab initio molecular dynamics (AIMD) methods. 37,38 They have been employed to study anharmonic renormalization of phonon frequencies as well as its effects on the phase stability, 37,39 the superconducting critical temperature, 25,[40][41][42][43] and the LTC of high-temperature phases. 35,44 For example, the self-consistent ab initio lattice dynamics (SCAILD) 33 and the stochastic self-consistent harmonic approximation (SSCHA) 34 are SCP-based approaches which can incorporate the effect of lattice anharmonicity at the mean-field level. In these methods, anharmonic phonon frequencies, or equivalently, effective harmonic force constants are obtained self-consistently by repeatedly calculating atomic forces in supercells with suitably chosen atomic configurations. These stochastic algorithms are quite useful since a calculation of fourth-order anharmonic force constants is unnecessary. More recently, another efficient implementation of the SCP theory was developed, 35 which employs anharmonic force constants calculated using the compressive sensing lattice dynamics method. 45 The temperature-dependent effective potential (TDEP) method 37,46 is an AIMD-based approach, in which effective harmonic and cubic force constants are extracted from displacement-force data sets calculated for atomic configurations sampled by AIMD. Although the AIMD-based approaches are not valid in the low-temperature region due to their inability to account for the zero-point motion, they should be accurate in the high-temperature region because anharmonic effects are fully included. As a result of these continuous efforts, an accurate first-principles 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 SCP-based approaches developed by the authors. 35 Using microscopic force constants estimated by supercell methods as inputs, 45,47,48 our method can calculate temperature-dependent phonon frequencies and LTC efficiently. We demonstrate the validity of our method by applying it to cubic SrTiO 3 , which is realized only at temperatures above 105 K. 49 In Sect. II, we present theoretical formulations and some numerical techniques that underlie our SCP approach. In Sect. 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 Sect. IV.

II. METHODOLOGY
A. Taylor expansion potential Let us start by expressing the potential energy of an interacting nuclear system U as a Taylor expansion with respect to the atomic displacement: . . . ; n κ n )u µ 1 ( 1 κ 1 ) · · · u µ n ( n κ n ). ( Here, u µ ( κ) is the displacement of atom κ in the th unit cell along the µ (= x, y, z) direction and Φ µ 1 ...µ n ( 1 κ 1 ; . . . ; n κ n ) is the nth-order derivative of the potential energy with respect to displacement, which is usually termed the interatomic force constant (IFC). In Eq. (2), the IFCs are determined at the ground state and therefore have no temperature dependence.
In the HA, all of the anharmonic terms U n (n > 2) are neglected. Then, the Hamiltonian of the system H 0 = T + U 2 can be represented as H 0 = 1 2 q j ω q j A q j A † q j with harmonic phonon frequency ω q j and associated displacement operator where b q j and b † q j are the annihilation and creation operators of a phonon with crystal momentum q and branch index j, respectively. The phonon frequency ω q j can be obtained from the dynamical matrix where M κ is the mass of atom κ and r( ) is the translational vector of the th unit cell. By diagonalizing the matrix D(q), one obtains the squared harmonic frequency and the polarization vector as D(q)e q j = ω 2 q j e q j . Since the harmonic IFCs are temperature-independent, the intrinsic temperature dependence of phonon frequencies and polarization vectors {ω q j , e q j } is absent within the HA. The atomic displacement u µ ( κ) can be represented in terms of A q j as follows: Here, N is the number of q points in the first Brillouin zone (BZ) and e µ (κ; q j) is a component of e q j . By substituting Eq. (4) into Eq. (2), we obtain Here and in the following, we use q for the shorthand notation of (q, j), satisfying q = (q, j) and −q = (−q, j). The term ∆(q) is 1 if q is a vector of the reciprocal lattice G and 0 otherwise. The coefficient V(q 1 ; . . . ; q n ) is defined as ..µ n (0κ 1 ; 2 κ 2 ; . . . ; n κ n )e i(q 2 ·r( 2 )+···+q n ·r( n )) . (7) Equation (5) is the nth-order potential energy represented in terms of the harmonic displacement operator A q .

B. Anharmonic perturbation theory
To describe the intrinsic phonon scattering processes and the 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 a perturbation H of the non-interacting Hamiltonian H 0 as Here, we omitted fifth-and higher-order terms since their contributions are much smaller than the cubic and quartic terms. Let G q (ω) be the one-phonon Green's function and G 0 q (ω) be that of the non-interacting system H 0 . Then, the following Dyson equation holds: with Σ q (ω) being the anharmonic self-energy, which can be estimated within a systematic diagrammatic approximation. If the self-energy correction is small and the condition ω q |Σ q (ω q )| is well satisfied, the phonon quasiparticle picture is still valid and the phonon frequency shift ∆ q and the linewidth Γ q are given as ∆ q = − 1 ReΣ q (ω q ) and Γ q = 1 ImΣ q (ω q ). Figure 1 shows the first-and second-order self-energy diagrams associated with the cubic and quartic anharmonic terms. The tadpole (a) and bubble (b) are the second-order diagrams resulting from the cubic terms and their explicit formulae are given as follows: 6 Here, n i = n(ω q i ) = 1/(e β ω q i − 1) is the Bose-Einstein distribution function and ω c = ω + i0 + with 0 + being a positive infinitesimal. In addition, the summation in Eq. (11) is restricted to the pairs (q 1 , q 2 ) satisfying the momentum conservation q 1 + q 2 = q + G. The tadpole diagram is real and therefore gives rise to only a frequency shift. In Eq. (10), we neglect the processes involving zone-center acoustic modes as the intermediate state 0 j 1 because of the singularity of the matrix element V(−q; q; 0 j 1 ) 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 quasi-harmonic approximation. 50 The bubble diagram provides the dominant contribution to the phonon linewidth Γ bubble q (ω q ), and the quasiparticle lifetime τ q can be estimated as τ q = 1/(2Γ bubble q (ω q )). The loop diagram (c) is the first-order contribution from the quartic terms, whose expression is Since Σ loop q is real, the loop diagram only gives a frequency shift. The expressions for the higher-order 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. [51][52][53] First-principles calculations of the bubble and tadpole diagrams have been performed to explain the temperature dependences of the Raman shift and linewidth of group IV and III-IV semiconductors, 7,8 . Recently, there have been many DFT-based estimations of the bubble diagram mostly for the thermal conductivity prediction described in Sect. II D. On the other hand, the calculation of the loop diagram has only been reported for relatively simple systems until very recently. [54][55][56] This is mainly because the calculation of Σ loop q relies on quartic IFCs, which are technically more difficult to obtain from DFT than harmonic and cubic IFCs. 43,57 In Sect. II E, we will introduce some recent technical developments for estimating the quartic terms.

C. Self-consistent phonon theory
When the anharmonic terms are comparable with the harmonic term, they need to be treated nonperturbatively. The SCP theory, originally developed by Hooton 32 and followed by more elegant formulations based on the many-body theory, 23 Here, H 0 is the effective harmonic Hamiltonian, which can be written as H 0 = 1 2 q Ω q A q A † q with the renormalized phonon frequency Ω q and the associated displacement operator A q . Then, the free energy of the original system F can be written as the cumulant expansion Here, F 0 = − 1 β log Z and X H 0 = Z −1 Tr(Xe −βH 0 ) with the partition function Z = Tr(e −βH 0 ), and X n c indicates the nthorder cumulant of the operator X. The first-order SCP theory only considers the first cumulant in Eq. (14), in which the following Gibbs-Bogoliubov inequality is satisfied: Since the variational principle holds in the first-order SCP theory, the solution can be found by minimizing the right-hand side of the above equation with respect to the adjustable parameters such as anharmonic frequencies {Ω q }, 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 ∂F 1 /∂Ω q = 0 with F 1 = F 0 + H − H 0 H 0 , the following SCP equation can be readily derived: By solving Eqs. (16) and (17) self-consistently, anharmonic phonon frequencies {Ω q } can be obtained. The equation is valid even if some of the harmonic phonon frequencies are unstable (ω 2 q < 0). These unstable modes are renormalized by the second term 2Ω q I q to give a stable phonon mode (Ω 2 q ≥ 0). Since the formulation assumes the existence of a welldefined Hamiltonian H 0 , an imaginary frequency (Ω 2 q < 0) is not allowed as a solution of the equation.
The same result has also been obtained on the basis of the Green's function approach. 35,60,61 Diagrammatically, the SCP theory corresponds to the inclusion of the infinite set of diagrams that can be generated from the loop diagram shown in Fig. 2(a). Therefore, the SCP theory automatically includes the contribution from the higher-order diagram shown in Fig. 1(e). Substituting ω q for Ω q 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 assumed that the polarization vectors are changed by anharmonic effects. Such an assumption is reasonable for simple systems containing a few atoms in the primitive cell 40 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 extended the SCP equation to include the off-diagonal elements of the loop self-energy and developed an efficient algorithm for solving the equation. 35 We applied the extended method to cubic SrTiO 3 and demonstrated the importance of the PM for the soft modes at high-symmetry q points. 35 The same implementation was also employed to study anharmonic effects in the superconducting hydrogen sulfides. 42 Moreover, it is possible to extend the theory to include the tadpole diagram. 60 The deterministic implementation based on the reciprocal formalism [Eqs. (16) and (17)] is efficient and applicable to relatively complex structures. Also, we can easily evaluate the finite-size effect by increasing the q 1 mesh in Eq. (17) based on the Fourier interpolation, which can be significant near a structural phase transition. 35,62 However, as already mentioned in Sect. II B, first-principles calculation of the quartic IFCs necessary for I q has been a technical challenge. To avoid the cumbersome estimation of the quartic IFCs, several stochastic implementations have been developed 33,34,36,43 and applied to study anharmonic effects in simple metals, 33,63 light-element superconductors, 25,41 transition metal dichalcogenides, 64 and simple perovskite materials. 44 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) and (17)] because they do not require the expansion of Eq. (1) nor the truncation of higher-order terms carried out in Eq. (8). However, since a large number of DFT calculations is generally needed to obtain convergence, the stochastic approaches are more computationally intensive than the deterministic one.

D. Lattice thermal conductivity
The LTC has been intensively studied as it plays an important role in optimizing the thermoelectric figure of merit ZT . Most of the theoretical predictions of LTC rely on either the Boltzmann transport equation (BTE) or the MD method. According to the formulation of BTE for phonons developed by Peierls, 65 the LTC is given as where V is the unit cell volume and v q = ∂ω q ∂q is the group velocity. The vector f q is a solution of the following linear equation: Here, Λ q q and Λ q qq (Λ q q q ) are the intrinsic transition rates of phonon q involving two and three phonons, respectively. For example, Λ q q includes the effect of phonon-isotope scattering and Λ q qq (Λ q q q ) is dominated by three-phonon scattering processes associated with the cubic lattice anharmonicity. To obtain the LTC of bulk solids, Eqs. (18) and (19) must be solved on a fine q 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. 9 Owing to the increased computational power and the development of efficient implementations, it is now feasible to solve Eq. (19) for relatively complex materials using iterative algorithms, 9,66,67 or even by using the direct method for simple systems. 20,68 To make the computation more feasible, the relaxationtime approximation (RTA) is often employed in the literature, whereby the LTC is simplified to In the RTA, the current vertex corrections are neglected and the transport relaxation time is approximated by the quasiparticle lifetime τ q . 69 Therefore, the RTA incorrectly treats the normal process of three-phonon scattering as resistive, which results in κ RTA < κ. The underestimation by the RTA is predominant in high-LTC materials such as diamond 70 as well as in a variety of two-dimensional materials. [10][11][12][13] Except for these cases, the RTA usually yields results similar to the full solution to the BTE. Hence, the RTA has successfully been employed to analyze and predict the LTC of various kinds of solids. [14][15][16][17][18][19][20][21][22] While the predictive accuracy of the conventional BTE [Eq. (18)] and BTE-RTA [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 from estimating the LTC of hightemperature phases because of the imaginary modes within the HA. To overcome this limitation, we made the following extension: 35 where Ω q is the SCP frequency,ṽ q = ∂Ω q ∂q is the group velocity of the renormalized phonon, andñ q = n(Ω q ). The lifetime of the renormalized phononτ q 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 {ω q , e q } are replaced with the SCP solution {Ω q , q }. We have successfully applied the new method to predict the LTC of cubic SrTiO 3 including its unusual temperature dependence, which will be discussed in detail in Sect. III. Other previous studies employed a similar approach for studying phonon transport properties in Bi 2 Te 3 , 71 PbTe, 72 and SrTiO 3 , 73 in which the effective harmonic and cubic IFCs were obtained by the temperature-dependent effective potential method. 37,46 These studies also indicated the importance of the renormalization of IFCs in the LTC of the studied materials.
Incorporating the four-phonon scattering process Λ q q q q into the BTE [Eq. (19)] is still a major computational challenge even within the RTA. In order to achieve this, the higherorder diagram shown in Fig. 1(d) must be calculated with a dense q mesh, which has only been reported for the group IV semiconductors based on empirical potentials. 74 More work is thus needed to develop a robust understanding of the fourphonon scattering effects on the LTC of various kinds of solids including low-dimensional systems. MD is another powerful tool for calculating the LTC of solids. Since the anharmonicity is fully included, the MDbased methods, such as the non-equilibrium MD 75,76 and Green-Kubo methods, 77 are in principle more accurate than the 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 the reliable prediction of the LTC with MD, it is necessary to employ a reasonably large supercell to sample long-wavelength phonons that make dominant contributions to the LTC. 48,78,79 In addition, the simulation length must be sufficiently long to reduce statistical errors. These requirements make the MD methods very costly, especially for high-LTC materials due to large mean free paths and long relaxation times of phonons. Therefore, straightforward applications of AIMD have only been attempted for a few materials having low LTC. 80,81 To make the computation feasible, classical force fields are commonly employed, but they can be a source of inaccuracy. To address this problem, much effort has been made to improve the accuracy of classical force fields or to develop a new force field by using DFT results as training data. [82][83][84][85]

E. First-principles 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: density-functional perturbation theory (DFPT) and the supercell method.
DFPT is one of the standard ab initio methods for calculating harmonic IFCs. 5,86,87 In DFPT, harmonic dynamical matrix and linear electron-phonon coupling coefficients for an arbitrary q-point can be obtained. Recently, an efficient implementation based on DFPT has also been developed for cubic IFCs. 88 While DFPT is very accurate and efficient since it does not involve expensive supercell calculations, its implementation is relatively complicated and its extension to quartic and higher-order 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 or inversely as 89 which are valid when the displacement u µ ( κ) is small, say u µ ( κ) ∼ 0.01 Å. F ν [ κ ; u µ ( κ)] is the atomic force acting on the κ th atom in the th cell along direction ν when the displacement of u µ ( κ) is made with all other atoms kept at the equilibrium positions. Equation (23) is defined for every single coefficient Φ µν ( κ; κ ), many of which are linearly dependent on each other because of the space group symmetry, the permutation symmetry, and the lattice periodicity. 47,48 For the convenience of the following discussion, let us introduce the column vectors u and F, which consist of 3s atomic displacements and forces in the supercell containing s atoms, respectively. We also denote the ith displacement pattern as u i and the associated force vector as F i . Then, denoting a column vector comprising n unique harmonic coefficients by Φ 2 , the set of linear equations [Eq. (23)] can be symbolically written as Here, A 2 is a 3sm × n matrix defined as , and m is the number of displacement patterns. Note that Eq. (24) is valid if the displacements are sufficiently small to neglect anharmonic effects in F . If we consider sufficient displacement patterns and make the matrix A 2 full rank, the harmonic IFCs can be obtained, in the leastsquares sense, as 89,90 where A + 2 is the pseudoinverse of the matrix A 2 , which can be readily obtained by singular value decomposition. A similar equation can be formed for the cubic terms as 90 in which ∆F 3 = F − A 2 Φ 2 is the cubic contribution to the atomic forces and A T 3 = [A T 3 (u 1 ), . . . , A T 3 (u m )] with A 3 (u) defined as A 3 (u) = −∂ 2 U 3 /∂Φ T 3 ∂u. Again, Eq. (26) assumes that the quartic and higher-order contributions to atomic forces F are negligible, which is valid when the displacement is small. To make the matrix A 3 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 shortrange, which is usually the case in covalent compounds. By contrast, long-range cubic interactions may be significant in ionic compounds, 91 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 linear equations for quartic and higher-order IFCs, but they are difficult to solve because of a large number of DFT calculations necessary to make the matrix A n (n ≥ 4) full rank. Moreover, an accurate DFT calculation of ∆F n is difficult because each ∆F n is very small in the small-displacement limit. While this numerical issue may be avoided by using a larger displacement, finding an appropriate value for the displacement is formidable as it depends on the atomic environment.
Performing AIMD simulation is a reasonable solution to this problem as it can generate physically relevant atomic configurations automatically. From the displacement-force data set obtained from the AIMD trajectory, we may estimate IFCs by solving the ordinary least-squares (OLS) problem 48 Here, Φ is the column vector defined as Given a sufficient number of AIMD snapshots, the unknown IFCs can be estimated asΦ = A + F . Using this approach, we successfully extracted anharmonic IFCs of Si and Mg 2 Si from AIMD trajectories. 48 Recently, Zhou et al. proposed a more sophisticated approach called the compressed sensing lattice dynamics (CSLD) method, 45 where they estimated IFCs using the least absolute shrinkage and selection operator (LASSO) 92 : Owing to the 1 -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 noise 93 and a penalty term helps to avoid the overfitting inherent in the OLS method. The coefficient λ is a hyperparameter that controls the tradeoff between the sparsity and accuracy of the model; all IFCs will be zero in the large-λ limit and λ = 0 corresponds to the OLS regression. An appropriate value of λ can be estimated, for instance, by cross-validation, as will be demonstrated in Sect. III. In Ref. 45, Zhou et al. extracted harmonic and anharmonic IFCs of Si, NaCl, and Cu 12 Sb 4 S 13 using Eq. (28) combined with AIMD simulation. On the basis of the IFCs estimated by LASSO, they calculated the LTC of these materials and obtained good agreement with experimental data. Using the same approach, we successfully extracted anharmonic IFCs of SrTiO 3 , 35 superconducting hydrogen sulfides, 42 and thermoelectric clathrates and applied them to SCP and BTE calculations. Also, compressed sensing has successfully been employed to develop interatomic potentials 94 and a cluster expansion model 95 based on DFT.

III. APPLICATIONS TO STRONTIUM TITANATE
In Ref. 35, we calculated anharmonic phonon properties of cubic SrTiO 3 (c-STO), which is stable above 105 K. All of the DFT calculations were performed using the Vienna Ab initio Simulation Package (vasp), 96,97 where the generalized gradient approximation for solids 98 (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. 99,100 Using quartic IFCs estimated by LASSO, we performed SCP simulations and obtained overall good agreement 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 Heyd-Scuseria-Ernzerhof hybrid functional (HSE06, hereafter HSE) 101 . Before performing phonon calculations, we optimized the lattice constant of c-STO using the HSE functional and obtained 3.900 Å, in good agreements with the experimental value of 3.905 Å (Ref. 102, 293 K) and the previous HSE result of 3.904 Å. 99 The computational conditions employed in this study were basically the same as those in the previous PBEsol study, 35 except that BZ integration was performed with the 8×8×8 Monkhorst-Pack k-point grid. All of the phonon calculations reported here were performed with a 2×2×2 cubic supercell (40 atoms) using the alamode package, 48,103 an open-source software developed for modeling lattice anharmonicity and thermal conductivity.

A. 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 2×2×2 supercell were considered. For quartic IFCs, we only considered onsite, two-body, and three-body IFCs within a cutoff length of 6 Å and neglected less important four-body terms. Also, only onsite and two-body terms within the same cutoff radius were considered for the fifth-and sixth-order 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. 47,48 The numbers of independent parameters were found to be 698, 2215, 43, and 125 for cubic, quartic, quintic, and sextic 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 about 5,000 DFT calculations. To determined these IFCs by LASSO, we prepared 40 displacement-force data sets using AIMD, whose detailed procedure is described in Ref. 35. After that, the data set was split into four smaller subsets and the hyperparameter λ in Eq. (28) was selected by four-fold cross-validation. 104 The solution path was obtained by using the coordinate descent method, 105 for which each column vector of the matrix A was standardized beforehand. Upon solving Eq. (28), harmonic IFCs were fixed to the values obtained by the OLS method. Figure 3 shows the result of the cross-validation obtained for PBEsol. Since very similar results were obtained for HSE, they are not shown here. The top panel shows the relative error of the atomic forces, which is defined as the square root of AΦ − F 2 2 / F 2 2 . 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 monotonically decreases, whereas the vali- This value of λ was selected as an optimal choice because it is expected to give the best prediction accuracy for 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 number of elements. The accuracy of the extracted IFCsΦ was carefully verified by applying them to the independent atomic configurations sampled by AIMD.

B. 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 off-diagonal elements of the loop diagram as in Ref. 35. The convergence of the anharmonic frequencies {Ω q } at the commensurate 2×2×2 qpoint grid was carefully checked with respect to the intermediate q 1 grid in Eq. (17) and found sufficiently converged results with 8×8×8 q 1 points. In addition, the non-analytic correction to the harmonic dynamical matrix was included to reproduce the LO-TO splitting around the Γ point. To this end, we employed the Born effective charges and the dielectric constant reported in Ref. 35 for PBEsol. Calculating these quantities using the HSE functional was very expensive. Therefore, as an approximation, we employed the Born effective charges obtained by PBEsol and the experimental dielectric constant of 5.18 108 in the following HSE calculations. measured at 300 K, 106,107 the theoretical curves at 300 K are highlighted as thick solid lines. In the harmonic phonon dispersion, unstable phonon modes (ω 2 q < 0) occur at points Γ (0, 0, 0) and R ( 1 2 , 1 2 , 1 2 ), which correspond to the ferroelectric (FE) and antiferrodistortive (AFD) modes, respectively. The harmonic phonon frequency of the FE mode is 101i cm −1 with HSE, whose absolute value is larger than the PBEsol result of 58i cm −1 , 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 76i cm −1 (PBEsol) to 25i cm −1 (HSE) as summarized in Table I. These tendencies are consistent with the previous DFT study of Wahl et al. 99 As can be seen in Fig. 4, the quartic anharmonicity generally increases the vibrational frequencies in c-STO, which is particularly predominant in the low-energy optical modes. The PBEsol result obtained here is almost identical to our previous result in Ref. 35. The SCP frequency of the FE mode Ω FE is 133 cm −1 at 300 K with PBEsol, which overestimates the INS data of 87.1±5.6 cm −1 (Ref. 49, 293 K) and the IR data of 89 cm −1 (Ref. 109, 300 K). By contrast, the frequency of the AFD mode at point R is underestimated with PBEsol. We have found that the underestimation at point R can be reduced by employing HSE as shown in the right panel of Fig. 4. However, Ω FE was improved only slightly by HSE and the SCP frequency of 128 cm −1 still overestimates the experimental values. In what follows, we discuss the origin of the disagreement by going beyond the SCP theory.

C. 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 reveal its effect quantitatively, we calculated the bubble self-energy shown in Fig. 2(b) and estimated the frequency shift of the FE mode from ∆ bubble q = − 1 ReΣ bubble q (Ω q ). The results of Ω q + ∆ bubble q are shown in Table I as "SCP+B" (B stands for bubble). The correction ∆ bubble q 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 −13 cm −1 with PBEsol and −19 cm −1 with HSE. For the AFD mode, the correction was very large, ∆ bubble q ≈ −29 cm −1 with PBEsol, and the corrected frequency Ω q + ∆ bubble q was almost zero. It is interesting to note that the frequency dependence of ReΣ bubble q (ω) was marginal in the low-frequency region (ω < 200 cm −1 ) and the difference between ∆ bubble q and − 1 ReΣ bubble q (0) was minor, especially when the SCP frequency was small. Therefore, it may be reasonable to make the static approximation ofΣ bubble q (ω) ≈Σ bubble q (0) in order to calculate the frequency shift due to the bubble diagram, as carried out in Ref. 110. With the correction from the bubble diagram, the theoretical values based on HSE are in better agreement with the experimental values. Although other theoretical and numerical aspects must be carefully investigated, our results indicate the importance of accurate theoretical and computational modeling of anharmonic effects for robust understanding of the lattice dynamics in severely anharmonic materials.
To make a straightforward comparison with the experimental data, we also calculated the spectral function A q (ω) = − 1 π ImG q (ω) using the following formula: Since A q (ω) has a direct connection with the dynamical structural factor measured in INS experiments, it should be more reasonable and accurate to estimate theoretical phonon fre- quencies from peak positions of Eq. (29). The estimated frequencies are also shown in Table I, which differ from the "SCP+B" results. This indicates the limited accuracy of the SCP+B results when the condition of Ω q |Σ q | is not well satisfied. Figure 5 shows the spectral function at 300 K calculated along the high-symmetry lines of the BZ. For comparison, we also show the SCP dispersion curves by dashed lines. It is evident from the figure that the AFD mode and high-energy optical modes above ∼ 400 cm −1 are strongly damped by the cubic anharmonicity. Moreover, the frequency shift caused by the bubble diagram is most significant in the low-lying optical modes around point Γ and along path R-M.

D. Anomalous thermal transport in SrTiO 3
Finally, we calculated the LTC of c-STO using Eq. (21) with 20×20×20 q-grid points. In Fig. 6, we compare our computational results with experimental data from Refs. 111 and 112. As can be seen in the figure, we obtained overall good agreement with the 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 is because the 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 of the LTC (κ ∝ T −α , α ∼ 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 κ ∝ T −1 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. 36, 72, and 80. We think the effect of the frequency renormalization by the quartic terms can be significant in other ultralow-LTC materials such as SnSe and clathrate, which will be a topic of future work.
It is interesting to note that in c-STO more than 70% of the total thermal conductivity is carried by high-frequency optical modes above 150 cm −1 , 35,73 as can be seen in the inset of Fig. 6. Here, κ(ω) is the thermal conductivity spectrum defined as κ(ω) = q κ q δ(ω − ω q ) with κ q being the contribution from each phonon mode q. 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 first-principles framework to study lattice dynamical properties of strongly anharmonic crystals. By combining an efficient implementation based on the self-consistent 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 lattice-dynamics wavefunctions, which is essential for predicting the lattice thermal conductivity (LTC) of anharmonic crystals from first principles.
We demonstrated the high predictive accuracy of the developed computational approach by carefully comparing softmode frequencies and LTC values of cubic SrTiO 3 with available experimental data, for which overall good agreement was obtained. In addition, we showed that accurate, albeit expensive, calculation based on the HSE hybrid functional was feasible within the present framework, where a marked 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 provide more insight into anharmonic effects in severely anharmonic materials, such as ferroelectric, thermoelectric, and photovoltaic materials, as well as in light-element superconductors where the conventional harmonic approximation and the perturbative treatments break down.