# Phonon anharmonicity and negative thermal expansion in SnSe

###### Abstract

The anharmonic phonon properties of SnSe in the Pnma phase were investigated with a combination of experiments and first-principles simulations. Using inelastic neutron scattering (INS) and nuclear resonant inelastic X-ray scattering (NRIXS), we have measured the phonon dispersions and density of states (DOS) and their temperature dependence, which revealed a strong, inhomogeneous shift and broadening of the spectrum on warming. First-principles simulations were performed to rationalize these measurements, and to explain the previously reported anisotropic thermal expansion, in particular the negative thermal expansion within the Sn-Se bilayers. Including the anisotropic strain dependence of the phonon free energy, in addition to the electronic ground state energy, is essential to reproduce the negative thermal expansion. From the phonon DOS obtained with INS and additional calorimetry measurements, we quantify the harmonic, dilational, and anharmonic components of the phonon entropy, heat capacity, and free energy. The origin of the anharmonic phonon thermodynamics is linked to the electronic structure.

###### pacs:

63.20.kg, 63.20.Ry, 63.20.dk, 65.40.De^{†}

^{†}preprint:

Notice: This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05- 00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

## I Introduction

Thermoelectric materials are of current interest for cost-effective, reliable power generation applications, either by utilizing natural sources or recovering man-made waste heat Snyder and Toberer (2008); Zebarjadi et al. (2012). In order to increase their heat-to-electricity conversion efficiency, thermoelectric materials need to have a low thermal conductivity, while maintaining a high electrical conductivity. A detailed understanding of phonons is necessary to rationalize both the thermodynamics and the thermal transport properties in thermoelectrics. Phonons are the main contributors to the entropy and heat capacity Wallace (1972, 2003); Grimvall (1999); Fultz (2010). In addition, phonons are the dominant heat carriers in semiconductors, and understanding deviations from harmonic lattice dynamics is necessary to account for the bulk lattice thermal conductivity, . From a fundamental standpoint, phonons are sensitive to the nature of chemical bonding Rhyee et al. (2009); Nielsen et al. (2013); Lee et al. (2014); Li et al. (2015a) and the temperature dependence of phonon spectral functions often provides important insights into the anharmonicity of the interatomic potential, directly affecting thermal resistivity through phonon-phonon scattering Delaire et al. (2011a); Shiga et al. (2012); Li et al. (2015a). The coupling between phonons and electronic structure also has an important effect on phonon frequencies (and linewidths) Grimvall (1999); Delaire et al. (2011b, 2008a, 2008b); Bansal et al. (2015), and electronic instabilities can produce large anharmonicity, helping to suppress and improve thermoelectric efficiency Li et al. (2015a); Lee et al. (2014); Hong and Delaire (2016).

Tin-selenide was recently reported to achieve an outstanding thermoelectric figure-of-merit, with Na-doped samples maintaining a high power factor over a broad range of temperatures Zhao et al. (2014); Sassi et al. (2014); Chen et al. (2014); Zhao et al. (2016). At high temperature, SnSe undergoes a continuous structural phase transition between the low-symmetry Pnma (below K) and the higher-symmery Cmcm phases (above ). This phase transition has been the subject of a number of early crystallographic and thermomechanical experimental investigations Chattopadhyay et al. (1986); Adouby et al. (1998); Wiedemeier and Csillag (1979); Chattopadhyay et al. (1984); von Schnering and Wiedemeier (1981); Feutelais et al. (1996); Sharma and Chang (1986); Balde et al. (1981); Dembovskii et al. (1963); Zhdanova (1961). More recently, the ultra-low lattice thermal conductivity of SnSe has attracted great interest for thermoelectric applications, and has been studied both experimentally and through ab-initio simulations Zhao et al. (2014); Li et al. (2015a); Carrete et al. (2014); Sassi et al. (2014); Ding et al. (2015); Loa et al. (2015, 2016); Chen et al. (2014); Guo et al. (2015); Kutorasinski et al. (2015); Wei et al. (2015); Sanchez et al. (2015); Li et al. (2015b); Han et al. (2015); Ge et al. (2015). The coupling between lattice distortion and electronic structure was recently shown to be responsible for the large anharmonicity, which persists in a broad range of temperatures, and leads to the low thermal conductivity Li et al. (2015a); Hong and Delaire (2016).

The crystal structure of SnSe is illustrated in Fig. 1. Importantly, x-ray and neutron diffraction studies showed that, in the Pnma phase, the axis, parallel to the direction of corrugation of Sn-Se bilayers below , exhibits a negative thermal expansion (NTE) coefficient Adouby et al. (1998); Wiedemeier and Csillag (1979); Chattopadhyay et al. (1986). The thermodynamic and atomistic origins of this effect in SnSe remain to be established. The NTE implies that the thermodynamic generalized Grüneisen tensor must have negative values for the diagonal element corresponding to the in-plane layer Grimvall (1999). Yet, this is in contrast to positive values predicted from quasiharmonic ab-initio simulations Zhao et al. (2014); Li et al. (2015a); Guo et al. (2015). Further, our inelastic neutron scattering measurements showed pronounced differences in the behaviors of phonons propagating along in- and out-of-plane crystallographic axes Li et al. (2015a). Since the coupling between the lattice distortion and the electronic structure is a strong source of anharmonicity Hong and Delaire (2016), it is important to understand its role in the unusual thermal expansion of SnSe.

Here, we investigate the origin of the NTE and anharmonic phonon thermodynamics in SnSe, focusing on the Pnma phase, using inelastic neutron scattering (INS) and nuclear resonant inelastic X-ray scattering (NRIXS) measurements, and first-principles simulations. Our measurements show a strong temperature dependence of the phonon frequencies, revealing a high degree of anharmonicity in the interatomic potential. We show how this anharmonicity is related to the structural phase transition. We model the thermal expansion and its anisotropic behavior, from first-principles, especially the lattice contraction parallel to the Sn-Se layers, and find good agreement with reported lattice parameter measurements. In addition, we quantify the harmonic, dilational, and anharmonic contributions of phonons to the entropy, internal energy and free energy.

