Thermal Conductivity of Monolayer MoSe2 and MoS2 - The Journal of

Subscriber access provided by La Trobe University Library

Article 2

2

Thermal Conductivity of Monolayer MoSe and MoS Yang Hong, Jingchao Zhang, and Xiao Cheng Zeng

J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.6b07262 • Publication Date (Web): 25 Oct 2016 Downloaded from http://pubs.acs.org on October 29, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Thermal Conductivity of Monolayer MoSe2 and MoS2 Yang Honga, Jingchao Zhangb*, Xiao Cheng Zenga* a

Department of Chemistry, University of Nebraska-Lincoln, Lincoln, NE 68588, USA

b

Holland Computing Center, University of Nebraska-Lincoln, Lincoln, NE, 68588, USA

ABSTRACT: Two dimensional (2D) MoSe2 and MoS2 monolayers, two prototype transition metal dichalcogenides (TMDCs) materials, have attracted growing interests as promising 2D semiconductors. In this work, thermal conductivity () of the monolayer MoSe2 is computed using large-scale classical non-equilibrium molecular dynamics (NEMD) simulations for the first time. The predicted  of monolayer MoSe2 with infinite length (or MoSe2 2D sheets) are 43.88 ± 1.33 and 41.63 ± 0.66 W/mK in armchair and zigzag direction, respectively. These simulation results are further confirmed by independent simulations using the Green-Kubo method (GKM) which yields computed  of 44.38 ± 2.08 and 44.63 ± 2.50 W/mK, respectively. For 2D MoS2 sheet, the computed  based on the NEMD method are 101.43 ± 1.13 and 110.30 ± 2.07 W/mK, respectively, in armchair and zigzag direction, whereas those based on the GKM are 102.32 ± 6.05 and 108.74 ± 6.68 W/mK, respectively. The predicted  values of MoS2 monolayer are twice larger than those of MoSe2 monolayer. Both types of 2D monolayers exhibit isotropic properties in thermal conduction. Effects of system dimensions, heat flux, and temperature on  are investigated comprehensively. The predicted  value increases monotonically with the system length but decreases with temperature.

* Corresponding authors: [email protected] (402-472-9894); [email protected] (402-472-6400) 1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 33

INTRODUCTION Graphene is one of the most widely studied 2D materials in recent years. Despite of its superb thermal properties1-2 and high intrinsic mobility3, its applications in nanoelectronic devices such as field effect transistors (FETs) are limited due to the zero band gap property.4-6 Although graphene nanotube7-8 and electrically gating bilayer graphene9 can enlarge its band gap up to 0.4 eV, its high mobility would lose concomitantly. To compensate this deficiency, researchers explored another group of 2D materials, categorized as the transition metal dichalcogenides (TMDCs)10-11, such as MoSe2 and MoS2 monolayers. It has been shown that mechanically exfoliated MoSe2 and MoS2 monolayers from bulk counterparts12-13 can exhibit a direct band gap of 1.55 eV12 and 1.8 eV14, respectively, rendering both monolayers highly promising for applications in FET and other optical devices.15-16 Indeed, MoS2-based FETs were reported to possess relatively low electron mobility of 200 cm2/Vs but extremely high switching on/off ratio (~108), compared to graphene.17 The bandgap of MoSe2 matches the optimum bandgap of singlejunction solar cells and photo eletrochemical devices, enabling its applications for energy conversion from sun light.12,

18

Moreover, Huang et al.19 discovered that the enhanced

photoluminescence of in-plane heterojunctions between monolayer MoSe2 and WSe2 allows their usages as in-plane transistors and diodes. The high thermally stable and desired direct-gap semiconducting properties along with many other unique physical and optical properties facilitate growing application interests of monolayer MoSe2 and MoS2 in sensors, saturable absorber of Q-switched Erbium-doped fiber laser, photocatalyst, and electroluminescence.20-25 As MoSe2 and MoS2 monolayers are promising candidates for the next generation of nanoelectronic materials, studies of their thermal conductivities () become timely and crucial. Some previous studies have reported thermal conductivity of monolayer MoSe2 and MoS2 at 2


Page 3 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


room temperature, both theoretically and experimentally. Previous molecular dynamics (MD) simulation studies often underestimated the  values of monolayer MoS2 within the range of 1.35~5.8 W/mK, due largely to the force field employed.26-27 The first-principles calculations and realistic experiments give relative higher  values. In a temperature-dependent Raman spectroscopy experiment, Yan et al.28 measured the  value of monolayer MoS2 at 34.5 ± 4 W/m·K. Using opto-thermal Raman technique, Zhang et al.29 found the  value of monolayer MoS2 and MoSe2 to be 84 ± 17 and 59 ± 18 W/m·K, respectively. Cai et al.30 obtained a  value of 23.2 W/mK for monolayer MoS2, based on density functional perturbation theory (DFPT) computation. By combining Peierls-Boltzmann transport equation (PBTE) and first-principles calculation, Li et al.31 predicted a  value of 83 W/mK for monolayer MoS2. Using a similar method, Gu et al.32 estimated the  value of monolayer MoS2 and MoSe2 to be 103 W/mK and 54 W/mK, respectively. Despites of numerous studies of thermal conductivity of monolayer MoS2, to our knowledge, the thermal properties of MoSe2 monolayer have not yet been studied comprehensively with MD approaches. In this work, thermal conductivities of monolayer MoSe2 sheet in armchair and zigzag directions are computed for the first time using non-equilibrium molecular dynamics (NEMD) approach33. Meanwhile, thermal conductivities of monolayer MoS2 sheet are also investigated for the purpose of comparison with previous studies and with those of MoSe2. Size effects on the thermal conductivity of monolayer MoSe2 and MoS2 are examined by analyzing computed thermal conductivities versus the system length and width. Temperature and energy dependences of the thermal conductivity for monolayer MoSe2 and MoS2 are also studied. Lastly, the thermal conductivities of monolayer MoSe2 and MoS2 sheets in both armchair and zigzag directions are

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 33

