Phonon Thermal Transport in Silicene/Graphene Heterobilayer Nanostructures: Effect of Interlayer Interactions.

Heterostructuring, as a promising route to optimize the physical properties of 2D materials, has attracted great attention from the academic community. In this paper, we investigated the room-temperature in-plane and cross-plane phonon thermal transport in silicene/graphene van der Waals (vdW) heterostructures using molecular dynamics simulations. Our simulation results demonstrated that heat current along the graphene layer is remarkably larger than that along the silicene layer, which suggests that graphene dominates the thermal transport in silicene/graphene heterostructures. The in-plane phonon thermal conductivity of the silicene/graphene heterostructures could be a compromise between monolayer graphene and monolayer silicene. Heterostructuring can remarkably reduce the in-plane thermal conductivity of the graphene layer but increase the in-plane thermal conductivity of the silicene layer in heterobilayers compared with the freestanding monolayer counterparts because of their different structures. We also simulated the interlayer interaction strength effect on the in-plane phonon thermal conductivity and cross-plane interfacial thermal resistance of silicene/graphene heterostructures. Total in-plane phonon thermal conductivity and interfacial thermal resistance both decrease with the increase in the interlayer interaction strength in the silicene/graphene heterobilayers. In addition, the calculated interfacial thermal resistance shows the effect of the thermal transport direction across the interface. This study provides a useful reference for the thermal management regulation of 2D vdW heterostructures.

Introduction

Graphene has attracted great attention from the scientific community because of its novel physical properties, such as high mechanical strength, superior thermal conductivity, and ultrahigh carrier mobility, which make graphene a promising candidate for a wide variety of revolutionary applications in micro-electronics.[1−3] However, its zero band gap greatly impedes its application in the semiconductor industry. Therefore, graphene is inevitably hybridized with other 2D materials when used to build field-effect transistors. Silicene, a hexagonal monolayer of silicon with a buckled structure, is of interest due to its complete compatibility with modern silicon technologies. Silicene has been fabricated using epitaxial growth on the surfaces of metal substrate and 2D materials.[4] The first silicene-based field-effect transistors to work at room temperature (RT) have been made recently.[5] Thus, silicene has promising application prospects in the future because of its excellent electrical properties. However, the obstacle that hinders the application of silicene is its low thermal conductivity and mechanical strength.[6,7] This may seriously limit the stability and performance of electronic devices. Graphene has an exceptional combination of electrical properties, high thermal conductivity, and mechanical strength; thus, these two monolayer materials can be stacked on top of each other to form an artificial heterostructure, which may tap the desirable properties of both monolayers. The atoms in different layers interact through weak van der Waals (vdW) forces. Recent studies showed that silicene on a graphene layer or between two graphene layers can exist stably at RT because the graphene layer in the heterostructure can provide mechanical support for the silicene, and the weak interlayer interaction maintains the excellent electronic properties of silicene and graphene.[8−10] Therefore, the silicene/graphene vdW heterostructure can be highly promising for applications in field-effect transistors, which would further require proper thermal management.

Many previous studies have reported the thermal conductivity of freestanding graphene and silicene.[11−17] Theoretical and experimental studies have shown that the extremely high thermal conductivity of freestanding graphene is in the range of 2000–5000 W/mK near RT.[2,18] Accurate first-principle calculations and machine learning potentials are used to study the mechanical properties and thermal conductivity of silicene nanostructures.[19,20] The RT thermal conductivity of silicene (∼10–50 W/mK[6,15−17]) is considerably smaller than that of graphene and bulk silicon, indicating that silicene has a greater advantage in thermoelectric materials. Recently, in-plane and cross-plane thermal conductivity of heterostructures were more under discussion because of their importance in the thermal management of nanodevices. In addition, high-intensity field-effect transistors are likely to create larger heat range. Thus, it is significant to choose the type of thermal interface material in such a system.

