Interfacial coupling effects on the thermal conductivity of few-layer graphene

Advanced nanotechnologies greatly promote the development of several devices with lengths in the nanometer order. However, one of the critical problems for such small electronic devices is the efficient management and dissipation of the generated excessive heat during their operation process. Two-dimensional materials have attracted a massive interest due to their unique and fascinating properties discovered over the past decades [1]. Graphene, the most famous member of this family with excellent electronic, optoelectronic, and thermal properties, can provide innovative solutions in the microelectronics field [2, 3]. Since the pioneering work of Novoselov et al in 2004 [4], the exploration of new materials, such as hexagonal boron nitride, molybdenum disulfide, and black phosphorus, is progressing constantly. It is therefore urgent to precisely evaluate the thermal transport properties of thin graphene films, in order to push the development of graphene-based electronic devices.

As reported in several experimental works [5–7], the single layer of graphene had a thermal conductivity of 5000 W m⁻¹ K⁻¹ [5], which highly depended on both temperature and flake size. The theoretical thermal conductivity obtained from the molecular dynamics (MD) simulations reached >10 000 W m⁻¹ K⁻¹ [8]. However, graphene can either be supported or embedded when it is integrated into a device in the majority of the current practical applications. Only in rare cases it is used in its suspended state. According to the experimental report from Seol et al [9], the thermal conductivity of supported graphene was measured to be only 600 W m⁻¹ K⁻¹, which was much lower than that of its suspended state. Jang et al [10] showed that the thermal conductivity of encased graphene was even lower than 160 W m⁻¹ K⁻¹. Thus, it is of pivotal importance to explore the thermal conductivity of graphene supported by a substrate and to understand the fundamental physics that controls the thermal transport. Several previous investigations [11–14] have shown that the coupling between graphene and a substrate plays an important role when graphene is in contact with the substrate. The damping of the flexural acoustic (ZA') phonons can lead to the reduction of the thermal conductivity [11, 12, 14]. One can attribute such unusual thermal transport phenomena of the graphene-substrate interface to the phonon coupling [15].

A few-layer graphene (FLG) sample confined between amorphous oxides exhibited an enhanced thermal conductivity upon the increase in its layer thickness [10]. This result presented a totally opposite correlation between the thermal conductivity and the layer thickness when compared with the results of suspended FLG [16, 17]. Moreover, it remains unclear whether the thermal conductivity reduction due to phonon scattering is dominated by the interface roughness or the coupling strength [11, 18, 19]. Besides, the substrates previously investigated are usually amorphous, and most of them are silicon dioxides. Crystalline silicon is commonly employed in microelectronic or photovoltaic industries since it plays an important role in determining the specific properties of silicon-based devices. In this work, the effects of the layer thickness, the substrate roughness, and the strength of graphene-substrate coupling and interlayer coupling on the thermal conductivity of graphene supported on a crystalline silicon substrate are investigated.

Non-equilibrium molecular dynamics (NEMD) simulations were performed to calculate the thermal conductivities of a series of supported and suspended FLG films. Figures 1(a), (b) shows the simulation model for the one single layer graphene (SLG) sample supported on a smooth silicon surface. The periodic boundary conditions were applied along the y-direction, whereas fixed boundaries were applied along the x-direction. No boundary conditions were defined for the out-of-plane direction (or z-axis). Figure 1(c) shows the atomic construction of graphene supported on a rough silicon substrate. In the case of a real silicon surface, the roughness is random and the contact size can vary over a broad range. To simplify the simulation process as well as the calculations, periodical grooves with a width of l = 0.54 nm and a spacing of l_p = 1.08 nm were designed on the surface of the silicon substrate. The groove depth was set to h = 0.54 nm.

Zoom In Zoom Out Reset image size

The graphene was initially placed on top of the crystalline silicon with a separation distance (d) of 0.31 nm, while the bottom atoms of the silicon substrate were fixed. After the system reached the equilibrium at 300 K, the equilibrium separation distance between the supported graphene and the substrate could be obtained. In order to obtain the equilibrium structures, the temperature of the system was adjusted by performing a constant temperature simulation with 8 × 10⁶ time steps. After relaxation,    two Langevin heat baths at 320 K and 280 K were applied to the regions of the samples located near the fixed walls to simulate the behavior of the heat source and sink, respectively. Different regions in the simulation system, including the fixed wall and the thermal bath (heat source/heat sink), are shown in figures 1(a), (b). In order to stabilize the heat flow under temperature gradient conditions, the simulations were performed with an additional 5 × 10⁶ time steps. In these conditions, the heat flux passing through the system is time independent and the temperature gradient can be clearly defined. At last, the thermal conductivity (κ) can be calculated as follows:

$$\begin{eqnarray}&&\kappa =-\displaystyle \frac{q}{{\rm{\nabla }}T}\end{eqnarray} \tag{ 1 }$$

where ∇T is the temperature gradient and q is the heat flux density across the FLG along the x-direction. The q value can be calculated from the heat source or the heat sink via the expression:

$$\begin{eqnarray}&&q=\displaystyle \frac{E}{{\rm{S}}\cdot t}\end{eqnarray} \tag{ 2 }$$

where S is the cross section area calculated depending on the number of layers n, the thickness of an SLG, and the width of the sample; E corresponds to the energy added into the heat source or removed from the heat sink; t is the effective simulation time. Only the heat flux density transported in the FLG region was recorded in the simulations for the supported FLG. Since the obtained heat flux values of the heat source and heat sink were almost identical, their average was used for the calculations. The temperature gradient was calculated according to the temperature of the FLG sample along the x-direction. The average temperature of a thin slice of 0.2 nm in length centered at the position x can be obtained by using the following expression:

$$\begin{eqnarray}&&T=\displaystyle \frac{\displaystyle {\sum }_{i=1}^{N}{m}_{i}{\overrightarrow{v}}_{i}^{2}}{3N{k}_{{\rm{B}}}}\end{eqnarray} \tag{ 3 }$$

where N is the total number of atoms in the slice, m_i is the mass of the atom i, ${\overrightarrow{v}}_{i}$ is the velocity of the atom i, and k_B is the Boltzmann constant. As a result, T corresponds to the instantaneous temperature at x for a defined time step.

The MD simulations were performed by using the LAMMPS package [20]. The Newton's equations of motion were solved by using the velocity-Verlet integral algorithm with a time step of 0.5 fs. The Stillinger-Weber interatomic potential was used to describe the atomic interaction between silicon-silicon in the substrate [21]. The optimized Tersoff potential was employed to model the interactions between the carbon atoms within the graphene layer [22]. The van der Waals (vdW) force was obtained by using the Lennard-Jones (LJ) potential to determine the interactions between atoms belonging to different graphene layers, as well as between the graphene sample and the substrate. The expression used is the following:

$$\begin{eqnarray}&&{V}_{ij}^{LJ}({r}_{ij})=4\chi {\varepsilon }_{ij}\left[{\left(\displaystyle \frac{{\sigma }_{ij}}{{r}_{ij}}\right)}^{12}-{\left(\displaystyle \frac{{\sigma }_{ij}}{{r}_{ij}}\right)}^{6}\right]\end{eqnarray} \tag{ 4 }$$

where, ε_ij and σ_ij are the LJ potential parameters, χ is a dimensionless proportional factor and assigned to unit by default, and r_ij is the interatomic distance. The LJ potential parameters of the C–C bond are identical to the ones reported in [23]. The LJ parameters for the C–Si interactions were obtained from [24]. The cut-off distance in the LJ potential was set to 2.5 times the value of σ_ij for both the C–C and C–Si interactions [11]. The neighbor list was dynamically updated every 5 fs in the simulations.

Due to the different periodic structures in graphene and the silicon crystal in the simulation system, the lattice constants were adjusted to adapt to the periodic boundary conditions along the y-direction. Here, the lattice constant of crystalline silicon, a, varies from 0.543 nm to 0.54 nm. The thermal conductivities of the silicon samples obtained by simulating the two lattice constants are similar, as shown in figure 2. Therefore, the calculation of the thermal conductivity of graphene is not influenced by the variation of the silicon lattice constants.

Zoom In Zoom Out Reset image size

The correlation between the thermal conductivity of the FLG structure and the number of layers at room temperature is shown in figure 3. The thickest graphene layer considered in this study has n = 20. This corresponds to a total number of 183808 atoms including the substrate atoms. The thermal conductivity of a suspended FLG sample was investigated at room temperature, and the results show that the thermal conductivity decreases as the number of layers increases. The thermal conductivity initially is 511.2 W m⁻¹ K⁻¹ in the case of a single layer of graphene. Then, it decreases to 497.6 W m⁻¹ K⁻¹ when 5 layers are present. This value further decreases to 491.3 W m⁻¹ K⁻¹ with 20 layers. The trend for the suspended FLG sample in figure 3 is similar to the results obtained in previous investigations [16, 25, 26]. However, for supported FLG samples on a smooth silicon, the thermal conductivity presents an increase as a function of the number of layers in a similar two-step. The increase rate of thermal conductivity decreases when the number of layers is up to 5. When the FLG is considered as a three-dimensional system, several interactions should be considered due to the dominant role of the damping of the ZA' phonons. These interactions need to be examined not only between graphene and the substrate but also among different graphene layers. The layer-to-layer coupling induces a thermal conductivity reduction in the suspended FLG sample upon an increase in its thickness [16, 25, 26]. For the supported graphene films, the coupling between the graphene and substrate also influence the thermal conductivity. However, as the number of layers increases, the effect of substrate coupling on phonon transport becomes weak, resulting in the thermal conductivity of supported graphene approaching that of bulk graphite.

Zoom In Zoom Out Reset image size

The thermal conductivity of a FLG sample supported on a rough crystalline silicon was calculated. Figure 3 shows that the thermal conductivity values of graphene samples supported on a rough silicon substrate are smaller than those on smooth silicon. With 5 layers, for instance, the rough condition induces a thermal conductivity of 465.1 W m⁻¹ K⁻¹, 2.5% lower than that calculated in smooth conditions. It seems that the contact surface would decrease when the graphene sample is supported on the rough substrate in comparison to the smooth substrate case, which may result in an increasing thermal conductivity. However, since the graphene is extremely flexible along the in-plane direction, the graphene layer at the bottom presents distinct wrinkles along the grooves, which in fact increases the interfacial coupling. This is also consistent with the results of Zhang et al [27], which show that the energy coupling is enhanced when graphene is in contact with a rough silicon substrate. They investigated the energy transport perpendicular to the interface between graphene and the silicon substrate and found that the rough substrate could increase the heat conduction. In our work, the energy transport parallel to the interface between graphene and the substrate was mostly investigated. The enhanced cross-plane interfacial coupling generally results in a reduced thermal conductivity along the in-plane direction. In addition, the deformation of the graphene surface close to the substrate increases the phonon scattering, resulting in a lower thermal conductivity of a graphene sample in contact with a rough surface when compared to a smooth one. It can also be noticed in figure 3 that the convergence of thermal conductivity with the number of graphene layers on rough silicon is a little slower than that of the suspended graphene and the graphene supported on smooth substrate. By analyzing the atomic structure of the model, we found that the fluctuations of the bottom layer of the supported sample on the rough substrate along the thickness direction for 9, 12 and 20 layers were 0.065 Å, 0.0438 Å, and 0.033 Å, respectively. This indicates that the thermal conductivity of supported graphene on rough substrate does still not saturate when the layer number is 20 due to the underneath patterns.

As mentioned, the coupling strength between graphene and the substrate, as well as the interlayer interactions within the multilayer graphene sample, are important factors in determining the thermal transport properties of the material. Here, the correlation between the interaction strength in a graphene film and its thermal conductivity was investigated. The interaction strengths of layer-layer and graphene-substrate were regulated via the dimensionless coupling strength parameter χ in the LJ potential (equation (4)), which in the 1–10 range. Figure 4 shows the coupling strength dependence of the thermal conductivity for both the suspended and supported FLG samples. Note that there are no interfacial vdW interactions within the suspended single layer graphene. Thus, the thermal conductivity of suspended graphene with n = 1 is not dependent on the coupling strength and is labelled by a solid line in figure 4(a). For a fixed value of χ, when increasing the layer number, the thermal conductivity of suspended graphene decreases, while the thermal conductivity of supported graphene increases, which is identical to those reported in figure 3. Upon the increase in the coupling energy, either the thermal conductivity of the suspended or of the supported graphene sample decreases, for all specimens with three different layer numbers. The results on two types of FLG samples are here compared: the suspended FLG specimen with 5 layers presents a decrease of only 6.1% in its thermal conductivity (from 497.6Wm⁻¹ K⁻¹ to 467.0Wm⁻¹ K⁻¹) as the value of χ increases from 1 to 10. In the case of the supported sample with the same number of layers, the decrease in amplitude reaches ∼19% on both a smooth and a rough substrate. This difference clearly indicates the impact of the graphene-substrate interaction to the heat transport inside the graphene films. Given that the increase of the coupling strength between the graphene layers and the substrate can significantly shorten the phonon lifetime of the ZA' modes [28, 29], one may expect that the thermal conductivity of graphene decreases once it is supported on the substrate. Besides, the influence of enhancing the coupling strength is more significant for an SLG sample on supported silicon than for suspended FLG specimens. With three coupling parameters, the thermal conductivity of an SLG on a smooth silicon substrate decreases from 461.4 to 403.5 and 303.5 W m⁻¹ K⁻¹, corresponding to a 34.2% reduction. When the SLG sample is supported on the rough substrate, its thermal conductivity even drops to 281.4 W m⁻¹ K⁻¹, a 36.2% decrease in the amplitude. As a single layer can be easily deformed, increasing the graphene-substrate interaction may result in a closer contact between the two surfaces, especially in the case of a rough substrate with a large number of grooves. This leads to an enhanced vdW force and shows significant effects on the thermal conductivity in graphene layers. It has to be mentioned that according to the previous studies [30, 31], the surface morphology can influence the interfacial coupling strength. In this work, we mainly focused on the relative changes of thermal conductivity among different structures, rather than the exact value of the thermal conductivity. In fact, it is difficult to obtain the thermal conductivity accurately using the MD simulation because of the inaccuracy of the classical potential function. Therefore, we did not consider the effect of surface morphology on interfacial coupling, but simply treated interface coupling as an adjustable parameter.

Zoom In Zoom Out Reset image size

Previous works have shown that the phonons are the main heat carries in graphene [32, 33]. According to the kinetical theory [34], the thermal conductivity can be determined by the phonon capacity, the group velocity, and the mean free path. In order to understand how the interlayer coupling strength influences both the phonon group velocity and scattering, the phonon dispersion of an AB stacked two-layer graphene model is calculated along the ΓK direction of the first Brillouin zone (FBZ) with different interlayer coupling strengths. The phonon dispersion can also be calculated using the same method and does not influence the main conclusion below. There are four atoms in the unit cell of a two-layer graphene. Thus, twelve phonon modes exist in FBZ for each wavevector, including eight in-plane phonon modes (labelled by black solid lines in figure 5) and four out-of-plane phonon modes (labelled by solid red lines in figure 5). It is shown that the in-plane phonon dispersion varies slightly with χ = 3 and 10, when compared to that of a normal interlayer coupling χ = 1. This indicates that the in-plane phonon contribution to the thermal conductivity is not influenced by the interlayer coupling effects. However, upon an increase in the interlayer coupling from χ = 1 to χ = 10, figures 5(b), (d), (f) shows that the frequency of the ZA' mode in the zone center of FBZ is significantly enhanced and passes from 2.14 to 6.78 THz. Contrarily, the frequency of the ZA' mode around the boundary of FBZ almost does not increase. On the one hand, the variation in the ZA' phonon dispersion as a function on an increasing interlayer coupling directly reduces the phonon group velocity and thus suppresses the thermal conductivity. On the other hand, the phonon band gap between the ZA' mode and the high optical phonon modes decreases when the ZA' phonon frequency is enhanced. This provides a higher number of scattering channels for the ZA' modes to satisfy the energy and momentum conservation [14, 17, 35]. Moreover, it enhances the phonon scattering and further cuts down the thermal conductivity.

Zoom In Zoom Out Reset image size

The thermal conductivities of both suspended and supported FLG samples have been investigated via NEMD simulations. The results show that the thermal conductivity of a FLG specimen supported on a crystalline silicon substrate is lower than that of a suspended one. Moreover, there are different trends of the thermal conductivity as a function of the number of layers for these two cases. The thermal conductivity of a suspended FLG sample decreases upon an increase in the layer number, whereas in supported FLG samples it takes an increase with the thickness of the specimen. Both values nearly converge above a thickness of 1.68 nm (5 layers). When a FLG is supported on a rough crystalline silicon, its thermal conductivity shows the same rising trend as a function of the number of layers, but with a lower magnitude than that recorded with a smooth substrate. In order to explain the effects of the supported substrate on the in-plane thermal conductivity, the thermal conductivities of both types of FLG on a smooth and a rough silicon surface are calculated and compared by modulating the coupling strength. The results show that the thermal conductivity of FLG decreases upon an increase in the coupling strength, no matter whether the graphene is suspended or supported. Further investigation on the effect of the interlayer coupling on the phonon dispersion relation indicates that increasing the interlayer coupling does not only reduce the phonon group velocity, but also increases the phonon scattering due to the reduction of the phonon band gap. This leads to a decrease in the thermal conductivity as a function of the interlayer coupling. These findings can be used to address several problems in the use of graphene for research and applications, as well as in designing and manufacturing graphene-based devices, when graphene is integrated and interfacial effects on the thermal transport have to be considered.

The authors would like to thank the financial support from the National Natural Science Foundation of China (51975120 and 51575109) and the Fundamental Research Funds for the Central Universities (2242019k1G011).