Heat conduction in graphene: experimental study and theoretical interpretation

We review the results of our experimental investigation of heat conduction in suspended graphene and offer a theoretical interpretation of its extremely high thermal conductivity. The direct measurements of the thermal conductivity of graphene were performed using a non-contact optical technique and special calibration procedure with bulk graphite. The measured values were in the range of ∼3000–5300 W mK−1 near room temperature and depended on the lateral dimensions of graphene flakes. We explain the enhanced thermal conductivity of graphene as compared to that of bulk graphite basal planes by the two-dimensional nature of heat conduction in graphene over the whole range of phonon frequencies. Our calculations show that the intrinsic Umklapp-limited thermal conductivity of graphene grows with the increasing dimensions of graphene flakes and can exceed that of bulk graphite when the flake size is on the order of a few micrometers. The detailed theory, which includes the phonon-mode-dependent Gruneisen parameter and takes into account phonon scattering on graphene edges and point defects, gives numerical results that are in excellent agreement with the measurements for suspended graphene. Superior thermal properties of graphene are beneficial for all proposed graphene device applications.

Abstract. We review the results of our experimental investigation of heat conduction in suspended graphene and offer a theoretical interpretation of its extremely high thermal conductivity. The direct measurements of the thermal conductivity of graphene were performed using a non-contact optical technique and special calibration procedure with bulk graphite. The measured values were in the range of ∼3000-5300 W mK −1 near room temperature and depended on the lateral dimensions of graphene flakes. We explain the enhanced thermal conductivity of graphene as compared to that of bulk graphite basal planes by the two-dimensional nature of heat conduction in graphene over the whole range of phonon frequencies. Our calculations show that the intrinsic Umklapplimited thermal conductivity of graphene grows with the increasing dimensions of graphene flakes and can exceed that of bulk graphite when the flake size is on the order of a few micrometers. The detailed theory, which includes the phonon-mode-dependent Gruneisen parameter and takes into account phonon scattering on graphene edges and point defects, gives numerical results that are in excellent agreement with the measurements for suspended graphene. Superior thermal properties of graphene are beneficial for all proposed graphene device applications.

Introduction
It was recently discovered experimentally that graphene [1,2] has an extremely high thermal conductivity [3,4]. The measurements reported by Balandin et al [3,4] were performed using the non-contact optical technique, where the local temperature rise due to laser heating was determined through the independently measured temperature coefficients of the graphene Raman G peak [5,6]. It was found that the near room-temperature (RT) thermal conductivity of partially suspended single-layer graphene is in the range K ∼ 3000-5300 W mK −1 depending on the graphene flake size. These experimental findings stimulated a body of theoretical work on the subject. Nika et al [7] performed a detailed numerical study of the lattice thermal conductivity of graphene using the phonon dispersion obtained by the valence-force field (VFF) method. The authors treated the three-phonon Umklapp scattering directly considering all phonon relaxation channels allowed by the energy and momentum conservation in the graphene two-dimensional (2-D) Brillouin zone (BZ) [7]. Jiang et al [8] calculated the thermal conductance of graphene in the pure ballistic limit obtaining a high value, which translates to a thermal conductivity in excess of ∼6600 W mK −1 . A higher thermal conductivity is expected for the ballistic regime when no scattering is included. Lan et al [9] determined the thermal conductivity of graphene nanoribbons by combining the tight-binding approach and the phonon non-equilibrium Green's function method. The authors found the thermal conductivity K = 3410 W mK −1 [9], which is clearly above the bulk graphite limit of 2000 W mK −1 and in agreement with the first experiments [3,4]. Strong edge and flake size effects were also revealed by the numerical data in line with the experiments [3,4].
In this paper, we review our experimental results, provide details of the data extraction procedure and offer a physical interpretation of heat conduction in graphene. Our theoretical model allowed us to explain the higher values of thermal conductivity in suspended graphene as compared to that of basal planes of bulk graphite. The rest of the paper is organized as follows. Section 2 describes the method developed for measuring the thermal conductivity of the suspended graphene. In section 3, we provide details of the measurements and calibration procedure. Section 4 outlines the formal theory of the thermal conductivity. It is an accurate approach but very challenging computationally. In order to treat the three-phonon Umklapp processes exactly, we propose the scattering diagram method. In section 5, we offer a simple model of the thermal conductivity of graphene. This model is less accurate but allows one to clearly see the physical reasons behind the higher thermal conductivity of graphene as compared 3 to that of basal planes of bulk graphite. It also explains the inherent importance of the lateral size effects on the thermal conductivity of graphene flakes.

