Thermal Transport between Graphene Sheets and SiC Substrate by Molecular-Dynamical Calculation

Using nonequilibrium molecular dynamics, we investigate the mechanisms of thermal transport across SiC/graphene sheets. In simulations, 3C-, 4H-, and 6H-SiC are considered separately. Interfacial thermal resistances between Bernal stacking graphene sheets and SiC (Sior C-terminated) are calculated at the ranges of 100K∼700K. The results indicate, whether Si-terminated or C-terminated interface, the interfacial thermal resistances of 4Hand 6H-SiC have similar trends over temperatures. Si-terminated interfacial thermal resistances of 3C-SiC are higher than those of 4Hand 6H-SiC in a wide temperature range from 100K to 600K. But, for C-rich interface, this range is reduced from 350K to 500K.


Introduction
For the atomically thin structure of graphene, it has attracted much attention due to potential applications to thermoelectric and photoelectric devices.Over the past decade, researchers mainly focus on the fabrication of graphene and its physical characteristics.The most common used methods of graphene fabrication are epitaxial growth on silicon carbide (SiC) and chemical vapor deposition (CVD).As to epitaxial growth method, Si atoms sublimate from a single-crystal SiC substrate and create large area graphene sheets.CVD method employs carbon source gas to react with the specific substrate chemically, so as to synthesize graphene sheets.In recent years, the techniques of graphene production have been improved dramatically.Tzalenchuk et al. [1] synthesized 50 um 2 graphene sheets by epitaxial growth method.Soon after, Bae et al. [2] made graphene films with 762 mm diagonal length onto Cu substrate.In 2011, Zhang et al. [3] first realized the low temperature growth of graphene onto several substrates such as SiC, SiO 2 , and Al 2 O 3 by plasma enhancement chemical vapor deposition (PECVD).Through solid-state molecular beam epitaxy (SSMBE), Tang et al. [4] made graphene films grown on Si(111) substrate and analyzed the effects of SiC buffer layer on stabilizing the surface of Si substrate.Janssen et al. [5][6][7] investigated the Hall resistance in epitaxial graphene grown on Si-terminated SiC and other graphene material systems.At present, although numerous studies have demonstrated the high thermal conductivities of graphene and graphene sheets along tangent direction, little work has been done to address the characteristics of graphene sheets along normal direction, let alone SiC/graphene composite films.

