Interfacial Thermal Conductance across Graphene / MoS 2 van der Waals Heterostructures

: The thermal conductivity and interface thermal conductance of graphene stacked MoS 2 (graphene / MoS 2 ) van der Waals heterostructure were studied by the ﬁrst principles and molecular dynamics (MD) simulations. Firstly, two di ﬀ erent heterostructures were established and optimized by VASP. Subsequently, we obtained the thermal conductivity ( K ) and interfacial thermal conductance ( G ) via MD simulations. The predicted K of monolayer graphene and monolayer MoS 2 reached 1458.7 W / m K and 55.27 W / m K, respectively. The thermal conductance across the graphene / MoS 2 interface was calculated to be 8.95 MW / m 2 K at 300 K. The G increases with temperature and the interface coupling strength. Finally, the phonon spectra and phonon density of state were obtained to analyze the changing mechanism of thermal conductivity and thermal conductance. the vDOS and phonon spectra of graphene and MoS 2 have been further analyzed for a better understanding of the a ﬀ ecting factors on G and K . The out-plane phonons in graphene transport most thermal energy. The increasing coupling strength X and temperature can improve thermal conductance.


Introduction
Due to high carrier mobility [1], high thermal conductivity [2], and ultra-high Young's modulus [3], graphene has been widely used in the fields of energy storage [4], optoelectronics [5], and sensors [6]. The MoS 2 as one kind of transition metal sulfides (TMDs), which overcomes graphene's zero-bandgap, has been used as a catalyst to generate hydrogen [7] and in photothermal conversion area [8]. The heterostructures is the two-dimensional (2D) layered materials stacked by van der Waals (vdW) forces. It can be integrated into complex devices applied in photodetectors [9], photocatalytic, molecular sieves, and electrodes [10]. Ren et al. [11] reported that the TMDs/BN (boron nitride) vertical heterostructures have a direct bandgap and excellent optical properties. It can be applied to photocatalytic, photovoltaic, and optical areas. Chen et al. [12] fabricated MoTe 2 /MoS 2 hetero-structures. It has a fast response and a broad wavelength range and can be applied to photodetection and on-chip logic circuits. Li et al. [13] prepared graphene/silicene heterostructures and proved weak van der Waals interactions between layers. The vertical graphene/silicene/Ru(Ruthenium) heterostructures show well Schottky rectification behavior.
The electrical and optical properties of graphene/MoS 2 heterostructures have been studied. Biroju et al. [14] prepared monolayer graphene with a single or a few layered MoS 2 van der Waals vertical heterostructures and studied the catalytic activity mechanism by the density functional theory (DFT). The heterostructures that graphene is over the MoS 2 showed higher photo-response and hydrogen evolution reaction activity. The work of Safeer et al. [15] showed the spin effect in the graphene/MoS 2 heterostructure and proved a superimposed spin-to-charge current conversion in the device based on it. Sun et al. [16] reported that the saturable absorption of graphene/MoS 2 hetero-structure is adjustable by changing the thickness of MoS 2 . Furthermore, the structure can be applied in solid-state lasers, which can generate ultrafast pulses.
The development of applications with graphene/MoS 2 vdW heterostructures requires the thermal management of these structures. Gong et al. [17] developed a model of monolayer graphene/Li/monolayer MoS 2 . The heterostructures' electrical conductivity is higher than that of the MoS 2 , and the thermal conductivity of graphene/monolayer MoS 2 with three graphene layers can maintain 85.5 W/m K at 100~500 K by molecular dynamics. Zhang et al. [18] simulated MoS 2 /graphene hybrid nanosheet (MGHN) structures and calculated thermal conductivity by molecular dynamics. In the simulation, heat baths are applied to both the MoS 2 and graphene layers. The results show that the graphene layer in the MGHN structure can transfer most heat. It has a very high thermal conductivity and can be tuned by the interlayer coupling strength, temperature, and contact area. Oh et al. [19] stacked the fingerprint-like randomly aligned graphene nanoribbons on the monolayer MoS 2 by block copolymer (BCP) lithography. The electrical conductivity and power factor of the hetero-structure are 902 S m −1 and 222.1 µW m −1 K −2 . The hetero-junction network has superior thermoelectric properties and may be applied in next-generation electronic and energy devices. Therefore, understanding of van der Waals heterostructures' thermal transport is essential for future use in optoelectronic and thermo-electronic applications. Previous researches investigated the thermal properties of graphene/MoS 2 with different lengths, layers, and strain by experiment and simulation. However, few focused on different initial cell structures and defect effects on heterostructure thermal conductivity.
In this work, two graphene/MoS 2 van der Waals heterostructures were established and investigated by the density functional theory (DFT). The bilayer nanoribbons were established, then the thermal conductivity and thermal conductance were obtained by molecular dynamics. We analyzed the thermal conductivity with different lengths, temperature, defect type, and defect concentration. The interfacial thermal conductance across the interface with different temperature and interface coupling strengths were calculated subsequently. We also obtained the phonon spectra to analyze defects and the coupling strength effects on the bilayer system thermal transports.