Experimental investigation of heat conduction in graphene
Since none of the conventional methods for the thermal conductivity measurement work for graphene, we developed our own non-contact optical approach. We took advantage of the fact that graphene has distinctive signatures in Raman spectra with a clear G peak and 2D band [10]- [13]. Moreover, we also found that the G peak of graphene's Raman spectra exhibits strong temperature dependence [5,6]. The latter means that the shift in the position of the G peak in response to laser heating can be used for measuring the local temperature rise. The correlation between the temperature rise and amount of power dissipated in graphene, for a sample with given geometry and proper heat sinks, can give the value of the thermal conductivity K (see the schematic of the experiment in figure 1(a)). Even a small amount of power dissipated in graphene can be sufficient for inducing a measurable shift in the G peak position due to the extremely small thickness of the material-one atomic layer. The suspended portion of graphene served several essential functions for (i) accurately determining the amount of power absorbed by graphene through the calibration procedure, (ii) forming 2-D in-plane heat front propagating toward the heat sinks and (iii) reducing the thermal coupling to the substrate through the increased micro and nanoscale corrugations (see figure 1(b)).
The first step in the measurements was determining the temperature coefficient χ G for the G peak. To accomplish this task the laser excitation power was kept at a minimal level and the temperature of the graphene flake was changed externally through the hot-cold cell [5,6]. The results of our measurement of the temperature coefficients were verified by measuring the coefficients for bulk graphite, which were in excellent agreement with the previously published data for graphite [5,6]. After the change in the temperature T has been correlated with the change in the ω peak position, the micro-Raman spectrometer can be used as a thermometer. During the measurement of the thermal conductivity, the excitation power is intentionally increased to induce the local heating. The local temperature rise is determined through the expression T = ω G /χ G . It is important to mention here that the measurement technique is steady state. Each data point in the thermal conductivity measurement, i.e. in recording the G peak position as a function of the excitation power, takes sufficient time (several minutes) for achieving the steady state. The energy deposited by the laser light to the electron gas in graphene is transferred to phonons very fast. The time constant for the energy transfer from the electrons to acoustic phonons in graphene is on the order of several picoseconds [14]- [16]. Thus, for the large graphene flakes utilized in our experiments (tens of microns) the changes in the induced hot spot due to the finite thermalization time are small and can be neglected. From the other side, our measurement time was small compared to the hours that are required in order to induce damage or surface modification in graphene by laser light [17].
The long graphene flakes for these measurements were produced using the standard technique of mechanical exfoliation of bulk Kish and highly oriented pyrolytic graphite (HOPG) [1,2]. The trenches were fabricated using reactive ion etching. The width of these trenches ranged from 1 to 5 µm with a nominal depth of 300 nm. In the first set of measurements, we selected graphene flakes of approximately rectangular shape connected to large graphitic pieces, which acted as heat sinks. The rectangular shape was selected in order to use a simple data extraction procedure based on the 1-D heat diffusion equation. These graphitic pieces  were at a distance of a few micrometers from the trench edges to ensure that the transport was at least partially diffusive and the phonon mean free path (MFP) was not limited just by the length of the flake. In the later measurements, we utilized well-defined massive metal heat sinks and an elaborate procedure for the thermal conductivity extraction based on the numerical solution of the heat diffusion equation. The single-layer graphene flakes were selected using micro Raman spectroscopy by checking the intensity ratio of G and 2D peaks and by 2D band deconvolution [10]- [13]. The combination of these two Raman techniques with atomic force microscopy (AFM) and scanning electron microscopy (SEM) allowed us to verify the number of atomic planes and flake uniformity with a high degree of accuracy. The 2D band deconvolution is currently considered to be the most robust, non-destructive and accurate technique for counting the number of atomic layers in graphene. One should note here that the thermal interface resistance between the flake and the heat sink is not such a big issue as in the case of carbon nanotubes (CNTs). The graphene flake is naturally attached to the graphitic piece and heat spreads over a large area.
The challenge in the measurement of the thermal conductivity with the described optical technique is in accurately determining the power absorbed in the graphene. Only a fraction P G of the laser light focused on the graphene flake will actually be dissipated in the graphene. Most of the light will be reflected back after the light travels through the flake to the trench bottom and is reflected back (see figures 1(a) and (b)). The power, which is measured by the detector placed at the position of the flake, is the total power P D , part of which goes into the graphene flake after two transmissions (incident pass and reflected pass) and the rest is lost in the silicon wafer P Si . It is now known that the fraction of the power absorbed by graphene is 2.3% per layer for light wavelength λ > 500 nm. Our measurements were performed at a smaller wavelength (λ = 488 nm) where the absorption is enhanced [18,19]. The stronger absorption at smaller wavelengths is explained by the effect of surface contaminations and bending of the suspended part of graphene. The near-field effects in the gap between the suspended graphene and trench bottom may also lead to additional absorption. Thus, it was important to determine the absorbed power in the specific conditions of our experiment. The power P G was measured through the calibration procedure with bulk graphite serving as a reference. The calibration and measurement of the absorbed power is illustrated in figure 2. It is based on comparison of the experimentally determined integrated Raman intensity for the G peak from singlelayer graphene and bulk graphite. In the next section, we give details of the derivation of the calibration formula for single-layer graphene. The measurements over the trenches with different depths excluded interference effects [12]. No electrical bias was used at the back gate or heat sinks to avoid the electrostatic effects on the G peak position. The RT position of the G peak was verified for all samples.

