Lateral and flexural thermal transport in stanene/2D-SiC van der Waals heterostructure

After the invention of graphene, other group-IV two-dimensional (2D) elements, including silicene [1], germanene and stanene, have also been experimentally synthesized [2, 3] in recent years. Among these 2D elements, stanene [4] has fascinated as an attractive element for applications in field effect transistors [5], quantum anomalous Hall insulators [6], topology-based electronics as well as magnetic and optical devices [7]. However, the zero band gap without spin orbiting coupling [8] and low thermal conductivity (∼3–11.6 ${\text{ }}W{\text{ }}{m^{ - {\text{ }}1}}{\text{ }}{K^{ - {\text{ }}1}}$) [9–12] impede its application in various nanoelectronic and spintronic devices. In contrast, another promising nanomaterial, monolayer SiC, has proven its prominence for not only introducing a wide bandgap of ∼2.52 eV but also upholding the robust structure like graphene [13, 14]. Further, it holds a large exciton binding energy of ∼2.0 eV [15, 16], larger in-plane stiffness [17], Bose Einestein condensate effect [18] and comparatively high thermal conductivity (∼313 ${\text{ }}W{\text{ }}{m^{ - {\text{ }}1}}{\text{ }}{K^{ - {\text{ }}1}}$) [19–21]. Being an ultrathin semiconductor and possessor of these unique properties, 2D-SiC could be a unique material complementary to stanene which manifests a zero bandgap and low thermal conductivity at its pristine form.

On the other hand, van der Waals heterostructures (vdWHs) consisting of two dissimilar 2D materials have recently attracted significant interest because of their exceptional electrical, thermal, and optical properties. In particular, extraordinary transportation of charge carriers [22], bandgap engineering [23], semiconductor band alignment [24], and new optical absorption [25] have been realized from the vdWHs of distinct 2D materials. Moreover, the electrical and optical properties of such vdWHs can be tuned in a number of ways, such as twisted angles of component layers [26], material selections in stacking [27], interlayer distance [28], imposing an external electric field and strain [29] and so on. The wide range of unique properties and their tunability have stimulated the application of these vdWHs in flexible field effect transistors [30], tunneling transistors [31], highly sensitive photodetectors [32], photovoltaics [33], light-emitting diodes [34], multilevel inverters [35] and magnetic information storage devices [36–38], for instance.

In this regard, a new vdWH consisting of stanene and 2D-SiC, has recently been proposed [28]. Ferdous et al investigated the different electronic properties of this hetero-bilayer by varying the stacking patterns and observed a wide bandgap of ∼160 meV at the K point. In fact, the bandgap which was explored in their work is the widest among the previously studied hetero-bilayers associated with stanene [39–41]. In addition, the high binding energy (>223 meV per Sn atom), extraordinary electronic stability and insignificant lattice mismatching (∼1%) between stanene and 2D-SiC could make this structure promising for experimental realization as well as practical applications. It has also been observed that 2D-SiC preserves numerous electrical properties, including a large gap 2D quantum spin Hall state, superconductivity and high carrier mobility of stanene in an excessive amount. These exceptional and extraordinary electronic properties of this hetero-bilayer would encourage its application in high-speed tunneling transistors and spintronic devices. Moreover, the bandgap of a stanene/2D-SiC heterostructure can be modulated by applying biaxial tensile strain as well as by varying the interlayer distance between 2D-SiC and stanene, making it a prominent candidate for flexible electronics and optoelectronic devices, which include photodetectors, photovoltaics and light-emitting devices.

However, along with the extraordinary electronic properties, it is necessary to ensure decent thermal properties of a nanomaterial for its use in nanoelectronics. In fact, thermal management has become a major issue in nanoelectronic and optoelectronic devices due to the ongoing trend of shrinking and high-power densification. As phonons are the main heat energy carrier in 2D materials, the density of phonons will be increased with an increasing current, which in turn leads to a strong electron–phonon scattering [42, 43]. Consequently, local hot spots are engendered within the nanodevices and if this excess thermal energy does not disperse proficiently, it could create the foundation of defects, degradation in performance or even destruction of the devices. Though the thermal conductivity of monolayer stanene and 2D-SiC has already been studied, to the best of our knowledge, the thermal transportation behavior of stanene/2D-SiC hetero-bilayer has not yet been studied. Moreover, the interface in a vdWH plays a significant role in the thermal transport. Prior studies revealed that the low thermal conductance at the hetero-bilayer interface usually hinders the efficient thermal transport and heat elimination. In addition, the thermal properties of constituents in a vdW heterostructure vary with respect to the adjacent layers with which they are supported. For instance, Mortazavi et al [44] evaluated the infinite length in-plane thermal conductivity of monolayer graphene as 3000 ± 100 ${\text{ }}W{\text{ }}{m^{ - {\text{ }}1}}{\text{ }}{K^{ - {\text{ }}1}}$ by using the molecular dynamics (MD) simulation. On the other hand, when graphene is used as a supporting layer in graphene/stanene [45], graphene/MoS₂ [46] and graphene/C₃N [43], its lateral thermal conductivity reduces to 685.4 ${\text{ }}W{\text{ }}{m^{ - {\text{ }}1}}{\text{ }}{K^{ - {\text{ }}1}}$, 2898 ${\text{ }}W{\text{ }}{m^{ - {\text{ }}1}}{\text{ }}{K^{ - {\text{ }}1}}$and 1475.20 ${\text{ }}W{\text{ }}{m^{ - {\text{ }}1}}{\text{ }}{K^{ - {\text{ }}1}}$, respectively. Therefore, investigating only the overall thermal conductivity of a vdWH does not provide a clear explanation of thermal transportation behavior. Hence, it has now become imperative to thoroughly investigate the in-plane and out-of-plane thermal conductivity of this heterostructure as well as its constituents before its mass application.