computed and confirmed by using the equilibrium molecular dynamics (EMD) Green-Kubo method (GKM).34 COMPUTATIONAL METHODS AND MODEL SYSTEMS All MD simulations are performed with the large-scale atomic/molecular massively parallel simulator (LAMMPS).35 The Stillinger-Weber (SW) potential developed by Kandemir et al.36 is used to describe the interatomic interactions within monolayer MoSe2 and MoS2. The formula of SW potential is divided into a 2-body bond stretching term and a 3-body bond bending term. To assure an accurate description of the 3-body bond bending term, the top layer Se/S atoms are treated as a different type of atoms from the bottom layer Se/S atoms. Note the SW potential was developed by fitting lattice constants, bond lengths, elastic constants, and vibrational properties of monolayer MoSe2 and MoS2 via the particle swarm optimization method.37 Thermal properties generated by this potential are in good agreement with results from pervious first-principle computations30-32 and experimental measurements28,

38-39

. Thus, it is credible to use the SW

potential to evaluate thermal conductivities of monolayer MoSe2 and MoS2. Atomic configuration schematics of monolayer MoSe2 and MoS2 are shown in Figure 1. Their thermal conductivities along the armchair (x) and zigzag (y) directions are computed by using the NEMD method, and the GKM (for the purpose of verification). In the NEMD simulations, periodic boundary condition is applied to the width direction (x) while free boundary conditions in the heat flux direction (y) and out-of-plane direction (z). Atoms in the outmost atomic chain at both ends are fixed. The next four chains of atoms at both ends are grouped together as heat-bath and heat-sink regions (see Figure 2a for a schematic plot of the NEMD simulation system with the heat flux in the armchair direction). The thicknesses of 4


Page 5 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


monolayers are half of their vertical lattice constant c. For monolayer MoSe2, the thickness is 6.469 Å40 and for MoS2 the thickness is 6.1475 Å.41 The trajectory of the system at each time step of 0.5 fs is recorded in the simulation. At the beginning of MD simulation, the monolayer is placed in the canonical ensemble (NVT) for 500 ps at 300 K and then is in the microcanonical ensemble (NVE) for another 500 ps. After thermal equilibrium simulation, the system remains in the NVE ensemble while the NEMD method is applied to compute the thermal conductivity for additional 2.5 ns, during which the same amount of kinetic energy is concurrently and constantly added to the heat bath region while subtracted from the heat sink region until the whole system reaches the steady state. Once the heat flux becomes steady, a stable temperature gradient along the heat flux direction (y) is generated.33 The slope of the temperature gradient is used to compute the thermal conductivity of monolayers via the Fourier’s law of heat conduction:

J y   y  A 

T , y

(1)

where Jy is heat flux in y direction; ∂T/∂y stands for the temperature change along the heat flux direction, i.e., the slope of temperature gradient, and A is cross section area. Addition and subtraction of kinetic energy at heat bath and heat sink regions are achieved by multiplying or dividing the velocity of atoms within these regions by a same factor. As the total energy of the system remains constant and the linear momentum at each time step is conserved, no external thermostating is needed. Snapshots of atomic velocities within the whole system are recorded before and after the heating/cooling process. Once system reaches steady state, the velocities of atoms within the system should follow the Maxwell–Boltzmann distribution  m  PM  4 v 2    2 k BT 

3/2

e



mv 2 2 k BT

,

(2)

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 33

where PM is the probability of an atom with velocity ν; m represents atomic mass, and kB is the Boltzmann constant. The simulated velocity distribution and the Maxwell–Boltzmann distribution are compared to confirm that the system has reached the steady state. During the last 500 ps of NEMD simulation, the system is uniformly divided into many slabs along the heat flux direction (y) and each slab contains at least one chain of atoms. The kinetic energies of these slabs are recorded and used to calculate their temperatures via the energy equipartition equation: N 1 3  E    m i 2  Nk BT , 2 1 2

(3)

where vi is the velocity of atom i and N is number of atoms within a slab. The calculated temperatures of each slab are averaged over 500 ps after the system reaches the steady state. The temperature gradient profile is produced by combining the average temperatures and positions of each slab. Then, the thermal conductivity can be evaluated based on the slope of temperature gradient curve as discussed above. In the GKM, the thermal conductivity is expressed in term of time-integration of heat current autocorrelation function (HCAF),34

 xy 

1 k BT 2V



J

x

(t )  J y (0)  dt ,

(5)

0

where V is the system volume; kB is Boltzmann constant; T is the temperature of system and Jx, Jy are the heat current along x and y directions. The angular brackets denote the average over time. The upper limit of HCAF integration time in the Eq. (5) is infinite, while the integration time is finite in MD simulations. Thus, as long as the upper limit of integral time is chosen longer than the time beyond which the current autocorrelations converge to zero, the use of Eq. 6


Page 7 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(5) would be meaningful. Once the system reaches equilibrium, the heat current of every 10 time steps is saved for 1.6×106 steps, giving the HCAF integration time of 800 ps. The MD simulation time lasts for 8 ns, ten times longer than the HCAF integration time in order to obtain an accurate statistical average. To reduce the influence of HCAF noise on thermal conductivity predictions, for each monolayer, its thermal conductivity is computed 10 times with different initial velocity seed. RESULTS AND DISCUSSION Heat flux of Jin = 9.74 × 107 W is added to the heat bath at each time step and the same amount Jout is subtracted from the heat sink simultaneously for 2.5 ns. Figure 2b is an example of temperature gradient along the heat flux direction for a 43.15 nm long armchair MoSe2 nanoribbon. The black dots represent MD temperatures and the red line stands for the linear fitting result based on Eq. (1). During the extremely fast energy exchange process in the NEMD simulation, kinetic energy and potential energy near the heating and cooling regions are nonequilibrium and phonon boundary scatterings are very rapid in these areas. Hence, the temperature changes in these regions are non-linear and should be removed from the linear fittings.42-45 To confirm that the system has reached the steady state before and after the heating/cooling process, atomic velocities of the whole system are extracted and compared with Maxwellian velocity distributions according to Eq. (2). As shown in Figure 3, the velocity distributions of 163.46 nm long armchair MoSe2 nanoribbon are mapped across the range from 0 to 1000 m/s. The solid lines stand for the theoretical Maxwellian velocity distributions and the dots represent MD simulation results. It is observed that both before (Figure 3a) and after (Figure 3b) the heating/cooling process, MD simulation results match Maxwellian velocity distributions

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 33

