Distinct anharmonic characteristics of phonon-driven lattice thermal conductivity and thermal expansion in bulk MoSe₂ and WSe₂

Abstract

Using inelastic neutron scattering and X-ray diffraction measurements, together with ab initio and machine-learning molecular dynamics simulations, we bring out the distinct nature of anharmonicity in the phonon spectra of MoSe₂ and WSe₂ relevant to thermal transport and thermal expansion behaviour. We show that the perturbation method, including 4^th-order force constants, is insufficient to capture the temperature-dependent explicit anharmonicity. The Green–Kubo method captures the explicit anharmonicity and reproduces the thermal conductivity (κ_l) with high fidelity. Our mode-resolved calculation reveals that the major contribution (∼90%) to κ_l is attributed to a small explicit anharmonicity of low-energy phonons. In contrast, these modes exhibit large positive Grüneisen parameters (implicit anharmonicity), causing the large thermal expansion of the material.

---

Introduction

Low dimensional materials such as transition-metal dichalcogenides (TMDCs) are known to exhibit interesting electronic, optical and thermodynamic properties and are promising materials for applications across thermoelectrics, transistors, optoelectronics, topological materials, solid lubrication, electrode materials for Mg and Zn batteries, energy storage and photovoltaic devices.^1–16 Many TMDCs, such as WSe₂, WS₂, MoSe₂ and MoS₂, are considered promising candidates for thermoelectric applications, and comprehensive studies have investigated their electronic, optical and transport properties.^2–6,17,18 In particular, WSe₂ has been found to exhibit ultralow thermal conductivity in single-layer form.¹⁹ These 2D materials are also promising candidates for flexible electronic devices and nano-electromechanical systems due to their ultralow weight, flexibility and high tensile strength.^7,8,20–22

Limited experimental and theoretical studies exist on the thermal conductivity of TMDCs, and they show discrepancies between different measurements and simulations.^10,23–26 P. Jiang et al. measured the anisotropic thermal conductivity of MX₂ (M = Mo, W and X = S, Se) using the time-domain thermoreflectance technique, from which they obtained in-plane (out-of-plane) conductivities ∼35 W m⁻¹ K⁻¹ (2.5 W m⁻¹ K⁻¹) and 40 W m⁻¹ K⁻¹ (2.2 W m⁻¹ K⁻¹), in MoSe₂ and WSe₂, respectively, at room temperature.²³ S. Kumar et al. used first-principles calculations and found that the three lowest optical modes for bulk MoSe₂ and WSe₂ contribute to the lattice thermal conductivity almost as much as the acoustic modes.¹⁰ The authors reported in-plane (out-of-plane) thermal conductivity values at a room temperature of ∼57 W m⁻¹ K⁻¹ (∼6 W m⁻¹ K⁻¹) and ∼48 W m⁻¹ K⁻¹ (∼6 W m⁻¹ K⁻¹), for MoSe₂ and WSe₂, respectively.¹⁰ Another theoretical investigation of thermal conductivity by Lindroth et al. on WSe₂ and MoSe₂ found in-plane and out-of-plane conductivities ∼33 W m⁻¹ K⁻¹ and 4 W m⁻¹ K⁻¹, respectively, for both compounds.²⁴ These calculations were performed with boundary scattering contribution considering a sample size of ∼1 μm. There are many other TMDCs, for which thermal conductivity measurements and simulations show significant variance in their values.^23,25–30 Hence, a comprehensive theoretical study exploring different techniques is informative to clarify the underlying physics of the thermal conductivity in these TMDCs and the possible source of error/variance.

Besides the specifics of TMDCs, a microscopic understanding of thermal properties is generally important to design materials for applications. The anharmonic characteristics of phonons are crucial in both thermal transport and thermal expansion behaviour.^31,32 For instance, Grüneisen parameters, which estimate phonon energy shifts due to anharmonicity upon expanding volume, are needed to account for thermal expansivity, while the phonon lifetimes, reflecting anharmonic phonon interactions control thermal conductivity.³³ Hence, predicting these quantities accurately and rationalizing their microscopic origins can help design materials with tailored properties. Several methods are widely used to compute the lattice thermal conductivity of solids: solving the phonon Boltzmann transport equation (BTE) within perturbation theory,^34–36 computing the heat-flux auto-correlation function with equilibrium molecular dynamics,³⁷ or computing the thermal gradient or heat flow using non-equilibrium molecular dynamics,^38–40 homogeneous non-equilibrium molecular dynamics,^41,42 atomic green function methods,^43,44etc. While the perturbation and Green–Kubo methods are among the most popular approaches, all these techniques have certain limitations and advantages.

The lowest-order perturbation theory of the linearized phonon BTE uses the second and third-order force-constants computed around an equilibrium configuration (at T = 0 K) and treats finite temperatures by considering the respective phonon occupation numbers. An extension of this method consists in using renormalized second and third-order force-constants at finite temperature, as obtained from ab initio molecular dynamics (AIMD) simulations, which can extend accuracy to a broader temperature range. However, even at low-T, these methods do not perform well for strongly anharmonic solids where high-order interactions are sizeable. Classical molecular dynamics simulations intrinsically capture all orders of anharmonicity in the interatomic potential, but do not capture the quantum effect at low temperatures. Hence, this approach tends to work best in moderate and high-temperature regimes.

