Graphene is currently a prime candidate material for micro-scale and nano-scale heat management applications Cahill et al. (2014) due to its exceptional thermal conductivity Huang et al. (2021); Chen et al. (2020); Pop et al. (2012); Nika and Balandin (2012); Balandin (2011); Li and Zhang (2013). By nanoscale patterning or shaping, the material properties can be tailored to suit particular demands. In the context of thermal transport, this pattering is done to locally modify the phonon dispersion relation Hu et al. (2016) and the corresponding research area, phononics, is sharply on the rise Maldovan (2013); Li et al. (2012). Phononics aims at engineering devices to control transport and detect phonons (heat carriers) in the same way electrons carry information in electronics. Just as the semiconductor diode paved the way for integrated electric circuitry, efforts to construct efficient acoustic and thermal diodes aim at a similar revolution in phononics Maldovan (2013); Zhu et al. (2018); Zhang et al. (2022).

In a thermal diode, or rectifier, the magnitude of the heat flux resulting from a temperature difference applied at the two terminals differ depending on which terminal is the hotter Li et al. (2012); Shih et al. (2015). An important figure of merit is the rectification ratio, R R , defined as the relative difference between these two heat fluxes (currents) or R R = ( J + − J − ) / J − in which J + and J − are left to right (forward) and right to left (reverse) heat fluxes, respectively. However, in contrast to acoustic diode proposals, which mostly rely on single-mode phononic analogs of photonic crystals, the implementation of thermal diodes is not as straightforward. This is because heat is carried by a wide spectrum of phonon modes.

The necessary conditions for a device to display rectification are nonlinearity of the inter-atomic potential and broken mirror symmetry. If this leads to a temperature dependent mismatch of the vibrational density of states (vDOS) between the left and right regions, thermal rectification may ensue. Achieving this condition is easy by attaching two pieces of bulk solid materials with different dependence of their thermal conductivity to the temperature. An example is a junction of bulk pristine crystal to an alloyed or amorphous variant of that in Kobayashi et al. (2009). However, inducing asymmetry and difference in vDOS in the same material needs extra measures. Most proposals relating to carbon nanotubes Wu and Li (2007); Ni et al. (2011), graphene, and other 2D materials have thus so far focused on asymmetry in (a) the number of layers, (b) internal geometry, (c) external or edge shape, (d) material (hetero structure), and (e) asymmetry due to defect i.e., pristine-defective junctions. Implementing heat rectifying devices using the above methods may be very hard to realize, even if it shows promising R R in theory. For a few examples of the above-mentioned devices see Yang et al. (2009); Hu et al. (2017); Hu et al., (2009); Xu et al. (2014); Wang et al. (2014).

Here we show that simply using two weakly connected graphene ribbon segments of unequal length can provide sufficient mismatch of vDOS for significant thermal rectification to emerge. We attribute this to the temperature dependency of the vDOS which emanates from the non-linearity of atomic potentials as well as the length difference which leads to different amount of mechanical strain applied to each section. This will result in different vDOS overlaps for thermal bias in the forward and the backwards directions, which in turn leads to a highly asymmetric flow of heat. Chief among our results is that the rectification ratio and its sign can be adjusted by changing the length of the short and long nano-ribbon sections on-the-fly using electric field-induced strain. An array of electrodes patterned under the free standing graphene controls the point of force application. The point of applying force (strain) acts as the third terminal of the device and causes breaking the symmetry of heat flow. Our proposed principle is akin to applying a normal magnetic field to a two-terminal two-dimensional electron gas (2DEG) device and breaking the time-reversal symmetry and causing non-reciprocity in electric current. The value of R R can reach up to 60% which is comparable or better than some of the proposals based on more elaborate nano-patterning of graphene by defects [See Malik and Fobelets (2022); Zhao et al. (2022)].

