Tuning graphene thermal modulator by rotating

ABSTRACT Exploring the thermal transport of graphene is significant for the application of its thermal properties. However, it is still a challenge to regulate the thermal conductivity of graphene interface. We study the interfacial thermal transport mechanism of the bilayer graphene by utilizing the molecular dynamics simulations. During the simulation, the interfacial thermal conductivity is regulated and controlled by lattice matching and tailoring. The lattice mismatched bilayer graphene model, combining the straining and torsion, can increase the interfacial thermal resistance (ITR) about 3.7 times. The variation trend of the ITR is explained by utilizing the vibrational spectra and the overlap factor. Besides, the thermal conductivity is proportional to the overlapping area. Our results show that the tailoring models can regularly control the thermal conductivity in a wide range by twisting the angle between upper and lower layers. These findings can provide a guideline for thermoelectric management and device design of thermal switch.


Introduction
Thermal transport across solid-solid interfaces has gained extensive attention because of the pervasive application in microelectronics such as energy conversion, cooling applications and so forth. The key problem is thermal resistance when heat pass through the dissimilar interface. Graphene, the two-dimensional material, hold a great promise in potential application in cooling electron devices due to its outstanding thermal [1][2][3][4] and electronic [5][6][7] properties. Over the past decade, considerable efforts have been devoted to studying and modulating interfacial thermal resistance (ITR) of multilayer graphene and graphene based interfacial systems, particularly in experiment [8,9] and theory [10][11][12][13].
Compared to the graphene-based systems, the graphene/graphene interface has a higher interfacial thermal conductivity [14][15][16][17][18]. For example, the thermal conductivity of graphene/graphene interface is 15 times that of graphene/MoS2 [14]. The ultralow ITR of graphene/graphene interface is attributed to the well-matched crystal lattice. For other graphene-based interface, the lattice pattern (also known as Moiré Pattern) has occurred at interface due to the difference of lattice constant between materials [19]. That will lead to the increase of phonon scattering and ITR at interface. In addition, when graphene sheet is twisted or shifted a certain angle or distance, the similar Moiré Pattern also occurred [20], which affects the thermal transport across graphene sheets. Cao et al. [21,22] found new electronic ground states in twisted bilayer graphene. For rotation angles of less than 1:05 � , the vertically stacked atomic regions can form narrow electron energy bands, and the electron 'correlation' effects are enhanced, which can generate a nonconducting state (a Mott insulator). When charge carriers are added to the graphene system in this state, it can be converted into a superconducting state. Nie et al. [23] studied the effect of interlayer twist angle on the ITR of six-layer graphene. Their results suggest that the ITR of six-layer graphene can increase to 1.8 times of the pristine graphene ITR (at 300 K) by modulation the twist angle between the two-layer graphene sheet. At this temperature, when the interlayer twist angle reaches 30°, the interfacial conductivity is enhanced to the minimum value about 54MW � m À 2 K À 1 . The enhanced ITR of the twisted bilayer graphene model was explained by the modification of the Brillouin zone and emergence of numerous folded phonon branches which increases the phonon Umklapp and normal scattering.
To extend its modulation range, the other ITR tuning methods (eg tailoring, strain engineering, etc.) should be introduced and combined with torsion. Chen et al. [13] reported that the ITR of few-layer graphene increases monotonically (50%-300%) by uniaxial cross-plane strain, the variation of the ITR is uncertain with compressive inplane strain and the effect of tensile in-plane on ITR can be neglect [13]. Gao et al. [24] demonstrated that the thermal conductivity of bilayer graphene kirigami heterostructure can be tuned by reasonable design of tailoring pattern and applying the external tensile strain. Therefore, the larger modulated range of ITR could be obtained by combining interlayer twist and tailoring. However, the combined effects of the interlayer twist and intralayer tailoring on ITR of bilayer graphene still remain poorly understood, which prohibits continuously growing development of aforementioned applications.
To our best knowledge, for interfacial thermal transport, the interlayer non-bonded interactions have a great influence on the accuracy of the ITR. In previous studies, when the interlayer non-bonded interactions are described by the 12-6 Lennard Jones (LJ) potential, the interface friction is underestimated. Therefore, it may lead to the underestimation of the ITR. Recently, a new Dihedral-angle-corrected Registry-dependent Interlayer Potential (DRIP) is developed, which can correctly reproduce the binding, sliding, and twisting energies and forces obtained from ab initio total-energy calculations based on density functional theory [25]. In this work, we investigate ITR of twisted bilayer graphene by using classical molecular dynamics (MD) simulations. Both DRIP and LJ potential are utilized to describe the interfacial interactions. Then, the effect of interlayer twist angle and strain on the ITR of bilayer graphene is investigated. In addition, a novel thermal switch with a high thermal tuning ratio can be realized by rational design of tailoring. The size effect is also investigated and the conclusion is extended to large scale theoretically. Our research provides a new method to regulate ITR and exhibit great potential applications in thermal management, microelectronics, and aerospace.

