Phonon quarticity induced by changes in phonon-tracked hybridization during lattice expansion and its stabilization of rutile TiO{}_{2}

# Phonon quarticity induced by changes in phonon-tracked hybridization during lattice expansion and its stabilization of rutile TiO$_2$

## Abstract

Although the rutile structure of TiO is stable at high temperatures, the conventional quasiharmonic approximation predicts that several acoustic phonons decrease anomalously to zero frequency with thermal expansion, incorrectly predicting a structural collapse at temperatures well below 1000 K. Inelastic neutron scattering was used to measure the temperature dependence of the phonon density of states (DOS) of rutile TiO from 300 to 1373 K. Surprisingly, these anomalous acoustic phonons were found to increase in frequency with temperature. First-principles calculations showed that with lattice expansion, the potentials for the anomalous acoustic phonons transform from quadratic to quartic, stabilizing the rutile phase at high temperatures. In these modes, the vibrational displacements of adjacent Ti and O atoms cause variations in hybridization of electrons of Ti and electrons of O atoms. With thermal expansion, the energy variation in this “phonon-tracked hybridization” flattens the bottom of the interatomic potential well between Ti and O atoms, and induces a quarticity in the phonon potential.

## I Introduction

Titanium dioxide (TiO) is of longstanding interest in physics, chemistry, surface science and materials science. It has numerous technological applications including pigments, optical coatings, catalysis, solar cells and gas sensors.Linsebigler et al. (1995); Diebold (2003); Du et al. (2009); Lu and Z. (2009) Rutile is the stable phase of TiO at temperatures below 1800 K, and the two other naturally-occurring phases of TiO, anatase and brookite, both convert to rutile upon heating.Isaak et al. (1998) The thermal properties of rutile TiO are central for applications that involve heat generation, thermal transport and temperature-driven phase transitions. Nevertheless, even after numerous experimental Smith et al. (2009); Zhang and Banfield (1998); Traylor et al. (1971); Samara and Peercy (1973); Burdett et al. (1987) and theoretical Lee et al. (1994); Vast et al. (2002); Muscat et al. (2002); Montanari and Harrison (2004); Labat et al. (2007); Mitev et al. (2010); Refson et al. (2013); Sikora (2005); Lan et al. (2012a) investigations, the stability of the rutile phase at high temperatures remains poorly understood.

The quasiharmonic approximation (QHA) is based on how phonon frequencies change with volume. In the QHA, all shifts of phonon frequencies from their low temperature values are the result of thermal expansion alone.Fultz (2010); Li et al. (2015) For a lattice expansion of 0.5%, the energies of transverse phonons decrease by 10% to 50% in the QHA, giving Grüneisen parameters as large as 100 for rutile TiO.Mitev et al. (2010); Refson et al. (2013) In the QHA, many transverse acoustic modes decrease to zero frequency at a lattice expansion of 0.75%, so the QHA predicts a collapse of the rutile structure at temperatures below 1000  K. The QHA does not account for all non-harmonic effects, however. Although the QHA accounts for some frequency shifts, the phonon modes are still assumed to be harmonic, non-interacting, and their energies depend only on the volume of the crystal.

In non-harmonic potentials, phonon-phonon interactions are responsible for pure anharmonicity that shortens phonon lifetimes and shifts phonon frequencies. 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 Ye et al. (1987) and the possible stabilization of bcc Fe-Ni alloys at conditions of the Earth’s core.Dubrovinsky et al. (2007) For PbTe and ScF, there are recent reports of anharmonicity being so large that both the QHA and anharmonic perturbation theory fail dramatically.Delaire et al. (2011); Li et al. (2011) These cases are suitable for ab initio molecular dynamics (AIMD) simulations, however, which naturally includes all orders of anharmonicity.de Koker (2009); Thomas et al. (2010); Lan et al. (2014) The AIMD method should be reliable when the electrons are near their ground states and the nuclear motions are classical.

In this work, we study phonons at high temperatures in rutile TiO by neutron inelastic scattering, AIMD, and other methods. We report a remarkable quartic anharmonicity that develops because the hybridization between electron states at Ti and O atoms changes dynamically during atom vibrations, and these changes are highly sensitive to lattice expansion. This “phonon-tracked hybridization” induces the phonon quarticity that is essential for stabilizing the rutile phase at high temperatures, and affects properties such as ferroelectricity and thermal transport.

## Ii Experiment

### ii.1 Inelastic neutron scattering

