Thermal conductivity of MoS₂ monolayers from molecular dynamics simulations

MD simulations for thermal conductivity require the use of empirical forcefields that can reproduce phonon populations, phonon and sound velocities and phonon mean free paths accurately. Similar to previous studies,^24,25,31,32 we parameterize a SWFF,³³ which models covalent interatomic interactions in the MoS₂ 2H crystal structure using 2-body and 3-body terms that mimic energy profiles for bond-stretching and bond-bending (i.e. bond-angle distortion) interactions. The total potential energy of the given system of N atoms located at [r₁, r₂, …, r_N] in the SWFF can be written as

where r_ij = |r_j−r_i|. The 2-body term, V₂, is defined as

The two-body term is defined by 3 optimizable parameters, A, B and γ. The 3-body term, V₃ around a central atom i is given by

where λ, γ₁ and γ₂ are optimizable parameters. The exponential decay of interatomic energies and the explicitly-defined cutoff distance for interatomic interactions makes the forcefield short-ranged and local, which makes for easier numerical implementation for large-scale MD simulations.

All SWFF parameters (A, B, γ, K, γ₁ and γ₂) are optimized for the MoS₂ crystal structure using the GULP program.³⁴ The remaining geometric parameters, namely, the three cutoff distances, r_cut, r₁ and r₂ and θ₀ are obtained directly from ab initio simulations (See Section V of supplementary material). θ₀ is the included angle between Mo-S-Mo and S-Mo-S triplets and r_cut for each pair of atoms is set equal to 1.5 times the nearest-neighbor distance between those types of atoms in the ideal 2H MoS₂ crystal structure. 2-body terms are parameterized for all type of atom-pairs: Mo-Mo, Mo-S and S-S. 3-body terms are parameterized only for the Mo-S-Mo and S-Mo-S. Other non-bonded triplets do not contribute to the potential energy of the crystal. The SWFF is optimized to reproduce three quantities, namely the lattice constant, in-plane elastic constants and phonon dispersion curves. During optimization, GULP uses the BFGS scheme to minimize the weighted residuals for these three observables.

where i loops over the lattice constant, the elastic constant and phonon dispersion and f^DFT is the DFT-calculated value for observable i and f^SW(x) is the value obtained from the Stillinger Weber forcefield with the parameter set x.

For the purpose of residual evolution and forcefield optimization, the phonon dispersion is discretized at 21 points for each band along the Γ-M-K-Γ direction. Further, for an accurate calculation of thermal conductivity, it is desirable that the empirical forcefield accurately reproduces the behavior of low-energy acoustic phonon modes which are responsible for the bulk of heat conduction at temperatures relevant to this study (T = 100 K – 400 K).¹⁷ Therefore, in our optimization scheme, the three acoustic branches have a larger weighting factor, w_i = 30 for the residual evaluation, compared to w_i = 1 for the high energy optical modes, which have negligible contributions to thermal transport.¹⁷ 

NEMD simulations are performed with the new SWFF (Tables I and II) using the LAMMPS molecular dynamics program.³⁵ The simulations are performed on a rectangular monolayer MoS₂ single-crystal of dimensions L × 2L with periodic boundary conditions along the x- and y-directions. Atomic positions and simulation cell parameters of the MoS₂ monolayer are initially relaxed using the SWFF until there no net forces on the atoms and no residual stress on the simulation cell. The system is then heated under the NPT ensemble to the desired target temperature where thermal conductivity is to be measured. Next, as defined in Figure 1, two thin (20Å wide) regions denoted ‘Hot’ and ‘Cold’, which span the entire width, L, of the simulation cell are defined at x = 3L/2 and x = L/2. To establish a thermal steady state, a fixed energy flux, $Q̇$⁠, is added to the kinetic energy of atoms in the ‘Hot’ region and an identical energy flux is removed the kinetic energy of atoms in the ‘Cold’ region. This method^36,37 allows us to impose arbitrarily small thermal fluxes, in contrast to the momentum-swapping method of reverse non-equilibrium molecular dynamics,³⁸ which can lead to large temperature gradients spanning ∼100 K and non-linear temperature profiles near the heat source and sink.

