# Nuclear quantum effect with pure anharmonicity and the anomalous thermal expansion of silicon

^{a}Department of Applied Physics and Materials Science, California Institute of Technology, Pasadena, CA 91125;^{b}Neutron Data Analysis and Visualization Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831;^{c}Department of Mechanical and Civil Engineering, California Institute of Technology, Pasadena, CA 91125;^{d}Instrument and Source Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831;^{e}Department of Mechanical Engineering, University of California, Riverside, CA 92521;^{f}Quantum Condensed Matter Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831

See allHide authors and affiliations

Edited by Alexi Maradudin, University of California, Irvine, CA and accepted by Editorial Board Member Zachary Fisk December 26, 2017 (received for review May 18, 2017)

## Significance

Silicon has a peculiar negative thermal expansion at low temperature. This behavior has been understood with a “quasiharmonic” theory where low-energy phonons decrease in frequency with volume contraction. We report inelastic neutron scattering measurements of phonon dispersions over a wide range of temperatures. These measurements cast doubt upon quasiharmonic theory, which predicts the wrong sign for most phonon shifts with temperature. Fully anharmonic ab initio calculations correctly predict the phonon shifts and thermal expansion. Crystal structure, anharmonicity, and nuclear quantum effects all play important roles in the thermal expansion of silicon, and a simple mechanical explanation is inappropriate. The quantum effect of nuclear vibrations is also expected to be important for thermophysical properties of many materials.

## Abstract

Despite the widespread use of silicon in modern technology, its peculiar thermal expansion is not well understood. Adapting harmonic phonons to the specific volume at temperature, the quasiharmonic approximation, has become accepted for simulating the thermal expansion, but has given ambiguous interpretations for microscopic mechanisms. To test atomistic mechanisms, we performed inelastic neutron scattering experiments from 100 K to 1,500 K on a single crystal of silicon to measure the changes in phonon frequencies. Our state-of-the-art ab initio calculations, which fully account for phonon anharmonicity and nuclear quantum effects, reproduced the measured shifts of individual phonons with temperature, whereas quasiharmonic shifts were mostly of the wrong sign. Surprisingly, the accepted quasiharmonic model was found to predict the thermal expansion owing to a large cancellation of contributions from individual phonons.

A quantized harmonic oscillator was Einstein’s seminal idea for understanding atom vibrations in solids. Better accuracy for crystalline solids is achieved when the vibrations are resolved into normal modes. Each normal mode is quantized, with a zero-point energy and excitations called phonons. However, harmonic models are limited to quadratic terms in the interatomic potential, and it is well known that higher-order terms are necessary to describe properties of real solids such as thermal conductivity and thermal expansivity. Despite this knowledge, the necessary and sufficient contributions to nonharmonic effects remain less clear. A popular approach is the quasiharmonic model (QH), which assumes harmonic oscillators, but with frequencies renormalized to account for the thermal expansion. In a QH, the energy of the phonon mode i changes with crystal volume, V. Changes to phonon energies are usually described by a mode Grüneisen parameter,

The cubic and quartic, and higher-order terms of the interatomic potential, cause the normal modes to interact and exchange energy. This is pure anharmonicity, where the energy of a phonon is altered by the presence of other phonons irrespective of the volume of the solid. Phonon anharmonicity is essential for thermal conductivity and other thermophysical properties. Anharmonic effects increase with larger thermal atomic displacements. Sometimes this causes a misperception that pure anharmonicity is important only at high temperatures, and quasiharmonic models may be valid at low and moderate temperatures owing to low phonon populations. However, the leading-order terms of both quasiharmonicity and anharmonicity are linear in temperature (4), so, if anharmonicity is important at high temperatures, it can have the same relative importance at low temperatures, too. Furthermore, at low temperatures, the “zero-point” energy gives atom displacements that allow a nuclear quantum effect to engage the high-energy phonon modes that are half-occupied.