DFT Calculations
The establishment of geometries and energy calculation are based on DFT performed by software VASP (Vienna ab-initio Simulation Package) [20]. Here we used the projector-augmented wave (PAW) method to describe the electron-ion interactions. GGA-PBE [21] exchange-correlation function was selected for describing the interactions, and DFT-D2 [22] method was applied to describe van der Waals forces between layers. We stacked a monolayer of graphene whose supercell is 5 × 5 × 1 and a monolayer of MoS 2 whose supercell is 4 × 4 × 1. The lattice constants of monolayer graphene and MoS 2 are shown in Table 1. As in the previous First-principles study, the parameters are the same and agree with the experimental data. Two structures are selected; for the first structure (Structure 1) shown in Figure 1a, the lattice constant of a H was taken as 12.26 Å [23], which is the optimized graphene lattice constant (a G ). The lattice constant of MoS 2 (a M ) and the interlayer spacing were taken as 3.12 Å and 3.37 Å, respectively. A vacuum layer of 25 Å is added to reduce the interaction between the periodic structures in the Z-direction. And, for the second structure (Structure 2), a G = 2.47 Å, and a M = 3.16 Å [24]. The lattice mismatch ratio of it was about 2.3%. The lattice constant of aH is 12.46 Å after optimized, and the distance between layers shown in Figure 1c is 3.32 Å. The cut-off energy of the plane wave was taken as 550 eV for structure 1, and then another is 500 eV. The electronic iterations convergence was determined as 1.00 × 10 −5 eV using the Normal (blocked Davidson) algorithm, and the structure relaxations finished when the convergence precision of each interatomic force reached 0.02 eV/Å. The 3 × 3 × 1 and 6 × 6 × 1 Monkhorst-Pack characteristic K points were adopted for geometries optimization. The building energy (Eb) per carbon atom between graphene and MoS 2 is calculated as: where, E G/M , E G , E M are energies of the heterostructure, isolated monolayer graphene, and monolayer MoS 2 . N C is the total number of carbon atoms in the cell (50 in this simulation). For structure 1, E b = −33.3 meV/atom; for structure 2, E b = −34.1 meV/atom. Due to the different exchange-correlation function and van der Waals forces correction, these values are higher than the previous study E b = −25.1 meV/atom [30]. That also can be found in the previous study of monolayer MoS 2 stacked on the graphene bilayer system [31].

