Thermal transport across graphene step junctions

The emergence of two-dimensional (2D) materials has brought new opportunities to explore fundamental physical properties and to exploit these materials for new applications [1]. As the first isolated 2D material [2–4] and due to its extraordinary transport properties [5], graphene has been extensively studied especially for electronic applications. However, the properties of graphene can be altered due to crystal imperfections which appear, for example, during graphene growth by chemical vapor deposition (CVD). One such type are grain boundaries (GBs) [6], i.e. line defects where two graphene grains (of the same thickness) are stitched together. Other defects are graphene junctions (GJs), i.e. the steps between regions with different number of graphene layers, such as monolayer-to-bilayer (1L–2L) junctions.

The properties of GBs are relatively well understood, having been measured electrically [7], thermally [8], and mechanically [9]. For example, GBs reduce the overall electrical [7] and thermal conductivity [10, 11] of graphene due to electron and phonon scattering, respectively. However, GJs have only recently attracted more interest with few experimental studies of their properties in electronics [12], optoelectronics [13], and as p-n junctions [14]. A theoretical study assigned thermal rectification properties to GJs [15], however this has not been examined experimentally. Other simulations also showed that heat transfer at GJs is non-trivial, because in the multilayer region different layers may have different temperatures [16]. Knowledge of heat flow across GJs is important not just fundamentally, but also for practical applications in terms of how they modify the overall thermal conductivity of graphene (as GBs do [10, 17]), or where GJs could act as phonon filters. As an example, electronic devices based on graphene and other 2D materials often contain GJs, but little is known about how their thermal resistance affects the overall device performance [18, 19].

Here, we investigate for the first time the temperature-dependent heat flow across GJs supported on SiO₂ substrates. Our experimental results combined with molecular and lattice dynamics (MD and LD) simulations indicate thermal decoupling between layers caused by a large thermal boundary resistance (TBR). Thus, we establish a microscopic understanding of thermal conduction across GJs and clarify their role in large-area thermal management applications of graphene.

Figure 1(a) illustrates the schematic of the device structure we used to measure the thermal conductance across GJs. Graphene used in this study (see Methods) is mechanically exfoliated onto a SiO₂/Si substrate (supplementary section 1 (stacks.iop.org/TDM/6/011005/mmedia)) and step junctions were identified by optical microscopy, atomic force microscopy (AFM), Raman spectroscopy (see Methods) and were finally confirmed by scanning electron microscopy (SEM) after all measurements were completed. During thermal measurements two parallel metal lines were used as the heater and thermometer, interchangeably [20, 21]. A thin layer of SiO₂ (~40 nm, electron-beam evaporated, see Methods) underneath the metal lines provided electrical isolation from the graphene. Figures 1(b) and (d) show SEM images of the two devices measured, which correspond to 1L–2L and 2L–4L (bilayer to four-layer) GJs, respectively. Figures 1(c) and (e) show Raman spectra obtained on each side of the GJ, determining the number of graphene layers. The Raman spectra do not show discernible D peaks even after patterning the metal lines, confirming relatively defect-free, crystalline graphene regions (supplementary section 2).

Zoom In Zoom Out Reset image size

We performed heat flow measurements from 50 K to 300 K on these GJ samples and on similar control samples without graphene. We also measured heat flow across the GJs in both directions by swapping the heater and sensor, to test for possible asymmetry in the heat flow as a consequence of phonon scattering at the junction, which would lead to thermal rectification for large temperature differentials [15]. The measurements are performed as follows. Current is forced into a metal line, which acts as a heater, while both metal lines are used to sense temperature, setting up a temperature gradient across the GJ. The metal lines are thermo-resistive elements, which allow us to convert measured changes of electrical resistance into variation of the temperature of the sensor, ΔT_S, and heater, ΔT_H, as a function of the heater power P_H (supplementary section 4). We calibrated both metal lines for each sample by monitoring the resistance over a slightly wider temperature range, from 40 K to 310 K, to determine the temperature coefficient of resistance (TCR) and quantify temperature variations (supplementary sections 6 and 7).