Interfacial thermal resistance (ITR) describes the finite transmission resistance when incident phonons propagate through the interface, which is also known as thermal boundary resistance at the interface in the heterojunction. The ITR is attributed to the differences in lattice properties and scattering by the interface, which is a significant factor in heat management in many advanced applications. Although graphene is a 2D material with high thermal conductivity, the ITR still plays a critical role in dictating the overall heat transport of the graphene hybrid systems.[21] Liu et al.[22] investigated the interfacial thermal conductance of a silicene/graphene heterobilayer using molecular dynamics (MD) simulations and found that as the interface coupling strength increases, the enhanced coupling between the phonons of graphene and silicene facilitates the heat transfer between graphene and silicene and thus brings about the decrease in the ITR. Recently, Han et al.[23] calculated in-plane and out-of-plane thermal conduction of the graphene/C3N heterolayer and evaluated the effects of system size, temperature, and interlayer coupling strength on the in-plane and cross-plane thermal conductivity. The interlayer weak interaction also plays an important role in the ITR between graphene and other 2D materials, such as graphene/hBN, graphene/stanene, graphene/MoS2, and graphene/phosphorene vdW heterostructures.[24−29] Particularly, for graphene hybrid monolayer and heterobilayer, the thermal conductance across the interface increases with temperature.[22,23] Surface engineering approaches, including atomic intercalation[30] and chemical functionalization,[22,31] are usually applied to tune interfacial thermal conductance. However, the effect of the interface on the in-plane thermal transport behavior of silicene/graphene vdW heterostructures is less studied in the past. In the present paper, we deposited the zigzag silicene nanoribbon (z-SR) on the zigzag graphene nanoribbon (z-GR) to model heterostructure z-SR/z-GR (Figure ). We calculated the phonon thermal transport of silicene nanoribbon (SR) layer and graphene nanoribbon (GR) layer in the vdW heterostructure, as well as those of the freestanding monolayer SR and monolayer GR, using MD simulations. The in-plane phonon thermal conductivity and cross-plane ITR of the SR/GR hetrobilayers were discussed. This study is expected to provide theoretical guidance for the regulation of thermal transport in heterostructure devices.

Figure 1

Schematic representation of a z-SR/z-GR heterostructure: (a) top view; (b) front view; and (c) aerial view. Yellow and gray balls represent the silicon atoms of the silicene layer and the carbon atoms of the graphene layer, respectively. The armchair edge is shown in the longitudinal direction.

Results and Discussion

In-Plane Thermal Conductivity of GR and SR Monolayers

We first calculated the values for freestanding GR and SR monolayers to benchmark our calculations with the literature. We test the width dependence of RT in-plane phonon thermal conductivity (κ) for z-SR and z-GR under periodic boundary conditions and free boundary conditions applied in the width direction, respectively. From Figure , we found that the boundary conditions used in the MD simulations largely affect the in-plane thermal conductivity. When the nanoribbon width increases from 2 to 12 nm, the calculated RT in-plane κ of SR and GR both increases for free boundary condition (FBC) applied in the width direction, whereas the calculated RT in-plane κ of SR and GR appears to be saturated for periodic boundary condition (PBC) in the width direction at the width greater than 4 nm. To prevent the finite-size effect in the one direction and thus to reduce the computing burden, we used PBC in the width direction for the calculations later in this paper.

Figure 2

Width dependence of RT in-plane κ for (a) z-SR and (b) z-GR under PBC and FBC in the width direction, respectively. Here, the length (L) between the heat source and the heat sink in nanoribbon is about 40 nm.

Figure a,b displays the dependence of in-plane κ on the length (L) between the heat source and sink for freestanding z-GR and z-SR with the width of about 6 nm, respectively. Given that the L of the studied monolayer z-GR increases from 20 to 1200 nm, additional long-wavelength modes are included, and the calculated κ increases from ∼859.4 ± 12.5 to ∼2485.1 ± 35.9 W/mK at 300 K, consistent with the previous MD result (∼3500 W/mK for L = 3000 nm).[32] The calculated RT κ value of the studied GR is much lower than that of the infinite graphene (above 3000 W/mK),[33] which is due to the remarkable phonon boundary scattering induced by the limited size of atomic models for ballistic transport. Indeed, the size of the studied GR is considerably smaller than the effective phonon mean-free-path (MFP) of pristine graphene (about 775 nm at 300 K[14]); thus, it is not possible to incorporate the contribution of phonon modes of which MFP is much larger than the simulated domain size. Nonetheless, the calculated RT κ of freestanding z-SR increases with the increase in L from 5 to 130 nm and becomes saturated when L is larger than 20 nm. This result suggests that the effective MFP of phonon (∼18 nm[15]) of silicene is comparable to the length between the heat source and sink considered here.

