First-Principles Lattice Dynamics Method for Strongly Anharmonic Crystals

TL;DR

This paper introduces a nonperturbative first-principles lattice dynamics method based on self-consistent phonon theory, enabling accurate analysis of strongly anharmonic crystals where traditional methods fail.

Contribution

The paper presents a novel computational approach that efficiently incorporates anharmonic effects nonperturbatively using recent numerical techniques and self-consistent phonon theory.

Findings

Accurately predicts phonon frequencies in cubic strontium titanate

Successfully calculates lattice thermal conductivity in strongly anharmonic materials

Demonstrates compatibility with hybrid exchange-correlation functionals

Abstract

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.

Tables

Table 1. Table 1: Phonon frequencies (cm -1 ) of the Γ 15 subscript Γ 15 \Gamma_{15} (FE) and R 25 subscript 𝑅 25 R_{25} (AFD) soft modes in c-STO calculated with different levels of approximation. Calculated anharmonic frequencies are evaluated at 300 K.

	PBEsol					HSE				Expt.
	HA	SCP	SCP + B	Peak in $A_q\left(\omega \right)$		HA	SCP	SCP + B	Peak in $A_q\left(\omega \right)$
${\Gamma }_{15}$ (FE)	58$i$ a)	133	120	120		101$i$	128	109	104	87.1$\pm$5.6 b),  89 c)
$R_{25}$ (AFD)	76$i$ a)	31	2	18		25$i$	66	53	52	48.7$\pm$1.8 d)

Equations

$$U-U_0=U_2+U_3+U_4+\dots ,$$

$$U_n=\frac{1}{n!}\sum _{\{\ell ,\kappa ,\mu \}}{\Phi }_{{\mu }_1\dots {\mu }_n}({\ell }_1{\kappa }_1;\dots ;{\ell }_n{\kappa }_n)u_{{\mu }_1}({\ell }_1{\kappa }_1)\cdots u_{{\mu }_n}({\ell }_n{\kappa }_n).$$

$$D_{\mu \nu }(\kappa {\kappa }^{\prime };\mathbf{q})=\frac{1}{\sqrt{M_{\kappa }{M}_{{\kappa }^{\prime }}}}\sum _{{\ell }^{\prime }}{\Phi }_{\mu \nu }(0\kappa ;{\ell }^{\prime }{\kappa }^{\prime })e^{i\mathbf{q}\cdot \mathbf{r}({\ell }^{\prime })},$$

$$u_{\mu }(\ell \kappa )=\frac{1}{\sqrt{N{M}_{\kappa }}}\sum _{\mathbf{q}j}\sqrt{\frac{\hslash }{2{\omega }_{\mathbf{q}j}}}{A}_{\mathbf{q}j}{e}_{\mu }(\kappa ;\mathbf{q}j)\exp [i\mathbf{q}\cdot \mathbf{r}(\ell )].$$

$$U_n=\frac{1}{n!}\sum _{\{\mathbf{q},j\}}\Delta ({\mathbf{q}}_1+\cdots +{\mathbf{q}}_n)V(q_1;\dots ;q_n)A_{q_1}\cdots A_{q_n}.$$

$$\displaystyle V(q_{1};\dots;q_{n})=\bigg{(}\frac{\hbar}{2}\bigg{)}^{\frac{n}{2}}\frac{\Phi(q_{1};\dots;q_{n})}{\sqrt{\omega_{q_{1}}\cdots\omega_{q_{n}}}},$$

$$\Phi (q_1;\dots ;q_n)$$

$$=N^{1-\frac{n}{2}}\sum _{\{\kappa ,\mu \}}(M_{{\kappa }_1}\cdots M_{{\kappa }_n}{)}^{-\frac{1}{2}}{e}_{{\mu }_1}({\kappa }_1;q_1)\cdots e_{{\mu }_n}({\kappa }_n;q_n)$$

$$\times \sum _{{\ell }_2,\dots ,{\ell }_n}{\Phi }_{{\mu }_1\dots {\mu }_n}(0{\kappa }_1;{\ell }_2{\kappa }_2;\dots ;{\ell }_n{\kappa }_n)e^{i({\mathbf{q}}_2\cdot \mathbf{r}({\ell }_2)+\cdots +{\mathbf{q}}_n\cdot \mathbf{r}({\ell }_n))}.$$