Molecular Dynamics Simulations
The simulation was performed by the large-scale atomic molecular massively parallel simulator (LAMMPS) package [32]. The optimized Tersoff potentials [33] was adopted to describe C-C interactions in graphene because it can reproduce the correct phonon dispersion and the density of states. We used Stillinger-Weber (SW) potential [34] and SW15 parameters [35], which can accurately describe the thermal and mechanical properties of monolayer MoS 2 . The 12-6 Lennard-Jones (LJ) potential was adopted to describe van der Waals interactions between graphene and MoS 2 atoms: where, x, ε, σ, and r are the interface interaction strength, energy parameter, distance parameter, and the distance between the atoms, respectively. In this study, the default setting of x is 1, ε Mo-C = 1.39 meV, σ Mo-C = 3.8 Å, ε S-C = 7.53 meV, σ S-C = 3.22 Å, and r = 1 nm [36]. The non-equilibrium molecular dynamics (NEMD) method was adopted to calculate the heterostructures' thermal conductivity shown in Figure 2. The length of the rectangle bilayer system is 21.56 nm, and its width is 4.98 nm. Periodic boundary conditions were applied in yand zdirections. In Figure 2a, the black region with a width of 2.56 Å contains fixed atoms. The left region is the heat source with a width of 1.53 nm, and the blue region is the heat sink with a width of 1.53 nm. Initially, the system equilibrated at a temperature of 300 K, pressure of 1 bar using Nosé-Hoover thermostat under the isothermal-isobaric (NPT) ensemble for 1 ns with the time step of 1 fs. Subsequently, the thermostat was removed and changed into the canonical (NVT) ensemble at a temperature of 300 K for 1 ns. Finally, the system was switched to the microcanonical ensemble for 5 ns. Simultaneously, the Langevin thermostat was applied to maintain heat source temperature at 310 K and heat sink temperature at 290 K. The system equilibrated after 1 ns, and the energy profile was obtained. According to Fourier's law, the thermal conductivity (K) is obtained by: where, J is the heat flux, A is the cross-sectional area equivalent to width multiplying thickness, and ∂T/∂X is the temperature gradient in the X direction. In this study, the thickness of the heterostructure is 9.53 Å, which is a sum of monolayer graphene's thickness of 3.44 Å and monolayer MoS 2 thickness of 6.09 Å. A system with different point defects shown as Figure 2b has also been calculated. The concentration of defect (C) is defined as: where, N D is the number of defect atoms in graphene or MoS 2 layers, N G denoted the total number of atoms in graphene or MoS 2 , which is 4000 and 3840 in this paper. The thermal conductance (G) across the graphene/MoS 2 interface has been obtained by the thermal relaxation method [37]. Due to its small size and high thermal conductivity, the lumped heat-capacity model can be applied to calculate the thermal conductance (G), where, C H is the effective constant volume heat capacity of the heterostructure. C H can be expressed as where, C G and C M are the effective constant volume heat capacity of monolayer graphene and MoS 2 , which can be written as where, C i,µ denoted the contribution of µ atom in i direction to the heat capacity. k B is Boltzmann constant, h is Planck constant, T is temperature in Kelvin, and g(ω) is phonon spectrum power at the frequency ω, which can be expressed as where, V j (t) is the velocity of the atom j at time t. τ is the relaxation time, which can be calculated by where, t 0 is origin time, ∆T is the temperature difference between graphene and MoS 2 . After the system finally reached equilibrium, the graphene was rapidly heated to 500 K by velocity rescaling. The MoS 2 layer was still kept at 300 K. The Langevin thermostats were used to maintain graphene at 500 K and MoS 2 at 300 K for 1 ns. Finally, the heat baths were removed, and the system reached equilibrium for another 1 ns.

