Is Thermal Conductivity of Graphene Divergent and Higher Than Diamond?

TL;DR

This study uses advanced first-principles calculations to determine that monolayer graphene's thermal conductivity converges at about 1300 W/(m·K) when four-phonon scattering is considered, and is lower than diamond's.

Contribution

It provides the first comprehensive first-principles analysis including four-phonon scattering, resolving the debate on graphene's thermal conductivity and its size dependence.

Findings

Thermal conductivity converges at 1300 W/(m·K) with four-phonon scattering.

Considering only three-phonon scattering leads to divergence.

Graphene's thermal conductivity is lower than diamond's.

Abstract

The thermal conductivity of monolayer graphene is an outstanding challenge with no consensus reached on its exact value and length convergence so far. We consider four-phonon scattering, phonon renormalization, and an exact solution to phonon Boltzmann transport equation (BTE) from first principles. Using this computational formalism with unprecedented sampling grid, we show that when four-phonon scattering is included the thermal conductivity is convergent with system size at a room temperature value of 1300 W/(m$\cdot$K), which is lower than that of diamond. On the contrary, considering three-phonon scattering only yields divergence with size due to the momentum-conserving normal processes of flexural phonons.

Full text

††preprint: APS/123-QED

Graphene is a subject of extensive research due to its unique properties and exotic transport phenomena observed [1, 2, 3, 4, 5, 6]. Among them, thermal conductivity ($\kappa$) of graphene, which is dominated by phonon transport, is one of the open questions that remains widely debated: its reported divergence with system size is controversial and its value higher than diamond and graphite is widely believed. Despite numerous experimental studies, the measured values at room temperature are scattered [7, 8, 9, 10, 11]. The difficulty in experiments is exemplified in the first measurement of graphene thermal conductivity using Raman thermometry [7]. Later, strong nonequilibrium phonon transport was revealed which causes additional challenges in extracting the thermal conductivity [12, 13]. In parallel, theoretical works also report a large spread of $\kappa$ between 800 to 3500 W/(m$\cdot$K) in atomistic simulations [14, 15, 16, 17, 18, 19].

To resolve the thermal conductivity of graphene, one must investigate two important yet elusive questions. One, the out-of-plane lattice vibration in graphene represented as the flexural phonons (ZA) [20] contributes significantly to $\kappa$ and its scattering was thought to be greatly suppressed by a selection rule for three-phonon (3ph) scattering [14]. But such selection rule does not apply to four-phonon (4ph) scattering [19], a higher-order mechanism found to be important in other solids [21, 22, 23, 24]. Two, except for the out-of-plane vibrations, graphene resembles the two-dimensional (2D) nonlinear lattice that was originated from the seminal Fermi-Pasta-Ulam-Tsingou (FPUT) model [25]. Heat transport in such perfect low-dimensional systems is non-Fourier [26] and $\kappa$ should diverge with system size as $\kappa_{\rm 2D}\propto\log(L)$ [27, 28]. It is a matter of debate whether real quasi-2D systems like graphene can recover diffusive regime when $L\to\infty$ and acquires a finite intrinsic thermal conductivity.

