Thermal rectification in asymmetric graphene nanoribbons

Introduction

A thermal diode conducts heat more readily in one direction than the other. The basic idea was formalised in nonlinear lattice models and has since expanded into thermal regulators, switches, and phononic devices [1,2]. For a useful nanoscale thermal diode, high rectification alone is not enough: the device must also retain adequate forward heat current, operate at realistic temperature bias, and tolerate fabrication disorder.

Graphene is an attractive but difficult platform for thermal rectification. It combines exceptional in-plane thermal conductivity, long phonon mean free paths, and strong sensitivity to edges, strain, isotope disorder, and substrate coupling [3,4,5,6,7,8,9]. Graphene nanoribbons add lateral confinement and edge scattering, which can lower thermal conductivity and reshape phonon spectra relative to extended graphene sheets [10,11,12]. These same features can create directional heat-flow asymmetry if the ribbon lacks inversion symmetry.

Early molecular-dynamics studies reported thermal rectification in graphene nanoribbons and asymmetric graphene ribbons, linking the effect to geometry-dependent phonon spectra and temperature-dependent mode overlap [14,15]. Later work explored thickness asymmetry, defects, hydrogenation, and ballistic multi-terminal junctions as additional routes to rectification [16,17,18,19,20]. The reported mechanisms are physically plausible, but the magnitude of rectification varies widely because simulations use different ribbon lengths, edge terminations, thermostat schemes, potentials, temperature biases, and definitions of rectification.

This article revisits asymmetric graphene nanoribbon rectification with a controlled simulation matrix. The aim is not to report a record rectification ratio. It is to separate geometric tapering, edge roughness, and chemical grading; quantify the conductance penalty that accompanies rectification; and identify which phonon bands are responsible for directional heat-current asymmetry under realistic nanoribbon dimensions.

Nanoribbon geometries

The base structures are suspended single-layer graphene nanoribbons with either armchair or zigzag long edges. Ribbons are hydrogen-passivated at open edges unless otherwise noted. Four length families were generated: 60, 100, 160, and 240 nm. The minimum width was varied from 3.2 to 8.5 nm, and the maximum width from 7.8 to 18.6 nm. Width-tapered ribbons used a linear edge taper over the central 80% of the device, with short rectangular contact regions at both ends to reduce thermostat artefacts.

Three asymmetry classes were studied. The first is geometric tapering: a narrow hot end and wide cold end define the forward direction, while reversing thermostats defines the reverse direction. The second adds edge roughness by randomly removing edge dimers with a prescribed correlation length and then repassivating dangling bonds. Roughness is applied either symmetrically, concentrated at the narrow end, or concentrated at the wide end. The third class uses chemical grading through partial hydrogenation of one side of the ribbon, ramped from 0 to 18% coverage along the transport direction.

The study intentionally avoids extreme geometries where the narrow section becomes a mechanically fragile molecular wire. The minimum neck is always wider than 3 nm, and no structure contains disconnected carbon islands. After geometry generation, each ribbon is relaxed at 0 K and then equilibrated at the target mean temperature for 2 ns. Structures with out-of-plane buckling above 0.9 nm after equilibration are discarded because their thermal current becomes dominated by geometric instability rather than by controlled rectification.

The selected geometry range overlaps the scale where experiments and simulations show strong size effects in graphene thermal transport [8,10,11,12]. It is still smaller than many practical lithographic devices, so the conclusions should be read as nanoribbon-scale mechanisms rather than wafer-scale thermal-management predictions.

Simulation protocol

Non-equilibrium molecular dynamics simulations were performed with fixed end slabs, Langevin reservoirs, and a Newtonian interior. The hot and cold reservoirs each occupied 6 nm of ribbon length. Heat current was measured from the net thermostat energy exchange after the initial transient. Each production run lasted 18 ns, divided into nine 2 ns blocks for uncertainty estimation. For selected cases, reverse non-equilibrium molecular dynamics was also performed as a check on thermostat sensitivity, following the velocity-swap logic introduced by Muller-Plathe [13].