In this manuscript, we calculated the thermal conductivity of MoSe₂ and WSe₂ by different techniques. We find that the perturbation methods are not satisfactory, and the machine-learning molecular dynamics (MLMD) simulations quantitatively well capture the higher-order anharmonicity. In addition, we measured the temperature dependence of the phonon density of states (DOS) from 313 K to 1098 K using inelastic neutron scattering (INS), which is used to benchmark our simulations and assess the strength of the anharmonicity. Furthermore, we also measured the thermal expansion by X-ray diffraction (XRD) and calculated the same within the quasi-harmonic approximation (QHA). By combining the experiments and simulations, we found that the contribution of the low-energy (low-E) phonons dominates the heat transport as well as the thermal expansion. These low-E phonons show the right magnitude of the phonon lifetimes and the Gruneisen parameters, which quantitatively reproduce the experimental thermal properties. The details about the sample synthesis, INS and XRD measurements, and state-of-the-art computational techniques of the temperature-dependent effective potential (TDEP) and MLMD are given in the ESI.†⁴⁵

Experimental

A polycrystalline MoSe₂ sample has been prepared from a stoichiometric ratio of high purity Mo (in powder form) and Se in small shots. To start with, the Mo powders and selenium shots were thoroughly ground in an agate mortar. The mixed powders were then cold-pressed into pellets. The pellets were then sealed under 10⁻⁶ mbar vacuum in a quartz ampoule. The quartz ampoule was then loaded into a resistive heating box-type furnace. The temperature of the furnace was gradually raised at a rate of 15 °C h⁻¹ up to 900 °C and kept at this temperature for about 48 h. Then, the furnace was cooled at a rate of about 20 °C h⁻¹ down to room temperature. We purchased the WSe₂ polycrystalline sample from Sigma Aldrich. Powder XRD was performed to characterize the samples and their phase purity. Temperature-dependent X-ray diffraction studies were carried out using an 18 kW rotating Cu anode-based powder diffractometer operating in the Bragg–Brentano focusing geometry with a curved crystal monochromator. Data were collected in the continuous scan mode at a scan speed of 1° per minute and step interval of 0.02°.

The inelastic neutron scattering measurements were carried out using a high-flux time-of-flight (IN4C) spectrometer at the Institut Laue Langevin (ILL), France, covering a wide range of scattering angles from 10° to 110°. Thermal neutrons of wavelength 2.4 Å (14.2 meV) were used for the measurements. The orientation-averaged scattering function, S(Q,E), was measured in the neutron energy gain mode with momentum transfer, Q, extending up to 7 Å⁻¹. The polycrystalline samples were loaded into a cylindrical niobium sample holder and mounted in a furnace. INS data were collected with samples heated to 313 K, 498 K, 673 K and 898 K. The neutron-weighted phonon DOS, g⁽ⁿ⁾(E) was extracted from the measured S(Q, E) using the incoherent one-phonon approximation as follows:

where A is a normalization constant and 2W(Q) is the Debye–Waller factor. The multiphonon contribution at each temperature was calculated using the Sjolander⁴⁶ formalism and subtracted.

In simulations, one can also compute g⁽ⁿ⁾(E) from the calculated partial phonon DOS, g_k(E), as follows:

where B is a normalization constant, m_k, b_k and g_k(E) are the mass, neutron scattering length and partial phonon DOS of the k^th type of atom, respectively. The weighting factors for W, Mo, and Se are 0.025, 0.060, and 0.105 in units of barns per amu, respectively. The neutron scattering lengths for various atoms can be found in ref. 47.

Computational details

Ab initio molecular dynamics

Ab initio molecular dynamics simulations were performed within the density functional theory (DFT) framework implemented in the Vienna Ab initio Simulation Package (VASP).^48,49 The simulations used the projector-augmented wave formalism⁵⁰ and the generalized gradient approximation exchange–correlation functional parameterization by Perdew, Burke, and Ernzerhof.^51,52 The DFT-D3 method with the Becke–Johnson damping function was used to account for the van der Waals dispersive interactions between WSe₂/MoSe₂ layers along the c-axis.^53,54 Under ambient conditions, WSe₂ and MoSe₂ crystallize in a hexagonal structure (space group no. 194, P6₃/mmc) that contains six atoms in the unit-cell (2 metal atoms and 4 Se) (Fig. S1†⁴⁵). A 4 × 4 × 1 supercell (96 atoms) was used for AIMD simulations. For DFT calculations, we used a 2 × 2 × 2 k-point mesh generated using the Monkhorst–Pack method⁵⁵ and a plane-wave energy cut-off of 800 eV. Before AIMD simulations, we fully relaxed the crystal structure with an energy and Hellmann–Feynman force convergence thresholds of 10⁻⁸ eV and 10⁻⁶ eV Å⁻¹, respectively, and obtained the optimized lattice constants a and c shown in Table S1.†⁴⁵ The total energy convergence criteria were set to 10⁻⁶ eV. The simulations were performed from 100 K to 1000 K in intervals of 100 K within the NVT framework, and the temperatures were controlled using a Nose–Hoover thermostat.⁵⁵ All the simulations were run for 10–15 ps with a time step of 2 fs.