In this work, thermal transports of stanene/2D-SiC have been investigated by classical MD simulation. Firstly, the reverse non-equilibrium molecular dynamics (RNEMD) method is introduced to characterize the lateral thermal conductivity (k). The length of the heterostructure is changed from 25 nm to 400 nm with the intention of elucidating the length dependence thermal conductivity. Afterwards, we evaluated the effect of temperature on 'k' by varying the temperature from 100 K to 700 K. Usually, the quantum effect on the thermal conductivity calculations is neglected by classical MD simulation at lower temperatures. The effect of quantum correction on in-plane thermal conductivity is also studied. On the other hand, the transient pump–probe method is introduced to calculate the flexural interfacial thermal resistance (R). The effects of temperature, system size (area) and coupling strength on 'R' are investigated to realize the thermal transportation at the interface. Moreover, various phonon modes are investigated to reveal the principal mechanisms of thermal transportation in the lateral and flexural direction of the hetero-bilayer.

A supercell was constructed to build the stanene/2D-SiC hetero-bilayer structure. Since the lattice constant of stanene (4.5 ${\mathop A\limits^ \circ }$) [28] is about 33% more than the lattice constant of 2D-SiC (3.095) [28], a lateral periodicity of 3 $\times$ 3 2D-SiC and a lateral periodicity of 2 $\times$ 2 stanene were considered. The lattice constant of the stanene/2D-SiC hetero-bilayer was taken as ${a_{{\rm{Sn}}/{\rm{SiC}}}} = 3{a_{{\rm{SiC}}}} = 9.285 \mathop A\limits^ \circ$ with a tensile strain of only ∼1% imposed on stanene to preserve the commensurability. Thus, the stanene/2D-SiC supercell contains nine carbon atoms, nine silicon atoms and eight Sn atoms, as presented in figure 1(a). To obtain the desired structural size, the supercell cell was replicated along the x and y directions. The RNEMD [47] method, where heat flux is enforced on the structure and the temperature gradient is obtained from the structure, was introduced to determine the lateral thermal conductivity of the stanene/2D-SiC heterstructure. The size of the system is an important issue in the MD simulation. We implemented periodic boundary conditions in all x, y and z directions with a view to avoiding size effects in the RNEMD simulation. Initially, the van der Waals distance between stanene and 2D-SiC was set as 3.3 $\mathop A\limits^ \circ$. The distance between these two monolayers will be attuned automatically after the relaxation of the hetero-bilayer. A vacuum space of 20 $\mathop A\limits^ \circ$ was introduced along the direction perpendicular to the surface of the heterostructure with the purpose of avoiding interaction of the sample with its image [28]. The optimized Tersoff potentials proposed by Cherukara et al [11] and Albe et al [48] were used to describe the interatomic interaction between the Sn–Sn bond in the stanene layer and the Si–C bond in the 2D-SiC layer, respectively. The van der Waals interaction between the stanene and 2D-SiC monolayer was defined by the conventional 12–6 Lennard Jones (LJ) potential. The equation of the LJ potential can be expressed as

$$\begin{equation}V\left( r \right) = 4{{\chi }}\varepsilon \left[ {{{\left( {\frac{\sigma }{r}} \right)}^{12}} - {\text{ }}{{\left( {\frac{\sigma }{r}} \right)}^6}} \right];r < {r_c}_,\end{equation} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

where $\varepsilon ,$ $r$ and $\sigma$ signify the energy parameter, interatomic distance and distance parameter, respectively. The coupling strength between two layers can be adjusted by changing the parameter ${{ \chi }}$. The universal force field [49] was used to calculate the LJ parameters, where ${\varepsilon _{Si - Sn}}{\text{ }}$ = 20.70 ${\text{ meV}}$ and ${\varepsilon _{C - Sn}}$ = 10.58 ${\text{ meV}}$; ${r_{Si - Sn}}$ = 4.34 $\mathop A\limits^ \circ$ and ${r_{C - Sn}}{\text{ }}$ = 4.11 $\mathop A\limits^ \circ$. Another important consideration in the LJ potential is ${r_c}$, which indicates the cutoff distance. It was set as 14 $\mathop A\limits^ \circ$ for all vdW interactions. A large-scale atomic/molecular massively parallel simulator [50] was employed to accomplish all the simulations in this work which has been broadly utilized to characterize the thermal transportation of diverse 2D nanostructures [51–53].