Modeling and methods
All MD simulations are performed by utilizing the large-scale atomic/molecular massively parallel simulator (LAMMPS) package. We designed several simulation models of ABstacked bilayer graphene to analysis that the influence of lattice mismatch and overlapping area on the ITR, as shown in Figure 1. The LJ potential and the DRIP are utilized to describe interlayer nonbonded interaction respectively, from which the effect of these two potentials on the ITR is investigated. The LJ potential has following form: Where r is the distance between two atoms. The parameters are set as σ = 0.00239 eV and ε = 0.3412 nm [26].
The DRIP functional form is: Where r −6 term models the attractive van der Waals (vdW) interactions; the exponential term is designed to capture the registry effect; f c is a cutoff function. The parameters of DRIP are obtained from Wen et al. [25].
And the Tersoff potential [27] is employed to describe the covalent interaction between carbon atoms in bilayer graphene, as shown in Figure 1(a,b). A time step of 0.5 fs is employed in the MD simulations. The nonperiodic boundaries are performed for all directions. The above parameters apply in Section 3.1.
Other simulations included strain engineering and tailoring, as shown in Figure 1(c-e). Tailored models add hydrogen atoms to passivate the tailoring edges' carbon atoms, as shown in Figure 1(e). Considered tailored models add hydrogen atoms and unifying all models' interaction potential make it easy to analyze the results, the AIREBO potential, which has been proved to be the most suitable potential for the system containing carbon and hydrogen atoms [28], is utilized to describe the intra-layer sp 2 interactions. While the LJ potential was utilized to model the van der Waals interaction among different layers. The time step is considered as 0.2 fs. And the boundaries are non-periodic in all directions.
The ITR of bilayer graphene is investigated by a transient heating technique, which mimics the experimental pump-probe approach. This approach has been successfully applied to compute the ITR of various heterostructures such as silicene/graphene and graphene/phosphorene [29,30]. To characterize the interfacial thermal resistance within the bilayer graphene, the system is initially placed into a canonical ensemble (NVT) for 10 ps at temperature 300 K. Temperature controls are applied to upper and lower monolayer graphene separately to avoid internal temperature differences. After the system reaches the steady state at a given temperature, the system converts to a micro-canonical ensemble (NVE). Then an ultrafast thermal impulse is imposed on the upper graphene layer. After this excitation, the temperature of upper monolayer graphene increases to 400 K, while the temperature of lower monolayer graphene remains at 300 K. In the following thermal relaxation process, temperature evolution in both upper and lower graphene monolayers is recorded. The total energy of each layer graphene is outputted at every time step and averaged every 100 steps to suppress data noise. Five independent simulations are respectively performed in different radii to eliminate errors caused by temperature fluctuations. The ITR can be calculated by the following equation: Where A is the interfacial area; T 1 and T 2 are the temperatures of upper and lower graphene monolayers respectively; E t is total energy of upper/lower layer graphene at time t. The specific calculation process is shown in Figure 2, and its design refers to Hong et al. [30]. The temperature evolution of the upper and lower monolayer graphene after introducing the thermal impulse is showed in Figure 2(a). After ~380 ps, the temperatures of both upper and lower monolayer graphene reach to 350 K. The profile of total energy evolution in upper monolayer graphene is also represented in Figure 2(a), which can be fit by the exponential function:E fit ¼ 0:0074 � exp À x=113:00 ð Þ þ 4:34. In addition, the relationship between the total energy (E t ) and the time integral of temperature difference ( ð ΔTdt) is obtained by integrating Eq. (1) (as shown in Figure 2(b)). It can be seen that the E t profile has a linear relationship with the ð ΔTdt. Then the E t profile can be divided into many segments, while the ITR can be determined by a linear fitting of the curve in each segment. As presented in Figure 2(c), the instant ITR vary slightly around the overall fitting results, and the average value is 0.143K � m 2 M À 1 W À 1 .
To get substantial insight into interfacial thermal transport mechanism, the vibrational density of states (VDOS) P ω ð Þ is calculated by taking the Fourier transform of the velocity autocorrelation function as: Moreover, to quantify the match between the vibration spectra of atoms at the interfaces, an overlap factor S is defined based on the correlation parameter [31], which is utilized to explore the insight of interfacial interactions: where P u ω ð Þ and P l ω ð Þdωare the phonon spectra at frequency ω of upper and lower layer graphene atoms, respectively.