Figure 3

Length (L) dependence of RT in-plane κ of (a) z-GR and (b) z-SR monolayers, and the inverse thermal conductivity (1/κ) of (c) z-GR and (d) z-SR nanoribbons vs its inverse length (1/L) at the ballistic-to-diffusive regime. z-SR width = 6.03 nm, thickness = 0.42 nm; z-GR width = 6.03 nm, thickness = 0.34 nm.

It is still challenging to conduct MD simulations for samples with tens of the micrometer scale. Generally, the thermal conductivity of the infinite-length nanoribbon (κ∞) and the effective MFP (Λ) were estimated by linearly fitting the data using the following equation[15,34]

It is noted that the length of the system used for fitting should be much larger than Λ in order to avoid an unphysical prediction.[15] As shown in Figure c,d, the inverse thermal conductivity (1/κ) of the nanoribbon is linear with its inverse length (1/L) at the ballistic-to-diffusive regime, which is confirmed in previous studies.[17,35] Here, at the RT κ∞ of silicene, the monolayer is estimated to be ∼14 W/mK by linear fitting with the least square method, in good agreement with the previous simulation (∼12 W/mK) reported by Zhang.[15] More accurate results are recently reported using first-principles methods (20–30 W/mK;[36] 15–30 W/mK[37]) and machine learning potentials (33.7 ± 0.6 W/mK;[38] 32.4 ± 2.9 W/mK[39]). The underestimation of the phonon thermal conductivity by the SW potential may be due to a significant underestimation of the group velocities of the phonons around the Γ point by the SW potential as compared to the DFT results.[38] Moreover, our estimated κ∞ of monolayer GR at 300 K is about 2940 W/mK, which is also close to previous theoretical and experimental results (2000–5000 W/mK[2,18,32,40,41]) for suspended monolayer graphene.

Next, we investigated the in-plane thermal transport in the z-SR/z-GR heterostructure. The calculated κ of the z-GR layer in the z-SR/z-GR heterostructure with the same length of 40 nm is about 849.2 ± 5.2 W/mK at 300 K, which is smaller than the κ value of the same-sized freestanding z-GR (1057.6 ± 12.3 W/mK at 300 K). However, the calculated κ of the z-SR layer in z-SR/z-GR heterostructure is 23.6 ± 0.6 W/mK at 300 K, which is larger than that of the freestanding z-SR layer (12.9 ± 0.3 W/mK at 300 K). Although the thermal transport through the monolayer silicene is inefficient, the supported GR layer efficiently diffuses the heat in the heterostructure. For example, under T = 300 K and ΔT = 60 K between the heat source and the heat sink, the major amount of heat energy in the heterostructure is transferred along the GR layer, with the heat flux along the GR layer being 20-fold higher than that along the SR layer (see Figure S1 in Supporting Information).

In-Plane Thermal Conductivity of SR/GR Heterostructures

Compared with monolayer z-GR with L = 40 nm (κ = 1057.6 ± 12.3 W/mK at 300 K), the apparent reduction in the κ of the z-SR/z-GR heterostructure (by ∼70%) is because of the following two reasons. As shown in Figure , the κ∞ of monolayer silicene is 2 orders of magnitude lower than that of monolayer graphene; therefore, the GR layer is the dominant thermal transport path in the heterostructure. In addition, MD simulation-predicted κ of 2D nanomaterial depends on the thickness (τ) for the nanosheet/nanofilm with the same width (i.e., κ ∝ 1/τ). In comparison, the κ of the z-SR/z-GR heterobilayer (∼328.2 ± 2.2 W/mK) is still about 25 times the κ of freestanding monolayer z-SR with L = 40 nm at 300 K. This result shows that the introduction of the graphene support can remarkably enhance the heat dissipation of silicene-based electronic devices.

