First-Principles Lattice Dynamics Method for Strongly Anharmonic Crystals

First-Principles Lattice Dynamics Method for Strongly Anharmonic Crystals

Terumasa Tadano International Center for Young Scientists (ICYS), National Institute for Materials Science, Tsukuba 305-0047, Japan Research and Services Division of Materials Data and Integrated System (MaDIS), National Institute for Materials Science, Tsukuba 305-0047, Japan    Shinji Tsuneyuki Department of Physics, The University of Tokyo, Tokyo 113-0033, Japan Institute for Solid State Physics, The University of Tokyo, Kashiwa 277-8581, Japan

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 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 phonon-mediated superconductors.Ziman (1960) With the advent of the increased computational power and developments of computational methods based on density-functional theory (DFT),Hohenberg and Kohn (1964); Kohn and Sham (1965) 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 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 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 volume-dependence 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). 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 lowest-order contribution associated with the third-order anharmonic terms is considered to describe the intrinsic phonon-phonon 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 zero-point 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 high-temperature phases of dielectric materials due to the existence of unstable phonon modes within the HA. This seriously limits the predictive power of the first-principles 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 DFT-based methods have recently been developed for treating anharmonic effects in solids nonperturbatively, which rely on the vibrational self-consistent-field theory,Bowman (1978); Monserrat et al. (2013) the self-consistent 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 high-temperature phases.Tadano and Tsuneyuki (2015); van Roekeghem et al. (2016b) For example, the self-consistent ab initio lattice dynamics (SCAILD)Souvatzis et al. (2008) and the stochastic self-consistent harmonic approximation (SSCHA)Errea et al. (2014) 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,Tadano and Tsuneyuki (2015) which employs anharmonic force constants calculated using the compressive sensing lattice dynamics method.Zhou et al. (2014) The temperature-dependent effective potential (TDEP) methodHellman et al. (2011); Hellman and Abrikosov (2013) 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. With these constant efforts, accurate first-principles modeling of lattice vibration is becoming practical even for severely anharmonic materials.

In this paper, we review the recent developments of first-principles lattice dynamics methods for strongly anharmonic materials, particularly focusing on one of the SCP-based 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 temperature-dependent 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 exchange-correlation 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:


Here, is the displacement of atom in the th unit cell along the direction and is the th-order 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


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:


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


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


Equation (5) is the th-order 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 non-interacting Hamiltonian as


Here, we omitted fifth- and higher-order terms since their contributions are much smaller than the cubic and quartic terms. Let be the one-phonon Green’s function and be that of the non-interacting system . Then, the following Dyson equation holds:


with being the anharmonic self-energy which can be estimated within a systematic diagrammatic approximation. If the self-energy correction is small and the condition is well-satisfied, the phonon quasiparticle picture is still valid and the phonon frequency shift and the linewidth are given as and .

Figure 1: Feynman diagrams of anharmonic self-energies up to the second-order. The solid lines represent the free phonon propagator and the small open circles are vertices involving three or four phonons. The diagrams (a), (b), and (c) are termed tadpole, bubble, and loop diagrams, respectively. The diagrams (a), (c), and (e) are frequency-independent and real, whereas diagrams (b) and (d) are complex and have frequency dependence.

Figure 1 shows the first- and second-order self-energy diagrams associated with the cubic and quartic anharmonic terms. The tadpole (a) and the bubble (b) are the second-order diagrams resulting from the cubic terms and their explicit formulae are given as follows:Maradudin and Fein (1962)


Here, is the Bose-Einstein 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 zone-center 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 quasi-harmonic 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 first-order contribution from the quartic terms, whose expression is


Since is real, the loop diagram solely causes the frequency shift. The expression 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.Tripathi and Pathak (1974); Monga and Pathak (1978); Procacci et al. (1992)

First-principles calculations of the bubble and tadpole diagrams have been performed for explaining the temperature dependence of Raman shift and linewidth of group IV and III-IV semiconductors,Narasimhan and Vanderbilt (1991); Debernardi et al. (1995). Recently, there are many DFT-based 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 Self-consistent phonon theory

When the anharmonic terms are comparable with the harmonic term, they need to be treated nonperturbatively. The self-consistent phonon (SCP) theory originally developed by Hooton,Hooton (1958) followed by more elegant formulations based on the many-body 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


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


Here, and with the partition function , and indicates the th-order cumulant of the operator . 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 , 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:


By solving Eqs. (16) and (17) self-consistently, 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 well-defined Hamiltonian , imaginary frequency () is not allowed as a solution of the equation.

