The thermal properties of two-dimensional (2D) materials, such as MoS2, are known to be affected by interactions with their environment, but this has primarily been studied only with SiO2 substrates. Here, we compare the thermal conductivity (TC) and thermal boundary conductance (TBC) of monolayer MoS2 on amorphous (a-) and crystalline (c-) SiO2, AlN, Al2O3, and h-BN monolayers using molecular dynamics. The room temperature, in-plane TC of MoS2 is ∼38 Wm−1 K−1 on amorphous substrates and up to ∼68 Wm−1 K−1 on crystalline substrates, with most of the difference due to substrate interactions with long-wavelength MoS2 phonons (<2 THz). An h-BN monolayer used as a buffer between MoS2 and the substrate causes the MoS2 TC to increase by up to 50%. Length-dependent calculations reveal TC size effects below ∼2 μm and show that the MoS2 TC is not substrate- but size-limited below ∼100 nm. We also find that the TBC of MoS2 with c-Al2O3 is over twice that with c-AlN despite a similar MoS2 TC on both, indicating that the TC and TBC could be tuned independently. Finally, we compare the thermal resistance of MoS2 transistors on all substrates and find that MoS2 TBC is the most important parameter for heat removal for long-channel (>150 nm) devices, while TBC and TC are equally important for short channels. This work provides important insights for electro-thermal applications of 2D materials on various substrates.

Two-dimensional (2D) semiconductors are being investigated for three-dimensional integrated circuits (3D-ICs), opto-electronics, and flexible electronics,1–3 in large part due to their relatively good electron and hole mobility in a sub-nanometer thin material.4 Thermal properties of 2D materials are also important in this context as the numerous interfaces and poor thermal conductivity (TC) substrates (e.g., SiO2 and polyimide) present in electronic applications lead to self-heating, which decreases transistor performance5 and reliability6 during operation.7 The ultrathin nature of 2D materials is known to lead to a dependence of their electrical mobility8 and TC9 on so-called remote phonons or impurities, i.e., those belonging to the materials below and above the 2D semiconductor. For example, our recent calculations9 have shown that monolayer MoS2 supported by amorphous SiO2 suffers decreased TC compared to the freely suspended material and that encasing the 2D material further lowers its TC.

However, the effects of insulators other than SiO2 on MoS2 thermal transport remain unknown, as does the impact of insulator crystallinity. Moreover, the thermal boundary conductance (TBC) of the interface between MoS2 and such insulators is not well understood, especially with respect to its dependence on substrate type. The TBC at 2D material van der Waals (vdW) interfaces could be a stronger bottleneck for heat dissipation in electronics because it is relatively low, of the order ∼15 MWm−2 K−1 at room temperature,6,10,11 which is equivalent to the thermal resistance of nearly 100 nm of SiO2 (i.e., the Kapitza length).

Here, we use atomistic molecular dynamics (MD) to calculate the TC and TBC of MoS2 supported by SiO2, AlN, and Al2O3, all technologically relevant insulators with a wide range of TCs themselves. We examine both amorphous and crystalline forms of the substrates to quantify their effect on the MoS2 thermal properties. In addition, because electrical properties of 2D materials are known to improve with a hexagonal boron nitride (h-BN) interfacial layer,12,13 we also study the thermal consequences of adding an h-BN layer between MoS2 and each substrate. The atomistic simulations enable a direct comparison of thermal properties between structures and provide insights into the frequency dependence of TC, the phonon mean free path, and length-dependent TC. Finally, we develop an analytical model to evaluate the effectiveness of heat removal from simple, back-gated MoS2 transistors of various geometries and compare the impact of different substrates by using our calculated TCs and TBCs for supported MoS2.

We use MD for all calculations of TC and TBC in this study. While the TC of crystalline materials can be described well by the Peierls–Boltzmann transport equation paired with density functional theory,14 these methods currently struggle to model systems with defects, interfaces, and amorphous materials.15 In contrast, MD can calculate the TC of spatially complex structures without any modifications or approximations.16–18 It also naturally incorporates all anharmonicities of a solid and does not make assumptions about the dynamics of the system.15,19 This makes MD ideal for realistic TBC calculations, especially for supported 2D materials where concepts of cross-plane phonon group velocities are not well defined.20 

All the results in this paper are calculated using the Graphics Processing Units Molecular Dynamics (GPUMD-v2.5.1) package.21–25 The atomic interactions within and between materials are as follows: MoS2 is modeled by the reactive empirical bond-order potential with a Lennard-Jones (LJ) addition (REBO-LJ),26–28 which has been shown to predict thermal properties accurately;29 SiO2 is modeled by the Tersoff potential30 parameterized by Munetoh et al.;31 AlN32 and Al2O3 (Ref. 33) are modeled by the Vashishta potential; h-BN is modeled by the Tersoff potential parameterized by Sevik et al.;34 and all inter-material, vdW interactions are modeled by the LJ potential with Lorentz–Berthelot mixing rules (see Sec. S1 in the supplementary material). Unless otherwise stated, the time step for each simulation is 0.5 fs and the temperature is 300 K.

Figures 1(a)1(f) display the MoS2-with-substrate combinations considered in this work, including both amorphous and crystalline substrates. The structures with an h-BN monolayer inserted between the MoS2 and substrate are shown in Fig. S1 of the supplementary material. For all simulations, we use a 4080 atom MoS2 sheet with an area of 10.8 × 11 nm2, which has been shown to be of sufficient size for TC calculations with periodic in-plane boundary conditions.9,29 We fix the lateral dimensions of the simulation cell and construct all substrates to fit into those dimensions. This ensures that MoS2 is strain-free regardless of the substrate, as a strain of 1% could change the TC by up to 5%.35 We create supported MoS2 structures by placing MoS2 on each substrate (with or without h-BN) and minimizing the energy of each heterostructure to establish the correct vdW distances between materials. The monolayer thickness of MoS2 is 6.15 Å36 and that of h-BN is 3.33 Å.37 The a-SiO2, a-AlN, a-Al2O3, c-SiO2 (quartz), c-AlN, and c-Al2O3 (sapphire) substrates are 2.7, 3.5, 2.9, 2.7, 2.95, and 2.5 nm thick, respectively, as seen in Fig. 1. For details about structure creation, including information on the surface orientation of each substrate, please see Sec. S2 in the supplementary material.

FIG. 1.

(a)–(f) Visualizations of the structures used for supported MoS2 thermal conductivity and thermal boundary conductance calculations. The area of each structure is 10.8 × 11 nm2, and the relative thicknesses can be determined from the scale bar in (a). The in-plane and out-of-plane nature of thermal conductivity and thermal boundary conductance is highlighted by the red cartoons in (b). (g) The thermal conductivity of monolayer MoS2 when suspended (blue bar, left) and supported by crystalline (c-) and amorphous (a-) AlN, Al2O3, and SiO2. The gray horizontal lines above the bars represent the thermal conductivity of MoS2 if a single layer of h-BN is inserted between the MoS2 and each substrate. All calculations are at 300K.

FIG. 1.

(a)–(f) Visualizations of the structures used for supported MoS2 thermal conductivity and thermal boundary conductance calculations. The area of each structure is 10.8 × 11 nm2, and the relative thicknesses can be determined from the scale bar in (a). The in-plane and out-of-plane nature of thermal conductivity and thermal boundary conductance is highlighted by the red cartoons in (b). (g) The thermal conductivity of monolayer MoS2 when suspended (blue bar, left) and supported by crystalline (c-) and amorphous (a-) AlN, Al2O3, and SiO2. The gray horizontal lines above the bars represent the thermal conductivity of MoS2 if a single layer of h-BN is inserted between the MoS2 and each substrate. All calculations are at 300K.

Close modal

For TC calculations along the MoS2, we use the homogeneous nonequilibrium MD (HNEMD) method,24,38 which has been shown to work well for supported 2D materials.9 This method enables a direct calculation of TC, but, in practice, we calculate the running average of the TC with24 

κ(t)=1t0tJ(τ)neTVFedτ,
(1)

