Thermal rectification in ultra-narrow hydrogen functionalized graphene: a non-equilibrium molecular dynamics study

In this study, the non-equilibrium molecular dynamics simulation (NEMD) has been used to evaluate the thermal properties, especially the rectification of ultra-narrow edge-functionalized graphene with hydrogen atoms. The system’s small width equals 4.91 Å (equivalent to two hexagonal rings). The dependence of the thermal rectification on the mean temperature, hydrogen concentration, and temperature difference between the two baths was investigated. Results reveal that the thermal rectification increases to 100% at 550 K by increasing the mean temperature. Also, it is disclosed that hydrogen concentration plays a vibrant role in thermal rectification. As a result of maximum phonon scattering at the interface, a thorough rectification is obtained in a half-fully hydrogenated system. As well, the effects of temperature difference of baths ΔT on thermal rectification has been calculated. As a result, the thermal rectification decreases even though the current heat increases with ΔT. Finally, the thermal resistance at the interface using a mismatching factor between the two-phonon density of states (DOS) on both sides of the interface has been explained.


Introduction
Graphene [1], a two-dimensional (2D) monolayer sheet with honeycomb form, has drawn more attention in the past decade due to its outstanding properties. Graphene shows excellent optical and electronic [2,3], mechanical [4], and thermal properties [5]. Besides graphene, novel 2D materials have been discovered or synthesized functionalized graphene for different purposes, such as carbon nitride [6], molybdenum disulfide [7], borophene [8], phosphorene [9], and graphene [10]. Functional atoms or groups can affect physical properties. For example, graphane was theoretically predicted [10] and synthesized [11,12] by hydrogenating both sides of the graphene in an alternating manner. In this situation, the sp 2 bond type is transformed to sp 3 and affects all properties. In addition, the number of hydrogen atoms on graphane can be changed to control properties. Moreover, we investigated the thermal conductivity of other systems and explored the parameters that affect the thermal conductivity [13][14][15].
The hybrid of 2D monolayers such as graphene and graphane (functionalized graphene with an arbitrary number of hydrogen atoms) is also enjoyable. Thermal rectification in such systems was investigated [16,17] recently. In thermal rectifier (thermal diode) systems, the heat current across the system depends significantly on the thermal gradient direction. As electrical diodes playing a vital role in electronics, thermal diodes have substantial implications for many thermal management applications, from nanoscale calorimeters to microelectronic processors to macroscopic refrigerators and energy-saving buildings [18]. Therefore, the heat current prefers one direction to another one. Rajabpour et al. [16] have studied the graphene and graphane junction to estimate thermal rectification. They obtained thermal rectification of about 20% through NEMD simulation. Moreover, they also examined thermal properties at the interface, which is responsible for phonon scattering, leading to thermal resistance (known as Kapitza resistance). They also showed that the Kapitza resistance depends on the temperature gradient, which causes the heat current to prefer one direction to the opposite.
Shavikloo et al. [14] studied thermal rectification in graphene and graphane with grain boundary interfaces between two sides. They examined thermal rectification by randomly removing hydrogen atoms from graphane. They found that the thermal rectification is less than 10% at room temperature. Furthermore, they investigated the effect of the hydrogen concentration on thermal rectification in the hybrid system. Removing hydrogen helps us tune the thermal properties so that they have succeeded in increasing thermal rectification by 17% by removing 40% more hydrogen atoms than fully hydrogenated graphane. In the two above studies, the system's width was about 100 Å, and periodic boundary conditions were applied to in-plane directions. The width and boundary condition can be influenced the thermal properties. Functionalization with hydrogen has been used in other systems such as single-walled carbon nanotube [19], nanoporous graphene [20], and pristine graphene [21] with particular patterns. Moreover, thermal rectification has also been seen in other systems, in graphene-C 3 N junction [22,23], graphene-C 3 B [24], graphene-pillared [25], radial rectification [26], telescopic silicon nanowire [27], and so on. Thermal rectification in ultra-narrow graphene with functionalized edge via NEMD simulation has been investigated. The dependency of thermal rectification on temperature and the temperature difference between baths and hydrogen concentration was studied. Additionally, the phonon density of states was obtained to explain thermal resistance at the interface.