Result and Discussion
3.1. Thermal Conductivity of the Heterostructure Figure 3 shows the temperature distribution and energy variation of structure 1, which is stacked in Figure 1a. The heterostructure is nanoribbons with a length (L) of 21.56 nm and width (W) of 4.98 nm. Avoid the heat bathes and boundary effects, the middle part of the system was selected for temperature fitting. The temperature gradients of graphene, MoS 2 , and graphene/MoS 2 bilayer by linear fitting are 0.01902, 0.04633, and 0.03255 K/Å. All the energies subtracted from the heat source and added to the heat sink increase proportionally with the time, as shown in Figure 3b. The energy slope of graphene, MoS 2 , and total bilayer system are 1482.52, 238.64, and 1721.15 eV/ns, respectively. According to Fourier's law mentioned above, the thermal conductivities of graphene in the bilayer (K GH ), MoS 2 in heterostructure (K MH ), and the heterostructure are 739.26, 27.6, and 181.03 W/m K, respectively. All the values are lower than 821.28 W/m K, which is the thermal conductivity of monolayer graphene (K G ). The result of MoS 2 is lower than 34.5 W/m K measured experimentally [38]. Because the phonon mean free path (MFP) of graphene is 775 nm [39], and the simulation model size is much smaller than it. The thermal conductivity of graphene is much lower than 3080-5150 W/m K with a large size [39]. Similar results also can be found in the previous study for the graphene/MoS 2 /graphene structure [40]. Compared with the free-standing monolayer graphene, the K of heterostructure is significantly decreased for two reasons. One is the bilayer's area, which is thicker than monolayer graphene. According to Fourier's law, the area of the system can affect thermal conductivity. Another is that the K of MoS 2 is much lower than graphene. In the vdW heterostructure, the thermal transport in MoS 2 is much slower than graphene. At the same heat flux, the bilayer system thermal transport is inefficient, so that the thermal conductivity of it is lower than monolayer graphene. For structure 2 stacked as Figure 1b, the K GH , K MH , and K H are 649.89, 21.985, and 154.74 W/m K, respectively. Table 2 presents the comparison of K and G of Graphene/MoS 2 obtained in this study with other structures. Since the thermal resistance is actually connected in parallel, the thermal conductivity of heterostructure is lower than that of the monolayer. The higher thermal conductivity materials can transfer more heat, so the K of the vertically stacked heterostructure is higher underlying the same substrate.  The K of heterostructure (116~355 W/m K) with a fixed width of 4.98 nm is lower than that of monolayer graphene (601~1189 W/m K) but higher than that of monolayer MoS 2 (19~36 W/m K), as shown in Figure 4a. Because the MFP of MoS 2 is about 18.1 nm, both K of monolayer MoS 2 and K of the MoS 2 in heterostructure change mildly, as shown in Figure 4a. According to the Casimir Limit [50], thermal conductivity values have a strong size dependence. The relation between K −1 and L −1 is depicted in Figure 4b, and the K ∞ of the system (infinite length) can obtain by: where, λ is MFP and L is the distance between the heat source and heat sink. K ∞ of the heterostructure is correlated to be 671.14 W/m K, and the MFP is 539.17 nm that between monolayer graphene and MoS 2 . K ∞ of MoS 2 is 55.27 W/m K that approaching the previous 39.8 W/m K, which is the thermal conductivity of armchair monolayer MoS 2 with a width fixed at 2 nm [51]. The K of monolayer graphene with width 4.98 nm and infinite-length is 1458.7 W/m K. It is close to 1500 W/m K [52] that is also obtained by applying the optimized Tersoff potentials.  Figure 5a presents the K of graphene, MoS 2 , and bilayer system at different temperatures. The initial model structure 1 (21.2 nm × 4.9 nm) was simulated. It is observed that all these three thermal conductivities decrease with the temperature. However, the slopes for different structures vary very much. K of graphene decreases very quickly, while K of MoS 2 only changes slightly. Because the Debye temperature of graphene is 2100 K [53] and that of MoS 2 is 263.2 K, the quantum correction of thermal conductivity is not required in this study. The frequency of phonon increases with temperature, and the phonon collisions increase at the same time. Subsequently, the MFP of graphene gets shorter, and the K of graphene declined obviously. Figure 5b reveals that the defect in the monolayer can affect the K of the bilayer system significantly. The simulation results are for structure 2 (21.5 nm × 4.98 nm), and the styles of single and double point defects are shown in Figure 2b. It clearly shows that the single point vacancy in graphene has the most apparent impact on the thermal conductivity of bilayer. The defects of MoS 2 have not many effects because the graphene layer plays a paramount role in the thermal transport of the bilayer. The K H decreases 68.7% when single vacancy concentration is 0.5%, while it only decreases by 39% when double vacancies rise to 0.52%. Because of their less stable two-coordinated atoms, the single vacancies reduce the thermal conductivity much more effectively [54].
The atoms do not follow the normal pattern of vibrations. Thus it causes a higher degree of scattering than double vacancies.  Figure 6a shows the temperature evolution of bilayer after Langevin thermostats were removed. By nonlinear fitting, ∆T = 203.16e −t/108.02 for the initial structure 1 stacked as Figure 1a shown. For structure 2, τ can be obtained as 105.04 ps. C H at temperature 300 K is 7.74 J·K −1 mol −1 that can be calculated from Figure 6c. So the thermal conductance of graphene/MoS 2 is 8.95 MW/m 2 ·K for model 1 and 9.2 MW/m 2 ·K for the other model. It is found that the initial model has little effect on thermal conductance. That is higher than previous investigations on the thermal conductance of the bilayer graphene/MoS 2 system, which is 5.8 MW/m 2 ·K [49]. Moreover it is lower than of graphene/MoS 2 (13.8 MW/m 2 ·K) [25] and higher than that of monolayer MoS 2 on the Au substrate (0.44 MW/m 2 ·K) [55].    It can be observed in Figure 8a that the overlapped region of graphene and MoS 2 are mostly in the frequency range from 1 THz up to 15 THz. The low-frequency phonons in graphene (out-plane) and MoS 2 play a significant role in thermal transport. These phonons can provide more channels for phonon transport, and lower frequency phonons transport most heat in the bilayer system. When the temperature rises, most phonons are excited, and Umklapp scattering is enhanced at the same time. More channels are provided, and the relaxation time τ is decreased, so the G increases with the temperature. When X increases, the efficiency of thermal transport is improved, and the G would be enhanced [49]. At the same time, the phonon in the MoS 2 layer can be much easier obtained from out-plane graphene to in-plane graphene, and the thermal transport between layers is enhanced. The interface thermal conductance would be increased as a consequence of it.