where T is the system temperature, V is the system volume, Fe is the driving force parameter [simplified to a scalar due to the isotropic in-plane thermal conductivity of MoS2 (Ref. 29)], and ⟨J(τ)⟩ne is the non-equilibrium heat current due to the driving force. This non-equilibrium heat current can be spectrally decomposed using the spectral heat current (SHC) method.24,25,39,40 The frequency-domain, spectral TC can be written as

κ(ω)=2K~(ω)TVFe,
(2)

where K~(ω) is the Fourier transform of the virial-velocity correlation as defined in Refs. 24 and 25. The SHC method is also more robust than other spectral TC methods when considering systems with more anharmonicity.25 Note that the in-plane TC can be decomposed into contributions from in-plane atomic motion (dominant in longitudinal and transverse acoustic phonons) and out-of-plane atomic motion (dominant in flexural, or ZA, acoustic phonons).23,24 Simulation and procedural details of the HNEMD and SHC methods can be found in Sec. S3 of the supplementary material.

The SHC method calculates the spectral TC if the structure has no temperature gradient and is large enough for diffusive thermal transport; however, it can also calculate the thermal conductance G(ω)23,41 during NEMD simulations (which have large temperature gradients) with ballistic thermal transport. The details of our G(ω) calculations and NEMD simulations are in Sec. S6 of the supplementary material. With both κ(ω) and G(ω), we also calculate size-related thermal properties, such as the spectral phonon mean free path (MFP), λ(ω) = κ(ω)/G(ω),24 and, subsequently, the length-dependent TC,24,41

κ(L)=0dω2πκ(ω)1+λ(ω)/L.
(3)

Here, L is the length between the two thermal baths.

In addition to the in-plane TC, the cross-plane thermal transport may also change drastically when MoS2 is supported by different substrates. To investigate this, we calculate the TBC between the MoS2 sheet and its substrate using the approach to the equilibrium MD (AEMD) technique.42–44 In these simulations, we allow previously heated MoS2 to cool into its substrate, and we track the temperature difference (ΔT = TMoS2 − Tsub) between the two over time t. Modeling each structure as a resistive-capacitive (RC) thermal circuit, we extract the TBC G using

ΔT(t)=ΔT0e(1CMoS2+1Csub)AGt,
(4)

where A is the contact area of the MoS2 with its substrate, ΔT0 is the initial temperature difference, and C is the heat capacity with subscripts denoting the respective material. Our AEMD setup and simulation details can be found in Sec. S7 of the supplementary material.

Based on the models we have chosen for our MD simulations, there are three areas to consider when interpreting the TC and TBC of supported MoS2: the vdW interactions, the mass differences, and the structure of each material. First, the vdW bonds, facilitated through the LJ potential, mediate any interaction between MoS2 and a substrate. (Note that we are not considering other interfacial effects such as chemical bonding,45 but our materials form otherwise atomically intimate contacts.) The strength of vdW bonds significantly influences the thermal properties of MoS2; an increase in vdW strength increases TBC and decreases TC.9,44,46 Our choice to use materials with elements next to each other on the periodic table (i.e., Al—Si, and N—O) results in LJ parameters that are similar across each of the substrates and, speaking to our second effect, similar elemental masses in each substrate. Thus, any significant changes in the TC of supported MoS2 likely do not come from differences in the strength of the vdW bonds or atomic masses in this study.

This leaves the atomic structure as the primary differentiator when interpreting the TC and TBC of supported MoS2 in this study. We know that amorphous materials have vastly different vibrational modes than their crystalline counterparts, resulting in vastly different interactions and phonon scattering events with MoS2. We also know that thermal properties of different crystalline substrates can vary greatly too, with the TC of c-SiO2 (Ref. 47) being over an order of magnitude smaller than that of c-AlN.48,49 The TCs of amorphous substrates tend to be similar, typically near 1–2 Wm−1 K−1.50–53 Which differences and similarities in the substrates matter to the thermal properties of MoS2 will be central to our discussion.

1. Total thermal conductivity

We first calculate the TC of freely suspended monolayer MoS2 at room temperature, which serves as a reference value for all other supported MoS2 structures. The TC of suspended MoS2 is 122.0 ± 1.45 Wm−1 K−1 with respective contributions from in-plane and out-of-plane atomic motion of 90.7 ± 1.3 Wm−1 K−1 and 31.3 ± 0.8 Wm−1 K−1. This value agrees well with recent suspended monolayer MoS2 TC measurements, which are around 100 Wm−1 K−1 (Refs. 11 and 54) and with previous calculations.9,29

We then calculate the TC of MoS2 supported by amorphous and crystalline SiO2, AlN, and Al2O3, as well as those same structures with an h-BN interlayer between MoS2 and the substrate. (See Figs. 1 and S1 in the supplementary material for structure visualizations, respectively.) The results of these calculations are shown in Fig. 1(g), and all TCs are listed in Table S2 of the supplementary material. We find that the TC of MoS2 supported by any amorphous substrate is approximately the same, in the range of ∼36–39 Wm−1 K−1. However, the TC of MoS2 supported by crystalline substrates can be notably different. When MoS2 is on c-AlN or c-Al2O3, its TC is in the range of ∼66–68 Wm−1 K−1, whereas MoS2 has a TC of only ∼42 Wm−1 K−1 when supported by c-SiO2 (quartz). Of these structures, only MoS2 on a-SiO2 has been studied experimentally, with a reported MoS2 TC of 63 ± 22 Wm−1 K−1,55 although a sputtered Ni capping layer is known to damage monolayer MoS2,45 likely affecting their TC measurement.

Next, we consider the effect of adding a single layer of h-BN between MoS2 and each substrate. The gray, horizontal lines in Fig. 1(g) show the new TC of MoS2 recalculated in these configurations. We find that adding an h-BN interlayer increases the TC of MoS2 in all cases. The largest increase of TC—by ∼20 Wm−1 K−1—is seen in the systems on SiO2 after the addition of the h-BN interlayer. For the other substrates, MoS2 sees a larger increase in TC when h-BN is placed on amorphous substrates (by ∼12 Wm−1 K−1) compared to on crystalline substrates (by ∼3.7–8.7 Wm−1 K−1). Interestingly, the h-BN interlayer prevents the substrates from reducing the in-plane contributions to TC, which results in the larger total TCs of MoS2. The out-of-plane contributions are approximately the same as without the h-BN interlayer or lower. Previous calculations suggested that the TC of graphene was higher when supported by an h-BN substrate than by an a-SiO2 substrate,56 indicating that other supported 2D materials could also have higher TCs if interfaced with h-BN.

2. Frequency dependence of thermal conductivity

By calculating the frequency-dependent (spectral) TC using the SHC method, we can better understand how each substrate affects MoS2. We first calculate the room temperature spectral TC of freely suspended MoS2, with the results shown in Fig. 2(a). Here, we only show the contributions from the acoustic modes (i.e., <8 THz) as optical mode contributions are negligible. (See Fig. S4 in the supplementary material for more details.) We find that the TC of suspended MoS2 is heavily influenced by low-frequency phonons, with 50% of the total TC contributions from frequencies below ∼2.1 THz. On the other hand, for substrate-supported MoS2, we find that those crucial, low-frequency contributions are the ones most greatly affected by the substrate, and the spectral TC contribution even decreases toward zero frequency as seen in Figs. 2(b)2(d). Part of this degradation is due to the various substrates almost completely suppressing TC contributions from out-of-plane atomic motion below ∼2 THz (i.e., long-wavelength flexural modes; see Fig. S4 in the supplementary material), which is where ∼50% of those contributions come from in suspended MoS2.

FIG. 2.

The frequency-dependent thermal conductivity κ(ω) when monolayer MoS2 is (a) suspended or supported by (b) SiO2, (c) AlN, and (d) Al2O3. The solid lines in (b)–(d) are for crystalline substrates, and the dotted lines are for amorphous substrates. The gray vertical dashed lines show the phonon frequency at which 50% of the total thermal conductivity has been contributed. All calculations are at 300K.