Machine learning molecular dynamics (MLMD)

The DEEP-MD code was used to generate a machine-learning surrogate force field based on a neural network representation.⁵⁶ Subsequent MLMD simulations were performed with this DEEP-MD potential using LAMMPS.⁵⁷ The generated force field was tested against several quantities also computed from AIMD (see Fig. S2†⁴⁵). The phonon spectral energy density (SED) was computed based on MLMD on a 20 × 20 × 10 supercell (24 000 atoms) up to 100 ps, which provides a Q- and E-resolution of 0.05 rlu and 0.05 meV, respectively.

Lattice thermal conductivity from MLMD

To compute the thermal conductivity within the Green–Kubo formalism, we performed MLMD simulations on a 24 × 24 × 5 supercell (17 280) up to 2 ns with 1 fs time steps at 100 K, 300 K, 500 K and 700 K for both WSe₂ and MoSe₂. The heat-flux vector, J(t) is given by:where ε_i and v_i are the energy and velocity of the i^th particle and F_ij is the force between the i^th and j^th atoms. The first and second terms in eqn (3) represent the convection and conduction terms of heat flux, respectively. Convection results from the mass transport of mobile (diffusive) species, while the second term describes the energy exchange due to interatomic interactions, e.g., lattice vibrations in crystalline solids. In non-diffusive solids, the convection term is insignificant. The total lattice thermal conductivity tensor (κ^uv_l) was calculated using the Green–Kubo formalism,³⁹ within which it is expressed as a time integral of the heat current autocorrelation function. In this approach, the lattice thermal conductivity tensor element κ^uv_l (u,v = x,y,z) is given by:where J_u(t) is the heat-flux component along the u direction at time t.

Lattice thermal conductivity simulation from perturbation theory

The lattice thermal conductivity (κ_l) using the perturbation method was calculated using 2^nd, 3^rd and 4^th order force constants using ShengBTE. These force-constants were evaluated at 0 K. To investigate the effect of temperature on force-constants and subsequently on κ_l, we have computed κ_l using 2^nd and 3^rd-order renormalized force constants (FCs) at 300 K implemented in TDEP.^58,59

The thermal conductivity, κ_l is expressed as:

where C_vλ, υ_gλ and τ_λ are the heat capacity, group velocity, and lifetime, which correspond to a specific phonon mode λ, respectively. λ consists of q and j, which are the phonon wave vector and branch number.

The C_vλ and υ_gλ can be determined from harmonic phonons, while τ_λ can be determined by calculating the phonon scattering rate due to three or higher order phonon–phonon scattering as well as from other scattering processes. Here, we only consider the phonon–phonon scattering. The harmonic interatomic force constants (IFCs) are calculated using a 4 × 4 × 1 supercell and a 2 × 2 × 2 Monkhorst–Pack grid within a finite displacement scheme, as implemented in the open-source software package Phonopy.⁶⁰ The third and fourth-order IFCs were calculated through third order and fourth order packages of ShengBTE, considering up to the third-nearest neighbours. The same supercell size and k-point mesh were used to calculate the third and fourth IFCs by the finite difference method. With these IFCs, the phonon frequencies and velocities are then obtained by diagonalizing the dynamical matrix, and the three-phonon scattering rates are calculated from the iterative solution of the BTE, while four-phonon scatterings were calculated within the relaxation time approximation. The thermal conductivity is solved on a 24 × 24 × 5 q-mesh, sufficient to integrate the BZ (Fig. S3†⁴⁵). The expressions for all involved scattering rates are described in ref. 34 and 36.

Spectral energy density (SED)

The phonon SED at wavevector q⃑ and energy E were obtained from MLMD trajectories using the following formalism:⁶¹where N is the number of unit cells in a supercell (N = N₁ × N₂ × N₃), and summation index α runs over Cartesian x, y, and z; index k runs over the number of particles in the unit cell. m_k and are the mass of the k^th atom and its equilibrium position in the n^th unit cell, and is the velocity of the k^th atom in the n^th unit cell at time t.

Results and discussion

X-ray diffraction, inelastic neutron scattering and phonon dynamics

We have performed the temperature-dependent X-ray diffraction measurements on polycrystalline samples of MoSe₂ and WSe₂ from 20 K to 280 K, as shown in Fig. 1. There are no impurity peaks in the XRD pattern, confirming the single-phase nature of the sample. The XRD peaks are indexed as per the data reported in Powder Diffraction File ICDD JCPDS-PDF 04-004-8782 and 00-038-1388 for hexagonal (space group P6₃/mmc) MoSe₂ and WSe₂, respectively.

Fig. 1 The temperature dependence of XRD patterns of MoSe₂ (ICDD JCPDS-PDF 04-004-8782) and WSe₂ (ICDD JCPDS-PDF 00-038-1388). All the peaks are indexed using hexagonal symmetry with space group P6₃/mmc.