Once the temperature difference between the metal lines is known as a function of heater power, the thermal conductance across the junction is obtained by processing the experimental data using a three-dimensional (3D) finite element model (FEM) [20, 22] (see Methods). In this simulation, the graphene channel region between heater and sensor is treated with an effective thickness h  =  0.34n nm, where n  =  2 in both devices because most of the two channels are covered by 2L graphene (see arrows in figures 1(b) and (d)). In other words, the FEM fits the graphene channel with an effective thermal conductivity, k, between heater and sensor. The effective channel thermal conductance is G  =  kh(W/L), where W and L are the graphene channel width and length.

The FEM shown in figure 1(f) accurately replicates the experimental setup taking into account: (i) all geometric dimensions of the metal lines, determined using SEM images (supplementary section 1); (ii) the thickness of the SiO₂ under the graphene from ellipsometry (supplementary section 3) and its temperature-dependent thermal conductivity from measurements of the control sample (supplementary section 5); (iii) the Si thermal conductivity for Si wafers with the same doping density [23] (supplementary section 3). The FEM also includes the effect of TBR at Si–SiO₂ interfaces [20] from the control sample, graphene–SiO₂ [24] and SiO₂-metal [25] interfaces, based on previous measurements of similar samples [20]. Figure 1(f) shows the simulated temperature distribution with current applied through the heater for the 1L–2L junction device. The thermal conductivity k of the graphene channel is varied in the simulation until ΔT_S and ΔT_H versus P_H modeling results match well with the experimental data.

We also measured a control sample without graphene in the channel to validate our method and to obtain the thermal properties of the parallel heat flow path through the contacts, the supporting SiO₂, the SiO₂–Si interface and the Si substrate (supplementary section 5). These thermal properties obtained after processing the experimental data with the FEM show good agreement with well-known data from literature [20, 26, 27] over the full temperature range. Consequently, these data were used as inputs for the FEM simulation of the GJ structures.

Figure 2 shows the experimental heater temperature rise (in red) and sensor temperature rise (in blue) normalized by the heater power, ΔT/P_H, as a function of temperature obtained for the two junctions studied, 1L–2L and 2L–4L. The heat flow was studied in both directions across the GJ to account for possible thermal rectification effects. The uncertainty of ΔT/P_H is ~0.5%–1%, which agrees well with our previous experiments that used similar metal lines [20].

Zoom In Zoom Out Reset image size

Figures 3(a) and (b) show the schematic of the 1L–2L and 2L–4L GJ samples. The rectangular colored sections on top illustrate the thermal conductivities, i.e. 1L (k₁), 2L (k₂) and 4L (k₄) for non-junction regions, while 1L–2L (k_1–2) and 2L–4L (k_2–4) represent the GJ channel region. These are used by the FEM to process the raw experimental data from figure 2, yielding the effective thermal conductivities shown in figures 3(c)–(f). While the effective thermal conductivity of the GJ channel (k_1–2 and k_2–4) is determined from the temperature gradient between heater and sensor, the thermal conductivity of 1L, 2L and 4L is mainly determined from temperature variations only at the heater surroundings. Although not the main topic of this study, these supported 2L and 4L graphene thermal conductivity estimates are among the first of their kind (others being discussed below).

Zoom In Zoom Out Reset image size