To characterize the in-plane thermal conductivity, the hetero-bilayer structure was first relaxed at the desired temperature and then the RNEMD method was applied to the structure. The schematic of the RNEMD setup on the stanene/2D-SiC heterostructure is depicted in figure 1(b). For the relaxation of the system, the structure was minimized initially by using the conjugated gradient algorithm. After energy minimization, an isothermal–isobaric (NPT) ensemble was performed for 0.05 ns with a time step of 0.5 fs to preserve the value of pressure and temperature constant at the desired value. After that, the structure was placed under a microcanonical (NVE) ensemble along with velocity rescaling for 0.1 ns in order to raise the temperature of the structure from the initial value to the preferred value. The temperature and pressure of the system were stabilized by applying an additional NPT integration for 0.25 ns. In order to check whether the system was relaxed or not, the structure was placed under the NVE ensemble for another 0.15 ns. If the temperature of the system does not change with time, it indicates that the system is equilibrated at that particular temperature. After the relaxation of the system, the whole structure was placed under the RNEMD simulation. For the application of RNEMD into the system, a 'hot slab' was introduced at a distance of L/4 of the nanosheet and a 'cold slab' was introduced at a distance of 3 L/4; where L is the total length of the structure. Then a small amount of heat $\Delta \varepsilon$ was imposed into the 'hot slab' at every 0.05 fs time step. On the other hand, the $\Delta \varepsilon$ was taken out from the 'cold slab' at every 0.05 fs time step. Consequently, the temperature was increased at the hot region and likewise the temperature at the cold region was decreased through a velocity rescaling procedure. The heat flux J which is responsible for this temperature variation at the hot slab and cold slab can be calculated by

$$\begin{equation}J = {\text{ }}\frac{{\Delta \varepsilon }}{{2A\Delta t}},\end{equation} \tag{ 2 }$$

where A and $\Delta t$ indicate the cross-sectional area and time step, respectively. The factor 1/2 is incorporated owing to the periodicity of the arrangement which signifies that the thermal energy can move in both directions. The velocity of every atom will be changed with the heat addition into the hot slab and removal from the cold slab. The temperature at any slab 'a' can be calculated by

$$\begin{equation}{T_a} = {\text{ }}\frac{1}{{3{N_a}{k_B}}}{\text{ }}\mathop \sum \limits_{b = 1}^{{N_a}} {m_b}{v_b}^2,\end{equation} \tag{ 3 }$$

where k_B, v_b and m_b represent the Boltzmann constant, velocity and mass of atom 'b' in the slab, respectively. N_a indicates the total number of atoms in slab 'a'. After imposing the heat flux, NVE integration was applied for 0.25 ns with a view to achieving the perfect transference of heat between the 'hot slab' and the 'cold slab'. In order to achieve a time-averaged temperature profile, NVE integration was applied for another 12 × 10⁶ time steps. The whole structure was split into 'n' slabs towards the x direction with the intention of calculating the temperature gradient accurately. The temperature gradient $\partial T/\partial x$ can be calculated easily from this time-averaged temperature profile of every slab. After obtaining the temperature gradient, the in-plane thermal conductivity k was evaluated based on Fourier's law of heat conduction, which is as follows:

$$\begin{equation}k = {\text{ }}\frac{J}{{\partial T/\partial x}} = {\text{ }}\frac{{\Delta \varepsilon }}{{2A\Delta t(\partial T/\partial x)}}\end{equation} \tag{ 4 }$$

The thickness of a 2D material is an important factor when calculating k. Normally, the thickness of a 2D material is calculated by taking the height of the monolayer [54, 55]. However, the thickness of a hetero-bilayer structure needs to be cautiously selected because there is a vdW distance between two monolayers [46, 56–58]. In this case, the total thickness of the hetero-bilayer is considered as the summation of the thicknesses of two individual monolayers. In our work, the thickness of 2D-SiC, d_2D-SiC was considered as 3.5 Å and the thickness of stanene, d_Sn was taken as 4.5 Å, which have been extensively used in previous works [11, 56]. Thus, the total thickness of the stanene/2D-SiC hetero-bilayer, d_Sn/2D-SiC will be the summation of d_2D-SiC and d_Sn, which is equal to 8 Å. In contrast, Wu et al [59] recently proposed that the monolayer thickness of different 2D materials should not be altered with the deviation of the materials. They also discussed that the thickness of graphene, which is 3.35 Å, should be used for the thermal conductivity calculation of all 2D materials. According to their work, the thickness of both stanene (d_Sn) and 2D-SiC (d_2D-SiC) is 3.35 Å. Thus, 6.7 Å will be the overall thickness of the hetero-bilayer (d_Sn/2D-SiC). In this work, we used both these selections of thickness and consequently calculated the in-plane thermal conductivity.

The phonon spectra of both stanene and 2D-SiC monolayers were calculated to elucidate the thermal transportation behavior in the stanene/2D-SiC bilayer. To estimate the phonon modes, the velocity autocorrelation function (VACF) was calculated by considering the initial velocity $v\left( 0 \right)$ of an atom and velocity of that particular atom $v\left( t \right)$ after a time interval t. The phonon density of states (PDOS) was evaluated by conducting the Fourier transform of the VACF. The equation of PDOS in terms of VACF can be represented as