Theory
As to theoretical study, researchers usually adopt analytical or numerical approaches to perform nanoscale heat transport simulations, involving molecular dynamics (MD) method, Monte Carlo (MC) method, Boltzmann transport equation, and Green's function.Compared with other methods, MD is much closer to physical reality, which defines the simulation region with atoms, supposed to the arrangement of chemical structures.Based on Newton's law, MD sets cut-off radius to build atom groups, in which atoms interact with each other by elastic forces.At the beginning of MD simulation, motion parameters of each atom should be initialized according to the imposed temperatures.After enough circulations, system reaches to an equilibrium and the temperature can be obtained by Boltzman equipartition theorem (BET).According to the different calculation principles, MD can be divided into equilibrium MD (EMD) and nonequilibrium MD (NEMD) methods.As to EMD, thermal conductivity may be obtained by Green-Kubo linear response under equilibrium state [8].For NEMD, through modulating the temperature of hot bath and cold bath, thermal flux is created and temperature gradient emerges in the system, and then thermal conductivity can be evaluated by Fourier's theorem [9,10].NEMD method is analogous to the experimental measurements of heat conduction and is preferred over the EMD-based approach to compute thermal conductivity of inhomogeneous material systems [11,12].Recently, low-sized nanoscale materials, such as core-shell nanowires, superlattice, graphene, and nanotube [13][14][15][16][17][18][19][20][21][22], were investigated by NEMD method.In this paper, we use a NEMD-based approach to study the thermal performance of SiC/graphene composite films.
The key points of NEMD method are the creation of thermal flux and the interaction potential.In order to create thermal flux, we may keep the temperature gap between the hot part and the cold part, which drives thermal flux.Another way is to swap the coldest atoms in the hot part with the hottest ones in cold part.From the computation viewpoint, the latter is more effective, because the number of atoms involved in the calculation is less than the former.Therefore, we utilize the former way to drive thermal flux in NEMD simulations.
In semiconductors, heat transfer mainly depends on the transmission of elastic vibration by interatomic force, which is the derivative of the interaction potential.Hence, whether the right interaction potential is utilized determines the correctness of the results.In this work, Tersoff [23,24], Airebo [2], and Lennard-Jones (LJ) functions are used together to construct interaction potential energy model.Tersoff potential as three-body interaction function is adopted to calculate the interatomic force of Si-C covalent bond in SiC material.Airebo and Lennard-Jones functions belong to twobody potential interaction.For a single sheet of graphene, the C-C short-range force can be solved by the former.The latter is used to evaluate C-C long-range force between the graphene sheets.
A NEMD ball and stick 3C-SiC film model in our study is depicted in Figure 1.Along x, y, and z direction, this model has three periodic boundary conditions.Hot zone is set in the middle and cold zone is set to the left part.Heat flux flows from the hot zone to the left side and right side.After equilibrium, temperature gradient may be created in this model.
In Figure 2, the results of 3C-SiC, graphene, and Si films by NEMD method are compared to confirm its validity.The cross section of this system is 6 × 8 unit cells, with a length of 140 unit cells at 100 K.We set the time step as 0.766 fs and the average temperature of this system equals 100 K.The processes of this simulation can be divided into two stages.The first stage is to relax atoms, and the next one is to establish stable temperature gradient.As shown in Figure 2, for 3C-SiC, graphene, and Si films, the fitting curve of temperature distribution shows a good linear dependence of cell serial number.Because of the interface scattering, there exist temperature jumps prominently near the hot bath and cold bath [11].The right inset figure in Figure 2 shows the heat flux, which equals the ratio of the total transferred kinetic energy to time and the cross-sectional area.
Exchanging kinetic energy of atoms in hot bath and cold bath, system reaches to an equilibrium state and thermal conductivity can be obtained from the following equation: where  denotes heat flux, ∇ is the gradient of temperature profile,     represents the area of cross section normal to thermal flux, () is the mass of atom in the cold bath (hot bath),  is the simulation time, and V ℎ (V  ) indicates the velocity of the coldest (hotest) atoms in hot bath (cold bath).We also calculate the in-plane thermal conductivities of monolayer graphene on SiO 2 substrate, as shown in Figure 3.
In this NEMD model, Tersoff, Airebo, and LJ potential functions referring to Zhao's work [25] are employed to calculate the chemical bonding interactions of C-C, Si-Si, Si-O, and O-O.In simulations, NEMD results are about 50% lower than the experimental values [26], which is perhaps due to quantum errors between realistic temperatures and molecular dynamics temperatures.Thus, the proximity between the experimental values and the NEMD results confirms the validity of the physical model we built.Interfacial thermal resistance  can be calculated from where Δ is the temperature drop at the interface of graphene/SiC.Figure 4 shows the schematic diagram of a graphene/SiC composite film and the lattice unit cell structures of 3C-, 4H-, and 6H-SiC polytypes.In simulations, we utilized AB (Bernal) stacking graphene sheets in place of AA, ABC stacking, because above 80% monocrystalline graphite belongs to AB stacking in reality [17,27].
The lattice parameters used in our work are listed in Table 1.a, b, and c are the lengths of the three axial-vectors.In Figure 5, we build a model of 4H-SiC/graphene sheets and calculate its temperature distribution, at the average temperature of 300 K. Hot bath locates at the middle of graphene sheets and cold bath is set at the right end of 4H-SiC film.The enlarged graph in Figure 5 shows a temperature drop at the interface due to the interfacial thermal resistance while on both sides of the interface, the fitting curve of temperature distribution is still linear.It is noteworthy that the profiles of temperature are not symmetrical and the values of Δ on the two sides of the hot bath are different, which is because the interfacial morphologies are not identical strictly.
The lattice mismatch of SiC/graphene system can strongly affect the special thermal properties and electronic characteristics.Correspondingly, it is necessary to investigate the influences of different interfacial morphologies on thermal transports.As shown in Figure 6, we investigate the averaged interfacial thermal resistances between Si-terminated SiC polytypes and graphene sheets.Referring to the work of Hass and so forth, we set the distance between the first graphene layer and Si layer of SiC as 2 ∼ 2.5 Å.From Figure 6, 4H-and 6H-SiC have similar convex increasing trends of interfacial thermal resistances.The interfacial thermal resistance of 3C-SiC initially increases with the increasing temperatures from 100 K to 300 K and then shows a downward bending above 300 K.The results indicate 3C-SiC has larger interfacial thermal resistances than those of 4H-and 6H-SiC below 600 K and has an obvious peak value about 2.5×10 −9 Km 2 /W at 300 K.Over 400 K, the values of 6H-SiC similarly conserve to a constant about 1.75×10 −9 Km 2 /W.Different from 3C-SiC and 6H-SiC, 4H-SiC keeps linear increasing trend from 100 K to 700 K.
In the CVD graphene growing method, we can also use C-rich surface to bond with C atoms.The interface of SiC/graphene becomes different from Si-rich interface as depicted in Figure 7. Compared with Figure 6, the maximum values of the three SiC polytypes are larger than those in Figure 7.The curve of 3C-SiC complies with approximately Gauss distribution and the peak value locates at about 400 K.The trend of 6H-SiC is a little changed; however the interfacial thermal resistances of 4H-SiC may not keep linear and show convex profile as 6H-SiC.We suggest the bond energy of C-C is stronger than that of Si-C bond; that is, C-C bond has larger potential energy.Therefore, C-rich surface has bigger interfacial thermal resistances than Si-rich surface.

