# Renormalized Phonon Microstructures at High Temperatures from First-Principles Calculations: Methodologies and Applications in Studying Strong Anharmonic Vibrations of Solids

Authors: Tian Lan; Zhaoyan Zhu
Journal: Advances in Condensed Matter Physics
 (2016)

---

## Abstract

While the vibrational thermodynamics of materials with small anharmonicity at low temperatures has been understood well based on the harmonic phonons approximation, at high temperatures, this understanding must accommodate how phonons interact with other phonons or with other excitations. To date the anharmonic lattice dynamics is poorly understood despite its great importance, and most studies still rely on the quasiharmonic approximations. We shall see that the phonon-phonon interactions give rise to interesting coupling problems and essentially modify the equilibrium and nonequilibrium properties of materials, for example, thermal expansion, thermodynamic stability, heat capacity, optical properties, thermal transport, and other nonlinear properties of materials. The review aims to introduce some recent developements of computational methodologies that are able to efficiently model the strong phonon anharmonicity based on quantum perturbation theory of many-body interactions and first-principles molecular dynamics simulations. The effective potential energy surface of renormalized phonons and structures of the phonon-phonon interaction channels can be derived from these interdependent methods, which provide both macroscopic and microscopic perspectives in analyzing the strong anharmonic phenomena while the traditional harmonic models fail dramatically. These models have been successfully performed in the studies on the temperature-dependent broadenings of Raman and neutron scattering spectra, high temperature phase stability, and negative thermal expansion of rutile and cuprite structures, for example.


---

## Body

## 1. Introduction

Today our understanding of the vibrational thermodynamics of materials at low temperatures is developing nicely, based on the harmonic model in which phonons are independent. At high temperatures, however, this understanding must accommodate how phonons interact with other phonons (so called anharmonic phonon-phonon interactions) or with electron or magnon excitations. These anharmonic processes and thermal excitations induce the frequency shifts and lifetime broadenings of the interacting quasiparticles and contribute to most thermal properties at high temperatures.The anharmonic phonon couplings and excitations are mostly as a result of the nonlinear terms in the atomic potentials. The anharmonicity can be much more prominent in some materials than others owing to their physical and chemical properties. For example, an open lattice structure tends to give rise to large anharmonicity in general because of the flexibility adapting to atomic vibrations especially for some rotational modes. Some particular responses of electronic band structures to the lattice variation may cause strong phonon anharmonicity too. Typical examples may include rigid-unit modes, rattling ions in structural cages, orbital driven ferroelectric-like lattice instability, charge disproportionation, valence fluctuations, and vibrational driven hybridization [1–6]. At low temperatures the harmonic part is dominant in most materials but examples of extremely anharmonic potentials are also seen. However, even in the former cases, with elevated temperatures, the contribution of anharmonicity is growing fast and can become substantial owing to its nonlinear nature. The anharmonicity renormalizes the vibrational quanta and serves to break the assumptions of harmonicity and independence in phonon dynamics and substantially changes the characteristics of phonons.These topics are rich and of great importance for the rational design and engineering of next-generation materials in energy based applications. For example, the anharmonic dynamics and other types of thermal excitations are the origin of most thermal energy transport processes and therefore greatly influence the performance of these materials in applications of harvesting, storing, and transporting energy. A reliable estimate of the anharmonic entropy is also crucial for synthesizing materials. For example, for metals and oxides, it seems that pure anharmonic contributions become large enough to affect phase stability at temperatures above half the melting temperature, which is the temperature range where materials are often processed or used [7]. Anomaly in thermal expansion is another prominent example. Recently, the large negative thermal expansion (NTE) of ScF3 and Ag2O was found to have strong dependence with these high temperature vibrational dynamical properties [8, 9].Modern inelastic scattering techniques with neutrons or photons are ideal for sorting these properties out. Analysis of the experimental data can generate vibrational spectra of the materials, that is, their phonon densities of states (DOS) and phonon or spin wave dispersions. We are developing the data reduction software to obtain the high quality data from inelastic neutron spectrometers [10–13]. With accurate phonon DOS and dispersion curves we can obtain the vibrational entropies of different materials. The understanding of the underlying reasons for differences in DOS curves and entropies then relies on the development of the fundamental theories and the computational methods.To date, most ab initio methods for calculating materials structures and properties have been based on density functional (DFT) methods and evaluating the internal energy,E, of materials at a temperature of zero Kelvin. For example, a harmonic or quasiharmonic model usually used to account for the vibrational thermodynamics at low temperatures and it is commonplace today to calculate harmonic phonons by methods based on DFT [7, 14, 15]. The quasiharmonic approximation (QHA) is based on how phonon frequencies change with volume and all shifts of phonon frequencies from their low temperature values are considered as a result of thermal expansion alone [7].In the QHA, the vibrational free energy can be minimized as a function of volume:(1)FV,T=E0+∫-∞+∞gωħω2+kBTln⁡1-e-ħω/kBTdω,where E0 is the energy calculated from the relaxed structure at T=0K. Thermodynamic properties are therefore derived from here [7, 16].Although the QHA accounts for some frequency shifts, the phonon modes are still assumed to be harmonic, noninteracting, and their energies depend only on the volume of the crystal. This can be adequate when the temperatures of service of the materials are low or when differences of chemical potentials are much larger than kT. Therefore, QHA has been found to be able to predict thermodynamic properties well consistent with the experiment results especially at high pressure [17, 18]. However, for most applications of materials in energy involving even modest temperature, this E alone is insufficient because the anharmonic vibrational dynamics and different types of thermal excitations play important roles and have significant thermodynamic effects at elevated temperatures.Phonon-phonon interactions are responsible for pure anharmonicity that shortens phonon lifetimes and shifts phonon frequencies especially at high temperatures. Anharmonicity competes with quasiharmonicity to alter the stability of phases at high temperatures, as has been shown, for example, with experiments and frozen phonon calculations on bcc Zr [19] and the possible stabilization of bcc Fe-Ni alloys at conditions of the earth’s core [20]. For PbTe, ScF3, and rutile TiO2, there are recent reports of anharmonicity being so large that both the QHA and anharmonic perturbation theory fail dramatically [6, 8, 21, 22].These cases are suitable for ab initio molecular dynamics (AIMD) simulations, which should be reliable when the electrons are near their ground states and the nuclear motions are classical. The big advantage of ab initio MD is that it can account for all effects of harmonic, anharmonic, and even some of the electron-phonon interactions. However, advanced postprocessing methodologies are requied to extract concrete information from these simulations. In the few examples where comparisons have been made with ab initio MD, agreement has been surprisingly good even for highly anharmonic materials [6, 22, 23]. Today, by validating these calculated results from inelastic scattering experiments with facilities such as the Spallation Neutron Source for neutrons, we can obtain sufficient details about phonon-phonon interactions, electron-phonon interactions, and other excitations at elevated temperatures.To understand the microscopic picture of these interactions, models of the effective vibrational potential energy surface and the fine structures of decay channels of phonon-phonon interactions have been proposed [23–26]. Based on the quantum perturbation theory of many-body interactions and first-principles molecular dynamics simulations, these methods are used to renormalize quasiharmonic phonons and to identify the three-phonon and four-phonon kinematics. We can assess the strengths of phonon-phonon interactions of different anharmonic orders or via different decay channels. These methods, with high computational efficiency, are promising directions to advance our understandings of nonharmonic lattice dynamics and thermal transport properties.In this review article, we discuss several first-principles computational techniques available recently which proved to be useful for assessing the anharmonic vibrational thermodynamics of solids. In particular, the computational details are discussed, followed by concrete examples that demonstrate how the applications of these interdependent methods can unveil interesting anharmonic properties of materials and their relationships with NTE, vibrational energy shift, and phase stability.