Effects of the interfacial interaction potentials
Multilayer graphene exhibits vdW interactions between layers. The interface friction coefficients are significantly different while different interaction potentials are utilized. Figure 3(a) [25] shows the relationship between the rotation and the out-of-plane force by utilizing different potentials. The DRIP results are validated correctly. However, the results of using the LJ potential shows that the force in out-of-plane direction is independent of the rotation. Thus we utilize both potentials to study and analyze on the ITR of bilayer graphene with twist angles.
The radius of the bilayer graphene is 60 Å with the fixed boundary of 5 Å in width. And then twist the upper monolayer graphene. Several twist angles ranging from 0° to 30° are selected to compare the ITR difference between these two potentials. The results in the range from 30° to 60° can be obtained according to the symmetry of the model structure. From the results shown in Figure 3(b), a similar variation trend can be observed between the ITR results of these two potentials, at the twist angle of 0°(60°) has the smallest ITR. Meanwhile the thermal conductivity is the largest. When the twist angle is 30°, the ITR reaches a maximum value. It can be noted that these two interaction potentials don't have a significant impact on the ITR, although their interface friction is obviously different. For example, as for bilayer graphene with the twist angle of 30°, the ITR value of using the DRIP is 0.160 K � m 2 M À 1 W À 1 , compared to the ITR result of utilizing the LJ potential is 0.148 K � m 2 M À 1 W À 1 . Therefore, we utilize LJ potential for the follow-up work to reduce the amount of calculation.

Effects of the lattice mismatch
Twist can lead to lattice mismatch, which sequentially has an effect on ITR. Straining the graphene also has a similar way of the lattice mismatch for affecting ITR. Therefore, the ITR of bilayer graphene can be regulated by combined utilization of interlayer twist and applying global tensile strain. According to the relationship between interlayer twist angle and the ITR of bilayer graphene, the maximum ITR is 0:062K � m 2 M À 1 W À 1 at the interlayer twist angle of 30°. Therefore, the global tensile strain applied to bilayer graphene with interlayer twist angle of 30° can obtain a larger ITR. Thus, in this part, the models' twist angle is 30°, and then straining the upper layer (as shown in Figure 4(a)) or both layers of graphene.
The global tensile strain of the model can be altered by changing the lattice constant of the graphene, which is defined as: ε ¼ l À l 0 ð Þ=l 0 , where ε denotes the global tensile strain, the l and l 0 are the lattice constants of the graphene for the unstrained and strained structure, respectively. The relationship between the ratio of ITR and the global tensile  [25]. The rotation include twisting and sliding. When the rotation changes from 0 � to � 60 � , the bilayer configuration accordingly changes from AB staking to AA staking. (b) The ITR results of bilayer graphene which respectively utilize DRIP potential and LJ potential to describe interfacial interaction.
strain can be seen in Figure 4  graphene can increase the lattice mismatch, which enhances the phonon scattering, thus the ITR increase.
In addition, to obtain an in-depth understanding of the strain dependent ITR, we investigated the VDOS at interfaces by analyzing the correlations of atomic vibrations. The in-plane and out-of-plane vibrational spectra of carbon atoms in the upper and lower graphene layers under tensile strains from 0% to 10% are represented in Figure 5(a,b). The in-plane and out-of-plane vibrational spectra of the upper and lower graphene layers completely overlap when the system is free of local strain. The higher frequency peak denotes the G-band, and it is located at 52.1 THz. When global tensile strain is applied on the upper layer graphene, the frequencies of in-plane phonon modes (especially those in the G-bands) undergo a red-shift under tension (see Figure 5(a), from 52.1 THz at strainfree to 28.7 THz at strain of 10%, which is labeled by a red arrow). This trend agrees well with previous studies [16].
Moreover, to quantitatively determine the match of VDOS between the upper and lower graphene layers, the overlap factor S is calculated (as shown in Figure 5(c)). In order to compare the variation trend of S more conveniently, all results in the figure are normalized with the strain-free result S 0 as the standard. For the overlap factor S 0 in the absence of tensile strain. Due to the red-shift of the in-plane G-band in upper layer graphene, the overlapped area of VDOS in upper and lower layers is reduced (Figure 5(a)). As the tensile strain increases to 10%, the overlap factor S decreases around 64.4%. Under the same situation, the out-of-plane overlap factor S only declines about 11.4%. Therefore, it can be summed up that the in-plane phonon coupling plays a main role on the ITR.