Samples were commercial TiO powder (Alfa Aesar, Ward Hill, MA) with a rutile phase fraction of at least 99.9%. The sample powder was loaded into a niobium sachet for neutron scattering measurements at high temperatures. Measurements were performed with the ARCS spectrometer at the Oak Ridge National Laboratory Abernathy et al. (2012), using incident neutron energies of 75 meV and 30 meV. A vacuum furnace with niobium radiation shields was used for measurements at 300, 673, 1073, and 1373 K. Furnace backgrounds were measured with empty sample containers at each temperature. Data reduction was performed with the software package DRCS.DRCS () The raw data of individual neutron detection events were first binned into scattering angle and energy transfer, , and normalized by the proton current on target. The data were corrected for furnace background and detector efficiency. The data were then rebinned into intensity, , where is the momentum transfer to the sample. This is averaged over all crystallographic directions, and is a reasonably well-balanced sampling of all phonons in the material.

### ii.2 Results

Surprisingly, temperature causes a global increase in energy of the entire TA branch. The lowest phonon dispersions formed by the TA branch at 14 meV at 300 K increase significantly in energy at 1373 K, as seen in the phonon spectra (Fig. 1). There are no dispersions that soften between 2 and 4 , and there are no signs of phonons that collapse to zero frequency. This anomalous thermal stiffening is also apparent in the phonon DOS curves, which were obtained after corrections for multiphonon and multiple scattering using the getdos package.Kresch et al. (2008) As shown in Fig. 2(a), peak 1 (centered at 14 meV) exhibits an unusual stiffening as large as 2.7 meV, and it remains sharp. The large stiffening of this DOS peak is even more apparent in measurements with higher energy resolution, using an incident energy of 30 meV (see Fig. 3(a)). On the other hand, as shown in Fig. 2(a) and Fig. 3(a), most other features of the phonon DOS above 20 meV undergo substantial softening and broadening with temperature, in good agreement with the individual phonons measured by Raman spectrometry at high temperatures.Samara and Peercy (1973); Lan et al. (2012a) It is well-known that rutile TiO has a ferroelectric soft mode A at the -point that stiffens with heating.Samara and Peercy (1973) but the global stiffening of the entire lower-lying TA branch is surprising.

Figure 3(b) presents the DOS calculated with the quasiharmonic model at 1373 K, showing that contrary to the experimental trend in Fig. 3(a), the first peak in the phonon DOS undergoes a large softening with temperature. Consistent with previous results,Mitev et al. (2010); Refson et al. (2013) our QHA calculations predict that modes in the TA branch soften through zero frequency with the thermal expansion of 1373 K, giving imaginary frequencies in the DOS (shown as negative in Fig. 3(b)) that would destabilize the rutile structure at high temperatures.

## Iii Calculations

### iii.1 First principles molecular dynamics simulations and techniques for calculating the vibrational spectra

First-principles calculations using the local density approximation (LDA) of density functional theory (DFT) Perdew and Zunger (1981) were performed with the VASP package.Kresse and Joubert (1999); Kresse and Furthmüller (1996) Projector augmented wave pseudopotentials and a plane wave basis set with an energy cutoff of 500 eV were used in all calculations. Previous work showed that for best accuracy, the Ti pseudopotential should treat the semicore 3 and 3 states as valence,Montanari and Harrison (2004); Mitev et al. (2010); Refson et al. (2013) and we took this approach with a similar LDA functional. Our calculated elastic properties, lattice dynamics, and dielectric properties derived from the optimized structure for 0 K, were in good agreement with results from experiment and from previous DFT calculations. Phonon dispersions and phonon DOS spectra were obtained in the quasiharmonic approximation by minimizing the vibrational free energy as a function of volume

 F(V,T)=E0++∞∫−∞g(ω)(ℏω2+kBTln(1−e−ℏω/kBT))dω, (1)

where is the energy calculated for the relaxed structure at K.

First-principles Born-Oppenheimer AIMD simulations for a supercell and a -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 ensemble with temperature control by a Nosé thermostat, then simulated as an ensemble for 20 ps with time steps of 1 fs. Good relaxations with residual pressures below 0.5 GPa were achieved in each calculation that accounted for thermal expansion. The lattice parameters and thermal expansivities determined by the simulations are presented in Table 1.

#### Fourier transformed velocity autocorrelation method