The forward heat-current direction is defined as heat flowing from the narrow end to the wide end for geometrically tapered ribbons. Rectification is reported as R = |J_forward/J_reverse| - 1, with positive R indicating larger heat current in the narrow-to-wide direction. We also report the forward conductance G_forward = J_forward/(A_eff Delta T), where A_eff uses the local ribbon width averaged over the central region and an effective graphene thickness of 0.335 nm. Because this area convention is arbitrary for monolayer graphene, comparisons within the same study are more meaningful than absolute W m^-1 K^-1 values.

Mean temperatures of 250, 300, 400, 500, and 600 K were used with imposed temperature differences of 20, 40, 80, 120, and 160 K. The range deliberately includes biases too large for many devices, because rectification can be nonlinear and may vanish near small bias. Heat-current convergence was checked by comparing block-averaged currents, left and right thermostat exchange, and spatial temperature profiles. Cases with temperature jumps exceeding 18% of the imposed bias at either contact were re-run with longer reservoirs.

Spectral energy density was computed from velocity autocorrelation in 20 spatial bins along each ribbon. The analysis separates in-plane longitudinal, in-plane transverse, and out-of-plane flexural contributions by projecting velocities onto local ribbon axes. Spectral overlap between hot-side and cold-side local densities of states is used as a qualitative indicator, not as a full Landauer transmission calculation. This restraint is important because classical molecular dynamics includes anharmonic scattering but not quantum heat capacity corrections at low temperature.

Thermal conductivity and finite-size checks

Rectification cannot be interpreted without checking baseline thermal transport. Rectangular armchair ribbons with width 8.5 nm showed length-dependent thermal conductivity increasing from 670 +/- 40 W m^-1 K^-1 at 60 nm to 1180 +/- 90 W m^-1 K^-1 at 240 nm at 300 K. Zigzag ribbons were 8-14% more conductive over the same range. These values are below suspended large-area graphene measurements, as expected for narrow ribbons with strong edge scattering, but they reproduce the known trend that graphene thermal conductivity depends strongly on length and boundary scattering [4,5,8,10].

Edge roughness reduced thermal conductivity substantially. A 0.5 nm root-mean-square roughness amplitude lowered conductance by 24% in armchair ribbons and 31% in zigzag ribbons at 160 nm length. This is consistent with theoretical studies showing that rough edges scatter phonons and suppress nanoribbon thermal conductivity [7,11,12]. The conductance penalty matters because a ribbon can show attractive rectification while conducting too little heat for a useful thermal-routing element.

The finite-size dependence of rectification was weaker than that of conductance but still significant. In purely tapered armchair ribbons with taper ratio 3:1 and Delta T = 120 K, R increased from 0.14 at 60 nm to 0.24 at 160 nm, then saturated within uncertainty at 240 nm. Short ribbons over-emphasised contact and ballistic contributions, while very long ribbons reduced asymmetry as diffusive scattering washed out the spectral mismatch. We therefore use the 160 nm family for most mechanism comparisons.

Thermostat choice changed absolute heat currents by up to 7% but changed rectification ratios by less than 0.03 for the main cases. The exception was the narrowest rough-edge ribbons, where local temperature profiles were noisy and contact jumps were large. Those cases are reported in the supplementary data but excluded from the headline comparison.

Geometric rectification

Purely width-tapered ribbons showed moderate rectification. At 400 K mean temperature and Delta T = 120 K, armchair ribbons tapered from 4.1 to 15.7 nm produced R = 0.24 +/- 0.03 and forward conductance of 0.61 nW K^-1. The corresponding zigzag ribbon gave R = 0.19 +/- 0.03 with slightly higher conductance. Rectification increased with taper ratio up to about 3.5:1 and then plateaued because the narrow neck became the dominant resistance in both directions.