Most theoretical studies do concern the above two questions but consensus has not been reached yet. Lindsay et al. employed the linearized phonon Boltzmann transport equation (BTE) and coupled it with empirical potential [14] or first principles [16] at 3ph level. They proved that the exact solution to BTE is necessary for graphene due to its strong momentum-conserving scattering (normal process) and the length-dependent thermal conductivity didn’t show convergence up to 50 $\mu$m. Their finite-size results at room temperature ($\sim$ 3000 W/(m$\cdot$K)) are in accordance with another theoretical study at 3ph level [17] reporting length convergence at around 1 mm, but in this study [17] the sampling grid of the Brillouin Zone was not checked for $\kappa$ convergence. Along this BTE approach, Feng and Ruan [19] computed 4ph scattering rates with empirical potential and showed strong 4ph effect in ZA mode and consequently, a great reduction in $\kappa$ to around 800 W/(m$\cdot$K). The sampling grid accessible then was relatively coarse. Gu et al. further incorporated phonon renormalization into this calculation with empirical potential and argued that ZA scattering rates are reduced by temperature modification [29]. However, the empirical interatomic potentials used in these studies are not accurate representations of experiments and the exact role of four-phonon scattering and the value of $\kappa$ remains elusive. Another mainstream approach is molecular dynamics (MD) with two different flavors in equilibrium (EMD) and nonequilibrium (NEMD) treatments. Regardless of classical statistics and empirical potential, MD approach is expected to capture all orders of anharmonicity. However, the EMD simulations show finite $\kappa$ around 2000 W/(m$\cdot$K) [15, 18] while the NEMD simulations show logarithm length dependence up to several microns [11, 30]. Interestingly, one MD study presents a saturation of $\kappa$ when $L$ is extended to 100 $\mu$m [31].

In this Letter, we revisit the thermal conductivity of graphene with first-principles-computed four-phonon scattering rates and the exact solution to BTE after careful check on convergence. In addition, to describe the temperature-dependent potential field we consider the phonon renormalization effect in both harmonic phonon dispersion [33, 34] and anharmonic lattice dynamics [35, 36]. The considerations of four-phonon scattering and phonon renormalization have successfully explained the Raman linewidth of suspended graphene in our prior work [37]. We note that the prediction in Ref. [37] is supported by experimental measurements [38, 12] which validate our approach to compute Raman-active phonon scattering rates, and in this work we extend the methodology to all phonon modes. The phonon BTE is exactly solved by an iterative scheme [39] incorporating both three-phonon and four-phonon scattering after we manage to store the iterative processes within one terabyte memory space accessible in modern supercomputer architecture. We compute from first principles that ZA mode has strong four-phonon scattering rates which are comparable to three-phonon counterpart. The sampling grid ($q-$mesh) is carefully checked for $\kappa$ convergence of infinitely long graphene sample and we see convergence for thermal conductivity including four-phonon scattering ($\kappa_{\rm 3ph+4ph}$) but not in thermal conductivity including only three-phonon scattering ($\kappa_{\rm 3ph}$). Finally, the temperature-dependent $\kappa$ is compared to diamond and found to be lower than diamond from 300 K to 800 K. At room temperature, we predict that graphene has an intrinsic thermal conductivity about 1300 W/(m$\cdot$K) and is finite. Our findings through rigorous first-principles calculations have implications in the thermal transport of quasi-low-dimensional systems, question the perception of graphene being a better heat conductor than diamond, and will motivate future experimental efforts.

We consider naturally occurring, monolayer graphene in our simulation with a vacuum space of 14 Å between periodic graphene layers. The first-principles calculations are based on density functional theory as implemented in the VASP package [40]. The phonon renormalization effect is included by a temperature-dependent effective potential method (TDEP) [34] to compute temperature-dependent phonon dispersions and interatomic force constants (IFCs) [36]. Phonon scattering rates summed up by Matthiessen’s rule [41] $\tau_{\lambda}^{-1}=\tau_{\lambda,\rm{3ph}}^{-1}+\tau_{\lambda,\rm{4ph}}^{-1}+\tau_{\lambda,\rm iso}^{-1}$ are then computed by our FourPhonon code [42], which is an extension module to ShengBTE package [43]. The exact solution of BTE is implemented in the same solver by a shared-memory parallel computing strategy [44].