well, demonstrating the system has reached the steady state before heating/cooling and data collection. In the NEMD method, the length of system along the heat flux direction is finite. When the length of system is shorter than the phonon mean free path (MFP) of the infinite system, the thermal conductivity varies with the system size.46 To analyze the size effect on the thermal conductivity of monolayer MoSe2 system in armchair and zigzag direction, two sets of models are considered with different length, while the width is the same at ~10 nm. Specifically, for the armchair direction  computations, the length varies from 10.64, 21.48, 43.15, 86.48, 163.46, 326.53, 433.16 to 519.83 nm. For the zigzag direction  computation, the length varies from 10.04, 20.25, 40.66, 81.48, 163.12, 326.40, 408.04 to 530.06 nm. To obtain more accurate thermal conductivity, five independent MD simulations are performed and the results are averaged. The error bar is calculated using the standard deviation. As summarized in Figure 4a, the computed  value shows a monotonically increasing trend with the length in both directions. For the armchair direction, with increasing system length, the  value increases from 1.83 to 24.12 W/mK. For zigzag direction with similar length changes, the  value rises from 2.04 to 24.89 W/mK, indicating the isotropy of thermal conduction for monolayer MoSe2. To minimize the size effect and predict the  value of monolayer MoSe2 sheet, the following linear function is used to fit the computed  values of finite length monolayer MoSe2 nanoribbons,34 1





1

l (  1) ,  L

(4)

where l is effective phonon mean free path and  is thermal conductivity of infinite length monolayer. The fitting results of 1/ with 1/L are given in Figure 5a. With 5 data points, the 8


Page 9 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


fitted  of monolayer MoSe2 sheet in armchair and zigzag direction are 42.55 and 40.97 W/mK, respectively. With 4 data points, the fitted  of monolayer MoSe2 sheet in armchair and zigzag direction are 45.21 and 42.28 W/mK. The difference between 4-point fitting and 5-point fitting are only 5.88% and 3.10%, individually. Thus, we can conclude that, at room temperature, thermal conductivities of monolayer MoSe2 sheet are 43.88 ± 1.33 and 41.63 ± 0.66 W/mK in armchair and zigzag direction, respectively. These values are quite close to the  value of 54 W/mK32 from the first-principles computation, and the experiment value of 59 ± 18 W/m·K29 at room temperature. The size effect of monolayer MoS2 is also studied. As shown in Figure 4a, the  values with length of 10.22, 20.62, 41.43, 83.04, 160.23, 299.85, 415.92 and 499.50 nm in armchair direction and with length of 9.96, 20.07, 40.30, 80.76, 161.68, 303.30, 404.45 and 485.37 nm in zigzag direction are obtained. For the armchair direction, with longer system length, the  value increases from 5.46 to 65.34 W/mK. For zigzag direction with similar length changes, the  value rises from 5.46 to 61.04 W/mK, again indicating the isotropy of thermal conduction for monolayer MoS2. With 5 data points, the fitted  values of monolayer MoS2 sheet in armchair and zigzag direction are 100.30 and 108.23 W/mK, respectively. With 4 data points, the fitted  values of infinitely long monolayer MoSe2 are 102.56 and 112.36 W/mK, respectively. The differences between 4-point fitting and 5-point fitting are 2.20% and 3.68%, respectively. Therefore, we can predict that the room temperature thermal conductivities of monolayer MoS2 sheet are 101.43 ± 1.13 and 110.30 ± 2.07 W/mK in armchair and zigzag direction, respectively. Both values are consistent with the first-principles predicted thermal conductivity of 103

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 33

W/mK32 and the experiment value of 84 ± 17 W/m·K29 for monolayer MoS2 at room temperature. Since periodic boundary condition is applied to the width direction of system, the system width should be considered as infinitely long and not affect the result of thermal conductivity. To confirm this prediction, we examine systems with 5.10, 10.04, 20.25, 30.45 and 40.66 nm width of MoSe2 sheet but the same length of 43.15 nm. As shown in Figure 4b, the thermal conductivity is converged with the system width. To lower the computational cost, we consider all systems of ~10 nm width with periodic boundary condition for thermal conductivity computation. The working temperature of monolayer MoSe2 and MoS2 can be different from the room temperature. To examine the temperature effect on the thermal conductivity, two ~160 nm long monolayer MoSe2 and MoS2 systems are subjected to 100, 200, 300, 400 and 500 K heat bathes separately. The computed  values are summarized in Figure 6a. For monolayer MoSe2,  values in both armchair and zigzag directions exhibit a monotonically decreasing trend with the increasing temperature. At high temperature, more phonons with higher frequency would be excited, and involved in the thermal transport. As a result, the entire phonon population is raised. Therefore, the Umklapp phonon scatterings become more active and can prominently limit the thermal conductivity of monolayer MoSe2 at high temperature.47 When the temperature of monolayer MoSe2 increases from 100 to 500 K, the  value is lowered by 43.66% and 44.37% in armchair and zigzag direction, separately. The  data are fitted with an inverse relationship with temperature (~1/T) as shown in Figure 6a. The fitting curves can well match the simulation results, indicating the Umklapp scattering is more dominant in this temperature range. The

10


Page 11 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


decreasing speed and trend of thermal conductivity in armchair direction are similar to those in zigzag direction, showing the isotropy of temperature effects for thermal conduction of monolayer MoSe2. Likewise, the  value of monolayer MoS2 in armchair and zigzag directions also decrease with temperature in a nearly similar fashion, suggesting that the temperature effect on monolayer MoS2 thermal conductivity is isotropic as well. With the temperature of monolayer MoS2 increasing from 100 to 500 K,  value is reduced by 48.70% and 46.86% in armchair and zigzag direction, respectively. As shown in Figure 6a, the 1/T fitting curves for MoS2 also perfectly match with the simulation results, indicating the Umklapp scattering is the main cause in this temperature range as well. The influence of heat flux on the predicted thermal conductivity is investigated. By altering the added/subtracted heat flux for the 43.15 nm long armchair MoSe2 system from 2.99×109 to 29.90×109 W/m2, the temperature gradient increases from 13 to 157 K. As presented in Figure 6b, the thermal conductivity is converged with heat flux. The change of added/subtracted heat flux amount to/from the system is associated with the change of temperature gradient; however, it has little influence on the thermal conductivity.