Finding the relative importances of quasiharmonicity and anharmonicity should be done by quantitative analysis, but, to date, the dominance of quasiharmonicity for silicon has been assumed, in part, because quasiharmonic models predict the thermal expansion with reasonable accuracy (5⇓⇓⇓⇓⇓–11). The quasiharmonic model predicts the anomalous negative thermal expansion of silicon from 10 K to 125 K and predicts the low thermal expansion up to the melting temperature (12⇓⇓⇓–16). The positive thermal expansion coefficients observed at moderate and high temperatures are anomalous in their own right—they are small compared with diamond and other materials with zincblende structures (14). Further validation of the quasiharmonic approximation was provided by measurements of the Raman mode and a few second-order Raman modes of silicon under pressure, which were accurately predicted by volume-dependent density functional theory (DFT) calculations at low temperature (17, 18). The negative Grüneisen parameters of the low-energy transverse acoustic (TA) modes have received considerable attention and have been attributed to the “openness” of the diamond cubic structure (16), the stability of angular forces (9), or entropy in general (8). Nevertheless, the precise role of the TA modes in thermal expansion remains unclear (7, 9). With increasing temperature, phonons are excited in higher-energy phonon branches, and their positive Grüneisen parameters are expected to cause the overall thermal expansion to change sign. Today, this quasiharmonic model is the workhorse for predicting thermal expansion.

“Nontrivial” phonon shifts that were not accounted for by thermal expansion were reported in an earlier experimental paper on phonon dispersions in silicon up to 300 K (3). The importance of pure anharmonicity in temperature-dependent phonon shifts at moderate and high temperatures was also found in work based on molecular dynamics, many-body perturbation theory, and ab initio calculations on silicon (19⇓⇓⇓⇓⇓⇓⇓–27). The uncertainty principle and quantum distributions of nuclear positions influence the exploration of atomic potential landscapes. The zero-point motion was shown to be important, but does not by itself reproduce the correct thermal expansion coefficients (28, 29). Temperature-dependent phonon shifts from pure phonon anharmonicity with zero-point energy could give a nuclear quantum effect that alters thermophysical properties. A more detailed study of the temperature dependence of phonons in silicon is therefore appropriate because very few modes were previously assessed (3, 23, 24, 30).

We report inelastic neutron scattering measurements of phonon dispersions of silicon above 300 K along with fully anharmonic ab initio calculations using the stochastically initialized temperature-dependent effective potential method (s-TDEP). This stochastic method samples and fits the phonon potential landscape the same way a Born–Oppenheimer molecular dynamics potential energy surface is fitted to a model Hamiltonian (17). This method can accurately describe highly anharmonic systems and includes higher-order contributions of the lattice dynamic Hamiltonian, which intrinsically includes the phonon–phonon interactions as well as the nuclear quantum effects (17, 31⇓⇓–34). These measurements are in conflict with the quasiharmonic theory, which predicts the wrong sign for phonon shifts with temperature. We show that the crystal structure, quasiharmonicity, pure anharmonicity, and nuclear quantum effects all play important roles in the thermal expansion of silicon. Methods for both the measurements and the calculations are described in *Materials and Methods* and *Supporting Information*.

Fig. 1 shows phonon dispersions as bright intensities. The dispersions at low temperatures are in excellent agreement with previous work that used triple-axis spectrometers (3, 30). With increasing temperature, the majority of phonon modes, including the low-energy TA modes, soften in proportion to their energy. This self-similar behavior of phonon softening was reported previously (25).

Results from calculations by the s-TDEP method (with anharmonicity and thermal expansion) and conventional quasiharmonic ab initio calculations (with no anharmonicity) are shown in Fig. 2. There are large discrepancies in the signs and magnitudes of phonon energy shifts between the two models. Most interestingly, Fig. 2 *B* and *C* shows that the s-TDEP calculations predict a reduction in phonon energy, a thermal “softening,” in the transverse modes (roughly <35 meV), whereas the quasiharmonic calculations predict an increase in phonon energy, “stiffening,” at 1,500 K [with negative Grüneisen parameters as reported previously (7⇓–9)].

We calculated the fractional shifts of energies, *q* points. Fig. 2*D* compares the density of fractional phonon shifts from quasiharmonic and anharmonic (s-TDEP) calculations. The density of fractional shifts, *E* from the s-TDEP method at 700 K. Compared with the quasiharmonic predictions for the TA modes (shown at the top of Fig. 2*D*), the anharmonic shifts are an order-of-magnitude larger, have opposite signs, and follow opposite thermal trends. Such large discrepancies allow for definitive experimental tests.

Individual phonon energies were obtained from constant *q* fits to the measured *S*(*q*,ε), as shown in *Supporting Information*. Fig. 2 *F*–*I* shows that the trends from the anharmonic s-TDEP calculations are in far better agreement with experiment than are the quasiharmonic trends. Thermal trends for individual phonons at the L, X, and K points (Fig. 2 *F*–*H*) are presented for their importance in the interpretation of quasiharmonic results (7). Another example for a phonon mode located away from a high-symmetry line is shown in Fig. 2*I*.