$$\begin{equation}F\left( \omega \right) = \frac{1}{{\sqrt {2\pi } }}\mathop \smallint \limits_{ - \infty }^\infty dt{e^{i\omega t{\text{ }}}}\frac{{{v\left( 0 \right) \cdot v\left( t \right)}}}{{{v\left( 0 \right) \cdot v\left( 0 \right)}}},\end{equation} \tag{ 5 }$$

where $F\left( \omega \right)$ represents the PDOS which is a function of frequency $\omega$. The phonon power spectral analysis provides a clear idea of thermal transport in any semiconductor as the phonon is the main heat energy carrier here.

The thermal resistance at the interface R is considered as the main parameter to investigate the flexural heat conduction proficiency in a hetero-bilayer. The transient pump-probe technique is an extensively applied technique to characterize R. Figure 1(c) portrays the arrangement for implementing the transient pump-probe scheme in the stanene/2D-SiC heterostructure. This method has previously been used to characterize the flexural thermal transportation in graphene/silicon [60], graphene/h-BN [61], silicene/SiO2 [62], phosphorene/silicon [53], graphene/phosphorene [63] and graphene/copper [64]. To measure the R between stanene and 2D-SiC, the structure was relaxed at a particular temperature by following the same steps which have been described at the in-plane thermal conductivity calculation section. After the relaxation of the system, an ultrafast heat impulse was applied to the stanene layer for 50 fs. As a result, the temperature of the stanene layer increases, whereas the temperature of the 2D-SiC layer remains constant. After the removal of the heat impulse, the temperature of stanene will be decreased gradually and the temperature of the 2D-SiC layer will be increased slowly as the surrounding environment is a vacuum. The heat transportation process will be continued as long as the system reaches the thermal equilibrium condition. The temperature progression of stanene (T_Sn) and 2D-SiC (T_2D-SiC) are recorded during the whole thermal equilibrium process. The total energy of stanene (E_t ) at every step is also recorded to calculate thermal resistance R, which is determined from the following formula

$$\begin{equation}\frac{{\partial {E_t}}}{{\partial t}} = {\text{ }}\frac{{{A_{R{\text{ }}}} \cdot \left( {{T_{Sn}} - {T_{2D - SiC}}} \right)}}{R},\end{equation} \tag{ 6 }$$

where A_R represents the contact area between stanene and 2D-SiC. The data points which are collected at every time step should be averaged after every 100-time steps in order to minimize data noise.

Recently, the thermal conductivity of 2D-SiC and the stanene monolayer has widely been studied [9, 11, 19]. These studies revealed that 2D-SiC offers higher thermal conductivity (almost twice in order) than stanene. Since a hetero-bilayer consisting of these two monolayers has been synthesized recently [28], it has become an important issue to determine the thermal behaviors of this heterostructure. In this work, RNEMD is performed to estimate the in-plane thermal conductivity of the individual stanene film, the 2D-SiC layer, and finally the hetero-bilayer. The calculated temperature profile utilizing the RNEMD simulation is depicted in figure 2. There are two regions at the temperature distribution profile, one is linear and the other one is non-linear which is well matched with previous MD studies [45, 65, 66]. The non-linear region is observed due to the rapid energy exchanges between potential and kinetic energies, which is close to the hot slab and cold slab. This non-linear zone is generally neglected at the time of calculating the temperature gradient. It is always preferred that the thermal conductivity be calculated by considering the temperature gradient attained from the linear region, which is colored in red for better understanding.

Zoom In Zoom Out Reset image size

At first, we calculated the length dependence in-plane thermal conductivity of the stanene/2D-SiC hetero-bilayer. The length towards the lateral direction is changed from 25 nm to 400 nm. The width of the system is taken to be fixed as 10 nm to reduce the size effect. The thermal transportation behavior of 2D-SiC, stanene, and the bilayer is calculated separately. It has been perceived that the value of heat flux needed to maintain a constant temperature variation at the heat source and heat sink is different for stanene and 2D-SiC. Stanene requires lower energy with respect to 2D-SiC to increase a certain temperature at the heat source, which indicates that the thermal transportation capability of 2D-SiC is greater than that of stanene. The thermal conductivities of stanene, 2D-SiC and the bilayer were evaluated for a structure sized (400 x 10) nm² to be 9.7 ${\text{ }}W{\text{ }}{m^{ - 1}}{K^{ - 1}}$, 132.91 $W{\text{ }}{m^{ - 1}}{K^{ - 1}}$ and 58.34 ${\text{ }}W{\text{ }}{m^{ - 1}}{K^{ - 1}}$, respectively.