Simulations based on the GKM. Lastly, the thermal conductivities of monolayer MoSe2 and MoS2 sheets at room temperature are also evaluated independently by using the GKM to confirm the NEMD results. The periodic boundary conditions are applied in both length and width directions, while z direction is still the free boundary condition. Therefore, the size effect of GKM is much smaller 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 33

than that of NEMD method.34 For the size effect analysis,  values of monolayer MoSe2 and MoS2 systems containing 5×5 to 60×60 unit cells are computed. As shown in Figure 7,  values of monolayer MoSe2 and MoS2 are both converged at 50×50 unit cells. Thus, the systems containing 50×50 unit cells, which are 34.02 nm and 32.76 nm long for monolayer MoSe2 and MoS2, respectively, are used for GKM thermal conductivity simulation. As shown in Figure 8a, the averaged HCAFs of monolayer MoSe2 decrease to 0 beyond ~150 ps and the overall  values of ten trails are converged after 250 ps in both armchair and zigzag directions (see Figures 9a and 9b). The  value of monolayer MoSe2 sheet is then calculated by averaging the values of  from 250 to 800 ps, which is 44.38 ± 2.08 and 44.63 ± 2.50 W/mK in armchair and zigzag direction, respectively. Both values are consistent with the NEMD results of 43.88 ± 1.33 and 41.63 ± 0.66 W/mK, and only differ by 1.13% and 6.72%, respectively. Note that both methods indicate isotropy of thermal conduction for the monolayer MoSe2. For monolayer MoS2, the averaged HCAFs decline to 0 after 200 ps (see Figure 8b) and the overall  values of ten trails are converged after 500 ps for both armchair and zigzag directions (see Figures 9c and 9d). The  value of monolayer MoS2 sheet is evaluated by taking average of overall  values from 500 to 800 ps, which is 102.32 ± 6.05 and 108.74 ± 6.68 W/mK in armchair and zigzag direction, respectively. Again, both values are very close to the NEMD results of 101.43 ± 1.13 and 110.30 ± 2.07 W/mK, and only differ by 0.87% and 1.41%, respectively. Although both methods show ~5% variance of  values in armchair and zigzag directions, in practical applications, this small difference can be ignored. So we can conclude that the thermal conductions of both monolayer MoSe2 and MoS2 are isotropic.

12


Page 13 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Comparison of Thermal Conductivities with Other 2D Monolayer Materials. In Table 1, we summarize the room temperature  values of several popular 2D monolayer materials that have been computed based on classical MD simulations, compared with those computed in this work. The corresponding experimental values and results from previous first-principle computations are also shown in Table 1. Only monolayer phosphorene shows stronger anisotropic  values while other materials are all isotropic. Suspended graphene sheet has the highest  among the 2D monolayers listed in Table 1, more than one order of magnitude higher than those of monolayer MoSe2 and MoS2. By nano-manufacturing hybrid structures of graphene and monolayer MoSe2/MoS2 or employing graphene as an interfacial material between monolayer MoSe2/MoS2 and substrates, one may be able to combine high thermal conductivity of graphene and direct-gap semiconducting property of monolayer MoSe2 and MoS2 for practical applications. CONCLUSION Using large-scale and independent NEMD and GKM simulations, thermal conductivities of monolayer MoSe2 are computed and compared with monolayer MoS2. Thermal conductivities of monolayer MoSe2 sheet predicted by NEMD are 43.88 ± 1.33 and 41.63 ± 0.66 W/mK in armchair and zigzag direction, respectively, very close to the GKM results of 44.38 ± 2.08 and 44.63 ± 2.50 W/mK. The variances between two simulation methods are 1.13% and 6.72%, correspondingly. Moreover, the computed thermal conductivities of armchair and zigzag monolayer MoS2 sheet based on NEMD are 101.43 ± 1.13 and 110.30 ± 2.07 W/mK, respectively, while those based on GKM are 102.32 ± 6.05 and 108.74 ± 6.68 W/mK, respectively. The differences between the two methods are only 0.87% and 1.41%, respectively. The thermal conductivity of monolayer MoS2 sheet is twice higher than that of monolayer MoSe2

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 33

sheet. Both 2D materials show isotropic properties of thermal conduction. For the NEMD simulations, effects of system dimensions, heat flux, and temperature on thermal conductivity are investigated comprehensively. The predicted thermal conductivity increases monotonically with the system length but decreases with temperature for both 2D materials.

ACKNOWLEDGEMENT Support of this work from Nebraska Center for Energy Sciences Research and Holland Computing Center of the University of Nebraska-Lincoln is greatly appreciated.

REFERENCES 1.

Jacoby, M. Graphene's Thermal Conductivity. Chem Eng News 2010, 88 (15), 5-5.

2.

Balandin, A. A.; Ghosh, S.; Bao, W. Z.; Calizo, I.; Teweldebrhan, D.; Miao, F.; Lau, C.

N. Superior thermal conductivity of single-layer graphene. Nano Letters 2008, 8 (3), 902-907. 3.

Bolotin, K. I.; Sikes, K.; Jiang, Z.; Klima, M.; Fudenberg, G.; Hone, J.; Kim, P.; Stormer,

H. Ultrahigh electron mobility in suspended graphene. Solid State Commun 2008, 146 (9), 351355. 4.

Shimizu, T.; Haruyama, J.; Marcano, D. C.; Kosinkin, D. V.; Tour, J. M.; Hirose, K.;

Suenaga, K. Large intrinsic energy bandgaps in annealed nanotube-derived graphene nanoribbons. Nat Nanotechnol 2011, 6 (1), 45-50. 5.