Computational method
In semiconductors, the contribution of electrons in the thermal transport is less than the phonon's contribution. Therefore, we ignore the electrons and compute phonon thermal current in our ultra-narrow system. This study implemented all NEMD simulations using a Large-scale Atomic/ Molecular Massively Parallel Simulator (LAMMPS) [28] as a molecular dynamics program. The desired structure is indicated in Fig. 1 with a length of 200 Å and a width of 4.91 Å. Both edges in the left half of the structure are fully functionalized with hydrogen atoms, while the right side has no hydrogen.
The hydrogen concentration can be different in some simulations due to investigating the effect of that on thermal rectification. Periodic boundary conditions along the x-direction, but free boundary conditions were assumed in other directions. Moreover, to model the interactions in this system between C-C and C-H, we employed the Adaptive Intermolecular Reactive Empirical Bond Order (AIREBO) [29] potential with a cutoff 3.0 (σ scale factor). This potential is a three-body type and is appropriate for modeling breaking, creating bonds, and describing hydrocarbon systems. The Velocity Verlet algorithm integrates the equation of motion of atoms with a time step of 0.2 fs. The time step was considered small because this system has hydrogen atoms that can oscillate with high frequency.
First, to relax the system and eliminate all extra stress, the NPT ensemble with Nose-Hoover [30] barostat and thermostat was applied to the system along x-direction so that the pressure tensor element P xx becomes zero at room temperature. In the mentioned situation, the system can expand or shrink if desired. After relaxing, as represented in Fig. 1, two slabs at two ends of the system were considered fixed atoms so that the particles inside of them could not move in during simulation time. Moreover, two other slabs near two fixed layers were selected as thermal baths (hot and cold sources). We used an NVT ensemble with a Nose-Hoover thermostat to keep the temperature fixed in thermal sources. The NVE ensemble was applied elsewhere.
Here, to estimate the thermal rectification, we use the equation below [20], where Q L→R ( Q R→L ) is the heat current from left to right (right to the left). We assume that in Eq. 1, Moreover, the phonon density of states (DOS) on both sides of the interface was calculated using the Fourier transform of the velocity autocorrelation function. Therefore, we choose two slabs on both sides of the interface and use the equation below [31,32], where s is the atom type and symbol < > indicates the velocity autocorrelation function. Parameters of m, k B, and T are the mass, Boltzmann constant, and system temperature. The is the angular frequency. (1)