FIG. 2.

The frequency-dependent thermal conductivity κ(ω) when monolayer MoS2 is (a) suspended or supported by (b) SiO2, (c) AlN, and (d) Al2O3. The solid lines in (b)–(d) are for crystalline substrates, and the dotted lines are for amorphous substrates. The gray vertical dashed lines show the phonon frequency at which 50% of the total thermal conductivity has been contributed. All calculations are at 300K.

Close modal

From earlier discussion, we know that MoS2 supported by the various amorphous substrates has a similar total TC. The dotted curves in Figs. 2(b)2(d) reveal that the frequency dependence of MoS2 TC on each amorphous substrate is similar as well. [A direct comparison of these curves is in Fig. S7(a) of the supplementary material.] To explain this, we turn to the vibrational density of states (VDOS) of each material, which represents the density of vibrational modes (or phonon modes, if in a crystal) at a given energy57 and is shown in Fig. 3. The overlap of VDOS between two materials has been tied to the transmission of phonons at an interface in previous TBC studies20,58–62 and can be used here to understand the interaction between the MoS2 and its substrate, with the presumption that a smaller overlap of VDOS with MoS2 will result in a smaller TC degradation in MoS2.

FIG. 3.

Vibrational density of states (VDOS) of crystalline (line) and amorphous (shaded area) (a) SiO2, (b) AlN, and (c) Al2O3 compared to MoS2 at 300K. The MoS2 curves are calculated when supported by the amorphous substrate corresponding to the panel.

FIG. 3.

Vibrational density of states (VDOS) of crystalline (line) and amorphous (shaded area) (a) SiO2, (b) AlN, and (c) Al2O3 compared to MoS2 at 300K. The MoS2 curves are calculated when supported by the amorphous substrate corresponding to the panel.

Close modal

Interestingly, despite the similar spectral TCs of MoS2 supported by amorphous substrates, the VDOS of the amorphous substrates is different. We see that the VDOS of a-SiO2 is larger at frequencies below 2 THz showing that, for the case of amorphous substrates, a larger overlap in VDOS with MoS2 does not result in a lower MoS2 TC. Similar TCs may partially be due to our choice of substrate materials as they have atomic species of similar mass, similar vdW bond strength to MoS2, and, because all three are amorphous, similar atomic structures. Amorphous materials have a diverse set of vibrational modes consisting of propagons (sinusoidal, phonon-like modes), diffusons (delocalized, non-sinusoidal modes), and locons (localized vibrations),63,64 which can possess enough momentum and energy variation in their propagating and localized vibrations65,66 to maximize phonon scattering in MoS2 (i.e., the remote phonon scattering is not substrate limited). From an interaction perspective, these amorphous substrates likely look identical to MoS2 despite their VDOS differences. This position is confirmed by vdW force spectrum calculations46 in Figs. S10(a), S10(c), and S10(e) of the supplementary material, which show that the forces each amorphous substrate exerts on MoS2 are comparable.

In contrast, the spectral TC of MoS2 on crystalline substrates, which is shown as solid lines in Figs. 2(b)2(d), is quite different among the materials. [A direct comparison of these spectral TC curves is in Fig. S7(b) of the supplementary material.] We find qualitative similarities between the spectral TC for MoS2 on c-SiO2 and c-AlN, with both having a peak spectral TC at ∼3 THz; however, spectral TC contributions for MoS2 on c-AlN are much larger, with its peak ∼2× larger than for MoS2 on c-SiO2. For MoS2 on c-Al2O3 (sapphire), we see that contributions to TC remain high, even down to 0 THz, similar to suspended MoS2. This suggests that lower-frequency, longer-wavelength phonons are not as severely affected when MoS2 is on c-Al2O3 as the other substrates. This observation is in alignment with our vdW force spectrum calculations, which show that c-Al2O3 exerts the smallest force on MoS2 in the range of 0–4 THz compared to all the other substrates. A more detailed look at these calculations can be found in Sec. S5 of the supplementary material.

Interestingly, the spectral TCs of MoS2 on c-SiO2 and a-SiO2 [Fig. 2(b)] are found to be comparable, with total TCs less than ∼4 Wm−1 K−1 apart. The VDOS of both substrates over 0–8 THz (where acoustic modes of MoS2 are located) are similar too [Fig. 3(a)], suggesting that crystalline substrates with a large VDOS overlap with MoS2 can also provide sufficient interaction (i.e., vibrational modes with proper energies and momenta to scatter MoS2 phonons) to significantly reduce the TC of supported MoS2. In contrast, we find that the VDOS of c-AlN and c-Al2O3 are ∼2–5 times smaller than of a-AlN and a-Al2O3 in the 0–4 THz range, respectively [see Figs. 3(b) and 3(c)]. The spectral TCs of c-AlN and c-Al2O3, in Figs. 2(c) and 2(d), differ the most from their amorphous counterparts within this frequency range, which suggests that, with a smaller VDOS overlap, the remote phonon scattering from the c-AlN and c-Al2O3 substrates is limited compared to a-AlN and a-Al2O3. However, we note that, despite their similar VDOS in the 0–4 THz range [see Fig. S8(b) in the supplementary material], c-AlN and c-Al2O3 yield very different spectral TCs, emphasizing that VDOS overlap alone cannot be used to fully characterize the changes in the TC of supported MoS2.

Finally, we consider the effects of an h-BN interlayer on the spectral TC of MoS2. Recall that the total TC of MoS2 always increases when an h-BN interlayer is used, as shown in Fig. 1(g). The spectral TC reveals that the increase is entirely from contributions of phonon frequencies up to 3 THz (see Fig. S6 in the supplementary material). To some degree, the h-BN layer “blocks” the low-frequency interactions between MoS2 and the substrate, and the shape of the h-BN-supported spectral TC curves become closer to resembling the suspended MoS2 spectral TC curve in Fig. 2(a). Note that the contributions from out-of-plane motion, shown in Fig. S5 of the supplementary material, do not improve with an h-BN interlayer and that all TC gains are strictly from in-plane atomic motion (i.e., not flexural modes).

Since h-BN is a stiffer, high-TC material,55,67,68 it has, on average, 1.5–2 times fewer phonon modes than MoS2 over the acoustic mode frequencies of MoS2 (see Fig. S9 in the supplementary material), with this difference being as high as 7.7× between 0 and 2 THz. These observations suggest that remote phonon scattering with h- BN should be severely limited and the TC of MoS2 should not be strongly affected by it. This is easily tested by calculating the TC of MoS2 in a heterostructure with only h-BN, where we find that the TC of MoS2 only drops 17% to ∼101 Wm−1 K−1 [see Fig. S5(a) and Table S2 in the supplementary material]. In structures with substrates, we expect the small number of h-BN phonon modes and the additional distance between MoS2 and the substrate to limit remote phonon scattering in MoS2.

3. Length dependence of thermal conductivity

Next, we calculate the frequency dependence of the phonon MFP as λ(ω) = κ(ω)/G(ω),24 where the frequency-dependent conductance G(ω) describes ballistic thermal transport in MoS2. The results of these calculations are shown in Figs. 4(a) and 4(b). If we examine the MFP of suspended MoS2 in Fig. 4(a), we see that it drops rapidly from a peak of 3.3 μm at the Γ point (0 THz) to tens of nanometers by 6 THz. Previous theoretical studies, which calculated a suspended MoS2 TC similar to ours, have also shown MFPs up to several micrometers.69,70 Other simulations that calculated peak MFPs in the tens of nanometers only reported TCs between 20 and 40 Wm−1 K−1.71–73 Measurements of bulk in-plane TC also support claims of MFPs on the order of micrometers.74 Note that it is the very long MFPs, in combination with high phonon group velocities at low phonon frequencies (i.e., ≤2 THz), that dominate the majority of in-plane thermal transport for suspended MoS2, as shown in the spectral TC of Fig. 2(a).

FIG. 4.