Figures 3(c) and (d) display the extracted thermal conductivity of the 1L, 2L and 4L graphene regions, as well as the effective thermal conductivity of the 1L–2L and 2L–4L junctions, in both directions of heat flow. The 1L, 2L and 4L thermal conductivities show similar values over the entire range of temperature. Their room temperature values are ~500 Wm⁻¹ K⁻¹ for 1L, ~400 Wm⁻¹ K⁻¹ for 2L, and ~450 Wm⁻¹ K⁻¹ for 4L graphene, respectively. These are consistent with earlier measurements by Seol et al [28], Sadeghi et al [29], and by Jang et al [21] who found the thermal conductivity of SiO₂-supported 1L, 2L and 4L graphene were ~580, ~600 and ~480 Wm⁻¹ K⁻¹ at room temperature, respectively. To obtain the various thermal conductivities from the FEM fitting, we used the same TBR between graphene and SiO₂ for all layers, following Chen et al [24], but there may be small differences in the TBR that could be behind this small variation. However, our results are in good agreement with values reported by Sadeghi et al [29], which show that the thermal conductivity of SiO₂-supported graphene few-layers remains very similar.

In comparison, figures 3(c) and (d) show that the effective thermal conductivity in the GJ regions, i.e. k_1–2 and k_2–4, is lower than in the graphene layers, i.e. k₁, k₂ and k₄. This difference becomes more evident as the temperature reduces from 300 K to 50 K. Figures 3(e) and (f) show the effective thermal conductance of the GJ regions, calculated by dividing the thermal conductivity with the metal line separation (see Methods). The thermal conductance for 1L–2L varies from 4.8  ±  1.1  ×  10⁷ W m⁻² K⁻¹ to 9.1  ±  1.2  ×  10⁸ W m⁻² K⁻¹ at 50 K and 300 K respectively, while for 2L–4L it varies from 6.1  ±  1.3  ×  10⁷ W m⁻² K⁻¹ to 7.7  ±  1.2  ×  10⁸ W m⁻² K⁻¹ at 50 K and 300 K respectively. Bae et al [20] explained that as we shorten the length of a graphene channel, quasi-ballistic phonon transport effects reduce its thermal conductivity, because the longest phonon mean free paths become limited by the length of the channel. In other words, the graphene thermal conductivity is length-dependent in this sub-micron regime. The thermal conductivity of our GJ samples is consistent with values reported by Bae et al [20] for length-dependent graphene without junctions. Additionally, that the thermal conductance of the 1L–2L and 2L–4L channels is almost identical for both heat flow directions, i.e. k_1–2  ≈  k_2–1 and k_2–4  ≈  k_4–2, indicates no measurable asymmetry in the heat flow or thermal rectification effects on supported graphene at the junction.

To explain the measured thermal conductance of the GJs in both heat flow directions we consider two possible scenarios. The first scenario consists of thermal decoupling between the top and bottom layers of graphene, which could be attributed to the presence of a large TBR between layers. The thermal decoupling between layers would cause the heat to flow only through one layer, i.e. the bottom one, which would result in similar conductance values as the work of Bae et al [20]. Moreover, the large TBR between layers would make phonon scattering at the junction negligible, which would support the idea of a non-asymmetry or thermal rectification effect. The second possible scenario would be a perfect coupling between the top and bottom graphene layers, i.e. very small TBR between layers, which would explain the similarity of the GJs thermal conductance with those shown by Bae et al [20]. However, under these circumstances, we would expect the junction to scatter phonons more efficiently, which might induce some thermal asymmetry across the junction.

To quantitatively understand the phonon physics at the GJ, we performed atomistic molecular dynamics (MD) simulations and LD calculations (see Methods). First, we evaluate a suspended 1L–2L junction by non-equilibrium molecular dynamics (NEMD) simulations [30] as shown in figure 4. The length of the MD models is up to 200 nm, which, although about half the size of the experimental device, still captures its essential physical properties. To produce a stationary heat current (J), the ends of the device are kept at 350 and 250 K (see Methods), respectively, by two Langevin thermostats. If the system displayed thermal rectification its thermal conductance, computed as G  =  J/ΔT, would differ if the heat current went from 1L  →  2L or from 2L  →  1L. Setting up the NEMD simulations we have two options to treat the bilayer side of the GJ: we can either apply the thermostat to both layers, as in [16], or treat only the top layer as a thermal bath. In the first case we find that the thermal conductance is near that of a single graphene layer, too large compared to the experiments (supplementary section 8). Thus, we focus our analysis on the second case. In fact, NEMD simulations show the thermal conductance of the device is the same, within the statistical uncertainty, regardless of the direction of the heat current. Hence, our simulations also confirm that this system does not display thermal rectification.