## Ii Sample Preparation

Single crystals of SnSe were synthesized from high purity Sn and Se (Alfa Aesar, 99.999%) in fused silica ampoules. After an initial reaction, a polycrystalline precursor was melted and subsequently crystallized while cooling at 0.4-0.5C/h; a 24 hr annealing at 830C occurred before cooling to room temperature. One of the single crystals was ground into a fine powder for neutron scattering measurements of the phonon density of states. Some properties of these samples were previously reported, and further details can be found in Ref. 10

## Iii Spectroscopy

### iii.1 Inelastic Neutron Scattering

We measured the phonon DOS of polycrystalline (mass g) using the time-of-flight wide angular-range chopper spectrometer (ARCS) at the Spallation Neutron Source (SNS) at Oak Ridge National Laboratory Stone et al. (2014). The powder sample was contained in a standard thin-walled aluminum can. Measurements at low ( K) and high temperature ( K) were performed using aluminum can inside a closed-cycle helium refrigerator and low-background resistive furnace. At low temperatures, we filled the sample chamber with low pressure of helium to facilitate cooling. We used two incident energies, meV and 55 meV, when combined it provided an entire phonon spectrum and high resolution datasets to distinguish between closely spaced phonon peaks. An oscillating radial collimator was employed to minimize scattering from the sample environment. We followed the similar steps in data normalization, reduction to - (momentum and energy transfer) space, background subtraction, multiphonon scattering, and removal of elastic peak as described in our previous work Bansal et al. (2015).

The neutron scattering cross-sections are different for different elements. The phonon DOS derived from INS measurements is weighted by the respective scattering cross-sections () and mass (m). The neutron-weighting factors for Sn, and Se are , and , respectively (in units of barns/amu). Due to the larger neutron-weighting factor for Se in comparison to Sn, phonon DOS is over-weighted by Se vibrational contributions. Accordingly, the neutron-weighted phonon DOS can be expressed as following:

(1) |

where , and are the partial densities of states of Sn, and Se. For comparison of the experimental phonon DOS with the simulations, the neutron-weighting factor can be applied to the partial phonon DOS.

### iii.2 Nuclear Resonant Inelastic X-ray Scattering

NRIXS experiments on single-crystal SnSe platelets (natural Sn) were performed at Sector 30 beamline at the Advanced Photon Source, Argonne National Laboratory. Since SnSe exhibit anisotropic behavior, we could extract the direction-projected phonon DOS by changing the orientation of the sample with respect to the incident x-ray beam. For the -axis projected phonon DOS, the sample was mounted such that the incident x-ray beam was parallel to -axis, while for the projected phonon DOS perpendicular to -axis, the incident beam was parallel to the b-c plane. The incident x-ray energy was 23.88 keV, the nuclear resonance energy of Sn. The energy bandpass of incident beam was reduced to 1.2 meV using a high resolution crystal monochromator Toellner et al. (2011). The inelastic signal collected is comprised mainly of 23.88 keV nuclear fluorescence signal, delayed in time relative to the prompt pulse by 20-130 nsec. The inelastic spectra were collected over an energy range of 160 meV with 0.25 meV steps. Repeat scans were taken and data were subsequently added. The Sn specific phonon DOS was extracted from the NRIXS data using the PHOENIX package Sturhahn (2000). The samples were single-crystal SnSe platelets (no isotopic enrichment), cut from INS crystals. The platelet samples were mounted with either transmission geometry (beam along –axis) or grazing geometry (beam nearly parallel to plane). Because of the crystalline anisotropy, the direction-projected partial DOS for each configuration provides information about different Sn motions ( or , respectively).

## Iv Modeling

### iv.1 Structural Relaxation and Phonon Dispersions

First-principles simulations were performed in the framework of density functional theory (DFT) as implemented in the Vienna Ab initio Simulation Package (VASP 5.3) Kresse and Hafner (1993); Kresse and Furthmller (1996a, b). A Monkhorst-Pack electronic k-point mesh was used for all of our simulations, with a plane-wave cut-off energy of 500 eV. The projector-augmented-wave potentials explicitly included four valence electrons for Sn (), and 6 for Se (). The lattice parameters and atomic positions were optimized until forces on all atoms were smaller than 1 meV Å. We carefully evaluated the accuracy of our phonon calculations, comparing the local-density approximation (LDA) exchange-correlation (XC) functional Perdew and Zunger (1981) and the generalized gradient approximation (GGA) with XC functional given by the Perdew-Burke-Ernzerhof (PBE) parametrization Perdew et al. (1996). The relaxed unit cell parameters are compared with the experimentally reported structure (Ref. 23) at 300 K in Table 1. The agreement between the LDA and PBE structures and the experimental structure is good, with, as expected, LDA underestimating the lattice constants and PBE overestimating them, by opposite amounts.

The phonon dispersions were calculated in the harmonic approximation, using the finite displacement approach as implemented in Phonopy A. Togo et al. (2008), with the atomic forces in the distorted supercells obtained with VASP. To construct the force-constant tensor using the finite displacement approach in Phonopy A. Togo et al. (2008), we computed 8 independent atomic displacements, each with atomic displacement amplitudes of one hundredth of lattice constants. The LDA phonon calculations used a supercell of the primitive cell, containing 256 atoms, while GGA calculations used up to supercells containing 600 atoms. We compare the phonon dispersion and phonon DOS of SnSe at 300 K with experimental phonon dispersions measured with INS in Fig. 2, and with the NRIXS and INS phonon DOS in Fig. 11 and 6, respectively. The experimental phonon dispersions in Fig. 2 were extracted from INS datasets originally reported by our group in Ref. 10. From these comparisons with INS measurements, one can see how the PBE simulations systematically underestimate the phonon frequencies in this material, while LDA simulations are much more accurate. Even with the larger supercell, the GGA phonon frequencies show a worse agreement with INS data than the LDA calculations. The mean phonon energy calculated from experimental phonon DOS at 300 K (corrected for neutron weighting), GGA simulations, and LDA simulations are 12.92, 12.15, and 12.90 meV, respectively. Additionally, the GGA phonon group velocities are on average 10% lower than the LDA values near the point. Thus, first-principles calculations of the lattice thermal conductivity () using GGA (Refs. 32; 33) are expected to carry a sizeable systematic error, since is proportional to , with a further dependence of phonon scattering rates on the dispersions, while LDA calculations are expected be more accurate Guo et al. (2015).