The frequency-dependent mean free path MFP(ω) of MoS2 when supported by (a) amorphous substrates and (b) crystalline substrates at 300K. Both (a) and (b) also show MFP(ω) for suspended MoS2. (c) The room temperature length-dependent thermal conductivity κ(L) for suspended and supported MoS2. The red, green, and orange dotted (solid) lines show MoS2 when supported by amorphous (crystalline) substrates, and the blue line shows suspended MoS2.

FIG. 4.

The frequency-dependent mean free path MFP(ω) of MoS2 when supported by (a) amorphous substrates and (b) crystalline substrates at 300K. Both (a) and (b) also show MFP(ω) for suspended MoS2. (c) The room temperature length-dependent thermal conductivity κ(L) for suspended and supported MoS2. The red, green, and orange dotted (solid) lines show MoS2 when supported by amorphous (crystalline) substrates, and the blue line shows suspended MoS2.

Close modal

When MoS2 is supported by amorphous substrates, we find that the MFPs [also shown in Fig. 4(a)] are much smaller than in suspended MoS2, with maximum MFPs in the range of a few hundred nanometers and significant deviations from suspended MoS2 below 4 THz. This shows, again, that thermal transport and MFPs of long-wavelength phonons are severely disrupted by the amorphous substrates. In Fig. 4(b), we see that MoS2 on c-SiO2 (quartz) has a similar spectral MFP curve to a-SiO2 in Fig. 4(a), and the c-AlN- and c-Al2O3-supported MoS2 starts to see MFPs near 1 μm below 3 THz. Finally, if we consider the MFPs for systems with h-BN interlayers (Fig. S14 in the supplementary material), we find that the MFPs for systems of both amorphous and crystalline substrates increase significantly, with the MFP distribution more closely resembling that of suspended MoS2.

With the spectral TC and MFP, we use Eq. (3) to calculate the length-dependent TC, as shown in Fig. 4(c). We find that the TC of suspended MoS2 is length-dependent in samples shorter than ∼10 μm and only converges to its bulk limit at longer length scales. Recent experimental work corroborates our results showing that the TC of MoS2 continues to increase up to a suspended MoS2 membrane diameter of ∼13 μm.11 In comparison, the TC of graphene, which is over an order of magnitude larger than that of MoS2,75,76 has been calculated to converge at length scales also an order of magnitude longer, around ∼100 μm.77–79 A special symmetry selection rule forbids certain phonon–phonon scattering events in graphene,76 enabling flexural modes with very long MFPs. As MoS2 is three atoms thick, it does not follow this rule,80 resulting in diffusive thermal transport at a much shorter length than graphene.

The length-dependent trends for the TC of MoS2 on other substrates [also shown in Fig. 4(c)] follow expectations based on their shorter MFPs. MoS2 supported by crystalline substrates converges to its final value on the order of a few micrometers, and MoS2 supported by amorphous substrates converges to its final value on the order of several hundred nanometers. There are some important implications of the length-dependent TC in samples shorter than ∼100 nm. In this range, the TC of MoS2 is size-limited no matter the substrate [including cases with h-BN interlayers, as seen in Fig. S14(c) of the supplementary material]. This suggests that some applications, such as patterned nanoscale transistors,7,81 may not see a significant TC benefit when choosing one substrate over another. In such devices, heat sinking may be dominated by the TBC with thermal pathways through the substrate, gate insulator, and gate (see Sec. III C of this section). However, the MoS2-substrate vdW bond strength and force may still play a role because a substrate that is more tightly coupled to MoS2 could reduce the TC a non-negligible amount.9,46 In addition, substrate-induced strain effects82 could also influence the TC of MoS2, although only on the order of 5%.35 

Because the MoS2-substrate TBC is a key property for heat removal in systems based on 2D materials, we run AEMD simulations and use Eq. (4) to extract the TBC for all structures near room temperature. These TBCs are shown with the TCs in Fig. 5, and additional details can be found in Sec. S7 and Table S2 of the supplementary material. We find that the TBC for MoS2 on a-SiO2 and c-SiO2 is similar, in the range of 22–23 MWm−2 K−1. These TBCs are comparable to previous measurements and within the range of experimental error.6,10,11,83 Based on the similar spectral TC, VDOS, and vdW force spectrum between the two SiO2 substrates, this outcome is expected. Given that bulk c-SiO2 (quartz) has a TC of 10.7 Wm−1 K−1 parallel to the c-axis (or 6.7 Wm−1 K−1 perpendicular to it),47 whereas a-SiO2 has a TC of ∼1.4 Wm−1 K−1 (Ref. 84) near room temperature, c-SiO2 is clearly preferred if trying to maximize heat removal in such circumstances.

FIG. 5.

A scatterplot showing both the thermal boundary conductance (TBC) and thermal conductivity (κ) of MoS2 supported by crystalline and amorphous SiO2 (red), AlN (green), and Al2O3 (orange). The black arrow emphasizes that a larger TBC and thermal conductivity will enable better heat removal from MoS2. The green and orange arrows simply connect related calculations. The dashed vertical blue line denotes the thermal conductivity of suspended MoS2.

FIG. 5.

A scatterplot showing both the thermal boundary conductance (TBC) and thermal conductivity (κ) of MoS2 supported by crystalline and amorphous SiO2 (red), AlN (green), and Al2O3 (orange). The black arrow emphasizes that a larger TBC and thermal conductivity will enable better heat removal from MoS2. The green and orange arrows simply connect related calculations. The dashed vertical blue line denotes the thermal conductivity of suspended MoS2.

Close modal

Similarly, we find that the TBC of MoS2 on c-Al2O3 is comparable to that on a-Al2O3, with both in the range of ∼32–34 MWm−2 K−1. This means that the sapphire substrate yields a high TBC and TC for MoS2, making it superior for heat removal. This is especially true because the TC of bulk c-Al2O3 [∼34 Wm−1 K−1 (Ref. 85)] is much higher than for a-Al2O3 [∼1.6 Wm−1 K−1 (Ref. 53)] at room temperature. Understanding why the TBC of MoS2 on c-Al2O3 is comparable to a-Al2O3 requires further investigation; however, we see in Fig. S10(f) of the supplementary material that c-Al2O3 exerts a large vdW force on MoS2 in the 4–8 THz range, which is not observed for other substrates. We hypothesize that this is responsible for the larger TBC at the MoS2/c-Al2O3 interface. Note that this large vdW force also corresponds with a large VDOS overlap between c-Al2O3 and MoS2 seen in Fig. 3(c) over the same frequency range. The a-Al2O3 substrate has a similar VDOS overlap but not a similar vdW force spectrum, suggesting that a crystalline structure is the differentiator.

In contrast to the other substrates, the TBCs of MoS2 on a-AlN and c-AlN are significantly different at ∼29 and ∼14 MWm−2 K−1, respectively. However, because the TC of a-AlN [∼1.7 Wm−1 K−1 (Ref. 52)] is significantly lower than the TC of bulk c-AlN [>200 Wm−1 K−1 (Ref. 48)], it is unclear which substrate is best for heat removal, a scenario we will investigate in Sec. III C. While previous experimental work agrees with our TBC calculation of MoS2 on c-AlN,10 we have difficulty explaining its smaller value compared to our TBC of MoS2 on c-Al2O3. In Fig. 3(b), we also see a large VDOS overlap between MoS2 and c-AlN in the 4–8 THz range; however, we do not see corresponding vdW forces like for MoS2 on c- Al2O3 [see Figs. S10(d) and S10(f) in the supplementary material]. This indicates that, for crystalline substrates, the VDOS overlap alone is insufficient to estimate the relative TBCs of supported MoS2, and more detailed structural information must be considered. Future work leveraging techniques, such as interface conductance modal analysis (ICMA),15,19 can provide the necessary structural and frequency-domain insights to fully understand these TBCs. Even without this analysis, the lack of a trend in Fig. 5 still shows that the interactions that affect the TC and TBC of substrate-supported MoS2 are not the same (i.e., there is not necessarily a trade-off between the two properties).