Temperature bias was crucial. At Delta T = 20 K, R was below 0.05 for all purely geometric cases. At 80 K, R reached 0.15-0.20 depending on edge orientation, and at 160 K the best geometric case reached 0.29. This nonlinearity is expected for a thermal diode: reversing the bias must change local temperature-dependent phonon populations enough to alter transmission [1,2,15]. It also means that large reported ratios at high bias should not be extrapolated to small-signal thermal management.

Spectral analysis shows that the forward direction preserves overlap between low-frequency flexural modes in the narrow and wide sections better than the reverse direction. In reverse bias, the wide hot end populates a broader set of transverse and edge-localised modes, many of which scatter at the narrowing transition. The bottleneck is therefore spectral and geometric. It is not simply that the narrow end has smaller area; the area resistance is present in both directions, while mode population and anharmonic scattering change with temperature profile.

A short abrupt taper produced larger apparent rectification than a long linear taper in 60 nm devices, but the effect disappeared at longer length and produced large contact jumps. We interpret the abrupt-taper result as a contact-scattering artefact rather than a robust diode mechanism. The more gradual tapers have lower peak R but more stable temperature profiles and better forward conductance.

Defects, roughness, and hydrogenation

Asymmetric roughness increased rectification when the rough segment was placed near the narrow end. For the 160 nm armchair ribbon, adding roughness to the narrow half increased R from 0.24 to 0.31 at Delta T = 120 K, but reduced forward conductance by 38%. Roughness on the wide half produced smaller rectification, R = 0.18, because both bias directions encountered strong scattering before or after the taper. Symmetric roughness lowered conductance without a consistent rectification gain.

Point-defect gradients produced a similar trade-off. A vacancy concentration ramp from 0 to 0.8% along the ribbon raised R to 0.34 +/- 0.04 but cut forward conductance nearly in half. Structural-defect studies have reported enhanced graphene rectification through asymmetric scattering, but our results show that the gain must be evaluated together with current loss [17,19]. For thermal regulation, a high R at vanishing conductance is not useful.

Hydrogenation gradients produced the largest rectification among the tested designs. A ramp from pristine graphene at the narrow end to 18% hydrogen coverage at the wide end gave R = 0.46 +/- 0.05 at 400 K and Delta T = 120 K. The mechanism differs from geometric tapering: hydrogenated regions suppress high-frequency in-plane modes and localise some edge modes, so reverse bias encounters a stronger spectral mismatch. This agrees qualitatively with earlier work on hydrogenation-controlled graphene nanoribbon rectification [18].

The conductance penalty for hydrogenation was severe. Forward conductance fell by 55-62% relative to the pristine tapered ribbon. The hydrogenated devices also showed larger sample-to-sample variation because the exact distribution of C-H sites changed local vibrational spectra. Hydrogenation is therefore a powerful mechanism for studying rectification physics, but not automatically the best engineering design.

Operating map

We summarise the design space with a rectification-conductance map. Pure tapering occupies the high-conductance, moderate-rectification region. Roughness and vacancies move devices toward higher rectification but lower conductance. Hydrogenation reaches the highest rectification but with the largest current penalty. Three-terminal or junction-based graphene concepts can introduce additional asymmetry, but they also introduce contact and fabrication complexity beyond the two-terminal nanoribbons studied here [20].

Mean temperature changes the ranking. At 250 K, rectification is lower because the classical MD model overpopulates some high-frequency modes but still shows reduced anharmonic redistribution; we therefore treat low-temperature numbers qualitatively. At 500-600 K, rectification increases for rough and hydrogenated ribbons because phonon-phonon scattering and mode mismatch become stronger, but edge reconstruction and hydrogen stability become more questionable. The most plausible operating window for the modelled suspended ribbons is 300-450 K.

The best practical design in our matrix is not the maximum-R hydrogenated ribbon. It is a moderately tapered armchair ribbon with controlled roughness at the narrow end, giving R = 0.31 and retaining enough forward conductance for local heat spreading. This design is still far from a macroscopic thermal diode. It would be better described as a nanoscale thermal regulator whose directional response could bias heat flow near hotspots or thermoelectric junctions.