LDA | GGA | Exp. Adouby et al. (1998) | |
---|---|---|---|

a (Å) | 11.309 (1.68) | 11.756 (-2.20) | 11.502 |

b (Å) | 4.119 (0.82) | 4.205 (-1.25) | 4.153 |

c (Å) | 4.300 (3.37) | 4.547 (-2.18) | 4.450 |

V (Å) | 200.302 (5.77) | 224.776 (-5.74) | 212.567 |

LDA | GGA | Exp. Adouby et al. (1998) | ||||
---|---|---|---|---|---|---|

x | z | x | z | x | z | |

Sn | 0.1169 | 0.0908 | 0.1206 | 0.1144 | 0.1208 | 0.1060 |

Se | 0.8580 | 0.4753 | 0.8551 | 0.4754 | 0.8551 | 0.4760 |

### iv.2 Negative Thermal Expansion

Non-cubic materials often exhibit anisotropic thermal expansion coefficients, reflecting the anisotropy in crystal structure and bonding. In the case of orthorombic SnSe, the anisotropy is particularly striking, with negative thermal expansion (NTE) of the axis, parallel to the direction of corrugation of Sn-Se layers (also corresponding to the direction of atomic motions across the phase transition), while the other in-plane axis () and the out-of-plane –axis show normal, positive expansion Wiedemeier and Csillag (1979); Chattopadhyay et al. (1986); Adouby et al. (1998). The relative change of lattice parameters with temperature reported by Adouby et al. in Ref. 23 are shown in Fig. 3. This distortion can also be analyzed as a rotation of SnSe tetrahedra, as illustrated in Fig. 1 and further discussed below. We proceed to rationalize this behavior based on first-principles thermodynamics and a microscopic analysis of bonding.

We start by noting that simple Grüneisen parameter calculations are insufficient to explain the anisotropy and NTE along one or more crystallographic axes. Instead, one needs to compute the 3-D total free energy surface, whose temperature dependence is dominated by the phonon entropy. In order to clearly make this point, we start by performing the computation of Grüneisen parameters. A generalized Grüneisen tensor for a given polarization index and wave vector is given by: Brugger (1965); Cantrell (1980)

(2) |

where, is the component of applied isothermal strain, and is the atomic vibration frequency. Summation over polarization index and wave vector, weighted by heat capacity at constant configuration , yields an average value of a generalized Grüneisen tensor. We have calculated the generalized Grüneisen tensor from our DFT simulations of phonons . In the case of , considering the symmetry of the second order tensor for the orthorhombic unit cell, non-diagonal components are exactly zero. Thus, the only isothermal strains required are , , and , corresponding to stretching or compressing the lattice along crystallographic axes , , and , respectively. Subsequently, the generalized Grüneisen tensor, is obtained by averaging over all polarizations and wave vectors, with each mode weighted by its heat capacity at constant configuration as follows:

(3) |

We have also calculated the Grüneisen tensor along high symmetry directions by averaging over polarizations only. Results are listed in Table 2. As can be seen from table 2, isothermal strains along , , and axes lead to positive Grüneisen parameters i.e., softening of phonon branches along all directions, except for isothermal strain along the axis, which produces a stiffening of phonon modes along and . However, summation over the entire Brillouin zone at constant configuration yields positive effective Grüneisen parameters for all cases of isothermal strain as noted in Eq. (3) and Table 2, which is at odds with the experimentally observed negative thermal expansion along the axis.

Isothermal strain | |||
---|---|---|---|

(–axis) | (–axis) | (–axis) | |

BZ integration | 1.44 | 1.35 | 0.64 |

3.26 | 1.34 | -1.01 | |

1.58 | 1.48 | 0.27 | |

1.27 | 1.35 | 0.69 | |

1.69 | 1.38 | 0.61 | |

1.66 | 1.43 | 0.88 | |

1.05 | 1.59 | 1.97 | |

1.66 | 1.40 | 0.85 | |

1.71 | 1.13 | -0.06 |

Alternatively, the relationship between the thermodynamic generalized Grüneisen tensor and the thermal expansion tensor, , can be expressed in terms of the isothermal bulk modulus , volume , and heat capacity at constant configuration , as following: Wallace (1972)

(4) |

The isothermal bulk modulus tensor, , can be determined either from the phonon group velocities along different crystallographic directions, or from DFT simulations by applying isothermal longitudinal and transverse strains. We have calculated with both approaches. The independent components of , constrained by the orthorhombic symmetry, are shown in Table 3. The phonon group velocities at 300 K were obtained from INS on single crystals, previously reported by our group in Ref. Li et al. (2015a). The tensor components not listed in table 3 are zero. As can be seen in table 3, the values for obtained from the two approaches are in fair agreement.

DFT (GPa) | Exp. (GPa) | |
---|---|---|

85.5 | 82.5 | |

89.1 | 97.6 | |

52.0 | ||

41.2 | ||

23.5 | ||

19.7 | ||

49.8 | 31.4 | |

22.9 | 17.9 | |

19.1 | 23.7 |

From the experimentally measured thermal expansion coefficient Wiedemeier and Csillag (1979), heat capacity (calculated from experimentally measured phonon DOS, see Fig. 6), volume Adouby et al. (1998), and bulk modulus (DFT values were used when experimental values were unavailable) at T = 300 K, the thermodynamic generalized Grüneisen tensor computed from Eq. (4) is:

(5) |