Zoom In Zoom Out Reset image size

An analysis of the temperature profile at stationary conditions (figure 4) shows that the top and bottom layers of the junction are thermally decoupled, and the main source of TBR is not the step at the junction, but rather the weak coupling between the two stacked graphene layers. Such weak coupling causes a larger temperature difference (ΔT ~ 34 K) between the top and bottom layer of the device, whereas the temperature discontinuity at the step of the junction is only ~5 K. Hence the main resistive process occurs at the interface between the overlapping layers, which is symmetric, thus explaining why no thermal asymmetry or rectification occurs. Even with a very large temperature difference at the two ends of the device (ΔT ~ 450 K), thermal rectification remains negligible (supplementary section 8).

Our experiments and simulations appear at odds with the NEMD results of Zhong et al [15]. In this work, the system is set up such that there is no thermal decoupling between layers in the thermal reservoir, and this effect is not probed in the non-thermostated junction. Hence these former simulations suggest an asymmetric phonon scattering at the junction that depends on the heat flow direction (thermal rectification effect). By comparing our simulations with theirs, we conclude that an apparent thermal rectification could be observed by sampling the system at non-stationary conditions, stemming from poor equilibration of the thermal baths. This is especially a problem for poorly ergodic systems such as graphene and carbon nanotubes [31].

While NEMD sheds light on the microscopic details of heat transport at the GJ, it does not allow a quantitative estimate of the conductance that can be compared to experiments. In fact, due to the classical nature of MD simulations, quantum effects are not taken into account. Considering that the Debye temperature of graphene exceeds 2000 K and experiments are carried out at room temperature and below, quantum effects are expected to play a major role in determining the conductance. Thus, we also calculated the thermal conductance of the 1L–2L junction, treated as an open system, using the elastic scattering kernel method (ESKM) [32]. ESKM is an LD approach equivalent to Green's functions [33], implemented in a scalable code that allows us to compute coherent phonon transport in systems of up to 10⁶ atoms [34]. Thus, we could calculate the thermal conductance of suspended and SiO₂-supported GJs with the same overlap length as in the experiments. LD calculations give the phonon transmission function $\mathcal{T}\left(\omega \right)$ for an open system with semi-infinite thermal reservoirs, resolved by mode frequency and polarization. The thermal conductance is then computed by the Landauer formula [35], integrating $\mathcal{T}\left(\omega \right)$ over all frequencies:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G=\frac{1}{2\pi}\int_{0}^{{{\omega}_{\max}}}{d\omega ~\hslash \omega \mathcal{T}\left(\omega \right)}\frac{\partial {{f}_{{\rm BE}}}(\omega ,T)}{\partial T},\nonumber \end{align} \tag{ 1 }$$

where T is the temperature and ${{f}_{{\rm BE}}}$ is the Bose–Einstein distribution function, accounting for the quantum population of phonons. In this approach we neglect anharmonic phonon-phonon scattering. This assumption is justified a posteriori by comparing the conductance of a suspended device with overlap length of 25 nm, computed by NEMD, G  =  1.16  ±  0.09  ×  10⁹ W m⁻² K⁻¹, with that obtained by LD using a classical phonon distribution function, G  =  0.92  ×  10⁹ W m⁻² K⁻¹. A ~20% difference between LD and NEMD calculations of G is acceptable, as it may stem not only from neglecting anharmonic scattering in LD, but also from the finite $\Delta T$ in NEMD.