The computed phonon dispersion and phonon scattering rates are presented in Fig. 1(a) and 1(b), respectively. Overall, graphene is quite rigid as shown by very small change of phonon frequency at finite temperatures (see Fig. 1(a)). The frequency of in-plane optical phonons decreases with increasing temperature, a signature that has been analyzed and experimentally verified in our prior work [37]. Contrary to the decreasing trend of frequency shift for in-plane optical phonons, we find that the flexural phonons are hardened with rising temperature (see Fig. 1(a) right panel). This is a result of the coupling between flexural phonons and in-plane degrees of freedom in free-standing graphene [20]. Renormalized flexural phonons do not have a strictly quadratic dispersion. With the renormalized phonon dispersion and temperature-dependent IFCs, one can then compute the spectral phonon scattering rates $\tau_{\lambda}^{-1}$ using Fermi’s golden rule, as shown in Fig. 1(b). Only the optical phonons and flexural phonons have comparable four-phonon scattering rates ($\tau_{\rm{4ph}}^{-1}$) with three-phonon counterpart ($\tau_{\rm{3ph}}^{-1}$), while the rest of phonon modes are dominated by three-phonon scattering. Optical phonons having large $\tau_{\rm{4ph}}^{-1}$ are understood as they easily satisfy the energy conservation in the recombination scattering events: $\lambda_{1}+\lambda_{2}\to\lambda_{3}+\lambda_{4}$ [45]. On a different ground, the reason for the flexural phonons lies in the selection rule for general quasi-2D systems [14] that only even numbers of flexural phonons can be involved in a phonon scattering event due to reflection symmetry. As a result, the phase space for three-phonon scattering is much smaller than that of four-phonon scattering [19]. Although the ZA mode was considered as the major heat carrier in graphene under the 3ph scattering picture [14, 16, 46], this statement may no longer be true after its strong $\tau_{\rm{4ph}}^{-1}$ is included, which would greatly reduce the predicted thermal conductivity.

We next solve for the intrinsic thermal conductivity by an iterative scheme [39] to account for the collective phonon excitations. This approach distinguishes the normal processes (N) that are momentum-conserving and Umklapp processes (U) that are resistive. Such a treatment is important for graphene since it has strong normal scattering process associated with its phonon hydrodynamics nature [47, 48]. The coupled equations are presented in the Supplymental Material [49]. For the BTE solution to generate meaningful results, one has to check the convergence of $\kappa$ with respect to the sampling grid in the Brillouin Zone ($q-$mesh). A convergence being reached means that an infinitely large system has finite thermal conductivity. Figure 2 shows our test of convergence when 3ph or 3ph+4ph is included in the calculations. Mesh size $N\times N\times 1$ is uniform in-plane. In all these calculations, the energy broadening factor is unity to ensure the accuracy within our computational power. We find that with four-phonon scattering, $\kappa_{\rm 3ph+4ph}$ converges around a mesh size of $N=40$ and further increasing the mesh to $N=52$ does not change the value of $\kappa$. In contrast, $\kappa_{\rm 3ph}$ does not reach convergence up to $N=240$ and scales logarithmly with mesh size $N$. Note that the prior work [17] used $N=128$ for 3ph. Our numerical results share similarities with a theoretical work on carbon nanotube (CNT) [50] where the authors estimated that third-order anharmonicty leads to a divergence of thermal conductivity due to vanishing scattering of long wavelength phonons. They empirically showed that higher order process can remove such divergence. While the study on CNT is worth further investigations [51, 52, 53], we show here that four-phonon scattering converges the thermal conductivity of graphene from first principles. To inspect the origin of this behavior, we decompose the N/U scattering for ZA mode and plot the ratio of two types of scattering events $\tau^{-1}_{N}/\tau^{-1}_{U}$ in the inset of Fig. 2. Apart from the dominant role of normal process for both scattering channels, this plot indicates that for 3ph the N/U ratio keeps rising for low-frequency flexural phonons. Given the fact that $\tau^{-1}_{N}$ decreases when $\omega\to 0$, we conclude that 3ph of long wavelength flexural phonons is almost entirely contributed by normal processes. The extreme case is the heat conduction in 2D nonlinear lattice where all scattering processes are momentum-conserving [28] and the thermal conductivity is logarithmic divergent. Thus, we explain the finite intrinsic thermal conductivity of graphene from two related arguments: one, four-phonon scattering provides an additional scattering channel; two, for long wavelength phonons, four-phonon scattering has considerable number of resistive U scattering events but three-phonon processes are nearly all N scattering events implying a divergence.