Thus, (Eq. (5)), and (Eq. (3)) are quite different. The Grüneisen parameter estimated from experimental data is negative along the axis, but the DFT prediction is positive. Importantly, while we do observe a substantial negative along by applying the strain along (cf Table 2), this alone is insufficient to produce negative along as a whole, once the integration is performed over the entire Brillouin zone. We note that this is not a failure of the DFT calculation per se, but rather arises because of the omission of couplings between different crystallographic axes.

The difficulty to capture the NTE along from the Grüneisen parameter prompted us to instead compute the three-dimensional free energy surface, on a 3D grid of isothermal strains. To construct the total 3-D free energy surface, we created a grid of lattice parameters with 1% expansion and compression of lattice parameter along each crystallographic axes with total of 27 grid points. The inherent advantage of this approach is that it captures the evolution of both electronic and vibrational free energy by varying lattice parameters individually (i.e., {}, {}, and so on), as well as coupled terms (i.e., {}, {}, and so on), and determine the total free energy minima. The total free energy is a sum of electronic free energy, , and vibrational free energy, , which can be expressed as following:

(6) |

and

(7) |

where electronic entropy, , and the Fermi-Dirac occupation function, , are given by:

(8) |

and

(9) |

respectively. Here, is ground state energy, is phonon energy of mode, is the energy at Fermi level, is Boltzmann’s constant, and is the number of atoms in the unit cell. We obtain the ground-state electronic energy, and the electronic density of states (eDOS) from DFT simulations for all points on the grid. The electronic entropy was determined from eq. (8), and (9). The vibrational free energy was obtained from eq. (7), using the phonon density of states (phonon DOS) calculated as described in section IV.1.

The measured thermal expansion between 300 K and 813 K along , , and is approximately 1.8%, 3.5%, and -3.1%, respectively Adouby et al. (1998). Additional points were added to extend our grid to larger expansions/compressions of 3% along , 5% along , and 5% along . However, phonons become unstable for extended grid points. Thus, we limit our procedure to the original 27-point grid, extended with the seven new points for which phonons are stable, shown in Fig. 3(b). To capture the NTE along over a wider range of temperature, we performed a quadratic extrapolation of the free energy surface (see Eq. (IV.2)). We note that the computational cost scales cubicly with the linear grid size, hence the practicality of finer grids is limited.

We determine the equilibrium lattice parameters, , , and by minimizing a second-order polynomial fit to the free-energy, , with respect to lattice parameters at ,

(10) |

The results are shown in Fig. 3(a). As can be seen on this figure, the experimental temperature dependence of lattice parameters is qualitatively captured, including the NTE along . This indicates an interplay between electronic and vibrational components of the free energy, which ultimately controls the minima of . The deviation between calculated and experimental lattice parameters at high is partially due to the fact that the theoretical lattice parameters at , Å, are smaller than the experimental values: Å, leading to shorter theoretical bond lengths, especially in the plane. In addition, strong anharmonic effects, not included in our quasiharmonic calculations of phonon spectra for strained cells, are seen to occur in the last K below in our heat capacity measurements (see below). One could also include anharmonic effects in the phonon calculations, for example with ab initio molecular dynamics or other methods incorporating anharmonic phonon renormalization occurring because of the thermal bath Li et al. (2015a); Zhao et al. (2014); Carrete et al. (2014). However, our goal is not to exactly match the measured thermal expansion coefficient, but rather to rationalize the origin of the anisotropy and NTE along , which can be explained with this simple and insightful approach.

The electronic ground state energy at K is minimum for the fully relaxed configuration (relaxed lattice parameters and atomic positions), , and increases with any strain. We note that, for the strained configurations, the minimum in the electronic ground state energy is related to the lattice parameters such that there are two possible scenarios: or . Hence, the electronic ground state energy alone can not uniquely define the equilibrium lattice structure. Thus, it is the vibrational free energy that brings the minimum of towards . The calculated phonon DOS for three different configurations is shown in Fig. 4. From this figure, we observe that leads to additional spectral weight at low energy in the phonon DOS (“softening”), while results in a stiffening. The softening in the first case increases the phonon entropy, and thus lowers . As further discussed below, the configuration leads to considerable softening of high energy phonons.

### iv.3 Structural Evolution and NTE

Cmcm | 2.75 | 0 | 3.05 | 89.85 | 4.32 | 4.31 | 4.31 | 4.31 |
---|---|---|---|---|---|---|---|---|

Pnma | 2.72 | 7.71 | 2.80 | 95.67 | 4.49 | 4.15 | 4.45 | 4.15 |

We now investigate the structural changes inside the unit cell in connection with the NTE. In particular, we show that the NTE of SnSe is closely related to the electronic instability and its coupling to the anharmonic vibrations of Sn atoms along . In a separete work, we showed the strength of bonds is considerably weakened in Pnma compared to Cmcm, as the chemical bond between Sn and Se in essentially breaks Hong and Delaire (2016). Thus, the apical bond and the two in-plane bonds can be seen to form triangular pyramids with Sn at the apex, as shown in Fig. 1. The phase transition from Cmcm to Pnma is driven by the electronic instability, which reduces the energy by lifting the degeneracy of four bonds in-plane in the Cmcm phase, favoring the off-centering of Sn along , and breaking a mirror plane symmetry. Hong and Delaire (2016) This distortion of the Cmcm phase overlaps most strongly with a strongly anharmonic soft mode at the zone boundary (optic mode at point in Cmcm), coupled to another anharmonic zone-center () optic mode, which together induce rotations of those pyramids as pseudo-rigid units Hong and Delaire (2016), and cause an increase in upon cooling.

In Fig. 1a, we introduce as the rotation angle of with respect to , and the angle between the bonds in the pyramid. For simplicity, we ignore the tilt of the plane formed by two bonds, away from plane. Within a close approximation, we have:

(11) |

(12) |

From Table 4 we can see that decreases from Cmcm to Pnma structure, and therefore, the expansion along can be seen to result from the rotation of pyramids (), which couples strongly to the mode ( almost remains constant). From the plots in Fig. 5, we can confirm that the NTE along is indeed induced by the rotation of pyramids, which is strongly coupled to the anharmonic modes (soft-modes in Cmcm phase generating zone-center soft-modes in Pnma phase) Hong and Delaire (2016). The strongly anharmonic character of these modes partly explains the shortcomings of the quasiharmonic free-energy calculations near the phase transition, seen in Fig. 3.