At steady state, this technique establishes a gradient in local temperature between L/2 and 3L/2 with an approximate and small temperature gradient of 10^-2 K/Å, which falls within the linear regime for Fourier law of heat conduction.²⁴ In these simulations, local (i.e. per-atom) temperature is defined in a classical sense based on the equation $12mvi2=32kBTi$⁠, where the brackets denote averages over time. Alternative definitions of local temperature based on quantification of phonon populations are infeasible due to the large phonon mean free path of MoS₂, which is larger than the size of the simulation cells used in this study.³⁹ 

NEMD simulations are performed for 24 ns within the NVE ensemble, and the temperature profile in the system is obtained from the last 12 ns of the simulation, which ensures that a steady state has been established. Thermal conductivity of the simulation cell is obtained from the Fourier law for heat conduction

where κ is the thermal conductivity of the MoS₂ monolayer and ∇T is the temperature gradient in the steady-state from the source to the sink due to the energy flux $Q̇$⁠. L and t are respectively the width and effective thickness of the MoS₂ monolayer and their product is the cross-sectional area through which heat conduction occurs. The factor ½ comes from the geometry of the simulation system, where heat conduction from the ‘Hot’ to ‘Cold’ region occurs along the +x and -x directions. It is important to note that thermal conductivity calculated from Fourier’s law of heat conduction requires an estimation of the effective thickness, t, through which the heat flux, $Q̇$⁠, is transported. There are several competing definitions for what constitutes the thickness of a two-dimensional crystal.⁴⁰ In this study, we have chosen the effective thickness of the MoS₂ monolayer to be equal to 6.16 Å, the experimentally calculated interlayer distance in bulk and multilayer MoS₂ samples,⁴¹ consistent with previous studies.^17,25

Thermal transport in NEMD occurs over two square regions of size L×L. The dimensions of this artificial system prevent the participation of larger mean-free-path phonons in the heat conduction process. Therefore, thermal conductivity values obtained from NEMD simulations show a strong dependence on the size of the simulated system, with larger systems resulting in larger calculated values of κ due to the participation of a larger subset of phonons in heat transfer. All NEMD simulations in this study are performed in the Casimir limit,⁴² i.e. in the regime where the system size is smaller than the mean free path, Λ, of acoustic phonons (estimated to be between 17.2 nm and 2200 nm for MoS₂^16,17). Therefore, to obtain the material-intrinsic (i.e. size-independent) thermal conductivity of the system, we follow the scaling scheme used by Oligschleger and Schon,⁴³ which invokes Matthiessen’s rule of independent scattering events,⁴⁴ to extrapolate κ to a (hypothetical) system of infinite size using the relation below.

where 1/l is the inverse scattering length due to the source-sink distance and 1/l₀ is the material-intrinsic inverse scattering length due to phonon-phonon scattering. The intrinsic size-independent thermal conductivity κ₀ of the material can be obtained at the limit 1/l → 0. It is important to note that such linear extrapolation schemes in the Casimir regime are extremely sensitive to the size of the systems used for extrapolation.⁴⁵ 

System size dependence of thermal conductivity is evaluated by calculating thermal conductivity values from NEMD simulations on crystals of size L×2L where L is one of 19, 30, 40, 60, 144, 200 and 300 nm. Figure 2(a) shows the dependence of the computed thermal conductivity at 300 K on simulation cell sizes, reflecting increasing thermal conductivity values for larger system sizes. Figure 2(b) demonstrates the system size scaling method based on Matthiessen’s rule and extracts an intercept value of $1/κMD300K=0.019WmK−1$⁠. The large system sizes used in this study up to 300 nm ensure that linear extrapolation leads to bulk thermal conductivity predictions that lie within 20% of those predicted by lattice dynamics.

This procedure is repeated at different temperatures, 100 K, 200 K, 300 K and 400 K to obtain κ_MD(T) (Figure 3(a)). These NEMD-calculated thermal conductivity values show a monotonic decrease with increasing temperatures due to greater phonon population and correspondingly greater phonon-phonon scattering rates (Figure 3(b)). This trend of decreasing thermal conductivity with increasing temperature is consistent with experimental measurements only at high temperatures above the Debye temperature, θ_D. One of the most significant drawbacks of NEMD-based calculations of thermal conductivity is the completely classical nature of the simulation, which is applicable only to solids above θ_D.³⁹ Therefore, calculated κ_MD values for T > θ_D ∼ 263 K⁴⁶ are realistic, but extrapolations of this method to low temperatures (∼100 K)²⁴ are not justified.

