Graphene (1, 2) is a single-atom-thick two-dimensional material, which has attracted immense interests in recent years because of its extensive fascinating properties. The ultra-high thermal conductivity of graphene (3–5) has raised the exciting prospect of using them for phononics devices (6–9). Moreover, rapidly increasing power density in electronic devices made efficient heat removal a crucial issue for progress in information technology. In addition to electronic and photonic devices, thermal management is also required for energy devices, such as solar cells. The increasing interest to the thermal conductivity of graphene was driven by the application for heat dissipation. It is a strong motivation to study thermal property of graphene and graphene nanoribbons (GNRs) from the viewpoints of fundamental physics and practical application.

Room temperature thermal conductivity of suspended single-layer graphene (SLG) exfoliated from highly oriented pyrolytic graphite was measured to reach a high value of 4800–5300 W/m-K [10]. This high thermal conductivity was confirmed in both exfoliated and chemical vapor deposition grown graphenes (11–13). However, thermal conductivity in fully supported graphene is much lower. The room temperature thermal conductivity of SLG supported on silicon dioxide substrate revealed a value about 600 W/m-K, which may be due to backside scatterings and the flexural phonons leakage into the substrate (14). Wang et al. measured thermal conductance of both suspended and supported few-layer graphene (FLG) by using typical thermal-bridge configuration, and the impacts of sample size and substrate have been reported (15).

There are rich physical phenomena about thermal property of GNRs. The presence of a substrate (16–18) and layer-layer interaction (19, 20) can lower the thermal conductivity, which was attributed to the graphene-substrate and interlayer coupling. Moreover, the effects of edge, roughness, strain, random defect and hydrogen termination have been investigated by different simulation methods (21–28). For instance, with only 10% isotopic doping, the thermal conductivity of GNR is reduced up to 50% (24). With 20% hydrogen coverage, thermal conductivity is only about 30% of that of pristine GNR (25). Moreover, compare to the random doped GNR, much stronger dependence of the thermal conductivity on the isotope concentration for the superlattice structures was observed, with 30% reduction in the thermal conductivity at 50% of the isotope concentration (26). Obviously, almost all of these impacts reduce the thermal conductivity of graphene. As the high thermal conductivity is one of the main advantages of graphene for thermal management device, all the above factors are “negative” impacts. A natural question immediately arises: can we enhance the extremely high thermal conductivity of GNR further? Obviously, exploring approach to enhance thermal conduction of GNR is of both fundamental significance and technological interest.

Moreover, in existing theoretical studies, thermal conductivity of GNR is calculated from the Fourier's law: J = −κ∇T, where J is defined as the heat energy transported along the GNR in unit time through unit cross sectional area, and ∇T is the temperature gradient. Although thermal conductivity of GNRs has been studied theoretically, only average heat flux is used in the calculation. So far it is not clear does heat flux distribute uniformly in the cross section plane of GNR? And on the cross section, which section is with higher capability for heat transfer? Here heat current is the rate of heat energy transfer through a given surface, the SI unit is watt [W]. Heat flux is the heat current per unit area, the SI unit is [W/m²]. In this letter, we will investigate numerically the position dependent local heat flux distribution in cross section plane of perfect GNRs. Then we propose to enhance the ultra-high thermal conduction of GNR by introducing a small gap at the center, i.e., construct a “comb” GNR structure. Our numerical results demonstrate that 0.48 nm gap in cross sectional area can induce 23% increase of average heat flux at room temperature. The surprising thermal conduction enhancement we discovered can be understood from the analysis of phonon scattering mechanism. Our study extends the understanding of heat transport characteristic in GNRs, and may inspire experimentalists to develop new approach to enhance thermal conduction of graphene ribbons.