The INS measurements were carried out at five temperatures: 313 K, 498 K, 673 K, 898 K, and 1098 K (Fig. 2 and S4†⁴⁵). The phonon spectra of MoSe₂ and WSe₂ extend up to 45 and 40 meV, respectively. The temperature-dependent spectra show several peaks at around 5, 14, 16, 19, 22, 26, and 36 meV in MoSe₂ and 5, 13, 15, 24, 27 and 32 meV in WSe₂. From our partial DOS calculations, we found that the entire spectral range (0–45 meV) is contributed by both Mo/W and Se dynamics. Interestingly, even at 1098 K, the low-E (<10 meV) phonon spectra do not show any significant change, while above 10 meV, the phonon peaks clearly soften. The low-E portion of the spectra is largely attributed to phonons from acoustic branches and the three lowest optical branches. The softening and broadening of phonon modes are attributed to lattice expansion (implicit anharmonicity) and phonon–phonon interactions (explicit anharmonicity). The implicit anharmonic shifts largely contribute to the thermal expansion, while the explicit anharmonicity dominates thermal transport. Hence, we have computed both anharmonic components in the entire BZ to identify the modes relevant to thermal expansion and lattice thermal conductivity.

Fig. 2 (a and d) Temperature dependence of the measured neutron-weighted PDOS g⁽ⁿ⁾(E) of MoSe₂ and WSe₂; (b and e) MLMD simulations of g⁽ⁿ⁾(E) of MoSe₂ and WSe₂; (c and f) the low-energy part of S(E) spectra of MoSe₂ and WSe₂ measured with INS, with the simulated S(E) shown in the inset.

We performed NVT-MLMD simulation to explicitly probe the temperature evolution of the neutron-weighted phonon spectra (Fig. 2(b) and (e)). We observe that the simulated low-E modes do not show significant change with temperature, while the high energy modes slightly soften at elevated temperatures as experimentally observed, although less in magnitude than in the measured spectra. The AIMD simulated phonon DOS as a function of temperatures are also shown in Fig. S7.†

To further probe the temperature behaviour of low-E modes, we plot the measured Q-integrated dynamical structure factor, S(E) (1.0 Å⁻¹ < Q < 4.0 Å⁻¹) in Fig. 2(c) and (f). This plot clearly shows that the modes barely shift upon warming but visibly broaden. We expect the explicit anharmonicity due to the phonon–phonon interactions and increased phonon amplitudes and the implicit anharmonicity due to the volume changes. Our detailed calculations show that the negligible shift of the low-E modes occurs due to a cancellation of the hardening due to the fourth and higher order explicit anharmonicity and the softening due to the implicit anharmonicity of these modes.

To investigate the effect of third-order and higher-order anharmonicity on phonon spectra, we computed the phonon spectral energy density (SED) along various high-symmetry directions in the Brillouin zone (BZ) using two different approaches. The first approach is based on the TDEP method, which computes the phonon SED using third-order force-constants only (referred to here as TDEP-SED), and the second approach is based on molecular dynamics trajectories (NVT-MLMD)⁶¹ and referred to here as MLMD-SED. The latter approach includes third and higher-order anharmonicity to phonon SED. Hence, comparing the phonon SED from the two methods helps us identify the region in the BZ dominated by fourth or higher-order anharmonicity. In Fig. 3, we plot TDEP-SED and MLMD-SED at 100 K and 500 K. TDEP-SED at 100 K shows sharp features in intensity profiles and does not change significantly at 500 K in both compounds. However, while MLMD-SED also shows sharp phonon peaks at 100 K, these become significantly broader at 500 K, making MLMD-SED look more diffuse compared to TDEP-SED.

Fig. 3 (Top panel) The TDEP-calculated phonon SED (including up to third-order anharmonicity) at 100 K and 500 K of MoSe₂ and WSe₂. (Bottom-panel) The MLMD-calculated phonon SED (including all orders of anharmonicity) at 100 K and 500 K of MoSe₂ and WSe₂. The red markers in WSe₂ SED are the phonon energies measured at 300 K by Cai et al.⁸⁰

Lattice thermal conductivity and anharmonic phonons

In the harmonic picture of lattice dynamics, the phonons are non-interacting and infinitely long-lived quasiparticles and hence the infinite lattice thermal conductivity. However, phonons in real systems get scattered by other phonons (intrinsic scattering processes), leading to a finite lifetime of phonons. Crystal defects and a finite sample size further reduce the phonon lifetimes, and the shortened mean free paths limit the phononic thermal conductivity (extrinsic scattering processes). In ideal single-crystals, the finite lattice thermal conductivity results from intrinsic phonon scattering processes, namely Umklapp and normal processes. Generally, normal processes do not directly contribute to the thermal resistance. Umklapp processes are the leading phonon scattering processes in heat conduction at high temperatures. Furthermore, intrinsic phonon scattering processes may involve three or more phonons. Three-phonon scattering processes include two phonons annihilating and creating a new phonon or one phonon splitting into two phonons. This is generally the dominant scattering channel in crystalline materials and leads to a 1/T dependence of κ_l. However, some materials such as BAs, perovskites, Heusler alloys and others^62–70 show limited scattering space for three-phonon processes and significantly higher scattering phase space for four-phonon processes, and as a result do not exhibit the typical 1/T dependence of κ_l, with a lower conductivity than three-phonon processes alone.

