# Phonon Transport in Large Scale Carbon-Based Disordered Materials: Implementation of an Efficient Order-N and Real Space Kubo Methodology

## Abstract

We have developed an efficient order- real space Kubo approach for the calculation of the phonon conductivity which outperforms state-of-the-art alternative implementations based on the Green’s function formalism. The method treats efficiently the time-dependent propagation of phonon wave packets in real space, and this dynamics is related to the calculation of the thermal conductance. Without loss of generality, we validate the accuracy of the method by comparing the calculated phonon mean free paths in disordered carbon nanotubes (isotope impurities) with other approaches, and further illustrate its upscalability by exploring the thermal conductance features in large width edge-disordered graphene nanoribbons (up to nm), which is out of the reach of more conventional techniques. We show that edge-disorder is the most important scattering mechanism for phonons in graphene nanoribbons with realistic sizes and thermal conductance can be reduced by a factor of 10.

###### pacs:

72.80.Vp,72.15.Rn,73.22.Pr## I Introduction

In recent years, the understanding of phonon transport in carbon-based materials such as carbon nanotubes (CNTs) (1) and graphene-based materials (2) has become particularly important, both for fundamental studies of coherent transport and also in view of novel applications. The thermal conductivity of suspended and supported single graphene layers has been found to be extremely high (3); (4) owing to micrometer long phonon mean free paths (MFP). Such high thermal conductivity of two-dimensional graphene increases its potential for faster nano-electronic devices with less energy dissipation (5). On the other hand, damaging a material like graphene could present interesting opportunities like achieving a high thermoelectric figure of merit and observation of Anderson localization of phonons. The question whether or not Anderson localization of acoustic phonons can be demonstrated unambiguously in disordered materials has been a long-standing problem in phonon physics, a phenomenon originating from the interference between multiple scattering paths was found ubiquitous in wave physics (6); (7). Besides, disordered CNT-based bundles have been found to exhibit a tendency towards a thermal insulating regime (8); (9). Also, isotope disorder was shown to strongly impact on the high-energy phonon modes, resulting in very low mean free paths, in marked contrast with the genuine robustness of ballistic conduction for low-energy phonon modes (10). As a result, the contribution of quantum localization effects of high energy modes was found to be important but completely masked when considering the thermal conductance of the disordered material (being an integrated quantity over the entire phonon spectrum). Graphene nanoribbons (GNRs) offer an alternative to carbon nanotubes. By constructing heterostructures from pristine graphene mixed with disordered (isotope impurities) GNRs (11), or selective coverage of hydrogen impurities (12), it is possible to tune the resulting phonon transport capability. Another important source of disorder is provided by unavoidable ribbon edge irregularities, absent in CNTs which have already been shown to significantly reduce thermal conductance in small width GNRs (13); (14).

Exploring such possibilities, however, demand the development of efficient computational approaches which are able to tackle large scale (and realistic) simulations of material of interest. The Green’s function (GF) methods are able to incorporate microscopic details of the system, but they require matrix inversions which unavoidably limit the accessible size of the simulated systems. For the case of GNRs, although the system length can be upscaled without difficulty thanks to the decimation procedure, the computational cost becomes prohibitive for lateral sizes above 10 nm. (10); (13)

In this Rapid Communication, using the real space Kubo (RSK) formalism we first demonstrate that a time-dependent phonon wave packet formalism can be connected to the calculation of the thermal conductance. After validating this numerical approach by comparing the obtained phonon MFPs in disordered CNTs (with isotope impurities) with previously computed ones by means of GF-based method (10), we apply the new algorithm to large width GNRs, and focus on the impact of edge-disorder profiles. Scaling properties of phonon MFP and temperature-dependent thermal conductance are calculated as a function of edge-disorder strength and for lateral ribbon sizes accessible to today’s state-of-the-art lithography. The broad generality of this method could offer a novel framework to explore other types of complex materials.

## Ii Computational phonon transport methodology