The length dependence of k is illustrated in figure 3. Two types of thicknesses, such as vdW thickness and unified thickness, were considered, as depicted in figures 3(a) and (b), respectively. Both of these figures reveal that the in-plane thermal conductivity increases with the increasing length. When the length of the hetero-bilayer increases from 30 nm to 400 nm, the thermal conductivity increases monotonically from 13.65 to 58.42 ${\text{ }}W{\text{ }}{m^{ - 1}}{K^{ - 1}}$. Initially, the thermal conductivity increases linearly with the increasing length, but after a certain length of the heterostructure, the increasing slope tends to be flat. The phonon mean free path (MFP) plays an important role in the thermal conductivity of a material. Normally, thermal transportation will be diffusive when the system length is greater than the MFPs of the phonons. In contrast, if the length of the nanostructure is shorter than the MFPs of the phonons, the thermal transport will be ballistic. Initially, the in-plane thermal conductivity in this work increases linearly with the increasing length due to the participation of more phonons whose MFPs are greater than the system length. Afterward, MFPs of the phonons become smaller with respect to the increasing system size which impede the linear increment of the thermal conductivity.

Zoom In Zoom Out Reset image size

Figure 3 also shows that the k value of stanene is much lower than the value of 2D-SiC. Recently, Che et al investigated various single-component systems and binary systems [67]. They reported that the chances of phonon–phonon scattering are dominant in a two-atomic system rather than a single-element system and concluded that the thermal conductivity of a binary system is 37% less than a single-component system. According to their research, the thermal conductivity of stanene should have been greater than that of 2D-SiC. However, these results are in constrast to the present work since the phonon frequency of stanene (∼6 THz) is very small with respect to that of 2D-SiC (∼42 THz), which is depicted in figure 4. At the same time, the high atomic mass and the buckling structure of stanene might also be responsible for this kind of discrepancy. Previously, Islam et al [19] calculated the thermal conductivity of 2D-SiC by using the same RNEMD method [47]. They predicted that the participation of low-frequency transverse acoustic (TA), longitudinal acoustic (LA) as well as the participation of flexural acoustic (ZA) phonons increases with the increasing sheet length. The observation is also compatible with our study and from figure 5, it is clear that the participation of various phonon modes increases with the size of the system. As a result, the value of the in-plane thermal conductivity of 2D-SiC as well as the stanene/2D-SiC hetero-bilayer rises with the growing sheet length.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also used a quantitative method to obtain the thermal conductivity of our hetero-bilayer. As the heat flux required to calculate k in the bilayer is the summation of the heat fluxes required to generate the temperature gradient at stanene and 2D-SiC, it is very common to calculate the in-plane thermal conductivity of the hetero-bilayer by using the following formula developed by Liu et al [46]

$$\begin{equation}{q_B} \cdot {\text{ }}{A_B} = {\text{ }}{q_{Sn}} \cdot {\text{ }}{A_{Sn}} + {\text{ }}{q_{2D - SiC}} \cdot {A_{2D - SiC}},\end{equation} \tag{ 7 }$$

where q = $- {\text{ }}k \cdot \Delta T$, q indicates the heat flux density. A denotes the area of the structure where $A = d \cdot w$. Here, w and d indicate the width and thickness of the structure, respectively. Since the temperature gradient and width are the same for stanene, 2D-SiC and the bilayer, equation (7) can be simplified as follows,

$$\begin{equation}{k_B} \cdot {d_B} = {\text{ }}{k_{Sn}} \cdot {d_{Sn}} + {\text{ }}{k_{2D - SiC}} \cdot {d_{2D - SiC}}\end{equation} \tag{ 8 }$$

As the values of ${\text{ }}{k_{Sn}}$, ${\text{ }}{d_{Sn}}$, ${k_{2D - SiC}}$, ${d_{2D - SiC}}$ and ${d_B}$ are known, it is pretty much straight-forward to calculate the thermal conductivity of the hetero-bilayer from equation (8). The estimated values are also depicted in figure 3, which coincide well with our calculated values from the MD simulations.

It is also clear from figure 3 that the slope of the thermal conductivity of stanene is not significant. In contrast, the thermal conductivities of 2D-SiC, as well as the hetero-bilayer, tend to increase exponentially until L > 64.17 nm. From these data, it is apparent that the size effect on the thermal conductivity is stronger for 2D-SiC than stanene. Therefore, it is important to calculate the infinite length thermal conductivity of 2D-SiC as well as the hetero-bilayer. The infinite length thermal conductivity of stanene, 2D-SiC, and the bilayer can be extracted by matching various dimensions to the MD outcomes. The correlation between the total length and the thermal conductivity can be derived by the combination of Matthiessen's rule and kinetic theory [68, 69]. The equation which is used to fit the MD results can be presented as:

$$\begin{equation}\frac{1}{k} = {\text{ }}\frac{1}{{{k_\infty }}}\left( {1 + {\text{ }}\frac{{{\text{ }}{L_{MFP}}}}{L}} \right),\end{equation} \tag{ 9 }$$