From INS measurements and MLMD simulations, we find that the phonons in both compounds do not change significantly at elevated temperatures; hence phase space and group velocity would not show any significant change. However, phonons substantially broaden on heating, corresponding to decreased lifetimes. The temperature dependence of thermal conductivity in these compounds is therefore primarily governed by phonon lifetimes. We first calculated the lattice thermal conductivity using third-order force constants at 0 K using ShengBTE (κ^Sheng3). In this approach, the phonon lifetime is only contributed by three-phonon scattering processes. In Fig. 4(a), (b) and S5,† we have shown the three-phonon scattering rates, life times and group velocities of MoSe₂ and WSe₂ at 300 K. The lattice thermal conductivity of MoSe₂ and WSe₂ was measured by Puqing Jiang et al.²³ using a time-domain thermoreflectance approach and found large anisotropy between the in-plane and out-of-plane lattice thermal conductivity. In Fig. 4(c)–(f), we have shown the calculated thermal conductivity from three-phonon processes in the ab-plane and along the c-axis and also shown the experimental values.²³ However, three-phonon scattering has shown a significant overestimation of thermal conductivity compared to experimental values (κ^expt). Including the four-phonon scattering rates improves this discrepancy, and we examined their impact on the thermal conductivity and phonon lifetimes using ShengBTE (κ^Sheng4). Here also, the fourth-order force constants were evaluated at 0 K. We found that four-phonon scattering resistance significantly reduces the thermal conductivities and phonon lifetimes but is still overestimated compared to experimental values (Fig. 4). This suggests that either higher-order processes dominate and/or that the third and fourth-order force constants at 0 K are insufficient. To further probe these aspects, we also used the renormalized second and third-order force constants to compute the thermal conductivity at respective temperatures using TDEP (κ^TDEP). Surprisingly, we find that the κ^TDEP at low-T shows lower thermal conductivity values than κ^Sheng4; however, they are nearly the same at higher temperatures (Fig. 4). This indicates that the 0 K anharmonic force constants cannot explain the measured thermal conductivity, and one also needs to consider the temperature dependence of higher-order force-constants. However, it is currently not feasible to evaluate the temperature dependence of fourth-order force constants and their effect on thermal conductivity due to TDEP limitations.

Fig. 4 (a and b) Calculated mode-resolved phonon scattering rates on a (24 × 24 × 7) q-point grid in the BZ in MoSe₂ and WSe₂ at 300 K using 3-phonon scattering, calculated from 2^nd and 3^rd-order force-constants at 0 K (ShengBTE-3ph), and renormalized 2^nd and 3^rd-order force-constants at 300 K (TDEP-3ph). The scattering rates are also calculated using 3- and 4-phonon scattering processes using 0 K 3^rd and 4^th order force-constants (ShengBTE 3ph + 4ph). (c and d) The calculated in-plane lattice thermal conductivity (κ_xy) and (e and f) out-of-plane lattice thermal conductivity (κ_z), using 3-phonon scattering from ShengBTE-3ph, TDEP-3ph and ShengBTE 3ph + 4ph methods and the Green–Kubo method based on equilibrium molecular dynamics simulation, includes all orders of anharmonicity. The calculated quantities are compared with available experimental measurements.²³ The numerical values of thermal conductivity from different methods (contribution of 4^th and higher order interactions) are given in Table S3.†

To further analyze the full impact of anharmonicity on κ, we performed MLMD simulations and evaluated the thermal conductivity within the Green–Kubo formalism, κ^GK. Strikingly, the κ^GK from MLMD reproduces the measured in-plane thermal conductivity (κ^expt_xy) very well and over a wide range of temperatures. However, the out-of-plane component (κ^expt_z) is somewhat overestimated. It is important to note that these TMDCs have a weak interplanar interaction and are governed by van der Waals forces. In the DFT framework, precisely capturing the van der Waals interactions is difficult and may lead to an overestimated thermal conductivity along the c-axis. Hence, the extensive simulation based on the perturbation method and MLMD Green–Kubo approach indicates that the importance of the temperature effect on higher-order force-constants is critical in understanding and predicting the thermal conductivity. To determine which phonon modes are most critical in heat transport in the ab-plane and along the c-axis, we computed the mode-resolved in-plane and out-of-plane lattice thermal conductivity as a function of phonon energy at 300 K (Fig. 5). Interestingly, more than 90% of lattice thermal conductivity is attributed to modes below 20 meV. By comparing the different methods, it appears that the higher-order scattering processes and temperature-renormalised force constants are critical for precisely evaluating these low-E phonon lifetimes and, hence, the thermal conductivity of TMDCs. The effect of higher-ordered scattering processes and their temperature dependence was also observed in strongly anharmonic superionics and 2-d materials, including WS₂ and graphene, where the impact of higher-order anharmonicity on thermal conductivity is found to be significant.^71–79