Finally, we calculate the TBC of MoS2 with each h-BN-capped substrate. These structures are thermally more complex, and we must consider two TBCs: one between MoS2 and h-BN and another between h-BN and each substrate. Additional details about the thermal circuit and TBC extractions can be found in Sec. S7 of the supplementary material. We find the limiting TBC to be at the MoS2/h-BN interface, similar to experimental work on a-SiO2,86 with TBCs between 18 and 26 MWm−2 K−1 across all substrates. Interestingly, when the h-BN layer is introduced, the vdW forces on MoS2 from c-Al2O3 in the 4–8 THz range vanish [see Fig. S11(g) in the supplementary material], and at the same time, the TBC drops by nearly 2× (Table S2 in the supplementary material). This supports our hypothesis that the forces from c-Al2O3 in that range of frequencies is responsible for the higher TBC between MoS2 and c-Al2O3 (∼34.3 MWm−2 K−1 in Fig. 5 and Table S2 in the supplementary material).

To investigate the impact of the TC and TBC of supported MoS2 on the temperature rise in a transistor due to self-heating, we develop an analytical model for the thermal resistance of a simple back-gated MoS2 transistor,7 shown in Fig. 6(a), and validate it through finite element method simulations (see Sec. S9 in the supplementary material). Underneath the monolayer MoS2 channel is a 50 nm-thick insulating film on a silicon substrate, which, as in many laboratory measurements, serves as the back gate. For consistency with our MD simulations, the insulating films will be limited to crystalline and amorphous SiO2, AlN, and Al2O3. The MoS2 lies between and underneath the two contacts, which are 100 nm thick and 500 nm long [i.e., extended in the x-direction in Fig. 6(a)]. In this section, we consider a wide device (10 μm in the y-direction) to focus on the length dependence of the peak device temperature rather than on effects arising from narrow channel widths. (For narrow channel widths,87 the lateral heat spreading beyond MoS2 in the y-direction, which increases the thermal footprint of the MoS2 transistor, is non-negligible, and the TC of MoS2 may be reduced by the narrow channel.) The peak thermal resistance is

Rth=ΔTmaxP,
(5)

where P is the power assumed to be uniformly dissipated in the MoS2 channel,6 consistent with a transistor operating in the linear regime, and ΔTmax is the resulting maximum temperature rise, which occurs at the center of the channel. The dependence of Rth on the transistor geometry and material properties has previously been modeled as87,88

Rth=1gL[1(coshL2LH+gLHRconsinhL2LH)1],
(6)

where LH = (κWt/g)1/2 is the thermal healing length along the MoS2, i.e., the characteristic distance over which the temperature drops by 1/e from the contacts. In addition, g is the thermal conductance per unit length from the channel to the bottom of the Si substrate. Note that a substrate with a small TC like polyimide2 can significantly reduce g, but, for the high-TC Si substrate (κsub = 150 Wm−1 K−1) considered in this section, the substrate contribution to Rth is rather small, except for long, wide devices where the substrate can account for about half of Rth. The contact thermal resistance is Rcon, and L, W, t, κ are the MoS2 channel length, width, thickness, and TC, respectively. The dimensions are labeled in Fig. 6(a). The TC of MoS2 affects Rth through its influence on Rcon and LH, while the TBC of MoS2 impacts g. The detailed descriptions of these quantities are given in Sec. S9 of the supplementary material.

FIG. 6.

(a) The back-gated transistor geometry with relevant properties defined for the analytical model (not to scale). Note that this geometry is similar to semiconductor-on-insulator (SOI) transistors, but with the semiconductor at its atomically thin limit with monolayer MoS2. (b) The thermal resistance of a back-gated MoS2 transistor, plotted as a function of channel length for different insulator materials, all on a silicon wafer substrate as is often the case in simple proof-of-concept experimental transistors. Solid and dashed lines correspond to the crystalline and amorphous phases of the insulators, respectively, with appropriate MoS2 TCs and TBCs used from our calculations. (c) The same data, plotted with a linear-scale vertical axis, to highlight the thermal resistance variation for short-channel transistors. The gray, vertical dashed lines in (b) and (c) denote the approximate transition between long- and short-channel regimes (x-direction).

FIG. 6.

(a) The back-gated transistor geometry with relevant properties defined for the analytical model (not to scale). Note that this geometry is similar to semiconductor-on-insulator (SOI) transistors, but with the semiconductor at its atomically thin limit with monolayer MoS2. (b) The thermal resistance of a back-gated MoS2 transistor, plotted as a function of channel length for different insulator materials, all on a silicon wafer substrate as is often the case in simple proof-of-concept experimental transistors. Solid and dashed lines correspond to the crystalline and amorphous phases of the insulators, respectively, with appropriate MoS2 TCs and TBCs used from our calculations. (c) The same data, plotted with a linear-scale vertical axis, to highlight the thermal resistance variation for short-channel transistors. The gray, vertical dashed lines in (b) and (c) denote the approximate transition between long- and short-channel regimes (x-direction).

Close modal

In the long-channel limit (L ≫ 3LH), the effect of the contacts is negligible and Rth depends purely on the heat transfer down into the supporting insulator and the Si substrate below it. In this case, we have Rth ≈ 1/(gL), which for a thin underlying insulator can be further approximated as

Rth1WL×TBCMoS2ins,
(7)

where W is the channel width and TBCMoS2−ins is the TBC between MoS2 and the underlying insulator. (In this limit, when only TBC is considered, 3LH ≈ 150 nm for all material interfaces in this work.) Equation (7) shows that the peak temperature is directly influenced by the TBC of MoS2. If the thermal resistance contribution of the underlying insulator is not negligible, a better estimate is given by

RthTBCMoS2-ins1+TBCins-sub1+tins/κinsWL,
(8)

where TBCins-sub is the TBC between the insulator and the Si substrate (which can usually be neglected when compared to other thermal resistances), κins is the (vertical) insulator thermal conductivity, and tins is the insulator thickness.

In the short-channel limit (L < 3LH), heat flow from MoS2 directly into the insulator underneath is negligible, and heat removal is facilitated primarily by the contacts. In this limit, Rth ≈ Rcon/2, where the factor of two is the result of the two contacts providing heat sinking paths in parallel.89 If the thermal resistance contributions of the underlying insulator and the contact metals are negligible, Rth can be simplified as

Rth12Wtκ(TBCMoS2ins+TBCMoS2met),
(9)

where TBCMoS2−met is the TBC from MoS2 to the contact metals. Equation (9) depends only on the thermal properties of MoS2 and its interfaces and can be viewed as the intrinsic thermal resistance of the MoS2 channel, without contributions from external insulating layers. This also shows that the TC of MoS2 and total TBC to the materials on either side are equally important for short-channel devices. Note that TBCMoS2−met is typically on the same order of magnitude as TBCMoS2−ins and cannot be neglected. In this work, we assume TBCMoS2−met = 20 MWm−2 K−1, which is typical of the interface with Au, Ti, and Al.74,90,91

Equation (9) represents a good approximation for the transistor geometry and materials considered here, as long as the contacts are longer than ∼500 nm in the x-direction. In general, this expression will hold if the underlying insulator is thin (≲100 nm) and the contact metals act as good heat sinks. Contacts are good heat sinks if they are moderately long (≥500 nm), extend beyond the width of the device to connect to other devices, or join larger interconnects through vias. For shorter contacts, Rth depends more on the TBC of MoS2 than the TC for all channel lengths [see Figs. S18(b) and S18(d) in the supplementary material]. In addition, TBCMoS2−met has a smaller impact on Rth than TBCMoS2−ins because short metal contacts are of limited utility in cooling the device. In our simple model, the contacts are assumed to be exactly as wide as the channel and not connected to an external heat sink. They merely act as lateral heat spreaders, ultimately assisting heat flow into the substrate.