Velocity trajectories were extracted from the MD simulation at each temperature, and were then transformed to the corresponding vibrational energy and/or momentum domain.de Koker (2009); Thomas et al. (2010); Lan et al. (2012a, 2014); the () Because the FTVAC method does not assume a form for the Hamiltonian, it is a robust tool for obtaining vibrational spectra from MD simulations, even with strong anharmonicity. The phonon DOS is given by

 g(\omegaup)=∑n,b∫e−i\omegaupt⟨→vn,b(t)→v0,0(0)⟩dt, (2)

where is an ensemble average, and is the velocity of the atom in the unit cell at time . Further projection of the phonon modes onto each 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.

#### 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.Maradudin and Fein (1962); Lan et al. (2012b) To distinguish the roles of cubic and quartic anharmonicity, the temperature-dependent effective potential (TDEP) method Hellman et al. (2011); Hellman and Abrikosov (2013) was used. In the TDEP method, an effective model Hamiltonian 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 Hellman et al. (2011)

 H = U0+12∑imp2i+12∑ijαβϕαβijuαiuβj (3) +13!∑ijkαβγψαβγijkuαiuβjuγk,

where and are second- and third-order force constants, is momentum, and is the Cartesian component of the displacement of atom . In the fitting, the “effective” harmonic force constants are renormalized by the quartic anharmonicity. The cubic anharmonicity, however, is largely accounted for by the third-order force constants , and can be understood in terms of the third-order phonon self-energy that causes linewidth broadening.Maradudin and Fein (1962); Li et al. (2014)

The resulting Hamiltonian was used to obtain the renormalized phonon dispersions (TDEP spectra) accounting for both the anharmonic shifts , and broadenings , of the mode . These are derived from the real and imaginary parts of the cubic self-energies , respectively.Maradudin and Fein (1962)

 Δ(→qj;Ω) = −18ℏ2∑→q1j1∑→q2j2∣∣V(→qj;→q1j1;→q2j2)∣∣2Δ(→q1+→q2−→q) (4) ×℘[n1+n2+1Ω+\omegaup1+\omegaup2−n1+n2+1Ω−\omegaup1−\omegaup2 +n1−n2Ω−\omegaup1+\omegaup2−n1−n2Ω+\omegaup1−\omegaup2] Γ(→qj;Ω) = 18πℏ2∑→q1j1∑→q2j2∣∣V(→qj;→q1j1;→q2j2)∣∣2Δ(→q1+→q2−→q) (5) ×[(n1+n2+1)δ(Ω−\omegaup1−\omegaup2) +2(n1−n2)δ(Ω+\omegaup1−\omegaup2)],

where is the renormalized phonon frequency and denotes the Cauchy principal part. The ’s are elements of the Fourier transformed third order force constants obtained in the TDEP method. The ensures conservation of momentum.

### iii.2 Crystal orbital Hamilton population method

The Crystal Orbital Hamilton Population Method (COHP) projects plane waves to local orbital basis functions and is used to partition the band-structure energy into the bonding, nonbonding and antibonding interactions between two atoms. COHP provides a semiquantitative interpretation of the net bonding characteristics.

In this work, COHP spectra were calculated with the package Local-Orbital Basis Suite Towards Electronic-Structure Reconstruction (LOBSTER).Dronskowski and Blchl (1993); Maintz et al. (2013) The phonon displacement patterns were generated in a supercell with a -point sampling according to the eigenvectors of the mode. Self-consistent wavefunctions were then computed and projected onto a local basis of Slater type orbitals with proper orthonormalization, giving the COHP spectra of the system. Since negative values of the COHP suggest bonding contributions, – COHP spectra are presented, following convention.

### iii.3 Results

As shown in Fig. 2(b), the neutron-weighted DOS spectra calculated with FTVAC method reproduced the thermal shifts and broadenings observed in the experimental DOS of Fig. 2(a). In particular, upon heating to 1373 K, an anomalous stiffening of peak 1 by 2.1 meV was obtained.

Figure 4 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. Especially for this TA branch, Fig. 4(b) shows an enormous discrepancy of phonon energies between the MD calculation and the QHA (orange dashed line) at 1373 K. The unstable phonon modes predicted by the QHA are fully stable in the AIMD simulations at high temperatures, however.

Using the same MD trajectories as for the FTVAC method, the calculated TDEP dispersions agree well with the FTVAC results as shown in Fig. 4(a)(b), and also in Table 2. (The full TDEP dispersions are shown in the Supplemental Material.sup ()) The cubic anharmonicity of rutile TiO is strong for phonons at energies above 25 meV, Lan et al. (2012a) causing broadening of the phonon DOS of Fig. 2 and broadening of the dispersions.sup () Nevertheless, at 1373 K the TA modes below 20 meV have only small linewidth broadenings. Furthermore, they are close in energy to those calculated if all are set to zero in Eq. 3, showing the dominance of quartic anharmonicity and the small cubic anharmonicity of the TA modes.