where ${L_{MFP}}$ indicates the infinite length MFP, ${k_\infty }$ indicates the infinite-length thermal conductivity, $L$ represents the length of the structure. The predicted values from MD simulations are fitted from the length of 90 nm–400 nm. The infinite-length thermal conductivity can be predicted by obtaining the relation between 1/k and 1/L. The relation between these parameters is depicted in figures 6(a) and (b). When the vdW thickness is used, the infinite length of the thermal conductivities of stanene, 2D-SiC, and the bilayer are 13.79 ${\text{ }}W{\text{ }}{m^{ - 1}}{K^{ - 1}}$, 151.75 ${\text{ }}W{\text{ }}{m^{ - 1}}{K^{ - 1}}$ and 66.67 ${\text{ }}W{\text{ }}{m^{ - 1}}{K^{ - 1}}$. In addition, the infinite length thermal conductivity of the bilayer changes to 79.62 $W{\text{ }}{m^{ - 1}}{K^{ - 1}}$ when a unified thickness is used. It is important to mention that the infinite length in-plane thermal conductivity of 2D-SiC is predicted in previous MD simulations as 313 ${\text{ }}W{\text{ }}{m^{ - 1}}{K^{ - 1}}$, which is greater than the value we have investigated in our work. The thermal conductivity of 2D-SiC reduces because of the lessening of the participation of ZA phonons. The reason behind this reduction of ZA phonons is the enhanced phonon scattering between the interlayers of the bilayer system [70, 71]. When 2D-SiC is used as a single monolayer, no interface coupling coefficient factor is introduced. Therefore, the in-plane and out-of-plane phonons in a free-standing 2D-SiC are well decoupled. In contrast, when 2D-SiC is used as a supporting layer with a stanene monolayer, various parameters (for example, rotation, translation, and reaction) of 2D-SiC as well as stanene are fragmented. Consequently, the phonon properties of the system will be altered by this stanene/2D-SiC interaction and a new scattering channel will be introduced, which in turn increases the interaction of flexural as well as in-plane phonons. This kind of phenomenon was also observed in previous studies of various heterostructures [43, 45, 46]. For instance, the infinite lengths of the thermal conductivity of graphene and C₃N monolayer were predicted as 3000 ± 100 ${\text{ }}W{\text{ }}{m^{ - 1}}{K^{ - 1}}$ and 820 ${\text{ }}W{\text{ }}{m^{ - 1}}{K^{ - 1}}$, respectively, whereas, the infinite length thermal conductivity of graphene and C₃N in a graphene/C₃N [43] heterostructure was predicted recently as 1475.20 ${\text{ }}W{\text{ }}{m^{ - 1}}{K^{ - 1}}$ and 775.19 ${\text{ }}W{\text{ }}{m^{ - 1}}{K^{ - 1}}$, respectively. Similarly, the infinite length thermal conductivity of graphene in a graphene/stanene [45] heterostructure was predicted as 685.4 ${\text{ }}W{\text{ }}{m^{ - 1}}{K^{ - 1}}{\text{ }}$, which is again not only less than the value of single-layer graphene but also less than the value of supported graphene in a graphene/C₃N hetero-bilayer. This is because the phonon MFP in stanene is shorter than that of C₃N, which ultimately increases the phonon scattering in the graphene/stanene hetero-bilayer and reduces the thermal conductivity of graphene and the overall bilayer as well.

Zoom In Zoom Out Reset image size

2D materials are predicted to be the key elements of future nanodevices and can often be placed under extreme temperature conditions. The high temperature can alter the thermal conductivity of these materials and can consequently cause the failure of, or permanent damage to, the nanoelectronic devices. Therefore, the variation of the thermal transportation of these materials with temperature is very important to determine. The thermal conductivity variation of the 2D-SiC, stanene and bilayer are depicted in figure 7. To predict the thermal transport behavior with temperature, a system of (10 x 100) nm² is taken into consideration. The green symbol in figure 7 represents the thermal conductivity of the hetero-bilayer predicted by equation (8). When the temperature increases from 100 K to 700 K, the thermal conductivity decreases (from 160.28 to 56.725 $W{\text{ }}{m^{ - 1}}{K^{ - 1}}$ for 2D-SiC, from 10.12 to 0.53 $W{\text{ }}{m^{ - 1}}{K^{ - 1}}$ for stanene, and from 70.11329 to 24.814 $W{\text{ }}{m^{ - 1}}{K^{ - 1}}$ for the bilayer). The decreasing path of stanene follows the 1/T law which coincides well with the previous studies for single-layer graphene, phosphorene and stanene [10, 72]. In contrast, it can be clearly observed from figure 7 that the value of the thermal conductivity of the bilayer decreases monotonically with the increasing temperature. The monotonic decreasing behavior of 2D-SiC and the bilayer can be described by the PDOS, which is presented in figure 8. The TA, LA and ZA mode phonons which occupy low frequency regions contribute the most to the transportation of thermal energy in binary 2D materials [19, 73]. The thermal conductivity decreases due to the reduction of low-frequency phonons with the increasing temperature. On the other hand, at a higher temperature, the contribution of high-frequency phonons becomes dominant and phonon–phonon Umklapp scattering occurs. Due to these phenomena at a higher temperature, reduced thermal conductivity and a slow decreasing rate have been observed in 2D-SiC and the bilayer as well.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The way in which 2D materials exist in practical application is different from our simulation environment. In real life, it is well known that the phonons undergo a freezing condition at low temperature [74, 75], which means they will be in a resting position until the surrounding is heated up to a sufficient temperature. The limit is called the Debye temperature. Therefore, the thermal conductivity calculated from MD is not valid before this Debye limit and it is desirable to predict thermal conductivity at low temperature by the quantum correction method. At first, normal mode-specific heat is calculated from the PDOS spectra of the hetero-bilayer. The calculation of specific heat can be represented by,