## V Anharmonic Phonon Thermodynamics

### v.1 Phonon DOS and Phonon Entropy

Much of the free-energy discussion above was limited to a quasiharmonic approximation. In order to estimate the anharmonic component of the vibrational entropy, we investigate the temperature dependence of the measured phonon DOS. Fig. 6a and 6b show the experimental phonon DOS measured at different temperatures for meV and 30 meV, respectively. The low and high energy phonon peaks soften by 0.3 and 1.0 meV, respectively as the sample temperature is increased from 10 to 750 K. One can observe a drastic broadening of the spectrum with increasing , in particular for optical modes above 13 meV. While the lowest energy zone-center TO mode has been shown to soften strongly with increasing in Pnma in our previous INS dispersion measurements (Ref. 10), that mode is confined to the zone-center and has a limited spectral weight in the DOS. We note that the neutron-weighted phonon DOS from harmonic DFT simulations agrees well with the total and direction-projected phonon DOS measured at 300 K (Fig. 11 and 6). The direction-projected experimental phonon DOS measurements (Fig. 11) provide substantial information in separating the two low energy Sn dominated phonon peaks at 6 and 11 meV. The low energy (6 meV) phonon peak include both out-of-plane (along -axis) and in-plane (b-c plane) Sn vibrations, while 11 meV phonon peak primarily originates from in-plane (b-c plane) Sn vibrations. Both Sn dominated phonon peaks show substantial broadening in the phonon spectrum with increasing temperature (Fig. 6).

With increasing temperature, it is the coupling between lattice distortion (especially the large Sn in-plane displacements) and the electronic structure that generates strong anharmonicity (leading to phonon-phonon interactions) Li et al. (2015a), primarily responsible for a significant phonon softening of high energy phonons and phonon broadening of entire phonon spectra seen in Fig. 6. This conclusion is further supported by our DFT simulations of NTE modeling where we identify that the minimization of electronic free energy alone is not sufficient to obtain the experimental lattice parameters with temperature, and including the change in vibrational free energy due to lattice distortions is necessary to reach equilibrium. As we have detailed in Ref. 17, the origin of the anharmonicity for Sn in-plane motions is the Jahn-Teller instability involving the Se p-states.

To quantify the anharmonic contribution to , we estimated a phonon broadening parameter, , by convolving the generalized (neutron-weighted) phonon DOS from DFT, , with a damped harmonic oscillator function, , Lovesey (1984)

(13) |

and matching the result with the measured phonon DOS. Here,

(14) |

and is a shift in phonon energy calculated by fitting the Lorentzian to the peaks in the measured phonon DOS. The resulting phonon DOS, , with temperature-dependent phonon broadening parameter, , is shown in Fig. 7. This model reproduces the temperature dependence of the experimental phonon DOS well, and enables us to assess the anharmonic . For this purpose, we use the non-neutron-weigthed version of .

The vibrational entropy can be expressed as Fultz (2010); Wallace (1972, 2003),

(15) |

where is the Bose-Einstein occupation factor for phonons, and is the number of atoms in the crystal.

The harmonic phonon entropy, , was calculated by convolving the harmonic phonon DOS from DFT with a damped anharmonic oscillator function at K. The resulting phonon DOS was substituted in Eq. (V.1) to evaluate at different temperatures. To evaluate the total entropy, , at different temperatures, the –dependent damped anharmonic oscillator function was convolved with the DFT phonon DOS, and substituted in Eq. (V.1). The vibrational entropy due to thermal expansion of the lattice is given by

(16) |

where is the volumetric thermal expansion coefficient, is the compressibility (inverse of bulk modulus), and and are specific heat capacity of material at constant pressure and temperature.

Using the total , harmonic , and dilational components of the vibrational entropy, the nonharmonic and anharmonic entropy were calculated as: , and , respectively. For the calculation of dilational entropy, we used the thermal expansion coefficient from Wiedemeier Wiedemeier and Csillag (1979), and temperature dependent compressibility from He et al. He et al. (2013). In addition, we have also estimated the total entropy from our experimental heat capacity measurements () using the following expression,

(17) |

The calculated harmonic, dilational, non-harmonic, anharmonic, and total vibrational entropy values are shown in Fig. 8.

As we can observe from Fig. 8, the main entropic contribution is the harmonic term. However, as temperature increases, the nonharmonic contribution becomes significant and cannot be ignored. The nonharmonic entropy is a summation of dilational and anharmonic entropy contributions. We have also calculated nonharmonic contributions to the vibrational internal energy and vibrational free energy following methods recently reported in Ref. 57. We should note that since tin-selenide is a indirect band gap semiconductor (band gap of 0.86 eV in Pnma phase at room temperature)Zhao et al. (2014), we could expect the contribution of electronic excitations to specific heat capacity/thermal expansion coefficient to be negligible at room temperature. The band gap in SnSe is temperature dependent, and reduces to eV in Cmcm phase Zhao et al. (2014). While the temperature dependence of band gap may affect our results slightly, we do not expect this effect to be significant. From these calculations, a similar trend is observed in vibrational internal energy (Fig. 9) and vibrational free energy (Fig. 10), with the anharmonic contribution rising as is approached. As we can observe, at high , thermal expansion (dilation) alone cannot account for nonharmonic contributions, and anharmonicity plays a large role near the instability.

It is interesting to note that in the experimental phonon DOS (Fig. 6), there is considerably larger softening, %, of high energy phonons ( meV) compared to the softening of low energy phonons, %, between and 750 K. As evident from Fig. 11a, low energy and high energy modes are dominated by and , respectively. We selectively investigate the change in low and high energy phonon energies by applying 1% isothermal strain along (3-D free energy surface minima direction as described earlier), and find that high energy phonons soften by % in comparison to % for low energy phonons.