Details of the measurements and data extraction technique
In this section, we provide details of the derivation of the formulae, which allow us to determine the fraction of light power absorbed in graphene. The intensity of Raman scattering in graphene 6 is given by where N is the number of scattering atoms on the illuminated surface of the cross-section A, I 0 is the laser intensity and σ G is the Raman scattering cross-section. Let us assume that an electromagnetic (EM) wave with the intensity I = I 0 is incident on a graphene layer. The electric field of the wave transmitted through the graphene layer can be written as where n 0 is the refractive index, κ is the extinction index and a G is the graphene thickness. Using the notations A 1 = 2πκa G /λ and A 2 = 2π n 0 a G /λ, we can rewrite the above equation as The power absorbed in the graphene monolayer is given by the expression where the averaging is carried over the wave period T = 2π/ω and α G = 4πκ/λ. The EM wave incident on the surface of silicon under the suspended graphene is given as where z 0 is the distance from the graphene monolayer to the silicon surface (trench bottom), which gives an extra phase factor of exp[2π iz 0 /λ]. The reflected wave from silicon will have a different amplitude E r and additional phase ϕ given as By definition, the measurable macroscopic reflection coefficient is a ratio of the intensities of the incident and reflected waves, i.e. R Si = I 3 /I 2 . Here, we can write the following: Thus, the reflection coefficient can be expressed as: The amplitude of the reflected wave on the bottom side of the suspended graphene monolayer is given as After transmitting through the graphene layer for the second time, i.e. from the bottom up, the electrical field of the EM wave will be written as As a result, the power absorbed in the graphene monolayer is given by In the above equations, we took into account that The total power absorbed in graphene after the EM wave passes through it twice can now be written as Expressing I 0 through P G we have Plugging equation (10) into equation (1), we obtain the following equation for the integrated intensity of Raman light reflected from graphene: Let us now consider the light absorption and Raman scattering in bulk graphite. The highly oriented crystalline graphite consists of atomic planes of graphene, which are bound by weak Van der Waals interaction. For light with energy above 0.5 eV graphite behaves essentially as a collection of independent graphene layers [20]. This is also confirmed by the results of measurements of the light absorption in graphene multilayers, which prove that for a wide range of wavelengths the absorption is constant per atomic layer [18,19]. In our analysis, we take into account that the absorption coefficient and Raman scattering cross-section defined per layer are the same for graphene and for the single atomic plane in bulk graphite. The electric field intensity created by the incident EM wave in the nth atomic layer is , where a B is the distance between atomic planes in bulk graphite. After the Raman scattering event the intensity in the scattered wave will be given by where δ n (t) is the phase, which appears after absorption or emission of a phonon. This phase is a random function of time varying from phonon to phonon. Summing the electric field intensities of the Raman scattered waves one can obtain the following: × cos(4πn 0 n a B /λ − ωt + 2πδ n ). In order to average over the random phases we recall that the cosines' values are between −1 and +1 with the average for the random variables being equal to zero for n = n . Therefore, one can write the following expression: Averaging over the whole time period, we obtain for the phase factor where ϕ 0 = 4πn 0 na B /λ + 2π δ n .