Results and discussion
Obtaining heat current across a system (in the x-direction), the accumulative energies added to and removed from the heat baths were calculated (see Fig. 2). Indeed, the slope of the linear fit to energy profile dE∕dt is the heat current. As indicated in Fig. 2, the input and output energies from the baths and their slope are equal. It shows that the heat current reached to steady state. The heat current depends significantly on the temperature gradient sign because there is a mass asymmetry across the system. It causes the phonons to prefer to move more in one direction than the opposite one, usually from thick to thin.
Moreover, we present the sensibility of the heat current and the thermal rectification of the mean temperature of the system. Therefore, we enhance the temperature from 200 to 650 K.
As shown in Fig. 3a, the heat current decreases with increasing the temperature. The behavior of the heat current in both temperature gradient biases arises from two phenomena. First, the phonon-phonon scattering will increase by increasing temperature, which is a negative factor for heat current. On the other hand, by increasing temperature, the high frequency (high energy) phonons will excite, leading to increased heat current (positive factor). Competition between these two factors specifies the behavior of the heat current. Here, the simulation results indicate that the phonon-phonon scattering has been dominant. Furthermore, the heat current along the positive direction is greater than the negative one. It is the origin of thermal rectification.
In Fig. 3b, the thermal rectification increases with increasing temperature. It should be noted that the decreasing rate of negative heat current than positive leads to increase thermal rectification. The magnitude of thermal rectification is more significant than 30% in this case. Figure 3b shows the estimated trend via two fitted lines in Fig. 3a. Two linear fits calculated the estimated trend in Fig. 3a. Consider L1 and L2 in Fig. 3a; for a given temperature, we calculated (L1-L2)/(L2)*100 for the estimated trend. The estimated trend remains valid as long as the heat flux behavior is linear. In analogy with other studies by Rajabpour [16], Shavikloo [17], Yousefi [20], Cartoixà [27], we get a higher thermal rectification. Therefore, we introduce a promising thermal rectifier.
Thermal rectification under various conditions, for example, different ΔT. This parameter can be changed under application circumstances. Therefore, we check its effect on heat current and thermal rectification.
According to Fourier's law ( J = − ∇T ), the heat flux is proportional to ΔT. Therefore, we expect that J increases with increasing ΔT. The NEMD result confirms this proportionality, and the heat current increases with increasing ΔT, as shown in Fig. 4a. For the positive direction, the slope of the fitted line is greater than the negative due to thermal conductivity depending on the temperature gradient sign. Although the difference between two heat currents gets larger, the heat current for the negative direction gets more significant too. Therefore, we obtain a decreasing behavior of thermal rectification as indicated in Fig. 4b. Thermal  One of the crucial factors that help us tune the thermal rectification is hydrogen concentration in half of the system. Figure 5a and b explored the effect of removed hydrogen atoms in the system in 5-95%. For any amount of removed hydrogen, the heat current in the positive direction is more significant than the negative. For 95%, due to the dilution of hydrogen atoms on the left side, the heat currents increase in both directions and close to each other because the interface is vanishing. On the other hand, the difference between the heat currents arises from thermal rectification. Figure 5b shows that the NEMD result indicates that the thermal rectification suppresses in the mentioned range by removing hydrogen. With a slight decrease in hydrogen (5%), the thermal rectification drops from 40 to 18%. With more removing the hydrogen atoms up to 50%, thermal rectification increases slightly. This behavior was reported already elsewhere [17]. This shows that the hydrogens are responsible for phonon scattering and thermal rectification. Thermal rectification beyond 65% decreases because the phonon scattering decreases at the interface in a negative direction, leading to an increase in heat current.
To understand the underlying mechanism of the thermal rectification, we should investigate the origin of the thermal resistance at the interface (known as Kapitza) [33,34]. Kapitza originates from phonon scattering at the interface. Therefore, a formalism called interface conductance modal analysis (ICMA) developed by Gordiz and Henry [35,36] applies to knowing the number of phonons is conductance or thermal resistance at the interface. The phonon density of states (DOS) on both sides of the interface in positive and negative temperature gradient signs specify mismatching (see Fig. 6). The mismatching indicates that some phonon modes are present on only one side. This is attributed to the Kapitza resistance at the interface, which depends on thermal gradient bias. Moreover, the mismatching factor between two DOSs on both sides of the interface can be calculated through the equation below [20], where D 1 and D 2 are the DOS on the left and the right side, respectively. For positive and negative temperature gradients, the mismatching factor is 0.01253 and 0.01332, respectively, for T = 300K and ΔT = 40K.

Conclusion
Using a series of NEMD simulations, we reported the thermal rectification in ultra-narrow hydrogen functionalized graphene. Thermal rectification is significantly high and even increases up to 100% in temperature 550 K. The dependency of thermal rectification on mean temperature, ΔT, and hydrogen atoms concentration was investigated. We found that the thermal rectification in the narrow structure in this work increases with increasing mean temperature while then decreases with increasing ΔT. Furthermore, we can tune the thermal rectification between 10-30% by removing some percent of hydrogen atoms. Finally, by calculating the phonon density of states on both sides of the interface, the underlying mechanism of the phonon scattering at the interface was shown.
Author contribution All authors have the same contribution.

Declarations
Ethics approval On behalf of all authors, the corresponding author approves.

Consent to participate
The corresponding author consents on behalf of all authors.

Consent for publication
On behalf of all authors, the corresponding author consents to publish under the subscription model. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Competing interests N/A.
Springer Nature or its licensor holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.