Thermal Conductance of Graphene/MoS 2 Heterostructure
The vDOS of MoS 2 with different defects shown in Figure 8b,c are obtained to explain their effect on the thermal conductivity in detail. The vDOS of MoS 2 with V S1 or V S2 point vacancy is very close to pristine MoS 2 . It can be observed in the enlarged box in Figure 8b that the vDOS of MoS 2 with Vs 2 is slightly less than that of V S1 . Therefore, the K H is little changed, and V S1 vacancy in MoS 2 has little influence on K H , as shown in Figure 5b. With the concentration of defect increasing, the DOS of phonon in MoS 2 is decreased, and the vibrational modes are suppressed. Hence the relaxation time and phonon MFP decrease at the same time. The K of the system decreases with the concentration of vacancy. It has less influence on K H than graphene owing to the relatively minor role played in thermal transports.
The phonon spectra shown in Figure 9 are obtained by the GULP [56]. The optimized Tersoff and SW potentials have been used to describe the interaction force. In Figure 9, the maximum wavenumber of monolayer graphene is around 50 TH Z , and that of MoS 2 is around 15 TH Z . This result has also been reflected in Figure 8, where the phonon density states of graphene are located in the high frequency around 50 TH Z , and that of MoS 2 cut off at around15 TH Z . Along Г to M direction in Figure 9, the TA, LA, and TO branches of graphene cut off at around 36 TH Z , and the ZO and ZA branch cut off at around 21 TH Z , corresponding to the PDOS of graphene shown in Figure 8a. The ZA, TA, and LA branches of MoS 2 cut off at around 5 TH Z along Г to M direction, so the enlarged area of PDOS shown in Figure 8b, c is selected from 4 TH Z to 6 TH Z . The results also show that the lower frequency phonons play a major role in thermal transport, which can also be found in Figure 8.

Conclusions
In this study, graphene/MoS 2 heterostructures have been established and optimized by density functional theory. A non-equilibrium molecular dynamics simulation method was applied to calculate the thermal conductivity and interfacial thermal conductance of the bilayer system. The results showed that K H increases with length, and the K of infinitely long graphene/MoS 2 (width = 4.98 nm) has been obtained to be 671.14 W/m K. The influencing factors of K H also have been analyzed carefully. It is revealed that K H increases with the temperature and the concentration of defect contained. The single vacancy in graphene can significantly decrease the K H . The thermal conductance of the heterostructure has been found to be 8.95 MW/m 2 ·K, which is much higher than that of MoS 2 on the substrate. Finally, the vDOS and phonon spectra of graphene and MoS 2 have been further analyzed for a better understanding of the affecting factors on G and K. The out-plane phonons in graphene transport most thermal energy. The increasing coupling strength X and temperature can improve thermal conductance.