8
As a result, for the intensity of the Raman scattered light from the stack of atomic planes, we can write where α B = 4π κ B /λ. Based on the above expression, the integrated Raman scattering intensity can be written as . Taking into account that part of the incident power is reflected from graphite and does not contribute to the Raman scattering intensity, one obtains where R B is the graphite reflection coefficient.
Combining equations (11) and (16) and introducing the ratio ς = I G / I Bulk , which has to be determined experimentally, we obtain the final expression for the power absorbed in graphene: The above expression allowed us to measure the heating power dissipated in the suspended portion of graphene under the specific conditions of the experiment. The ratio ς = I G / I Bulk , measured for graphene and reference bulk graphite for the G peak and the same frequency interval, stays nearly constant over the excitation power range used in the experiments. The reflection coefficients for silicon and bulk graphite are tabulated quantities but can also be measured directly while P D is the actual reading of the power detector taken at the position of the sample. Since the microscopic in-plane Raman cross-sections and absorption coefficients are same for graphene and bulk graphite, the square bracket term in equation (17) is close to one. With equation (17) serving as the calibration law for converting P D to P G , the measurement of the thermal conductivity of suspended graphene reduces to measuring the Raman shift ω G as a function of the heating power P D determined by the detector (see figure 3). The measured slope ω/ P D , ratio of the integrated intensities ζ and the temperature coefficient χ G give the value of the thermal conductivity of graphene. Initially, the thermal data extraction was accomplished using the simple 1-D model, and later improved by utilizing the numeric solution of the heat diffusion equation for a given flake shape. It has been found that the thermal conductivity for the single-layer partially suspended graphene is in the range of ≈ 3000-5300 W mK −1 depending on the width (lateral size) of the graphene flakes. This value is very high compared to other carbon materials [21]- [26]. Table 1 summarizes the thermal conductivities of graphene, CNTs and bulk graphite. One can see that there is a wide data scatter for the reported values of the thermal conductivity of CNTs. The conventionally accepted values for CNTs are K ≈ 3000-3500 W mK −1 . Thus, graphene can outperform CNTs as a heat conductor. Owing to its planar geometry, graphene may have potential for lateral heat spreading. It still remains unclear as to how the thermal conductivity of graphene will be affected when a graphene layer is embedded inside the device structure or how the heat conduction ability will evolve by increasing the number of atomic layers.