Zhang, S. J.; Lin, S. S.; Li, X. Q.; Liu, X. Y.; Wu, H. A.; Xu, W. L.; Wang, P.; Wu, Z. Q.;

Zhong, H. K.; Xu, Z. J. Opening the band gap of graphene through silicon doping for the

14


Page 15 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


improved performance of graphene/GaAs heterojunction solar cells. Nanoscale 2016, 8 (1), 226232. 6.

Han, S. J.; Jenkins, K. A.; Garcia, A. V.; Franklin, A. D.; Bol, A. A.; Haensch, W. High-

Frequency Graphene Voltage Amplifier. Nano Letters 2011, 11 (9), 3690-3693. 7.

Ming-Wei, L.; Cheng, L.; Yiyang, Z.; Hyeun Joong, Y.; Mark Ming-Cheng, C.; Luis, A.

A.; Nicholas, K.; Noppi, W.; Zhixian, Z. Room-temperature high on/off ratio in suspended graphene nanoribbon field-effect transistors. Nanotechnology 2011, 22 (26), 265201. 8.

Han, M. Y.; Özyilmaz, B.; Zhang, Y.; Kim, P. Energy Band-Gap Engineering of

Graphene Nanoribbons. Physical Review Letters 2007, 98 (20), 206805. 9.

Zhang, Y.; Tang, T.-T.; Girit, C.; Hao, Z.; Martin, M. C.; Zettl, A.; Crommie, M. F.;

Shen, Y. R.; Wang, F. Direct observation of a widely tunable bandgap in bilayer graphene. Nature 2009, 459 (7248), 820-823. 10.

Wilson, J. A.; Yoffe, A. D. The transition metal dichalcogenides discussion and

interpretation of the observed optical, electrical and structural properties. Advances in Physics 1969, 18 (73), 193-335. 11.

Mattheiss, L. F. Band Structures of Transition-Metal-Dichalcogenide Layer Compounds.

Physical Review B 1973, 8 (8), 3719-3740. 12.

Tongay, S.; Zhou, J.; Ataca, C.; Lo, K.; Matthews, T. S.; Li, J.; Grossman, J. C.; Wu, J.

Thermally Driven Crossover from Indirect toward Direct Bandgap in 2D Semiconductors: MoSe2 versus MoS2. Nano Letters 2012, 12 (11), 5576-5580. 13.

Novoselov, K. S.; Jiang, D.; Schedin, F.; Booth, T. J.; Khotkevich, V. V.; Morozov, S. V.;

Geim, A. K. Two-dimensional atomic crystals. Proceedings of the National Academy of Sciences of the United States of America 2005, 102 (30), 10451-10453.

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

14.

Page 16 of 33

Mak, K. F.; Lee, C.; Hone, J.; Shan, J.; Heinz, T. F. Atomically Thin MoS2: A New

Direct-Gap Semiconductor. Physical Review Letters 2010, 105 (13), 136805. 15.

Wang, Q. H.; Kalantar-Zadeh, K.; Kis, A.; Coleman, J. N.; Strano, M. S. Electronics and

optoelectronics of two-dimensional transition metal dichalcogenides. Nat Nano 2012, 7 (11), 699-712. 16.

Chhowalla, M.; Shin, H. S.; Eda, G.; Li, L.-J.; Loh, K. P.; Zhang, H. The chemistry of

two-dimensional layered transition metal dichalcogenide nanosheets. Nat Chem 2013, 5 (4), 263275. 17.

RadisavljevicB; RadenovicA; BrivioJ; GiacomettiV; KisA. Single-layer MoS2 transistors.

Nat Nano 2011, 6 (3), 147-150. 18.

Mirhosseini, H.; Kiss, J.; Roma, G.; Felser, C. Reducing the Schottky barrier height at the

MoSe2/Mo(110) interface in thin-film solar cells: Insights from first-principles calculations. Thin Solid Films 2016, 606, 143-147. 19.

Huang, C.; Wu, S.; Sanchez, A. M.; Peters, J. J. P.; Beanland, R.; Ross, J. S.; Rivera, P.;

Yao, W.; Cobden, D. H.; Xu, X. Lateral heterojunctions within monolayer MoSe2–WSe2 semiconductors. Nat Mater 2014, 13 (12), 1096-1101. 20.

Chu, H. P.; Lei, W. Y.; Liu, X. J.; Qu, J. H.; Li, J. L.; Zhu, G.; Niu, L. Y.; Pan, L. K.

MoSe2 visible-light photocatalyst for organic pollutant degradation and Cr(VI) reduction. J Mater Sci-Mater El 2016, 27 (5), 5483-5489. 21.

Li, H.; Yin, Z.; He, Q.; Li, H.; Huang, X.; Lu, G.; Fam, D. W. H.; Tok, A. I. Y.; Zhang,

Q.; Zhang, H. Fabrication of Single- and Multilayer MoS2 Film-Based Field-Effect Transistors for Sensing NO at Room Temperature. Small 2012, 8 (1), 63-67.

16


Page 17 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


22.

Late, D. J.; Doneux, T.; Bougouma, M. Single-layer MoSe2 based NH3 gas sensor.

Applied Physics Letters 2014, 105 (23). 23.

Tsai, M.-L.; Su, S.-H.; Chang, J.-K.; Tsai, D.-S.; Chen, C.-H.; Wu, C.-I.; Li, L.-J.; Chen,

L.-J.; He, J.-H. Monolayer MoS2 Heterojunction Solar Cells. Acs Nano 2014, 8 (8), 8317-8322. 24.

Sundaram, R. S.; Engel, M.; Lombardo, A.; Krupke, R.; Ferrari, A. C.; Avouris, P.;

Steiner, M. Electroluminescence in Single Layer MoS2. Nano Letters 2013, 13 (4), 1416-1421. 25.

Ahmad, H.; Suthaskumar, M.; Tiu, Z. C.; Zarei, A.; Harun, S. W. Q-switched Erbium-

doped fiber laser using MoSe2 as saturable absorber. Opt Laser Technol 2016, 79, 20-23. 26.