In our simulations, classical molecular dynamics (MD) method is adopted, by using the LAMMPS simulation package. (29) The advantage of the MD method is that it does not rely on any thermodynamic-limit assumptions and is thus applicable for the study of realistic nanoscale systems. The carbon-carbon interactions are described by using adaptive intermolecular reactive empirical bond order (AIREBO) potential (30), which depends not only on the distance between atoms, but also on the bounding environment around atoms, so they implicitly contain many-body information. The accuracy of this potential has been demonstrated in MD simulations in many carbon-based systems, such as carbon nanotubes and graphene for studying thermal as well as mechanical properties. (25) Figure 1 shows the structure of the armchair graphene nanoribbon used in this work. The length is 9.9 nm and width various from 0.97 to 7.3 nm. Following the conventional notation (21), N-GNR denotes a GNR with N carbon-dimer lines across the ribbon width. Before non-equilibrium MD simulation, the canonical ensemble MD simulation with Nosé–Hoover heat bath (31, 32) first runs for 10⁷ steps to equilibrate the whole system at room temperature. After structure relaxation, fixed boundary condition is used at the two ends of the length (x) direction (Figure 1). Next to the fixed boundary, Nosé–Hoover heat baths with different temperature are applied to the two ends of x direction to establish a temperature gradient along the longitudinal direction. At each end, 4-layers of Nosé–Hoover heat bath are used to avoid any artificial effect induced by the localized edge modes. (33) Free boundary condition is used in the width (y) direction and out-of-plane (z) direction. The equations of motion for the atom i in heat baths are: (1)p→˙i=f→i−γp→i;γ˙=1Q(p→i·p→imi−NkBT);Q=NkBTτ2

Here the subscript i runs over all the atoms in the thermostat, m_i and p→i are the mass and momentum of the i-th atom, f→i is the total force acting on the i-th atom, Nis the total number of degrees of freedom of the atoms in the thermostat, k_B is the Boltzmann constant, γ and τ are the dynamic parameter and relaxation time of the thermostat. Tis the desired temperature of the thermostats. T_i(t) is the instant temperature of atom i calculated from the kinetic energy Ti(t)=mi3kB∑αvα2(t), where αdenotes Cartesian coordinate and v_α (t) is the time-dependent velocity. Due to the high Debye temperature in graphene, quantum correction is applied to MD calculated temperature (19). The velocity Verlet algorithm is used to integrate the differential equations of motion. The non-equilibrium MD simulations are performed long enough (10⁸ time steps, time step is 0.5 fs) to allow the system to reach the non-equilibrium steady state where the temperature gradient is well established and the heat current passing through the system is time independent.

To explore the spacial distribution of local heat flux, we refer to the different atom lines in the GNR width direction as I = 0, 1 … n lines, with N = 2 × n + 1. Here “0” line corresponds to the central dimer of GNR, and “n” line is the edge dimer. Figure 2 shows the local heat flux distribution along the cross section of 21-GNR. Using the local heat flux along the center of GNR (line 0) as reference, we show the normalized local heat flux here. It is interesting the local heat flux along the center of GNR is obviously larger than those close to the GNR boundary. For local heat flux at the boundary, it is only about 35% of that at the center. This demonstrates heat flux distribution across GNR cross section is not uniform. In non-equilibrium steady state, the temperature in the same layer (atoms in the same layer means they have the same x coordinate) is the same, so different dimer line is with the same temperature gradient. Thus in GNR, the thermal transport capability of the boundary atoms is much lower than those at the centers. This can be understood from the energy distribution of phonon localized modes. It has been shown that in GNR (34), the intensity of localized modes is almost zero in the center, while with finite value at the boundary. Compared with two-dimensional sheet, there are edge scattering-induced localized phonon modes in the low frequency region in nanoribbons. The edge-localized phonons are responsible for the simulation results that the local heat flux along the center of GNR is larger than those close to the boundary.

This position dependent local heat flux also provides idea in manipulating thermal conduction by other factors, such as random doping. Thermal conduction is obviously less sensitive to edge doping than doping at the center, simply because the local heat flux at the edge is much lower than that at the center. To control thermal conduction, phononic engineering at the center of GNR is more efficient.

To avoid artifact caused by different empirical potential used, we also calculated the heat flux distribution with Tersoff potential (35). As shown in Figure 2, the results by Tersoff potential coincide with those by AIREBO potential indicating that our results are independent of the empirical potential used. In the next, all the results are calculated by using AIREBO potential.

Next we increase the width of the GNR, change N from 9 to 61. As shown in Figure 3A, for all the width range, the local heat flux in the GNR is confined to the central region. Thus the observed spacial dependent local heat flux is a general characteristic in GNRs. Figure 3B shows the ratio of the edge heat flux to the average heat flux decreases with the GNR width. With width increases, the phonon boundary scattering is suppressed, leads to increases in average heat flux across the GNR. Thus compare to the average heat flux, the contribution from edge to total heat flux decreases.

Furthermore, it is interesting to observe the central heat flux decreases with GNR width. With the increase of width, more and more phonons are excited, which results in the increase of thermal conductivity. Moreover, when the width increases, thermal conductivity increases as a consequence of suppressed phonon scattering at edge. On the other hand, the increase of phonon modes will also increase phonon-phonon scattering that in turn will increase the thermal resistance, thus suppress the energy transport. The thermal conduction is determined by these effects that compete with each other. In narrow GNR, the central atoms are distanced from the edge, and the less phonon-phonon scattering leads to the higher central local heat flux than that in wide GNRs.

