Effect of grain boundaries on thermal transport in bi-layer graphene nano-ribbons

Graphene and its derivatives, such as few-layer graphene and carbon nanotubes, are regarded as promising materials for various next-generation applications due to their high thermal^1–3 and electrical⁴ conductivities, as well as excellent mechanical stability.⁵ These applications include high-performance electronics,⁶ sensors,⁷ thermal management,^8,9 energy storage,¹⁰ and quantum computing.¹¹ The 2D nature of graphene and its compatibility with other materials make it highly suitable for heterostructures,^12,13 including the usage as a capping layer on metalization materials at the break-end-of-line¹⁴ to prevent electromigration and rapid heat dissipation. In addition, stacking multiple layers of graphene sheets provides further opportunities to tune their properties.^15,16 One prominent example is the discovery of superconductivity in bi-layer graphene with a twist angle of 1.1°.¹⁶ Therefore, understanding the role of interlayer interactions in layered 2D structures, including bi-layer graphene, is crucial.

Despite its outstanding properties and promising applications, defects and grain boundaries (GBs) commonly occur during graphene growth and processing.^17–20 These defects and GBs are known to significantly impact graphene’s performance in electronics, especially when used for heat dissipation.²¹ Randomly distributed defects typically lower the thermal conductivity, κ, as defect concentration increases,^22–24 following the classical theories of Klemens.²⁵ Conversely, GBs induce an anisotropic thermal transport behavior. While κ can be effectively reduced both across^26–28 and along²⁹ GBs, the reduction mechanisms differ.²⁹ In graphene, these GBs are commonly formed by pentagon–heptagon (5–7) rings.^30–32 It has been previously studied that κ along such GBs in graphene can be reduced by at least 60% due to a combination of specular and anharmonic scattering.²⁹ To date, most theoretical investigations on the effect of defects and GBs on thermal transport have focused on monolayer graphene, despite the versatile applications of multilayer graphene. Experimental characterization has revealed unique defect profiles in different layers of graphene with two or more layers,³³ similar to those observed in multilayer transition-metal dichalcogenides.^34,35 This can lead to various topological features due to different GB-induced symmetries.³⁵ However, the impact of monolayer GBs on the overall thermal transport in multilayer graphene remains unclear.

In this study, we investigated the thermal transport along a string of periodic 5–7 GBs in one of the two layers of bi-layer graphene nanoribbons (GNRs) as an example to understand the impact of monolayer defects and GBs on overall in-plane thermal properties. Using non-equilibrium molecular dynamics (NEMD) simulations, we found that while the GBs in one of the bi-layer GNR sheets decrease κ, the adjacent pristine layer mitigates this reduction. In addition, we observed a slight κ reduction in the pristine layer of the bi-layer defective GNR (d-GNR) due to stronger flexural phonon mode scattering, which is offset by enhanced flexural acoustic modes acting as the main heat carriers. In contrast, in the defective layer, the major phonon scattering events occur in the longitudinal acoustic phonon modes. These different mechanisms for reducing κ in different layers provide insights into controlling thermal transport in graphene structures and other 2D materials.

The setup for calculating κ along the 5–7 GB in the bi-layer GNR using NEMD is presented in Fig. 1(a). To prevent the GNR from rolling and sliding, two unit cells at each end are fixed. The hot and cold baths are formed by five adjacent unit cells with temperatures of T_h = 310 K and T_c = 290 K, respectively. The GB, consisting of a string of 5–7 defect rings, is periodically aligned with two rows of hexagonal unit cells in between in one of the two GNR sheets. This configuration is referred to as a bi-layer defective GNR (d-GNR). We prepared three widths (Ws) of bi-layer d-GNRs: 3.3, 7.3, and 10.3 nm. As a reference, bi-layer pristine GNRs (p-GNRs) of equal widths are also created with a similar NEMD setup, as shown in Fig. 1(a), maintaining an average interlayer spacing of 3.4 Å, consistent with the atomic thickness of graphene.

Before fixing the two ends, a conjugate gradient relaxation was conducted at 0 K to relax the structure. Unlike the monolayer GNR, where a mirror symmetry is preserved about the GB,²⁹ the relaxed bi-layer GNR reveals that the GB breaks the symmetry. In particular, on the left side of the GB, the two graphene layers maintain an AA stacking pattern, while the right side of the GB presents a moiré-like stacking. The ripples in the defective layer caused by 5–7 rings create a larger average interlayer spacing of up to 3.93 Å, as shown in Figs. 1(b)–1(d). Due to the van der Waals (vdW) interaction between the two GNR sheets, the pristine layer also exhibits a bending conformation directly underneath the 5–7 defects. Consequently, both layers maintain similar bending and rippling profiles despite the absence of defects in the pristine layer. It is also interesting to note that in both sheets on the left side of the GB, less rippling is observed since the stacking is mostly maintained as AA stacking, while on the right side of the GB, where a moiré-like pattern is formed, the sheets are periodically rippled along the GB direction, i.e., the y-axis, slightly increasing the average interlayer distance.