Formal theory of the phonon heat conduction in graphene
In this section, we outline the formal theory of heat conduction in graphene. The heat flux along a graphene atomic plane can be calculated according to the expression given in [27]. Multiplied by the system volume the 'total heat flux' can be written as where v(s, q)hω s ( q) is the energy carried by one phonon, v(s, q) = dω s /dq is the phonon group velocity and N (ω, q) = N 0 (ω, q) + n(ω, q) is the number of phonons in the flux. Here N 0 is the Bose-Einstein distribution function and n = −τ tot ( v∇T )∂ N 0 /∂ T is the non-equilibrium part of the phonon distribution function N , where τ tot is the total phonon relaxation time and T is the absolute temperature.
Comparing the microscopic expression with the macroscopic definition of the thermal conductivity we obtain the following expression for the thermal conductivity tensor: Here L x = d is the sample width (graphene flake width), L y is the sample length and h = a G = 0.35 nm is the thickness of graphene. The diagonal element of the thermal conductivity tensor, which corresponds to the phonon flux along the temperature gradient, is given by Finally, making a transition from the summation to the integration and taking into account the 2D density of phonon states, we obtain the following expression for the scalar thermal conductivity: The detailed description of the theoretical formalism for the phonon heat conduction in graphene was recently reported by some of us elsewhere [7]. Our theoretical approach utilized an original phase-diagram technique to account for all three-phonon Umklapp scattering channels allowed by the energy and momentum conservation. We consider two types of the three-phonon Umklapp scattering processes [27]- [29]. The first type is the scattering when a phonon with the wave vector q(ω) absorbs another phonon from the heat flux with the wave vector q (ω ), i.e. the phonon leaves the state q. For this type of scattering process, the momentum and energy conservation laws are written as where b i , i = 1, 2, 3, is one of the vectors of the reciprocal lattice (see figure 4(a)). The processes of the second type are those that occur when the phonons q of the heat flux decay into two phonons with the wave vectors q and q leaving the state q, or, alternatively, when two phonons q (ω ) and q (ω ) merge together forming a phonon with the wave vector q(ω), which corresponds to the phonon coming to the state q(ω). The conservation laws for this type are given by To find all the possible three-phonon processes, we used a fine mesh q j = ( j − 1) q ( j = 1, . . . , 1001) with the step q = q max /1000 ≈ 0.015 nm −1 . For each phonon mode (q i , s), we found all pairs of the phonon modes ( q , s ) and ( q , s ) such that the conditions of equations (24) and (25) are met. As a result, we constructed in ( q )-space the phase diagrams for all allowed three-phonon transitions [7]. A sample phase diagram is depicted in figure 4(b) for different (q i , s) and ( q , s ). Using the general expression for a matrix element of the three-phonon interaction [27]- [29] and taking into account all relevant phonon branches and their dispersion as well as all unit vectors of the reciprocal lattice b 1 . . . b 6 , directed from the point to the centers of the neighboring unit cells (see figure 4(a)), we obtain the following for the Umklapp scattering rates: 1 Here q l and q ⊥ are the components of the vector q parallel or perpendicular to the lines defined by equations (24) and (25), correspondingly, γ s ( q) is the mode-dependent Gruneisen parameter, which is determined for each phonon wave vector and polarization branch and ρ is the surface mass density. In equation (26), the upper signs correspond to the processes of the first type, whereas the lower signs correspond to those of the second type. The integrals for q l , q ⊥ are taken along and are perpendicular to the curve segments, correspondingly, where the conditions of equations (24) and (25) are met (see figure 4(b)). Integrating along q ⊥ , we obtain the line integral The integration in equation (27) is carried out along the curve segments l, as shown in figure 4(b). The integration in the expression for the scattering rate is a nontrivial task. Previously, in the calculation of the thermal conductivity of conventional materials, it was done 12 rather roughly with many simplifying assumptions. Most authors assumed that the phonon with the wave vector q is on the surface of the BZ so that for any q, the following relation is satisfied: The frequency ω(q ) in the equations similar to our equation is replaced with ω(q ) = ω T , where ω T is the maximum value of the transverse acoustic phonons. Next, for each q on the BZ surface, a circle with the radius r C was introduced within which (for a given q) the conditions of equations (24) and (25) are satisfied. In such an approach, the integration over q is reduced to multiplication by the circle area πr 2 C . In order to capture the specifics of heat conduction in a 2D system such as graphene we avoided these simplifying assumptions, and performed the integration in equation (27) accurately. The latter was made possible via the scattering phase diagram technique developed by us.
As the first step, we found the points in BZ for which the condition ω( q) + ω( q ) = ω( q ) is satisfied. To do this, we fixed q and numerically examined all pairs of q and q belonging to BZ, and selected such points in 2D reciprocal space ( q , ω s ( q )) and ( q , ω s ( q )) where the conditions of equations (24) and (25) are valid for a given q. Such points form the curve segments in the scattering diagrams (see figure 4(b)). These curves are nothing but the lines along which the integral in equation (27) is taken. This procedure was repeated for all phonon polarization branches and all values of q. As a result, we found numerically, for all possible q and for all polarization branches TA, LA, ZA, TO, LO, ZO, the lines in BZ along which the conditions of equations (24) and (25) are valid. The integral in equation (27) was taken numerically along these lines. The contributions from the separate lines (segments in the phase diagrams shown in figure 4(b)) were summed to obtain the total Umklapp scattering rate. The outlined procedure constitutes the essence of our diagram technique.
The combined scattering rate in both types of the three-phonon Umklapp processes for a phonon in the state (s, q) can be calculated as 1 τ U (s, q) = 1 τ I U (s, q) .
One should note here that for the small phonon wave vectors (long wavelength), q → 0, the Umklapp-limited phonon life-time τ U → ∞. For this reason, the calculation of the intrinsic thermal conductivity with only Umklapp scattering is not possible without an arbitrary truncation procedure. To avoid non-physical assumptions about the limits of the integration, we accurately include the phonon scattering on boundaries. In the case of graphene, the boundary scattering term corresponds to scattering from the rough edges of graphene flakes. No scattering happens from the top and bottom sides of graphene flake since it is only one atomic layer thick and the phonon flux is parallel to the graphene plane. We evaluate the rough edges' scattering using the Ziman equation [30]: Here p is the specularity parameter, which depends on the roughness of the graphene edges and L is the width of the graphene flake. The total phonon relaxation rate is given by It is important to emphasize that the thermal conductivity of a 2-D system such as graphene cannot be determined without the restriction on the phonon MFP in the long wavelength limit. The phonon scattering on edges restricts the MFP in the formal theory of thermal conductivity. In this sense the thermal conductivity limited by the Umklapp and boundary scattering can be considered as an intrinsic property of a graphene flake of a particular size. The extrinsic effects, which reduce the thermal conductivity, such as phonon scattering on defects, impurities and grain boundaries are not included into the consideration. The RT intrinsic thermal conductivity of graphene as a function of the specularity parameter is shown in figure 5. The specularity parameter is defined by the edge roughness: p = 0 corresponds to the perfectly diffusive phonon scattering from the edges while p = 1 corresponds to the perfectly specular scattering, which does not contribute to the thermal resistance. The results were calculated for different sizes (widths) of the graphene flakes. The experimental data point after Balandin et al [3,4] is also shown for comparison.