Device integration remains difficult. Substrate coupling can damp flexural modes and reduce the spectral mismatch responsible for rectification. Contacts can dominate thermal resistance in short ribbons. Lithographic edge roughness may be less controlled than the correlated roughness used here. These effects tend to reduce both conductance and reproducibility, so suspended idealised simulations should be treated as upper-bound mechanism studies.

Comparison with previous mechanisms

The geometric rectification ratios reported here are close to earlier molecular-dynamics studies of asymmetric graphene ribbons when similar temperature biases and lengths are used [14,15]. Differences arise mainly from thermostat length, edge passivation, and how the effective cross-section is defined. Thickness-asymmetric ribbons and chemically modified ribbons can give larger R because they introduce stronger spectral asymmetry than width taper alone [16,18].

The results also explain why graphene is both promising and frustrating for thermal diodes. High intrinsic thermal conductivity provides large heat current, but long phonon mean free paths and broad spectra can reduce asymmetry unless the structure strongly filters modes. Defects and functionalisation provide filtering but lower conductance. This is the central engineering compromise for graphene rectifiers.

Our spectral-overlap analysis supports a mode-selective interpretation. Low-frequency flexural modes carry a large fraction of heat in pristine ribbons and are sensitive to boundary conditions, length, and substrate coupling [8,9]. Edge-localised and high-frequency in-plane modes contribute more strongly in narrow or chemically modified regions. Rectification appears when reversing the thermal gradient changes how these mode families are populated and scattered along the asymmetric profile.

The same interpretation cautions against over-reading a single R value. A quoted rectification ratio without temperature bias, mean temperature, length, edge condition, and forward conductance is not enough to compare devices. We recommend reporting an operating map rather than one headline number.

Limitations

The main limitation is the use of classical molecular dynamics. Quantum heat-capacity effects are important below room temperature, and empirical potentials approximate anharmonicity and edge chemistry imperfectly. The low-temperature trends should therefore be read qualitatively. Ab initio force constants or machine-learned potentials could improve spectral accuracy, but at substantially higher computational cost.

The second limitation is idealised geometry. The ribbons are suspended, clean, and contacted by ideal reservoirs. Real devices have polymer residue, substrate coupling, metal contacts, strain gradients, oxidation, and roughness statistics imposed by fabrication rather than by a random generator. These effects can reduce flexural-mode contributions and change rectification. Experimental validation would require suspended devices with well-characterised edge profiles and independent thermal-contact measurements.

The third limitation is timescale. Production runs of tens of nanoseconds are long for atomistic simulation but short compared with rare defect migration or chemical reconstruction. Hydrogenated ribbons at high temperature may evolve chemically in ways not captured by fixed bonding topology. We therefore use hydrogenation as a controlled spectral perturbation, not as a final device prescription.

Finally, the model does not include electron heat conduction. For undoped graphene at the simulated dimensions, phonons dominate heat transport, but electrostatic gating, metallic contacts, or narrow-band electronic states could change the balance. Coupled electron-phonon thermal rectification remains outside the present scope.

Conclusion

Asymmetric graphene nanoribbons can rectify heat flow through a combination of geometry-dependent phonon spectra, edge scattering, and chemically induced mode filtering. Pure width tapering provides moderate rectification while retaining useful conductance. Edge roughness, vacancies, and hydrogenation increase rectification but impose large conductance penalties. The best designs therefore depend on the intended function: maximum heat-current asymmetry, maximum forward heat flow, or robust operation under fabrication variability.

The study reinforces a simple design rule for graphene thermal diodes: report rectification and conductance together. High R in a strongly damaged ribbon may be less useful than moderate R in a conductive tapered ribbon. Future work should combine atomistic simulation with phonon transmission calculations, substrate-coupled device models, and experimental suspended nanoribbon measurements to determine whether the mechanisms identified here survive realistic contacts and fabrication disorder.

Data and code availability