## 2. Renormalized Phonon Spectra from Molecular Dynamics Simulations

### 2.1. Molecular Dynamics Simulation and Fourier Transformed Velocity Autocorrelation Method

Pure anharmonicity contributes to phonon-phonon interactions that shorten phonon lifetimes and shift phonon frequencies. The vibrational energy spectra of these renormalized phonons can be produced by velocity trajectories extracted from the MD simulation at each temperature. It is based on nonequilibrium statistics initiated by Green and Kubo [27, 28]. In essence, the FTVAC method transforms the vibrational representation from the time and spatial domain to the corresponding energy and/or momentum domain [29–31]. Because the FTVAC method does not assume a form for the Hamiltonian, it is a robust tool for obtaining vibrational spectra of renormalized phonons from MD simulations, even with strong anharmonicity. In the FTVAC model, the phonon DOS is given by(2)gω=∑n,b∫e-iωtv→n,btv→0,00dt,where 〈〉 is an ensemble average and v→n,b(t) is the velocity of the atom b in the unit cell n at time t. Further projection of the phonon modes onto each k point in the Brillouin zone was performed by computing the phonon power spectrum with the FTVAC method, with a resolution determined by the size of the supercell in the simulation(3)gq→,ω=∫dte-iωt∑n,beiq→·R→nv→n,btv→0,00,where R→n is the equilibrium position of the cell n and q→ is the phonon wavevector. Equation (3) is both a time and space Fourier transform and gives the frequency and lifetime of each phonon mode.


### 2.2. Temperature-Dependent Effective Potential Method

In general, the cubic phonon anharmonicity contributes to both the phonon energy shift and the lifetime broadening, whereas the quartic anharmonicity contributes only to the phonon energy shift [26, 32]. To distinguish the roles of cubic and quartic anharmonicity, the Temperature-Dependent Effective Potential (TDEP) method [6, 24] is used. In the TDEP method, an effective Hamiltonian model is used to sample the potential energy surface, not at the equilibrium positions of atoms, but at the most probable positions for a given temperature in an MD simulation [24](4)H=U0+12∑impi2+12∑ijαβϕijαβuiαujβ+13!∑ijkαβγψijkαβγuiαujβukγ,where ϕij and ψijk are second- and third-order force constants, p is momentum, and uiα is the Cartesian component α of the displacement of atom i. In the fitting, the “effective” harmonic force constants ϕij are renormalized by the quartic anharmonicity. The cubic anharmonicity, however, is largely accounted for by the third-order force constants ψijk and can be understood in terms of the third-order phonon self-energy that causes linewidth broadening [32].The above Hamiltonian was used to obtain the renormalized phonon dispersions (TDEP spectra) accounting for both the anharmonic shiftsΔ and broadenings Γ of the mode q→j. These are derived from the real and imaginary parts of the cubic self-energies Σ(3), respectively [32](5)Δq→j;Ω=-18ħ2∑q→1j1∑q→2j2Vq→j;q→1j1;q→2j22Δq→1+q→2-q→℘n1+n2+1Ω+ω1+ω2-n1+n2+1Ω-ω1-ω2+n1-n2Ω-ω1+ω2-n1-n2Ω+ω1-ω2,Γq→j;Ω=18πħ2∑q→1j1∑q→2j2Vq→j;q→1j1;q→2j22Δq→1+q→2-q→n1+n2+1δΩ-ω1-ω2+2n1-n2δΩ+ω1-ω2,where Ω is the renormalized phonon frequency and ℘ denotes the Cauchy principal part. The V(·)’s are elements of the Fourier transformed third-order force constants ψijk obtained in the TDEP method. Δ(q→1+q→2-q→) ensures conservation of momentum.


### 2.3. Example: The Quartic Phonons and Their Stabilization of Rutile Phase of TiO2 at High Temperatures