$$H=H_0+H^{\prime }\approx H_0+U_3+U_4.$$

$$[{\mathbf{G}}_q(\omega ){]}^{-1}=[{\mathbf{G}}_q^0(\omega ){]}^{-1}-{\mathbf{\Sigma }}_q(\omega ),$$

$${\Sigma }_q^{\mathrm{tadpole}}=-\frac{1}{\hslash }\sum _{q_2,j_1=\mathrm{TO}}V(-q;q;0j_1)V(0j_1;q_2;-q_2)\frac{2n_2+1}{{\omega }_{0j_1}},$$

$${\Sigma }_q^{\mathrm{bubble}}(\omega )=\frac{1}{2\hslash }\text{sideset}{}^{\prime }{\sum }_{q_1,q_2,s=\pm 1}|V(-q;q_1;q_2){|}^2$$

$$\times \left[\frac{(n_1+n_2+1)}{s{\omega }_c+{\omega }_{q_1}+{\omega }_{q_2}}-\frac{(n_1-n_2)}{s{\omega }_c+{\omega }_{q_1}-{\omega }_{q_2}}\right].$$

$${\Sigma }_q^{\mathrm{loop}}=-\sum _{q_1}V(q;-q;q_1;-q_1)\frac{2n_1+1}{2}.$$

$$H={\mathcal{H}}_0+(H_0-{\mathcal{H}}_0+U_3+U_4)={\mathcal{H}}_0+{\mathcal{H}}^{\prime }.$$

$$F={\mathcal{F}}_0+{⟨{\mathcal{H}}^{\prime }⟩}_{{\mathcal{H}}_0,\mathrm{c}}+\frac{1}{2}{⟨({\mathcal{H}}^{\prime }{)}^2⟩}_{{\mathcal{H}}_0,\mathrm{c}}+\cdots .$$

$$F\le {\mathcal{F}}_0+{⟨{\mathcal{H}}^{\prime }⟩}_{{\mathcal{H}}_0,\mathrm{c}}={\mathcal{F}}_0+{⟨H-{\mathcal{H}}_0⟩}_{{\mathcal{H}}_0}.$$

$${\Omega }_q^2={\omega }_q^2+2{\Omega }_q{I}_q,$$

$$I_q=\sum _{q_1}\frac{\hslash \Phi (q;-q;q_1;-q_1)}{4{\Omega }_q{\Omega }_{q_1}}\frac{[2n({\Omega }_{q_1})+1]}{2}.$$

$$\kappa =\frac{\hslash }{NV{k}_{\mathrm{B}}T}\sum _q{\omega }_q{\mathbf{v}}_q\otimes {\mathbf{f}}_q{n}_q(n_q+1),$$

$${\beta }^{-1}{\mathbf{v}}_q\left(\frac{\partial n_q}{\partial T}\right)=\sum _{q^{\prime }}({\mathbf{f}}_q-{\mathbf{f}}_{q^{\prime }}){\Lambda }_q^{q^{\prime }}$$

$$\displaystyle+\sum_{q^{\prime},q^{\prime\prime}}\bigg{[}(\bm{f}_{q}+\bm{f}_{q^{\prime}}-\bm{f}_{q^{\prime\prime}})\Lambda_{qq^{\prime}}^{q^{\prime\prime}}+\frac{1}{2}(\bm{f}_{q}-\bm{f}_{q^{\prime}}-\bm{f}_{q^{\prime\prime}})\Lambda_{q}^{q^{\prime}q^{\prime\prime}}\bigg{]}.$$

$${\kappa }_{\mathrm{RTA}}=\frac{{\hslash }^2}{NV{k}_{\mathrm{B}}{T}^2}\sum _q{\omega }_q^2{\mathbf{v}}_q\otimes {\mathbf{v}}_q{n}_q(n_q+1){\tau }_q.$$

$${\tilde{\kappa }}_{\mathrm{RTA}}=\frac{{\hslash }^2}{NV{k}_{\mathrm{B}}{T}^2}\sum _q{\Omega }_q^2{\tilde{\mathbf{v}}}_q\otimes {\tilde{\mathbf{v}}}_q{\tilde{n}}_q({\tilde{n}}_q+1){\tilde{\tau }}_q,$$