The values of Grüneisen parameter along high-symmetry directions are reported in Table 5. The large values of Grüneisen parameter across the various high symmetry directions in reflect the pronounced anharmonicity, which is also shown quantitatively in the anharmonic vibrational entropy, internal energy and free energy calculations, while the significant variation in values across different directions is attributed to the structural anisotropy in this material.

Isothermal strain along: | |||

All modes | modes 13 meV | modes 13 meV | |

5.51 | 6.49 | 4.49 | |

2.79 | 2.92 | 2.66 | |

1.91 | 1.39 | 2.46 | |

2.46 | 2.08 | 2.86 | |

2.20 | 1.79 | 2.63 | |

0.63 | -0.39 | 1.68 | |

2.21 | 1.83 | 2.61 | |

2.90 | 2.59 | 3.23 |

### v.2 Partial Phonon DOS and Thermal Displacement Parameters

Fig 11-a shows the partial phonon DOS of tin and selenium calculated from DFT simulations for , the ground state equilibrium lattice parameters for fully relaxed structure with Sn and Se atoms predominantly occupying the low and high energy phonon spectrum, respectively. To enable the direct comparison with experimentally measured NRIXS Sn partial phonon DOS parallel and perpendicular to -axis, we have also calculated the projected phonon DOS of Sn as shown in Fig 11-b, and 11-c. The good agreement between experiment and simulations further validates the accuracy of DFT-LDA. Additionally, from the partial phonon DOS and f-factor (Lamb-Mössbauer factor) measured from NRIXS Sturhahn and A.Chumakov (1999); Sturhahn (2000), the mean square thermal displacement of Sn can be calculated. The f-factors corresponding to partial phonon DOS parallel and perpendicular to -axis are 0.13700.0098 and 0.14870.0095, respectively. We should note that, the exact sample orientation for in-plane ( a) was not known; however, the in-plane NRIXS measurement agrees well with the simulation for the -axis and we believe this was the likely orientation.

Furthermore, we have calculated the mean square thermal displacement parameter, , from DFT phonon spectra, as:

(18) |

where, denotes thermal displacement direction, is mass of atom at location , is phonon frequency of phonon branch at wave vector , is the mean Bose-Einstein occupation factor, and is phonon wave vector. The results for K, listed in Table 6, are consistent with values measured by Chattopadhyay et al. at room temperature with neutron diffraction Chattopadhyay et al. (1986). The thermal displacement parameters for and have similar magnitudes at 300 K, thus indicating that the average bond stiffness at each site is comparable. Indeed, our DFT simulations give a small increase ( 12%) in the trace of Se on-site force-constant matrix in comparison to Sn. The ratio of average frequency for Se and Sn vibrations calculated from DFT is 1.38, close to the effect of mass ratio, . It is worth emphasizing, that while the magnitude of thermal displacements of Se atoms reported by Chattopadhyay et al. increases linearly with temperature, off-centered Sn atoms show a superlinear increase for in-plane thermal displacements near the phase transition Chattopadhyay et al. (1986); Adouby et al. (1998), which reflects the strong anharmonicity in this regime.

(in Å) | This work (DFT) | This work (NRIXS) | Ref. Chattopadhyay et al. (1986) | This work (DFT) | Ref. Chattopadhyay et al. (1986) |
---|---|---|---|---|---|

0.0165 | 0.0136 | 0.0151 | 0.0143 | 0.0130 | |

0.0148 | 0.0130 | 0.0130 | 0.0118 | 0.0107 | |

0.0170 | – | 0.0177 | 0.0126 | 0.0143 |

### v.3 Heat capacity

To further investigate the nature of the phase transition and anharmonic effects in SnSe, we measured the heat capacity. The measurements were performed with a Netzsch DSC 404C differential scanning calorimeter, with the sample loaded inside a Pt crucible, under an ultra-pure Ar purge gas cycled through a Ti gettering furnace. The scans were performed after careful evacuation and purging of the sample chamber. The heating and cooling rates were 20 K/min and 20 K/min, respectively. A sapphire standard and empty-crucible baseline measurements were performed in identical conditions. The heat capacity curves measured during heating and cooling are shown in Fig. 12. The heat capacity exhibits a lambda shape in the vicinity of the phase transition, akin to the classic case of liquid helium Buckingham and Fairbank (1961). This behavior is generally indicative of a second-order phase transition, in agreement with the nearly continuous evolution of structural parameters observed with diffraction Adouby et al. (1998); Wiedemeier and Csillag (1979); Chattopadhyay et al. (1986). The -polarized lowest-energy transverse optic soft-mode was also shown to continuously condense across in our previous INS study Li et al. (2015a). The phase transition temperature, K, obtained from our heat capacity measurements is in good agreement with values (802–813 K) reported in the literature Adouby et al. (1998); Wiedemeier and Csillag (1979); Chattopadhyay et al. (1986). A slight hysteresis in of about 4 K was observed in our DSC. Some of the hysteresis may possibly be caused by some Se evaporation at high , oxidation, or temperature lags in the instrument. However, it is also possible for the transition to exhibit a partially first-order character, since there is some reported evidence for a small latent heat Feutelais et al. (1996); Sharma and Chang (1986); Balde et al. (1981). We point out that the values of from our DSC measurements (Fig. 12), while in excellent agreement with those of Sassi et al. Sassi et al. (2014), are significantly larger than the linear estimate from laser flash measurements reported in Ref. 18, which also misses the phase transition behavior. For example, at 700 K, our measurement is 16% larger than the linear estimate of Zhao et al. Zhao et al. (2014), and at 790 K the discrepancy approximately reaches a factor of two. However, the discrepancy is more minimal above the phase transition, in the Cmcm phase.

We calculated the harmonic and dilational heat capacity from INS phonon DOS and thermal expansion coefficient measurements. The harmonic ()and dilational () heat capacity are given by,

(19) |

(20) |