Additional s-TDEP calculations of densities of thermal shifts suggest why the quasiharmonic theory has been so apparently successful. Calculations were performed for volumes that were 1% larger and 1% smaller than the 0 K harmonic volume calculated for Fig. 2*A*, and the results are shown on the left and right sides of Fig. 3 for the TA modes (Fig. 3 *A*–*C*) and all phonon modes (Fig. 3 *D*–*F*). For all three volumes, at low temperatures, there is a wide spread in the thermal phonon shifts, both stiffening and softening. At low temperatures, the average thermal shift from anharmonicity at a fixed volume is, surprisingly, nearly zero. At fixed volume, the shifts of all quasiharmonic phonons are zero, of course, so the two methods agree on the average. This is seen in Fig. 3 *A*–*C* for the TA modes and in Fig. 3 *D*–*F* for all modes. Nevertheless, the average phonon energies of TA modes from the s-TDEP method show an ordinary softening with increased volume and temperature, inconsistent with the negative Grüneisen parameters from quasiharmonic calculations. At high temperatures, Fig. 3 *D*–*F* shows that all of the modes tend to soften at similar rates. Differences in vibrational entropies from the s-TDEP and quasiharmonic methods were calculated using equations in *Supporting Information*. The difference in entropies *G*–*I*. For all volumes, the differences are negligible up to 125 K but increase at higher temperatures (Fig. 3).