Liu, X.; Zhang, G.; Pei, Q.-X.; Zhang, Y.-W. Phonon thermal conductivity of monolayer

MoS2 sheet and nanoribbons. Applied Physics Letters 2013, 103 (13), 133113. 27.

Jiang, J.-W.; Park, H. S.; Rabczuk, T. Molecular dynamics simulations of single-layer

molybdenum disulphide (MoS2): Stillinger-Weber parametrization, mechanical properties, and thermal conductivity. J Appl Phys 2013, 114 (6), 064307. 28.

Yan, R.; Simpson, J. R.; Bertolazzi, S.; Brivio, J.; Watson, M.; Wu, X.; Kis, A.; Luo, T.;

Hight Walker, A. R.; Xing, H. G. Thermal Conductivity of Monolayer Molybdenum Disulfide Obtained from Temperature-Dependent Raman Spectroscopy. Acs Nano 2014, 8 (1), 986-993. 29.

Zhang, X.; Sun, D.; Li, Y.; Lee, G.-H.; Cui, X.; Chenet, D.; You, Y.; Heinz, T. F.; Hone,

J. C. Measurement of Lateral and Interfacial Thermal Conductivity of Single- and Bilayer MoS2 and MoSe2 Using Refined Optothermal Raman Technique. Acs Appl Mater Inter 2015, 7 (46), 25923-25929. 30.

Cai, Y.; Lan, J.; Zhang, G.; Zhang, Y.-W. Lattice vibrational modes and phonon thermal

conductivity of monolayer MoS2. Physical Review B 2014, 89 (3), 035438.

17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

31.

Page 18 of 33

Li, W.; Carrete, J.; Mingo, N. Thermal conductivity and phonon linewidths of monolayer

MoS2 from first principles. Applied Physics Letters 2013, 103 (25), 253103. 32.

Gu, X.; Yang, R. Phonon transport in single-layer transition metal dichalcogenides: A

first-principles study. Applied Physics Letters 2014, 105 (13), 131903. 33.

Müller-Plathe, F. A simple nonequilibrium molecular dynamics method for calculating

the thermal conductivity. The Journal of Chemical Physics 1997, 106 (14), 6082-6085. 34.

Schelling, P. K.; Phillpot, S. R.; Keblinski, P. Comparison of atomic-level simulation

methods for computing thermal conductivity. Physical Review B 2002, 65 (14), 144306. 35.

Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. Journal of

Computational Physics 1995, 117 (1), 1-19. 36.

Ali, K.; Haluk, Y.; Alper, K.; Tahir, Ç.; Cem, S. Thermal transport properties of MoS2

and MoSe2 monolayers. Nanotechnology 2016, 27 (5), 055703. 37.

Kennedy, J.; Eberhart, R. In Particle swarm optimization, Neural Networks, 1995.

Proceedings., IEEE International Conference on, Nov/Dec 1995; 1995; pp 1942-1948 vol.4. 38.

Jo, I.; Pettes, M. T.; Ou, E.; Wu, W.; Shi, L. Basal-plane thermal conductivity of few-

layer molybdenum disulfide. Applied Physics Letters 2014, 104 (20), 201902. 39.

Sahoo, S.; Gaur, A. P. S.; Ahmadi, M.; Guinel, M. J. F.; Katiyar, R. S. Temperature-

Dependent Raman Studies and Thermal Conductivity of Few-Layer MoS2. The Journal of Physical Chemistry C 2013, 117 (17), 9042-9047. 40.

Böker, T.; Severin, R.; Müller, A.; Janowitz, C.; Manzke, R.; Voß, D.; Krüger, P.; Mazur,

A.; Pollmann, J. Band structure of MoS2, MoSe2, and α-MoTe2: Angle-resolved photoelectron spectroscopy and ab initio calculations. Physical Review B 2001, 64 (23), 235305.

18


Page 19 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


41.

Swanson H E, G. N. T. a. U. G. M. Standard X-ray Diffraction Powder Patterns NBS

Circular 1955, 5, 1953-1960. 42.

Zhong, H.; Lukes, J. R. Interfacial thermal resistance between carbon nanotubes:

Molecular dynamics simulations and analytical thermal modeling. Physical Review B 2006, 74 (12), 125403. 43.

Landry, E. S.; McGaughey, A. J. H. Thermal boundary resistance predictions from

molecular dynamics simulations and theoretical calculations. Physical Review B 2009, 80 (16), 165304. 44.

Zhang, J.; Wang, Y.; Wang, X. Rough contact is not always bad for interfacial energy

coupling. Nanoscale 2013, 5 (23), 11598-11603. 45.

Hong, Y.; Li, L.; Zeng, X. C.; Zhang, J. Tuning thermal contact conductance at graphene-

copper interface via surface nanoengineering. Nanoscale 2015, 7 (14), 6286-6294. 46.

Sellan, D. P.; Landry, E. S.; Turney, J. E.; McGaughey, A. J. H.; Amon, C. H. Size

effects in molecular dynamics thermal conductivity predictions. Physical Review B 2010, 81 (21), 214305. 47.

Liu, B.; Baimova, J. A.; Reddy, C. D.; Law, A. W.-K.; Dmitriev, S. V.; Wu, H.; Zhou, K.

Interfacial Thermal Conductance of a Silicene/Graphene Bilayer Heterostructure and the Effect of Hydrogenation. Acs Appl Mater Inter 2014, 6 (20), 18180-18188. 48.

Hong, Y.; Zhang, J.; Huang, X.; Zeng, X. C. Thermal conductivity of a two-dimensional

phosphorene sheet: a comparative study with graphene. Nanoscale 2015, 7 (44), 18716-18724. 49.

Ng, T. Y.; Yeo, J. J.; Liu, Z. S. A molecular dynamics study of the thermal conductivity

of graphene nanoribbons containing dispersed Stone-Thrower-Wales defects. Carbon 2012, 50 (13), 4887-4893.

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

50.

Page 20 of 33

Zhang, J.; Huang, X.; Yue, Y.; Wang, J.; Wang, X. Dynamic response of graphene to