Below θ_D, the quantization of vibrational energy becomes the major source of error in NEMD-calculated thermal conductivity. Specifically, lattice thermal conductivity of crystals is linearly proportional to the lattice specific heat, C_v.

where v and Λ are the group velocity and mean free path of phonons respectively. The specific heat in classical NEMD simulations is constant at $CvMD=3NkB$⁠, at all temperatures, which corresponds to the high temperature limit of experimental specific heat values, $Cvq(T)$⁠. Therefore, the non-classical system-level correction to the MD thermal conductivity, κ_MD(T), can be obtained by substituting the classical value of specific heat in κ_MD(T) with the experimentally realistic value of $Cvq(T)$⁠.

This empirical system-level correction accounts for the relative occupation of different vibrational energy levels and provides a good approximation of κ_MD at temperatures T < θ_D. More rigorous, phonon-mode-specific corrections have been proposed for crystalline systems, which additionally correct for different relaxation times and phonon line-widths at different temperatures,^47,48 but these corrections are not attempted in our work.

The quantized vibrational energy of the crystal is given by the equation

where G(ω) is the phonon density of states and BE(ω, T) is the temperature-dependent Bose-Einstein distribution given by $BEω,T=expħωkBT−1−1$⁠. The lattice specific heat is then given by the equation

Figure 4(a) shows the computed C_v, which approaches the classical value of 3k_B at high temperatures. Figure 4(b) compares κ_MD(T) with the quantum-corrected κ_corr(T) showing significantly reduced values below θ_D.

Phonon scattering in experimental MoS₂ samples come from four sources.

• Phonon-phonon scattering, which controls the temperature-dependence of thermal conductivity

• Scattering by isotopes

• Scattering by point defects, primarily sulfur vacancies⁴⁹ 

• Scattering by crystal edges, boundaries and terminations

For the case of material-intrinsic thermal conductivity (i.e. for a defect-free, infinitely large crystal), the last two terms vanish, leaving phonon-phonon and phonon-isotope scattering as the primary mechanisms controlling thermal conductivity. In order to quantify the impact of phonon scattering by cation and anion isotopes, we compare the thermal conductivity of a crystal containing an explicit distribution of atomic isotopes at their natural concentration, with a model crystal where the mass of all atoms of a given species are set equal to its abundance-weighted average of atomic masses of all isotopes. Table S4 lists the isotopes of Molybdenum and Sulfur along with their naturally occurring abundances.^50–52 As Figure 5 demonstrates, the lack of phonon scattering by discrete isotopes leads to a nearly 45% greater value for thermal conductivity of the model crystal with homogeneous atomic masses. This behavior is qualitatively similar to that observed in thermal conductivity experiments on diamond and uranium, where isotopically pure samples show significantly larger thermal conductivity.^53,54 However, the magnitude of thermal conductivity suppression due to isotopes is significantly larger in the case of MoS₂ due primarily to the polydispersity of isotopic distribution in Molybdenum and Sulfur compared to Carbon and the relatively larger normalized isotope mass distribution in Molybdenum compared to Uranium.¹⁵ This effect of isotopic distribution is more pronounced at low temperatures, where the other scattering mechanism, phonon-phonon scattering, is less effective.

In conclusion, we have performed NEMD simulations with a newly-parameterized empirical Stillinger Weber forcefield to calculate the thermal conductivity of suspended monolayer MoS₂ as a function of temperature in the range of 100 K – 400 K. NEMD-calculated thermal conductivity values were corrected to account for the suppression of specific heat at low temperatures as well as phonon scattering by explicit isotopes in the crystal. The reported thermal conductivity values will serve as a reference for monolayer MoS₂ systems.

See supplementary material for details about forcefield parameterization and calculation of vibrational properties.

This work was supported as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award Number DE-SC0014607. Simulations were performed at the Argonne Leadership Computing Facility under the DOE INCITE program and at the Center for High Performance Computing at the University of Southern California.