Although the rutile structure ofTiO2 is known to be stable at high temperatures, the QHA predicts that several acoustic phonons decrease anomalously to zero frequency with thermal expansion, incorrectly predicting a structural collapse at temperatures well below 1000 K [33, 34].Inelastic neutron scattering was used to measure the temperature dependence of the phonon density of states (DOS) of rutileTiO2 from 300 to 1373 K. Surprisingly, as shown in Figure 1, these anomalous acoustic phonons centered at 14 meV (peak 1 in the Figure) were found to increase in frequency with temperature instead [6]. The green curves in Figure 1 present the corresponding peak 1 of DOS spectra calculated with the quasiharmonic model at 300 K and 1373 K. It shows that QHA calculations predict that the translational acoustic (TA) branch softens dramatically to zero frequency with the thermal expansion of 1373 K, giving imaginary frequencies in the DOS that would destabilize the rutile structure at high temperatures.Figure 1
Neutron weighted phonon DOS of rutileTiO2 from experimental measurements at temperatures from 300 to 1373 K with an incident energy of 30 meV (black). The dashed spectrum corresponds to the experimental result at 300 K, shifted vertically for comparison at each temperature. Simulation results of peak 1 of phonon DOS from the QHA at 300 K and 1373 K are shown in green and compared with experimental spectra.Our QHA calculation is consistent with previous computational results [33, 34], but in obvious contrary to the experimetal observation. The dramatic failure of QHA suggests the existence of strong anharmonicity that could intrinsically alter the harmonic characteristics of these acoustic phonons subject to the phonon renormalization. It turns out that the FTVAC and TDEP methods based on the first-principles MD simulations are reliable for this investigation.For MD, first-principles calculations using the local density approximation (LDA) of density functional theory (DFT) were performed with the VASP package [35, 36]. First-principles Born-Oppenheimer AIMD simulations for a 2×2×4 supercell and a 2×2×1k point sampling were performed to thermally excite phonons to the target temperatures of 300 and 1373 K. For each temperature, the system was first equilibrated for 3 ps as an NVT ensemble with temperature control by a Nosé thermostat and then simulated as an NVE ensemble for 20 ps with time steps of 1 fs. Fine relaxations with residual pressures below 0.5 GPa were achieved in each calculation that accounted for thermal expansion. With the simulated atomic trajectories derived from MD simulations, the FTVAC and TDEP methods were able to reproduce the renormalized phonon dispersions and effective vibrational potential surface as detailed in Sections 2.1 and 2.2.Figure2 shows the vibrational energies of the TA branch, calculated by the FTVAC method with AIMD trajectories. From 300 to 1373 K, the TA branch increases in energy by an average of about 2.1 meV. For this TA branch, Figure 2(b) shows an enormous discrepancy of phonon energies between the FTVAC calculation and the QHA (orange dashed line) at 1373 K. Apparently, the unstable phonon modes predicted by the QHA are fully stable in the AIMD simulations at high temperatures.Figure 2
Diffuse curves are TDEP phonon dispersions of rutileTiO2 below 25 meV at (a) 300 K and (b) 1373 K, compared with the results from the FTVAC method (red circles). The white curves are phonon dispersions for the quasiharmonicity plus quartic anharmonicity calculated with all ψijk set to zero in (4). In (b), the dispersions are compared to the quasiharmonic dispersions (orange dashed curve) and the single quartic oscillator model (orange triangles). (c) Frozen phonon potential (black) of TA mode at R point with q=(0.5,0,0.5) at 1373 K, showing the harmonic component (red) and quartic component (blue). The low temperature potential surface is also shown (dashed black). Reprinted figure with permission from Lan et al. [6]. Copyright by the American Physical Society.
(a) (b) (c)Using the same MD trajectories as for the FTVAC method, the calculated TDEP dispersions agree well with the FTVAC results as shown in Figures2(a) and 2(b). At 1373 K the TA modes below 20 meV have only small linewidth broadenings, suggesting the small cubic anharmonicity. Furthermore, they are close in energy to those calculated if all ψijk are set to zero in (4), showing the dominance of quartic anharmonicity and the small cubic anharmonicity of the TA modes.For more details about the anomalous anharmonicity of the TA modes, the frozen phonon method was adopted to calculate the potential energy surfaces for specific phonons, as a typical example presented in Figure2(c). It should be emphasized that the frozen phonon method is still based upon the harmonic theory since this method isolates a particular phonon mode and tries to explore the anharmonic potential landscape defined by its own harmonic eigenvectors. The FTVAC and TDEP methods, however, are able to fully account for the couplings of all phonon modes and therefore reproduce the potential energy surface renormalized by the interactions of all other phonons. Nevertheless, the frozen phonon method is a convenient way to evaluate the strength of anharmonicity, although the potential energy surface which is derived from the frozen phonon method is not the real anharmonic potential in general and cannot be used to investigate the phonon-phonon interactions.We shall see that, at 300 K, the frozen phonon potential energy of the TA mode at theR point is nearly quadratic, with a small quartic part. With the lattice expansion characteristic of 1373 K, the potential energy curve transforms to being nearly quartic. In fact, for all modes in the TA branch that were evaluated by the frozen phonon method, the potential energy surface develops a quartic form with lattice expansion. For a quantum quartic oscillator, the vibrational frequency stiffens with temperature owing to the increasing spread between the energy levels [8, 37]. We assessed a high temperature behavior by assigning Boltzmann factors to the different oscillator levels derived from frozen phonon potentials, giving the energies of the quartic TA modes at 1373 K. As shown in Figure 2(b), they are reasonably close to the FTVAC and TDEP results. The similarity of potential surfaces derived from the TDEP method and the frozen phonon method suggests that the isolated frozen phonon potential does not reshape itself much in response to the renormalization as a result of the interactions of other phonons. This is unusual and indicates a global phonon quarticity of acoustic branch in rutile.Based on FTVAC and TDEP, further investigation reveals that the positive and negative displacements of these acoustic modes show a significant accumulation of charge in the Ti–O bond of shorter distance and a depletion in the bond of longer distance. For rutileTiO2, the hybridization follows the atom displacements in thermal phonons instead of a displacive phase transition, and this “phonon-tracked hybridization” changes with lattice parameter [6]. With lattice expansion, the short-range repulsion is weaker, and hybridization favors electrons between the shorter Ti–O pairs in the phonon displacement pattern. The hybridization in the Ti–O bond is very sensitive to interatomic distance, much as has been noticed in the ferroelectric distortion of BaTiO3 [38]. It not only provides a source of extreme phonon anharmonicity, but also provides thermodynamic stability for rutile TiO2. It may occur in other transition metal oxides that show unusual changes of properties with lattice parameter or with structure. Besides altering thermodynamic phase stability, properties such as ferroelectricity and thermal transport will be affected directly.


## 2.1. Molecular Dynamics Simulation and Fourier Transformed Velocity Autocorrelation Method

Pure anharmonicity contributes to phonon-phonon interactions that shorten phonon lifetimes and shift phonon frequencies. The vibrational energy spectra of these renormalized phonons can be produced by velocity trajectories extracted from the MD simulation at each temperature. It is based on nonequilibrium statistics initiated by Green and Kubo [27, 28]. In essence, the FTVAC method transforms the vibrational representation from the time and spatial domain to the corresponding energy and/or momentum domain [29–31]. Because the FTVAC method does not assume a form for the Hamiltonian, it is a robust tool for obtaining vibrational spectra of renormalized phonons from MD simulations, even with strong anharmonicity. In the FTVAC model, the phonon DOS is given by(2)gω=∑n,b∫e-iωtv→n,btv→0,00dt,where 〈〉 is an ensemble average and v→n,b(t) is the velocity of the atom b in the unit cell n at time t. Further projection of the phonon modes onto each k point in the Brillouin zone was performed by computing the phonon power spectrum with the FTVAC method, with a resolution determined by the size of the supercell in the simulation(3)gq→,ω=∫dte-iωt∑n,beiq→·R→nv→n,btv→0,00,where R→n is the equilibrium position of the cell n and q→ is the phonon wavevector. Equation (3) is both a time and space Fourier transform and gives the frequency and lifetime of each phonon mode.


## 2.2. Temperature-Dependent Effective Potential Method

In general, the cubic phonon anharmonicity contributes to both the phonon energy shift and the lifetime broadening, whereas the quartic anharmonicity contributes only to the phonon energy shift [26, 32]. To distinguish the roles of cubic and quartic anharmonicity, the Temperature-Dependent Effective Potential (TDEP) method [6, 24] is used. In the TDEP method, an effective Hamiltonian model is used to sample the potential energy surface, not at the equilibrium positions of atoms, but at the most probable positions for a given temperature in an MD simulation [24](4)H=U0+12∑impi2+12∑ijαβϕijαβuiαujβ+13!∑ijkαβγψijkαβγuiαujβukγ,where ϕij and ψijk are second- and third-order force constants, p is momentum, and uiα is the Cartesian component α of the displacement of atom i. In the fitting, the “effective” harmonic force constants ϕij are renormalized by the quartic anharmonicity. The cubic anharmonicity, however, is largely accounted for by the third-order force constants ψijk and can be understood in terms of the third-order phonon self-energy that causes linewidth broadening [32].The above Hamiltonian was used to obtain the renormalized phonon dispersions (TDEP spectra) accounting for both the anharmonic shiftsΔ and broadenings Γ of the mode q→j. These are derived from the real and imaginary parts of the cubic self-energies Σ(3), respectively [32](5)Δq→j;Ω=-18ħ2∑q→1j1∑q→2j2Vq→j;q→1j1;q→2j22Δq→1+q→2-q→℘n1+n2+1Ω+ω1+ω2-n1+n2+1Ω-ω1-ω2+n1-n2Ω-ω1+ω2-n1-n2Ω+ω1-ω2,Γq→j;Ω=18πħ2∑q→1j1∑q→2j2Vq→j;q→1j1;q→2j22Δq→1+q→2-q→n1+n2+1δΩ-ω1-ω2+2n1-n2δΩ+ω1-ω2,where Ω is the renormalized phonon frequency and ℘ denotes the Cauchy principal part. The V(·)’s are elements of the Fourier transformed third-order force constants ψijk obtained in the TDEP method. Δ(q→1+q→2-q→) ensures conservation of momentum.