thermal impulse. Physical Review B 2011, 84 (23), 235416. 51.

Li, C. J.; Li, G.; Zhao, H. J. A Molecular Dynamics Study about Two Way Tuning of

Thermal Conductivity in Graphene: Strain and Doping. 2013 13th Ieee Conference on Nanotechnology (Ieee-Nano) 2013, 913-916. 52.

Hu, J. N.; Ruan, X. L.; Chen, Y. P. Thermal Conductivity and Thermal Rectification in

Graphene Nanoribbons: A Molecular Dynamics Study. Nano Letters 2009, 9 (7), 2730-2735. 53.

Evans, W. J.; Hu, L.; Keblinski, P. Thermal conductivity of graphene ribbons from

equilibrium molecular dynamics: Effect of ribbon width, edge roughness, and hydrogen termination. Applied Physics Letters 2010, 96 (20). 54.

Lee, J.-U.; Yoon, D.; Kim, H.; Lee, S. W.; Cheong, H. Thermal conductivity of

suspended pristine graphene measured by Raman spectroscopy. Physical Review B 2011, 83 (8), 081419. 55.

Kong, B. D.; Paul, S.; Nardelli, M. B.; Kim, K. W. First-principles analysis of lattice

thermal conductivity in monolayer and bilayer graphene. Physical Review B 2009, 80 (3), 033406. 56.

Tabarraei, A. Thermal conductivity of monolayer hexagonal boron nitride nanoribbons.

Comp Mater Sci 2015, 108, 66-71. 57.

Mortazavi, B.; Rémond, Y. Investigation of tensile response and thermal conductivity of

boron-nitride nanosheets using molecular dynamics simulations. Physica E: Low-dimensional Systems and Nanostructures 2012, 44 (9), 1846-1852. 58.

Lindsay, L.; Broido, D. A. Enhanced thermal conductivity and isotope effect in single-

layer hexagonal boron nitride. Physical Review B 2011, 84 (15), 155421.

20


Page 21 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


59.

Wang, C. R.; Guo, J.; Dong, L.; Aiyiti, A.; Xu, X. F.; Li, B. W. Superior thermal

conductivity in suspended bilayer hexagonal boron nitride. Sci Rep-Uk 2016, 6, 25334. 60.

Jo, I.; Pettes, M. T.; Kim, J.; Watanabe, K.; Taniguchi, T.; Yao, Z.; Shi, L. Thermal

Conductivity and Phonon Transport in Suspended Few-Layer Hexagonal Boron Nitride. Nano Letters 2013, 13 (2), 550-554. 61.

Pei, Q.-X.; Zhang, Y.-W.; Sha, Z.-D.; Shenoy, V. B. Tuning the thermal conductivity of

silicene with tensile strain and isotopic doping: A molecular dynamics study. J Appl Phys 2013, 114 (3), 033526. 62.

Zhang, X.; Xie, H.; Hu, M.; Bao, H.; Yue, S.; Qin, G.; Su, G. Thermal conductivity of

silicene calculated using an optimized Stillinger-Weber potential. Physical Review B 2014, 89 (5), 054310. 63.

Guo, Y.; Zhou, S.; Bai, Y. Z.; Zhao, J. J. Tunable Thermal Conductivity of Silicene by

Germanium Doping. J Supercond Nov Magn 2016, 29 (3), 717-720. 64.

Xie, H.; Hu, M.; Bao, H. Thermal conductivity of silicene from first-principles. Applied

Physics Letters 2014, 104 (13), 131906. 65.

Cherukara, M. J.; Narayanan, B.; Kinaci, A.; Sasikumar, K.; Gray, S. K.; Chan, M. K. Y.;

Sankaranarayanan, S. K. R. S. Ab Initio-Based Bond Order Potential to Investigate Low Thermal Conductivity of Stanene Nanostructures. The Journal of Physical Chemistry Letters 2016, 37523759. 66.

Kuang, Y. D.; Lindsay, L.; Shi, S. Q.; Zheng, G. P. Tensile strains give rise to strong size

effects for thermal conductivities of silicene, germanene and stanene. Nanoscale 2016, 8 (6), 3760-3767.

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

67.

Page 22 of 33

Peng, B.; Zhang, H.; Shao, H.; Xu, Y.; Zhang, X.; Zhu, H. Low lattice thermal

conductivity of stanene. Sci Rep-Uk 2016, 6, 20225. 68.

Zhang, Y. Y.; Pei, Q. X.; Jiang, J. W.; Wei, N.; Zhang, Y. W. Thermal conductivities of

single- and multi-layer phosphorene: a molecular dynamics study. Nanoscale 2016, 8 (1), 483491. 69.

Xu, W.; Zhu, L.; Cai, Y.; Zhang, G.; Li, B. Direction dependent thermal conductivity of

monolayer phosphorene: Parameterization of Stillinger-Weber potential and molecular dynamics study. J Appl Phys 2015, 117 (21), 214308. 70.

Zhu, L. Y.; Zhang, G.; Li, B. W. Coexistence of size-dependent and size-independent

thermal conductivities in phosphorene. Physical Review B 2014, 90 (21). 71.

Jain, A.; McGaughey, A. J. H. Strongly anisotropic in-plane thermal transport in single-

layer black phosphorene. Sci Rep-Uk 2015, 5, 8501.

22


Page 23 of 33

(a)

(d)

zigzag

TABLE AND FIGURES

zigzag

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


MoSe2

MoS2

x y

z (b)

top

bottom (c)

x

top

bottom

armchair

y

z (e)

armchair

top

bottom (f)

top

bottom

Figure 1. (a) Top view of the monolayer MoSe2 structure. Zigzag boundary is along the x direction and armchair boundary is along the y direction. (b) Front view of monolayer MoSe2 in the x direction. The top and bottom Se atoms are treated as two different atomic types in the SW potential to generate an accurate description of interatomic interactions. (c) Side view of monolayer MoSe2 in the y direction. (d) Top view of the monolayer MoS2 structure. (e) Front view of monolayer MoS2 in the x direction. The top and bottom S atoms are treated as two atomic types as well. (f) Side view of monolayer MoS2 in the y direction.