The reported values of R R vary among different works. For example, in graphene-based devices this changes from 6% Yousefi et al. (2020a) to 2,580% Yang et al. (2018). Partly, this depends on how the rectification ratio is defined, what is the base temperature, and how much is the applied temperature gradient or bias across the junction. Also some definition of R R are based on the difference of thermal conductivity Δ κ = κ + − κ − rather than the heat flux difference ( Δ J = J + − J − ) . Roberts and Walker (2011), Malik and Fobelets (2022), and Zhao et al. (2022) have reviewed the theoretical and experimental works on heat rectification in solid and 2D materials and listed the achieved R R in detailed tables. While some numbers in these tables show very promising values, it must be remembered, that these predictions assume nano-fabrication with high tolerances, and may thus be sensitive to structural changes on the nanometer scale. Also the reported R R values are not adjustable after the fabrication of the device. Furthermore, we here prefer to use a more conservative definition of R R defined as | J + − J − | / m a x ( J + , J − ) . This makes R R insensitive to direction and confines it to an interval between zero and unity, where R R = 0 means no rectifying behavior and R R = 1 implies completely unidirectional transport.

The remainder of this paper is organized as follows. In the next section the molecular dynamics calculations using LAMMPS package and post-processing steps are explained. To demonstrate that thermal rectification can be achieved without detailed nano-patterning, we consider a simple thermal junction composed of a short and a long graphene nano-ribbon connected via a thermal weak link, and review the basic physical explanation of the origin of this effect in terms of temperature dependent vDOS-overlap. Although we can attribute the effect to the behavior of the overlap of the vDOS, we show that the conventional measure of the overlap is here too simple to quantitatively explain the large observed rectification. Thereafter, we present a simplified system made adjustable using an external force due to applied DC voltage. This showcase the role of the third terminal and the mechanical strain which causes symmetry breaking in analogy to the normal magnetic field applied to 2DEG layer. The success of every proposed scheme to achieve heat rectification depends on how much strain is induced in each section of the heat diode. We proved that strain enhances the Umklapp process or down-conversion of phonons and this modifies the low-frequency portion of the vDOS Shiri and Isacsson (2019).

The molecular dynamics simulation is performed using LAMMPS package Plimpton (1995). The graphene nano ribbons are of rectangular shape in zigzag direction and their widths and lengths are (22 nm and 22 nm) and (10 nm and 100 nm), respectively along x and y directions. The size of the simulation box along the z direction is ± 50 Å to avoid interaction of the graphene ribbon with its periodic replica along the z direction. The ribbon is built with ATOMSK code Hirel (2015) starting from a unit cell made of four carbon atoms cut from the upper half of a carbon hexagonal ring. The metal unit is chosen. The molecular dynamics simulation flow is divided into the following sequential steps.

• energy minimization The energy of the ribbon is minimized under periodic boundary conditions in all directions using Conjugate Gradient (CG) algorithm with box/relax command. The Tersoff three-body potential by Lindsay and Broido (2010) which is modified for heat transport in graphene is used. The integration time step is 0.5 fs.

For visualizing purposes e.g., viewing the output geometry at each time step, The OVITO package provided by Stukowski (2010) is used. The vibrational density of states of atoms within a region is calculated from velocity auto correlation function (VACF) which is implemented in LAMMPS. The VACF for each region is recorded in the output file for each time step (e.g., d t = 0.5 fs). For each component of the atomic velocity there is one column of VACF data which is accessible using [1], [2] and [3] suffixes and [4] is used for the total value of VACF i.e., ⟨ v ( t ) v ( 0 ) ⟩ = ⟨ v x ( t ) v x ( 0 ) ⟩ + ⟨ v y ( t ) v y ( 0 ) ⟩ + ⟨ v z ( t ) v z ( 0 ) ⟩ . Fourier transforming the VACF in MATLAB returns the vibrational DOS spanning from zero to frequencies around 60 THz for the left and the right side respectively. Finally, the overlap S is defined by Lan and Li (2006) as S = ∫ D O S L ω D O S R ω d ω ∫ D O S L ω d ω ∫ D O S R ω d ω (1) where ω is phonon frequency and L and R correspond to Left and Right sections of the thermal junction, respectively. To simulate the effect of force it is possible to apply the force (max 0.05 eV/Å) every few time steps to the chosen area using addforce command and this divides the 10 nm × 100 nm nano ribbon into long and short sections. Alternatively it is possible to add electric charge of one electron to the atoms of the chosen area and apply electric field along z direction using efield command with a value between 0.01–0.05 V/Å.

