Quantification of lattice thermal conductivity of two-dimensional semiconductors like MoS_{2} is necessary for the design of electronic and thermoelectric devices, but direct experimental measurements on free-standing samples is challenging. Molecular dynamics simulations, with appropriate corrections, can provide a reference value for thermal conductivity for these material systems. Here, we construct a new empirical forcefield of the Stillinger-Weber form, parameterized to phonon dispersion relations, lattice constants and elastic moduli and we use it to compute a material-intrinsic thermal conductivity of 38.1 W/m-K at room temperature and estimate a maximum thermal conductivity of 85.4 W/m-K at T = 200 K. We also identify that phonon scattering by the large isotopic mass distribution of Mo and S contributes a significant correction (>45%) to the thermal conductivity at low temperatures.

In-plane thermal transport is an important consideration for the design of nanoelectronic and thermoelectric devices made from layered and two-dimensional (2D) materials which possess some of the lowest out-of-plane thermal boundary conductance for solid-solid interfaces.^{1} Lattice thermal conductivity also determines the thermal profile of growing two-dimensional crystals in exothermic synthesis techniques like chemical vapor deposition and is an important metric for the fabrication of energy harvesting devices based on the thermoelectric effect.^{2–4} Despite this importance, definitive experimental measurements of fundamental thermal properties of these materials like lattice thermal conductivity have not been conclusive. The primary experimental challenges, as detailed in recent reviews,^{2,5} stem from difficulties in accurately measuring temperature gradients on atomic length scales in suspended monolayer samples. Further, thermal conductivity of suspended monolayers and sub-nm-thin materials is greatly quenched (by up to 85%) when these monolayers are supported on a substrate or are investigated using contact-probes.^{6,7} Even state-of-the art non-contact techniques for thermal conductivity measurement like opto-thermal Raman are undermined by several implicit assumptions about optical absorption and equilibrium of phonon modes^{8} and are affected by heat losses to the environment^{9} leading to large errors in reported values of thermal conductivity.^{10–12} Further, such techniques involving localized laser heating can only be used to investigate room-temperature and high-temperature thermal conductivity and cannot readily access low-temperature regimes.^{13} In contrast, atomistic MD simulations are an attractive method for estimating temperature-dependent thermal conductivity of suspended monolayer systems directly from nm-scale temperature gradients. Specifically, classical non-equilibrium molecular dynamics (NEMD) provides a computationally inexpensive way of computing lattice thermal conductivity for monolayer materials, which includes the full anharmonicity of interatomic forces within a steady-state (i.e. dynamic and non-equilibrium) simulation. Classical NEMD simulations offers several advantages over the calculation of lattice thermal conductivity by *ab initio* atomistic simulations. Solutions to the phonon Boltzmann transport equation from *ab initio* simulations involve strong assumptions about phonon scattering rates and lifetimes as well as boundary conditions for phonon scattering.^{14–19} These assumptions are reflected in the large scatter in the reported values of calculated thermal conductivity for 2D materials like MoS_{2}^{16,19–25} (Table S1). Further, thermal conductivity values from classical NEMD simulations are consistent with values extracted from classical equilibrium MD (EMD) simulations using Green-Kubo techniques,^{26} but have a smaller magnitude of systematic errors^{25} (Table S2). However, classical MD simulations ignore the quantization of vibrational energy levels, which is responsible for sub-classical specific heat at low temperatures and correspondingly suppressed thermal conductivity. MD simulations in literature also fail to account for additional phonon scattering introduced by the naturally-occurring distribution of isotopic masses. Here, we report quantum-corrected and isotope-corrected thermal conductivity values for a representative two-dimensional crystal, MoS_{2}, calculated from NEMD simulations using a new forcefield parameterized to reproduce experimental lattice constants, elastic moduli and *ab initio* calculated phonon dispersion curves. MoS_{2} is an ideal material for thermal conductivity calculations since despite its widespread use in nano-electronic applications, there is a large scatter in reported values of experimentally determined thermal conductivity ranging from 34 W/m-K to 84 W/m-K at 300 K.^{27–30}

## PARAMETERIZATION OF STILLINGER WEBER FORCEFIELD (SWFF)

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,^{33} which models covalent interatomic interactions in the MoS_{2} 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**_{1}, **r**_{2}, …, **r**_{N}] in the SWFF can be written as

where *r*_{ij} = |*r*_{j}−*r*_{i}|. The 2-body term, V_{2}, is defined as

The two-body term is defined by 3 optimizable parameters, *A*, *B* and *γ*. The 3-body term, V_{3} around a central atom *i* is given by