The above arguments can further be seen in the spectral $\kappa$ as shown in Fig. 3. The 3ph and 3ph+4ph pictures have distinct spectral trend for low-frequency phonons. While the 3ph+4ph calculation can saturate the spectral $\kappa$ when $\omega\to 0$, the 3ph only calculation shows a up-soaring spectral trend. The integration of the shaded area in Fig. 3 should give the total $\kappa$ and this spectral analysis indicates that the 3ph case has a singularity at near-zero frequency.

Finally, we present the calculated thermal conductivity of graphene as a function of temperature and some comparisons to other carbon allotropes in Fig. 4. Since our 3ph calculations are not converged within the accessible mesh size, we consider a $10~{}\mu$m diameter graphene sample and add a boundary scattering term in our simulation [16, 17]: $\tau^{-1}_{\lambda,b}=\left|\mathbf{v}_{\lambda}\right|/L$, where $\mathbf{v}_{\lambda}$ is the group velocity of phonon mode $\lambda$. Our results on finite-size $\kappa_{L}$ at 3ph level are consistent with prior first-principles studies [16, 17]. The comparison in Fig. 4(a) shows that incorporating 4ph in the calculation reduces $\kappa_{\rm 3ph}$ nearly by half. Nevertheless, our results of $\kappa_{\rm 3ph+4ph}$ are higher than previous simulation with empirical potential at moderate mesh size [19]. At room temperature, we predict that $\kappa_{\rm 3ph+4ph}=1298$ W/(m$\cdot$K), and the thermal conductivity of graphene is lower than that of diamond from 300 K to 800 K. The first-principles results for diamond (black lines in Fig. 4(a)) are from Ref. [54] where similar methodologies were applied including 3ph/4ph scattering, phonon renormalization and iterative solution to BTE. The dashed black line is their result at 3ph level while solid black line at 3ph+4ph level. Note that the thermal conductivity of diamond is well received, with good agreements between simulations and experiments [55, 57, 60, 56]. Contrary to the situation of diamond, no consensus has been reached on the measured thermal conductivity value of graphene so far. While our finding challenges a popular perception that graphene is a better heat conductor than diamond, the numerical results are self-consistent with advancements in theories and computational power. Under current computational formalism, a plausible reason for this observation is that diamond does not have strong four-phonon scattering for acoustic modes [22, 54] but flexural phonons in graphene do have large four-phonon scattering rates originated from its 2D nature. In this sense, the reduction of $\kappa$ from diamond to graphene is understandable. Another consequence of our results is that the predicted graphene thermal conductivity is lower than experimentally reported $\kappa$ of graphite [58, 59] (Fig. 4(b)). Theoretical study on graphite considering both four-phonon scattering and phonon renormalization is not seen in literature yet. Future studies could be focusing on the layer-dependent transport behavior of multilayer graphene and graphite. In Fig. 4(b), we also cautiously compare our results to a recent Raman measurement on suspended graphene sample [12] that is based on apparent phonon temperature. We expect such comparison may motivate further experimental efforts to fully resolve the intriguing finding here. A tentative explanation might be that in bulk phase graphite the four-phonon scattering of ZA mode is suppressed by interlayer interactions.

In summary, we conduct a first-principles study on the thermal conductivity of monolayer graphene. Our calculations include four-phonon scattering, phonon renormalization effect and the exact solution to BTE, all of which are state-of-the-art computational formalism. Our methods can reveal both phonon properties, their detailed scatterings and eventually the thermal conductivity of graphene over a wide range of temperature. Part of the computed results have been supported by recent experiments. Our results provide a strong computational evidence of $\kappa$ convergence to date. We expect that our study may inspire further experimental efforts on graphene and theoretical understanding on general low-dimensional systems.