## 2.3. Example: The Quartic Phonons and Their Stabilization of Rutile Phase of TiO2 at High Temperatures

Although the rutile structure ofTiO2 is known to be stable at high temperatures, the QHA predicts that several acoustic phonons decrease anomalously to zero frequency with thermal expansion, incorrectly predicting a structural collapse at temperatures well below 1000 K [33, 34].Inelastic neutron scattering was used to measure the temperature dependence of the phonon density of states (DOS) of rutileTiO2 from 300 to 1373 K. Surprisingly, as shown in Figure 1, these anomalous acoustic phonons centered at 14 meV (peak 1 in the Figure) were found to increase in frequency with temperature instead [6]. The green curves in Figure 1 present the corresponding peak 1 of DOS spectra calculated with the quasiharmonic model at 300 K and 1373 K. It shows that QHA calculations predict that the translational acoustic (TA) branch softens dramatically to zero frequency with the thermal expansion of 1373 K, giving imaginary frequencies in the DOS that would destabilize the rutile structure at high temperatures.Figure 1
Neutron weighted phonon DOS of rutileTiO2 from experimental measurements at temperatures from 300 to 1373 K with an incident energy of 30 meV (black). The dashed spectrum corresponds to the experimental result at 300 K, shifted vertically for comparison at each temperature. Simulation results of peak 1 of phonon DOS from the QHA at 300 K and 1373 K are shown in green and compared with experimental spectra.Our QHA calculation is consistent with previous computational results [33, 34], but in obvious contrary to the experimetal observation. The dramatic failure of QHA suggests the existence of strong anharmonicity that could intrinsically alter the harmonic characteristics of these acoustic phonons subject to the phonon renormalization. It turns out that the FTVAC and TDEP methods based on the first-principles MD simulations are reliable for this investigation.For MD, first-principles calculations using the local density approximation (LDA) of density functional theory (DFT) were performed with the VASP package [35, 36]. First-principles Born-Oppenheimer AIMD simulations for a 2×2×4 supercell and a 2×2×1k point sampling were performed to thermally excite phonons to the target temperatures of 300 and 1373 K. For each temperature, the system was first equilibrated for 3 ps as an NVT ensemble with temperature control by a Nosé thermostat and then simulated as an NVE ensemble for 20 ps with time steps of 1 fs. Fine relaxations with residual pressures below 0.5 GPa were achieved in each calculation that accounted for thermal expansion. With the simulated atomic trajectories derived from MD simulations, the FTVAC and TDEP methods were able to reproduce the renormalized phonon dispersions and effective vibrational potential surface as detailed in Sections 2.1 and 2.2.Figure2 shows the vibrational energies of the TA branch, calculated by the FTVAC method with AIMD trajectories. From 300 to 1373 K, the TA branch increases in energy by an average of about 2.1 meV. For this TA branch, Figure 2(b) shows an enormous discrepancy of phonon energies between the FTVAC calculation and the QHA (orange dashed line) at 1373 K. Apparently, the unstable phonon modes predicted by the QHA are fully stable in the AIMD simulations at high temperatures.Figure 2
Diffuse curves are TDEP phonon dispersions of rutileTiO2 below 25 meV at (a) 300 K and (b) 1373 K, compared with the results from the FTVAC method (red circles). The white curves are phonon dispersions for the quasiharmonicity plus quartic anharmonicity calculated with all ψijk set to zero in (4). In (b), the dispersions are compared to the quasiharmonic dispersions (orange dashed curve) and the single quartic oscillator model (orange triangles). (c) Frozen phonon potential (black) of TA mode at R point with q=(0.5,0,0.5) at 1373 K, showing the harmonic component (red) and quartic component (blue). The low temperature potential surface is also shown (dashed black). Reprinted figure with permission from Lan et al. [6]. Copyright by the American Physical Society.
(a) (b) (c)Using the same MD trajectories as for the FTVAC method, the calculated TDEP dispersions agree well with the FTVAC results as shown in Figures2(a) and 2(b). At 1373 K the TA modes below 20 meV have only small linewidth broadenings, suggesting the small cubic anharmonicity. Furthermore, they are close in energy to those calculated if all ψijk are set to zero in (4), showing the dominance of quartic anharmonicity and the small cubic anharmonicity of the TA modes.For more details about the anomalous anharmonicity of the TA modes, the frozen phonon method was adopted to calculate the potential energy surfaces for specific phonons, as a typical example presented in Figure2(c). It should be emphasized that the frozen phonon method is still based upon the harmonic theory since this method isolates a particular phonon mode and tries to explore the anharmonic potential landscape defined by its own harmonic eigenvectors. The FTVAC and TDEP methods, however, are able to fully account for the couplings of all phonon modes and therefore reproduce the potential energy surface renormalized by the interactions of all other phonons. Nevertheless, the frozen phonon method is a convenient way to evaluate the strength of anharmonicity, although the potential energy surface which is derived from the frozen phonon method is not the real anharmonic potential in general and cannot be used to investigate the phonon-phonon interactions.We shall see that, at 300 K, the frozen phonon potential energy of the TA mode at theR point is nearly quadratic, with a small quartic part. With the lattice expansion characteristic of 1373 K, the potential energy curve transforms to being nearly quartic. In fact, for all modes in the TA branch that were evaluated by the frozen phonon method, the potential energy surface develops a quartic form with lattice expansion. For a quantum quartic oscillator, the vibrational frequency stiffens with temperature owing to the increasing spread between the energy levels [8, 37]. We assessed a high temperature behavior by assigning Boltzmann factors to the different oscillator levels derived from frozen phonon potentials, giving the energies of the quartic TA modes at 1373 K. As shown in Figure 2(b), they are reasonably close to the FTVAC and TDEP results. The similarity of potential surfaces derived from the TDEP method and the frozen phonon method suggests that the isolated frozen phonon potential does not reshape itself much in response to the renormalization as a result of the interactions of other phonons. This is unusual and indicates a global phonon quarticity of acoustic branch in rutile.Based on FTVAC and TDEP, further investigation reveals that the positive and negative displacements of these acoustic modes show a significant accumulation of charge in the Ti–O bond of shorter distance and a depletion in the bond of longer distance. For rutileTiO2, the hybridization follows the atom displacements in thermal phonons instead of a displacive phase transition, and this “phonon-tracked hybridization” changes with lattice parameter [6]. With lattice expansion, the short-range repulsion is weaker, and hybridization favors electrons between the shorter Ti–O pairs in the phonon displacement pattern. The hybridization in the Ti–O bond is very sensitive to interatomic distance, much as has been noticed in the ferroelectric distortion of BaTiO3 [38]. It not only provides a source of extreme phonon anharmonicity, but also provides thermodynamic stability for rutile TiO2. It may occur in other transition metal oxides that show unusual changes of properties with lattice parameter or with structure. Besides altering thermodynamic phase stability, properties such as ferroelectricity and thermal transport will be affected directly.