Because the temperature difference between the SR and GR layers is much smaller under the steady-state, interlayer heat transfer can be negligible compared to the thermal energy passing through each layer. The temperature gradients have minor differences among the SR layer, GR layer, and SR/GR bilayer, as shown in Figure S2. Since the heat flux of the bilayer heterostructure is equal to the summation of the heat flux of the silicene monolayer and graphene, the in-plane κ of the bilayer heterostructure can be approximated as follows[42,43]

It is worth noting that the composite rule of the mixture, as shown in eq , is only applicable when the two monolayers have strong distinct in-plane κ values. Here, the thickness of the bilayer, which is the sum of the thickness of the two monolayers, is about 2.24 times the thickness of the GR layer, that is, τbilayer = 2.24 τGR. Therefore, the total κ of the bilayer is expected to compromise between the thermal conductivity of the two monolayers.

The phonon density of states (PDOS) may reveal the variations in the internal vibration mode of a system under the impact of external factors. We calculated the PDOS using the Fourier transformation of velocity auto-correlation function in MD simulations. Figure shows the total value and components of the PDOS of the SR and GR layers in the SR/GR heterostructure. The phonon frequencies of the GR layer extend up to 50 THz, whereas those in the SR layer range from 0 THz to ∼15 THz. The large mismatch in PDOS leads to phonon Umklapp scattering,[44] which markedly lowers the total phonon thermal conductivity of the SR/GR heterostructure compared with that of GR.

Figure 4

PDOS of z-SR and z-GR layers in the z-SR/z-GR heterostructure. Subscripts x, y, and z indicate the PDOS component along the coordinate direction, and the subscript all shows the sum of components.

Interfacial Interaction Effects on Phonon Thermal Transport in SR/GR Heterostructures

Next, we investigated the effect of the interlayer interaction strength on the in-plane and cross-plane thermal conductivity of SR/GR heterostructures. The non-bonding vdW interaction is described by U(r) = χV(r), where V(r) is the standard Lennard-Jones (LJ) potential and χ is a scaling factor used to adjust the strength of the interaction. Figure shows the in-plane κ of the z-SR and z-GR layers in the SR/GR heterostructure as a function of χ, together with the total in-plane κ of SR/GR heterostructures. Interestingly, our modeled z-GR in the heterostructure with 40 nm length and about 6 nm width shows about 20% lower in-plane phonon κ at RT than that of similar-sized freestanding z-GR, whereas the room-temperature in-plane κ of z-SR in the heterostructure increases by about 83% compared with that of similar-sized freestanding z-SR (Figures and 5). Obviously, the effect of interfacial interaction on buckled SR differs from that of flat GR. As χ increases from 1 to 30, the in-plane κ of the GR layer decreases monotonically, whereas the in-plane κ of the SR layer increases initially and then reduces slowly after χ = 10 in the SR/GR heterostructures. The special behavior of silicene can be attributed to the initial buckling structure. The introduction of the interlayer interaction, which is equivalent to the strain induced in the in-plane direction, will affect the structures of SR and GR. For the z-SR/z-GR heterostructure, the interfacial distance between bilayers decreases with the increase in χ (Figure ). Thus, the enhanced vdW interaction at larger χ may lower the in-plane thermal conductivity of the graphene layer because of the suppressed contribution of the phonon flexural mode (ZA).[45] However, at small χ, the SR structure becomes less buckled because of bond rotation, which leads to increased in-plane stiffness and thus an increase in the thermal conductivity of the SR layer. Moreover, at large χ, Si–Si bonds are severely stretched and distorted, which decreases the in-plane stiffness and thus thermal conductivity.[46] Because the RT thermal conductivity of the silicene layer is 2 orders of magnitude smaller than that of the graphene layer, the in-plane κ of bilayer follows mostly the trend of the graphene monolayer with a variation of the χ compared to the silicene monolayer. We also found that the in-plane κ of the heterolayers predicted theoretically using eq is in good agreement with the results achieved by MD calculations (Figure c).

Figure 5

Effect of interfacial interaction strength (χ) on the in-plane κ of (a) z-GR layer and (b) z-SR layer in the heterostructure as compared to the (c) z-SR/z-GR bilayer. The thickness of the bilayer comes from the addition of the thickness of two monolayers. The length and width of the heterostructures are approximately 40 and 6 nm, respectively.