We consider one atomic layer thickness of graphene, i.e., 3.4 Å, in the calculation of A and 6.8 Å when considering the overall heat flux of bi-layer GNRs.

To compare κ in the defective and pristine layers of the bi-layer d-GNR, we calculated heat fluxes and temperature profiles for both layers. As shown in the example of a 3-nm wide and 50-nm long bi-layer d-GNR in Fig. 2, the defective layer conducts ∼50% less heat than the pristine layer despite having about a 27% higher temperature gradient. This comparison between the two GNR sheets within the same bi-layer d-GNR highlights the significant impedance of in-plane thermal transport due to GBs. This impedance can be attributed to strong in-plane phonon localization²⁹ induced by the GB, a mechanism similar to that reported in our previous studies of κ along screw dislocations in nanowires,^36,37 which is unaccounted for in classical Klemens theory.²¹ 

To understand the effect of GBs in one of the bi-layer d-GNRs on the overall thermal properties along the GB direction, we calculated the κ of bi-layer p-GNRs of the same widths. To capture both ballistic and diffusive transport, we calculated the κ of GNRs with lengths ranging from 50 to 700 nm, as shown in Fig. 3(a). As expected, due to the significant reduction in κ of the defective layer, the overall κ of the bi-layer d-GNR is consistently lower than that of its pristine counterpart across all lengths. Moreover, for all three widths of bi-layer GNRs, the κ reduction remains around 40% and 44% in the ballistic and diffusive regimes, respectively, which is less severe than the 60% reduction induced by the same type of GBs in monolayer graphene.²⁹ 

Although the pristine layer in the d-GNR mitigates the reduction in κ, it is unclear whether this effect is solely due to the pristine layer averaging out the thermal transport impedance of the defective layer or if interlayer coupling also plays a role. As shown in Figs. 1(b)–1(d), the pristine layer slightly deforms along with the defective layer due to vdW interaction. To investigate this, we compared the κ of the pristine and defective layers in bi-layer d-GNRs with that of p-GNRs, as presented in Fig. 3(b). Interestingly, both layers in the d-GNR exhibit a reduction in κ compared to those in p-GNRs, with a 20% reduction in the pristine layer and a 60% reduction in the defective layer. Thus, the overall κ reduction in bi-layer d-GNRs comprises the effect from both the GB in the defective layer and the interlayer coupling.

To understand the mechanism of the κ reduction in the bi-layer d-GNR, we calculated the total phonon density of states (PDOS) for these structures [Fig. 4(a)]. We observed a significant decrease and broadening of the G-mode peak around 49 THz in the d-GNR, as described using the re-optimized Tersoff potential,³⁸ indicating strain in the C–C bonds.³⁹ This broadening of the G-peak facilitates more Umklapp scattering by providing additional acoustic–acoustic–optical scattering channels for the heat-carrying acoustic modes.⁴⁰ 

A closer examination of different phonon modes in each individual layer of p- and d-GNRs presented in Figs. 4(b)–4(d) reveals that the most significant changes occur in the defective layer of d-GNRs. This layer shows a lowered and broadened G-peak at around 49 THz, with additional modes appearing around 52 THz, particularly in the two in-plane modes [Figs. 4(c) and 4(d)], indicating compression of the C–C bonds due to the presence of 5–7 rings. Changes in the pristine layer of d-GNRs around the G-peak are less pronounced, suggesting a minimal variation in bond lengths. Therefore, much of the aforementioned phonon scattering is attributed to the strain induced by the 5–7 rings in the defective layer of d-GNRs.