$$\begin{equation}{C_v} = {\rm{ }}\frac{{3N{K_B}{\rm{ }}\int_0^\infty {\frac{{{u^2}{e^u}}}{{{{\left( {{e^u} - 1} \right)}^2}}}.F\left( \omega \right).{\rm{ }}d\omega } {\rm{ }}}}{{\int_0^\infty {F\left( \omega \right).{\rm{ }}d\omega } }}{\rm{ }},\end{equation} \tag{ 10 }$$

where C_v represents the specific heat, $F\left( \omega \right)$ represents the phonon spectrum, ${K_B}{\text{ }}$ indicates the Boltzmann constant and u is equal to $\frac{{\hbar \omega }}{{{K_B}T}}$. The value of specific heat is presented in figure 9, where it is clear that the value of C_v increases with the temperature, and when the phonons are fully excited at the higher temperature, specific heat becomes saturated. Finally, the quantum corrected thermal conductivity, k_QC of the stanene/2D-SiC hetero-bilayer is calculated as,

$$\begin{equation}{k_{QC}} = {\text{ }}\frac{{{C_v}}}{{3N{K_B}}} \times {k_{MD}},\end{equation} \tag{ 11 }$$

Zoom In Zoom Out Reset image size

where ${k_{MD}}$ represents the value of thermal conductivity obtained from the MD simulation. Figure 10 elucidates the difference between the ${k_{QC}}{\text{ }}$ and ${k_{MD}}$. At first, the in-plane thermal conductivity increases with the temperature due to the increasing of excited phonons. No sooner had the system reached the Debye temperature limit than the thermal conductivity started decreasing with temperature. It is because of the reduction of ZA mode phonons at a higher temperature due to the increasing amount of phonon–phonon scattering [73, 76].

Zoom In Zoom Out Reset image size

Despite having a large thermal property along the lateral direction, most of the 2D heterostructures show low thermal transportation along the flexural direction. The in-plane thermal conductivity is much higher owing to the covalent sp² bonding between atoms. On the other hand, there are no sp² bondings towards the flexural direction of the hetero-bilayer. A weak vdW interaction is present merely to join the adjacent layers in a heterostructure. To characterize the flexural thermal resistance at the interface, a structure with a size of (100 x 10) nm² is taken. The periodic boundary condition is used along the x, y, and z direction to eradicate the size effect on the structure. After the relaxation of the structure at 300 K, a heat impulse of 2.66 $\times$ ${10^{15}}$ W ${m^{ - 2}}$ is enforced to the stanene monolayer for 50 fs. This high amount of thermal impulse is imposed to ensure a large temperature difference between stanene and the 2D-SiC layer. A large temperature difference is required to obtain a smooth and relaxed temperature and energy evolution graphs. The predicted value of interfacial resistance totally depends on these temperature and energy evolution profiles of the material. However, owing to this heat impulse, the temperature of stanene increases rapidly to ∼750 K after 50 fs while the temperature of 2D-SiC remains constant at ∼300 K. Afterwards, heat is removed and the structure reaches a stable condition at ∼400 K. The temperature of stanene (${T_{sn}}$), the energy of stanene (E_t ) as well as the temperature of 2D-SiC (${T_{2D - SiC}}$) are documented for another 5 × 10⁵ time steps. Figure 11 depicts the temperature evolution of 2D-SiC and stanene which are plotted against the right y-axis as well as the energy evolution of stanene which is plotted against the left y-axis. By using the integral form of equation (6), the thermal resistance can be easily calculated from these temperature and energy progression profiles. The integral formula of equation (6) can be represented by:

$$\begin{equation}{E_t} = {\text{ }}{E_0} + {\text{ }}\frac{{{A_R}}}{R}{\text{ }}\mathop \smallint \limits_0^t \left( {{T_{sn}} - {\text{ }}{T_{2D - SiC}}} \right) \cdot dt,\end{equation} \tag{ 12 }$$

Zoom In Zoom Out Reset image size

where the initial energy is denoted by ${E_0}$. The predicted value of R from the transient pump-probe technique is 3.0622 $\times$ ${10^{ - 7}}$ K.m² W⁻¹ for the structure area of 1000 nm².

We discussed in an earlier section that the thermal interface materials are used universally to manufacture FETs and other nanoelectronic devices. As they are getting smaller day by day, it has become a challenging task to manage the thermal energy inside these devices. At the same time, contact pressure variation caused by the condensed arrangement of these layered materials can affect the thermal conductivity. Thus, to design an effective nanodevice, it is necessary to understand the characteristics of interfacial resistance with different sizes, temperatures as well as various coupling factors.