## Iv Discussion

### iv.1 Phonon quarticity of TA branch

For more details about the anomalous anharmonicity of the TA modes, the frozen phonon method was used to calculate potential energy surfaces for specific phonons, as in Fig. 4(c). At 300 K the potential energy of the TA mode at the -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 (see Supplemental Material sup ()). Even more interesting behavior occurs along the directions of , and , where the potential energy surface develops a Landau-type form, with negative curvature at zero displacement, but a positive quartic shape at large displacements. (The height of the energy barrier between the two minima is lower than , however.) For a quantum quartic oscillator, the vibrational frequency stiffens with temperature owing to the increasing spread between the energy levels.Dorey and Tateo (1999); Li et al. (2011) 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 Fig. 4(b), they are reasonably close to the FTVAC and TDEP results.

### iv.2 Phonon-tracked hybridization and its change with thermal expansion

The patterns of atom displacements in the anomalous modes at , R, and along Z-, -M, and M-A were identified, and those for the R point in Fig. 5(a) are typical. In these anomalous modes, the O atoms were approximately stationary, and each O atom has a Ti neighbor that moves towards it and another Ti neighbor that moves away from it by approximately the same amount. These modes have “ring” patterns in which displacements of Ti atoms rotate a structural unit, and all the O atoms see approximately the same change in their Ti neighborhood. In the positive and negative displacements of these modes, the O atoms show an accumulation of charge in the Ti–O bond of shorter distance and a depletion in the bond of longer distance, as indicated by a much higher value of the electron localization function (ELF) Silvi and Savin (1994) for the short bond shown in Fig. 5(b).

We calculated the “bond-weighted” electron DOS by partitioning the band structure energy into bonding and no-bonding contributions and obtaining the crystal orbital Hamilton population (COHP) spectrum Dronskowski and Blchl (1993); Maintz et al. (2013) of rutile with different lattice parameters in a supercell. Figure 5(c) shows the COHP spectrum of the bonds formed by the Ti   and O  orbitals between 5.5 and 3.5 eV below the Fermi energy. With lattice expansion (of 1% or 2%), these bonding states become less favorable, and their COHP spread narrows. Also shown in color in Fig. 5(b) is the COHP with a frozen phonon mode at the R-point having 0.14  displacements of Ti atoms. On the scale of thermal energies, the broadening effect from the phonon changes considerably with lattice expansion. For the equilibrium lattice parameter (labeled “0%”), bonding states at both the top and bottom of the Ti   – O  band are shifted upwards by the phonon. With larger lattice parameter, however, the changes take place mostly in the states at the bottom of the band, which broadens the COHP spectra and helps to lower the bonding energy. Additional bonding states also appear near –3 eV. The overall effect is that the potential energy is relatively insensitive to small displacements, and the phonon potential has a much flatter bottom with 1% lattice expansion as shown in Fig. 4(c). For a 2% linear expansion, a Landau-type phonon potential evolves, with local minima away from zero displacement owing to the spread of states to the lowest part of the band.

The effect can be understood by adapting a hard sphere model for the chain of -Ti-O-Ti-O- atoms along the direction “A” or around the “ring” units in the phonon displacement pattern in Fig. 5(a). With an increase in lattice parameter, the longer Ti–O bond makes a smaller contribution to the interatomic force during its vibrational cycle. The shorter Ti–O bond gives a stronger hybridization of Ti   – O  orbitals as the Ti atom moves closer to the O atom. The hybridization serves to offset the energy of short-range repulsion. 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 “ring” patterns of the phonons play an important role in increasing the degree of hybridization as they complete the electron back-donation cycles from the O to the Ti atoms. For example, a 3% decrease of Bader effective charge (+2.22 at equilibrium) was found for the Ti atoms with 0.14   displacments in the “ring” patterns, which is comparable to the charge decrease of the Ti atoms during the ferroelectric transition of rutile.Montanari and Harrison (2004) However, if Ti atoms along the direction “A” were locked down at their equilibrium positions so the ring motion is broken, the resulting decrease of the effective charge dropped by 50% to 70% in the “ring” patterns. The potential was found to rise, and was largely quadratic even at 1% or 2% expansion, as shown in Fig. 6.