$${\Phi }_{\mu \nu }(\ell \kappa ;{\ell }^{\prime }{\kappa }^{\prime })=\frac{{\partial }^{2}U}{\partial u_{\mu }(\ell \kappa )\partial u_{\nu }({\ell }^{\prime }{\kappa }^{\prime })}\approx -\frac{F_{\nu }[{\ell }^{\prime }{\kappa }^{\prime };u_{\mu }(\ell \kappa )]}{u_{\mu }(\ell \kappa )},$$

$$F_{\nu }[{\ell }^{\prime }{\kappa }^{\prime };u_{\mu }(\ell \kappa )]=-{\Phi }_{\mu \nu }(\ell \kappa ;{\ell }^{\prime }{\kappa }^{\prime })u_{\mu }(\ell \kappa ),$$

$${\mathbb{A}}_2{\mathbf{\Phi }}_2=\mathcal{F}.$$

$${\mathbf{\Phi }}_2={\mathbb{A}}_2^{+}\mathcal{F},$$

$${\mathbb{A}}_3{\mathbf{\Phi }}_3=\Delta {\mathcal{F}}_3,$$

$$\tilde{\mathbf{\Phi }}=\underset{\mathbf{\Phi }}{\arg \min }\Vert \mathbb{A}\mathbf{\Phi }-\mathcal{F}{\Vert }_2^2.$$

$$\tilde{\mathbf{\Phi }}=\underset{\mathbf{\Phi }}{\arg \min }\Vert \mathbb{A}\mathbf{\Phi }-\mathcal{F}{\Vert }_2^2+\lambda \Vert \mathbf{\Phi }{\Vert }_1.$$

$$A_q(\omega )\propto \frac{4{\Omega }_q^2{\Gamma }_q^{\mathrm{bubble}}(\omega )}{\{{\omega }^2-{\Omega }_q^2-2{\Omega }_q{\Delta }_q^{\mathrm{bubble}}(\omega ){\}}^2+\{2{\Omega }_q{\Gamma }_q^{\mathrm{bubble}}(\omega ){\}}^2}.$$

Full text

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.Ziman (1960) With the advent of increased computational power and the development of computational methods based on density functional theory (DFT),Hohenberg and Kohn (1964); Kohn and Sham (1965) 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.Wallace (1972) The HA is in many cases valid and useful for obtaining phonon dispersion curves and discussing phase stability based on the vibrational free energy. Baroni et al. (2001) 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). 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, 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 heliumHorner (1967) and 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. 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 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. 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 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 SrTiO3, which is realized only at temperatures above 105 K.Yamada and Shirane (1969) 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 exchange-correlation functional is demonstrated. Finally, we make concluding remarks in Sect. IV.

II Methodology

II.1 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:

[TABLE]

Here, $u_{\mu}(\ell\kappa)$ is the displacement of atom $\kappa$ in the $\ell$th unit cell along the $\mu~{}(=x,y,z)$ direction and $\Phi_{\mu_{1}\dots\mu_{n}}(\ell_{1}\kappa_{1};\dots;\ell_{n}\kappa_{n})$ is the $n$th-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}=\frac{1}{2}\sum_{\bm{q}j}\hbar\omega_{\bm{q}j}A_{\bm{q}j}A_{\bm{q}j}^{\dagger}$ with harmonic phonon frequency $\omega_{\bm{q}j}$ and associated displacement operator $A_{\bm{q}j}=b_{\bm{q}j}+b_{-\bm{q}j}^{\dagger}$, where $b_{\bm{q}j}$ and $b_{\bm{q}j}^{\dagger}$ are the annihilation and creation operators of a phonon with crystal momentum $\bm{q}$ and branch index $j$, respectively. The phonon frequency $\omega_{\bm{q}j}$ can be obtained from the dynamical matrix

[TABLE]

where $M_{\kappa}$ is the mass of atom $\kappa$ and $\bm{r}(\ell)$ is the translational vector of the $\ell$th unit cell. By diagonalizing the matrix $D(\bm{q})$, one obtains the squared harmonic frequency and the polarization vector as $D(\bm{q})\bm{e}_{\bm{q}j}$ $=\omega_{\bm{q}j}^{2}\bm{e}_{\bm{q}j}$. Since the harmonic IFCs are temperature-independent, the intrinsic temperature dependence of phonon frequencies and polarization vectors $\{\omega_{\bm{q}j},\bm{e}_{\bm{q}j}\}$ is absent within the HA. The atomic displacement $u_{\mu}(\ell\kappa)$ can be represented in terms of $A_{\bm{q}j}$ as follows:

[TABLE]

Here, $N$ is the number of $\bm{q}$ points in the first Brillouin zone (BZ) and $e_{\mu}(\kappa;\bm{q}j)$ is a component of $\bm{e}_{\bm{q}j}$. By substituting Eq. (4) into Eq. (2), we obtain

[TABLE]

Here and in the following, we use $q$ for the shorthand notation of $(\bm{q},j)$, satisfying $q=(\bm{q},j)$ and $-q=(-\bm{q},j)$. The term $\Delta(\bm{q})$ is 1 if $\bm{q}$ is a vector of the reciprocal lattice $\bm{G}$ and 0 otherwise. The coefficient $V(q_{1};\dots;q_{n})$ is defined as

[TABLE]

Equation (5) is the $n$th-order potential energy represented in terms of the harmonic displacement operator $A_{q}$.

II.2 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^{\prime}$ of the non-interacting Hamiltonian $H_{0}$ as

[TABLE]

Here, we omitted fifth- and higher-order terms since their contributions are much smaller than the cubic and quartic terms. Let $\bm{G}_{q}(\omega)$ be the one-phonon Green’s function and $\bm{G}_{q}^{0}(\omega)$ be that of the non-interacting system $H_{0}$. Then, the following Dyson equation holds:

[TABLE]

with $\bm{\Sigma}_{q}(\omega)$ being the anharmonic self-energy, which can be estimated within a systematic diagrammatic approximation. If the self-energy correction is small and the condition $\hbar\omega_{q}\gg|\Sigma_{q}(\omega_{q})|$ is well satisfied, the phonon quasiparticle picture is still valid and the phonon frequency shift $\Delta_{q}$ and the linewidth $\Gamma_{q}$ are given as $\Delta_{q}=-\frac{1}{\hbar}\mathrm{Re}\Sigma_{q}(\omega_{q})$ and $\Gamma_{q}=\frac{1}{\hbar}\mathrm{Im}\Sigma_{q}(\omega_{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:Maradudin and Fein (1962)

[TABLE]

Here, $n_{i}=n(\omega_{q_{i}})=1/(e^{\beta\hbar\omega_{q_{i}}}-1)$ is the Bose-Einstein distribution function and $\omega_{c}=\omega+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 $\bm{q}_{1}+\bm{q}_{2}=\bm{q}+\bm{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 $\bm{0}j_{1}$ because of the singularity of the matrix element $V(-q;q;\bm{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.Paulatto et al. (2015) The bubble diagram provides the dominant contribution to the phonon linewidth $\Gamma_{q}^{\mathrm{bubble}}(\omega_{q})$, and the quasiparticle lifetime $\tau_{q}$ can be estimated as $\tau_{q}=1/(2\Gamma_{q}^{\mathrm{bubble}}(\omega_{q}))$.

The loop diagram (c) is the first-order contribution from the quartic terms, whose expression is

[TABLE]

Since $\Sigma_{q}^{\mathrm{loop}}$ 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.Tripathi and Pathak (1974); Monga and Pathak (1978); Procacci et al. (1992)

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,Narasimhan and Vanderbilt (1991); Debernardi et al. (1995). Recently, there have been many DFT-based estimations of the bubble diagram mostly for the thermal conductivity prediction described in Sect. 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 $\Sigma_{q}^{\mathrm{loop}}$ relies on quartic IFCs, which are technically more difficult to obtain from DFT than harmonic and cubic IFCs.Rousseau and Bergara (2010); Errea (2016) In Sect. 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 SCP theory, originally developed by HootonHooton (1958) and followed by more elegant formulations based on the many-body theory,Koehler (1966); Horner (1967); Werthamer (1970) is a nonperturbative approaches that can treat the anharmonic renormalization of phonon frequencies. To derive the SCP equation, we rewrite the Hamiltonian [Eq. (8)] as

[TABLE]

Here, $\mathscr{H}_{0}$ is the effective harmonic Hamiltonian, which can be written as $\mathscr{H}_{0}=\frac{1}{2}\sum_{q}\hbar\Omega_{q}\mathcal{A}_{q}\mathcal{A}_{q}^{\dagger}$ with the renormalized phonon frequency $\Omega_{q}$ and the associated displacement operator $\mathcal{A}_{q}$. Then, the free energy of the original system $F$ can be written as the cumulant expansion

[TABLE]

Here, $\mathscr{F}_{0}=-\frac{1}{\beta}\log{Z}$ and $\braket{X}_{\mathscr{H}_{0}}=Z^{-1}\mathrm{Tr}{(Xe^{-\beta\mathscr{H}_{0}})}$ with the partition function $Z=\mathrm{Tr}{(e^{-\beta\mathscr{H}_{0}})}$, and $\braket{X^{n}}_{\mathrm{c}}$ indicates the $n$th-order 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:

[TABLE]

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 $\{\Omega_{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 $\partial F_{1}/\partial\Omega_{q}=0$ with $F_{1}=\mathscr{F}_{0}+\braket{H-\mathscr{H}_{0}}_{\mathscr{H}_{0}}$, the following SCP equation can be readily derived:

[TABLE]

By solving Eqs. (16) and (17) self-consistently, anharmonic phonon frequencies $\{\Omega_{q}\}$ can be obtained. The equation is valid even if some of the harmonic phonon frequencies are unstable ($\omega_{q}^{2}<0$). These unstable modes are renormalized by the second term $2\Omega_{q}I_{q}$ to give a stable phonon mode ($\Omega_{q}^{2}\geq 0$). Since the formulation assumes the existence of a well-defined Hamiltonian $\mathscr{H}_{0}$, an imaginary frequency ($\Omega_{q}^{2}<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. Hermes and Hirata (2013, 2014); Tadano and Tsuneyuki (2015) 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 $\omega_{q}$ for $\Omega_{q}$ in the second term of Eq. (16) gives the correct perturbation limit, where we have $\Omega_{q}=(\omega_{q}^{2}-2\omega_{q}\Sigma_{q}^{\mathrm{loop}}/\hbar)^{\frac{1}{2}}\approx\omega_{q}-\frac{1}{\hbar}\Sigma_{q}^{\mathrm{loop}}=\omega_{q}+\Delta_{q}^{\mathrm{loop}}$.

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 cellErrea 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 extended 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 applied the extended method to cubic SrTiO3 and demonstrated the importance of the PM for the soft modes at high-symmetry $\bm{q}$ 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 evaluate the finite-size effect by increasing the $\bm{q}_{1}$ mesh in Eq. (17) based on the Fourier interpolation, which can be significant near a structural phase transition.Cowley (1996); Tadano and Tsuneyuki (2015) However, as already mentioned in Sect. II.2, 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 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 dichalcogenides,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) 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.

II.4 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,Peierls (1955) the LTC is given as

[TABLE]

where $V$ is the unit cell volume and $\bm{v}_{q}=\frac{\partial\omega_{q}}{\partial\bm{q}}$ is the group velocity. The vector $\bm{f}_{q}$ is a solution of the following linear equation:

[TABLE]

Here, $\Lambda_{q}^{q^{\prime}}$ and $\Lambda_{qq^{\prime}}^{q^{\prime\prime}}(\Lambda_{q}^{q^{\prime}q^{\prime\prime}})$ are the intrinsic transition rates of phonon $q$ involving two and three phonons, respectively. For example, $\Lambda_{q}^{q^{\prime}}$ includes the effect of phonon-isotope scattering and $\Lambda_{qq^{\prime}}^{q^{\prime\prime}}(\Lambda_{q}^{q^{\prime}q^{\prime\prime}})$ 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 $\kappa$ values that agree well with experimental data.Broido et al. (2007) 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,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 the LTC is simplified to

[TABLE]

In the RTA, the current vertex corrections are neglected and the transport relaxation time is approximated by the quasiparticle lifetime $\tau_{q}$.Sun and Allen (2010) Therefore, the RTA incorrectly treats the normal process of three-phonon scattering as resistive, which results in $\kappa_{\mathrm{RTA}}<\kappa$. The underestimation by the 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, 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.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 from estimating the LTC of high-temperature phases because of the imaginary modes within the HA. To overcome this limitation, we made the following extension:Tadano and Tsuneyuki (2015)

[TABLE]

where $\Omega_{q}$ is the SCP frequency, $\tilde{\bm{v}}_{q}=\frac{\partial\Omega_{q}}{\partial\bm{q}}$ is the group velocity of the renormalized phonon, and $\tilde{n}_{q}=n(\Omega_{q})$. The lifetime of the renormalized phonon $\tilde{\tau}_{q}$ 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 $\{\omega_{q},\bm{e}_{q}\}$ are replaced with the SCP solution $\{\Omega_{q},\bm{\epsilon}_{q}\}$. We have successfully applied the new method to predict the LTC of cubic SrTiO3 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 Bi2Te3,Hellman and Broido (2014) PbTe,Romero et al. (2015) and SrTiO3,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 indicated the importance of the renormalization of IFCs in the LTC of the studied materials.

Incorporating the four-phonon scattering process $\Lambda_{q}^{q^{\prime}q^{\prime\prime}q^{\prime\prime\prime}}$ into the BTE [Eq. (19)] is still a major computational challenge even within the RTA. In order to achieve this, the higher-order diagram shown in Fig. 1(d) must be calculated with a dense $\bm{q}$ 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 a robust understanding of the four-phonon 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 MD-based methods, such as the non-equilibrium MDEvans (1982); Müller-Plathe (1997) and Green-Kubo methods,Kubo et al. (1957) 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. Schelling et al. (2002); He et al. (2012); Tadano et al. (2014) 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.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, 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.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 methods 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 $\bm{q}$-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 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 first-order numerical derivative of the atomic forces as

[TABLE]

or inversely asParlinski et al. (1997)

[TABLE]

which are valid when the displacement $u_{\mu}(\ell\kappa)$ is small, say $u_{\mu}(\ell\kappa)\sim$ 0.01 Å. $F_{\nu}[\ell^{\prime}\kappa^{\prime};u_{\mu}(\ell\kappa)]$ is the atomic force acting on the $\kappa^{\prime}$th atom in the $\ell^{\prime}$th cell along direction $\nu$ when the displacement of $u_{\mu}(\ell\kappa)$ is made with all other atoms kept at the equilibrium positions. Equation (23) is defined for every single coefficient $\Phi_{\mu\nu}(\ell\kappa;\ell^{\prime}\kappa^{\prime})$, many of which are linearly dependent on 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 $\bm{u}$ and $\bm{F}$, which consist of $3s$ atomic displacements and forces in the supercell containing $s$ atoms, respectively. We also denote the $i$th displacement pattern as $\bm{u}_{i}$ and the associated force vector as $\bm{F}_{i}$. Then, denoting a column vector comprising $n$ unique harmonic coefficients by $\bm{\Phi}_{2}$, the set of linear equations [Eq. (23)] can be symbolically written as

[TABLE]

Here, $\mathbb{A}_{2}$ is a $3sm\times n$ matrix defined as $\mathbb{A}_{2}^{T}=[A_{2}^{T}(\bm{u}_{1}),\dots,A_{2}^{T}(\bm{u}_{m})]$, $A_{2}(\bm{u})=-\partial^{2}U_{2}/\partial\bm{\Phi}_{2}^{T}\partial\bm{u}$ is a $3s\times n$ matrix, $\bm{\mathcal{F}}^{T}=[\bm{F}_{1}^{T},\dots,\bm{F}_{m}^{T}]$, 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 $\bm{\mathcal{F}}$. If we consider sufficient displacement patterns and make the matrix $\mathbb{A}_{2}$ full rank, the harmonic IFCs can be obtained, in the least-squares sense, asParlinski et al. (1997); Chaput et al. (2011)

[TABLE]

where $\mathbb{A}^{+}_{2}$ is the pseudoinverse of the matrix $\mathbb{A}_{2}$, which can be readily obtained by singular value decomposition. A similar equation can be formed for the cubic terms asChaput et al. (2011)

[TABLE]

in which $\Delta\bm{\mathcal{F}}_{3}=\bm{\mathcal{F}}-\mathbb{A}_{2}\bm{\Phi}_{2}$ is the cubic contribution to the atomic forces and $\mathbb{A}_{3}^{T}=[A_{3}^{T}(\bm{u}_{1}),\dots,A_{3}^{T}(\bm{u}_{m})]$ with $A_{3}(\bm{u})$ defined as $A_{3}(\bm{u})=-\partial^{2}U_{3}/\partial\bm{\Phi}_{3}^{T}\partial\bm{u}$. Again, Eq. (26) assumes that the quartic and higher-order contributions to atomic forces $\bm{\mathcal{F}}$ are negligible, which is valid when the displacement is small. To make the matrix $\mathbb{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 short-range, which is usually the case in covalent compounds. By contrast, long-range cubic interactions 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 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 $\mathbb{A}_{n}(n\geq 4)$ full rank. Moreover, an accurate DFT calculation of $\Delta\bm{\mathcal{F}}_{n}$ is difficult because each $\Delta\bm{\mathcal{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) problemTadano et al. (2014)

[TABLE]

Here, $\bm{\Phi}$ is the column vector defined as $\bm{\Phi}^{T}=[\bm{\Phi}_{2}^{T},\bm{\Phi}_{3}^{T},\bm{\Phi}_{4}^{T},\dots]$ and $\mathbb{A}^{T}=[A^{T}(\bm{u}_{1}),\dots,A^{T}(\bm{u}_{m})]$ with $A(\bm{u})=-\partial^{2}U/\partial\bm{\Phi}^{T}\partial\bm{u}$. Given a sufficient number of AIMD snapshots, the unknown IFCs can be estimated as $\tilde{\bm{\Phi}}=\mathbb{A}^{+}\bm{\mathcal{F}}$. Using this approach, we successfully extracted anharmonic IFCs of Si and Mg2Si 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 using the least absolute shrinkage and selection operator (LASSO)Tibshirani (1996):

[TABLE]

Owing to the $\ell_{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 noiseCandès et al. (2006) and a penalty term helps to avoid the overfitting inherent in the OLS method. The coefficient $\lambda$ is a hyperparameter that controls the trade-off between the sparsity and accuracy of the model; all IFCs will be zero in the large-$\lambda$ limit and $\lambda=0$ corresponds to the OLS regression. An appropriate value of $\lambda$ can be estimated, for instance, by cross-validation, as will be demonstrated in Sect. III.

In Ref. Zhou et al., 2014, Zhou et al. extracted harmonic and anharmonic IFCs of Si, NaCl, and Cu12Sb4S13 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 SrTiO3,Tadano and Tsuneyuki (2015) superconducting hydrogen sulfides,Sano et al. (2016) and thermoelectric clathrates and applied them to SCP and BTE calculations. Also, compressed sensing has successfully been employed to develop interatomic potentialsSeko et al. (2014); *Seko:2015gq and a 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 SrTiO3 (c-STO), which is stable above 105 K. All of the DFT calculations were performed using the Vienna Ab initio Simulation Package (vasp),Kresse 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 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)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 were basically the same as those in the previous PBEsol study,Tadano and Tsuneyuki (2015) except that BZ integration was performed with the 8$\times$8$\times$8 Monkhorst-Pack $\bm{k}$-point grid. All of the phonon calculations reported here were performed with a 2$\times$2$\times$2 cubic supercell (40 atoms) using the alamode package,Tadano ; Tadano et al. (2014) an open-source software developed for modeling lattice anharmonicity and thermal conductivity.

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 2$\times$2$\times$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.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 sextic terms, respectively, which in total form a parameter vector $\bm{\Phi}$ of length 3018. Calculating all elements of $\bm{\Phi}$ 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. Tadano and Tsuneyuki, 2015. After that, the data set was split into four smaller subsets and the hyperparameter $\lambda$ in Eq. (28) was selected by 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 $\mathbb{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 $\|\mathbb{A}\tilde{\bm{\Phi}}-\bm{\mathcal{F}}\|^{2}_{2}/\|\bm{\mathcal{F}}\|_{2}^{2}$. In the large -$\lambda$ region, the difference between the training error and the validation error is marginal. With decreasing $\lambda$, the difference becomes more predominant; the training error monotonically decreases, whereas the validation curve has a minimum around $\lambda=2.3\times 10^{-4}$, indicated by the dotted vertical line in the figure. This value of $\lambda$ 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 $\lambda$, we obtained 2024 nonzero coefficients, which is about 67% of the total number of elements. The accuracy of the extracted IFCs $\tilde{\bm{\Phi}}$ 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 off-diagonal elements of the loop diagram as in Ref. Tadano and Tsuneyuki, 2015. The convergence of the anharmonic frequencies $\{\Omega_{q}\}$ at the commensurate 2$\times$2$\times$2 $\bm{q}$-point grid was carefully checked with respect to the intermediate $\bm{q}_{1}$ grid in Eq. (17) and found sufficiently converged results with 8$\times$8$\times$8 $\bm{q}_{1}$ points. In addition, the non-analytic correction to the harmonic dynamical matrix was included to reproduce the LO-TO splitting around the $\Gamma$ point. To this end, we employed the Born effective charges and the dielectric constant reported in Ref. Tadano and Tsuneyuki, 2015 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.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 (0) 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 ($\omega_{q}^{2}<0$) occur at points $\Gamma\;(0,0,0)$ and $\mathrm{R}\;(\frac{1}{2},\frac{1}{2},\frac{1}{2})$, which correspond to the ferroelectric (FE) and antiferrodistortive (AFD) modes, respectively. The harmonic phonon frequency of the FE mode is 101$i$ cm*-1* with HSE, whose absolute value is larger than the PBEsol result of 58$i$ 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 76$i$ cm*-1* (PBEsol) to 25$i$ cm*-1* (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 to our previous result in Ref. Tadano and Tsuneyuki, 2015. The SCP frequency of the FE mode $\Omega_{\mathrm{FE}}$ is 133 cm*-1* at 300 K with PBEsol, which overestimates the INS data of 87.1$\pm$5.6 cm*-1* (Ref. Yamada and Shirane, 1969, 293 K) and the IR data of 89 cm*-1* (Ref. Servoin et al., 1980, 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, $\Omega_{\mathrm{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.

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 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 $\Delta_{q}^{\mathrm{bubble}}=-\frac{1}{\hbar}\mathrm{Re}\tilde{\Sigma}_{q}^{\mathrm{bubble}}(\Omega_{q})$. The results of $\Omega_{q}+\Delta_{q}^{\mathrm{bubble}}$ are shown in Table 1 as “SCP+B” (B stands for bubble). The correction $\Delta_{q}^{\mathrm{bubble}}$ 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, $\Delta_{q}^{\mathrm{bubble}}\approx-29$ cm*-1* with PBEsol, and the corrected frequency $\Omega_{q}+\Delta_{q}^{\mathrm{bubble}}$ was almost zero. It is interesting to note that the frequency dependence of $\mathrm{Re}\tilde{\Sigma}_{q}^{\mathrm{bubble}}(\omega)$ was marginal in the low-frequency region ($\omega<200$ cm*-1*) and the difference between $\Delta_{q}^{\mathrm{bubble}}$ and $-\frac{1}{\hbar}\mathrm{Re}\tilde{\Sigma}_{q}^{\mathrm{bubble}}(0)$ was minor, especially when the SCP frequency was small. Therefore, it may be reasonable to make the static approximation of $\tilde{\Sigma}_{q}^{\mathrm{bubble}}(\omega)\approx\tilde{\Sigma}_{q}^{\mathrm{bubble}}(0)$ in order to calculate the frequency shift due to the bubble diagram, as carried out in Ref. Shukla and Cowley, 1998. 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}(\omega)=-\frac{1}{\pi}\mathrm{Im}G_{q}(\omega)$ using the following formula:

[TABLE]

Since $A_{q}(\omega)$ 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 results when the condition of $\hbar\Omega_{q}\gg|\Sigma_{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 $\sim 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 $\Gamma$ and along path R–M.

III.4 Anomalous thermal transport in SrTiO3

Finally, we calculated the LTC of c-STO using Eq. (21) with 20$\times$20$\times$20 $\bm{q}$-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 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 ($\kappa\propto T^{-\alpha}$, $\alpha\sim$ 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 $\kappa\propto 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. 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*-1*,Feng et al. (2015); Tadano and Tsuneyuki (2015) as can be seen in the inset of Fig. 6. Here, $\kappa(\omega)$ is the thermal conductivity spectrum defined as $\kappa(\omega)=\sum_{q}\kappa_{q}\delta(\omega-\omega_{q})$ with $\kappa_{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 soft-mode frequencies and LTC values of cubic SrTiO3 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.

Acknowledgements

This study was 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 was carried out using the facilities of the Supercomputer Center, Institute for Solid State Physics, The University of Tokyo.