Design of thermal switch
In the previous simulations, the results show that the lattice mismatch can alter the ITR at a small variation. Whether the lattice mismatch is caused by torsion or strain, the thermal conductivity of the graphene interface is not reduced to near 0, and its regulation is difficult to achieve linear changes. Nevertheless, our previous research results show that the thermal conductivity of the tailored graphene can theoretically be reduced to 0, becoming a thermal insulation material.
Firstly, we present a conceptual design of bilayer by assembling tailored graphene, which can seem as a thermal switch, to control ITR. The radius of the tailored models is 100 Å with 5 Å fixed boundary, and the inradius which is fixed region is 10 Å. Then tailor two sectors of the bilayer graphene. Each tailored angle of the sector is 105°.
When the graphene interface lattice mismatch is not considered, the ITR can be controlled by twisting the tailored graphene. Its mechanism of control is regulating the overlapping area of the upper and lower graphene layers, as shown in Figure 6. The overlapping area reaches the minimum value, when continually twisting the tailored graphene. At this state only the fixed regions are overlapped, and this state can be considered as the off-state of the thermal switch, as shown in Figure 6(a). The other state of the thermal switch is the on-state when the maximum value of the overlapping area is reached, as shown in Figure 6 (d). The process from off-state to on-state is shown in Figure 6(a-d). The variation of the ITR in this process may have a linear regulation which can be accurately controlled, because the overlapping area can be linear regulated, as shown in Figure 6(e). In order to make the thermal conductivity formula better match the simulation results, the thermal conductivity formula can be split into two parts, namely the graphene part and the boundary area part. During the process of twisting the tailored graphene model, the ITR can be described as this equation where G free and G edge are the interfacial thermal conductivity of the free region and fixed region, A 0 is the area of the graphene without tailoring,A free and A edge are the area of free region and fixed region.
Considered the design of the tailored graphene model may have an impact on the results, another model is designed to analysis, as shown in Figure 7(a). The fixed boundary is same as the previous tailored model, the number of the tailored sectors is increased to 3. For keeping the total tailored area the same as the previous tailored graphene model, each tailored sector's angle is 70°. Meanwhile, the size effect is also considered. These two kinds of tailored models with the radii of 60 Å, 80 Å, 100 Å, and 120 Å are utilized to analyze this effect. The ITR results of the tailored graphene model shown in Figure 7(b), the results demonstrate that the difference of the tailored graphene models doesn't have a significant effect on the ITR when the models are in same size and have the same overlapping area.
Then enlarging the range of model size, only the tailored graphene model with two tailored sectors is utilized for simulation to further elaborate the size effect. The model radii of 70 Å,140 Å, and 160 Å are added, the result is shown in Figure 7(b). As the size increases, the linear relationship between the interfacial thermal conductivity and the area of graphene becomes more obvious. The relationship between them is: where A is the overlapping area of the tailored graphene model. As the radius increases, the off-state interfacial thermal conductivity decreases and the onstate ITR increases within a narrow range. The ITR ratio of on-off state is shown in Figure 7 (c). The size of model plays a dominant role in the maximum variable range of ITR. When the system is large enough, the fixed edge area can be ignored. The thermal switch can offer maximum thermal insulation in the off state.

Conclusion
We have performed the simulations to study the effect of the lattice mismatch and the overlapping area on interfacial thermal resistance of bilayer graphene. In this work, we first utilize the DRIP and LJ potentials, respectively, to describe the interfacial interaction to calculate the ITR. There is no obvious difference in the ITR calculated by them. It shows that the interface friction has negligible effect on the ITR. Then, we analyze the effect of the lattice mismatch caused by the torsion combined strain on ITR. The VDOS and overlap factor show that the in-plane phonon coupling has a dominant effect on the ITR. Finally, we design a thermal switch model. Accurately control the interfacial thermal conductivity by controlling the overlapping area. And it is found that the thermal switch model has obvious size effect. When the size of this model increases to a certain extent, it can be considered that the off-state can be insulated. Our investigation about the interfacial thermal transport across the bilayer graphene provides a fundamental knowledge to the design of the thermal switch which can linearly regulate the interfacial thermal conductivity.

Disclosure statement
The authors declare no competing financial interest.