Fig. 5 (a and b) The calculated in-plane (κ_xy) and (c and d) out-of-plane (κ_z) lattice thermal conductivity contributed from individual modes in the BZ (on a 24 × 24 × 5 grid) in MoSe₂ and WSe₂ using (i) 3-phonon processes (lowest order perturbation method) at 0 K (Sheng-BTE 3ph), (ii) 3 and 4 phonon processes at 0 K (Sheng-BTE 3 + 4ph) and (iii) 3-phonon processes with a renormalized force constant at 300 K (TDEP 3ph). (e and f) The calculated mode Grüneisen parameters in MoSe₂ and WSe₂ using pressure dependence of phonon frequencies; the inset shows the fractional change in volume with temperature calculated within the QHA and compared with our X-ray measurements.

Lattice thermal expansion behaviour and phonon anharmonicity

Thermal expansion is also associated with anharmonicity on the potential energy surface. By knowing the phonon frequency shifts with temperature and volume, one can estimate how the material would behave at a given temperature. The phonon frequency shifts due to changes in the vibrational amplitudes with temperature are usually much smaller than the frequency shifts due to the change in volume with temperature.³² We have calculated the phonon frequencies at different volumes using lattice dynamics within the quasi-harmonic approximation (QHA).

In the quasi-harmonic approximation, the volume thermal expansion coefficient⁴² of a crystalline material is given by the following relation: . Here, is the mode Grüneisen parameter, which is a measure of the pressure response of the phonon energy to determine the nature of the expansion behaviour. For instance, the phonon modes with negative values of Γ tend to contract the lattice, while the positive ones would expand the lattice. C_Vi(T) is the specific heat contribution of the i^th phonon mode (of energy E_i) at temperature T, while B and V stand for the bulk modulus and volume of the compound, respectively. The calculated elastic constants and bulk modulus are given in Table S2.†⁴⁵ The pressure dependence of the phonon spectra within QHA has been calculated in the entire Brillouin zone to allow for the calculation of the energy dependence of the Grüneisen parameter Γ(E) and further processed to obtain the thermal expansion coefficient α_V(T).

The volume dependence of the phonon spectra within the QHA has been calculated in the entire Brillouin zone to allow for the calculation of the mode-wise Grüneisen parameter Γ(q_j,E) (Fig. 4(e) and (f)) and to obtain the thermal expansion coefficient α_V(T). We have calculated the mode-wise thermal expansion contribution on a 20 × 20 × 10 q-point grid in the BZ at 300 K, as shown in Fig. S6.†⁴⁵ It is evident that the most significant contribution to thermal expansion comes from low-E modes. We have performed temperature-dependent X-ray diffraction measurements on MoSe₂ and WSe₂ from 50 K to 300 K. The temperature dependence of fractional change in volume with respect to 50 K is shown in Fig. 5(e) and (f), which exhibits positive expansion. The large implicit anharmonicity of the low-E modes leads to large thermal expansion. It is interesting to note that the same low-E modes exhibit rather small explicit anharmonicity and contribute most to thermal conductivity.

Conclusions

We combined inelastic neutron scattering measurements and ab initio and machine-learning molecular dynamics simulations to unravel the phonon anharmonicity and their role in thermal transport and thermal expansion in MoSe₂ and WSe₂. We find that low-energy phonons play a crucial role in determining both the thermal-expansion behaviour and heat transport in these materials. Conventional perturbation methods are insufficient to capture the temperature-dependent anharmonicity of these phonon modes, particularly the low-E modes in strongly anharmonic systems. Machine-learning molecular dynamics simulations pave the way to completely treat the anharmonicity of phonons by performing temperature-dependent large-scale molecular dynamics simulations. These results establish the importance of anharmonic low-energy modes in TMDCs, which drive the thermal and transport properties.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The use of the ANUPAM super-computing facility at BARC is acknowledged. SLC acknowledges the financial support of the Indian National Science Academy for the INSA Senior Scientist position. OD acknowledges funding from the U.S. Department of Energy (DOE) BES MSED, under Award No. DE-SC0019978.

References