Conclusion
By chemical method, monolayer graphene and graphene sheets can be synthesized on SiC substrates.Many researchers have studied electrical and thermal properties of graphene; however, few works involve the thermal transport of SiC/graphene sheets.In this work, we apply a dedicated NEMD model to investigate interfacial thermal resistances of SiC/graphene sheets.Three polytypes of SiC (3C-, 4H-, and 6H-SiC) are set as substrates to bond with graphene sheets in simulation.Si-rich and C-rich interfaces between SiC/graphene are also analyzed at the temperatures from 100 K to 700 K. From the results of Si-rich interface, the interfacial thermal resistance of 3C-SiC is larger than those of 4H-and 6H-SiC below 600 K. From 100 K to 700 K, the curve of 6H-SiC is always greater than that of 4H-SiC.By fitting the calculated datum, the peak interfacial thermal resistances of 3C-, 4H-, and 6H-SiC are 2.3 × 10 −9 Km 2 /W (at 354 K), 1.74 × 10 −9 Km 2 /W (at 571 K), and 1.16 × 10 −9 Km 2 /W (at 457 K).As to C-rich interface in Figure 7, the trends of 4Hand 6H-SiC are not changed greatly, but the curve of 3C-SiC approximates to Gauss distribution and its maximum reaches to 2.95 × 10 −9 Km 2 /W (at 480 K).Especially, the values of 3C-SiC become lower than those of 3C-SiC and 4H-SiC over 600 K.
However, it should be acknowledged that surface reconstructions of SiC such as (3 × 3), ( √ 3 × √ 3)30 ∘ , and (6 √ 3 × 6 √ 3)30 ∘ are not considered in this work, which needs to be further investigated in the near future.We expect to give a theory reference to the applications of SiC/graphene films in the field of micro/nanothermoelectric devices.

Figure 1 :
Figure 1: Ball and stick model of 3C-SiC film in this work.

Figure 2 :
Figure 2: Temperature distributions of 3C-SiC, graphene sheets, and Si across cell serial number.(a) The enlarged figure.(b) The thermal flux.