where *λ*, *γ*_{1} and *γ*_{2} 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*, *γ*_{1} and *γ*_{2}) are optimized for the MoS_{2} crystal structure using the GULP program.^{34} The remaining geometric parameters, namely, the three cutoff distances, *r*_{cut}, *r*_{1} and *r*_{2} and *θ*_{0} are obtained directly from *ab initio* simulations (See Section V of supplementary material). *θ*_{0} 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_{2} 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).^{17} 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.^{17}

### Non-equilibrium molecular dynamics simulations

NEMD simulations are performed with the new SWFF (Tables I and II) using the LAMMPS molecular dynamics program.^{35} The simulations are performed on a rectangular monolayer MoS_{2} single-crystal of dimensions *L* × 2*L* with periodic boundary conditions along the *x*- and *y*-directions. Atomic positions and simulation cell parameters of the MoS_{2} 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* = 3*L*/2 and *x* = *L*/2. To establish a thermal steady state, a fixed energy flux, $Q\u0307$, 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,^{38} which can lead to large temperature gradients spanning ∼100 K and non-linear temperature profiles near the heat source and sink.

Interacting Atoms . | A . | ρ . | B . | r_{min}
. | r_{max}
. |
---|---|---|---|---|---|

S-S | 1.117 | 0.369 | 38.564 | 0.00 | 4.3463 |

Mo-S | 7.098 | 0.680 | 9.971 | 0.00 | 3.2033 |

Mo-Mo | 3.777 | 0.653 | 26.272 | 0.00 | 4.3463 |

Interacting Atoms . | A . | ρ . | B . | r_{min}
. | r_{max}
. |
---|---|---|---|---|---|

S-S | 1.117 | 0.369 | 38.564 | 0.00 | 4.3463 |

Mo-S | 7.098 | 0.680 | 9.971 | 0.00 | 3.2033 |

Mo-Mo | 3.777 | 0.653 | 26.272 | 0.00 | 4.3463 |

Triplet . | λ . | θ_{0}
. | γ_{1} = γ_{2}
. | r_{12}^{max}
. | r_{13}^{max}
. | r_{23}^{max}
. |
---|---|---|---|---|---|---|

S-Mo-S | 7.663 | 82.4795 | 0.870 | 3.2034 | 3.2034 | 4.3463 |

Mo-S-Mo | 30.230 | 82.4795 | 3.482 | 3.2034 | 3.2034 | 4.3463 |

Triplet . | λ . | θ_{0}
. | γ_{1} = γ_{2}
. | r_{12}^{max}
. | r_{13}^{max}
. | r_{23}^{max}
. |
---|---|---|---|---|---|---|

S-Mo-S | 7.663 | 82.4795 | 0.870 | 3.2034 | 3.2034 | 4.3463 |

Mo-S-Mo | 30.230 | 82.4795 | 3.482 | 3.2034 | 3.2034 | 4.3463 |

At steady state, this technique establishes a gradient in local temperature between *L*/2 and 3*L*/2 with an approximate and small temperature gradient of 10^{-2} K/Å, which falls within the linear regime for Fourier law of heat conduction.^{24} 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_{2}, which is larger than the size of the simulation cells used in this study.^{39}

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_{2} monolayer and ∇*T* is the temperature gradient in the steady-state from the source to the sink due to the energy flux $Q\u0307$. *L* and *t* are respectively the width and effective thickness of the MoS_{2} 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\u0307$, is transported. There are several competing definitions for what constitutes the thickness of a two-dimensional crystal.^{40} In this study, we have chosen the effective thickness of the MoS_{2} monolayer to be equal to 6.16 Å, the experimentally calculated interlayer distance in bulk and multilayer MoS_{2} samples,^{41} consistent with previous studies.^{17,25}

### System size dependence in NEMD simulations

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,^{42} 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_{2}^{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,^{43} which invokes Matthiessen’s rule of independent scattering events,^{44} 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*_{0} is the material-intrinsic inverse scattering length due to phonon-phonon scattering. The intrinsic size-independent thermal conductivity *κ*_{0} 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.^{45}

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/\kappa MD300K=0.019WmK\u22121$. 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.

### Temperature dependence of calculated thermal conductivity

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}.^{39} Therefore, calculated *κ*_{MD} values for *T* > *θ*_{D} ∼ 263 K^{46} are realistic, but extrapolations of this method to low temperatures (∼100 K)^{24} 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.

### Non-classical specific heat calculations

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\omega ,T=exp\u0127\omega kBT\u22121\u22121$. 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 by isotopes

Phonon scattering in experimental MoS_{2} 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

^{49}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_{2} 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.^{15} 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_{2} 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_{2} systems.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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.