Figure 5 displays the thermal conductance of the 1L–2L GJ calculated by LD as a function of the length of the bilayer part (a) and of the temperature (b), compared to experimental data. To assess the effect of the substrate in the experimental device, we consider models of GJ both suspended and supported on a SiO₂ substrate. The geometry of the suspended model is the same as the one used in NEMD (figure 4). G is independent of the length of the monolayer part of the device, as in this approach it conducts heat ballistically. G is normalized by the width of the GJ and by a nominal thickness of the bilayer part of 0.67 nm, which is the same convention used in processing the experimental data.

Zoom In Zoom Out Reset image size

The agreement between modeling and experiments is excellent at low temperature in figure 5. In the experimental device at higher temperature, heat transfer is still mainly dictated by the TBR between the two graphene layers, but the thermal bath also affects the bottom layer in the bilayer part of the device, thus making the conductance larger than that predicted by the model. The thermal conductance of the device increases with the interlayer overlapping surface area, which is determined by the length of the bilayer part (figure 5(a)). However, G does not grow linearly with the overlap surface and tends to saturate with the overlapping length. The conductance limit of this device is indeed dictated by the ballistic limit of a single graphene sheet [36].

The interaction with the SiO₂ substrate reduces the overall conductance of the device by about 30% at room temperature. In order to achieve quantitative agreement between theory and experiments, it is important to consider the conductance reduction in the model for supported structures. The temperature dependence of the GJ thermal conductance can be almost entirely ascribed to the quantum population of the phonon modes. In fact, figure 5(b) shows that theory and experiments display an excellent agreement at low temperature, while systematic deviations appear at T  >  200 K, allowing us to pinpoint the effect of anharmonic scattering, which is not taken into account in the calculations. Resolving the transmission function by mode polarization shows that only out-of-plane modes contribute to heat transport across the interlayer junction, consistent with another recent study [37] (supplementary section 9). We also observe that the interaction with the substrate causes an offset of the out-of-plane modes of the bottom graphene layer with respect to those of the top layer, thus hampering the transmission function even further and reducing the conductance of the device.

In conclusion, we have experimentally measured, for the first time, the temperature-dependent heat flow across GJs, i.e. 1L–2L and 2L–4L GJs, supported on SiO₂ substrates. MD and LD simulations were used to analyze the GJ thermal transport. The simulations show that the top and bottom layers of the junction are only weakly thermally coupled, and the main source of TBR is not the step at the junction, but rather the weak coupling between the two layers in bilayer graphene. The interaction with the substrate was observed to have a significant effect to achieve good agreement between the theory and experiments. In fact, the values obtained for the experimental and theoretical thermal conductance of supported GJs showed excellent agreement at low temperature (T  <  200 K), whose dependence can be almost entirely ascribed to the quantum population of the phonon modes. The deviations observed above 200 K, allowed us to quantify the effect of anharmonic scattering. Additionally, the thermal decoupling observed between layers suppresses the possibility of thermal rectification in GJs. Our findings shed new light on thermal transport across GJs, revealing thermal decoupling between layers that is behind the large TBR observed. These results also imply that the presence of GJs in large-area (e.g. CVD-grown) graphene should not affect the overall thermal conductivity of the material, unlike GB defects. Thus, the thermal properties of CVD-grown graphene are not expected to be affected by the presence of small bilayer islands, because most heat will be carried in the bottom layer.

4.1. Methods

We acknowledge the Stanford Nanofabrication Facility (SNF) and Stanford Nano Shared Facilities (SNSF) for enabling device fabrication and measurements. We acknowledge Woosung Park for his help setting the cryostat. This work was supported by the National Science Foundation Engineering Research Center for Power Optimization of Electro Thermal Systems (POETS) with cooperative agreement EEC-1449548, by the National Science Foundation (NSF) EFRI 2-DARE grant 1542883, the Air Force Office of Scientific Research (AFOSR) grant FA9550-14-1-0251, the National Research Foundation of Korea grant SRC2016R1A5A1008184, and in part by the Stanford SystemX Alliance.

Supplementary Information is available in the online version of the paper.

The authors declare no competing financial interests.