## 3. Fine Structures of Phonon-Phonon Interaction Channels

### 3.1. Computational Methodology

Anharmonicity tensors describe the coupling strengths for phonon-phonon interactions, but a prerequisite is that the phonons in these processes satisfy the kinematical conditions of conservation of energy and momentum as presented in (5). In the phonon-phonon interaction functional, an anharmonicity tensor element for an s-phonon process can be expressed as [39](6)Vj;q→1j1;…;q→s-1js-1=12s!ħ2Ns/2NΔq→1+⋯+q→s-1ωj0ω1⋯ωs-11/2Cj;q→1j1;…;q→s-1js-1,where Δ(q→1+⋯+q→s-1) enforces momentum conservation and the C(·)’s, elements of the s-phonon anharmonic tensor, are expected to be slowly varying functions of their arguments.If the anharmonicity tensor or its average does not vary significantly for different phonon processes, the coupling factor and the kinematic factor are approximately separable in (6). The separation of the anharmonic coupling and the kinematics has been used with success in many studies including recent reports on rutile TiO2 and SnO2 [25, 26]. We consider the term C(j;q→1j1;…;q→s-1js-1) to be a constant of the Raman mode j and use it as a fitting parameter. Although C(j;q→1j1;q→2j2) and C(j;j;q→1j1;-q→1j1) change with j1 and j2, an average over modes, 〈C(·)〉=∑1,2‍C(j;q→1j1;q→2j2)/∑1,2‍1, is needed by the fitting, where 1 and 2 under the summation symbol represent q→iji. We define the cubic and quartic fitting parameters as(7a)Cj3=Cj;q→1j1;q→2j2,(7b)Cj4=Cj;j;q→1j1;-q→1j1.To the leading order of cubic and quartic anharmonicity, the broadening of the Raman peaks is2Γ(3)(j;Ω). The frequency shift of the Raman peaks is ΔQ+Δ(3)+Δ(3′)+Δ(4), where the quasiharmonic part is denoted by ΔQ. These quantities can be written as functions of D(Ω,ω1,ω2) and P(Ω,ω1,ω2), weighted by average anharmonic coupling strengths [25, 26, 40](8a)Γ3j;Ω=πħ64ωj0Cj32∑q→1,j1∑q→2,j2ω1ω2DΩ,ω1,ω2=ωj0Cj32DωΩ,(8b)Δ3j;Ω=-ħ64ωj0Cj32∑q→1,j1∑q→2,j2ω1ω2PΩ,ω1,ω2=ωj0Cj32PωΩ,(8c)Δ3′j=-ħ16Nωj0Cj32∑q→2j2ωj2q→2nq→2j2+12,(8d)Δ4j=ħ8Nωj0Cj4∑q→1j1ωj1q→1nq→1j1+12,where Dω(Ω) and Pω(Ω) are functionals of D(Ω,ω1,ω2) and P(Ω,ω1,ω2) weighted by the kinematics of anharmonic phonon coupling. Dω(Ω) is the so called two-phonon density of states (TDOS) spectra, which characterize the size of the phonon-phonon interaction channels. Δ(3′) is an additional low-order cubic term that corresponds to instantaneous three-phonon processes [32]. It is nonzero for crystals having atoms without inversion symmetry, as in the case for the oxygen atom motions in the A1g mode of rutile. It is much smaller than other contributions, however, owing to symmetry restrictions.


### 3.2. Example  1: Anomalous Temperature-Dependent Broadening ofB2g Mode of Rutile SnO2

Recently, strong anharmonicity of rutileSnO2 was discovered by the high-resolution Raman spectrometer, and the anomaly of the temperature-dependent broadening of B2g mode is a prominent example [26]. As shown in Figure 3(b), at high temperatures above 500 K, the broadening of the B2g mode of rutile SnO2 shows an unusual concave downward curvature, while the other two Raman modes Eg and A1g broaden linearly. Moreoever, at low temperatures, the B2g mode has a much larger linewidth than the other two modes. The linewidth of the B2g mode extrapolated to 0 K is approximately 8 cm-1, whereas the linewidths of the Eg and A1g modes extrapolate to less than 2 cm-1.Figure 3
Temperature dependence of (a) frequency shifts and (b) linewidths as FWHMs, of the Raman modesEg, A1g, and B2g of rutile SnO2. The solid and open symbols represent the experimental data from powder and single crystal samples, respectively. Solid curves are the theoretical fittings with a full calculation of the kinematics of three-phonon and four-phonon processes. The dashed curve was calculated without considering the frequency dependence of D0↓(Ω) and the number of decay channels, at elevated temperatures. Reprinted figure with permission from Lan et al. [26]. Copyright by the American Physical Society.
(a) (b)The TDOS functionDω(Ω) in Figure 4 shows large variations with Ω which explains a trend in the thermal broadening of Figure 3(b). Owing to the high frequency of the B2g mode, at the temperatures of this study its phonon-phonon anharmonicity comes mostly from downconversion processes, as shown in Figure 4. Ignoring the small upconversion contribution, at high temperatures(9)DωΩ,T∝T∑q→1,j1∑q→2,j2δΩ-ω1-ω2≡TD0↓Ω,where D0↓(Ω) is the number of two-phonon downconversion channels. Therefore, the linewidth broadening resulting from the cubic anharmonicity is proportional to the temperature, modulated by C(Ω). Because Ω undergoes a shift with temperature, D0↓(Ω) is an implicit function of T. Usually the line broadening is linear in T because D0↓(Ω) does not vary much with temperature, as in the cases of modes Eg and A1g. However, the B2g mode at 774 cm-1 lies on a steep gradient of Dω(Ω,T) in Figure 4. On the other hand, the B2g mode undergoes a significant shift of frequency by more than 25 cm-1 with temperature as shown in Figure 3(a), which moves the downconversion channels down the gradient of Dω(Ω,T). These cooperative effects as a result narrow the capacity of these decay channels. Therefore, its broadening is significantly less than linear in T, exhibiting an unusual concave downwards shape. For comparison, the dashed line in Figure 3(b) was calculated without considering the frequency dependence of D0↓(Ω) at elevated temperatures, and it deviates substantially from the experimental trend.Figure 4
Two-phonon density of statesDω(Ω) of rutile SnO2 for 0 K and 800 K. The arrowheads mark the positions of the three Raman modes, Eg, A1g, and B2g, respectively. The upconversion and downconversion contributions at 800 K are shown in green and black dashed curves, respectively. There is no upconversion process at 0 K. The phonon DOS spectra are presented along with the TDOS but take a negative sign and shaded blue to make the view of TDOS and DOS separated.As discssed, the unusual temperature dependence of the linewidth of theB2g mode originates from the sharp peak in Dω(Ω,T) centered at 800 cm-1. We shall see that this feature in the TDOS originates with the phonon DOS of SnO2 which has a wide band gap between 360 cm-1 and 450 cm-1 associated with the mass difference of Sn and O atoms, as shown in Figure 4. The shape of the TDOS can be understood as the convolution of the phonon DOS with itself. With two approximately equal regions above and below the gap, the result is a peak at 800 cm-1, with steep slopes on both sides. For comparison, although the TDOS of rutile TiO2 is shaped as a broad peak [25], it does not have the sharp features of Figure 4 because the mass difference between Ti and O atoms does not cause a band gap in the phonon DOS of TiO2.The fine structures of phonon-phonon interaction channels can also explain the physics of large linewidth ofB2g mode at low temperature limit following the same approach as detailed above. In the low temperature limit, upconversion processes are prohibited because n→0. The peak linewidth extrapolated to 0 K is determined entirely by downconversion processes, quantified by the downconversion TDOS. The B2g mode has a significant broadening because its frequency is near a peak in the downconversion TDOS, as seen in Figure 4. On the other hand, the Eg and A1g modes are not broadened at low temperatures because their frequencies are at low values of the TDOS. The study on the structures of phonon-phonon interaction channels successfully marks the “hot” downconversion decay activities even at low temperature limit and accounts for the significant difference of linewidths between B2g and the other two Raman modes at low temperatures.


