I. INTRODUCTION
Graphene, according to its unique thermal and mechanical properties, Reference Balandin, Ghosh, Bao, Calizo, Teweldebrhan, Miao and Lau1 is a promising material for applications such as heat spreading material in very large-scale integration, Reference Grosse, Bae, Lian, Pop and King2 thermal interface materials, Reference Wang, Gong, Pei and Xu3,Reference Shahil and Balandin4 and so on. In experiments, the room-temperature thermal conductivity κ was measured up to 5300 W/mK for a monolayer graphene by inspecting the dependence of the frequency of Raman G peak on the excitation laser power. Reference Balandin, Ghosh, Bao, Calizo, Teweldebrhan, Miao and Lau1 This value is much higher than that of the best bulk thermal conductor, diamond (with κ in the range of 1000–2200 W/mK at room temperature), Reference Shinde and Goela5 and even higher than carbon nanotubes with κ = 3500 W/mK at room temperature. Reference Pop, Mann, Wang, Goodson and Dai6 However, large-scale single-crystalline graphene is very difficult to be produced. Reference Yakobson and Ding7 Although 30-inch continuous graphene films can be developed by the chemical vapor deposition (CVD) technique, Reference Bae, Kim, Lee, Xu, Park, Zheng, Balakrishnan, Lei, Ri Kim, Song, Kim, Kim, Ozyilmaz, Ahn, Hong and Iijima8 the polycrystalline nature Reference Huang, Ruiz-Vargas, van der Zande, Whitney, Levendorf, Kevek, Garg, Alden, Hustedt, Zhu, Park, McEuen and Muller9 of these films poses a question as to how the thermal transport properties of graphene is modified by the presence of lattice imperfections, especially the grain boundaries (GBs).
The thermal conductivity of graphene at room temperature was predicted theoretically with values ranging from a few hundreds to 8000 W/mK, depending on the interatomic potentials and simulation methods adopted. Reference Nika, Askerov and Balandin10–Reference Evans, Hu and Keblinski14 On the other hand, experimental measurements suggested values approximately from 2000 to 5000 W/mK, Reference Balandin and Nika15,Reference Balandin16 where the material actually features a polycrystalline nature with grain size of a few to hundreds of micrometers. Reference Balandin, Ghosh, Bao, Calizo, Teweldebrhan, Miao and Lau1,Reference Huang, Ruiz-Vargas, van der Zande, Whitney, Levendorf, Kevek, Garg, Alden, Hustedt, Zhu, Park, McEuen and Muller9,Reference Geng, Wu, Guo, Huang, Xue, Chen, Yu, Jiang, Hu and Liu17 It is not yet clear how the polygrain nature could modify the intrinsic thermal transport properties of graphene. In the literature, the GB was known to play a key role in determining the κ for semiconductors. Superlattice structures could reduce the thermal conductivity of the material to values below its alloy limit in the mechanism of minimum thermal conductivity, where the combined role of the particle and wave nature of phonon arises. Reference Garg and Chen18 For example, κ of the GaAs/AlAs superlattice was reported to be reduced by about one order of magnitude compared to the bulk conductivity of GaAs, Reference Capinski and Maris19 and this significant reduction was observed experimentally and also analyzed by several models. Reference Venkatasubramanian20,Reference Chen21 With this evidence in mind, it is thus interesting and of practical importance to assess the effect of GBs on thermal conduction in graphene.
Several studies have explored the thermal transport properties of GBs in graphene. For example, Bagri et al. Reference Bagri, Kim, Ruoff and Shenoy22 studied the thermal transport across twin GBs in polycrystalline graphene using molecular dynamics (MD) simulations and found that the boundary conductance is significantly higher than that of other thermoelectric interfaces. Cao and Qu Reference Cao and Qu23 have studied the temperature, grain size, and alignment dependence of κ for graphene with parallel GBs. A more recent study explored the effect of grain size, or dislocation density in the GB, more quantitatively, and showed that κ is not significantly reduced at room temperature if the grain size reaches several micrometers. Reference Serov, Ong and Pop24 However, all these studies were focused on graphene with parallel GBs, which obviously deviated from the structure of polycrystalline graphene identified in experiments. Reference Huang, Ruiz-Vargas, van der Zande, Whitney, Levendorf, Kevek, Garg, Alden, Hustedt, Zhu, Park, McEuen and Muller9 In this work, we will study κ of polycrystalline graphene consisting of hexagonal grains as a more realistic model, and the results are expected to be feasible to compare with experimental measurements. Reference Geng, Wu, Guo, Huang, Xue, Chen, Yu, Jiang, Hu and Liu17,Reference Kotakoski and Meyer25–Reference Wu, Fan, Speller, Creeth, Sadowski, He, Robertson, Allen and Warner28
MD simulations are usually limited by the system size it could model. To extend our discussion to the material level, an effective medium theory (EMT) where the polycrystalline graphene is considered as a two-phase composite was also used to predict thermal transport properties up to the micrometer scale. Reference Hasselman and Johnson29,Reference Nan, Birringer, Clarke and Gleiter30 This approach was widely used to identify thermal conductivities of material interfaces and nanocomposites. For example, Nan et al. Reference Nan, Birringer, Clarke and Gleiter30,Reference Hao, Fang and Xu31 introduced a so-called Kapitza radius for the interfacial thermal resistance by applying the Maxwell–Garnett EMT to predict thermal properties of nanocomposites, which works for both 2D and 3D composites. Hao et al. Reference Hao, Fang and Xu31,Reference Hao, Fang and Xu32 studied the effective thermal conductivity of graphene with defects, as well as crystalline Si–Ge nanocomposites using EMT by considering the interface between Si and Ge as an individual phase. We followed this approach in this work to analyze our results.
The paper is arranged as follows. After Sec. I, which presents the background and motivations, Sec. II introduces the atomic structure of polycrystalline graphene and the methodological details we used in MD simulations. Then in Sec. III, the dependence of κ for polycrystalline graphene on the grain size, alignment, and temperature is described and analyzed. In Sec. IV, the determination of κ for polycrystalline and the major factors affecting κ are discussed, as well as the temperature dependence, through the microscopic mechanisms of thermal conduction. Finally, main conclusions are summarized in Sec. V.
II. MODELS AND METHODS
A. The atomic structures of polycrystalline graphene
Recently, films of uniform hexagonal graphene flakes were synthesized on liquid copper surfaces, which held great promises in growing large-scale graphene films with uniform and well-controlled grain size and shape. Reference Geng, Wu, Guo, Huang, Xue, Chen, Yu, Jiang, Hu and Liu17,Reference Wu, Fan, Speller, Creeth, Sadowski, He, Robertson, Allen and Warner28 In this work, we created model structures of polycrystalline graphene from hexagonal grains following this concept, instead of randomly patching graphene flakes with irregular shapes into a 2D membrane. Reference Huang, Ruiz-Vargas, van der Zande, Whitney, Levendorf, Kevek, Garg, Alden, Hustedt, Zhu, Park, McEuen and Muller9,Reference Kotakoski and Meyer25 These graphene hexagonal flakes were defined by its grain size L and orientation angle θ, as denoted in Fig. 1. In our work, the grain size L is defined as the edge length of a hexagonal grain, and the orientation angle θ was defined as the angle between zigzag direction and X-axis in our coordinate system. Considering the hexagonal symmetry of the graphene lattice, a modulated value of θ between 0 and 30° was used, i.e., θ and 60° − θ are equivalent. Thus, a whole polycrystalline model can be denoted by a triple array of angles (α, β, and γ) that defines the orientation alignment between neighboring grains. To describe the orientation of the whole polycrystalline material, we defined a relative orientation angle (ROA) θr = (|α − β| + |α − γ| + |β − γ|)/3. We can denote the misorientation angle of the two crystalline dominions for a GB in graphene, Reference Yazyev and Louie33 and here the ROA is the mean value of the three orientational angles in polycrystalline graphene, so we can use ROA to describe the grain alignment of polycrystalline graphene. Successive procedures to optimize the GBs were carried out further by adding, removing carbon atoms, and adjusting local bonding topologies, to ensure the relative energy of each carbon atom to the value in bulk graphene less than 1 eV, corresponding to an energetically well-defined dislocation array in the GB.
(a) The atomic structure of polycrystalline graphene where hexagonal grains of graphene with size L and different orientations are patched into a super cell. Periodic boundary conditions are applied in both length and width directions. The lattice orientation of grains is denoted as θ, which is (0°, 15°, 30°) in this illustration. (b) The GBs are built up by pentagons, hexagons, and heptagons, and the typical atomic structure and intersection are plotted.
In principle, this approach could be used to create polycrystalline graphene with any grain size, shape, and ROA. But here, without loss of the generality, we focused on six polycrystalline graphene structures with hexagonal grains. Their angle triplet (α, β, and γ) are (0°, 5°, 10°), (0°, 5°, 15°), (0°, 5°, 20°), (0°, 10°, 20°), (0°, 10°, 25°), (0°, 15°, 30°) as denoted in Fig. 1(a), and the corresponding ROAs are 5, 10, 13.33, 13.33, 16.67, and 20°, respectively. Although not investigated in this work, other combinations of grain orientation angles are expected not to change the physics under discussion here qualitatively. The grain size effect for polycrystalline graphene was characterized by considering models with L = 1, 2, 3, and 5 nm. It should be noted that graphene with grain size smaller than 1 nm that we investigated is not discussed in this paper, as the high density of defects therein lead to diminishing difference between the grain and GB, and thus ambiguous definitions of their sizes.
B. Computational methods
To compute the thermal conductivity κ from MD simulations, we used the Green–Kubo method based on the linear response theory, Reference Zhang, Lee, Fonseca, Borders and Cho12,Reference Hao, Fang and Xu31 which applies for systems in thermal equilibrium where heat flux fluctuates around zero. κ could thus be expressed as an integration of the heat flux operator multiplied by a prefactor
$${{\rm{\kappa }}_{xy}} = {1 \over {V{k_{\rm{B}}}{T^{\rm{2}}}}}\int\limits_0^{{{\rm{\tau }}_{\rm{c}}}} {\left\langle {{{\bf{J}}_x}({\rm{\tau }}) \cdot {{\bf{J}}_y}(0)} \right\rangle {\rm{d\tau }}\quad {\rm{,}}}$$
where T is the temperature of system, k B is the Boltzmann constant, and V is the system volume that is defined here as the area of polycrystalline graphene multiplied by its nominal thickness (the inter layer distance of graphite, 3.4 Å). The upper limit of time integration τc needs to be long enough, so that the current–current correlation function decays to zero. Reference Schelling, Phillpot and Keblinski34 J x and J y are the heat current operators in the x and y directions, and the angular bracket represents the ensemble average, Reference Hardy35 namely the heat flux autocorrelation function (HFACF). The heat flux J of the system was computed from the expression J = (Σ i e i v i − Σ i S i v i )/V, where e i , v i , and S i are the total energy, velocity vector, and stress tensor of each atom i, respectively. We first integrated HFACF with an integration time τ following Eq. (1) and then obtained the relation between κ and τ. The decorrelation time for the heat flux is typically on the scale 10 ps for our models, and thus converged results for κ could be extracted when τ > τc in the simulations. The advantage of Green–Kubo method, compared to the nonequilibrium molecular dynamics (NEMD) simulations with thermal gradient built up across the material, is that it does not require additional perturbation to the equilibrium that may give rise to nonlinearity and has weaker dependence on the size of a simulated system. Due to the hexagonal symmetry of graphene lattice, the thermal conductivity of polycrystalline graphene is isotropic. We thus evaluate κ as the mean value of κ xx and κ yy that may differ in the finite system under simulation though, where the maximum difference between κ xx and κ yy in our simulations is less than 0.1κ.
All simulations were performed using the large-scale atomic/molecular massively parallel simulator. Reference Plimpton36 Periodic boundary conditions were applied to a 2D supercell of polycrystalline graphene, and to minimize the size effect, we used almost the same supercell size in the MD simulations (18.0 × 31.20 nm2 for polycrystalline graphene with L = 1, 2, 3 nm, 15.0 × 26.00 nm2 for polycrystalline graphene with L = 5 nm, and 17.9 × 31.20 nm2 for pristine graphene). The optimized Tersoff potential, Reference Lindsay and Broido37 which has been demonstrated to provide an excellent prediction for the phonon dispersion graphene in comparison with experimental measurements, was used for the interatomic interactions between carbon atoms. This method was widely used to predict thermal and mechanical properties of graphene and related materials.
In the MD simulations, the atomic structures of polycrystalline graphene were firstly equilibrated at ambient condition (temperature T = 300 K, where the quantum correction is negligible Reference Savin, Kosevich and Cantarero38 ) under a Nosé–Hoover thermostat for 200 ps, in which the time step is 0.2 fs. Then the structure was further equilibrated in a microcanonical ensemble for 50 ps before the thermal conductivity was calculated. In the equilibrium Green–Kubo calculations, the system was simulated in a microcanonical ensemble as well. The atomic positions and velocities were collected along the simulations to evaluate the heat flux and its autocorrelation functions. The thermal conductivity was finally obtained by following the Green–Kubo formula and averaging over eight independently sampled runs for each structure.
III. RESULTS
A. The dependence of κ on grain size and orientation
We first explored the grain size dependence at T = 300 K. Data points plotted in Fig. 2(a) are our calculation results for κ of polycrystalline graphene with grain size ranging from 1 to 5 nm, but with the same orientation (0°, 5°, 10°). The results show that κ increases for larger size of grains. The value for κ is 60.76 W/mK for L = 1 nm and 197.69 W/mK for L = 5 nm. The data can be well fitted into a linear function as κ = 26.67 + 35.24L (nm) W/mK. For comparison, the thermal conductivity of graphene was calculated to be 892.68 W/mK following the same simulation procedure. Thus, the incorporation of nanosized grains reduces the thermal conductivity of graphene remarkably, by amount up to 93% (for L = 1 nm). Further calculations for models with different ROAs yield similar results.
(a) Thermal conductivity of polycrystalline graphene as a function of the grain size L for the mismatch angle triplet (0°, 5°, 10°). (b) Thermal conductivity of polycrystalline graphene as a function of the ROA θr for different grain size, which are 5, 10, 13.33, 13.33, 16.67, and 20° from left to right, respectively.
In experiments, GBs with different orientational mismatch were identified, Reference Huang, Ruiz-Vargas, van der Zande, Whitney, Levendorf, Kevek, Garg, Alden, Hustedt, Zhu, Park, McEuen and Muller9 which corresponds to a certain density level of dislocations or defects in the GB. Considering that the GB as a linear array of edge dislocations, the density of consecutive dislocations is thus given by the Frank formula, ρ = 1/d = 2sin(θ/2)/|b|, where b is the Burgers vector and d is the distance between two neighboring dislocations. Reference Song, Artyukhov, Yakobson and Xu27 To quantify the effect of orientation angles, we calculate values of κ for polycrystalline graphene with different ROAs. The results are shown in Fig. 2(b), which suggest that κ decreases as the ROA θ r increases, and the trend is the same for different grain sizes. For L = 5 nm, κ is 197.69 W/mK when the angle triplet is (0°, 5°, 10°) and 126.09 W/mK for the angle triplet of (0°, 15°, 30°), which is reduced by 36% as θr increases from 5 to 20°.
B. Heat flux distribution and the effective width of the GB
We know that effective thermal conductivity of the polycrystalline graphene κ increases with the grain size L and decreases with the ROA θr, which is closely related to the atomic structure of polycrystalline graphene. The structure of a GB in our model comprises a series of topological defects, e.g., pentagon–heptagon (5|7) pairs, and the GB tends to break into an array of dislocation as shown in Fig. 1(b). The properties of polycrystalline graphene thus will rely on the characteristics of the GB (size, the density of defects, etc.) The role of dislocation in modifying mechanical and thermal properties was well identified for bulk materials, Reference Hull and Bacon39,Reference Ziman40 as material discontinuity in crystals. From the continuum theory where the dislocation is seen as a line in 3D or a point in 2D, one can define the effective width of the dislocation as Reference Ziman40
$$w = {{\rm{\gamma }}^2}{{\rm{\bf{b}}}^2}{{\rm{{\bf{q}}}} \mathord{\left/ {\vphantom {q 2}} \right. \kern-\nulldelimiterspace} 2}\∼\,{\rm{\pi }}{{\rm{\gamma }}^2}\,{{{\rm{\bf{b}}^2}} \mathord{\left/ {\vphantom {{{b^2}} {\rm{\lambda }}}} \right. \kern-\nulldelimiterspace} {\rm{\lambda }}}\quad ,$$
where the length of the Burgers vector b is usually about one lattice constant, q and λ are the wave vector and wave length of the phonon, respectively, and the parameter γ is a dimensionless parameter that describes the relation between dilatation and temperature when the temperature is low. Considering the dislocation core, i.e., the region where the atomic displacements are so large that elasticity theory cannot be adequately applied, one can follow the theory of Rayleigh scattering and estimate the scattering width as Reference Ziman40
$$w = a{\left( {{{\Delta D} \mathord{\left/ {\vphantom {{\Delta D} D}} \right. \kern-\nulldelimiterspace} D}} \right)^2}{\left( {\rm{\bf{q}}}a} \right)^3}\quad ,$$
where a is the radius of the scattering region and where the material has a constantly shifted density (from D to D + ΔD). Although the assumptions in Eqs. (2) and (3) are very crude, it reveals some basic aspects such as the effective width of dislocation is finite and is quite small for long-wave phonons.
If a crystal contains a GB of dislocations, a propagating phonon attempting to cross the GB may find itself in a region where it can no longer proceed in the same direction with the same velocity and energy. That is to say, just as a single dislocation has a finite scattering width, there must also be a bounded region with an effective width w. Inside this region, the crystal structure cannot be treated as a perfect one. To find more details about this bounded region as well as the definition of w, and to obtain more insights into the grain size and orientational dependence of thermal conductivity for polycrystalline graphene, we plot the spatial distribution of heat flux within the whole polycrystal in Fig. 3, which displays the heat flux of polycrystalline graphene for the specified angle triplet (0°, 5°, 10°). This was done by calculating atomic heat flux in NEMD simulations that was averaged over 10 ps in the steady state. The atomic heat flux is defined from the expression J i = e i v i − S i v i . The results clearly identified strong scattering of the heat flux in the region adjacent to the GBs. However, inside the grains and away from the GBs, the spatial distribution becomes close to that in a single crystal of graphene. When the grain size increases, the portion of scattered region becomes reduced. Thus, we can consider the polycrystalline graphene as a composite consisting of two phases—the pristine grains and GBs that feature quite different material properties. Furthermore, from the heat flux distribution, we could conclude visually that the scattering region spans in few lattice constants, and quantitative determination of w will be introduced below with respect to a predictable model for thermal conductivity of polycrystalline graphene. The interfacial thermal conductivity for such a GB is ∼1500 MW/m2K. This conductivity is higher than that value (10–500 MW/m2K) quantified for typical nonbonding interfaces, such as those between graphene and the supporting substrate or matrix material where it is embedded. Reference Xu and Buehler41,Reference Chang, Nair and Buehler42
Spatial distribution of the heat flux in pristine and polycrystalline graphene for the specified angle triplet (0°, 5°, 10°). (a) Pristine graphene. (b) Polycrystalline graphene with single GB across the periodic boundary, and the heat flux is parallel to the GB. (c–f) Polycrystalline graphene with grain size of 1, 2, 3, and 5 nm, respectively. The length of scale bar is 0.49, 0.15, 0.14, 0.20, 0.31, and 0.43 nm, respectively.
C. EMT for the prediction of κ
A validated effective model could be used to predict the thermal transport properties of polycrystalline graphene with the grain size much larger than the one our MD simulations are capable to capture. The key to establish such a model is to well define the model and fit parameters therein. In our polycrystalline graphene models, there are three GBs with different misorientation angles for each sample, and thus the effective width w and thermal conductivity κGB for the GB phase should be different with respect to these GBs. However, we noted in the aforementioned simulation results [Fig. 2(b)] and previous reports Reference Cao and Qu23 that there is only a small difference in κGB, e.g., κGB for a zigzag GB with orientation angle 21.7° is only 15% less than that of a zigzag GB with orientation angle 5.5°. Moreover, by assumption, these three GBs have the same volume fraction and distribute homogeneously in our polycrystalline graphene model although the density of dislocations is different, which can be seen from Fig. 1. Thus for simplicity, we divided the polycrystalline graphene into two different phases, i.e., the bulk phase with the properties as pristine graphene, and a GB phase with no microstructures specified. In other words, we considered the material as 2D composite with GBs embedded in bulk graphene as inclusions. The overall thermal resistance of the composite, κ, could thus be predicted as Reference Hao, Fang and Xu31
$${1 \mathord{\left/ {\vphantom {1 {\rm{\kappa }}}} \right. \kern-\nulldelimiterspace} {\rm{\kappa }}} = {1 \mathord{\left/ {\vphantom {1 {{{\rm{\kappa }}_{{\rm{GB}}}}}}} \right. \kern-\nulldelimiterspace} {{{\rm{\kappa }}_{{\rm{GB}}}}}} + {1 \mathord{\left/ {\vphantom {1 {{{\rm{\kappa }}_{\rm{G}}}}}} \right. \kern-\nulldelimiterspace} {{{\rm{\kappa }}_{\rm{G}}}}}\quad ,$$
where κ, κGB, and κG are the predicted effective thermal conductivities of the composite, the GB, and the bulk phase (pristine graphene), respectively, and f GB is the volume fraction of the GB phase. In EMT for composites, inclusions affect overall thermal transport behaviors through phonon scattering, Reference Nan, Birringer, Clarke and Gleiter30 with the same mechanism as here for polycrystalline graphene with GBs in two dimension. To use Eq. (4) for quantitative prediction, we need to first determine the effective width of the GB so f GB could be defined, through f GB = (2√3wL − w 2)/3L 2.
But how could we determine w from the MD simulation results? As mentioned previously, we defined w ≈ 0.7 nm visually by inspecting the spatial distribution of heat flux and measuring the size of scattering region in Fig. 3. On the other hand, the value thus determined could be verified by substituting the value and f GB it defines into Eq. (4) to see if the simulation results for κ could be well fitted. In statistics, the coefficient of determination, Reference Carpenter43 denoted as R 2, indicates how well data points fit a line or curve. In our fitting, R 2 is up to 0.99, so we can say that Eq. (4) fits the simulation results very well. From now on, without loss of accuracy, we define the effective width w constantly as 0.7 nm in our EMT model. Fitting the simulation results in Fig. 4(a) gives
$${\rm{\kappa }} = {{719.97} \mathord{\left/ {\vphantom {{719.97} {\left( {1 + 16.81{f_{{\rm{GB}}}}} \right)}}} \right. \kern-\nulldelimiterspace} {\left( {1 + 16.81{f_{{\rm{GB}}}}} \right)}}{{\rm{W}} \mathord{\left/ {\vphantom {{\rm{W}} {{\rm{mK}}\quad }}} \right. \kern-\nulldelimiterspace} {{\rm{mK}}\quad }}.$$
(a) Thermal conductivity of polycrystalline graphene as a function of the volume fraction of GB phase f GB for the specified angle triplet (0°, 5°, 10°). We fitted the results by using the effective medium model and the result shows κ = 719.97/(1 + 16.81f GB) W/mK. (b) The ratio κ/κ∞ as a function of the grain size for the specified angle triplet (0°, 5°, 10°), κ∞ is the value when the grain size is infinity. The black line is our prediction and the red dashed line is data taken from Serov et al. Reference Serov, Ong and Pop24
Based on Eq. (5), we fitted the parameters in the EMT model. We can then predict κ for polycrystalline graphene with arbitrary grain size, as the physics is not changed. The result shows that as L increases, κ increases and converges to a limit κ = 719.97 W/mK that corresponds to the value for pristine graphene. The L-dependence of κ for polycrystalline graphene with the specified angle triplet (0°, 5°, 10°) is plotted in Fig. 4(b), as well as the calculation results by Serov et al. Reference Serov, Ong and Pop24 for polycrystalline graphene with parallel GBs. The values were both normalized by the thermal conductivity of pristine graphene, i.e., the limit at infinitely large grain size. The comparison shows that although the trend is in general the same, the predicted effective thermal conductivity for our polycrystalline model with hexagonal grains is higher than the superlattice-like model explored by Serov. Reference Serov, Ong and Pop24 The difference may thus be attributed to the contrastive structure (the density of dislocations, orientation angle, etc.) and scattering of phonons propagating across the GBs perpendicularly in their work. Our results for the more realistic polycrystalline models suggest that at room temperature, the thermal conductivity of polycrystalline graphene with grain sizes larger than several hundred of nanometers shows no significant reduction in comparison with the pristine single crystal.
D. The temperature dependence of κ
The temperature dependence of κ for graphene and polycrystalline graphene calculated in our studies is shown in Fig. 5. For comparison, we also plot κ measured for suspended CVD-grown graphene at different temperatures (by Cai et al. Reference Cai, Moore, Zhu, Li, Chen, Shi and Ruoff44 ) and calculated for polycrystals with parallel GBs (by Cao and Qu Reference Cao and Qu23 ) in Fig. 5. For all the results, we found that κ monotonically decreases with the temperature that agrees well with other theoretical studies, for example Ref. Reference Wei, Wu, Yin, Shi, Yang and Dresselhaus45. The results in Fig. 5 suggest a temperature dependence of κ ∼ T −α for T from 200 to 1000 K for all the cases, and we can define an exponent α to capture this scaling behavior. κ for pristine graphene calculated in this work shows a scaling behavior of ∼T −0.96, while κ for experimental measurements is ∼T −1.03. There is thus ∼7% difference in the exponents α determined from our simulations compared with experimental measurements. The value of α is reduced when the polycrystalline nature is introduced. For the graphene with parallel GBs, Reference Cao and Qu23 α ∼ 0.79, and for our polycrystalline graphene models with hexagonal grains, the exponents α vary from 0.22 to 0.32 with increase in grain size from 1 to 5 nm, which are both much smaller in comparison with α ∼ 0.96 for the pristine graphene. From the scaling relation κ ∼ T −α, we conclude that the temperature dependence of κ is significantly reduced for polycrystalline graphene, which is evidenced in both the reduced exponent α compared to the pristine graphene and enhanced α at large grain size.
Thermal conductivity of pristine and polycrystalline graphene plotted with the temperature and grain size, where the polycrystalline graphene features a specified angle triplet (0°, 5°, 10°). Experimental and simulation results for the temperature dependence found in the literature Reference Cao and Qu23,Reference Cai, Moore, Zhu, Li, Chen, Shi and Ruoff44 are also plotted here for comparison.
IV. DISCUSSION
A. Phonon scattering at GB and the effective width
Considering the graphene with an armchair GB with orientation angle 15.4° as shown in Fig. 6(a), the propagating phonon will experience reflection and refraction when it transports across the GB, and its group velocity will have a mismatch between the two crystalline domains on each side of the GB. This mismatch suggests additional phonon scattering sources in addition to the lattice distortion and disorder at the GB. To identify the scattering effect on the thermal transport quantitatively, we calculate κ for the graphene model with this infinite GB (across the periodic boundary) using the Green–Kubo method described above. The results show that κ⊥ = 122.03 W/mK and κ|| = 249.68 W/mK, where κ⊥ and κ|| are thermal conductivities of graphene with GBs when the direction is perpendicular to GB and when the direction is parallel to GB, respectively. To further analyze these results, we divided the polycrystalline graphene with single GB into two phases, i.e., the bulk graphene and GB phases. For thermal transport perpendicular or parallel to the GB, we assumed that the effective width of GB w is w ⊥ and w ||, and the thermal conductivity of GB as κGB⊥ and κGB||, respectively. By applying Eq. (4) we obtain
$$({{1 \mathord{\left/ {\vphantom {1 {\rm{\kappa }}}} \right. \kern-\nulldelimiterspace} {\rm{\kappa }}}_{{\rm{GB}} \bot }} - {1 \mathord{\left/ {\vphantom {1 {{{\rm{\kappa }}_{\rm{G}}}}}} \right. \kern-\nulldelimiterspace} {{{\rm{\kappa }}_{\rm{G}}}}}){w_ \bot } = 0.07\,{\rm{nm}} \times {\rm{mK/W}}\quad {\rm{,}}$$
$$({{1 \mathord{\left/ {\vphantom {1 {\rm{\kappa }}}} \right. \kern-\nulldelimiterspace} {\rm{\kappa }}}_{{\rm{GB}}{\parallel} }} - {1 \mathord{\left/ {\vphantom {1 {{{\rm{\kappa }}_{\rm{G}}}}}} \right. \kern-\nulldelimiterspace} {{{\rm{\kappa }}_{\rm{G}}}}}){w_{\parallel} } = 0.03\,{\rm{nm}} \times {\rm{mK/W}}\quad .$$
(a) The atomic structure of graphene with parallel GBs. The mismatch angle is 15.4°. (b) The spatial distribution of strain field in the graphene with an armchair GB, calculated using the data from our MD simulations. (c) The shear strain field induced by a dislocation calculated using Eq. (8), Reference Ziman40 plotted along lines 1 and 2 as denoted in panel (b).
As seen from the above relation, the definition of effective width w couples with thermal conductivities of the GB, which makes the determination of these two parameters from simulation results nontrivial. To solve this problem, we can define that GB has an intrinsic thermal conductivity κGB and assume it as the mean value of κGB⊥ and κGB||, i.e., κGB = (κGB|| + κGB⊥)/2, This crude treatment, however, is convenient in estimating the effective properties of a realistic polycrystalline graphene with two-dimensional patches of grains. For a single GB, the effective width w should be defined uniquely and thus w ⊥ should equal to w || . From the assumption above and calculations following Eqs. (6) and (7), we have w = 0.59 nm, which conveys the thickness of a single dislocation. In our previous discussion based on the spatial distribution of heat flux in the polycrystalline graphene, the effective width was found to be 0.7 nm, which is close to the value w = 0.59 nm extracted here by considering the fact that both GBs consist of one array of dislocations. The agreement suggests that the assumption we made is reasonable and the definition of GB as an individual phase with well-defined structural and thermal transport properties is feasible in discussing the effective thermal conductivity of polycrystalline graphene.
To analyze our results through the EMT, we have introduced the effective width w to calculate the volume fraction of different phases. Naturally, a problem arises that how to determine the effective width of a GB. In the preceding section, we have obtained an effective width from the heat flux field in MD simulations, but now we discuss the definition of effective width directly from the strain or stress field of graphene with GBs. Ziman Reference Ziman40 argued that for a GB of orientation angle θ tends to break into an array of dislocations, spaced by d = b/2sin(θ/2) = b/θ. The effect of such a GB can be found in the same way as for a single dislocation. It is known that, a GB, spaced by d apart in the y direction, gives rise to a shear strain at the point (x, y)
$${{\rm{\varepsilon }}_{xy}}{ \,= }{{\rm{{\bf{b}}}} \over {{2\pi (1} - {\rm{\mu }}{)}}}{Re}\left[ {{{{{\pi }^{2}}x} \over {{d^{2}}}}{{\rm{cosec}}}{{{\rm{h}}}^{2}}{{{\pi (}x{ + }iy{)}} \over d}} \right]\quad ,$$
where μ is the Poisson's ratio, and x, y are measured from the dislocation core, as indicated in Fig. 6(c). At a large distance from the GB, this shear strain can be approximately written as a function of x
$${{\rm{\varepsilon }}_{xy}}{ \,=\, }{{\rm{{\bf{b}}}} \over {{2\pi (1} - {\rm{\mu }}{)}}}{{{4}{{\pi }^{2}}x} \over {{d^{2}}}}{{e}^{ - {2\pi }x{/}d}}\quad ,$$
which decreases exponentially from the GB and can be neglected for x >> d. It should be noted that Eq. (9) was originally derived for edge dislocations. Although the out-of-plane displacement of carbon atoms in graphene could release the in-plane stress and strain, recent discussions on the mechanical properties of polycrystalline graphene validates its applicability in qualitative predictions. Reference Yakobson and Ding7,Reference Song, Artyukhov, Yakobson and Xu27,Reference Wei, Wu, Yin, Shi, Yang and Dresselhaus45 So we can measure the spatial distribution of shear strain in experiments or simulations, namely ε xy . As described above, we can define a specified value w, when x is larger than half of this value w the shear strain will be small enough to be neglected. Then we can consider the specified value w as our effective width of GB. This approach can be used to determine the effective width of a GB without referring to thermal transport measurements. The value of w thus determined could be further verified by comparing it with the value obtained from MD simulations that was discussed in the previous sections. Here we define the characteristic strain of graphene as ε i = [(r 1 + r 2 + r 3) − r 0]/r 0, where ε i is the strain of the atom i, r 1, r 2, and r 3 are the length of the three bonds with atom i, and r 0 = 1.42 nm is the carbon–carbon bond length in pristine graphene. For the graphene model with parallel GBs as shown in Fig. 6(a), we calculated the spatial distribution of strain field from our MD simulation results and plotted it in Fig. 6(b). By assuming that the strain can be neglected when it is smaller than 0.05, we measured the effective scope from the distribution of strain and obtain that the value is in the range from 0.6 to 1.0 nm, which agrees to the effective width 0.7 and 0.59 nm as extracted in the previous sections.
The definition of effective width w also depends on the density of 5|7 dislocations, which is related to the orientation angle according to the Frank formula. Reference Song, Artyukhov, Yakobson and Xu27 In the polycrystalline graphene we constructed, there are three GBs featuring different density levels of 5|7 dislocations. By neglecting the difference in κGB, Reference Cao and Qu23 these three GBs were homogenized in our EMT model. Moreover, from the results of κ⊥ = 122.03 W/mK and κ|| = 249.68 W/mK or the armchair GB with orientation angle 15.4°, we know that w should also rely on the angle between the direction of heat flux and GBs, namely χ. In our polycrystalline graphene model, χ is different for GBs at different positions as the heat flux was applied along the X-direction. However, no significant difference was found for the spatial distribution of heat flux (as shown in Fig. 3), so we excluded the effect of χ-dependence when analyzing κ.
B. The factors that affect thermal transport in polycrystalline graphene
To see what factors play key roles in defining the thermal transport properties of the polycrystalline graphene, we extend the discussion with analysis based on the Boltzmann transport equation (BTE), which was widely used to solve thermal transfer problems in material at the carrier level by neglecting phonon coherence. The general form of BTE is
$$\left( {{\partial \over {\partial t}} + {\bf{v}} \cdot \nabla } \right)f({\bf{r}},{\bf{p}},t) = {\left( {{{\partial f} \over {\partial t}}} \right)_{{\rm{scattering}}}}\quad ,$$
where r is the position, p is the momentum, t is time, v is the particle velocity, and f(r, p, t) is the single particle distribution function in the phase space at given time. The factors that affect the thermal transport properties can be discussed based on this equation. To solve the BTE [Eq. (10)], one can use the single-mode relaxation time approximation, Reference Srivastava and Kresin46 which assumes that each phonon mode is characterized by its own relax time that is independent of all other phonons in the system. With this assumption, the scattering term in Eq. (10) can be approximated as
$${\left( {{{\partial f} \over {\partial t}}} \right)_{{\rm{scattering}}}} = {{ - {f^{\rm{'}}}} \over {\rm{\tau }}}\quad ,$$
where f′ is the fluctuation of the distribution function at equilibrium and τ is the relaxation time. For the phonon scattering by 5|7 dislocations, the relaxation time is related to the defect density ρ, the group velocity v, and the effective width of dislocation w as Reference Kim and Majumdar47
$${\rm{\tau }}\∼{1 \mathord{\left/ {\vphantom {1 {{\rm{\rho }}vw}}} \right. \kern-\nulldelimiterspace} {{\rm{\rho }}vw}}\quad .$$
We can further denote the phonon-specific heat as c(ω) and the mean free path as l(ω). The thermal conductivity Reference Berman48 can thus be expressed as
$${\rm{\kappa }} = {1 \over 3}\mathop\sum\limits_{\rm{\omega }} {c({\rm{\omega }})vl({\rm{\omega }})} = {1 \over 3}\mathop\sum\limits_{\rm{\omega }} {c({\rm{\omega }}){v^2}{\rm{\tau }}({\rm{\omega }})} \quad .$$
So according to Eqs. (12) and (13), the material with a higher defect density has the shorter relaxation time, and thus the lower thermal conductivity. For symmetric GBs in graphene, the density of 5|7 dislocation is determined by the orientation angle following the Frank formula. Reference Song, Artyukhov, Yakobson and Xu27 The GB with a larger orientation angle contains a higher 5|7 dislocation density. For the polycrystalline graphene with hexagonal grains as explored in this work, the density of defects increases as the grain size is reduced or the orientation angle increases. That is to say, the polycrystalline graphene with a smaller grain size or a larger orientation angle has the lower thermal conductivity, which agrees well with the results in Fig. 2. For other types of GBs, for example those with asymmetric GB or other defects, the scattering of heat flux is expected to reduce the thermal conductivity of polycrystalline graphene in a similar way as symmetric GBs we have discussed here.
C. The temperature-dependence of effective thermal conductivity
First, we consider the temperature dependence of κ for the pristine graphene. From the theory of phonon thermal transport we know that Reference Srivastava and Kresin46 the group velocity v is determined by the phonon dispersion that does not change with temperature and the mean free path l scales with the temperature T as 1/T, Reference Berman48 and the specific heat c should depend on the temperature T, namely c(T), scales as T 3 at low temperature and approaches a constant when the temperature is close to the Debye temperature T D (i.e., 1045 K for graphene). Reference Politano, Borca, Minniti, Hinarejos, de Parga, Farias and Miranda49,Reference Pop, Varshney and Roy50 Following Eq. (13) we can see that κ ∼ T −1 if the temperature dependence of c is weak, which is very close to our simulation results κ ∼ T −0.96 and Cai's Reference Cai, Moore, Zhu, Li, Chen, Shi and Ruoff44 experimental measurements κ ∼ T −1.03.
Then for the polycrystalline graphene that is treated as a composite material here, we can assume that only two mechanisms contribute to the phonon relaxation or life time obtained from our MD simulations, and those are (i) Umklapp scattering inside the grain (τU) and (ii) the interfacial scattering at the GBs (τGB). If we assume that these two processes are independent, we can invoke Matthiessen's rule and write the inverse lifetime of the phonon mode (k, ω) a sum of contributions from each mechanism
$${1 \mathord{\left/ {\vphantom {1 {\rm{\tau }}}} \right. \kern-\nulldelimiterspace} {\rm{\tau }}}(k,{\rm{\omega }}) = {1 \mathord{\left/ {\vphantom {1 {{{\rm{\tau }}_{\rm{U}}}}}} \right. \kern-\nulldelimiterspace} {{{\rm{\tau }}_{\rm{U}}}}}(k,{\rm{\omega }}) + {1 \mathord{\left/ {\vphantom {1 {{{\rm{\tau }}_{{\rm{GB}}}}(k,{\rm{\omega }})}}} \right. \kern-\nulldelimiterspace} {{{\rm{\tau }}_{{\rm{GB}}}}(k,{\rm{\omega }})}}\quad ,$$
where k and ω are the wave vector and frequency of the phonon, respectively. As temperature increases, the Umklapp scattering is enhanced as higher energy phonons are thermally populated Reference Kim, Shi, Majumdar and McEuen51 and τU is reduced. Reference Nika, Pokatilov, Askerov and Balandin52 In contrast, the GB scattering is weakened and τGB increases. Reference Cao and Qu23 Then as predicted by Eq. (13), κ increases with the relaxation time. From the simulation results plotted in Fig. 5, we can see that κ monotonically decreases with T, so we conclude that the Umklapp scattering plays a more important role than the GB scattering. Moreover, as the grain size decreases, the defect density increases and the GB scattering could make more contribution to the total relaxation time, which is, however, still lower than that from the Umklapp scattering. That is to say, the relaxation time τ of polycrystalline graphene with a smaller grain size is reduced less than that of polycrystalline graphene with a larger grain size as temperature increases, and the temperature dependence of κ is weakened as the grain size decreases.
V. CONCLUSIONS
In summary, we calculated the effective thermal conductivity κ of polycrystalline graphene using equilibrium MD simulations. It was found that the effective κ increases with the grain size, and the results showed a linear dependence when the grain size is a few nanometers. On the other hand, κ is reduced as the mismatch angle or dislocation density of the GB increases. We also found that κ is reduced when temperature increase and the temperature dependence of κ is much weakened when the grain size decreases. The underlying phonon scattering mechanisms were discussed through phonon scattering theories, with respect to the grain size, orientation, and temperature dependence. To predict κ with large grain size, we constructed an effective medium model, where polycrystalline graphene was considered as a two-phase composite. Based on the effective model we predicted κ for grain size when grain size ranges up to 10 μm. The results showed that the thermal conductivity of polycrystalline graphene converges as the grain size increases to a few micrometers, in the typical range of grain sizes of CVD-grown graphene. These understandings and the formulated effective model with predefined GBs also shed light on engineering thermal transport processes in a material with rationally designed defective microstructures.
ACKNOWLEDGMENTS
This work was supported by the National Natural Science Foundation of China through Grant Nos. 11222217 and 11002079, Tsinghua University Initiative Scientific Research Program 2011Z02174. The computation was performed on the Explorer 100 cluster system of Tsinghua National Laboratory for Information Science and Technology.