23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

Qout

Qin

Heat flux

x z

Page 24 of 33

y

heating

cooling

fixed

(b)

Figure 2. (a) Schematic plot to illustrate the NEMD simulation method. The outmost chains of black atoms are fixed as denoted in black. Periodic boundary condition is employed in the x (width) direction. Free boundary condition is applied to the y (heat flux) and z (out-of-plane) directions. Red and blue atoms indicate heat bath and heat sink region, respectively. (b) Temperature gradient along the heat flux direction of 43.15 nm long armchair monolayer MoSe2 nanoribbon. Red line stands for the linear fitting result of MD temperature for each atom (black dots). Atomic configuration after the heat flux reached the steady state is shown in the inset.

24


Page 25 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Atomic velocity distributions of 163.46 nm long armchair MoSe2 nanoribbon. (a) Velocity distribution after 500 ps NVT and 500 ps NVE simulations. (b) Velocity distribution after 2 ns NEMD simulations. MD simulation results well match the Maxwellian velocity distributions, indicating system has reached steady state before data collection.

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 33

(a)

(b)

Figure 4. Each data point represents an average of results from five independent MD simulations. The error bar is calculated using the standard deviation. (a) Length dependence of thermal conductivities for monolayer MoSe2 and MoS2 nanoribbons in armchair and zigzag directions. (b) Width dependence of thermal conductivity for 43.15 nm long armchair MoSe2 nanoribbon. Results converged at 5.90 W/m·K, indicating thermal conductivity is not affected by the system width due to the periodic boundary condition.

26


Page 27 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

(b)

Figure 5. Each data point represents an average of results from five independent MD simulations. The error bar is calculated using the standard deviation. (a) Linear fitting of 1/κ and 1/L for armchair and zigzag MoSe2 nanoribbons based on 5 data points and 4 data points, respectively. (b) Linear fitting of 1/κ and 1/L for armchair and zigzag MoS2 nanoribbons based on 5 data points and 4 data points, respectively.

27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 33

(a)

(b)

Figure 6. (a) Temperature dependence of thermal conductivity for ~160 nm long monolayer MoSe2 and MoS2 nanoribbons in armchair and zigzag directions. Fitting results by the inverse relationship with temperature ( ~ 1/T) are plotted with solid lines. (b) Energy dependence of thermal conductivity for 43.15 nm armchair MoSe2 nanoribbons. Results converged at 6.01 W/m·K, indicating thermal conductivity is independent of heat flux.

28


Page 29 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5×5

10×10 15×15 20×20 30×30 50×50 60×60 Unit cell

Figure 7. Unit-cell size dependence of thermal conductivities for monolayer MoSe2 and MoS2 sheets in armchair and zigzag directions in GKM simulations. Averaged thermal conductivities over armchair and zigzag directions are shown in solid lines, with vertical bars to give a clear view of size convergence. Thermal conductivities of monolayer MoSe2 and MoS2 are both converged at 50×50 unit cells.

29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 33

(a)

(b)

(c)

(d)

Figure 8. Time dependence of HCAFs for monolayer MoSe2 and MoS2 sheets in armchair and zigzag directions in GKM simulation. The HCAFs decay to zero at the end of 800 ps.

30


Page 31 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

(b)

(c)

(d)

Figure 9. Time dependence of thermal conductivities for monolayer MoSe2 and MoS2 sheets in armchair and zigzag directions in the GKM simulation. Dashed lines stand for the thermal conductivities of ten independent trials and the thick solid lines represent the overall thermal conductivities averaged over ten individual trials, and the vertical dotted line indicates where thermal conductivity starts to converge and from which the data are used for the prediction of final thermal conductivity of monolayers.

31



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 33

Table 1. A summary of previous classical MD simulations for computing thermal conductivities of popular 2D monolayer materials at 300 K, compared with first-principle results and experimental values. 2D Material

Reference

Method

MoSe2 MoSe2 MoSe2 MoS2 MoS2 MoS2 graphene graphene graphene graphene graphene graphene graphene graphene h-BN h-BN h-BN h-BN h-BN silicene silicene silicene silicene silicene stanene stanene stanene phosphorene phosphorene phosphorene phosphorene phosphorene

current work Zhang et al.29 Gu et al.32 current work Zhang et al.29 Gu et al.32 Hong et al.48 Ng et al.49 Zhang et al.50 Li et al.51 Hu et al.52 Evans et al.53 Lee et al.54 Kong et al.55 Tabarraei et al.56 Mortazavi et al.57 Lindsay et al.58 Wang et al.59 Jo et al.60 Pei et al.61 Zhang et al.62 Zhang et al.62 Guo et al.63 Xie et al.64 Cherukara et al.65 Kuang et al.66 Peng et al.67 Hong et al.48 Zhang et al.68 Xu et al.69 Zhu et al.70 Jain et al.71

Classical MD simulation Optothermal Raman technique First-principles Classical MD simulation Optothermal Raman technique First-principles Classical MD simulation Classical MD simulation Classical MD simulation Classical MD simulation Classical MD simulation Classical MD simulation Raman scattering spectroscopy First-principles Classical MD simulation Classical MD simulation First-principles Microbridge Microbridge Classical MD simulation Classical MD simulation Classical MD simulation First-principles First-principles Classical MD simulation First-principles First-principles Classical MD simulation Classical MD simulation Classical MD simulation First-principles First-principles

κ (W/m·K) armchair zigzag 42.31 41.96 59 ± 18 54 103.00 106.18 84 ± 17 103 1008.5 1086.9 53.6 78.0 317 ~100 ~2000 ~6000 ~1800 2200 277.78 588.24 80 ∼520 484 ~360 41 8.64 11.77 32 9.4 6.4 ± 0.3 2.8 ± 0.2 1.25 11.6 63.6 110.7 9.891 42.553 33 152.7 24.3 83.5 36 110

32


Page 33 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC GRAPHIC

Heat flux

33


Thermal Conductivity of Monolayer MoSe2 and MoS2 - The Journal of

Recommend Documents