### 3.3. Example  2: Negative Thermal Expansion of Cuprite Ag2O

Silver oxide (Ag2O) with the cuprite structure has attracted much interest after the discovery of its extraordinarily large negative thermal expansion (NTE) [41, 42], which exceeds -1×10-5 K-1 and occurs over a wide range of temperature from 40 K to its decomposition temperature near 500 K.A rigid-unit modes (RUMs) model of NTE considers tetrahedra ofAg4O around each O atom that bend at the Ag atoms linking the O atoms in adjacent tetrahedra. RUMs account for counteracting rotations of all such tetrahedra [1, 2]. Locally, the O–Ag bond length does not contract, but bending of the O–Ag–O links pulls the O atoms together, leading to NTE. These RUMs tend to have low frequencies owing to the large mass of the unit and hence are excited at low temperatures. This model correlated the NTE with quasiharmonic approximation and should explain the main physics at low temperatures. However, as shown in Figure 5, at temperatures above 250 K, there is a second part of the NTE behavior of cuprite Ag2O that is apparently beyond the predictions of quasiharmonic theory [9].Figure 5
Temperature dependence of lattice parameter ofAg2O from experimental data in [41], quasiharmonic calculations, and MD calculations, expressed as the relative changes with respect to their 40 K values, that is, a(T)/a(40K)-1.First-principles Born-Oppenheimer molecular dynamics simulations were performed for a3×3×3 supercell with temperature control by a Nosé thermostat. The simulated temperatures included 40, 100, 200, 300, and 400 K. For each temperature, the system was first equilibrated for 3 ps and then simulated for 18 ps with a time step of 3 fs. The system was fully relaxed at each temperature, with convergence of the pressure within 1 kbar.As shown in Figure5, the MD simulation predicts the NTE very accurately. The temperature dependence of this NTE behavior follows the Planck occupancy factor for phonon modes above 50 meV, corresponding to the O-dominated band of optical frequencies. In the QHA these modes above 50 meV do not contribute to the NTE. These modes are highly anharmonic as shown by their large broadenings and shifts.Simlar to rutileSnO2, because of the large mass difference between Ag and O atoms, the O-dominated phonon modes are well separated from the Ag-dominated modes. Partial phonon DOS analysis showed that the Ag-dominated modes had similar energies, forming the peak of the phonon DOS below 20 meV (peak 1), whereas the O-dominated modes had energies above 40 meV (peak 2 and 3) as shown in the DOS spectrum in Figure 6.Figure 6
The TDOS spectra ofAg2O, D(ω), at 40 K (dashed) and 400 K (solid). The downconversion and upconversion contributions are presented separately as black and green curves, respectively. The phonon DOS spectra are presented along with the TDOS but take a negative sign and shaded blue to make the view of TDOS and DOS separated.For cubic anharmonicity, as discussed, the two-phonon DOS (TDOS) is the spectral quantity parameterizing the number of phonon-phonon interaction channels available to a phonon. ForAg2O with the cuprite structure, the peaks in the TDOS overlap well with the peaks in the phonon DOS, as illustrated in Figure 6. The phonon DOS is analogous to a reservior used as the supply of phonons while TDOS indicates the gates that can direct them. We can therefore understand that most of the phonons have many possible interactions with other phonons, which contributes to the large anharmonicity of Ag2O with the cuprite structure and short lifetimes (or large linewidths). Although the lifetimes of most phonon modes in Ag2O with the cuprite structure are small and similar, the origins of these lifetime broadenings are intrinsically different. For peak 2 of the phonon DOS, the anharmonicity is largely from the upconversion processes: O ↦ O − Ag, while for peak 3 it is from the downconversion processes: O ↦ O + Ag. The anharmonicity of peak 1 is more complicated. It involves both upconversion and downconversion processes of Ag-dominated modes.Owing to explicit anharmonicity from phonon-phonon interactions, the thermodynamic properties ofAg2O with the cuprite structure cannot be understood as a sum of contributions from independent normal modes. The frequency of an anharmonic phonon depends on the level of excitation of other modes. At high temperatures, large vibrational amplitudes increase the anharmonic coupling of modes, and this increases the correlations between the motions of the Ag and O atoms, as shown by perturbation theory. Couplings in perturbation theory have phase coherence, so the coupling between Ag- and O-dominated modes at higher energies, as seen in the peak of the TDOS, causes correlations between the motions of Ag and O atoms. The ab initio MD simulations show that anharmonic interactions allow the structure to become more compact with increasing vibrational amplitude. The mutual motions of the O and Ag atoms cause higher density as the atoms fill space more coherently as a result of anharmonic interactions. The relationship between macroscopic thermal expansion and microscopic anharmonic coupling of different modes requies more quantitative investigation on, for example, a particular decay channel that contributes most to such effect. It should be noted that the large difference in atomic radii of Ag and O may substantially contribute to this effect in collaboration with the coherent anharmonic interaction. For cuprite Cu2O, which has less of a difference in atomic radii, the thermal expansion is much less anomalous although the phonon anharmonicity is found to be equally large.


## 3.1. Computational Methodology