To investigate bulk quantum phonon transport in disordered materials, the use of the Kubo formalism turns out to be the most natural and computationally efficient one. It has already been used for investigating thermal transport in disordered binary alloys or nanocrystralline silicon (15); (16). Similarly, many efforts have been made to explore time-dependent features of propagating phonon (or polaron) wave packets in complex quantum systems (17). Here a novel real space implementation of the Kubo formula for phonon propagation is given, establishing a direct computational bridge between phonon dynamics and the thermal conductance. In comparison to the previous methods, we extract the dynamical information from time evolution of the wave packet (16) and simulate Hamiltonian dynamics in place of Newtonian equations of motion (18), so that only one initial condition is required, i.e. the initial atomic displacements, and one does not need to calculate the atomic velocities in time. Additionally, by using the Lanczos technique which avoids any matrix inversion, a considerable computational efficiency gain is obtained, allowing the study of very large scale materials. The starting Hamiltonian which describes the phonon spectrum, taking only the harmonic interactions into account, is of the form

(1) |

where and are the displacement and momentum operators for the th degree of freedom, is the corresponding mass, and is the force constant tensor. Based on the linear response theory, the phonon conductivity can be defined as (16). is the component of the energy flux operator , and it can be expressed as , where is the velocity operator and is the equilibrium position of the atom to which the th degree of freedom belongs. After some algebra, becomes

(2) |

where is the Bose distribution function, and is the mass normalized dynamical matrix, denotes the trace of the matrix, and is the Dirac-delta operator. . The dynamical matrix can be obtained given the atomistic configuration of the system and it defines a Newtonian equation of motion which is of second order in time. Using the fact that the operators and have the same energy spectrum () and defining , one can write the thermal conductance of a one-dimensional system as

(3) |

The thermal conductance can also be derived from the Landauer formalism (19) or the nonequilibrium Green’s function approach (20) as

(4) |

with being the phonon transmission function. Comparing the two formulas, the transmission function is defined as

(5) |

which has the same form as the electron transmission function derived from the Kubo-Greenwood formula (21),

(6) |

This equivalence allows us to implement an order algorithm related to the Lanczos method, which has been successfully applied to electron conduction in complex materials (21). The Lanczos method has been previously employed in the past for calculating the vibrational density of states of disordered systems (22). The starting point of the RSK method is the transmission function as expressed in Eq. (6). One can rewrite in terms of the diffusion coefficient,

(7) |

where the diffusion coefficient has the form

(8) |

where is the position operator in the Heisenberg picture. In the diffusive regime , so that the transmission function reduces to

(9) |

where in the ballistic regime , with being the average group velocity over states with frequency . Since the number of channels is , the phonon mean free path can be approximated by . By computing , one can thus deduce and and . Note that the numerator in Eq. (8) can be rewritten as , and approximated by to the first order in perturbation theory, with , and . The trace can be efficiently calculated through an average over a few random phase states of atomic displacements as , being the number of degrees of freedom and the bra-ket corresponds to the local density of states (LDOS) associated with the vector which is calculated by using Chebyshev expansion of . LDOS is obtained by using the Lanczos method, and is also calculated by averaging the LDOS of through the Lanczos method.

## Iii Results

As a test case, we first consider a CNT(7,0) with 10.7% C impurities, which was studied in Ref. (10) using the GF method. From the saturation values of the time-dependent diffusion coefficients, we obtain the phonon MFP using the RSK scheme which compares very well with those obtained from the GF method (23).