We consider first a simple junction geometry as shown in Figures 1A, B, consisting of a 22 nm × 22 nm graphene sheet clamped at two ends, interrupted by a thermal weak link made from an equidistantly-spaced array of pinned circular regions. Using molecular dynamics (MD), the average heat flux flowing through the structure when the two clamped ends are at different temperatures was recorded.

Prior to imposing pinned sections on the graphene nano-ribbon, the thermal conductivity of the graphene sheet is calculated using Non-Equilibrium Molecular Dynamics (NEMD) Hu et al. (2016) for calibration purposes. Immediately after and before the left and right clamped ends, hot and cold regions are kept at 320 K and 280 K, respectively. Calculating the heatflux J , the thermal conductivity κ is determined from J i = − κ i j ∂ T ∂ x j , where i and j are coordinate indices as the heat conductivity is in general a tensor.

Values of κ x x = 570 W/m K and κ x x = 360 W/m K are obtained using modified Tersoff Lindsay and Broido (2010) and CH-AIREBO Stuart et al. (2000) potentials, respectively. These values are in close agreement with those obtained using MD for graphene nano-ribbons of similar sizes Hu et al. (2010); Guo et al. (2015). As the length of the ribbon is smaller than the intrinsic phonon mean free path (MFP) in graphene, l M F P = 600 nm Pop et al. (2012), a difference between thermal conductivity of a 22 nm × 22 nm in our model and a 20 nm × 20 nm nano ribbon in Hu et al. (2010) and Guo et al. (2015) is expected. The higher thermal conductivity calculated by using the Tersoff potential is attributed to the stiffer Carbon-Carbon bonding in the Tersoff potential as opposed to the AIREBO potential. Stiffer bonds gives a higher group velocity for acoustic phonons and, consequently, a larger thermal conductivity (Chen et al. (2020)) since k = 1 3 c V v g l M F P , where c V , v g , and l M F P are specific heat, group velocity, and mean free path of phonons in a bulk isotropic solid. Finally, it should be pointed out that classical MD neglects quantum corrections, which can be rather large for graphene due to its large Debye temperature Tewary and Yang (2009); Jiang et al. (2010). However, we are here interested in demonstrating a qualitative phenomenon and derive a ratio between heat-fluxes. Hence, for our purposes, classical MD suffices.

Note that the absolute value of thermal conductivity does not matter for study of heat rectification. What matters is the stress-strain relation and how soon or late the applied force leads to asymmetry in the heat flux. Therefore, we speculate that using CH-AIREBO (soft) potential instead of Tersoff (hard) potential leads to down-shift of frequency contents in VDOS, and it should not change the results qualitatively.

Figure 1C shows the result of simulating the phononic junction of Figures 1A, B where the terminals where kept at temperatures T L = 350 K, T R = 250 K and vice versa. A clear rectifying behavior with heat fluxes J ± obeying J − > J + and a rectification ratio of ≈ 41 % is observed. This rectification ratio corresponds to the 160% ratio of reverse to forward heat fluxes, i.e., J − / J + = 1.6 . Changing the pinning diameters to 20 Å gives similar results. Results for similar simulations at other base temperatures are shown in Table 1. Interestingly, we find values which are strongly temperature dependent and non-monotonic with increasing temperature.