A macroscopic elastic response to this phonon can also be identified with the assistance of Fig. 5(a). In equilibrium, the apex angles of the rhombuses are all equal, but with the rotation by , the vertical stretching of rhombuses along the line A is , and the contraction is , where is the semi-angle of the rhombus. For small , a Taylor expansion gives a net vertical (or horizontal) distortion of (or ). The distortions are proportional to , while the atom displacements in this TA mode are proportional to . A strain energy that goes as the square of this distortion is consistent with a quartic potential.

The hybridization in the Ti–O bond is very sensitive to interatomic distance, much as has been noticed in the ferroelectric distortion of BaTiO.Cohen (1992) For rutile TiO, however, the hybridization follows the atom displacements in thermal phonons (instead of a displacive phase transition), and this “phonon-tracked hybridization” changes with lattice parameter. It provides a source of extreme phonon anharmonicity, but also provides thermodynamic stability for rutile TiO. It may occur in other transition metal oxides that show unusual changes of properties with lattice parameter or with structure, and such materials may be tunable with composition or pressure to control this effect. Besides altering thermodynamic phase stability, properties such as ferroelectricity and thermal transport will be affected directly.

## V Conclusion

Several inter-dependent methods proved useful for assessing the anharmonic phonon dynamics of rutile TiO. Ab initio methods of molecular dynamics, temperature-dependent effective potential, and frozen phonons gave accurate accounts of the thermal stiffening of the TA phonons that was measured by inelastic neutron scattering on rutile TiO. These methods avoided the structural instability predicted by the quasiharmonic approximation (QHA). With thermal expansion, many acoustic modes were unstable in the QHA, as observed previously. From AIMD simulations and frozen phonon calculations, it was found that when the rutile structure undergoes thermal expansion, a hybridization between Ti and O atoms gives low-energy electron states when the Ti and O atoms are in closest proximity during their vibrations. This phonon-tracked hybridization flattens the bottom of the potential of the anomalous phonon modes, giving quartic potentials that stabilize the rutile structure at high temperatures.

###### Acknowledgements.
Research at the SNS at the Oak Ridge National Laboratory was sponsored by the Scientific User Facilities Division, BES, DOE. This work was supported by the DOE Office of Science, BES, under contract DE-FG02-03ER46055.

Supplemental Materials:

Phonon quarticity induced by changes in phonon-tracked hybridization during lattice expansion and its stabilization of rutile TiO

### Footnotes

1. Experimental data from Ref. [Traylor et al., 1971]
2. First principles calculations from Ref. [Sikora, 2005]
3. First principles calculations from Ref. [Mitev et al., 2010]
4. First principles calculations from Ref. [Mitev et al., 2010] with ambient temperature but 1% isotropic lattice expansion

### References