Figure 6

Dependence of interlayer distance (d) on interfacial interaction strength (χ) in the z-SR/z-GR heterostructure.

ITR (R) or interfacial thermal conductance (G) describes the cross-plane thermal transport across the interface of the heterostructure. As shown in Figure , for χ = 1 at 300 K, the ITR R = 1/G = (13.1 ± 0.7) × 10–2 m2 K/MW when the heat flow is from the GR layer to the SR layer and R = (10.4 ± 0.7) × 10–2 m2 K/MW when the heat flow is from the SR layer to the GR layer, indicating a significant effect of the thermal transport direction on the interfacial thermal conductance of the bilayer.[26] We note that the calculated R also depends on the strength of the vdW coupling between the layers, consistent with previous studies of graphene-based vdW heterostructured interfaces.[23,25−27] The ITR of silicene/graphene bilayer decreases with the increase in χ, and it shows an insignificant dependence on the heat flow direction at larger interlayer coupling strength.

Figure 7

ITR (R) in the z-SR/z-GR heterostructures with different interfacial interaction strength (χ).

Conclusions

MD simulations were performed to investigate the in-plane and cross-plane phonon thermal transport in silicene/graphene heterostructures. The total in-plane phonon thermal conductivity of the heterostructure is expected to compromise between the in-plane phonon thermal conductivity of two monolayers. The graphene layer is the dominant thermal transport path in the heterostructure; thus, the introduction of graphene support can remarkably enhance the heat dissipation of silicene/graphene heterostructure devices. We also simulated the effects of interlayer interaction strength on the in-plane and cross-plane phonon thermal conductivity in silicene/graphene heterostructures. When the interlayer interaction increases from 1 to 30 times, the in-plane phonon thermal conductivity of GR is remarkably decreased, but the in-plane phonon thermal conductivity of SR increases initially and then decreases slightly. However, the total in-plane phonon thermal conductivity of the heterostructure decreased with the increase in interlayer interaction strength. In addition, the calculated ITR depends not only on the strength of the vdW coupling between the layers but also on the thermal transport direction across the interface. Our study is helpful to understand the phonon thermal transport in graphene/silicene heterostructures and provides a useful reference for thermal regulation in 2D vdW heterostructures.

Computational Methods

The MD simulations were performed using the large-scale atomic/molecular massively parallel simulator (LAMMPS) package.[47] The Tersoff potential parameterized by Lindsay and Broido[32] was used to model the interactions between C–C bonds in the graphene layer, which can accurately describe phonon dispersion of graphene and was previously widely adopted to study the thermal and mechanical properties of graphene.[18,24,48] The Si–Si interactions in the silicene layer were described using the optimized Stillinger–Weber (SW2) potential,[15] which can provide a good description of the buckling atomic configuration of silicene.[49] The newly parameterized Stillinger–Weber (SW2) potential can well reproduce the mechanical and thermal properties of silicene nanostructures,[15,31,49] consistent with the first-principles results and recently developed machine learning interatomic potentials.[16,36−39,50] The non-bonding atomic interactions (i.e., vdW interactions) between the silicene and graphene layers were described by the standard LJ potentialwhere r represents the interatomic distance, ε is the energy parameter (ε = 0.008909 eV), and σ is the distance parameter (σ = 0.3326 nm). These parameters were previously used to model the thermal boundary conductance between SiO2 and carbon nanotube[51] and the mechanical properties of graphene/silicene/graphene heterostructures.[10] The cutoff distance of the LJ potential was set as 1.0 nm, which is about threefold of σ.

Structure Relaxation

The velocity–Verlet algorithm was used to integrate Newton’s equations of motion with a time step of 0.5 fs. The silicon atoms in SR were initially displaced in a low-buckled 2D structure with a buckling distance of 0.042 nm. The average interlayer distance d was initially set to 0.36 nm, according to previous first-principles calculations.[22,52] The lattice mismatch between SR and GR in the heterostructure was less than 5%. We apply the periodic boundary condition along the in-plane x and y directions to eliminate the edge effect. The model structures were equilibrated to obtain a steady thermal state with a successive NPT (constant particle number, constant pressure, and constant temperature) ensemble for 0.5 ns and NVT (constant particle number, constant volume, and constant temperature) ensemble for 0.5 ns by the Nosé–Hoover approach with a background temperature of 300 K. Afterward, we relaxed the system with the NVE (constant particles, volume, and no thermostat) ensemble for another 5 ns.