Figures 6(b) and 6(c) show the dependence of Rth on the transistor channel length L for different underlying insulator materials. Note the distinct long-channel regime, where Rth is inversely proportional to L (except for very long channels where the substrate thermal resistance, which is proportional to L1/2, becomes significant), and the short-channel regime, where Rth is relatively independent of L. The channel length that roughly separates these two regimes is 3LH ≈ 150 nm. For long channels, thanks to a higher TBC with MoS2 and a relatively high cross-plane thermal conductivity (18 Wm−1 K−1),92 the c-Al2O3 substrate achieves the lowest Rth. On the other hand, despite its high thermal conductivity (35 Wm−1 K−1),48,93 the c-AlN substrate results in a high Rth, exceeded only by a-SiO2 due to its low TC (1.4 Wm−1 K−1),84 in accordance with Eq. (8). At short channels, consistent with Eq. (9), the c-Al2O3 substrate once again yields the lowest Rth due to the high TC and TBC of MoS2 when supported by it, followed by the c-AlN substrate, which achieves the highest MoS2 TC. Note that the TCs listed for the crystalline insulators account for size constraints (i.e., not bulk) and anisotropy (e.g., for quartz and sapphire).

We note that we have only considered back-gated MoS2 transistors in this section, which is the simpler, more commonly encountered geometry in laboratory experiments. Top-gated transistors will be the subject of future work, as top gates have been found to play a more important role in heat spreading when the substrate TC is extremely low, such as is the case with polyimide.2 

We investigated the dependence of the TC of MoS2 on crystalline and amorphous SiO2, AlN, and Al2O3 substrates, showing that the TC of MoS2 is larger when supported by crystalline substrates than amorphous and that amorphous substrates yield approximately the same TC for MoS2 regardless of material. For any substrate, adding even a single h-BN interlayer improves the TC of MoS2. The degradation of MoS2 TC is primarily due to substrate interactions with long-wavelength, low-frequency MoS2 phonons (i.e., <2 THz), which a majority of suspended MoS2 contributions arise from. The in-plane phonon mean free paths for MoS2 can be very long (i.e., >1 μm), and the diffusive thermal transport regime is not reached until length scales greater than a micrometer. Additionally, at lengths below 100 nm, the TC of MoS2 becomes mostly limited by system size (i.e., MoS2 length and width) and not substrate interactions. Finally, we found the TBC of MoS2 to behave notably different for different crystalline substrates, i.e., the TBC with c-Al2O3 being over double that with c-AlN, despite yielding similar MoS2 TCs. This suggests that the mechanisms controlling the TC and TBC of MoS2 are different and that both could be optimized on a single substrate.

Finally, using an analytical model to determine the thermal resistance of a back-gated MoS2 transistor, we evaluated the impact of our calculated TBCs and TCs of MoS2 when supported by the insulating substrates considered in this study. Of these substrates, c-Al2O3 leads to the smallest temperature rise in an MoS2 transistor. MoS2 transistors on amorphous substrates will have larger temperature rises than on crystalline substrates as the TC of MoS2 is worse and the low substrate TCs contribute significantly to the total device resistance. For transistors with long channels (L > 3LH ≈ 150 nm) or short contacts (<500 nm), the TBC between MoS2 and the underlying substrate is the dominant mechanism for heat removal. For short-channel devices (L < 3LH ≈ 150 nm) with longer contacts (>500 nm), the TC of MoS2 and total TBC with the contacts and the substrate are equally important for heat removal. Overall, these TBC and TC results, as well as the analytical thermal model, provide important insights to evaluate heating effects during the operation of electronic and optical devices based on MoS2 and similar 2D materials.

See the supplementary material for details on Lennard-Jones parameters (Sec. S1), structure creation (Sec. S2), homogeneous nonequilibrium simulations (Sec. S3), vibrational density of states calculations (Sec. S4), van der Waals force spectrum calculations (Sec. S5), nonequilibrium simulations (Sec. S6), approach to equilibrium calculations (Sec. S7), final calculated thermal properties (Sec. S8), and analytic thermal model (Sec. S9).

Some of the computing for this project was performed on the Sherlock cluster at Stanford University. We would like to thank Stanford University and the Stanford Research Computing Center (SRCC) for providing computational resources and support that contributed to these results. This work was also partially supported by the Stanford SystemX Alliance and by ASCENT, one of the six centers in JUMP, a Semiconductor Research Corporation (SRC) program sponsored by DARPA. A.J.G. also acknowledges support from the Achievement Rewards for College Scientists (ARCS) Northern California Chapter.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
M. M. S.
Aly
,
M.
Gao
,
G.
Hills
,
C.-S.
Lee
,
G.
Pitner
,
M. M.
Shulaker
,
T. F.
Wu
,
M.
Asheghi
,
J.
Bokor
,
F.
Franchetti
 et al,
Computer
48
,
24
(
2015
).
2.
A.
Daus
,
S.
Vaziri
,
V.
Chen
,
Ç
Köroğlu
,
R. W.
Grady
,
C. S.
Bailey
,
H. R.
Lee
,
K.
Schauble
,
K.
Brenner
, and
E.
Pop
,
Nat. Electron.
4
,
495
501
(
2021
).
3.
X.
Li
,
L.
Tao
,
Z.
Chen
,
H.
Fang
,
X.
Li
,
X.
Wang
,
J.-B.
Xu
, and
H.
Zhu
,
Appl. Phys. Rev.
4
,
021306
(
2017
).
4.
C. D.
English
,
G.
Shine
,
V. E.
Dorgan
,
K. C.
Saraswat
, and
E.
Pop
,
Nano Lett.
16
,
3824
3830
(
2016
).
5.
C. J.
McClellan
,
E.
Yalon
,
K. K. H.
Smithe
,
S. V.
Suryavanshi
, and
E.
Pop
,
ACS Nano
15
,
1587
1596
(
2021
).
6.
E.
Yalon
,
C. J.
McClellan
,
K. K. H.
Smithe
,
M.
Muñoz Rojo
,
R. L.
Xu
,
S. V.
Suryavanshi
,
A. J.
Gabourie
,
C. M.
Neumann
,
F.
Xiong
,
A. B.
Farimani
 et al,
Nano Lett.
17
,
3429
3433
(
2017
).
7.
S.
Das
,
A.
Sebastian
,
E.
Pop
,
C. J.
McClellan
,
A. D.
Franklin
,
T.
Grasser
,
T.
Knobloch
,
Y.
Illarionov
,
A. V.
Penumatcha
,
J.
Appenzeller
 et al,
Nat. Electron.
4
,
786
799
(
2021
).
8.
I. M.
Datye
,
A. J.
Gabourie
,
C. D.
English
,
K. K. H.
Smithe
,
C. J.
McClellan
,
N. C.
Wang
, and
E.
Pop
,
2D Mater.
6
,
011004
(
2018
).
9.
A. J.
Gabourie
,
S. V.
Suryavanshi
,
A. B.
Farimani
, and
E.
Pop
,
2D Mater.
8
,
011001
(
2021
).
10.
E.
Yalon
,
Ö. B.
Aslan
,
K. K. H.
Smithe
,
C. J.
McClellan
,
S. V.
Suryavanshi
,
F.
Xiong
,
A.
Sood
,
C. M.
Neumann
,
X.
Xu
,
K. E.
Goodson
 et al,
ACS Appl. Mater. Interfaces
9
,
43013
43020
(
2017
).
11.
Y.
Yu
,
T.
Minhaj
,
L.
Huang
,
Y.
Yu
, and
L.
Cao
,
Phys. Rev. Appl.
13
,
034059
(
2020
).
12.
G.-H.
Lee
,
Y.-J.
Yu
,
X.
Cui
,
N.
Petrone
,
C.-H.
Lee
,
M. S.
Choi
,
D.-Y.
Lee
,
C.
Lee
,
W. J.
Yoo
,
K.
Watanabe
 et al,
ACS Nano
7
,
7931
7936
(
2013
).
13.
G.-H.
Lee
,
X.
Cui
,
Y. D.
Kim
,
G.
Arefe
,
X.
Zhang
,
C.-H.
Lee
,
F.
Ye
,
K.
Watanabe
,
T.
Taniguchi
,
P.
Kim
 et al,