Next, we consider zigzag graphene nanoribbons (ZGNR) of different widths with edge-disorder. We use the fourth nearest neighbor force constants for building the dynamical matrices (24). The ribbon widths are defined with the number of zigzag chains , 40, and 80, and the relative amount of edge defects (removed carbon atoms at the edges) is chosen to be 10%, and additionally 15% for (Fig. 1). The inset of Fig. 2 displays the evolution of the wave packet dynamics for different frequencies for ZGNR(80) with 10% edge-disorder.. The linear increase in at indicates ballistic transport at relatively short distances, whereas the decrease in is a signature of localization for this particular frequency. The fact that decreases slowly suggests the localization length is large. The saturation of to a maximum value characterizes diffusive transport (c.f. Eq. 9). In Fig. 2, the decay of with decreasing the ribbon width is shown at a fixed disorder strength. This behavior can be rationalized with the fact that the scattering rate decreases with increasing width, a behavior previously derived for electron transport in both disordered CNTs and GNRs (25). One notes that, for low frequency modes, the MFP are several hundreds of nanometers, and due to large values of the possibility to observe any onset of Anderson localization is jeopardized in the thermal conductance, as previously discussed for small diameter disordered carbon nanotubes (10). We then focus on ZGNR(80), and obtain the transmission according to , assuming a diffusive regime, disregarding quantum interference effects and assuming minimum contact resistances. The resulting frequency-dependent transmission function for different ribbon lengths is plotted in Fig. 3, while the thermal conductance using Eq. (4) is shown in Fig. 4. The downscaling of directly impacts on which is found to be reduced by one order of magnitude for 2 m (at room temperature) compared with the ballistic case. Note that the calculation time directly determines the largest that we can reach, the same corresponds to different evolution time for different frequencies, since . It takes longer calculation time for low frequency phonons to reach . Here for , has not been reached within the finite calculation time. To compute the contribution of these modes to , a linear extrapolation for the transmission is used and the results are compared with those obtained from the GF method. Our analysis shows that this approximation causes an error of less than 1.5% for the thermal conductance at room temperature. The conductivity of edge-disordered GNRs is obtained using , being the cross section area of the ribbons which is taken as the interplane distance of graphite layers. At room temperature, WmK and WmK with edge-disorder 10% and 15%, respectively for 2 m. Comparing these values with the experimental values for suspended(3) (3000-5000 WmK) and supported(4) (600 WmK) graphene, we conclude that edge-disorder is a very important source of scatterings not only in ultra-narrow ribbons, but also in GNRs as wide as at least 20 nm. A comparison of our results with that of Ref. (13) suggests that the effect of edge-disorder is reduced when increasing the width from 2 nm to 20 nm although it still plays an important role for large lateral sizes. The ratio of edge atoms affected from edge reconstructions to the total number of atoms decay inversely with the width of the GNR (14), therefore edge profile disorder predominates over reconstruction effects with increased ribbon width.

## Iv Conclusion

Conclusion.-A new real space and order approach to compute the phonon wave packet propagation and the thermal conductance has been reported. Diffusion coefficients and MFPs can be extracted directly from the wave packet propagation. Its computational accuracy and efficiency were demonstrated on disordered carbon nanotubes and large width graphene nanoribbons respectively. A strong impact of smooth edge-disorder on the thermal conductance was found. Unlike edge reconstruction (14), edge-disorder strongly suppresses thermal conduction not only for ultra-narrow GNRs, but also for realistically large ribbons. Phonons in edge-disordered GNRs, being scattering so effectively, pinpoints towards good thermoelectric properties of large width GNRs. The applicability of this new method goes far beyond quasi-one dimensional systems and the area of carbon-based materials (nanotubes, graphene, etc.) studied here, and could be applied without difficulty to a wide range of other materials, including Boron-nitride-based materials (26) or silicon-based materials (nanowires, superlattices, etc.) (27).

Acknowledgements.- This work was supported by the priority program Nanostructured Thermoelectrics (SPP-1386) of the German Research Foundation (DFG) (Contract No. CU 44/11-1), the cluster of excellence of the Free State of Saxony ECEMP - European Center for Emerging Materials and Processes Dresden (Project A2), and the European Social Funds (ESF) in Saxony (research group InnovaSens), the ANR/P3N2009 GRAPHENE_NANOSIM (Project No. ANR-09-NANO-016-01), and the Alexander von Humboldt Foundation. We acknowledge support from the WCU (World Class University) program sponsored by the South Korean Ministry of Education, Science, and Technology Program, Project no. R31-2008-000-10100-0. The authors are thankful to N. Mingo for fruitful discussions and W. L. thanks CAS-MPG joint doctoral program. The Center for Information Services and High Performance Computing (ZIH) at the TU-Dresden is acknowledged.