respectively Wallace (1972); Fultz (2010). To calculate the harmonic heat capacity from Eq. (19), we have used the non-neutron-weigthed version of phonon DOS at 10 K, as described earlier for the entropy calculations. The comparison of experimental measurements with harmonic and dilational heat capacity is shown in Fig. 12. For K, the quasi-harmonic (harmonic+dilational) heat capacity accounts well for the measured . However, as the phase transition is approached, the difference between the two curves increases significantly, reflecting the growing contribution of anharmonicity. A similar behavior was also seen in entropy, internal energy and free energy (Fig. 8, 9, and 10).

## Vi Conclusion

We have investigated the anharmonic thermodynamics and negative thermal expansion in SnSe with a combination of INS, NRIXS, calorimetry measurements, and first-principles simulations. We identified a pronounced contribution of anharmonicity at high temperature, especially within K of the structural phase transition. The NTE along the –axis, the direction of corrugation of Sn-Se bilayers in the Pnma phase, can be qualitatively accounted for with a quasiharmonic free energy minimization, but deviations arising from strong anharmonicity are evident at high temperatures. The structural distortion accompanying the NTE can be traced to rotations of SnSe tetrahedra, which overlap with the strongly anharmonic soft-modes condensing at the phase transition.

## Vii Acknowledgements

We would like to acknowledge Ayman H. Said for technical help with NRIXS measurements at APS sector 30. We also thank Amr Mohammed for help in preparing DSC samples and measurements. Data analysis, modeling, and phonon simulations (D.B, J.H.) were supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, through the Office of Science Early Career Award grant of O.D. Neutron scattering measurements were supported as part of the S3TEC EFRC, an Energy Frontier Research Center funded by the US Department of Energy, Office of Science, Basic Energy Sciences under Award # DE-SC0001299 (C.W.L, O.D.). Sample synthesis (A.F.M.) was supported by the US Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division. The use of Oak Ridge National Laboratory’s Spallation Neutron Source was sponsored by the Scientific User Facilities Division, Office of Basic Energy Sciences, US Department of Energy. Use of the APS was supported by DOE-BES under Contract No. DE-AC02-06CH11357. Theoretical calculations were performed using resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the US Department of Energy under contract no.DE-AC02-05CH11231. This research used resources of the Oak Ridge Leadership Computing Facility, which is supported by the Office of Science of the U.S. DOE.

## References