1. S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev and A. Kis, Nat. Rev. Mater., 2017, 2, 17033.
2. B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti and A. Kis, Nat. Nanotechnol., 2011, 6, 147–150.
3. Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman and M. S. Strano, Nat. Nanotechnol., 2012, 7, 699–712.
4. X. Hong, J. Kim, S.-F. Shi, Y. Zhang, C. Jin, Y. Sun, S. Tongay, J. Wu, Y. Zhang and F. Wang, Nat. Nanotechnol., 2014, 9, 682–686.
5. M. Massicotte, P. Schmidt, F. Vialla, K. G. Schädler, A. Reserbat-Plantey, K. Watanabe, T. Taniguchi, K. J. Tielrooij and F. H. L. Koppens, Nat. Nanotechnol., 2016, 11, 42–46.
6. H. Liang, Z. Cao, F. Ming, W. Zhang, D. H. Anjum, Y. Cui, L. Cavallo and H. N. Alshareef, Nano Lett., 2019, 19, 3199–3206.
7. K. Nassiri Nazif, A. Daus, J. Hong, N. Lee, S. Vaziri, A. Kumar, F. Nitta, M. E. Chen, S. Kananian, R. Islam, K.-H. Kim, J.-H. Park, A. S. Y. Poon, M. L. Brongersma, E. Pop and K. C. Saraswat, Nat. Commun., 2021, 12, 7034.
8. D. B. Velusamy, R. H. Kim, S. Cha, J. Huh, R. Khazaeinezhad, S. H. Kassani, G. Song, S. M. Cho, S. H. Cho, I. Hwang, J. Lee, K. Oh, H. Choi and C. Park, Nat. Commun., 2015, 6, 8063.
9. W. Huang, H. Da and G. Liang, J. Appl. Phys., 2013, 113, 104304.
10. S. Kumar and U. Schwingenschlögl, Chem. Mater., 2015, 27, 1278–1284.
11. H.-P. Komsa, J. Kotakoski, S. Kurasch, O. Lehtinen, U. Kaiser and A. V. Krasheninnikov, Phys. Rev. Lett., 2012, 109, 035503.
12. G. Wang, C. Robert, M. M. Glazov, F. Cadiz, E. Courtade, T. Amand, D. Lagarde, T. Taniguchi, K. Watanabe, B. Urbaszek and X. Marie, Phys. Rev. Lett., 2017, 119, 047401.
13. F. Wu, T. Lovorn, E. Tutuc, I. Martin and A. H. MacDonald, Phys. Rev. Lett., 2019, 122, 086402.
14. D. Li, H. Shan, C. Rupprecht, H. Knopf, K. Watanabe, T. Taniguchi, Y. Qin, S. Tongay, M. Nuß, S. Schröder, F. Eilenberger, S. Höfling, C. Schneider and T. Brixner, Phys. Rev. Lett., 2022, 128, 087401.
15. Z. Lin, Y. Liu, Z. Wang, S. Xu, S. Chen, W. Duan and B. Monserrat, Phys. Rev. Lett., 2022, 129, 027401.
16. Y.-M. Wu, Z. Wu and H. Yao, Phys. Rev. Lett., 2023, 130, 126001.
17. H. Guo, T. Yang, P. Tao, Y. Wang and Z. Zhang, J. Appl. Phys., 2013, 113, 013709.
18. D. Sarkar, W. Liu, X. Xie, A. C. Anselmo, S. Mitragotri and K. Banerjee, ACS Nano, 2014, 8, 3992–4003.
19. C. Chiritescu, D. G. Cahill, N. Nguyen, D. Johnson, A. Bodapati, P. Keblinski and P. Zschack, Science, 2007, 315, 351–353.
20. S. Manzeli, D. Dumcenco, G. Migliato Marega and A. Kis, Nat. Commun., 2019, 10, 4831.
21. L. Gao, Small, 2017, 13, 1603994.
22. W. Choi, N. Choudhary, G. H. Han, J. Park, D. Akinwande and Y. H. Lee, Mater. Today, 2017, 20, 116–130.
23. P. Jiang, X. Qian, X. Gu and R. Yang, Adv. Mater., 2017, 29, 1701068.
24. D. O. Lindroth and P. Erhart, Phys. Rev. B, 2016, 94, 115205.
25. J. J. Bae, H. Y. Jeong, G. H. Han, J. Kim, H. Kim, M. S. Kim, B. H. Moon, S. C. Lim and Y. H. Lee, Nanoscale, 2017, 9, 2541–2547.
26. X. Gu, B. Li and R. Yang, J. Appl. Phys., 2016, 119, 085106.
27. R. Yan, J. R. Simpson, S. Bertolazzi, J. Brivio, M. Watson, X. Wu, A. Kis, T. Luo, A. R. Hight Walker and H. G. Xing, ACS Nano, 2014, 8, 986–993.
28. J. Judek, A. P. Gertych, M. Świniarski, A. Łapińska, A. Dużyńska and M. Zdrojek, Sci. Rep., 2015, 5, 12422.
29. X. Zhang, D. Sun, Y. Li, G.-H. Lee, X. Cui, D. Chenet, Y. You, T. F. Heinz and J. C. Hone, ACS Appl. Mater. Interfaces, 2015, 7, 25923–25929.
30. I. Jo, M. T. Pettes, E. Ou, W. Wu and L. Shi, Appl. Phys. Lett., 2014, 104, 201902.
31. V. A. Drebushchak, J. Therm. Anal. Calorim., 2020, 142, 1097–1113.
32. R. Mittal, M. K. Gupta and S. L. Chaplot, Prog. Mater. Sci., 2018, 92, 360–445.
33. X. Qian, J. Zhou and G. Chen, Nat. Mater., 2021, 20, 1188–1202.
34. Z. Han, X. Yang, W. Li, T. Feng and X. Ruan, Comput. Phys. Commun., 2022, 270, 108179.
35. A. Togo, L. Chaput and I. Tanaka, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 91, 094306.
36. W. Li, J. Carrete, N. A. Katcho and N. Mingo, Comput. Phys. Commun., 2014, 185, 1747–1758.
37. A. J. McGaughey and M. Kaviany, Adv. Heat Transfer, 2006, 39, 169–255.
38. M. S. Green, J. Chem. Phys., 1954, 22, 398–413.
39. R. Kubo, Rep. Prog. Phys., 1966, 29, 255.
40. R. Kubo, M. Yokota and S. Nakajima, J. Phys. Soc. Jpn., 1957, 12, 1203–1211.
41. F. Müller-Plathe, J. Chem. Phys., 1997, 106, 6082–6085.
42. Z. Fan, H. Dong, A. Harju and T. Ala-Nissila, Phys. Rev. B, 2019, 99, 064308.
43. J.-S. Wang, J. Wang and N. Zeng, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 74, 033408.
44. N. Mingo and L. Yang, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 68, 245406.
45. See the ESI† for force-field validation, sample preparation, numerical data of thermal conductivities, phonon DOS, group velocity and lifetime calculations.
46. A. Sjolander, Ark. Fys., 1958, 14, 315.
47. Neutron Scattering, ed. D. Price and K. Skold, Academic Press, Orlando, 1986, vol. 35.
48. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
49. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
50. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
51. J. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1998, 80, 891.
52. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
53. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
54. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
55. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Solid State, 1976, 13, 5188.
56. H. Wang, L. Zhang, J. Han and E. Weinan, Comput. Phys. Commun., 2018, 228, 178–184.
57. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
58. O. Hellman, P. Steneteg, I. A. Abrikosov and S. I. Simak, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 104111.
59. O. Hellman and I. A. Abrikosov, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 144301.
60. A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1–5.
61. J. A. Thomas, J. E. Turney, R. M. Iutzi, C. H. Amon and A. J. McGaughey, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 081411.
62. T. Yue, Y. Zhao, J. Ni, S. Meng and Z. Dai, npj Comput. Mater., 2023, 9, 17.
63. Y. Xiao, Y. Zhao, J. Ni, S. Meng and Z. Dai, Mater. Today Commun., 2023, 35, 105450.
64. X. Song, J. Wang, Y. Zhao, J. Ni, S. Meng and Z. Dai, Phys. Lett. A, 2022, 456, 128550.
65. Y. Zhao, S. Zeng, G. Li, C. Lian, Z. Dai, S. Meng and J. Ni, Phys. Rev. B, 2021, 104, 224304.
66. C. Zhang, J. Sun, Y. Shen, W. Kang and Q. Wang, J. Phys. Chem. Lett., 2022, 13, 5734–5741.
67. T. Feng, L. Lindsay and X. Ruan, Phys. Rev. B, 2017, 96, 161201.
68. K.-C. Zhang, C. Shen, H.-B. Zhang, Y.-F. Li and Y. Liu, Phys. Rev. B, 2022, 106, 235202.
69. N. K. Ravichandran and D. Broido, Phys. Rev. X, 2020, 10, 021063.
70. F. Tian, B. Song, X. Chen, N. K. Ravichandran, Y. Lv, K. Chen, S. Sullivan, J. Kim, Y. Zhou, T.-H. Liu, M. Goni, Z. Ding, J. Sun, G. A. G. Udalamatta Gamage, H. Sun, H. Ziyaee, S. Huyan, L. Deng, J. Zhou, A. J. Schmidt, S. Chen, C.-W. Chu, P. Y. Huang, D. Broido, L. Shi, G. Chen and Z. Ren, Science, 2018, 361, 582–585.
71. A. She, Y. Zhao, J. Ni, S. Meng and Z. Dai, Int. J. Heat Mass Transfer, 2023, 209, 124132.
72. G. Sun, J. Ma, C. Liu, Z. Xiang, D. Xu, T.-H. Liu and X. Luo, Int. J. Heat Mass Transfer, 2023, 215, 124475.
73. G.-H. Liu, Z.-X. Xie, P.-Z. Jia, X.-J. Wu and X.-K. Chen, Diamond Relat. Mater., 2023, 137, 110116.
74. G. Zhang, S. Dong, C. Yang, D. Han, G. Xin and X. Wang, Appl. Phys. Lett., 2023, 123, 052205.
75. J. Tang, G. Li, Q. Wang, J. Zheng, L. Cheng and R. Guo, Int. J. Heat Mass Transfer, 2023, 207, 124011.
76. J. Han, Q. Zeng, K. Chen, X. Yu and J. Dai, Nanomaterials, 2023, 13, 1576.
77. Q. Ren, M. K. Gupta, M. Jin, J. Ding, J. Wu, Z. Chen, S. Lin, O. Fabelo, J. A. Rodríguez-Velamazán, M. Kofu, K. Nakajima, M. Wolf, F. Zhu, J. Wang, Z. Cheng, G. Wang, X. Tong, Y. Pei, O. Delaire and J. Ma, Nat. Mater., 2023, 22, 999–1006.
78. M. K. Gupta, J. Ding, D. Bansal, D. L. Abernathy, G. Ehlers, N. C. Osti, W. G. Zeier and O. Delaire, Adv. Energy Mater., 2022, 12, 2200596.
79. M. K. Gupta, J. Ding, N. C. Osti, D. L. Abernathy, W. Arnold, H. Wang, Z. Hood and O. Delaire, Energy Environ. Sci., 2021, 14, 6554–6563.
80. Q. Cai, B. Wei, Q. Sun, A. H. Said and C. Li, Mater. Today Phys., 2022, 28, 100856.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3ta03830k

---