ACS Nano
9
,
7019
7026
(
2015
).
14.
L.
Lindsay
,
C.
Hua
,
X. L.
Ruan
, and
S.
Lee
,
Mater. Today Phys.
7
,
106
120
(
2018
).
15.
H. R.
Seyf
,
K.
Gordiz
,
F.
DeAngelis
, and
A.
Henry
,
J. Appl. Phys.
125
,
081101
(
2019
).
16.
W.
Lv
and
A.
Henry
,
New J. Phys.
18
,
013028
(
2016
).
17.
B.
Mortazavi
,
M.
Pötschke
, and
G.
Cuniberti
,
Nanoscale
6
,
3344
3352
(
2014
).
18.
H.
Babaei
,
A. J. H.
McGaughey
, and
C. E.
Wilmer
,
Chem. Sci.
8
,
583
589
(
2017
).
19.
K.
Gordiz
and
A.
Henry
,
New J. Phys.
17
,
103002
(
2015
).
20.
C.
Monachon
,
L.
Weber
, and
C.
Dames
,
Annu. Rev. Mater. Res.
46
,
433
463
(
2016
).
21.
Z.
Fan
,
W.
Chen
,
V.
Vierimaa
, and
A.
Harju
,
Comput. Phys. Commun.
218
,
10
16
(
2017
).
22.
Z.
Fan
,
L. F. C.
Pereira
,
H.-Q.
Wang
,
J.-C.
Zheng
,
D.
Donadio
, and
A.
Harju
,
Phys. Rev. B
92
,
094301
(
2015
).
23.
Z.
Fan
,
L. F. C.
Pereira
,
P.
Hirvonen
,
M. M.
Ervasti
,
K. R.
Elder
,
D.
Donadio
,
T.
Ala-Nissila
, and
A.
Harju
,
Phys. Rev. B
95
,
144309
(
2017
).
24.
Z.
Fan
,
H.
Dong
,
A.
Harju
, and
T.
Ala-Nissila
,
Phys. Rev. B
99
,
064308
(
2019
).
25.
A. J.
Gabourie
,
Z.
Fan
,
T.
Ala-Nissila
, and
E.
Pop
,
Phys. Rev. B
103
,
205421
(
2021
).
26.
T.
Liang
,
S. R.
Phillpot
, and
S. B.
Sinnott
,
Phys. Rev. B
79
,
245110
(
2009
).
27.
T.
Liang
,
S. R.
Phillpot
, and
S. B.
Sinnott
,
Phys. Rev. B
85
,
199903
(
2012
).
28.
J. A.
Stewart
and
D. E.
Spearot
,
Modell. Simul. Mater. Sci. Eng.
21
,
045003
(
2013
).
29.
K.
Xu
,
A. J.
Gabourie
,
A.
Hashemi
,
Z.
Fan
,
N.
Wei
,
A.
Barati Farimani
,
H.-P.
Komsa
,
A. V.
Krasheninnikov
,
E.
Pop
, and
T.
Ala-Nissila
,
Phys. Rev. B
99
,
054303
(
2019
).
30.
J.
Tersoff
,
Phys. Rev. B
39
,
5566
5568
(
1989
).
31.
S.
Munetoh
,
T.
Motooka
,
K.
Moriguchi
, and
A.
Shintani
,
Comput. Mater. Sci.
39
,
334
339
(
2007
).
32.
P.
Vashishta
,
R. K.
Kalia
,
A.
Nakano
,
J. P.
Rino
, and
C. f. A. Computing
, and
Simulations
,
J. Appl. Phys.
109
,
033514
(
2011
).
33.
P.
Vashishta
,
R. K.
Kalia
,
A.
Nakano
, and
J. P.
Rino
,
J. Appl. Phys.
103
,
083504
(
2008
).
34.
C.
Sevik
,
A.
Kinaci
,
J. B.
Haskins
, and
T.
Çağın
,
Phys. Rev. B
84
,
085409
(
2011
).
35.
L.
Zhu
,
T.
Zhang
,
Z.
Sun
,
J.
Li
,
G.
Chen
, and
S. A.
Yang
,
Nanotechnology
26
,
465707
(
2015
).
36.
R. G.
Dickinson
and
L.
Pauling
,
J. Am. Chem. Soc.
45
,
1466
1471
(
1923
).
37.
R. S.
Pease
,
Nature
165
,
722
723
(
1950
).
38.
D. J.
Evans
,
Phys. Lett. A
91
,
457
460
(
1982
).
39.
K.
Sääskilahti
,
J.
Oksanen
,
S.
Volz
, and
J.
Tulkki
,
Phys. Rev. B
91
,
115426
(
2015
).
40.
K.
Sääskilahti
,
J.
Oksanen
,
J.
Tulkki
, and
S.
Volz
,
Phys. Rev. B
90
,
134312
(
2014
).
41.
Z.
Li
,
S.
Xiong
,
C.
Sievers
,
Y.
Hu
,
Z.
Fan
,
N.
Wei
,
H.
Bao
,
S.
Chen
,
D.
Donadio
, and
T.
Ala-Nissila
,
J. Chem. Phys.
151
,
234105
(
2019
).
42.
C.
Nyby
,
A.
Sood
,
P.
Zalden
,
A. J.
Gabourie
,
P.
Muscher
,
D.
Rhodes
,
E.
Mannebach
,
J.
Corbett
,
A.
Mehta
,
E.
Pop
 et al,
Adv. Funct. Mater.
30
,
2002282
(
2020
).
43.
Z.-Y.
Ong
and
E.
Pop
,
Phys. Rev. B
81
,
155408
(
2010
).
44.
S. V.
Suryavanshi
,
A. J.
Gabourie
,
A.
Barati Farimani
, and
E.
Pop
,
J. Appl. Phys.
126
,
055107
(
2019
).
45.
K.
Schauble
,
D.
Zakhidov
,
E.
Yalon
,
S.
Deshmukh
,
R. W.
Grady
,
K. A.
Cooley
,
C. J.
McClellan
,
S.
Vaziri
,
D.
Passarello
,
S. E.
Mohney
 et al,
ACS Nano
14
,
14798
14808
(
2020
).
46.
L.
Zhang
,
Y.
Zhong
,
X.
Qian
,
Q.
Song
,
J.
Zhou
,
L.
Li
,
L.
Guo
,
G.
Chen
, and
E. N.
Wang
,
ACS Appl. Mater. Interfaces
13
,
46055
46064
(
2021
).
47.
R. W.
Powell
,
C. Y.
Ho
, and
P. E.
Liley
,
Thermal Conductivity of Selected Materials
(
U.S. Department of Commerce, National Bureau of Standards; for sale by the Superintendent of Documents, U.S. Government Printing Office
,
Washington
,
DC
,
1966
).
48.
R. L.
Xu
,
M. M.
Rojo
,
S. M.
Islam
,
A.
Sood
,
B.
Vareskic
,
A.
Katre
,
N.
Mingo
,
K. E.
Goodson
,
H. G.
Xing
,
D.
Jena
 et al,
J. Appl. Phys.
126
,
185105
(
2019
).
49.
G. A.
Slack
,
R. A.
Tanzilli
,
R. O.
Pohl
, and
J. W.
Vandersande
,
J. Phys. Chem. Solids
48
,
641
647
(
1987
).
50.
T.
Yamane
,
N.
Nagai
,
S.
Katayama
, and
M.
Todoki
,
J. Appl. Phys.
91
,
9772
(
2002
).
51.
K. T.
Regner
,
D. P.
Sellan
,
Z.
Su
,
C. H.
Amon
,
A. J. H.
McGaughey
, and
J. A.
Malen
,
Nat. Commun.
4
,
1640
(
2013
).
52.
Y.
Zhao
,
C.
Zhu
,
S.
Wang
,
J. Z.
Tian
,
D. J.
Yang
,
C. K.
Chen
,
H.
Cheng
, and
P.
Hing
,
J. Appl. Phys.
96
,
4563
4568
(
2004
).
53.
I.
Stark
,
M.
Stordeur
, and
F.
Syrowatka
,
Thin Solid Films
226
,
185
190
(
1993
).
54.
X.
Yang
,
X.
Zheng
,
T.
Zhang
,
H.
Chen
, and
M.
Liu
,
Int. J. Heat Mass Transfer
170
,
121013
(
2021
).
55.
M.
Rahman
,
M.
Shahzadeh
, and
S.
Pisana
,
J. Appl. Phys.
126
,
205103
(
2019
).
56.
Z.
Zhang
,
S.
Hu
,
J.
Chen
, and
B.
Li
,
Nanotechnology
28
,
225704
(
2017
).
57.
S.-T.
Lin
,
M.
Blanco
, and
W. A.
Goddard
 III