- Snyder and Toberer (2008) G. Snyder and E. Toberer, Nature Materials 7, 105 (2008).
- Zebarjadi et al. (2012) M. Zebarjadi, K. Esfarjani, M. Dresselhaus, Z. Ren, and G. Chen, Energy Environmental Science 5, 5147 (2012).
- Wallace (1972) D. Wallace, Thermodynamics of Crystals (John Wiley & Sons Inc., 1972).
- Wallace (2003) D. Wallace, Statistical Physics of Crystals and Liquids: A Guide to Highly Accurate Equations of State (World Scientific Pub Co Inc., Singapore, 2003).
- Grimvall (1999) G. Grimvall, Thermophysical Properties of Materials (North Holland; 1 edition, Amsterdam, The Netherlands, 1999).
- Fultz (2010) B. Fultz, Progress in Materials Science 55, 247 (2010).
- Rhyee et al. (2009) J.-S. Rhyee, K.-H. Lee, S.-M. Lee, E. Cho, S.-I. Kim, E. Lee, Y. Kwon, J. Shim, and G. Kotliar, Nature 459, 965 (2009).
- Nielsen et al. (2013) M. Nielsen, V. Ozolins, and J. Heremans, Energy Environmental Science 6, 570 (2013).
- Lee et al. (2014) S. Lee, K. Esfarjani, T. Luo, J. Zhou, Z. Tian, and G. Chen, Nature Communications 5, 3525 (2014).
- Li et al. (2015a) C. Li, J. Hong, A. May, D. Bansal, J. Ma, T. Hong, S. Chi, G. Ehlers, and O. Delaire, Nature Physics 11, 1063 (2015a).
- Delaire et al. (2011a) O. Delaire, J. Ma, K. Marty, A. May, M. McGuire, M.-H. Du, D. Singh, A. Podlesnyak, G. Ehlers, M. Lumsden, and B. Sales, Nature Materials 8, 614 (2011a).
- Shiga et al. (2012) T. Shiga, J. Shiomi, J. Ma, O. Delaire, T. Radzynski, A. Lusakowski, K. Esfarjani, and G. Chen, Physical Review B 85, 155203 (2012).
- Delaire et al. (2011b) O. Delaire, K. Marty, M. Stone, P. Kent, M. Lucas, D. Abernathy, D. Mandrus, and B. Sales, Proceedings of the National Academy of Sciences 108, 4725 (2011b).
- Delaire et al. (2008a) O. Delaire, M. S. Lucas, J. A. Muñoz, M. Kresch, and B. Fultz, Physical Review Letters 101, 105504 (2008a).
- Delaire et al. (2008b) O. Delaire, M. Kresch, J. A. Muñoz, M. S. Lucas, J. Y. Y. Lin, and B. Fultz, Physical Review B 77, 214112 (2008b).
- Bansal et al. (2015) D. Bansal, C. Li, A. Said, D. Abernathy, J.-Q. Yan, and O. Delaire, Physical Review B 92, 214301 (2015).
- Hong and Delaire (2016) J. Hong and O. Delaire, arXiv , 1604.07077 (2016).
- Zhao et al. (2014) L.-D. Zhao, S.-H. Lo, Y. Zhang, H. Sun, G. Tan, C. Uher, C. Wolverton, V. Dravid, and M. Kanatzidis, Nature 508, 373 (2014).
- Sassi et al. (2014) S. Sassi, C. Candolfi, J. Vaney, V. Ohorodniichuk, P. Masschelein, A. Dauscher, and B. Lenoir, Applied Physics Letters 104, 212105 (2014).
- Chen et al. (2014) C.-L. Chen, H. Wang, Y.-Y. Chen, T. Day, and G. J. Snyder, J. Mater. Chem. A 2, 11171 (2014).
- Zhao et al. (2016) L.-D. Zhao, G. Tan, S. Hao, J. He, Y. Pei, H. Chi, H. Wang, S. Gong, H. Xu, V. P. Dravid, C. Uher, G. J. Snyder, C. Wolverton, and M. G. Kanatzidis, Science 351, 141 (2016).
- Chattopadhyay et al. (1986) T. Chattopadhyay, J. Pannetier, and H. Von-Schnering, Journal of Physics and Chemistry of Solids 47, 879 (1986).
- Adouby et al. (1998) K. Adouby, C. Perez-Vicente, J. C. Jumas, R. Fourcade, and A. A. Toure, Z. Kristallogr 213, 343 (1998).
- Wiedemeier and Csillag (1979) H. Wiedemeier and F. Csillag, Zeitschrift für Kristallographie 149, 17 (1979).
- Chattopadhyay et al. (1984) T. Chattopadhyay, A. Werner, H. Von-Schnering, and J. Pannetier, Revue de Physique Appliquee 19, 807 (1984).
- von Schnering and Wiedemeier (1981) H. von Schnering and H. Wiedemeier, Zeitschrift für Kristallographie 156, 143 (1981).
- Feutelais et al. (1996) Y. Feutelais, M. Majid, B. Legendre, and S. G. Fries, Journal of Phase Equilibria 17, 40 (1996).
- Sharma and Chang (1986) R. C. Sharma and Y. A. Chang, Bulletin of Alloy Phase Diagrams 7, 68 (1986).
- Balde et al. (1981) L. Balde, B. Legendre, C. Souleau, and P. Khodadad, Journal of Less-Common Metals 80, 45 (1981).
- Dembovskii et al. (1963) S. A. Dembovskii, B. N. Egorov, A. S. Pashinkin, and Y. A. Polyakov, Russian Journal of Inorganic Chemistry 8, 530 (1963).
- Zhdanova (1961) V. V. Zhdanova, Soviet Physics Solid State 3, 1174 (1961).
- Carrete et al. (2014) J. Carrete, N. Mingo, and S. Curtarolo, Applied Physics Letters 105, 101907 (2014).
- Ding et al. (2015) G. Ding, G. Gao, and K. Yao, Scientific reports 5, 9567 (2015).
- Loa et al. (2015) I. Loa, R. J. Husband, R. A. Downie, S. R. Popuri, and J.-W. G. Bos, J. Phys.: Condens. Matter 27, 072202 (2015).
- Loa et al. (2016) I. Loa, R. J. Husband, R. A. Downie, S. R. Popuri, and J.-W. G. Bos, J. Mater. Chem. C 4, 1685 (2016).
- Guo et al. (2015) R. Guo, X. Wang, Y. Kuang, and B. Huang, Physical Review B 92, 115202 (2015).
- Kutorasinski et al. (2015) K. Kutorasinski, B. Wiendlocha, S. Kaprzyk, and J. Tobola, Physical Review B 91, 205201 (2015).
- Wei et al. (2015) T.-R. Wei, C.-F. Wu, X. Zhang, Q. Tan, L. Sun, Y. Pan, and J.-F. Li, Phys. Chem. Chem. Phys. 17, 30102 (2015).
- Sanchez et al. (2015) F. S. Sanchez, M. Gharsallah, N. M. Nemes, F. J. Mompean, J. L. Martínez, and J. A. Alonso, Applied Physics Letters 106, 083902 (2015).
- Li et al. (2015b) Y. Li, X. Shi, D. Ren, J. Chen, and L. Chen, Energies 8, 6275 (2015b).
- Han et al. (2015) Y. M. Han, J. Zhao, M. Zhou, X. X. Jiang, H. Q. Leng, and L. F. Li, J. Mater. Chem. A 3, 4555 (2015).
- Ge et al. (2015) Z. H. Ge, K. Wei, H. Lewis, J. Martin, and G. S. Nolas, Journal of Solid State Chemistry 225, 354 (2015).
- Stone et al. (2014) M. Stone, J. Niedziela, D. Abernathy, L. DeBeer-Schmitt, G. Ehlers, O. Garlea, G. Granroth, M. Graves-Brook, A. Kolesnikov, A. Podlesnyak, and B. Winn, Review of Scientific Instruments 85, 045113 (2014).
- Toellner et al. (2011) T. S. Toellner, A. Alatas, and A. H. Said, Journal of Synchrotron Radiation 18, 605 (2011).
- Sturhahn (2000) W. Sturhahn, Hyperfine Interactions 125, 149 (2000).
- Kresse and Hafner (1993) G. Kresse and J. Hafner, Physical Review B 47, 558 (1993).
- Kresse and Furthmller (1996a) G. Kresse and J. Furthmller, Physical Review B 54, 11169 (1996a).
- Kresse and Furthmller (1996b) G. Kresse and J. Furthmller, Computational Materials Science 6, 15 (1996b).
- Perdew and Zunger (1981) J. Perdew and A. Zunger, Physical Review B 23, 5048 (1981).
- Perdew et al. (1996) J. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
- A. Togo et al. (2008) A. A. Togo, F. Oba, and I. Tanaka, Physical Review B 78, 134106 (2008).
- Brugger (1965) K. Brugger, Physical Review 137, 17 (1965).
- Cantrell (1980) J. Cantrell, Physical Review B 21, 4191 (1980).
- Sturhahn and A.Chumakov (1999) W. Sturhahn and A.Chumakov, Hyperfine Interactions 123, 809 (1999).
- Lovesey (1984) S. Lovesey, Theory of neutron scattering from condensed matter (Clarendon Press, Oxford, 1984).
- He et al. (2013) X. He, H. Shen, W. Wang, Z. Wang, B. Zhang, and X. Li, Journal of Alloys and Compounds 556, 86 (2013).
- Bansal et al. (2016) D. Bansal, A. J. Aref, G. F. Dargush, and O. Delaire, Journal of Physics: Condensed Matter 28, 385201 (2016).
- Buckingham and Fairbank (1961) M. J. Buckingham and W. M. Fairbank, Chapter 3: The nature of -transition in liquid helium – Progress in Low Temperature Physics, Vol. 3 (North Holland, Amsterdam, 1961) pp. 80–112.