### References

- J. C. Charlier, X. Blase, S. Roche, Rev. Mod. Phys. 79, 677 (2007).
- A. Cresti et al., Nano Res. 1, 361 (2008).
- A. A. Balandin et al., Nano Lett. 8, 902 (2008).
- J. H. Seol et al., Science 328, 213 (2010).
- E. Pop, Nano Res. 3, 147 (2010); M.-H. Bae, Z.-Y. Ong, D. Estrada and E. Pop, arXiv:1004.0287v2, Nano Lett. (to be published).
- P. W. Anderson, Phys. Rev. 109, 1492 (1958).
- A. Lagendijk, B. van Tiggelen, D. Wiersma, Phys. Today 62 (8), 24 (2009).
- J. Che, T. Çağın, W. A. Goddard, Nanotechnology 11, 65 (2000).
- R. S. Prasher et al., Phys. Rev. Lett. 102, 105901 (2009).
- I. Savic, N. Mingo and D. A. Stewart, Phys. Rev. Lett 101, 165502 (2008)
- N. Mingo et al., Phys. Rev. B 81, 045408 (2010).
- X. Ni et al., App. Phys. Lett. 95, 192114 (2009).
- H. Sevinçli and G. Cuniberti, Phys. Rev. B 81, 113401 (2010).
- J. Lan, J.-S. Wang, C. K. Gan and S. K. Chin, Phys. Rev. B 79, 115401 (2009).
- A. Alam and A. Mookerjee, Phys. Rev. B 72, 214207 (2005); A. Bodapati et al., Phys. Rev. B 74, 245207 (2006).
- P. B. Allen and J. L. Feldman, Phys. Rev. Lett. 62, 645 (1989); P. B. Allen and J. L. Feldman, Phys. Rev. B 48, 12581 (1993).
- G. G. Naumis, F. Salazar, C. Wang, Phil. Mag. 86, 1043 (2006); N. Kondo, T. Yamamoto, K. Watanabe, Jap. J. Appl. Phys. 45, L963 (2006); S. Monturet and N. Lorente, Phys. Rev. B 78, 035445 (2008).
- Y. L. Loh, S. N. Taraskin, and S. R. Elliot, Phys. Rev. E 63, 056706 (2001).
- L. G. C. Rego and G. Kirczenow, Phys. Rev. Lett 81, 232 (1998); N. Mingo, Phys. Rev. B 74, 125402 (2006); N. Mingo and L. Yang, Phys. Rev. B 68, 245406 (2003).
- A. Özpineci and S. Ciraci, Phys. Rev. B 63, 125415 (2001); T. Yamamoto and K. Watanabe, Phys. Rev. Lett 96, 255503 (2006).
- F. Triozon et al., Phys. Rev. B 69, 121410(R) (2004); A. Lherbier et al., Phys. Rev. Lett. 100, 036803 (2008).
- P. M. Lam, W. Bao, Z. Zheng, Z. Phys. B 59, 63 (1985).
- H. Sevinçli, W. Li, S. Roche and G. Cuniberti, to be published elsewhere.
- R. Saito, G. Dresselhaus, M. S. Dresselhaus, Physical Properties of Carbon Nanotubes, Imperial College Press, London (1998); J. Zimmermann, P. Pavone, G. Cuniberti, Phys. Rev. B 78, 045410 (2008).
- D. A. Areshkin, D. Gunlycke, C.T. White, Nano Lett. 7, 204 (2007);
- I. Savic, D. A. Stewart, N. Mingo, Phys. Rev. B 78, 235434 (2008); D. A. Stewart, I. Savic, N. Mingo, Nano Lett. 9, 81 (2009); L. Ci et al., Nature Mat. 9, 430 (2010).
- T. Markussen et al., Phys. Rev. Lett. 99, 076803 (2007); T. Markussen, A.-P. Jauho, M. Brandbyge, Phys. Rev. Lett. 103, 055502 (2009); D. Donadio and G. Galli, Phys. Rev. Lett. 102, 195901 (2009).