To understand the mechanism of rectification in this system, it is instructive to first recall the well studied model of a nonlinear one dimensional chain of oscillators connected by an interface layer (or interface bond) Li et al. (2004); Terraneo et al. (2002); Lepri et al. (1997); Yang et al. (2007) which is shown in Figure 2A. If the elastic properties on the two sides of the junction differ, by using for instance different spring constants and different strength of nonlinear potential, as a result the corresponding density of vibrational modes will differ. In the absence of inelastic scattering, the vDOS must have an overlap for phonons to cross the junction. Following the discussion of Benenti et al. (2016) and Li et al. (2004) we consider a system which is composed of two sections modeled as atoms connected by springs, and a Hamiltonian H = H L + H R . The left and right sections are weakly linearly coupled with a single link k l i n k yielding a total Hamiltonian H = H L + H R + H i n t . The necessary nonlinearity is introduced through adding a Frenkel-Kontorova (FK) co-sinusoidal potential to each section as shown by Li et al. (2004) which results in H L , R = ∑ i P i 2 2 m i + 1 2 k L , R x i − x i + 1 − a 2 − V L , R 2 π 2 cos 2 π x i .

Here V L , R are the depths of the FK potential, x i and a are atom positions and the lattice constant, respectively.

Rectifying behavior requires a combination of nonlinearity and broken symmetry. Symmetry breaking is obtained by setting V L > V R and k L > k R . This means that the atoms in the left section are connected by stiffer springs and are exposed to a deeper FK potential as opposed to the right section. The FK-potential is a softening potential, i.e., high amplitude vibrations have lower frequencies than low amplitude ones. Hence, at low temperatures they oscillate with higher frequencies than at high temperatures. With increased temperatures the vDOS moves toward lower frequencies [see Figure 2B]. The same trend exists for the right section of the junction however with different rate as well as different frequency. As springs on the right are softer than those on the left, the vDOS on the left is predominantly covering lower frequencies. This explains why for T R > T L the vDOS overlap is smaller than when T R < T L [see Figure 2B].

To enhance the rectification, it is favorable to use hardening and softening potentials on either sides. An example is to replace the FK-potential in the right with a Fermi-Pasta-Ulam (FPU) type system. The nonlinearity in the FPU- β model is an extra anharmonic term like β x 4 in the Hamiltonian with β > 0 . Due to this term, the behavior of the vDOS with temperature is opposite on the right and left sides. Hence, for an FPU- β potential Lan and Li (2006) the increase of temperature shifts the spectrum to higher frequencies according to ω ∝ ( β T ) 1 / 4 . This means that vDOS of the left and right sections of the thermal junction move in opposite directions with increasing temperatures and, consequently, enhance the thermal rectification [see Figure 2C].

We now return to the 2D-system of Figure 1. We consider the area with pinned circles as the weak link between the left and right sections creating the symmetry breaking. To test if the rectification mechanism is similar to the 1D-systems we have studied the vDOS D L , R ( ω ) corresponding to the left and right sections of the nano ribbon at different temperatures (Figure 3).

First we inspect the low-frequency part of D ( ω ) , as shown in Figure 3 for the two sections at T = 300 K in panels (a) and (b). As anticipated, and in contrast to the case of spring-mass models where the spring constant determines the bottom of the band, here the lateral confinement will determine the band bottom. Indeed, regardless of temperature, the left section (short) is gapped from 0 GHz to 58 GHz, while the right (long) section supports pronounced low frequency modes from 10 GHz and upwards. Hence, we can identify the longer (right) side with the soft-spring part and the short side (left) with the stiff-spring part.

Figure 3C shows how the acoustic parts of the spectra shifts with changing temperature for both sides. For these low frequency portions of the spectra, increasing the temperature causes blue-shift to higher frequencies, which resembles the way in which FPU- β model behaves Lan and Li (2006). The nonlinearities responsible for this comes partly from non-linearity of the Tersoff many-body potential Lindsay and Broido (2010) itself. Partly, it can be attributed to the negative thermal contraction of graphene which stiffens the acoustic ZA-modes. This leads to a situation akin to the FK-FK situation in Figure 2B where both spectra shift in the same direction but by unequal amounts. Qualitatively we indeed find agreement that J − > J + for all investigated cases listed in Table 1.