Anharmonicity tensors describe the coupling strengths for phonon-phonon interactions, but a prerequisite is that the phonons in these processes satisfy the kinematical conditions of conservation of energy and momentum as presented in (5). In the phonon-phonon interaction functional, an anharmonicity tensor element for an s-phonon process can be expressed as [39](6)Vj;q→1j1;…;q→s-1js-1=12s!ħ2Ns/2NΔq→1+⋯+q→s-1ωj0ω1⋯ωs-11/2Cj;q→1j1;…;q→s-1js-1,where Δ(q→1+⋯+q→s-1) enforces momentum conservation and the C(·)’s, elements of the s-phonon anharmonic tensor, are expected to be slowly varying functions of their arguments.If the anharmonicity tensor or its average does not vary significantly for different phonon processes, the coupling factor and the kinematic factor are approximately separable in (6). The separation of the anharmonic coupling and the kinematics has been used with success in many studies including recent reports on rutile TiO2 and SnO2 [25, 26]. We consider the term C(j;q→1j1;…;q→s-1js-1) to be a constant of the Raman mode j and use it as a fitting parameter. Although C(j;q→1j1;q→2j2) and C(j;j;q→1j1;-q→1j1) change with j1 and j2, an average over modes, 〈C(·)〉=∑1,2‍C(j;q→1j1;q→2j2)/∑1,2‍1, is needed by the fitting, where 1 and 2 under the summation symbol represent q→iji. We define the cubic and quartic fitting parameters as(7a)Cj3=Cj;q→1j1;q→2j2,(7b)Cj4=Cj;j;q→1j1;-q→1j1.To the leading order of cubic and quartic anharmonicity, the broadening of the Raman peaks is2Γ(3)(j;Ω). The frequency shift of the Raman peaks is ΔQ+Δ(3)+Δ(3′)+Δ(4), where the quasiharmonic part is denoted by ΔQ. These quantities can be written as functions of D(Ω,ω1,ω2) and P(Ω,ω1,ω2), weighted by average anharmonic coupling strengths [25, 26, 40](8a)Γ3j;Ω=πħ64ωj0Cj32∑q→1,j1∑q→2,j2ω1ω2DΩ,ω1,ω2=ωj0Cj32DωΩ,(8b)Δ3j;Ω=-ħ64ωj0Cj32∑q→1,j1∑q→2,j2ω1ω2PΩ,ω1,ω2=ωj0Cj32PωΩ,(8c)Δ3′j=-ħ16Nωj0Cj32∑q→2j2ωj2q→2nq→2j2+12,(8d)Δ4j=ħ8Nωj0Cj4∑q→1j1ωj1q→1nq→1j1+12,where Dω(Ω) and Pω(Ω) are functionals of D(Ω,ω1,ω2) and P(Ω,ω1,ω2) weighted by the kinematics of anharmonic phonon coupling. Dω(Ω) is the so called two-phonon density of states (TDOS) spectra, which characterize the size of the phonon-phonon interaction channels. Δ(3′) is an additional low-order cubic term that corresponds to instantaneous three-phonon processes [32]. It is nonzero for crystals having atoms without inversion symmetry, as in the case for the oxygen atom motions in the A1g mode of rutile. It is much smaller than other contributions, however, owing to symmetry restrictions.


## 3.2. Example  1: Anomalous Temperature-Dependent Broadening ofB2g Mode of Rutile SnO2

Recently, strong anharmonicity of rutileSnO2 was discovered by the high-resolution Raman spectrometer, and the anomaly of the temperature-dependent broadening of B2g mode is a prominent example [26]. As shown in Figure 3(b), at high temperatures above 500 K, the broadening of the B2g mode of rutile SnO2 shows an unusual concave downward curvature, while the other two Raman modes Eg and A1g broaden linearly. Moreoever, at low temperatures, the B2g mode has a much larger linewidth than the other two modes. The linewidth of the B2g mode extrapolated to 0 K is approximately 8 cm-1, whereas the linewidths of the Eg and A1g modes extrapolate to less than 2 cm-1.Figure 3
Temperature dependence of (a) frequency shifts and (b) linewidths as FWHMs, of the Raman modesEg, A1g, and B2g of rutile SnO2. The solid and open symbols represent the experimental data from powder and single crystal samples, respectively. Solid curves are the theoretical fittings with a full calculation of the kinematics of three-phonon and four-phonon processes. The dashed curve was calculated without considering the frequency dependence of D0↓(Ω) and the number of decay channels, at elevated temperatures. Reprinted figure with permission from Lan et al. [26]. Copyright by the American Physical Society.
(a) (b)The TDOS functionDω(Ω) in Figure 4 shows large variations with Ω which explains a trend in the thermal broadening of Figure 3(b). Owing to the high frequency of the B2g mode, at the temperatures of this study its phonon-phonon anharmonicity comes mostly from downconversion processes, as shown in Figure 4. Ignoring the small upconversion contribution, at high temperatures(9)DωΩ,T∝T∑q→1,j1∑q→2,j2δΩ-ω1-ω2≡TD0↓Ω,where D0↓(Ω) is the number of two-phonon downconversion channels. Therefore, the linewidth broadening resulting from the cubic anharmonicity is proportional to the temperature, modulated by C(Ω). Because Ω undergoes a shift with temperature, D0↓(Ω) is an implicit function of T. Usually the line broadening is linear in T because D0↓(Ω) does not vary much with temperature, as in the cases of modes Eg and A1g. However, the B2g mode at 774 cm-1 lies on a steep gradient of Dω(Ω,T) in Figure 4. On the other hand, the B2g mode undergoes a significant shift of frequency by more than 25 cm-1 with temperature as shown in Figure 3(a), which moves the downconversion channels down the gradient of Dω(Ω,T). These cooperative effects as a result narrow the capacity of these decay channels. Therefore, its broadening is significantly less than linear in T, exhibiting an unusual concave downwards shape. For comparison, the dashed line in Figure 3(b) was calculated without considering the frequency dependence of D0↓(Ω) at elevated temperatures, and it deviates substantially from the experimental trend.Figure 4
Two-phonon density of statesDω(Ω) of rutile SnO2 for 0 K and 800 K. The arrowheads mark the positions of the three Raman modes, Eg, A1g, and B2g, respectively. The upconversion and downconversion contributions at 800 K are shown in green and black dashed curves, respectively. There is no upconversion process at 0 K. The phonon DOS spectra are presented along with the TDOS but take a negative sign and shaded blue to make the view of TDOS and DOS separated.As discssed, the unusual temperature dependence of the linewidth of theB2g mode originates from the sharp peak in Dω(Ω,T) centered at 800 cm-1. We shall see that this feature in the TDOS originates with the phonon DOS of SnO2 which has a wide band gap between 360 cm-1 and 450 cm-1 associated with the mass difference of Sn and O atoms, as shown in Figure 4. The shape of the TDOS can be understood as the convolution of the phonon DOS with itself. With two approximately equal regions above and below the gap, the result is a peak at 800 cm-1, with steep slopes on both sides. For comparison, although the TDOS of rutile TiO2 is shaped as a broad peak [25], it does not have the sharp features of Figure 4 because the mass difference between Ti and O atoms does not cause a band gap in the phonon DOS of TiO2.The fine structures of phonon-phonon interaction channels can also explain the physics of large linewidth ofB2g mode at low temperature limit following the same approach as detailed above. In the low temperature limit, upconversion processes are prohibited because n→0. The peak linewidth extrapolated to 0 K is determined entirely by downconversion processes, quantified by the downconversion TDOS. The B2g mode has a significant broadening because its frequency is near a peak in the downconversion TDOS, as seen in Figure 4. On the other hand, the Eg and A1g modes are not broadened at low temperatures because their frequencies are at low values of the TDOS. The study on the structures of phonon-phonon interaction channels successfully marks the “hot” downconversion decay activities even at low temperature limit and accounts for the significant difference of linewidths between B2g and the other two Raman modes at low temperatures.