Due to the high local heat flux at center, it is interesting to find the total heat current in 21-GNR is 0.07 (in arb. unit), which is higher than half of total heat current (0.11) in 41-GNR. It is obvious that the total heat current in two 21-GNRs is larger than that of one 41-GNR, although they have the similar entire width. The total heat current in two 21-GNRs is about 27% higher than that in one 41-GNR. The enhancement of total heat current by a factor of 27% can promote performance as heat dissipation devices. As in steady state, different dimer line is with the same temperature gradient, the enhancement in total heat current corresponds to increase in thermal conductivity. This provides an idea to enhance the thermal conduction of GNRs. For a GNR consists two narrow branches, its total heat current is hoped to be higher than that in a whole GNR with the same entire width. This idea is verified in 21-GNR and 41-GNR. Inspired by this characteristic, next we explore how to increase total heat current (thermal conduction) of GNR by constructing some spectral structures.

To directly confirm the above mechanism to increase heat current, we “cut” a whole GNR by introducing a small gap at the center, i.e., construct a “comb” GNR structure. A typical comb GNR and the corresponding perfect GNR are shown in Figure 4. A 0.48 nm gap is used to eliminate the interaction between these two branches. The length of the whole GNR is 20 nm, and its width changes from 4.3 to 9.2 nm. The width of each branch in the comb GNR is W', with W = 2 × W' + Δ. Then the heat current in comb GNR is compared with that in whole GNR. Here we introduce the thermal enhancement ratio R to describe quantitatively the heat flux enhancement efficiency, R=JJ0, where J₀ is the average heat flux in the original GNR, and J is the average heat flux in the comb GNR.

As shown in Figure 5, there is obvious enhancement in average heat flux by introducing a small gap at the center of GNR. For comb GNR, the enhancement in average heat flux is due to less phonon-phonon scattering. On the other hand, the small gap will also increase phonon-boundary scattering that in turn will increase the thermal resistance, thus suppress the energy transport. The thermal conduction is determined by these effects that compete with each other. In narrow GNR branch, boundary scattering is the dominant process to limit thermal conduction. At wide range phonon-phonon scattering is the limiting mechanism. As a result, we conclude that the average heat flux enhancement arise primarily due to the phonon-phonon scattering in the narrow GNR is less than the effect in the wide ones. Moreover, this enhancement ratio depends on GNR width remarkably. We change the width of GNR and find the enhancement ratio reaches maximum when width W = 6.3 nm. For the comb GNR, with width increases, more phonons are excited in each branch, and the strong phonon-phonon scattering leads to weaken the enhancement ratio. On the other side, with width decreases, the impact of edge scattering increases, and also decreases the enhancement ratio. Thus there is an optimal width to maximize this enhancement effect. For a GNR with width of 6.3 nm, a 0.48 nm wide gap can enhance its average heat flux of 23%. As the temperature difference of the comb GNR is the same as that in the original GNR, it is reasonable to predict corresponding 23% enhancement in the thermal conductivity. One should note here that comb GNR has a smaller cross sectional area than the original GNR does, therefore the increase in total heat current (heat flux times cross sectional area) is about 16%. Recently, several experiments (36, 37) report the possibilities of cutting GNRs with controlled edge structures. Thus the comb GNR structure proposed in our work can be realized experimentally.

In summary, for the first time, we have demonstrated that Graphene nanoribbons with small gap at the center are promising high thermal conductivity materials. The proposed “comb” structure can increase average heat flux of GNR and the optimal GNR structure parameters are discussed. For GNR with width of 6.3 nm, a 0.48 nm-wide gap in the center can enhance the average heat flux (thermal conductivity) up to 23% at room temperature, and the corresponding increase in total heat current is 16%. The reported impacts can be well understood from the analysis of phonon edge scattering. Furthermore, we have also studied the spacial distribution of local heat flux in the cross section of GNR. Due to the edge scattering, the local heat flux at the edge is only 35% of that at the center of GNR. The ultra-high thermal conductivity of suspended graphene has raised the exciting prospect that the high efficiency heat dissipation device can be realized with graphene. Combined with the availability of experimental fabrication technology, our results shed light on heat conduction enhancement in graphene nanoribbons and may favor graphene applications in thermal management.