Abstract
The thermal conductivity (TC) of isolated graphene with different concentrations of isotope (C13) is studied with equilibrium molecular dynamics method at 300 K. In the limit of pure C12 or C13 graphene, TC of graphene in zigzag and armchair directions are ~630 W/mK and ~1000W/mK, respectively. We find that the TC of graphene can be maximally reduced by ~80%, in both armchair and zigzag directions, when a random distribution of C12 and C13 is assumed at different doping concentrations. Therefore, our simulation results suggest an effective way to tune the TC of graphene without changing its atomic and electronic structure, thus yielding a promising application for nanoelectronics and thermoelectricity of graphene-based nano device.
1. Introduction
Since it was fabricated in 2004 [1], graphene, a monolayer of sp2-bonded network of carbon atoms, has attracted much attention for its unique electronic properties [2]. Meanwhile, both recent theory and experiment studies [3, 4] have revealed that isolated graphene has shown an unusual high thermal transport capability, which is of great importance in thermal management of nanoelectronics. Different from conventional metallic materials, thermal energy carriers for graphene are mainly in the form of phonon vibrations [3–5], and phonon contribution to TC is approximately 50 times larger than electron contribution at room temperature [3], which suggests that electron thermal transport is negligible in our case. Furthermore, the electronic contribution to TC would be independent of isotope effect as is electronically identical to . For these reasons, we will study only the phonon contribution to TC in this paper.
In addition to utilizing its high TC, another possible application of graphene has been investigated for thermoelectric energy conversion [6], where low TC but high electric conductivity for graphene is required for obtaining high thermoelectric efficiency. Such efficiency is expressed as , where σ is the electric conductivity, S is the Seebeck coefficient, and is the thermal conductivity. To achieve high ZT for graphene, a general scheme is to minimize while keeping and S less changed. One practical method is to dope the graphene with the stable isotope since electronic structure of graphene is unchanged in this doping. In a pure crystal, one without defects or dislocations, phonon scattering in the presence of different isotopes has been strongly correlated with changes in thermal conductivity. Similar to the observation for isotope effects on TC of Ge, diamond, and boron nitride nanotubes [7, 8], it is interesting to check how effectively isotope-doping method will reduce TC of Graphene.
Modeling of thermal transport can be achieved by using Boltzmann transport equation (BTE) [9, 10], or molecular dynamics (MD) simulations [11]. In BTE method, the single mode relaxation time (SMRT) approximation is a commonly used technique involving the assignment of the relaxation time to different phonon scattering mechanisms. The relaxation time can be either fitted to the experimental TC value [12] or determined from MD simulations [10]. In MD simulation approach, TC can be predicted from either nonequilibrium MD [13], where a temperature is applied across the simulation cell, or equilibrium MD [14], where the so-called Green-Kubo method is used to compute TC from heat current fluctuations. One advantage for MD method is that there is no assumption needed for phonon interactions. As long as the phonon dispersion and anharmonicity of the potential are accurate, MD method provides a robust way to accurately compute thermal transport.
In Section 2, we introduce Green-Kubo method and describe the simulation procedures to compute TC. In Section 3, we show the isotope effects on thermal transport of graphene, and discuss the reason for TC reduction and its possible application for thermoelectric application. Section 4 presents a summary and conclusions.
2. Green-Kubo MD Simulation Method
The Green-Kubo formula [14] derived from linear response theory can express thermal conductivity tensor in terms of equilibrium heat current-current autocorrelation in the form where is the system volume defined as the area of graphene multiplied with van der Waals thickness (3.35 Angstrom), kB is the Boltzmann constant, T is the system temperature, and τm is the time required to be longer than the time for current-current correlations to decay to zero [11]. is denoted as the heat flux in α or β direction, its expression is commonly defined as where , , and are the position, velocity and site energy of atom , respectively. is the potential energy at site . Then, we use Hardy’s definition [15], where is the force exerted on atom i by atom j. One merit of Hardy’s definition is that it is independent of pair-wise or many-body potential formulas. Based on the above equations, heat current at each step is recorded on disk as a quantity defined by atom positions, velocities, and interatomic forces, which can be extracted from atomistic simulations. The last step is to compute TC tensor with the known heat current in (1).
In the equilibrium MD simulations, we used the second generation REBO carbon potential for its accuracy in describing bond strength and anharmonicity of carbon materials [16]. A unit cell of graphene, with a periodic boundary condition as shown in Figure 1, has been thermalized to 300 K for 200 ps with Berendson thermostats. Afterwards, the heat current of graphene is recorded every 2 fs in nine microcanonical ensemble simulations with uncorrelated initial conditions. Each microcanonical simulation needs to run up to 10 ns to obtain converged TC value. Compared with nonequilibrium molecular dynamics (NEMD), the Green-Kubo method is indeed less computationally efficient; however, it is free of troubles such as finite size and boundary effects, which are commonly unavoidable in NEMD.
3. Results and Discussion
We computed the TC of mass defect-free graphene as 630 W/mk and 1000 W/mk in armchair and zigzag direction, respectively. This is lower than the reported experimental data [4, 5]. The reason has to do with the discrepancies on phonon dispersion between experiment measurement and theoretical data predicted with original REBO potential [18]. Although the absolute TC of graphene is overall underestimated in our simulation, the relative difference of TC caused by isotope mass defect is still physically meaningful, and the normalized TC reduction seen in Figure 4 properly reflects the isotope effect on TC of graphene.
One primary test done before studying isotope effects on thermal transport is the convergence test for graphene with different unit cell periodic boundary lengths. As shown in Figure 2, we have tested three cases with the boundary length of 1.6, 3.4, and 5.1 nm respectively, and the converged result suggests that the majority contribution from different phonon vibration modes are well included in our simulations. Then, to be efficient, we choose the 1.6 nm unit cell as the structure for the following simulations.
The first thing we learned from the simulations is that TC of pure or graphene in zigzag direction is 58% greater than that in armchair direction. This can be conceptually explained from the acoustic phonon dispersion and Grüneisen parameter of graphene calculated with REBO potential as shown in Figure 3. In the SMRT approximation, the BTE method can express TC as a summation of each phonon mode contribution, and where and represent phonon modes and wave momentum, and is the specific heat of the phonons, which is constant in the classical case. is the group velocity for mode , is the relaxation time, and T is the temperature. From Figure 3(a), we notice that both out-of-plane acoustic (ZA) and transverse acoustic (TA) modes in zigzag direction have apparently larger group velocity than that in armchair direction, longitudinal acoustic (LA) modes group velocity in both directions is quite similar. Then, from Figure 3(b), the Grüneisen parameter, which is plotted as a function of phonon mode and momentum , shows a smaller difference in zigzag and armchair directions for each vibration mode. As the Grüneisen parameters for LA, TA, and ZA are similar, it is reasonable to assume that the relaxation times for phonon vibration in zigzag and armchair direction are the same. Based on the above analysis, we conclude that the difference in the TC between zigzag and armchair directions mainly comes from their different group velocities.
(a)
(b)
(a)
(b)
In the isotope effect study, we generated a wide range of graphene samples with different C13 concentrations randomly distributed, since it is a more realistic possible configuration after the synthesis of graphene. Recently, Mingo et al. have proposed a possible method to generate isotope clusters in graphene, and theoretically demonstrated TC reduction using nonequilibrium Green’s function method [19]. Figure 4 shows the normalized graphene TC values as a function of concentration in armchair and zigzag directions, and the maximum TC reduction for graphene can be made between 25% and 75% of the atomic concentrations.
The explanation of the almost “parabolic” shape of TC in Figure 4 can be found in an earlier important relation derived by Klemens [20]. As it is commonly accepted, the longer the phonon mean free path is, the larger TC would be. Klemens found that, for isotope scattering, the mean free path is proportional to , where and is the temperature. Ci and Mi represent the concentration and the mass of the constituent isotope atoms, respectively. Thus, the mean free path directly has to do with the g factor, which is the mass variation of isotope atoms. In our simulation, factor reaches its maximum for around 50% of C13, so there we have the minimum phonon mean free path and the minimum TC. As atomic concentration approaches to zero or 100%, becomes smaller, phonon free path is larger, and TC increases. It will not get to infinity because of the Umklapp and other phonon scattering mechanisms.
Among various methods that modulate the TC of graphene, the isotope-doping method provides an efficient way to improve thermoelectric efficiency, since isotope atoms do not change the electronic structure of graphene, the electric conductivity, and the Seebeck coefficient remains the same after the doping. As our simulations suggest, in armchair and zigzag directions can be reduced by up to 80%, which means an improvement of the ZT by a factor of 5 since ZT is inversely proportional to TC. Recently, two important simulation studies were carried out on ZT of armchair [21] and zigzag [22] graphene nanoribbons (GNRs). In both works, the authors looked for defect and disordered structures that minimize TC while keep electric conductivity less changed, thus maximizing the ZT. Compared with this, the isotope-doping method appears to be another attractive method that would promote ZT even further for GNRs when it is combined with those methods [21, 22]. We believe our results can motivate new experiments to check the efficiency of isotope-doping method for thermoelectric application.
4. Conclusion
In summary, we have studied isotope effect on TC of isolated graphene with equilibrium MD methods. Our simulation results suggest that TC of graphene can be effectively reduced by up to 80% in armchair and zigzag directions for isotope concentrations as low as 25%. The phenomenon that mass defect can reduce TC is explained with the relation between phonon mean free path and mass variation of the isotope mixtures [20]. Finally, our study shows that doping graphene with carbon isotope C13 could be a practical way to reduce TC without changing its electric property, thus promoting thermoelectric coefficient.
Acknowledgments
The authors are grateful to support from Lockheed Martin. G. Lee acknowledges support from SWAN, and A. F. Fonseca acknowledges partial support from the Brazilian agency CNPq. We thank Rodney Ruoff for suggesting to study isotope effect in graphene. The authors also appreciate discussions with Davide Donadio and Giulia Galli.