A quasiharmonic model with negative Grüneisen parameters gives a physically incorrect explanation of thermal expansion, although some of its predictions of average properties are preserved by gross cancellations of errors. As described in *Supporting Information*, zero-point energy (**S6** proves essential for an anharmonic model to predict the negative thermal expansion of silicon (Fig. 4). Nuclear quantum effects give nonzero anharmonic couplings between all phonons, even modes of higher energy that are not excited thermally at low temperature. These anharmonic couplings alter the self-energies of the lower-energy phonons that are excited at low temperatures, altering the volume dependence of the free energy. Calculated coefficients of linear thermal expansion are in excellent agreement with experiments (Fig. 5). Not only are quantum effects essential at lower temperatures, but differences persist up to melting temperatures. Varying the zero-point motion from changes in nuclear mass allows for an interesting engineering opportunity, too (29, 35⇓–37).

Measurements of the phonon dispersions of single-crystal silicon from 100 K to 1,500 K showed thermal shifts that contradict the trends predicted by the widely accepted QH, even at low temperatures. Pure phonon anharmonicity, i.e., phonon–phonon interactions, dominate the phonons in silicon from low to high temperatures, altering the effective interatomic potential and causing both positive and negative shifts of phonon energies. At low temperatures, the zero-point quantum occupancies of high-energy vibrational modes alter the energies of low-energy modes through anharmonic coupling. This nuclear quantum effect with anharmonicity (and quasiharmonicty) is the essential cause of the negative thermal expansion of silicon. The crystal structure, anharmonicity, and nuclear quantum effects of silicon all play important roles in the thermal expansion of silicon, and could be essential in other technologically important materials.

## Materials and Methods

For a detailed description of methods, see *Supporting Information*.

The experiments used a high-purity single crystal of silicon (mass ≈ 28.5 g) with ⟨110⟩ orientation, machined into a tube for optimal neutron scattering properties. The sample was rotated in a furnace on a direct geometry time-of-flight inelastic neutron scattering spectrometer called wide-angular range chopper spectrometer (ARCS) (38) at the Spallation Neutron Source at Oak Ridge National Laboratory. For each temperature, the 4D *S*(*q*,ε) data were reduced and multiphonon scattering was subtracted to give all phonon dispersions in the irreducible wedge of the first Brillouin zone. The multiphonon scattering produces a relatively smooth background between the phonon dispersions and was determined to produce the majority of the background intensity (*Supporting Information*) (39). Our “folding” technique of summing all of the *S*(*q*,ε) data (from >100 Brillouin zones) into an irreducible wedge increases the signal strength, suppresses polarization effects that alter intensities in some Brillouin zones (39), and averages out any possible effects of “anharmonic interference” (40).

All ab initio calculations were performed with the Vienna Ab Initio Simulation Package (VASP) (41⇓⇓⇓⇓⇓–47). An s-TDEP (17, 31, 33, 48) was implemented to obtain phonon shifts with temperature, including intrinsic phonon anharmonicities and nuclear quantum effects. Quasiharmonic calculations were also conducted as described previously (25).

## Acknowledgments

The authors thank F. H. Saadi, A. Swaminathan, I. Papusha, and Y. Ding for assisting in sample preparation and discussions. Research at Oak Ridge National Laboratory’s Spallation Neutron Source (SNS) was sponsored by the Scientific User Facilities Division, Basic Energy Sciences (BES), Department of Energy (DOE). This work used resources from National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the US Department of Energy under Contract DE-AC02-05CH11231. Support from the Swedish Research Council Program 637-2013-7296 is also gratefully acknowledged. Supercomputer resources were provided by the Swedish National Infrastructure for Computing. This work was supported by the DOE Office of Science, BES, under Contract DE-FG02-03ER46055.

## Footnotes

↵

^{1}D.S.K. and O.H. contributed equally to this work.- ↵
^{2}To whom correspondence may be addressed. Email: dennis.s.kim{at}icloud.com or btf{at}caltech.edu.

Author contributions: D.S.K., H.L.S., and B.F. designed research; D.S.K., O.H., H.L.S., J.L.N., C.W.L., D.L.A., and B.F. performed research; D.S.K., O.H., J.H., N.S., J.L.N., C.W.L., D.L.A., and B.F. contributed new reagents/analytic tools; D.S.K., O.H., J.H., H.L.S., J.Y.Y.L., N.S., J.L.N., C.W.L., D.L.A., and B.F. analyzed data; and D.S.K., O.H., J.H., H.L.S., J.Y.Y.L., N.S., J.L.N., C.W.L., D.L.A., and B.F. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission. A.M. is a guest editor invited by the Editorial Board.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1707745115/-/DCSupplemental.

Published under the PNAS license.

## References

- ↵
- Fultz B

- ↵
- Fultz B

- ↵
- Nilsson G,
- Nelin G

- ↵
- ↵
- ↵
- ↵
- Wei S,
- Li C,
- Chou MY

- ↵
- Liu ZK,
- Wang Y,
- Shang S

- ↵
- Xu CH,
- Wang CZ,
- Chan CT,
- Ho KM

- ↵
- ↵
- Rignanese GM,
- Michenaud JP,
- Gonze X

- ↵
- Middelmann T,
- Walkov A,
- Bartl G,
- Schödel R

- ↵
- ↵
- ↵
- ↵
- Shah JS,
- Strauman ME

- ↵
- Hellman O,
- Steneteg P,
- Abrikosov IA,
- Simak SI

- ↵
- Weinstein BA,
- Piermarini GJ

- ↵
- ↵
- Narasimhan S,
- Vanderbilt D

- ↵
- Balkanski M,
- Wallis RF,
- Haro E

- ↵
- Wang CZ,
- Chan CT,
- Ho KM

- ↵
- Menéndez J,
- Cardona M

- ↵
- Tsu R,
- Hernandez JG

- ↵
- Kim DS, et al.

- ↵
- Debernardi A

- ↵
- Lang G, et al.

- ↵
- Allen PB

- ↵
- Herrero CP,
- Ramírez R

- ↵
- Brockhouse BN

- ↵
- Klein ML,
- Horton GK

- ↵
- Shulumba N,
- Hellman O,
- Minnich AJ

- ↵
- Errea I,
- Calandra M,
- Mauri F

- ↵
- Wallace DC

- ↵
- Herrero CP,
- Ramírez R,
- Cardona M

- ↵
- Herrero CP

- ↵
- Noya JC,
- Herrero CP,
- Ramírez R

- ↵
- ↵
- Squires GL

- ↵
- Glyde HR

- ↵
- Kresse G,
- Hafner J

- ↵
- Kresse G,
- Hafner J

- ↵
- ↵
- Kresse G,
- Furthmüller J

- ↵
- ↵
- Mattsson AE,
- Armiento R

- ↵
- Armiento R,
- Mattsson AE

- ↵
- Hellman O,
- Abrikosov IA

- Stone MB,
- Niedziela JL,
- Loguillo MJ,
- Overbay MA,
- Abernathy DL

- Mantid (2013) Manipulation and analysis toolkit for instrument data. Mantid Project. Available at https://dx.doi.org/10.5286/SOFTWARE/MANTID. Accessed March 2, 2014.
- Blochl PE

- Gonze X,
- Lee C

- Yao Y,
- Kanai Y

- Dove MT

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Applied Physical Sciences