Although giving little contribution to the thermal transport, the high-frequency part of D L , R ( ω ) , corresponding to the optical phonon branches, is shown in Figure 3D. Here the short section supports high frequency modes up to 50.7 THz even at 50 K. However, by increasing the temperature from 300 K to 400 K we observe a gradual shift to lower frequencies e.g., 50 THz. The same trend is true for the long (right) section of the ribbon although the modes have generally lower frequency spanning from 50 THz to 50.4 THz. This is consistent with the observed red shift of the G-peak in Raman measurements, and suggests that increased temperature causes softening of the optical modes and reduction of group velocity due to tensile strain build-up Hong et al. (2016); Wei et al. (2011).

As a measure of the degree of overlap, Li et al. (2005) defined a phenomenological quantity S as the vDOS overlap which was shown in Equation 1 and relates to the rectification ratio by comparing J − / J + to S − / S + . To quantitatively support the above mentioned arguments we have calculated the ratio S − / S + according to Equation 1 for the simulations presented in Table 1. We obtain D L , R from using the velocity auto correlation function (see section 2). Qualitatively we find, in agreement with the simulated heat fluxes, that the overlap is always larger when T R > T L . This is reflected in Table 1 as S − / S + > 1 . However, we find no other quantitative correlations between S − / S + and R R . This is hardly surprising, considering the inclusion and thermal shift of optical modes in the calculation as well as the presence of localized modes in the spectra which may contribute to S but not to J . Hence, we conclude that the overlap S alone is not sufficient to quantify the rectifying behavior.

Based on the above method of dividing a nano-ribbon into a short and long section, we here propose a simpler design with two distinct advantages to other proposals. Instead of partitioning the ribbon in two sections by pinning, we consider applying a force to the nano-ribbon in an asymmetric fashion. The applied force could emanate from sources like mechanical actuation [e.g., piezoelectricity, Atomic Force Microscopy (AFM) tip Hod and Scuseria (2009), Scanning Tunneling Microscopy (STM) tip Neek-Amal et al. (2014)] or electrical force due to an electric field from a local back-gate. The first advantage is that no explicit nano-patterning is needed. The second, and more important advantage, is that changing the position and magnitude of the applied force allow for in situ reconfiguration and control of the rectification.

Specifically we envision a setup here as shown in Figure 4 which is similar to the experimental setup already used to modulate the electron transport properties in carbon nanotubes by Benyamini et al. (2013). Applying electric field to the selected electrode can divide the nano-ribbon into a short and long section. Here a proof of concept molecular dynamics simulation is performed to demonstrate the rectification as a result of force-induced asymmetry.

Figure 5A, shows a simulation set-up of a 10 nm × 100 nm nano-ribbon, in which a force is locally applied along the z -direction to a narrow strip causing deflection. The asymmetric application of the force divides the strip into a long and short section with unequal tensile stresses. The magnitude of forces used in the simulations are of the order nN (nano Newton), characteristic for nano-electromechanical devices. The diode sign in Figure 5 shows the heat rectification direction.

The recorded values of average heat fluxes are shown in Table 2 and they suggest rectification ratios of 35% and 42% for temperature/gradient values of 300   K ± 100   K and 300   K ± 50   K , respectively. If the position of force is changed and it is moved to the right side as in Figure 5B, the same values as above are obtained for the rectification ratios however, with the reversed sign. As for the junction rectifier, we also changed the base temperature to 200 K and 400 K respectively. As can be seen from Table 2, the value of the rectification ratio changes. In case of T L , R = 400   K ± 50   K and T L , R = 200   K ± 50   K , the R R values are 60% and 22%, respectively. The enhanced R R at higher temperature is due to stiffening of acoustic phonon branches which leads to more strain and more VACF asymmetry.