## 3.3. Example  2: Negative Thermal Expansion of Cuprite Ag2O

Silver oxide (Ag2O) with the cuprite structure has attracted much interest after the discovery of its extraordinarily large negative thermal expansion (NTE) [41, 42], which exceeds -1×10-5 K-1 and occurs over a wide range of temperature from 40 K to its decomposition temperature near 500 K.A rigid-unit modes (RUMs) model of NTE considers tetrahedra ofAg4O around each O atom that bend at the Ag atoms linking the O atoms in adjacent tetrahedra. RUMs account for counteracting rotations of all such tetrahedra [1, 2]. Locally, the O–Ag bond length does not contract, but bending of the O–Ag–O links pulls the O atoms together, leading to NTE. These RUMs tend to have low frequencies owing to the large mass of the unit and hence are excited at low temperatures. This model correlated the NTE with quasiharmonic approximation and should explain the main physics at low temperatures. However, as shown in Figure 5, at temperatures above 250 K, there is a second part of the NTE behavior of cuprite Ag2O that is apparently beyond the predictions of quasiharmonic theory [9].Figure 5
Temperature dependence of lattice parameter ofAg2O from experimental data in [41], quasiharmonic calculations, and MD calculations, expressed as the relative changes with respect to their 40 K values, that is, a(T)/a(40K)-1.First-principles Born-Oppenheimer molecular dynamics simulations were performed for a3×3×3 supercell with temperature control by a Nosé thermostat. The simulated temperatures included 40, 100, 200, 300, and 400 K. For each temperature, the system was first equilibrated for 3 ps and then simulated for 18 ps with a time step of 3 fs. The system was fully relaxed at each temperature, with convergence of the pressure within 1 kbar.As shown in Figure5, the MD simulation predicts the NTE very accurately. The temperature dependence of this NTE behavior follows the Planck occupancy factor for phonon modes above 50 meV, corresponding to the O-dominated band of optical frequencies. In the QHA these modes above 50 meV do not contribute to the NTE. These modes are highly anharmonic as shown by their large broadenings and shifts.Simlar to rutileSnO2, because of the large mass difference between Ag and O atoms, the O-dominated phonon modes are well separated from the Ag-dominated modes. Partial phonon DOS analysis showed that the Ag-dominated modes had similar energies, forming the peak of the phonon DOS below 20 meV (peak 1), whereas the O-dominated modes had energies above 40 meV (peak 2 and 3) as shown in the DOS spectrum in Figure 6.Figure 6
The TDOS spectra ofAg2O, D(ω), at 40 K (dashed) and 400 K (solid). The downconversion and upconversion contributions are presented separately as black and green curves, respectively. The phonon DOS spectra are presented along with the TDOS but take a negative sign and shaded blue to make the view of TDOS and DOS separated.For cubic anharmonicity, as discussed, the two-phonon DOS (TDOS) is the spectral quantity parameterizing the number of phonon-phonon interaction channels available to a phonon. ForAg2O with the cuprite structure, the peaks in the TDOS overlap well with the peaks in the phonon DOS, as illustrated in Figure 6. The phonon DOS is analogous to a reservior used as the supply of phonons while TDOS indicates the gates that can direct them. We can therefore understand that most of the phonons have many possible interactions with other phonons, which contributes to the large anharmonicity of Ag2O with the cuprite structure and short lifetimes (or large linewidths). Although the lifetimes of most phonon modes in Ag2O with the cuprite structure are small and similar, the origins of these lifetime broadenings are intrinsically different. For peak 2 of the phonon DOS, the anharmonicity is largely from the upconversion processes: O ↦ O − Ag, while for peak 3 it is from the downconversion processes: O ↦ O + Ag. The anharmonicity of peak 1 is more complicated. It involves both upconversion and downconversion processes of Ag-dominated modes.Owing to explicit anharmonicity from phonon-phonon interactions, the thermodynamic properties ofAg2O with the cuprite structure cannot be understood as a sum of contributions from independent normal modes. The frequency of an anharmonic phonon depends on the level of excitation of other modes. At high temperatures, large vibrational amplitudes increase the anharmonic coupling of modes, and this increases the correlations between the motions of the Ag and O atoms, as shown by perturbation theory. Couplings in perturbation theory have phase coherence, so the coupling between Ag- and O-dominated modes at higher energies, as seen in the peak of the TDOS, causes correlations between the motions of Ag and O atoms. The ab initio MD simulations show that anharmonic interactions allow the structure to become more compact with increasing vibrational amplitude. The mutual motions of the O and Ag atoms cause higher density as the atoms fill space more coherently as a result of anharmonic interactions. The relationship between macroscopic thermal expansion and microscopic anharmonic coupling of different modes requies more quantitative investigation on, for example, a particular decay channel that contributes most to such effect. It should be noted that the large difference in atomic radii of Ag and O may substantially contribute to this effect in collaboration with the coherent anharmonic interaction. For cuprite Cu2O, which has less of a difference in atomic radii, the thermal expansion is much less anomalous although the phonon anharmonicity is found to be equally large.


## 4. Conclusions

Today, our understanding about the vibrational thermodynamics of materials at low temperatures is broad and deep because it is based on the harmonic model in which phonons are independent, avoiding issues of anharmonic lattice dynamics and phonon-phonon interactions. However, in most cases, the failure of the harmonic theory also arises from the assumption of independent phonons which becomes increasingly inaccurate at high temperatures. We have seen that the anharmonic renormalization of the vibrational quanta could break the independent and harmonic phonon assumptions and substantially change their characteristics in a coherent way.The complexity of phonon anharmonicity arises because we need to consider how phonons interact with other phonons or with other excitations, which is an example of notorious many-body interaction problem. To this end, progress of anharmonic phonon theories and advanced computational methodologies discussed in this review provide both the macroscopic and the microscopic perspectives of phonon anharmonicty and its relationship with thermodynamic properties, such as high temperature vibrational energy distribution, thermal stability, and anomalous thermal expansion.Based upon quantum perturbation theory and first-principles molecular dynamics simulations, microstructures of renormalized phonons can be established. In particular, the effective potential energy surface as well as the anharmonic energy spectra of renormalized phonons can be built up by analyzing the atomic vibration trajectories derived from MD. Incorporated with frozen phonon method, we have seen that these methods give rise to the anharmonic potential landscape fully accounting for the phonon-phonon interactions. The evolution of the potential surface with temperature as well as behaviors of different orders of anharmoncity can be derived in detail as well. On the other hand, the fine structure of the phonon-phonon interaction channels by means of two-phonon density of states provides unique insight into the microscopic picture of pathways of phonon renormalization processes. With phonon density of states, these fine structures give crucial information about how phonons interact with each other and the distribution of three-phonon and four-phonon kinematics in these interactions.



---

*Source: 2714592-2016-10-19.xml*