All nanoribbon geometries, random edge seeds, hydrogenation masks, LAMMPS input files, thermostat parameters, raw heat-current logs, block-averaged statistics, spectral analysis scripts, and plotting notebooks are provided in the supplementary archive. Analysis used Python 3.11, NumPy 1.26, SciPy 1.11, MDAnalysis 2.6, OVITO Pro 3.10, and LAMMPS stable release 2 August 2023.

References

1. Li, B., Wang, L. & Casati, G. Thermal diode: rectification of heat flux. Phys. Rev. Lett. 93, 184301 (2004).
2. Wehmeyer, G., Yabuki, T., Monachon, C., Wu, J. & Dames, C. Thermal diodes, regulators, and switches: physical mechanisms and potential applications. Appl. Phys. Rev. 4, 041304 (2017).
3. Geim, A. K. & Novoselov, K. S. The rise of graphene. Nat. Mater. 6, 183-191 (2007).
4. Balandin, A. A. et al. Superior thermal conductivity of single-layer graphene. Nano Lett. 8, 902-907 (2008).
5. Ghosh, S. et al. Extremely high thermal conductivity of graphene: prospects for thermal management applications in nanoelectronic circuits. Appl. Phys. Lett. 92, 151911 (2008).
6. Pop, E., Varshney, V. & Roy, A. K. Thermal properties of graphene: fundamentals and applications. MRS Bull. 37, 1273-1281 (2012).
7. Nika, D. L., Pokatilov, E. P., Askerov, A. S. & Balandin, A. A. Phonon thermal conduction in graphene: role of Umklapp and edge roughness scattering. Phys. Rev. B 79, 155413 (2009).
8. Xu, X. et al. Length-dependent thermal conductivity in suspended single-layer graphene. Nat. Commun. 5, 3689 (2014).
9. Fugallo, G., Cepellotti, A., Paulatto, L., Lazzeri, M., Marzari, N. & Mauri, F. Thermal conductivity of graphene and graphite: collective excitations and mean free paths. Nano Lett. 14, 6109-6114 (2014).
10. Guo, Z., Zhang, D. & Gong, X.-G. Thermal conductivity of graphene nanoribbons. Appl. Phys. Lett. 95, 163103 (2009).
11. Aksamija, Z. & Knezevic, I. Lattice thermal conductivity of graphene nanoribbons: anisotropy and edge roughness scattering. Appl. Phys. Lett. 98, 141919 (2011).
12. Savin, A. V., Kivshar, Y. S. & Hu, B. Suppression of thermal conductivity in graphene nanoribbons with rough edges. Phys. Rev. B 82, 195422 (2010).
13. Muller-Plathe, F. A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity. J. Chem. Phys. 106, 6082-6085 (1997).
14. Hu, J., Ruan, X. & Chen, Y. P. Thermal conductivity and thermal rectification in graphene nanoribbons: a molecular dynamics study. Nano Lett. 9, 2730-2735 (2009).
15. Yang, N., Zhang, G. & Li, B. Thermal rectification in asymmetric graphene ribbons. Appl. Phys. Lett. 95, 033107 (2009).
16. Zhong, W.-R., Huang, W.-H., Deng, X.-R. & Ai, B.-Q. Thermal rectification in thickness-asymmetric graphene nanoribbons. Appl. Phys. Lett. 99, 193104 (2011).
17. Wang, Y., Chen, S. & Ruan, X. Tunable thermal rectification in graphene nanoribbons through defect engineering: a molecular dynamics study. Appl. Phys. Lett. 100, 163101 (2012).
18. Melis, C., Barbarino, G. & Colombo, L. Exploiting hydrogenation for thermal rectification in graphene nanoribbons. Phys. Rev. B 92, 245408 (2015).
19. Nobakht, A. Y. et al. Thermal rectification via asymmetric structural defects in graphene. Carbon 132, 565-572 (2018).
20. Ouyang, T., Chen, Y., Xie, Y., Wei, X. L., Yang, K., Yang, P. & Zhong, J. Ballistic thermal rectification in asymmetric three-terminal graphene nanojunctions. Phys. Rev. B 82, 245403 (2010).