In the bi-layer d-GNRs, the flexural (Z-) modes in the pristine layer are slightly sharpened, with a 16% increase in the PDOS peak around 13 THz for the flexural acoustic (ZA) mode and a 13% increase around 26 THz for the flexural optical (ZO) mode. These changes in the Z-modes in the pristine layer of d-GNRs could be due to an enlarged interlayer spacing caused by GBs in the defective layer. This may result in a competition between the κ enhancement due to an increased number of ZA modes, which are the main contributors to thermal transport in graphene,⁴¹ and the ZA–ZO phonon scattering, leading to a slight reduction in κ of the pristine layer in the bi-layer d-GNRs. In addition, there is a significant increase (about 56%) in the transverse acoustic (TA) modes around 13 THz in the pristine layer of d-GNRs compared to p-GNRs, although its effect on the overall κ is not as significant as that of the Z-modes.⁴² Little difference is observed in the longitudinal (L-) modes between the pristine layers of d-GNRs and p-GNRs.

Therefore, we can infer that although κ reduction occurs in both layers of the d-GNRs, different mechanisms are responsible for this reduction. In the defective layer, the GBs introduce strongly localized strains that severely scatter phonons, especially the in-plane modes (T- and L-modes). However, the slight κ reduction in the pristine layer is mainly due to ZA–ZO phonon scattering with its reduction being mitigated by the enhanced ZA modes due to increased interlayer spacing.

In our study, using molecular dynamics simulations, we found that thermal transport in bi-layer GNRs containing a GB in one of the layers is significantly hindered. However, the impedance is not as strong as in monolayer graphene²⁹ due to the presence of the pristine layer in the bi-layer GNR. We observed a reduction in κ in both defective and pristine layers of the d-GNR, each due to different mechanisms. By analyzing the phonon density of states (PDOS) of different phonon modes in separate layers, we determined that the reduction in κ in the defective layer is mainly caused by anharmonicity induced by the strains from the 5–7 defects, resulting in strong in-plane acoustic and optical phonon scattering. In contrast, the reduction in κ in the pristine layer is due to enhanced flexural acoustic and optical phonon scattering as a result of increased interlayer spacing. However, the reduction is partly offset by the enhanced flexural acoustic modes, which are the primary heat carriers. Notably, none of these reduction mechanisms have been accounted for by Klemens theory.²¹ Our findings can potentially be applied to engineer thermal transport in a wider range of graphene-based heterostructures¹² or vdW heterostructures.⁴³ Moreover, our results can suggest more engineering pathways for tuning phonon–electron/exciton coupling in exotic electronic materials.

The construction of graphene nanoribbons (GNRs) with grain boundary (GB) defects began with creating a GNR unit cell in VMD and Avogadro. The pentagon–heptagon (5–7) ring was introduced by selectively breaking a bond and adding an atom near the center of the structure, with the pentagon ring on top of the heptagon ring along the y-axis. The structure is then relaxed in Avogadro⁴⁴ with the universal force field. After constructing the defect-containing unit cell, the structure was systematically translated along the y-axis to create ribbons with lengths ranging from 50 to 70 nm and widths of 3.3, 7.3, and 10.3 nm for comparative analysis across three studies.

The non-equilibrium molecular dynamics (NEMD) simulations were carried out using the LAMMPS package (version August 2, 2024)⁴⁵ with a time step of 0.5 fs. We chose the optimized Tersoff potential³⁸ to describe the C–C interactions and the Lennard-Jones (L-J) potential with the energy parameter ɛ = 2.84 meV and the distance parameter σ = 3.4 Å to describe the interlayer vdW interactions. Such a combination of parameters has been proven to be able to accurately represent the configuration of collapsed carbon nanotubes in our previous study.⁴⁶ The cutoff distance was set to 10.2 Å throughout the simulation.⁴⁷ Since different periodicities exist in the pristine and defective layers of bi-layer d-GNRs, the initial energy minimization with the conjugate gradient method to obtain relaxed configurations shown in Fig. 1 is performed with a periodic box much larger than the dimensions of the bi-layer GNRs. We then fixed the two ends for thermal calculations.

For thermal calculations, the system was initially heated up to 300 K with a Nosé–Hoover thermostat for 1 ns (2 × 10⁶ time steps) followed by 0.5 ns of equilibration with the microcanonical ensemble. The two reservoirs were then rescaled at every time step to maintain the two heat baths at T_h = 300 K and T_c = 290 K, respectively, to enable NEMD. The steady state was reached after 1 ns or longer, depending on the sample size. We take the last 2 ns of the steady state run for the calculation of the heat fluxes and temperature profiles as presented in Fig. 2.

This work was partially supported by the Semiconductor Research Corporation (Grant No. SRC-3181.001). T.B. also acknowledges the Materials Science Program Graduate Fellowship at the University of Vermont for partially supporting this project.

The authors have no conflicts to disclose.

Temitope Boriwaye: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (supporting). Jihong Ma: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Resources (lead); Software (lead); Supervision (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (lead).

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