The variation of R with various structural sizes is depicted in figure 12. The value of R increases until the area of the structure reaches around 250 nm². It is important to mention that the value of thermal resistance becomes constant after this certain area. All the values are predicted at 300 K and by taking the coupling factor of ${{\chi }} = 1$. The variation of R with various temperature conditions is revealed in figure 13. It is perceived that the thermal resistance stays almost constant after a size of 250 nm². Thus, to realize the variation of R with temperature, a structure of 300 nm² is used. The temperature is varied from 100 K to 600 K. For ${{ \chi }} = 1$, the interfacial resistance between stanene and 2D-SiC is decreased by 70.49% from 3.83 $\times$ ${10^{ - 7}}$ K.m² W⁻¹ to 1.13 $\times$ ${10^{ - 7}}$ K.m² W⁻¹. The phonon modes are not properly excited at low temperature but with the increasing temperature they become excited, and high-frequency phonons are converted to a large volume of multifold low-frequency phonons. Normally, low-frequency phonons along the out-of-plane direction are the main heat carrier. These increased low-frequency phonons increase the probability of inelastic scattering among the phonon modes. Moreover, an earlier study reported that inelastic scattering plays a vital role in transporting heat energy. As a result, phonon thermal transport efficiency increases at higher temperatures which leads to a lowered thermal resistance for a vdWH. Another important factor of increasing thermal conductance with the increasing temperature along the flexural direction is anharmonicity inside the materials. The number of anharmonic scattering increases at a higher temperature, which increases the probability of heat energy rearrangement for various low-frequency phonons [77]. These kinds of phenomena at higher temperatures can lead to a high heat transfer rate across the out of plane direction of any hetero-bilayer.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also studied the cross-plane thermal resistance with respect to the coupling strength between stanene and the 2D-SiC monolayer. From figure 14, the R value decreases with the increasing coupling strength. The value of R is predicted by taking the coupling factor of ${{\chi }} = 0.5,{\text{ }}1,{\text{ }}2,{\text{ }}3,{\text{ }}4{\text{ and }}5$. When the coupling factor is set to ${{\chi }} = 0.5$, the value of R is 4.23 $\times$ ${10^{ - 7}}$ K.m² W⁻¹ and this value is reduced by 50.118% when the coupling factor is set to ${{\chi }} = 5$. There are two reasons for this kind of reduction of thermal resistance. (i) With the increasing coupling strength between stanene and 2D-SiC, the contact pressure between them will be increased. Consequently, the phonon coupling at the interface toughens, which directly increases the chance of heat transfer across the interface and also diminishes the interfacial resistance. (ii) The coupling between the lateral and flexural phonons of 2D-SiC is increased with the increasing coupling strength. This occurs because the Sn atoms of the stanene monolayer are now occupied as a scattering center of 2D-SiC which indirectly smooths the thermal conductivity through the interface. The trend of the decrement of thermal resistance is incongruous to the previous studies of MoSe₂/MoSe₂ [78], graphene/phosphorene [63], phosphorene/silicon [53] and graphene/C₃N [43].

Zoom In Zoom Out Reset image size

In summary, the lateral and flexural thermal conductivity of stanene/2D-SiC vdWH were explored via the RNEMD simulation and transient pump-probe approach. The effects of system dimension, interlayer coupling strength and temperature on the in-plane thermal conductivity (k) as well as the flexural thermal resistance (R) were studied. The phonon modes in standalone stanene and 2D-SiC as well as stanene/2D-SiC were calculated with a view to elucidating the underlying mechanism of thermal transport in stanene/2D-SiC vdWH. The predicted infinite length in-plane thermal conductivity of the stanene/2D-SiC bilayer was 66.67 ${\text{ }}W{\text{ }}{m^{ - 1}}{K^{ - 1}}$ by using the actual thickness of the system, whereas this value changes to 79.62 ${\text{ }}W{\text{ }}{m^{ - 1}}{K^{ - 1}}$ when a unified thickness was used. When the temperature increases from 100 K to 700 K, the in-plane thermal conductivity decreases from 70.11329 to 24.814 ${\text{ }}W{\text{ }}{m^{ - 1}}{K^{ - 1}}$. The PDOS revealed that the contribution of high-frequency phonons becomes dominant at a higher temperature, which in turns increases phonon–phonon Umklapp scattering and hence reduces the thermal conductivity. The thermal conductivity with quantum correction at low temperature exhibited a rising characteristic up to the Debye limit as a result of the contribution of ground state phonons in specific heat. The computed interfacial thermal resistance between stanene and 2D-SiC was 3.0622 $\times$ ${10^{ - 7}}{\text{ }}$ K.m² W⁻¹, which is close to other 2D vdWHs. Interface coupling strength and temperature can both lessen the R of the heterostructure. With the rising temperature from 100 K to 600 K, the R between stanene and 2D-SiC is lessened by 70.49% from 3.83 $\times$ ${10^{ - 7}}$ K.m² W⁻¹ to 1.13 $\times$ ${10^{ - 7}}$ K.m² W⁻¹. Likewise, the value of R is reduced by 50.118% when the coupling factor is increased from ${{\chi }} = 0.5$ to ${{\chi }} = 5$. These findings provide a deep understanding into thermal transport in the stanene/2D-SiC vdWH, which might be helpful for the superior design of nanoelectronic devices using this heterostructure.

There are no conflicts to declare.

The data that supports the findings of this study are available within the article.