In-Plane Phonon Thermal Conductivity Calculations

Following the equilibration, the non-equilibration MD method was employed for the simulations of the in-plane phonon thermal conductivity. 0.4 nm long regions at both ends in the longitudinal direction were set as rigid walls. 5 nm long regions next to the rigid walls were set as a heat source and sink, respectively, to establish a temperature gradient along the longitudinal direction. The temperature gradient was realized using a Langevin thermostat to maintain the temperatures of the heat sink and heat source at 270 and 330 K, respectively. When non-equilibration MD simulation runs for long enough time (5–10 ns, depending on the system size), the non-equilibrium steady-state temperature distribution in the system can be obtained. The sum of added energy and subtracted energy is equal to zero, and thus, total energy is conserved. The heat flux J is computed as the energy transported per unit time across the unit area, which can be recorded by the energy (E) injection/extraction rate in the heat source/sink (see Figure S1 in Supporting Information)where dE/dt is the thermal power across a cross-sectional area A.

For the 1D nanoribbon, the thermal conductivity κ of a finite system with length L is usually calculated by Fourier’s law of heat conductionwhere J is the heat flux in the longitudinal direction, and dT/dx is the temperature gradient. In this scheme, dT/dx is obtained from the slope of the so-called linear region of the temperature profile (see Figure S2 in Supporting Information), ignoring the nonlinear parts of the temperature profile near the thermal baths. This practice assumes that transport is diffusive, that is, in accordance with Fourier’s law that predicts a linear temperature profile at steady-state conditions. However, thermal transport at the nanoscale, especially in materials with high thermal conductivity such as graphene, is almost ballistic with a length-dependent κ(L).[53] Recent study[54] suggests that the temperature difference between the hot and cold thermostats should not exclude any local nonlinear regions of the temperature profile in the ballistic and ballistic-to-diffusive regimes, and one should directly calculate the thermal conductance from the temperature difference (ΔT) between the hot and cold thermostats and convert it into the thermal conductivity by multiplying it with the system length (L), that is

The test study shows the length dependence of κ predicted by eqs and 6, respectively, shows a similar trend, and is consistent with each other near diffusive regimes (see Figure S3 in the Supporting Information). Considering our simulation systems with tens of the nanometer length; here, we used eq to calculate κ of finite-sized GR and SR monolayers and its heterostructures.

Interfacial Thermal Resistance Calculations

With MD simulations, the thermal relaxation method is performed to calculate the ITR.[22,23,28,31] The thermal relaxation simulations mimic the experimental pump–probe approach.[55] In the experiment, while the pump laser applies a transient heat pulse to the sample, the probe laser detects the evolution of temperature. Compared to the traditional non-equilibration MD approach, thermal relaxation simulations could not only get accurate results but also save calculation cost.[56]

After the bilayer is equilibrated at 300 K, one of the layers (e.g., GR layer) is rapidly heated to a higher temperature (500 K) by velocity rescaling for 50 fs. Then, the thermostats are removed, and the layers are relaxed to the thermal equilibrium. During the thermal relaxation process, the heat is transferred from GR to SR in the form of kinetic energy and reaches a new equilibrium. Figure S4 in Supporting Information shows temperature changes of GR (TGR) and SR (TSR) and energy evolution of GR (EGR) after introducing the 50-fs thermal impulse. The energy variation of the GR layer with time can be described by[28]where Ainterfacial is the interfacial area, TGR and TSR represent the temperature of GR and SR layer, respectively, and ∂E/∂t denotes the variation of GR energy with time. Then, the ITR R can be calculated fromwhere E0 is the initial energy of the GR layer. The linear relationship betweenE and is shown in Figure S5 in the Supporting Information. It should be noted that the initial data (100 ps) are excluded from the linear fitting process because of the strong energy disturbance from thermal impulse to the system at the beginning of the thermal relaxing process. We performed six independent simulations to suppress the errors.