,
J. Chem. Phys.
119
,
11792
11805
(
2003
).
58.
C. J.
Foss
and
Z.
Aksamija
,
2D Mater.
6
,
025019
(
2019
).
59.
G. C.
Correa
,
C. J.
Foss
, and
Z.
Aksamija
,
Nanotechnology
28
,
135402
(
2017
).
60.
C. J.
Foss
and
Z.
Aksamija
,
Nanotechnology
32
(40),
405206
(
2021
).
61.
T.
Ishibe
,
R.
Okuhata
,
T.
Kaneko
,
M.
Yoshiya
,
S.
Nakashima
,
A.
Ishida
, and
Y.
Nakamura
,
Commun. Phys.
4
,
153
(
2021
).
62.
C.
Hua
,
X.
Chen
,
N. K.
Ravichandran
, and
A. J.
Minnich
,
Phys. Rev. B
95
,
205423
(
2017
).
63.
H. R.
Seyf
and
A.
Henry
,
J. Appl. Phys.
120
,
025101
(
2016
).
64.
P. B.
Allen
,
J. L.
Feldman
,
J.
Fabian
, and
F.
Wooten
,
Philos. Mag. B
79
,
1715
1731
(
1999
).
65.
J.
Moon
,
B.
Latour
, and
A. J.
Minnich
,
Phys. Rev. B
97
,
024201
(
2018
).
66.
F.
DeAngelis
,
M. G.
Muraleedharan
,
J.
Moon
,
H. R.
Seyf
,
A. J.
Minnich
,
A. J. H.
McGaughey
, and
A.
Henry
,
Nanoscale Microscale Thermophys. Eng.
23
,
81
116
(
2019
).
67.
Q.
Cai
,
D.
Scullion
,
W.
Gan
,
A.
Falin
,
S.
Zhang
,
K.
Watanabe
,
T.
Taniguchi
,
Y.
Chen
,
E. J. G.
Santos
, and
L. H.
Li
,
Sci. Adv.
5
,
eaav0129
(
2019
).
68.
H.
Ying
,
A.
Moore
,
J.
Cui
,
Y.
Liu
,
D.
Li
,
S.
Han
,
Y.
Yao
,
Z.
Wang
,
L.
Wang
, and
S.
Chen
,
2D Mater.
7
,
015031
(
2020
).
69.
A. N.
Gandi
and
U.
Schwingenschlögl
,
Europhys. Lett.
113
,
36002
(
2016
).
70.
B.
Peng
,
H.
Zhang
,
H.
Shao
,
Y.
Xu
,
X.
Zhang
, and
H.
Zhu
,
Ann. Phys.
528
,
504
511
(
2016
).
71.
X.
Wei
,
Y.
Wang
,
Y.
Shen
,
G.
Xie
,
H.
Xiao
,
J.
Zhong
, and
G.
Zhang
,
Appl. Phys. Lett.
105
,
103902
(
2014
).
72.
D.
Saha
and
S.
Mahapatra
,
Physica E: Low Dimens. Syst. Nanostruct.
83
,
455
460
(
2016
).
73.
Y.
Cai
,
J.
Lan
,
G.
Zhang
, and
Y.-W.
Zhang
,
Phys. Rev. B
89
,
035438
(
2014
).
74.
J.
Liu
,
G.-M.
Choi
, and
D. G.
Cahill
,
J. Appl. Phys.
116
,
233107
(
2014
).
75.
E.
Pop
,
V.
Varshney
, and
A. K.
Roy
,
MRS Bull.
37
,
1273
1281
(
2012
).
76.
L.
Lindsay
,
D. A.
Broido
, and
N.
Mingo
,
Phys. Rev. B
82
,
115427
(
2010
).
77.
Y.
Kuang
,
L.
Lindsay
,
S.
Shi
,
X.
Wang
, and
B.
Huang
,
Int. J. Heat Mass Transfer
101
,
772
778
(
2016
).
78.
G.
Fugallo
,
A.
Cepellotti
,
L.
Paulatto
,
M.
Lazzeri
,
N.
Marzari
, and
F.
Mauri
,
Nano Lett.
14
,
6109
6114
(
2014
).
79.
G.
Barbarino
,
C.
Melis
, and
L.
Colombo
,
Phys. Rev. B
91
,
035416
(
2015
).
80.
X.
Gu
,
Y.
Wei
,
X.
Yin
,
B.
Li
, and
R.
Yang
,
Rev. Mod. Phys.
90
,
041002
(
2018
).
81.
C. D.
English
,
K. K. H.
Smithe
,
R. L.
Xu
, and
E.
Pop
, in
2016 IEEE International Electron Devices Meeting (IEDM)
(
IEEE
,
2016
), p.
5.6.1
.
82.
G. H.
Ahn
,
M.
Amani
,
H.
Rasool
,
D.-H.
Lien
,
J. P.
Mastandrea
,
J. W.
Ager
, III
,
M.
Dubey
,
D. C.
Chrzan
,
A. M.
Minor
, and
A.
Javey
,
Nat. Commun.
8
,
608
(
2017
).
83.
S.
Vaziri
,
E.
Yalon
,
M.
Muñoz Rojo
,
S. V.
Suryavanshi
,
H.
Zhang
,
C. J.
McClellan
,
C. S.
Bailey
,
K. K. H.
Smithe
,
A. J.
Gabourie
,
V.
Chen
 et al,
Sci. Adv.
5
,
eaax1325
(
2019
).
84.
D. G.
Cahill
and
T. H.
Allen
,
Appl. Phys. Lett.
65
,
309
311
(
1994
).
85.
D. G.
Cahill
,
S.-M.
Lee
, and
T. I.
Selinder
,
J. Appl. Phys.
83
,
5783
5786
(
1998
).
86.
Y.
Liu
,
Z.-Y.
Ong
,
J.
Wu
,
Y.
Zhao
,
K.
Watanabe
,
T.
Taniguchi
,
D.
Chi
,
G.
Zhang
,
J. T. L.
Thong
,
C.-W.
Qiu
 et al,
Sci. Rep.
7
,
43886
(
2017
).
87.
A. D.
Liao
,
J. Z.
Wu
,
X.
Wang
,
K.
Tahy
,
D.
Jena
,
H.
Dai
, and
E.
Pop
,
Phys. Rev. Lett.
106
,
256801
(
2011
).
88.
M. J.
Mleczko
,
R. L.
Xu
,
K.
Okabe
,
H.-H.
Kuo
,
I. R.
Fisher
,
H.-S. P.
Wong
,
Y.
Nishi
, and
E.
Pop
,
ACS Nano
10
,
7507
7514
(
2016
).
90.
M.
Goni
,
J.
Yang
, and
A. J.
Schmidt
,
Nano Res.
11
,
2173
2180
(
2018
).
91.
K. M.
Freedy
,
D. H.
Olson
,
P. E.
Hopkins
, and
S. J.
McDonnell
,
Phys. Rev. Mater.
3
,
104001
(
2019
).
92.
B.
Dongre
,
J.
Carrete
,
N.
Mingo
, and
G. K. H.
Madsen
,
MRS Commun.
8
,
1119
1123
(
2018
).
93.
C.
Duquenne
,
M.-P.
Besland
,
P. Y.
Tessier
,
E.
Gautron
,
Y.
Scudeller
, and
D.
Averty
,
J. Phys. D: Appl. Phys.
45
,
015301
(
2012
).

Supplementary Material