Simple model of the phonon thermal conductivity of graphene
The calculation of the thermal conductivity within the formal theory outlined in the previous section is an accurate but rather complicated procedure. In this section, we describe a simpler model of the thermal conductivity of graphene. This model originates from Klemens' approach to the thermal conductivity analysis of bulk graphite and graphene [25,26]. We substantially modified it by using a more general expression for the thermal conductivity, introducing two Gruneisen parameters γ s obtained independently for each of the heat conducting phonon polarization branches s, and keeping separate the velocities and cut-off frequencies for each phonon branch. These changes allowed us to reflect the specifics of the phonon dispersion in graphene. The effective parameters γ s are computed by averaging the phononmode-dependent γ s ( q) for all relevant phonons (here q is the phonon wave vector). The phonon branches, which carry heat, are longitudinal acoustic (LA) and transverse acoustic (TA). The out-of-plane transverse acoustic phonons (ZA) do not make contributions to heat conduction due to their low group velocity and high γ s ( q). There is a clear difference in the heat transport in basal planes of bulk graphite and in single-layer graphene [25,26]. In the former, the heat transport is approximately 2-D only till some low-bound cut-off frequency ω C .
Below ω C there appears strong coupling with the cross-plane phonon modes and heat starts to propagate in all directions, which reduces the contributions of these low-energy modes to heat transport along basal planes to negligible values. In bulk graphite, there is a physically reasonable reference point for the onset of the cross-plane coupling, which is the ZO phonon branch near ∼4 THz observed in the spectrum of bulk graphite. The presence of the ZO branch and corresponding ω C allows one to avoid the logarithmic divergence in the Umklapplimited thermal conductivity integral and calculate it without considering other scattering mechanisms.
The physics of heat conduction is principally different in graphene where the phonon transport is purely 2-D all the way to zero-phonon frequency ω(q = 0) = 0. There is no onset of the cross-plane heat transport at the long-wavelength limit in the system, which consists of only one atomic plane. There is no ZO branch in the phonon dispersion of graphene (see figure 6). Thus, the cut-off frequency for Umklapp processes cannot be introduced by analogy with bulk graphite. Using an expression for the three-phonon Umklapp scattering from [25,26] but introducing separate lifetimes for LA and TA phonons, we have where s = TA or LA, υ s is the average phonon velocity for a given branch, T is the absolute temperature, k B is the Boltzmann constant, ω s,max is the maximum cut-off frequency for a given branch and M is the mass of the unit cell. To determine γ s , we averaged γ s (q) obtained from the accurate phonon dispersion calculated by using the VFF method [7] and the ab initio theory [31]. Substituting τ tot = τ K U,s in equation (23), one can obtain the following formula for intrinsic thermal conductivity in graphene: The above equation can be used to calculate the thermal conductivity with the actual dependence of the phonon frequency ω s (q) and the phonon velocity dω s (q)/dq on the phonon wavenumber. To simplify the model, we can use the liner dispersion ω s (q) = υ s q and rewrite it as Substituting equation (31) into equation (33) and performing integration we obtain where F(ω s,min , ω s,max ) = hωs,max/kBT In the above equation, ξ =hω/k B T , and the upper cut-off frequencies ω s,max are defined from the actual phonon dispersion in graphene, calculated by using the VFF model [7] and presented in figure 6: ω LA,max = 241 THz, ω TA,max = 180 THz. The low-bound cut-off frequencies ω s,min for each s are determined from the condition that the phonon MFP cannot exceed the physical size L of the flake, given as The integrand in equation (35) can be further simplified near RT whenhω s,max > k B T , and it can be expressed as The obtained equations (34) and (37) constitute a simple analytical model for calculation of the thermal conductivity of the graphene layer, which retains such important features of graphene phonon spectra as different υ s and γ s for LA and TA branches. The model also reflects the 2D nature of heat transport in graphene all the way down to zero-phonon frequency. Our equation (31) reduces to Klemens' formula for graphene [26] in the limit ξ → 0 (hω k B T ) and additional assumption of the same γ s and υ s for LA and TA phonons. Using equations (32) and (34)-(37) of the simple model we calculated the Umklapp-limited intrinsic thermal conductivity of graphene as a function of temperature. The results are shown in figure 7 for two different widths of graphene flakes. The Gruneisen parameters used in this calculation, γ LA = 1.8 and γ TA = 0.75, were obtained by averaging of γ s (q) [32]. The thermal conductivity calculated with the exact ω s (q) and the phonon velocity dω s (q)/dq (shown by the dashed curves) is slightly smaller (by 5-10%) than those obtained with the simplified equations (34)-(37) due to the decrease of dω s (q)/dq with increasing q. The difference between thermal conductivities increases with increasing temperature as the conditionhω s,max > k B T breaks down. One should note a very different temperature dependence obtained within this model for the relatively small graphene flakes (10 µm) as compared to that for large flakes (100 µm). In the very narrow graphene flakes and nanoribbons, the thermal conductivity increases with temperature, which is related to the size (edge) effect on the phonon MFP.
In figure 8, we present the calculated RT thermal conductivity as a function of the flake width. The dashed curves are obtained using equation (23) with the total relaxation rate 1/τ tot (s, q) = 1/τ K U,s (ω s (q)) + 1/τ B (s, q), explicitly taking into account the phonon scattering on graphene edges. The solid curve is obtained from the simplified approach (equations (31)-(37))  [3,4] is also shown for comparison.  [3,4] is also shown for comparison.
in which the edge scattering is not directly included. Instead, the phonon MFP is limited by the physical size of the flake. The dashed curves are plotted for different values of the specularity parameter p. Using both approaches we obtained similar dependences: the intrinsic thermal conductivity of graphene grows with the increasing linear size of the graphene flake L. This is a manifestation of the two-dimensional nature of the phonon transport in graphene. In real experimental conductions, the thermal conductivity will also be limited by extrinsic factors (defects, impurities, grain size, etc), which prevent the growth of the thermal conductivity for very large flakes. With decreasing specularity parameter p (more diffuse scattering), the thermal conductivity calculated with the boundary scattering term approaches the result obtained with the simple model. This is because in the simple model we neglect phonons with frequency ω < ω min,s by completely restricting the phonon MFP to the lateral sizes of the flake. The latter corresponds to the perfectly diffusive scattering case ( p = 0).
The described models are in good agreement with the measurements and shed light on the heat conduction properties of graphene. They also help to understand the extremely large range of thermal conductivity values in carbon materials, which vary from some of the lowest in amorphous carbon, nanocrystalline diamond and diamond-like carbon [33]- [36] to the highest in graphene and CNTs (see table 1). The simple model can be readily incorporated into simulation software for the analysis of heat conduction in graphene layers and graphene devices [37,38]. The superior thermal conductivity of graphene is beneficial for all of its proposed device applications, e.g. in low-noise transistors, sensors and interconnects [39,40]. The application of graphene or few-layer graphene for lateral heat spreading would require special heat sinks and graphene's integration with suitable electrical insulators, e.g. synthetic diamonds. The thermal conductivity of graphene embedded within SiO 2 layers can be reduced due to the interface effects. It also decreases with temperature due to stronger Umklapp scattering at elevated temperatures [41]. A recent study of energy dissipation in graphene field-effect transistors [41] revealed the temperature dependence of the thermal conductivity of graphene in line with our theoretical prediction [7]. One should note however that our measurements were performed for partially suspended graphene with lateral heat sinks while the study reported in [41] is for graphene on SiO 2 with heat sinks at the substrate bottom, which led to a different heat front and its dissipation through the oxide layer.

Conclusions
We reviewed the results of our experimental investigation of heat conduction in suspended graphene and presented their detailed theoretical interpretation. We explain the enhanced thermal conductivity of graphene as compared to that of bulk graphite basal planes by the 2D nature of heat conduction in graphene over the whole range of phonon frequencies. The intrinsic Umklapp-limited thermal conductivity of graphene grows with the increasing linear dimensions of graphene flakes and can exceed that of bulk graphite when the flake size is on the order of a few micrometers. The superior thermal properties of graphene are beneficial for all proposed graphene device applications.