Figure 2: (a) Diagrammatic representation of the SCP propagator where the loop diagram is considered. The double line represents the SCP propagator and the solid line is that of free phonon. (b) Bubble diagram defined based on the SCP propagator.

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 higher-order 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 off-diagonal elements of the loop self-energy 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 high-symmetry 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, first-principles 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 higher-order 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 figure-of-merit . 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


where is the unit cell volume and is the group velocity. The vector is a solution of the following linear equation:


Here, and are intrinsic transition rates of phonon involving two and three phonons, respectively. For example, includes the effect of phonon-isotope scatterings and is dominated by three-phonon 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 relaxation-time approximation (RTA) is often employed in the literature, whereby LTC is simplified as


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 three-phonon scattering as resistive, which results in . The underestimation by RTA is predominant in high-LTC materials such as diamond Ward et al. (2009) as well as in a variety of two-dimensional 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 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 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)


where is the SCP frequency, is the group velocity of the renormalized phonon, and . The lifetime of the renormalized phonon associated with the three-phonon 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 temperature-dependent 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 four-phonon scattering process into BTE [Eq. (19)] is still a big computational challenge even within RTA. In order to achieve this, the higher-order 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 four-phonon scattering effects on LTC of various kinds of solids including low-dimensional systems.

MD is another powerful tool for calculating LTC of solids. Since the anharmonicity is fully included, the MD-based methods, such as the non-equilibrium MDEvans (1982); Müller-Plathe (1997) and the Green-Kubo 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 long-wavelength 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 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.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 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 method for calculating harmonic IFCs.Giannozzi et al. (1991); Gonze (1997); Baroni et al. (2001) In DFPT, harmonic dynamical matrix and linear electron-phonon 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 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 first-order numerical derivative of the atomic forces as


or inversely asParlinski et al. (1997)


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


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 least-squares sense, asParlinski et al. (1997); Chaput et al. (2011)


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)


in which is the cubic contribution to the atomic forces and with defined as . Again, Eq. (26) assumes the quartic and higher-order 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 short-ranged, which is usually the case in covalent compounds. By contrast, long-range 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 higher-order 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 displacement-force data set obtained from the AIMD trajectory, we may estimate IFCs by solving the ordinary least-squares (OLS) problemTadano et al. (2014)


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):


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 trade-off 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 cross-validation 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 (c-STO), 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 Heyd-Scuseria-Ernzerhof hybrid functional (HSE06, hereafter HSE)Heyd et al. (2003); *Heyd:2006dc. 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. 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 Monkhorst-Pack -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, two-body, and three-body IFCs inside a cutoff length of 6 Å and neglected less important four-body terms. Also, only onsite and two-body terms inside 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.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 displacement-force 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 four-fold cross-validation.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: (Color online) (a) Relative errors in atomic forces and (b) the number of nonzero coefficients as a function of the hyperparameter obtained for PBEsol. The dotted vertical line indicates the selected value of that minimizes the cross-validation score.

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 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

Figure 4: (Color online) Temperature-dependent anharmonic phonon dispersion of cubic SrTiO calculated with the PBEsol (left) and HSE (right) exchange-correlation functionals. The gray solid lines in the left (right) figure represent the SCP solutions at different temperatures ranging from 200 K (0 K) to 1000 K. The dotted lines are harmonic lattice dynamics results and the red open circles are experimental values at 300 K adapted from Refs. Stirling, 1972 and Cowley et al., 1969. Theoretical curves at 300 K are highlighted as thick solid lines for comparison with experimental data.

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 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 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. 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 c-STO 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 c-STO, which is particularly predominant in the low-energy 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)

Table 1: Phonon frequencies (cm) of the (FE) and (AFD) soft modes in c-STO calculated with different levels of approximation. Calculated anharmonic frequencies are evaluated at 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 self-energy 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 low-frequency 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:


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 high-symmetry 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 high-energy optical modes above cm are strongly damped by the cubic anharmonicity. Moreover, the frequency shift by the bubble diagram is most significant in the low-lying optical modes around the point and along the R-M path.

Figure 5: (Color online) Phonon spectral function of cubic SrTiO at 300 K calculated with the HSE functional. The SCP solution at 300 K is also shown by white dashed lines for comparison.

iii.4 Anomalous thermal transport in SrTiO

Finally, we calculated LTC of c-STO 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.

Figure 6: (Color online) Temperature dependence of calculated LTC in cubic SrTiO compared with experimental values adapted from Refs. Muta et al., 2005; Popuri et al., 2014. The dashed line shows LTC values calculated without changing the SCP frequencies from those at 200 K, which approximately follow . (Inset) Thermal conductivity spectra at 300 K and 900 K obtained with the HSE functional.

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 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,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 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 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 light-element superconductors where the conventional harmonic approximation and the perturbative treatments break down.


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.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description