1. A. Linsebigler, G. Lu,  and J. T. Yates, Chem. Rev. 95, 735 (1995).
2. U. Diebold, Surf. Sci. Rep. 48, 53 (2003).
3. Y. Du, N. A. Deskins, Z. Zhang, Z. Dohnálek, M. Dupuis,  and I. Lyubinetsky, Phys. Rev. Lett. 102, 096102 (2009).
4. C. Lu and C. Z., Sens. Act. B 140, 109 (2009).
5. D. G. Isaak, J. D. Carnes, O. L. Anderson, H. Cynn,  and E. Hake, Phys. Chem. Minerals 26, 31 (1998).
6. S. J. Smith, R. Stevens, S. Liu, G. Li, A. Navrotsky, B.-G. J.,  and W. B. F., Am. Miner. 94, 236 (2009).
7. H. Zhang and J. F. Banfield, J. Mater. Chem. 8, 2073 (1998).
8. J. G. Traylor, H. G. Smith, R. M. Nicklow,  and M. K. Wilkinson, Phys. Rev. B 3, 3457 (1971).
9. G. A. Samara and P. S. Peercy, Phys. Rev. B 7, 1131 (1973).
10. J. K. Burdett, T. Hughbanks, G. J. Miller, J. W. Richardson, Jr.,  and J. V. Smith, J. Am. Chem. Soc. 109, 3639 (1987).
11. C. Lee, P. Ghosez,  and X. Gonze, Phys. Rev. B 50, 13379 (1994).
12. N. Vast, L. Reining, V. Olevano, P. Schattschneider,  and B. Jouffrey, Phys. Rev. Lett. 88, 037601 (2002).
13. J. Muscat, V. Swamy,  and H. N. M., Phys. Rev. B 65, 224112 (2002).
14. B. Montanari and N. M. Harrison, J. Phys.: Condens. Matter 16, 273 (2004).
15. F. Labat, P. Baranek, C. Domain, C. Minot,  and C. Adamo, J. Chem. Phys. 126, 154703 (2007).
16. P. D. Mitev, K. Hermansson, B. Montanari,  and K. Refson, Phys. Rev. B 81, 134303 (2010).
17. K. Refson, B. Montanari, P. D. Mitev, K. Hermansson,  and N. M. Harrison, Phys. Rev. B 88, 136101 (2013).
18. R. Sikora, J. Phys. Chem. Solids 66, 1069 (2005).
19. T. Lan, X. Tang,  and B. Fultz, Phys. Rev. B 85, 094305 (2012a).
20. B. Fultz, Prog. Mater. Sci. 55, 247 (2010).
21. C. W. Li, H. L. Smith, T. Lan, J. L. Niedziela, J. A. Muñoz, J. B. Keith, L. Mauger, D. L. Abernathy,  and B. Fultz, Phys. Rev. B 91, 144302 (2015).
22. Y. Ye, Y. Chen, K. Ho, B. N. Harmon,  and P. Lindgard, Phys. Rev. Lett. 58, 1769 (1987).
23. L. Dubrovinsky, N. Dubrovinskaia, O. Narygina, I. Kantor, A. Kuznetzov, V. B. Prakapenka, L. Vitos, B. Johansson, A. S. Mikhaylushkin, S. I. Simak,  and I. A. Abrikosov, Science 316, 1880 (2007).
24. O. Delaire, J. Ma, K. Marty, A. F. May, M. A. McGuire, M.-H. Du, D. J. Singh, A. Podlersnyak, G. Ehlers, M. D. Lumsden,  and B. C. Sales, Nature Materials 10, 614 (2011).
25. C. W. Li, X. Tang, J. A. Muñoz, J. B. Keith, S. J. Tracy, D. L. Abernathy,  and B. Fultz, Phys. Rev. Lett. 107, 195504 (2011).
26. N. de Koker, Phys. Rev. Lett. 103, 125902 (2009).
27. J. A. Thomas, J. E. Turney, R. M. Iutzi, C. H. Amon,  and A. J. H. McGaughey, Phys. Rev. B 81, 081411 (2010).
28. T. Lan, C. W. Li, J. L. Niedziela, H. Smith, D. L. Abernathy, G. R. Rossman,  and B. Fultz, Phys. Rev. B 89, 054306 (2014).
29. D. L. Abernathy, M. B. Stone, M. J. Loguillo, M. S. Lucas, O. Delaire, X. Tang, J. Y. Y. Lin,  and B. Fultz, Rev. Sci. Instum. 83, 015114 (2012).
30. DRCS, http://danse.us/trac/DrChops .
31. M. Kresch, M. Lucas, O. Delaire, J. Y. Lin,  and B. Fultz, Phys. Rev. B 77, 024301 (2008).
32. J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
33. G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
34. G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
35. T. Lan, Studies of phonon anharmonicity in solids, Ph.D. thesis, California Institute of Technology (2014) .
36. A. A. Maradudin and A. E. Fein, Phys. Rev. 128, 2589 (1962).
37. T. Lan, C. W. Li,  and B. Fultz, Phys. Rev. B 86, 134302 (2012b).
38. O. Hellman, I. A. Abrikosov,  and S. I. Simak, Phys. Rev. B 84, 180301 (2011).
39. O. Hellman and I. A. Abrikosov, Phys. Rev. B 88, 144301 (2013).
40. C. W. Li, O. Hellman, J. Ma, A. F. May, H. B. Cao, X. Chen, A. D. Christianson, G. Ehlers, D. J. Singh, B. C. Sales,  and O. Delaire, Phys. Rev. Lett. 112, 175501 (2014).
41. R. Dronskowski and P. E. Blchl, J. Phys. Chem. 97, 8617 (1993).
42. S. Maintz, V. L. Deringer, A. L. Tchougreeff,  and R. Dronskowski, J. Comput. Chem. 34, 2557 (2013).
43. See Supplemental Material at [URL will be inserted by publisher] .
44. P. Dorey and R. Tateo, J. Phys. A 32, L419 (1999).
45. B. Silvi and A. Savin, Nature 371, 683 (1994).
46. R. E. Cohen, Nature 358, 136 (1992).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
100422

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