To verify that the rectification is indeed induced from asymmetric length of the nanoribbon section, we also applied the force symmetrically, dividing the ribbon into two equal partitions. The resulting flux differences for the imposed temperature differences of Δ T = 50   K and Δ T = 100   K corresponded to R R = 6 % and R R = 1 %, respectively. This corroborates that the large asymmetry in the values of heat flux indeed stems from the strain-induced asymmetry in the nanoribbon. The significance of this work lies in its flexibility in changing the value and sign of heat rectification without resorting to a newly nano-patterned and/or functionalized graphene. It is also possible to generalize the proposed principle to other two-dimensional materials using not only the electric field but also the magnetic field with proper electrode design e.g., electrodes shaped as current loops under the free standing 2D material.

In the simulation setup the graphene nanoribbon is periodic along the width direction, which means that changing the width of the nanoribbon does not affect the results. In a real device, the change of width leads to a change of thermal conductivity because phonon branches contribute more to heat conduction. However, in this study, the focus is on effect of asymmetry to rectification, and the absolute value of the thermal conductivity does not matter. However, the change of length needs careful study. For a short nanoribbon, the applied force (stress) leads to larger strain value as both asymmetric parts suffer more length change. For a longer nanoribbon, the same force leads to less deflection and less strain as a result. This means that the same rectification ratio requires more applied force (voltage). With respect to the edge type, it was shown that the zigzag edge along the length of the nanoribbon has more stretchability than the armchair edge. For strain values under 20% the stress-strain plot of graphene is the same for both types of edges, and the difference appears for strain values above 30% in which zigzag shows more robustness. Thus, for a device proposed in our work, if the strain is below 20% both edge types perform equally in terms of heat rectification.

Based on molecular dynamics studies of non equilibrium heat transport in graphene nano-ribbons we showed that the thermal rectification ratio can be simply adjusted both in value and sign using a force to divide the graphene into strained long and short sections. Starting from a simple structure in which long and short sections are weakly linked by pinned areas, we showed how this partition lead to the necessary temperature dependent asymmetric vibrational density of states required for thermal rectification.

The importance of this observation and the proposal of Figure 4 is reflected in the simplicity in implementing an in situ tunable heat diode, as exemplified by using an array of local back gates which can selectively induce asymmetric stress in a suspended graphene membrane. Our findings show that the values of heat flux and the rectification ratio can be controlled on-the-fly while the device is operating. Consequently, to achieve rectification it is not necessary to fabricate elaborate asymmetric shapes, i.e., patterning arrays of holes, defects, inducing grain boundaries or the requirement of making a hetero-structure Hong et al. (2016) of pristine graphene with other kinds of graphene [multi-layer Zhang and Zhang (2011); Xu et al. (2014), defected Wang et al. (2012), patterned Yang et al. (2009); Hu et al. (2009); Hu et al., (2017); Wang et al. (2014)], etc. Thus, the combination of simplicity of implementation, and tunability of rectification ratio and sign, opens up for more controlled experiments and developments in the field of 2D-phononics.

Fabricating of gate arrays is a standard process for study of spin-based qubits in semiconducting nanowires and it is CMOS-compatible and is scalable to two dimensions. It was also used to study electron transport in quantum dots formed in a long carbon nanotube [See Benyamini et al. (2013)]. This process is much simpler than creating sophisticated patterns in graphene e.g., arrays of hexagonal holes, tapering the graphene with atom precision, or making graphene multilayer using CNT pillars as proposed in other theoretical works [See Yousefi et al. (2020b); Malik and Fobelets (2022); Zhao et al. (2022)].

In our study the applied voltage (force) is static although the place of the electrode can be changed in situ, however after the things are settled the same static mechanism explained here is valid. We did not study the role of dynamic (time-dependent) forces for example, applied sinusoidal voltage (force) on the heart rectification. As graphene is a nonlinear system, application of an external tone leads to mixing of different phonon modes and generation of different harmonics in VDOS. De Alba et al. (2016) et al. have shown that by applying sinusoidal force using back gate to a graphene disk, it is possible to parametrically amplify or de-amplify a vibrational mode. A recent theoretical study was brought to our attention in which the rectification ratio of 120% is achieved by application of strain gradient to a graphene nanoribbon Sai Kavuri and Sathian (2024).