Phonon anharmonicity and negative thermal expansion in SnSe

Phonon anharmonicity and negative thermal expansion in SnSe

Dipanshu Bansal, Jiawang Hong, Chen W. Li, Andrew F. May, Wallace Porter, Michael Y. Hu, Douglas L. Abernathy, and Olivier Delaire bansald@ornl.gov olivier.delaire@duke.edu Materials Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA
Advanced Photon Source, Argonne National Laboratory, Argonne, Illinois 60439, USA
Quantum Condensed Matter Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA
Mechanical Engineering and Materials Science, Duke University, Durham, North Carolina 27708, USA
July 12, 2019
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.

Figure 1: Crystal structure of SnSe in Pnma (a) and Cmcm (c) phases, illustrating the double bilayer structure and Pnma distortion, and corresponding Sn-Se bonding state (b,d). Sn atoms are in grey and Se atoms in green. and are labels for bonds in out-of-plane and in-plane directions. The orientation of crystal axes in the Cmcm structure is chosen to match the Pnma phase, in order to facilitate the comparison. The rotation of pyramids (shaded in (a)) and the angle between two bonds in the Pnma structure are also shown in (a). The conventional unit cell edges are indicated by black square box.

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
Table 1: Lattice constants and internal atomic coordinates for Pnma phase from DFT simulations compared with experimental data at 300 K. Internal coordinates are normalized to the lattice constants. The number in parentheses are percentage deviation from the experimental data reported in Ref. 23. Details are provided in the text.
Figure 2: Phonon dispersions of SnSe from DFT simulations with either LDA or GGA exchange-correlation functionals (relaxed unit cells), compared with experimental INS data at 300 K.

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
Table 2: Anisotropic Grüneisen parameter () of tin-selenide calculated by expanding and compressing the lattice along , , and axes. The Brillouin-zone integrated values are averages of over the entire Monkhorst-Pack q-point mesh, while the values along specific segments of the Brillouin zone (e.g. ) are averages over polarizations only. Further details are provided in the text.

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
Table 3: Fourth order isothermal bulk modulus tensor calculated from DFT simulations and experimental inelastic neutron scattering (INS) data. Experimental INS data is at T = 300 K. More details in text.

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.

Figure 3: (a) Lattice parameters calculated (open markers) from minimization of the free energy surface at different temperatures, compared with experimental data (solid markers) from Adouby et al. Adouby et al. (1998). (b) Grid of isothermal strains computed with DFT to sample the free-energy surface.

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.

Figure 4: Phonon density of states calculated from DFT for three different configurations of lattice parameters. The ground-state equilibrium lattice parameters are . The other two configurations are and .

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
Table 4: Rotation of pyramidal induces NTE along . Bond length and lattice constant in . All data is from Ref. Adouby et al. (1998) except for the calculated value and . The small discrepancy between calculated and results from the ignorance of tilt of two  bonds plane.

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.

Figure 5: Temperature dependence of bond length (a) and the rotation angle of triangular pyramids (b) extracted from structure data in Ref. 23. (c) Solids dots is the lattice parameters in Ref. 23, the curves are derived from Eq. 11 and 12. (d) The lattice change with temperature, the change is divided by the tilt of bond () and rotation of bond (), according to Eq. 11, is the lattice constant at 813 K from data in Ref. 23.

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 .

Figure 6: Neutron-weighted phonon DOS of measured with inelastic neutron scattering at different temperatures for incident neutron energies of a) 55 meV and b) 30 meV, compared with neutron weighted and experimental resolution convoluted DFT simulations.
Figure 7: Generalized phonon DOS of calculated with DFT, weighted with experimental neutron cross-sections. The phonon DOS is convolved with a damped harmonic oscillator function, as detailed in the text.

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.

Figure 8: Components of vibrational entropy, , of at T = 10, 150, 300, 450, 600, and 750 K. Lower panel shows the markers of total () and harmonic (red line) vibrational entropy contribution along with values calculated from experimental heat capacity data (cyan line), while upper panel displays the markers for nonharmonic (+), dilational (magenta line), and anharmonic entropy (). More details in the text.
Figure 9: Components of vibrational internal energy, , of at T = 10, 150, 300, 450, 600, and 750 K. Lower panel shows the markers of total () and harmonic (red line) vibrational internal energy contribution, while upper panel displays the markers for nonharmonic (+), dilational (magenta line), and anharmonic internal energy (). More details in the text.
Figure 10: Components of vibrational free energy, , of at T = 10, 150, 300, 450, 600, and 750 K. Lower panel shows the markers of total () and harmonic (red line) vibrational free energy contribution, while upper panel displays the markers for nonharmonic (+), dilational (magenta line), and anharmonic free energy (). More details in the text.

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
Table 5: Grüneisen parameter () of tin-selenide calculated by applying 1% isothermal strain along . Average along particular direction such as represents the values of averaged over q-point along that direction weighted by specific heat capacity. More details in text.

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.

Figure 11: a) Partial phonon DOS of tin and selenium calculated from DFT simulations for , the ground state equilibrium lattice parameters for fully relaxed structure, b) projected phonon DOS of Sn parallel to axis, and c) projected phonon DOS of Sn parallel to and axis compared with partial phonon DOS of Sn from NRIXS measurements at room temperature.
(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
Table 6: Mean square thermal displacement parameter of tin and selenium calculated from DFT (LDA) simulations at T = 300 K compared with experimental neutron diffraction and NRIXS data.

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. 89, and 10).

Figure 12: Heat capacity of SnSe measured with DSC on single-crystal samples (both heating and cooling), and its harmonic, dilational, and nonharmonic components, indicating the strong anharmonicity around the phase transition. The green dots are data measured on polycrystalline samples by Sassi et al. Sassi et al. (2014). The cyan triangles are data measured on single crystal samples by Zhao et al. Zhao et al. (2014)

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.
Comments 0
Request Comment
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 minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
101618
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
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