Science, Faculty of
Physics and Astronomy, Department of
DSpace
UBCV
Bruns, Daniel K.W.
Doctor of Philosophy - PhD
University of British Columbia
Heat conduction phenomena of carbon nanotubes (CNTs) have attracted great interest both from the viewpoint of engineering applications and fundamental science, but the thermal conductivity of individual single-walled CNTs remains a rather controversial topic. Starting from an empirical, realistic atomic interaction potential, we study the lattice thermal conductivity (LTC) of single-walled CNTs by employing two approaches: quantum mechanical calculations of three-phonon scattering rates in the framework of the Peierls-Boltzmann transport theory (PBTT) and classical molecular dynamics (MD) simulations. First, we compare the system-size and temperature dependence of the LTC determined from an iterative solution of the linearized phonon transport equation in the framework of the PBTT and from a nonequilibrium MD (NEMD) approach. At room temperature, qualitatively similar trends for the tube-length dependence are found in the limit of short tubes, where an extensive regime of ballistic heat transport prevailing in CNTs of lengths L ≲ 1 µm is independently confirmed. In the limit of long tubes, the PBTT-derived LTC diverges. Using PBTT and equilibrium MD approaches, we perform numerical calculations of acoustic phonon lifetimes to clarify the source of divergence. NEMD-derived temperature dependencies obtained for micrometer-long CNTs and temperatures T ≤ 800 K confirm the 1/T behavior of the LTC at moderately high temperatures. Next, we revisit the tube-length dependence of the LTC by use of the relaxation time approximation in the PBTT. Through a combination of numerics and analytical considerations, we derive exact asymptotic scaling laws of the LTC. In particular, we demonstrate the importance of tensile lattice strain, previously overlooked in the long-standing dispute over tube-length convergence vs divergence of the LTC. Namely, it is proved that, in the long-tube limit, the relaxation time approximation yields a finite value for stress-free but an infinite value of the LTC in any stretched tube configuration. Lastly, we pursue a matrix inversion approach to solve the linearized phonon transport equation in the framework of the PBTT. Here, it is shown that violations of acoustic sum rules cause spurious convergent behaviors of a length-dependent LTC.
Theory and molecular simulation of lattice vibrational heattransport in carbon nanotubesbyDaniel K.W. BrunsB.Sc., University of Rostock, 2014M.Sc., University of Rostock, 2016a thesis submitted in partial fulfillmentof the requirements for the degree ofDoctor of Philosophyinthe faculty of graduate and postdoctoral studies(Physics)The University of British Columbia(Vancouver)July 2022© Daniel K.W. Bruns, 2022The following individuals certify that they have read, and recommend to the Faculty ofGraduate and Postdoctoral Studies for acceptance, the thesis entitled:Theory and molecular simulation of lattice vibrational heat transport incarbon nanotubessubmitted by Daniel K.W. Bruns in partial fulfillment of the requirements for the degreeof Doctor of Philosophy in Physics.Examining Committee:Jo¨rg Rottler, Professor, Department of Physics & Astronomy, UBCCo-supervisorAlireza Nojeh, Professor, Department of Electrical & Computer Engineering, UBCCo-supervisorGeorge Sawatzky, Professor Emeritus, Department of Physics & Astronomy, UBCSupervisory Committee MemberMona Berciu, Professor, Department of Physics & Astronomy, UBCUniversity ExaminerChristoph Ortner, Professor, Department of Mathematics, UBCUniversity ExaminerPhilip B. Allen, Professor Emeritus, Department of Physics and Astronomy, Stony BrookUniversityExternal ExaminerAdditional Supervisory Committee Members:A. Srikantha Phani, Professor, Department of Mechanical Engineering, UBCSupervisory Committee MemberiiAbstractHeat conduction phenomena of carbon nanotubes (CNTs) have attracted great interest bothfrom the viewpoint of engineering applications and fundamental science, but the thermalconductivity of individual single-walled CNTs remains a rather controversial topic. Startingfrom an empirical, realistic atomic interaction potential, we study the lattice thermal con-ductivity (LTC) of single-walled CNTs by employing two approaches: quantum mechanicalcalculations of three-phonon scattering rates in the framework of the Peierls-Boltzmanntransport theory (PBTT) and classical molecular dynamics (MD) simulations.First, we compare the system-size and temperature dependence of the LTC determinedfrom an iterative solution of the linearized phonon transport equation in the frameworkof the PBTT and from a nonequilibrium MD (NEMD) approach. At room temperature,qualitatively similar trends for the tube-length dependence are found in the limit of shorttubes, where an extensive regime of ballistic heat transport prevailing in CNTs of lengthsL . 1µm is independently confirmed. In the limit of long tubes, the PBTT-derived LTCdiverges. Using PBTT and equilibrium MD approaches, we perform numerical calculationsof acoustic phonon lifetimes to clarify the source of divergence. NEMD-derived temperaturedependencies obtained for micrometer-long CNTs and temperatures T ≤ 800 K confirm the1/T behavior of the LTC at moderately high temperatures.Next, we revisit the tube-length dependence of the LTC by use of the relaxation timeapproximation in the PBTT. Through a combination of numerics and analytical consider-ations, we derive exact asymptotic scaling laws of the LTC. In particular, we demonstratethe importance of tensile lattice strain, previously overlooked in the long-standing disputeover tube-length convergence vs divergence of the LTC. Namely, it is proved that, in thelong-tube limit, the relaxation time approximation yields a finite value for stress-free butan infinite value of the LTC in any stretched tube configuration.Lastly, we pursue a matrix inversion approach to solve the linearized phonon transportequation in the framework of the PBTT. Here, it is shown that violations of acoustic sumrules cause spurious convergent behaviors of a length-dependent LTC.iiiLay SummaryCarbon nanotubes are technologically important materials that can be regarded as low-dimensional as only the tube length is principally unconstrained, while the tube diameteris constrained to nanometer length scales. Their reduced dimensionality may promoteunique heat transport properties, but accurately probing such properties at nanometer scalesconstitutes a formidable challenge for experiments. Here, we investigate heat transport inindividual carbon nanotubes using two atomic-scale modeling approaches. One approachconsists in tackling a long-established but hard-to-solve equation of statistical physics. Theother approach is based on molecular dynamics simulation and enables direct monitoring ofheat conducting tubes in space and time. Our results show that the notion of a material-intrinsic thermal conductivity does not apply to carbon nanotubes, whose conductivitycoefficient instead critically depends on the tube length, even at macroscopic length scales.ivPrefaceA version of Chapter 4 has been published asDaniel Bruns, Alireza Nojeh, A. Srikantha Phani, and Jo¨rg Rottler, “Heat trans-port in carbon nanotubes: Length dependence of phononic conductivity fromthe Boltzmann transport equation and molecular dynamics”, Phys. Rev. B 101,195408 (2020).A version of Chapter 5 has been published asDaniel Bruns, Alireza Nojeh, A. Srikantha Phani, and Jo¨rg Rottler, “Nanotubeheat conductors under tensile strain: Reducing the three-phonon scatteringstrength of acoustic phonons”, Phys. Rev. B 104, 075440 (2021).A version of Chapter 6 has been published asDaniel Bruns, Alireza Nojeh, A. Srikantha Phani, and Jo¨rg Rottler, “Commenton ”Ultrahigh Convergent Thermal Conductivity of Carbon Nanotubes fromComprehensive Atomistic Modeling””, Phys. Rev. Lett. 128, 259601 (2022).I performed all simulations and calculations in these publications and wrote the firstdraft of the manuscripts. The preparation of final drafts was assisted by A. Nojeh, A. S. Phani,and J. Rottler.The following computer codes were used. All molecular dynamics simulations were con-ducted with LAMMPS [1]. Calculations in Chapters 4 and 5 were performed with ShengBTE [2].I implemented modifications to ShengBTE to account for four acoustic modes of carbon nan-otubes and wrote additional subroutines to consider finite-length effects and to perform achannel-by-channel decomposition of three-phonon scattering rates. Calculations in Chap-ter 6 were carried out with kALDo [3].vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Atomistic modeling of lattice vibrational heat transport . . . . . . . . . . . 31.3 Carbon nanotubes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.3.1 Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.3.2 The phonon conductivity and its controversies . . . . . . . . . . . . 71.4 Empirical approach for sp2-hybridized carbon . . . . . . . . . . . . . . . . . 141.5 Outline and scope of this work . . . . . . . . . . . . . . . . . . . . . . . . . 152 The Peierls-Boltzmann phonon transport theory applied to carbon nan-otubes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.1 The harmonic lattice approximation . . . . . . . . . . . . . . . . . . . . . . 172.1.1 Phonons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.1.2 Heat capacity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.2 The linearized Peierls-Boltzmann transport equation . . . . . . . . . . . . . 232.2.1 Ballistic limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26vi2.2.2 Anharmonic phonon scattering . . . . . . . . . . . . . . . . . . . . . 272.2.3 Relaxation time approximation . . . . . . . . . . . . . . . . . . . . . 292.2.4 Iterative solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303 Simulating phonon transport in carbon nanotubes . . . . . . . . . . . . . 323.1 The nonequilibrium MD method . . . . . . . . . . . . . . . . . . . . . . . . 333.1.1 Simulation setup and temperature profiles . . . . . . . . . . . . . . . 343.2 Phonon properties from equilibrium MD simulations . . . . . . . . . . . . . 354 Comparison of calculation and simulation based approaches: Effects ofsystem size and temperature . . . . . . . . . . . . . . . . . . . . . . . . . . 394.1 System-size dependence of the thermal conductivity from the PBTT andNEMD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 404.1.1 Comparison with earlier modeling efforts of κ(L,D) . . . . . . . . . 434.1.2 PBTT calculations of κ(L→∞) . . . . . . . . . . . . . . . . . . . . 444.1.3 PBTT calculations of acoustic phonon lifetimes . . . . . . . . . . . . 464.1.4 EMD simulations of acoustic phonon lifetimes . . . . . . . . . . . . . 494.2 Temperature dependence of the thermal conductivity from the PBTT andNEMD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 524.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 545 Low-frequency three-phonon scattering rates and the role of tensilelattice strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 555.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 555.2 Continuum models and anharmonic scattering channels . . . . . . . . . . . 575.3 Channel-by-channel decomposition for stress-free and stretched configurations 595.3.1 Longitudinal modes . . . . . . . . . . . . . . . . . . . . . . . . . . . 595.3.2 Flexural modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 615.3.3 Twisting modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 625.4 Thermal conductivity under the relaxation time approximation . . . . . . . 655.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 685.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 706 Revisiting phonon transport beyond the relaxation time approximation 726.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 726.2 The collision matrix and the necessity of umklapp scattering . . . . . . . . 736.3 Link to the iterative solution . . . . . . . . . . . . . . . . . . . . . . . . . . 746.4 Calculations of κ(L→∞) via direct matrix inversion . . . . . . . . . . . . . 75vii6.4.1 The null space due to energy conservation . . . . . . . . . . . . . . . 766.4.2 Accuracy of third order force constants . . . . . . . . . . . . . . . . 766.4.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 786.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 807 Conclusions and future work . . . . . . . . . . . . . . . . . . . . . . . . . . 817.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 817.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85A Functional form of the Tersoff potential . . . . . . . . . . . . . . . . . . . 97B Spectral densities of acoustic phonons from EMD . . . . . . . . . . . . . 99C Asymptotic analysis of three-phonon scattering channels . . . . . . . . . 102C.1 Decay into two flexural modes . . . . . . . . . . . . . . . . . . . . . . . . . . 103C.2 Coalescence of two flexural modes . . . . . . . . . . . . . . . . . . . . . . . 104C.3 One longitudinal and two twisting modes . . . . . . . . . . . . . . . . . . . 105C.4 Absorption by optical phonons . . . . . . . . . . . . . . . . . . . . . . . . . 105viiiList of TablesTable 1.1 Structural parameters of achiral nanotubes . . . . . . . . . . . . . . . . 7Table 1.2 Review of modeling predictions regarding the existence of a diffusiveheat transport regime in pristine CNTs . . . . . . . . . . . . . . . . . . 12Table A.1 Parametrization of the Tersoff potential for graphitic carbon . . . . . . 98ixList of FiguresFigure 1.1 Chiral indices and structure of single-walled carbon nanotubes . . . . . 6Figure 1.2 Schematic of a heat conducting single-walled carbon nanotube . . . . . 8Figure 1.3 Review of qualitative trends of CNT thermal conductivity as a functionof tube diameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10Figure 1.4 Review of the temperature dependence of CNT thermal conductivity . 13Figure 2.1 Phonon dispersion branches of armchair CNTs . . . . . . . . . . . . . 21Figure 2.2 Displacement patterns of acoustic phonons . . . . . . . . . . . . . . . . 22Figure 2.3 Specific heat capacity under the harmonic lattice approximation . . . . 23Figure 2.1 Schematic of black body boundary conditions for a finite-length CNT . 25Figure 3.1 NEMD simulation setup for calculating the thermal conductivity of CNTs 33Figure 3.2 Steady state temperature profiles from NEMD . . . . . . . . . . . . . . 34Figure 3.1 Finite-temperature phonon dispersion from EMD . . . . . . . . . . . . 37Figure 4.1 System-size dependence of the CNT thermal conductivity from PBTTcalculations and NEMD simulations . . . . . . . . . . . . . . . . . . . . 41Figure 4.2 PBTT calculation results for infinitely long tubes . . . . . . . . . . . . 45Figure 4.3 Lifetimes of acoustic phonons from PBTT calculations . . . . . . . . . 48Figure 4.4 Selection of spectral densities of acoustic phonons from EMD simula-tions of a (10,10) CNT . . . . . . . . . . . . . . . . . . . . . . . . . . . 50Figure 4.5 Lifetimes of longitudinal and twisting modes of CNTs as obtained fromEMD simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51Figure 4.1 Thermal conductivity vs temperature for a (7,7) CNT of varying tubelength . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52Figure 5.1 Acoustic phonon dispersion with varying tensile strain amplitude . . . 58Figure 5.1 Strain-dependent anharmonic scattering rates of LA modes and decom-position into individual three-phonon scattering channels . . . . . . . . 60xFigure 5.2 Strain-dependent anharmonic scattering rates of FA modes and decom-position into individual three-phonon scattering channels . . . . . . . . 62Figure 5.3 Strain-dependent anharmonic scattering rates of TW modes and decom-position into individual three-phonon scattering channels . . . . . . . . 63Figure 5.4 Discontinuity in the strain dependence of scattering rates of long-wavelengthTW modes due to optical intrabranch transitions . . . . . . . . . . . . 65Figure 5.1 Lattice thermal conductivity κRTA for a (4,4) CNT as obtained fromstrain-dependent three-phonon scattering rates . . . . . . . . . . . . . 67Figure 5.2 Relative importance of summed κΣ vs integrated κLW thermal conduc-tivity contributions of acoustic modes as a function of tube length andstrain amplitude . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69Figure 6.1 Lattice-dynamical predictions of low-frequency phonon lifetimes as afunction of anharmonic force constant accuracy . . . . . . . . . . . . . 77Figure 6.2 Room temperature thermal conductivity predictions for an infinite-lengthCNT vs wave number grid resolution . . . . . . . . . . . . . . . . . . . 79Figure B.1 Spectral densities of acoustic phonons (full wave number range) fromEMD simulations of a (3,3) CNT . . . . . . . . . . . . . . . . . . . . . 100Figure B.2 Spectral densities of acoustic phonons (full wave number range) fromEMD simulations of a (4,4) CNT . . . . . . . . . . . . . . . . . . . . . 100Figure B.3 Spectral densities of acoustic phonons (full wave number range) fromEMD simulations of a (7,7) CNT . . . . . . . . . . . . . . . . . . . . . 101Figure B.4 Spectral densities of acoustic phonons (full wave number range) fromEMD simulations of a (10,10) CNT . . . . . . . . . . . . . . . . . . . . 101xiAcknowledgementsThis thesis marks the end of a long journey. A journey that I could not have undertaken ifnot for the support of many great people.First and foremost, I would like to express my gratitude to my academic supervisorsProf. Jo¨rg Rottler, Prof. Alireza Nojeh and Prof. Srikantha Phani. Starting out as a novicein the field of computational condensed matter physics, they gave me the opportunity toventure into new scientific territory and to build my expertise in the field of nanoscalethermal transport. Their guidance and constructive criticism helped me enormously tohone my skills as a researcher, for which I am truly grateful. Acting as my main advisor, Iowe a special thanks to Prof. Jo¨rg Rottler. Knowing that he would always make time whenI needed advice was immensely reassuring.I wish to thank all past and current members of the Rottler group. Discussing variousresearch projects, either in a group format or individually, has been a constant sourceof learning for me. I would like to particularly thank my fellow graduate student DerekFujimoto, who would always be willing to share some valuable hints related to softwareand coding best practices that helped me a lot in my day-to-day work. Many colleagues atUBC became friends outside of work. I would like to thank them for making my time inVancouver very enjoyable and for creating plenty of great memories.For their continuous encouragement, I would like to thank my family and friends backhome. I am especially grateful to my parents, Bernd and Sabine, for supporting me in mystudies abroad and for always being there for me, even from across the miles.Last but not least, I would like to express my deepest gratitude to my wife, Mariana.Her love, understanding and companionship have helped me, above all, to complete thisjourney.xiiChapter 1Introduction1.1 MotivationSeveral technologies of contemporary interest rely to a large extent on materials that eitherblock or move heat with great efficiency [4]. For example, demand is driven by energy con-version and storage devices, such as lithium ion batteries with fast charging capabilities forelectric vehicles, or by microscale computer chips with ever-increasing power densities forportable electronic devices. In recent years, deliberately structuring materials at a micro-scopic scale has become a promising avenue for degrading or improving a material’s abilityto transfer heat [5, 6]. Graphite-based compounds composed of carbon nanotubes (CNTs)in particular have emerged as attractive materials whose nanometer-sized building blocks,seamless cylinders of sp2-hybridized carbon, can be intentionally assembled so that boththermally insulating and highly conductive CNT-based nanostructures can be engineered [7,8]. CNT materials consisting of single and multi-walled CNTs exist in a large variety ofconfigurations with varying degrees of packing densities and nanotube alignment, including,for instance, randomly oriented sponges [9], buckypapers [10], aligned CNT fibers [11], andforests [12]. Experimentally measured thermal conductivities at room temperature rangefrom 0.1 Wm−1K−1 in randomly oriented CNT networks [13] to 766 Wm−1K−1 in high-density, aligned buckypapers [14]. Due to their versatile thermal performance, CNT-basedmaterials have shown prospects for technological use in rather different applications. Amongthese are thermoelectric [15] and thermal management (heat spreading) applications [16],where CNT materials are able to meet low and high thermal conductivity demands, re-spectively. Therefore, from the viewpoint of engineering technologies, it is indispensableto develop a clear understanding of thermal transport at the level of individual single-walled CNTs. Naturally, such an understanding lays the foundation for nanostructuringapproaches, helping the design of CNT-based materials with tailored thermal performance1from the bottom up.Bound together by strong covalent atom-atom interactions, CNTs constitute one-dimen-sional (1D) crystal lattices of cylindrically arranged carbon atoms. Rigid interatomic bonds,a relatively small atomic mass, and crystalline order promote an extremely efficient latticevibrational heat transport along individual tubes. In fact, thermal conductivity values as-signed to the parent materials of CNTs, two-dimensional graphene sheets, are commonlydeemed to be among the highest of all known materials [17], with experimental estimatestypically on the order of 1000s of Wm−1K−1 (see, e.g, Ref. [18] and references therein). Per-forming experimental measurements of heat transfer with nanometer-scale resolution is, ofcourse, a tremendous challenge and, in consequence, our knowledge of lattice heat transportin individual CNTs relies to a great extent on theoretical modeling studies. The first suchtheoretical study dates to the year 2000, when Berber et al. [19], by means of moleculardynamics (MD) simulations, predicted an intrinsic CNT conductivity of 6600 Wm−1K−1at room temperature. Since then, many modeling efforts have been undertaken based onever more refined MD simulations or through numerical calculations in the framework ofthe Peierls-Boltzmann transport theory (PBTT), but some aspects of lattice heat transportin CNTs remain controversial, even on a mere qualitative level to this day. In particular,qualitatively different PBTT and MD-derived predictions can be found in the literaturewith regard to size and temperature effects of CNT heat transport. A review of the currentliterature will be presented in Section 1.3.2.Arguably the most prominent controversy relates to the tube-length dependence of thelattice thermal conductivity. For truly 1D models of monoatomic chains, convincing analyt-ical [20–23] and numerical [24, 25] evidence has been presented for a violation of Fourier’slaw j = κ∆T/L, in the sense that a scaling relation between heat flux j and chain length Lof the form j ∼ L−1+α with 0 < α < 1 can be expected for L→∞, provided that the tem-perature bias ∆T between the chain ends is held constant. This then implies a power-lawdivergence of the lattice thermal conductivity κ ∼ Lα. Although CNTs are often suggestedas a real-world test ground for such low-dimensional heat transport anomalies [26, 27], it isfar from obvious whether 1D model predictions are actually applicable to quasi-1D crystalswhose atoms vibrate in three dimensions. The two main-stream approaches for modelinglattice heat transport in realistic crystals, the PBTT and MD, have thus far struggled toprovide a conclusive answer as to the existence of an intrinsically upper bounded transportcoefficient κ in the long tube limit. As a matter of fact, during the last two decades, thenumber of conflicting MD and PBTT-derived predictions in this respect has continuallygrown (see Table 1.2 below).Given the rather incoherent state of the current literature, the central theme of thisthesis is to re-evaluate the predictions that can be obtained from the two main-stream2modeling approaches of lattice vibrational heat transport when applied to single-walledCNTs. Focusing on PBTT and MD approaches, we seek to establish a better understandingof the reliability of these two state-of-the-art methods at various size and temperature scales.Modeling intricacies related to the reduced dimensionality of CNTs, that have received littleattention previously, shall be discussed in detail.The remainder of this introductory chapter is organized as follows. Section 1.2 providesa high-level perspective of the two modeling approaches, the PBTT and MD, which arepursued throughout this work, with special attention devoted to their conceptual advantagesand disadvantages. Single-walled CNTs, their structural classification and our present-day knowledge of their heat transport properties, are reviewed in Section 1.3. Both thePBTT and MD are considered as parameter-free methods, meaning in principle the onlyprerequisite for predicting lattice heat transfer rates is a model describing the interatomicforces. For the CNTs studied in this thesis, we resort to an empirical force field. Furtherdetails and incentives in regards to this choice are given in Section 1.4. Lastly, in Section 1.5,we state more precisely the outline of this thesis, including scopes and objectives of the mainresearch chapters.1.2 Atomistic modeling of lattice vibrational heat transportIn heat conducting, defect-free crystalline insulators, resistance to the flow of thermal en-ergy originates from two different sources: from interactions with two thermal reservoirskept at different temperatures and from intrinsic lattice anharmonicity. A crystal whoseinteratomic potential energy is purely quadratic in the atomic displacements, forming aharmonic lattice, gives rise to certain normal modes of vibration which, in the absence ofexternal perturbations, do not interact with each other to the extent that per-mode vibra-tional energies remain strictly conserved. The quanta of vibrational energies associated withthese normal modes, termed phonons, are subject to the laws of quantum mechanics. In thepresence of two thermal reservoirs, exact analytical predictions of phononic heat transferrates and temperature profiles can be obtained in the case of harmonic crystals, see, forexample, Ref. [28]. In the case of realistic crystals, however, the challenge of modelingvibrational heat transport consists in accounting for the effects of lattice anharmonicity.Atomic displacements contributing with higher than quadratic order to the interatomic po-tential energy induce phonon-phonon interactions and enable an exchange of energy fromone normal mode to another.To tackle the many-body problem of interacting phonons and related heat conductionphenomena in realistic crystals, two modern computational approaches exist that rely onrather different approximations: the Peierls-Boltzmann phonon transport theory (PBTT)and molecular dynamics simulation (MD) [29]. The PBTT adopts the notion of a kinetic3phonon gas whose transport characteristics in the presence of a temperature gradient arethen described by a Boltzmann transport equation. Here, the quantum nature of phononsis naturally embodied. At the same time, explicit assumptions on the type of phononscattering mechanisms have to be made. Importantly, the PBTT is built upon the notionof a weakly interacting phonon gas, meaning there is no choice but to rely on a perturbativetreatment of anharmonic phonon-phonon interactions. Given the complexity involved inthe interactions of phonons with numerous frequencies, it is standard practice to carry outPBTT calculations within lowest-order anharmonic perturbation theory. MD simulations,on the other hand, obviate the need for explicit assumptions as to the types of phononscattering mechanisms. Here, the dynamics of vibrational normal modes (or phonons) canbe monitored without a truncation of anharmonic effects. In particular, crystal-specific heattransport properties can be inferred by either imposing an actual heat flux within a so-callednonequilibrium MD (NEMD) approach or, alternatively, within a so-called equilibrium MD(EMD) approach, by considering a system at thermal equilibrium and applying some formof fluctuation-dissipation theory. A drawback of MD, however, is the lack of quantumfeatures. Average per-mode vibrational energies follow from classical equipartition statisticsand are therefore systematically overestimated within MD. Such overestimation tends to beparticularly problematic in systems with constituents of relatively low atomic mass. Thesesystems exhibit high-frequency normal mode oscillations which, according to the laws ofquantum mechanics, should be “frozen out” at moderate temperatures.In a sense, an all-encompassing computational framework of vibrational heat transportis nonexistent. Both the PBTT and MD face different limitations and shed light only oncertain aspects of lattice vibrational heat transport. In general, when these two comple-mentary approaches are brought together, a broader range of transport phenomena can becaptured.1.3 Carbon nanotubesAs single-walled CNTs and their lattice vibrational degrees of freedom are the central objectof interest in this thesis, Section 1.3.1 begins by introducing the lattice structure of these1D carbon allotropes. In Section 1.3.2, we then formulate the fundamental problem ofpinpointing CNT thermal conductivities. Here, a review of earlier related research studiesis provided, whereby the existing inconsistencies across different modeling methods areespecially highlighted.41.3.1 StructureIn CNTs, the arrangement of carbon atoms is determined by the way in which a two-dimensional hexagonal graphene lattice is rolled up to form a tube-like hollow cylinder.The number and orientation of hexagonal carbon rings around the cylinder axis are definedby the chiral vector Ch, which acts as the roll-up vector connecting two lattice points ofthe unrolled graphene sheet,Ch = na1 +ma2 ≡ (n,m). (1.3.1)Here, a1 and a2 are the two-dimensional unit vectors of graphene whose lattice parameteramounts to a0 = |a1| = |a2| ≈ 2.46 A˚ [30]. The 1D tube-like lattice of CNTs can then beuniquely specified in terms of the two integers 0 < n and 0 ≤ m ≤ n, called chirality indices.Starting from the honeycomb graphene lattice, the construction of three different nanotubeswith chiral vectors (8, 0), (6, 2) and (5, 5) is exemplified in Figure 1.1. The shown tubes eachgive an example for one of the three different types of nanotubes: zigzag tubes with chiralityindices (n, 0), armchair tubes with (n, n), and chiral tubes with (n,m) and 0 < m < n. Dueto the presence of reflection symmetries, zigzag and armchair nanotubes are also referred toas achiral tubes. As an aside, the chirality (n,m) can be shown to determine the electronictype of CNTs [30]. For example, armchair CNTs are metallic without exception, whereaszigzag types are metallic only if n is a multiple of 3.The diameter of CNTs is given byD =|Ch|pi=a0pi√n2 + nm+m2. (1.3.2)Of importance for a study of lattice vibrations in CNTs is the lattice periodicity along thetube axis, which can be found by drawing the smallest vector T perpendicular to Ch thatconnects two lattice points of the unrolled honeycomb lattice, as is illustrated in Figure 1.1for the case of the (5, 5) tube. The translational vector readsT =n+ 2mga1 − 2n+mga2, (1.3.3)where g stands for the greatest common divisor of (n+2m) and (2n+m). The translationallattice parameter of CNTs follows asa = |T | =√3(n2 + nm+m2)ga0. (1.3.4)A related quantity, that is of interest for a treatment of lattice vibrations in CNTs, isthe number N of atoms per translational unit cell. By dividing the area of one unrolled5Figure 1.1: Chiral indices and structure of single-walled carbon nanotubes. A nan-otube is constructed by means of the chiral vector Ch, which forms the tube’s circum-ference and which specifies structural parameters, such as diameter and lattice periodalong the tube axis. Shown by dashed arrows are the chiral vectors of (8, 0) (red),(6, 2) (dark yellow), and (5, 5) (cyan) nanotubes. Corresponding cylindrical atomicstructures and lattice constants (solid arrows) are sketched in the bottom. Adaptedfrom Ref. [31].6Table 1.1: Structural parameters of achiral nanotubes. a0 denotes the lattice constantof the graphene honeycomb lattice.Diameter Lattice constant Atoms per unit cellD/a0 a/a0 NArmchair (n, n)√3n/pi 1 4nZigzag (n, 0) n/pi√3 4ntranslational unit cell by the area of a six-membered carbon ring and accounting for twoatoms per ring, one findsN = 2 |Ch × T ||a1 × a2| =4(n2 + nm+m2)g. (1.3.5)Throughout this thesis, we restrict our attention to modeling achiral CNTs. For agiven diameter, these tubes have the smallest possible translational unit cell, which leadsto computationally less demanding phonon transport calculations. Simplified expressionsfor achiral tubes resulting from Equations (1.3.1)–(1.3.5) are summarized in Table 1.1.1.3.2 The phonon conductivity and its controversiesHeat transport in solids is described phenomenologically on the basis of Fourier’s law, whichdictates a linear relation between the rate of heat transfer and a temperature gradient,thereby providing the definition of the thermal conductivity κ. In the case of 1D CNTs, byreferring to the situation sketched in Figure 1.2, Fourier’s law can be written in terms ofscalar quantities, readingjz = −κ dTdz. (1.3.6)Here, jz = Jz/Ac stands for the heat flux density corresponding to the thermal energytransferred along the tube axis per unit time Jz through the cross-sectional area Ac, anddT/dz is a temperature gradient equalling the ratio ∆T/L, where ∆T is some temperaturebias between two reservoirs and L is the tube length. While the choice of the tube’s cross-sectional area seems ambiguous, we choose here an annular area Ac = piDδ of thicknessδ = 3.4 A˚ for all subsequently reported values of κ. This choice is traditionally justified bynoting that flat graphene sheets pack in space with interlayer separation ≈ 3.4 A˚ to formgraphite [19]. It further provides a meaningful normalization of the thermal conductance,an extensive material property defined and related to κ as G = Jz/∆T = κAc/L, whichallows one to study the diameter dependence of CNT heat transport.Although most previous efforts aimed at quantifying κ of CNTs relied on computa-tional modeling, some important experimental findings do exist that should be mentioned.7Figure 1.2: Schematic of a heat conducting single-walled carbon nanotube, showingkey geometric parameters. (a) The tube is aligned along the z-axis and a heat fluxJz is induced by two thermal reservoirs of length Lth separated by a distance L. (b)Conventional definition of the tube’s cross-sectional area Ac entering into Fourier’s law,see Equation (1.3.6).Shedding some light on transport phenomena in individual CNTs at varying temperaturesT , first measurements by Hone et al. were performed on crystalline ropes of CNTs [32](8 K ≤ T ≤ 350 K) and on magnetically aligned CNT films [33] (10 K ≤ T ≤ 400 K) tostudy both electrical σe and thermal conductivities κ of these bulk-grown CNT materials.Here, the experimenters observed that the electron contribution κe to the total thermalconductivity κ as inferred from measured σe values and from the Widemann-Franz law fellshort of measured κ values by two orders of magnitude, indicating that phonons (i.e. lat-tice vibrations) are the dominant heat carriers in CNTs, at least as long as T & 10 K.More recently, this observation was experimentally confirmed by Lian et al. [34]. Theycarried out transport measurements on chirality-sorted CNT films and found that heattransfer rates remained insensitive to the chirality dependent electronic type of constituentCNTs (metallic vs semi-conducting), further suggesting phonon dominated heat transportin CNTs with κe κ. The relatively weak electronic contribution κe can be supportedby an argument1 considering phononic and electronic dispersion relations of CNTs whoseband separations are on the order of meV’s and eV’s, respectively [35, 36]: even in metallicCNTs, thermally populated phonon branches greatly outnumber thermally populated elec-tronic conduction bands, assuming no contributions to transport from electronic states be-low the Fermi level. Developing this argument further, phonons and electrons are expectedto contribute equally to the total thermal conductivity only in metallic CNTs at ultra-lowtemperatures T < 5 K [35]. A few experimental studies exist that attempted to pinpointintrinsic κ values of individual CNTs. Pioneering measurements by Yu et al. [37] and Popet al. [38] indicated values on the order of 3000 Wm−1K−1 for tube lengths L ≈ 3µm atroom temperature, consistent with the more recent measurement of Liu et al. [39]. Lee et1This argument assumes idealized ballistic (i.e. scattering-free) heat transport in 1D. As a word of caution,such a simplified picture of phonon and electron-mediated heat transport in CNTs may not be applicableover the whole temperature range, especially so at high T as particle-particle interactions become important.8al. [40] managed to extract values of up to 8000 Wm−1K−1 from their room temperaturemeasurements on CNTs with lengths up to L ≈ 1 mm. Despite these trend-setting reportshinting at extremely efficient phonon-mediated heat transport in CNTs, reliable experimen-tal measurements of intrinsic κ values continue to pose a major challenge. Microfabricatedmeasurements setups inevitably introduce thermal contact resistances and pathways forparasitic heat flows, which hampers the separation of intrinsic heat transport in CNTs fromexternal factors. This issue is aggravated by unavoidable defects in real-world carbon lat-tices. More so, in past experiments on individual CNTs, exact chiralities as well as tubediameters D, and consequently Ac, could not be given with absolute certainty [37–40].By offering perfect control over both atomic structures and environmental factors, com-putational modeling approaches of lattice vibrational transport, either based on the Peierls-Boltzmann phonon transport theory (PBTT) or molecular dynamics (MD) simulation, ap-pear as a natural choice for probing intrinsic CNT heat transfer rates. In reviewing andcomparing some earlier efforts towards computing κ, for a start, it is important to note thatmodeling predictions of exact κ values should not be taken too seriously. This is so, be-cause realistic modeling of heat transport critically depends on an accurate representationof interatomic forces. Existing PBTT and MD studies on CNT heat transport, for the mostpart, have relied on various types of empirical interaction potentials, such as Tersoff [41],Brenner [42], Brenner-II [43], and AIREBO [44] potentials, each of which produces differ-ent atomic force-displacement relationships with varying degrees of nonlinearity, impactingabsolute κ predictions. In fact, for CNTs, a fourfold variability in absolute κ predictionsunder varying interaction potentials was reported [45]. Notwithstanding some doubts aboutquantitative predictions, finding general trends of CNT heat transport makes for the truepower of PBTT and MD approaches. In particular, on a systematic basis, these approachesallow for an investigation of fundamental geometric size and temperature dependencies ofthermal transport, κ = κ(D,L, T ). Surprisingly, in the case of CNTs, a remarkably largenumber of conflicting modeling predictions can be found in the literature, some of whichare reviewed in the following.Diameter dependenceGoing from two-dimensional graphene sheets to CNTs with successively smaller diameters,the question to what extent thermal transport is affected by the curvature of the carbonplane forming the CNT cylinder surface has been addressed previously through PBTTcalculations [47, 48] and nonequilibrium MD (NEMD) simulations [46, 49]. Figure 1.3 showsa compilation of diameter dependent κ predictions normalized against respective limitingvalues obtained for graphene, κGraphene = κ(D →∞). By applying different computationalapproaches to armchair CNTs, Thomas et al. [46] and Yue et al. [47] reported κ ≥ κGraphene9Figure 1.3: Qualitative trends of the CNT thermal conductivity κ at T = 300 K as afunction of tube diameter D as predicted by nonequilibrium MD (NEMD) simulationsand PBTT calculations, assuming different tube lengths L. Predictions by Thomas etal. [46] and Yue et al. [47] are for armchair CNTs and imply a saturation of κ withtube length (“L = ∞”). Calculations by Lindsay et al. [48] contain both chiral andachiral tubes. Simulations by Cao and Qu [49] consider armchair and zigzag tubes ofvarying tube lengths.and a monotonic increase of κ with decreasing tube diameter D. The exact opposite trendwith κ ≤ κGraphene was inferred from NEMD simulations by Cao and Qu [49], who consideredboth armchair and zigzag tubes. Performing PBTT calculations and including arbitrarytube chiralities, Lindsay et al. [48] found κ ≤ κGraphene and a nonmonotonic diameterdependence with a minimum at D ≈ 1.5 nm. Whether the choice of different tube lengthsL helps, in part, to explain these contradicting trends remains an open question.Length dependenceThe simplest description of phonon-mediated heat transport in solids resorts to kinetic gastheory [50, 51], yielding a phonon conductivity expressible asκ ∼ Cv¯Λ. (1.3.7)Here, C is the volumetric specific heat, v¯ is an average sound velocity, and Λ is a phononmean free path capturing the average distance a heat carrying phonon can travel beforeundergoing a particle-like phonon-phonon collision. Based on this notion, one anticipatesthat κ becomes a length dependent quantity at small length scales where L < Λ, withκ(L) exhibiting an increasing trend until at least L ≈ Λ. As part of their experimental10measurements of the CNT thermal conductivity, Hone et al. [32] and Pop et al. [38] gaveroom temperature estimates of Λ = 0.5 − 1.5µm and Λ = 0.5µm, respectively. Hence,for nanometer-sized tube lengths, a sensitive dependence of κ on L can be safely assumed,which is intuitively plausible under the assumption of a micrometer-long intrinsic phononmean free path Λ. However, entering micrometer length scales and beyond, no consensushas been reached about the macroscopic limit, κ(L → ∞), of CNT thermal transport.Due to the low dimensionality of CNTs, it is widely debated whether intrinsic anharmonicthermal resistance alone suffices to establish a finite limiting value κ(L→∞) <∞, markinga diffusive heat transport regime as in ordinary bulk materials, or not, implying an ever-increasing κ as is suggested by observations of superdiffusive heat transport in truly 1Dlattice models [20–25].Table 1.2 shows a selection of previous modeling efforts that arrived at contradictingconclusions concerning the long-tube limit of CNT thermal transport. NEMD studies [49,52, 54, 58], with one exception [46], predicted a power-law divergence κ ∼ Lα with 0 < α < 1for L → ∞. Transport simulation techniques based on equilibrium molecular dynamics(EMD) [55] and homogeneous nonequilibrium molecular dynamics (HNEMD) [57], on theother hand, gave evidence for the existence of a diffusive heat transport regime with κ <∞.Likewise, a finite value was found by means of PBTT calculations [47, 53, 56, 59], someof which [53, 56, 59] further indicated that extremely long tubes with L ' 100µm may berequired to observe a saturation of κ. Within the PBTT formalism, conflicting statements,however, were made regarding the role of higher-order phonon-phonon scattering. Mingoand Broido [53] as well as Lindsay et al. [56] provided numerical evidence that treatingthermal resistance within lowest-order anharmonic perturbation theory causes κ to growwithout bound, thus pointing out the necessity of higher order phonon-phonon scatteringin establishing a diffusive transport regime. This necessity could be stronger in 1D than inhigher dimensions. An inclusion of beyond-lowest-order anharmonic scattering effects, onthe other hand, was not found to be necessary in the more recent works of Yue et al. [47]and Barbalinardo et al. [59].Temperature dependenceThe temperature dependence of κ is expected to show an increasing-to-decreasing trendas the lattice temperature increases which can be rationalized based on Equation (1.3.7)by considering the competing temperature dependencies of the heat capacity C and of thephonon mean free path Λ [50, 51]. A compilation of earlier modeling efforts towards describ-ing κ(T ) of CNTs is shown in Figure 1.4. Inspired by their experimental measurements, Popet al. [38] proposed an empirical function κ(T, L) that goes as ∼ T 2 at low temperatures,peaks at T ≈ 300 K for L ' 1µm, transitions to ∼ T−1 at intermediate temperatures,11Table 1.2: Modeling predictions regarding the existence of a diffusive heat transportregime in pristine CNTs. For diffusive heat transport, κ is expected to saturate to afinite value for sufficiently long tubes, κ(L→∞) <∞. In the case of superdiffusiveheat transport, a power-law divergence, κ ∼ Lα with 0 < α < 1, becomes apparent.Reference Method L(nm) Saturation of κ(L→∞)?Maruyama(2003) [52]NEMD 10− 400 No. κ ∼ Lα with diameter dependentexponent 0.11 ≤ α ≤ 0.27.Mingo and Broido(2005) [53]PBTT 10− 109 Yes. Requires beyond lowest-orderanharmonic phonon scattering andL ' 100µm at T = 316 K.Zhang and Li(2005) [54]NEMD 2− 100 No. κ ∼ Lα with diameter and tempera-ture dependent exponent 0.12 ≤ α ≤ 0.4.Donadio and Galli(2007) [55]EMD 0.5− 200 Yes. In absence of a temperature gra-dient, MD simulations with L ' 20 nmallow one to deduce κ(L→∞) <∞.Lindsay et al.(2009) [56]PBTT 1− 109 Yes. Requires beyond lowest-orderanharmonic phonon scattering andL ' 1 cm at T = 300 K.Thomas et al.(2010) [46]NEMD 200− 1400 Yes. κ(L) = k∞ [1− exp(L/Lc)] withk∞ and Lc increasing with decreasing di-ameter.Cao and Qu(2012) [49]NEMD 100− 2400 No. κ ∼ Lα with α ≈ 0.21 independentof diameter.Yue et al.(2015) [47]PBTT ∞ Yes. Lowest-order anharmonic phononscattering suffices to ensure a saturation.Fan et al.(2019) [57]HNEMD 100 Yes. In absence of a temperature gra-dient, MD simulations with L ' 100 nmallow one to deduce κ(L→∞) <∞.Ray and Limmer(2019) [58]NEMD 5− 1000 No. κ ∼ Lα with α = 0.5.Barbalinardo et al.(2021) [59]PBTT 10− 108 Yes. Lowest-order anharmonic phononscattering and L ' 1 mm suffices to en-sure a saturation.and levels off as ∼ T−2 at high temperatures due to beyond-lowest-order phonon-phononscattering processes. Not captured by this particular function, a Debye model for a one-dimensional phonon gas at low temperatures when applied to a finite-size system2 yieldsκ(T ) ∼ T , which, however, might quickly transition to a superlinear, possibly quadratic,scaling with increasing temperature as high-frequency lattice vibrations of CNTs begin tocontribute to C(T ) [36]. A more problematic assumption incorporated into the empirical2It is assumed that the finite system size imposes an upper bound on the phonon mean free path Λ,rendering Λ temperature independent at low T .12Figure 1.4: Double-logarithmic plot of the temperature-dependent thermal conductiv-ity of CNTs as given by PBTT calculations [53, 60] and by an empirical model [38] thatwas suggested to fit experimentally measured data for a tube of diameter D ≈ 1.7 nm.CNTs of chirality (10,0) and (10,10) were considered in Ref. [53] and in Ref. [60],respectively.function of Ref. [38] concerns the tube length dependence of κ. In crystalline solids, resistive3phonon-phonon scattering processes are generally assumed to exhibit an exponential scalingat low temperatures, implying Λ ∼ exp(T0/T ) with some material-specific positive T0 [51].In the case of pristine CNTs, such behavior and its impact on the tube length dependenceof κ are not included in the model of Ref. [38], which instead assumes a length saturationof κ at L ≈ 10µm irrespectively of temperature. Using empirical interatomic interactionpotentials, PBTT calculations of κ(T ) were performed at different levels of sophistication byMingo and Broido [53] and by Lindsay and Broido [60]. As shown in Figure 1.4, assumingL = 3µm, both these studies agree on a peak location of κ(T ) at T ≈ 150 K, but makedifferent predictions about the high temperature scaling of κ(T ). As opposed to Ref. [60]which predicts κ(T ) ∼ T−1, the authors of Ref. [53] considered approximately the effects ofhigher order phonon-phonon scattering processes, predicting κ ∼ T−2.It is important to note here that the nonmonotonic behavior of κ(T ), leading to a max-imum value at finite temperatures, is a quantum feature whose description defies standardMD simulation approaches. In the case of CNTs, some attempts towards correcting classi-cal MD predictions of κ were made [61], but the applicability of such quantum corrections3The term “resistive” refers to momentum-destroying umklapp scattering processes. These are to bedistinguished from momentum-conserving normal processes which, in the absence of umklapp processes, donot cause resistance to heat flow, as will be further explained in Chapter 6.13remains generally debated [62]. Still, since no a priori assumptions about phonon-phononscattering have to be made, it is the high temperature regime where the predictive power ofMD simulations unfolds. To the best of our knowledge, however, a systematic investigationof the high temperature dependence of κ, in the range up to T = 1000 K, has not yet beenperformed by means of MD simulations.1.4 Empirical approach for sp2-hybridized carbonAtomistic modeling of lattice vibrational heat transport requires knowing a material’s latticeenergy U as a function of atomic coordinates. While such calculations are, in principle,within the reach of density functional theory, there are two major pitfalls of parameter-free ab initio approaches when modeling heat transport of individual CNTs. Within MD,probing high-κ as opposed to low-κ materials creates an increased overall computationaloverhead, as converged κ predictions tend to require larger simulation time and lengthscales [63]. Due to the computational cost, existing ab initio MD studies that examinedthermal transport in a variety of different materials were limited to a few hundred atomswith time-evolved trajectories not exceeding the picosecond scale [64–67]. Hence, for highlyconductive materials like CNTs, such efforts might be computationally prohibitive. In thecontext of the PBTT, another intricacy of ab initio calculations relates to the approximationof lattice anharmonicity. Anharmonic force constants (meaning spatial derivatives of U ofthird or higher order), which are an essential input to the PBTT method, are typicallyestimated by finite differences when derived from ab initio calculations [68]. Here, a supercellis considered and, for feasible calculations, small but finite force constants linking atomsoutside a predefined cutoff radius are truncated. Approximating force constants in thisway inevitably leads to a violation of translational and rotational invariances of the crystallattice, which in turn necessitates a careful symmetrization of ab initio force constants [69].Both these hurdles can be alleviated by employing an empirical interatomic interactionpotential, U ≡ U emp, which is the method of choice to model CNT heat transport inthis thesis. In particular, we resort to the bond-order potential developed by Tersoff forcovalently bonded carbon [70, 71] and adopt its optimized parametrization given by Lindsayand Broido [72]. The latter authors re-parametrized Tersoff’s original parameter set tomore accurately reproduce experimentally measured in-plane phonon-dispersion data forgraphite.4 As pointed out in Ref. [72], this new optimized Tersoff potential provides amore realistic description of thermal transport in sp2-hybridized carbon lattices such asCNTs. The functional form of U emp is described in Appendix A. When equipped with4It is generally expected that experimentally measured in-plane phonon dispersion relations of graphiteare almost identical to their counterparts in individual graphene layers, given that interlayer van der Waalsinteractions are relatively weak in comparison to sp2-hybridized carbon-carbon bonds, see Ref. [73].14modern high performance computing resources, optimized implementations of the Tersoffpotential enable MD simulations with hundreds of thousands of atoms over nanosecondtime scales, see, for example, Ref. [74]. Moreover, in the context of the PBTT, having acontinuously differentiable function for U , with well-defined short-range atomic interactionradius, significantly simplifies obtaining anharmonic force constants, in that, for one, theoverall number of nonzero force constants is significantly reduced and, more so, in thatcorrection procedures as detailed in Ref. [69] become redundant. On the downside, itneeds to be mentioned that the optimized Tersoff potential [72], like most other empiricalpotentials in general, was fitted with the objective of reproducing realistic harmonic latticeproperties, in this case of sp2-bonded carbon lattices, with no particular emphasis given toanharmonic terms influencing thermal resistance.1.5 Outline and scope of this workIn this thesis, we consider defect-free single-walled CNTs in the adiabatic approximation ofBorn and Oppenheimer [75], focusing exclusively on nuclear degrees of freedom. More specif-ically, for all calculations and simulations presented in this work, the nuclei are assumedto interact through the optimized Tersoff potential (see Section 1.4), which constitutes ourchoice for the Born-Oppenheimer energy surface. Regarding nonequilibrium heat transportscenarios, we further take it for granted that CNTs operate in a steady-state linear responseregime. That means a linear relationship between heat flux j (response) and temperaturebias ∆T (perturbation) as per Fourier’s law is assumed to hold true. More specifically,referring to Figure 1.2, we proceed on the assumption that the temperature along the tubeaxis can be described by an average T plus a small correction δT (z) leading to a constanttemperature gradient outside of the thermal reservoirs. As will become clear in Chapters 2and 3, approximating δT (z) to be strictly linear simplifies a treatment of finite-tube-lengtheffects in that effective temperature-dependent conductivity coefficients κ(T ) for varyingtube lengths L can be defined globally. We note that such an approximation can be relaxedby generalizing Fourier’s law to a nonlocal relation, jz(z) =∫dz′κ(z−z′)dT (z′)/dz′ [76, 77],which in turn allows one to account for possible nonlinear spatial variations of δT (z) [78].Applying nonlocal theories of phonon-mediated heat to realistic materials is an active areaof research [79–81] beyond the scope of this thesis.Our calculations in the framework of the PBTT revolve around solving the linearizedPeierls-Boltzmann transport equation for phonons, taking into account the full spectrum oflowest-order anharmonic phonon-phonon interactions (i.e three-phonon processes). Chap-ter 2 outlines the theoretical background of the PBTT. MD simulations are carried out intwo different settings. Direct simulations of heat conducting CNTs are performed withinan NEMD velocity-swapping approach. Individual normal mode vibrations of CNTs are15studied by the use of EMD simulations. Both MD techniques are explained in Chapter 3.Chapters 4–6 present our main findings. Their outline is as follows. In Chapter 4, weconcentrate on free-standing CNTs and conduct a side-by-side comparison of numericallyobtained thermal conductivity coefficients κ using PBTT calculation and MD simulationbased approaches. The length, diameter, and temperature dependence of κ is discussed.Here, we put special emphasis on the tube-length dependence of thermal transport byconsidering numerically the related issue of weakly damped long-wavelength phonon modesin both PBTT and MD modeling frameworks.In Chapter 5, we aim to provide deeper insights into PBTT-derived conductivity predic-tions. While standard PBTT modeling approaches rely heavily on numerics, we investigateto what extent numerical calculations can be supported by analytical considerations. Thischapter is inspired by an earlier work of Bonini et al. [82], in which the authors provideda detailed analysis of the issue of long waves in heat conducting graphene sheets. In thecase of these two-dimensional crystals, Bonini et al. further included and demonstrated theimportance of tensile strain effects. We do so here for 1D nanotube heat conductors.The purpose of Chapter 6 is to comment on a recently reported break-through in themodeling of low-dimensional heat transport. This concerns the work of Barbalinardo etal. [59] who devised a new numerical routine to calculate the CNT thermal conductivityfrom a linearized Peierls-Boltzmann equation with three-phonon interactions, arriving atthe conclusion that normal (i.e. diffusive) heat transport, κ <∞, in ultra-long free-standingtubes can be deduced. A source of systematic calculation errors in the framework of thePBTT is discussed.In Chapter 7, we conclude and give a final assessment of PBTT and MD-based modelingapproaches. In addition, a perspective on open problems and future research directions isoffered.16Chapter 2The Peierls-Boltzmann phonontransport theory applied to carbonnanotubesThe phonon picture of crystalline solids as well as the Peierls-Boltzmann phonon transporttheory (PBTT) are treated in several classic textbooks [50, 83, 84]. In this chapter, webriefly review the theory of phonons and the PBTT specifically tailored to CNTs. We beginby introducing the harmonic lattice approximation in Section 2.1, which we subsequentlyapply to showcase phonon dispersion relations and heat capacities of CNTs. In doingso, we establish some nomenclature and some elementary lattice-dynamical properties ofCNTs that are essential for heat transport calculations in the framework of the PBTT. InSection 2.2, we focus on the fundamentals of the PBTT. In particular, we follow the work ofMingo and Broido [53] to derive a linearized Boltzmann transport equation for finite-lengthCNTs. We then describe how the so-obtained transport equation enables a calculation ofthe CNT thermal conductivity under different approximations, thus providing a detailedaccount for the numerical results presented in later chapters.2.1 The harmonic lattice approximationThe starting point for modeling lattice vibrational heat transport of carbon nanotubes(CNTs) is the lattice HamiltonianHlatt = K ({r˙}) + U ({r})=∑l,p12mCr˙2l,p + U(r1,1, r1,2, . . . , r1,N , r2,1, . . . , rN,N ).(2.1.1)17Here, mC is the atomic mass of carbon, and instantaneous atomic velocities entering intothe kinetic energy K, r˙l,p = r˙l,p(t), as well as positions entering into the potential energyU , rl,p = rl,p(t), are labeled by two indices l = 1, . . . , N and p = 1, . . . ,N running overtranslational unit cells and atomic sites within each unit cell, respectively. In solids, a sta-tionary equilibrium position r0(l, p) can be assigned to each atom. Dynamical lattice prop-erties are then more readily expressed in terms of time-dependent displacement variablesul,p(t) = rl,p(t)−r0(l, p). With these new variables, the potential energy in Equation (2.1.1)can be expanded in a Taylor series in powers of ul,p = (uxl,p, uyl,p, uzl,p)T,U ({u}) = U0 + 12∑l1,p1,l2,p2∑α1,α2Φα1,α2l1,p1,l2,p2uα1l1,p1uα2l2,p2+∑s=31s!∑l1,p1,...,ls,ps∑α1,...,αsΦα1,...,αsl1,p1,...,ls,psuα1l1,p1· · ·uαsls,ps ,(2.1.2)with sth-order force constants defined asΦα1,...,αsl1,p1,...,ls,ps =∂sU∂uα1l1,p1 · · · ∂uαsls,ps∣∣∣∣∣uα1l1,p1=···=uαsls,ps=0. (2.1.3)In Equation (2.1.2), U0 is the cohesive lattice energy at equilibrium which has no impact onlattice-dynamical properties and will be dropped in the following. We further note that thefirst order term of the expansion is zero since we expand the potential around its minimumenergy configuration.Assuming atomic displacements to be small and truncating the Taylor expansion (2.1.2)after the second order term constitutes the so-called harmonic lattice approximation. Underthis approximation, the motion of each individual carbon atom can be decomposed into aset of noninteracting waves that travel along the tube axis with wave number k,ul,p(t) =∑k,νep,k,ν ei[k la−ωk,νt], (2.1.4)where a denotes the translational lattice constant of CNTs. The discrete periodicity alongthe tube axis brings with it a lower bound for possible wavelengths, λmin = 2a, so that|k| ≤ 2pi/λmin = pi/a. At the same time, by applying periodic boundary conditions to a tubeconsisting of N unit cells, only N distinct values of k are commensurate with ul,p = ul+N,p,implying a discretizationk =2piNank, with integers |nk| ≤ N/2. (2.1.5)18An additional index ν in Equation (2.1.4) relates to possible wave polarizations and extendsover ν = 1, . . . , 3N , matching the total number of atomic degrees of freedom 3NN for aCNT embedded in three dimensional space. Wave polarization vectors ep,k,ν and frequenciesωk,ν in Equation (2.1.4) are determined by solving the equations of motion of the harmoniclattice,mCu¨αl,p = −∑l′,p′∑α′Φα,α′l,p,l′,p′uα′l′,p′ . (2.1.6)By plugging an individual traveling wave as defined in Equation (2.1.4) into Equation (2.1.6),one arrives at the following matrix equationω2k,ν ep,k,ν =∑p′Ωp,p′(k) · ep′,k,ν , (2.1.7)where the details of atom-atom interactions are encoded in the dynamical matrixΩα,α′p,p′ (k) =1mC∑l′Φα,α′l,p,l′,p′ eik (l′−l)a, (2.1.8)which can be calculated with respect to an arbitrary unit cell labeled by l due to thetranslational invariance of the lattice potential energy. By diagonalizing the 3N × 3Nmatrix Ω(k), we obtain for each value of k a discrete set of 3N real-valued frequencies ων(k)labeled by ν, with complex-valued wave polarization vectors, ep,ν(k) = e∗p,ν(−k), definingwave-specific atomic displacement patterns. A phonon mode is generally associated with aneigenpair[ωk,ν , (e1,k,ν , . . . , eN ,k,ν)T], and finding the complete spectrum of phonon modesfor a given CNT is an essential preliminary calculation step in the PBTT formalism.Before we proceed to discuss actual data of phonon spectra in the case of CNTs, it isinstructive to introduce some finite temperature properties of the harmonic lattice. Thisis conveniently done by treating phonon modes as quasi-particles to which the tools ofquantum statistics can be applied. To do so explicitly, we may write atomic displacementsand velocities in terms of so-called normal mode coordinates Qk,ν and Pk,ν ,ul,p(t) =1√N mC∑k,νe˜p,k,ν eik laQk,ν(t),u˙l,p(t) =1√N mC∑k,νe˜p,k,ν eik la Pk,ν(t),(2.1.9)with polarization vectors inferred from Equation (2.1.7) and with tildes signifying normal-ization to unit length,∑p |e˜p,k,ν |2 = 1. Expressing atomic displacements and velocitiesin terms of normal mode coordinates, the lattice Hamiltonian (2.1.1) takes on a simplified19form,Hlatt = Hharm +Hanh=∑k,ν12Pk,νP−k,ν +12ω2k,νQk,νQ−k,ν +Hanh.(2.1.10)Comprising the kinetic energy K and quadratic terms in ul,p, the harmonic lattice Hamil-tonian Hharm transforms into a sum over noninteracting harmonic oscillators, representingindependent phonon modes, while anharmonic terms of cubic or higher order in ul,p aresubsumed into Hanh. Focusing on Hharm, the quantum, particle-like character of phononmodes is now formally introduced by following the standard second quantization procedure,where one writesQˆk,ν =√~2ωk,ν(aˆk,ν + aˆ†−k,ν),Pˆk,ν = i√~ωk,ν2(aˆ†−k,ν − aˆk,ν),(2.1.11)with bosonic annihilation and creation operators aˆk,ν and aˆ†k,ν , respectively. Using Equa-tion (2.1.11), the harmonic lattice Hamiltonian becomesHˆharm =∑k,ν~ωk,ν(aˆ†k,ν aˆk,ν +12), (2.1.12)which is the desired form that allows one to derive thermodynamic functions of the harmoniclattice.In the context of lattice heat transport, we quote the most important ensemble averages.The Bose-Einstein distribution function,n0k,ν(T ) = 〈aˆ†k,ν aˆk,ν〉 =(e~ωk,νkBT − 1)−1, (2.1.13)determines how much energy at lattice temperature T is stored in each phonon mode,Ek,ν(T ) = ~ωk,ν(n0k,ν +12). (2.1.14)Deriving Equation (2.1.14) with respect to T yields the per-mode contributions to the totallattice heat capacity,ck,ν(T ) =∂Ek,ν∂T=~2ω2k,νkBT 2n0k,ν(n0k,ν + 1). (2.1.15)Classical high temperature limits can be derived from Equations (2.1.13)–(2.1.15) under the20Figure 2.1: Phonon dispersion branches of armchair CNTs (n, n) with diameters equalto 5.43 A˚, 9.49 A˚, and 13.56 A˚ (from left to right). Colored and grey lines show thefour acoustic and the 12n− 4 optical branches, respectively, some of which are doublydegenerate. Shown by a grey dashed line is the radial breathing mode which hybridizeswith the longitudinal acoustic mode, causing an avoided crossing of respective phononbranches.assumption ~ωk,ν kBT .2.1.1 PhononsTo demonstrate both some general and some diameter specific features of phonon spectrain CNTs, Figure 2.1 shows the dispersion branches of (4, 4), (7, 7) and (10, 10) armchairtubes, as obtained from a numerical diagonalization of the dynamical matrix (2.1.8). Onlythe region 0 ≤ k ≤ pi/a is plotted since ων(k) = ων(−k). Given the stiffness of sp2-hybridized carbon-carbon bonds and the relatively low atomic mass of carbon, phonons inCNTs, as compared to other materials, reach relatively high frequencies of up to ≈ 50 THz,amounting to quasi-particle energies of up to ≈ 200 meV. Signifying the room temperaturethermal energy of ≈ 25.8 meV, the dot-dashed line in Figure 2.1 indicates that only a21Figure 2.2: Displacement patterns of a (10,10) CNT induced by the three differenttypes of acoustic phonon polarizations.small portion of low-frequency phonons allow for the classical approximation ~ων,k kBTat T = 300 K, which implies that lattice-dynamical properties of CNTs are governed byquantum effects even at moderate to high temperatures. Crossing dispersion branches νand ν ′ at the boundary k = pi/a exhibiting opposite, nonzero slope, ω′ν = −ω′ν′ , are aconsequence of screw-axis symmetry operations [85, 86], of which there are generally fewerin smaller diameter CNTs [31].In the low-frequency regime, four acoustic phonon branches exist with ων(|k| → 0) = 0.These are one longitudinal acoustic (LA) branch, one twisting acoustic (TW) branch, andtwo degenerate flexural acoustic (FA) branches. Inferred from phonon eigenvectors definedin Equation (2.1.7), corresponding displacement patterns of acoustic phonons in CNTsare visualized in Figure 2.2. LA and TW branches exhibit a linear dispersion relation,ωLA/TW ∼ k for |k| → 0, whereas FA branches disperse quadratically away from k =0, ωFA ∼ k2. Comparing the FA mode dispersion among the three different tubes inFigure 2.1, it can be seen that the prefactor of the quadratic dispersion increases withincreasing diameter, which is indicative of an increasing bending stiffness in larger diametertubes.2.1.2 Heat capacityEquipped with a list of allowed phonon frequencies, the total lattice heat capacity of CNTsat constant volume can be straightforwardly calculated by summing over the per-modeheat capacities according to Equation (2.1.15), CV =∑k,ν ck,ν . In Figure 2.3, we show thespecific lattice heat capacity, CmV = CV /NNmC, of CNTs with varying diameter resultingfrom a harmonic lattice calculation. A noticeable diameter dependence of CmV becomesapparent only at ultra-low temperatures, where thermal energies may excite a larger numberof low-energy FA modes in smaller diameter tubes. It can further be noted from Figure 2.3that, unless T > 2000 K, the classical Dulong-Petit approximation of ck,ν = kB, leading toCm,clV = 3kB/mC, largely overestimates the true lattice heat capacity of CNTs. For instance,focusing on the CNT (4,4), we find that the ratio CmV /Cm,clV amounts to 0.31, 0.83 and 0.95at temperatures of 300 K, 1000 K, and 2000 K, respectively.22Figure 2.3: Specific heat capacity of CNTs under the harmonic lattice approximation.The dot-dashed line indicates the classical Dulong-Petit limit Cm,clV = 3kB/mC. In theinset, the low-temperature region is plotted on log-log scale.2.2 The linearized Peierls-Boltzmann transport equationTaking into account anharmonic phonon-phonon interactions, Peierls [87] gave a theoreti-cal treatment of phonon-mediated transport in terms of a Boltzmann transport equation,providing a formalism to connect microscopic atom-atom interactions with a material’smacroscopic lattice heat transport coefficient κ. The basic idea of the Peierls-Boltzmannphonon transport theory (PBTT) is to derive the nonequilibrium phonon distribution func-tions nk,ν in the presence of a temperature gradient ∇T from a balance equation that, inits most general form [50, 83], readsn˙k,ν(r, t) = −vk,ν · ∇nk,ν(r, t) + n˙k,ν(r, t)∣∣∣scatt. (2.2.1)Here, the net change per unit time n˙k,ν(r, t) of phonons with wave vector k and polarizationν, located in space at r, is given by a drift contribution of phonons traveling with groupvelocity vk,ν = ∂ωk,ν/∂k along a temperature gradient and by a scattering contributionwhich, in the case of pristine, insulating crystals, originates from anharmonic phonon-phonon collisions. Since the unknown functions nk,ν for each individual phonon mode arecoupled in a nontrivial manner through the scattering term, some simplifying assumptionsare needed in order to deal with Equation (2.2.1). For an infinitely extended, spatiallyinvariant solid, at some average lattice temperature T , one typically considers an expansion23about the Bose-Einstein distribution function n0[X] =(eX − 1)−1 of the form [83, 84]:nk,ν = n0 [~ωk,ν/kBT −Ψk,ν ]≈ n0k,ν −∂n0 [X]∂X∣∣∣X=~ωk,ν/kBTΨk,ν= n0k,ν + n0k,ν(n0k,ν + 1)Ψk,ν ,(2.2.2)where the unknowns Ψk,ν , assumed to be linear in the temperature gradient, measuredeviations from n0[~ωk,ν/kBT ] = n0k,ν . In this way, under steady state conditions n˙k,ν = 0,Equation (2.2.1) can be transformed into a system of linear equations for Ψk,ν , whosesolution relates to a macroscopic heat flux asj =1NV∑k,ν~ωk,νvk,νnk,ν=1NV∑k,ν~ωk,νvk,νn0k,ν(n0k,ν + 1)Ψk,ν ,(2.2.3)with NV standing for the volume of a solid with unit cell volume V [50, 84].As pointed out by Mingo and Broido [36, 53], applying the PBTT to finite-length CNTsunder steady state conditions, a linearized transport equation has to be formulated so as toaccount for boundary conditions at the tube ends. This can be achieved by resorting to theexpansion (2.2.2) and by allowing for a spatial variation of Ψk,ν = Ψk,ν(z) along the tubeaxis. Considering a heat conducting CNT of length L as depicted in Figure 1.2, subject to asmall temperature bias ∆T = Th− Tc T = (Th + Tc)/2, black body boundary conditionscan be expressed asnσ→(0) = n0 [~ωσ/kBTh] ≈ n0σ +~ωσkBT 2n0σ(n0σ + 1)∆T2, (2.2.4a)nσ←(L) = n0 [~ωσ/kBTc] ≈ n0σ −~ωσkBT 2n0σ(n0σ + 1)∆T2, (2.2.4b)where σ ≡ (k, ν) is introduced for brevity in notation, and arrows indicate phonons with op-posite velocities traveling to the right (→) and to the left (←). By referring to the situationsketched in the top panel of Figure 2.1, a simplified expression for the drift contribution in24Figure 2.1: Schematic of black body boundary conditions in the case of a CNT withlength L. Red solid, blue dashed, and green solid lines indicate weak, intermediate andstrong phonon-phonon scattering, respectively. See text for explanations. Top panel:Linear approximation of the spatial dependence of nσ(z). Bottom panel: Spatiallyresolved average temperature in the nonequilibrium steady state.Equation (2.2.1) can now be obtained through a linear approximation of nσ(z), setting1dnσ→(z)dz≈ nσ→(L/2)− nσ→(0)L/2, (2.2.5a)dnσ←(z)dz≈ nσ←(L)− nσ←(L/2)L/2. (2.2.5b)Combining Equations (2.2.2), (2.2.4) and (2.2.5) and using the fact that vk,ν = −v−k,ν , wearrive at the steady state transport equation for finite-length CNTs [36, 53],vσ~ωσkBT 2∆TL=2|vσ|LΨσ(L/2)− 1n0σ(n0σ + 1)n˙σ∣∣∣scatt. (2.2.6)For infinitely extended systems, the drift term is often simplified via the so-called local1The assumption of linearity of nσ(z) can be relaxed by considering a “full” solution of a space-dependentlinearized Peierls-Boltzmann transport equation [79–81], which comes at the price of significantly increasedcomputational complexity.25thermal equilibrium assumption [50],dnσ(z)dz≈ ddzn0[~ωσkBT (z)]≈ − ~ωσkBT 2n0σ(n0σ + 1)∆TL, (2.2.7)which can be checked to reproduce Equation (2.2.6) in the limit L → ∞, where the firstterm on the right-hand side vanishes. In such a scenario, as illustrated by green linesin Figure 2.1, we expect phonon-phonon interactions to be sufficiently strong such thatphonons traveling between the reservoirs thermalize to a local equilibrium temperature,giving rise to a continuous temperature profile T (z) along the tube axis. If, on the otherhand, phonons scatter only weakly among each other, we may use Equation (2.2.6) tocalculate heat transfer rates in the ballistic limit.2.2.1 Ballistic limitFor L → 0, the phonon-phonon scattering term in Equation (2.2.6) can be assumed negli-gible, and we haveΨσ =~ωσkBT 2vσ|vσ|∆T2, (2.2.8)such that, according to Equation (2.2.3), a heat flux density follows asjz =1NaAc∑σ~2ω2σkBT 2n0σ(n0σ + 1)|vσ|2∆T=1NaAc∑σcσ|vσ|2∆T,(2.2.9)where cσ is the per-mode heat capacity as defined in Equation (2.1.15). This expressiondescribes a ballistic transport regime, in that phonon travel without hindrance along thetube axis, thus violating the notion of a local thermal equilibrium as depicted by red linesin Figure 2.1. For fully ballistic transport, the thermal conductance G = jzAc/∆T of CNTsbecomes tube-length independent. Comparing Equation (2.2.9) with Fourier’s law (1.3.6),we may define a ballistic, tube-length dependent thermal conductivity asκbal =LNaAc∑σcσ|vσ|2. (2.2.10)For a given tube length L and lattice temperature T , this form provides a quantum upper-bound to phonon-mediated heat transport in CNTs. Numerical values of heat transfer ratesof CNTs derived from Equation (2.2.10) will be discussed in Chapter 4. We note in passingthat Equations (2.2.9) and (2.2.10) can equivalently be derived in the Landauer formalismfor quantum transport [88, 89].262.2.2 Anharmonic phonon scatteringIn order to incorporate thermal resistance due to anharmonic phonon scattering into thePBTT, one considers an expansion of scattering processes of the formn˙σ∣∣∣scatt=∑σ2,σ312(P(3)σ2+σ3→σ − P(3)σ→σ2+σ3)+ P(3)σ3→σ+σ2 − P(3)σ+σ2→σ3+∑σ2,σ3,σ416(P(4)σ2+σ3+σ4→σ − P(4)σ→σ2+σ3+σ4)+12(P(4)σ2+σ3→σ+σ4− P(4)σ+σ4→σ2+σ3)+12(P(4)σ2→σ+σ3+σ4 − P(4)σ+σ3+σ4→σ2)+∑σ2,σ3,σ4,σ5. . .(2.2.11)where the P(s)i→f are collision rates with the involvement of s phonons, and the prefactors of1/2 and 1/6, etc., serve to prevent double counting [84, 90]. Individual collision rates aredetermined perturbatively via Fermi’s golden rule asP(s)i→f =2pi~∣∣∣〈f |Hˆ(s)anh|i〉∣∣∣2 δ(Ei − Ef ), (2.2.12)where Hˆ(s)anh is the sth-order anharmonic lattice Hamiltonian which is inferred from theTaylor series (2.1.2) and written in terms of phonon creation and annihilation operators,Hˆ(s)anh ∝(aˆk1,ν1 + aˆ†−k1,ν1). . .(aˆks,νs + aˆ†−ks,νs). (2.2.13)Evaluating the scattering term as per Equation (2.2.11) constitutes one of the computationalbottlenecks of the PBTT. In this thesis, we restrict ourselves to lowest-order anharmonicperturbation theory, which means that we retain only three-phonon collision rates in Equa-tion (2.2.11). In this case, one finds (see Refs. [50, 83]):n˙σ∣∣∣scatt= −∑σ2,σ312[nσ(nσ2 + 1)(nσ3 + 1)− (nσ + 1)nσ2nσ3 ]L−σ,σ2,σ3+ [nσnσ2(nσ3 + 1)− (nσ + 1)(nσ2 + 1)nσ3 ]L+σ,σ2,σ3 .(2.2.14)Here, the L±’s are the intrinsic, temperature independent three-phonon absorption (super-script +) and decay (superscript −) rates, which, for a given phonon triplet (σ, σ2, σ3),depend on frequencies (ωσ, ωσ2 , ωσ3), eigenvectors (e˜σ, e˜σ2 , e˜σ3), and third order force con-27stants Φα1,α2,α3l1,p1,l2,p2,l3,p3 , reading:L±σ,σ2,σ3 =pi~4Nm3C∣∣V ±σ,σ2,σ3∣∣2 δ(ωσ ± ωσ2 − ωσ3)ωσωσ2ωσ3 ∆(k ± k2 − k3), (2.2.15)with three-phonon coupling coefficients, V ±σ,σ2,σ3 = Vν,ν2,ν3(k,±k2,−k3), given byVν,ν2,ν3(k, k2, k3) =∑p1,l2,p2,l3,p3∑α1,α2,α3Φα1,α2,α3l1,p1,l2,p2,l3,p3× e˜α1p1,k,ν e˜α2p2,k2,ν2e˜α3p3,k3,ν3 eik2 l2aeik3 l3a,(2.2.16)where, just like in Equation (2.1.8), an arbitrary unit cell labeled by l1 may be chosenas reference, due to the translational invariance of the lattice potential energy. Whilethe conservation of energy of three-phonon collisions is strictly enforced through the deltafunction in Equation (2.2.15), phonon wave numbers are conserved only modulo 2pi/a, asreflected by the function∆(k) ={1 if k = 2pi j/a, with j = 0,±10 otherwise.(2.2.17)Collisions with j = 0 and j = ±1 are referred to as normal and umklapp processes, respec-tively. The distinction of normal and umklapp collisions plays a special role in the PBTTand will be discussed in greater detail in Chapter 6.As mentioned in Section 2.2, carrying out actual transport calculations based on Equa-tion (2.2.1) becomes possible by linearizing the nonequilibrium distribution functions nσabout the equilibrium Bose-Einstein distributions n0σ. By plugging the expansion (2.2.2)into Equation (2.2.14) and by dropping terms of higher than linear order in Ψσ, the scat-tering term writesn˙σ∣∣∣scatt= −∑σ2,σ312[Ψσ −Ψσ2 −Ψσ3 ] n0σ(n0σ2 + 1) (n0σ3 + 1)L−σ,σ2,σ3+ [Ψσ + Ψσ2 −Ψσ3 ]n0σn0σ2(n0σ3 + 1)L+σ,σ2,σ3 .(2.2.18)Finally, to achieve a linear relation between heat flux and temperature gradient, one definesΨσ = Fσ∆TL, (2.2.19)which, when combined with Equations (2.2.6) and (2.2.18), yields a computable form of thetransport equation for finite-length CNT under the inclusion of anharmonic three-phonon28scatterings, reading:vσ~ωσkBT 2=2|vσ|LFσ +∑σ2,σ312[Fσ − Fσ2 − Fσ3 ] Γ−σ,σ2,σ3+ [Fσ + Fσ2 − Fσ3 ] Γ+σ,σ2,σ3 ,(2.2.20)where the Γ±’s represent finite-temperature scattering amplitudes, defined asΓ±σ,σ2,σ3 ={n0σ2 − n0σ3n0σ2 + n0σ3 + 1}L±σ,σ2,σ3 , (2.2.21)with the top and bottom row in curly brackets denoting absorption (superscript +) anddecay (superscript −) amplitudes, respectively. Once the spectrum of nonvanishing scatter-ing amplitudes, as determined by phonon dispersion relations and three-phonon couplingcoefficients V ±σ,σ2,σ3 , is known, Equation (2.2.20) may be solved for the unknown coefficientsFσ at different levels of accuracy. Solutions are commonly obtained within a so-called relax-ation time approximation [84] or, more exactly, through an iterative solution approach [2,91, 92].2.2.3 Relaxation time approximationArising from anharmonic phonon collisions in a steady state heat transport scenario, Equa-tion (2.2.20) describes a complex interdependence of nonequilibrium phonon distributionfunctions as measured by the coefficients Fσ. A major simplification of the phonon trans-port picture can be achieved by assuming that each individual phonon relaxes back toits equilibrium distribution with all other phonon modes simultaneously residing in ther-mal equilibrium, thus ignoring nonequilibrium excitations of scattering partners. This as-sumption constitutes the relaxation time approximation (RTA). Mathematically, the RTAcorresponds to setting Fσ2 = Fσ3 = 0 in Equation (2.2.20). In this case, we find thatFRTAσ =~ωσkBT 2vσ(2|vσ|L+∑σ2,σ312Γ−σ,σ2,σ3 + Γ+σ,σ2,σ3)−1, (2.2.22)which yields a lattice thermal conductivity for CNTs of the formκRTA =1NaAc∑σcσv2στ0σ=1NaAc∑σcσv2σ(1τbsσ+1τ3phσ)−1,(2.2.23)29with boundary (τbsσ )−1 and anharmonic three-phonon scattering rates (τ3phσ )−1 defined as1τbsσ=2|vσ|L,1τ3phσ=∑σ2,σ312Γ−σ,σ2,σ3 + Γ+σ,σ2,σ3 . (2.2.24)Applying the RTA, the form (2.2.23) yields a refined kinetic-gas-theory expression (seeEquation (1.3.7)), in that individual phonon heat capacities cσ, group velocities vσ = ∂kωσ,and mean free paths Λσ = vστ0σ are fully taken into account. Numerical calculations ofheat transfer rates of CNTs inferred from Equation (2.2.23) will be presented in Chapters 4and 5.As an aside, we note that anharmonic scattering rates (τ3phσ )−1 as per Equation (2.2.24)can be equivalently found in a classic paper by Maradudin and Fein [93], where (τ3phσ )−1 =−2 Im [Σσ(ωσ)] with Im [Σσ(ω)] denoting the imaginary part of the phonon self-energy dueto lowest-order anharmonic interactions. For a system in thermal equilibrium, −Im [Σσ(ω)]becomes experimentally observable as line broadening of phonon frequencies [93].2.2.4 Iterative solutionIn order to allow for a coupling of nonequilibrium phonons, an iterative solution approachaimed at Equation (2.2.20) was first proposed by Omini and Sparavigna [91]. Rearrangingterms in Equation (2.2.20), we can writeFσ =~ωσkBT 2τ0σvσ + τ0σ∑σ2,σ312[Fσ2 + Fσ3 ] Γ−σ,σ2,σ3+ [Fσ3 − Fσ2 ] Γ+σ,σ2,σ3 ,(2.2.25)with 1/τ0σ = 1/τbsσ + 1/τ3phσ defined as per Equation (2.2.24). This expression represents aself-consistency equation for the coefficients Fσ. A solution may be approached successivelyby performing an iteration as follows:F (i+1)σ = F(0)σ + τ0σ∑σ2,σ312[F (i)σ2 + F(i)σ3]Γ−σ,σ2,σ3 +[F (i)σ3 − F (i)σ2]Γ+σ,σ2,σ3 , (2.2.26)with starting condition F (0) = ~ωστ0σvσ/kBT 2. Terminating the iteration at the zeroth-orderstep is tantamount to the RTA, F(0)σ = FRTAσ , see Equation (2.2.22). For some predefinedtolerance value ε, the iteration (2.2.26) may be carried out until |F (i+1)σ − F (i)σ | < ε for all30phonon modes. Finally, the lattice thermal conductivity of CNTs is obtained asκiter =1NaAc∑σ~ωσvσn0σ(n0σ + 1)Fσ=1NaAc∑σcσv2στeffσ .(2.2.27)Here, in the second line, effective quasi-lifetimes are introduced, reading τ effσ = FσkBT2/vσ~ωσ,which allows one to recover an expression that is reminiscent of the RTA, see Equa-tion (2.2.23). However, since the τ effσ ’s reflect an interplay of out-of-equilibrium phonons,situations with τ effσ < 0 are generally possible such that the notion of a lifetime should notbe taken too seriously.The iterative method, as compared to the RTA, is computationally more expensive inso far as an evaluation of Equation (2.2.26) demands storing a large number of transitionamplitudes Γ± explicitly in memory. Despite the computational cost, going beyond theRTA bears physical significance. A key result of the PBTT due to Peierls [87] is thenecessity of umklapp scattering processes in giving rise to a nonvanishing thermal resistanceof realistic crystal lattices. In other words, normal scattering processes alone are incapableof dissipating a heat flux. Not distinguishing between normal and umklapp scatteringprocesses, the RTA, in contrast to beyond-RTA solution methods, fails to capture thisimportant result of the PBTT. While discussing numerical data of both κRTA and κiter inthe case of CNTs in Chapters 4 and 5, we will first take this deficiency of the RTA forgranted without further explanation. This point will then be revisited in Chapter 6, wherewe clarify, in the context of the PBTT, the necessity of distinguishing normal and umklappprocesses in mathematical terms.31Chapter 3Simulating phonon transport incarbon nanotubesIn this thesis, we conduct two different types of classical molecular dynamics (MD) simu-lations to study phonon transport in CNTs: nonequilibrium MD (NEMD) and equilibriumMD (EMD) simulations. Under the former simulation approach, no preconceived approxi-mations of phonon transport are needed and, in the high temperature limit, exact thermalconductivity predictions can be obtained at the system-level. On the downside, in NEMDsimulations, typically all information at the level of individual phonons is lost. If, on theother hand, one accepts a description of phonon transport at the level of an RTA (see Sec-tion 2.2.3), then a microscopic, phonon-resolved perspective of thermal transport can beachieved by means of EMD simulations.To give a background for our simulation data presented in Chapter 4, we outline herethe fundamental principles of our chosen NEMD and EMD methods and show how theseare applied specifically to CNTs. In particular, our NEMD approach will be detailed inSection 3.1. Here, we draw attention to some ambiguity in the evaluation of nonequilibriumtemperature profiles, that was only recently conclusively clarified in the literature by thework of Li et al. [94]. Focusing on phonon-level EMD simulations in Section 3.2, we point outhow the RTA description of phonon transport links to normal mode correlation functions.In this context, the approach of Thomas et al. [95] to investigate CNT heat transport viaEMD-derived atomic velocities will be discussed. Inspired by Ref. [95], we motivate a refinedEMD-based sampling technique of atomic velocities, that we use in Chapter 4 to offer newphonon-resolved insights into the system-size dependence of the CNT thermal conductivity.32Figure 3.1: NEMD simulation setup for calculating the thermal conductivity of CNTs.A temperature profile T (z1), . . . , T (zM ) is measured across M bins. Eight consecutivebins with coordinates z1, . . . , z8 and zM/2+1, . . . , zM/2+8 are designated as cold and hotregion, respectively. A quasi-ring-like geometry is enforced through periodic boundaryconditions (PBC) along the tube axis.3.1 The nonequilibrium MD methodThe nonequilibrium molecular dynamics (NEMD) method for the calculation of the phononconductivity κ closely resembles a real laboratory experiment, where one directly probesthe cause-effect relation between a temperature gradient and a heat flow as per Fourier’slaw (1.3.6) in order to measure κ. In NEMD simulations of stationary heat transport, atime-averaged temperature gradient is more readily obtained than a time-averaged heatflux, especially so in the case of many-body potentials [96–98]. Therefore, most state-of-the-art NEMD protocols [99–102], which prescribe various ways of driving a system intoa nonequilibrium state, revolve around imposing on a system a heat flux, that is knownexactly, and around measuring the resulting response in terms of a temperature gradient,that is subject to statistical noise. Focusing on a CNT as depicted in Figure 1.2, theNEMD-derived thermal conductivity is then calculated asκNEMD =jz〈∆T 〉t/L, (3.1.1)with 〈∆T 〉t = 〈Th− Tc〉t denoting a temperature difference between hot and cold reservoirsaveraged over a steady state simulation run time t. Here, any ambiguities in the definitionof the temperature gradient arising from possible temperature jumps at the reservoirs areexplicitly circumvented by setting −dT/dz = ∆T/L. This practice ensures a descriptionthat is consistent with ballistic transport limits (cf. Section 2.2.1), see also Li et al. [94]and further discussion below.33Figure 3.2: Simulated steady state temperature profiles of a (10,10) CNT for varyingtube lengths L.3.1.1 Simulation setup and temperature profilesIn this thesis, to simulate CNTs in a nonequilibrium steady state enabling a calculation ofκNEMD, we follow the velocity swapping protocol proposed by Mu¨ller-Plathe [100]. Here,two chunks of a CNT are designated as cold and hot region, and heat transport along thetube is induced by continuously swapping the velocity vector of the fastest atom in the coldregion vc with the velocity vector of the slowest atom in the hot region vh. In this way,by setting a constant velocity swap frequency, a steady state heat transport regime can beestablished. It is then straightforward to tally the amount of kinetic energy transferred overa simulation run time t, yielding a heat flux density for CNTs expressible asjz =∑swapsmC(v2c − v2h)/2Act. (3.1.2)A schematic illustration of our NEMD simulations is given in Figure 3.1. Cold and hotchunks of a CNT whose atoms are subject to velocity swapping are chosen with a width ofLth = 20 nm, which, for all conducted simulations, ensures a sufficiently broad velocity dis-tribution in the hot and cold slab so that |vc| > |vh| for all swaps. We apply equally spacedbins along the tube axis in order to monitor local tube temperatures T (z1), . . . , T (zM ). Thebin width is chosen as ∆z = 2.5 nm. The so-obtained resolution is sufficiently fine-grainedso that possible nonlinear features of T (z) in the vicinity of the thermal reservoirs can beresolved. Periodic boundary conditions along the tube axis serve to reduce boundary effects.Using Equations (3.1.1) and (3.1.2), our NEMD-derived thermal conductivity of CNTs iscalculated asκNEMD =∑swapsmC(v2h − v2c)/22Act〈∆T 〉t/L , (3.1.3)where an extra factor of 1/2 accounts for heat flow in two directions.While we defer the presentation of numerical results based on Equation (3.1.3) to Chap-34ter 4, a comment regarding some data of time-averaged temperature profiles is appropriate.In Figure 3.2, we show typical profiles obtained for CNTs as a result of velocity swapping.In particular in the case of CNTs of submicron lengths, it can be seen that temperatures offreely evolving tube parts, outside of hot and cold reservoirs, suggest |dT (z)/dz| ≤ ∆T/L.This observation is in accord with ballistic heat transport as described by the PBTT, seeSection 2.2.1, where flat temperature profiles in the limit of short tube lengths were at-tributed to weak phonon-phonon scattering. It is therefore evident that values of ∆T/Lshould be used in NEMD heat transport calculations. Contrary to this practice, by fitting|dT (z)/dz|, one may wrongly deduce that κNEMD(L→ 0) becomes unphysically large.3.2 Phonon properties from equilibrium MD simulationsAs pointed out in the context of the PBTT and the RTA, see Section 2.2.3, a phonon-resolveddescription of heat transport in terms of anharmonic lifetimes τ3phσ connects directly to thephonon self-energy, which in turn links to normal mode correlation functions in an equilib-rium state [103, 104]. Quantum-mechanically, such equilibrium correlation functions, whichcapture temperature dependent shifts and broadening of phonon frequencies in the presenceof anharmonic phonon-phonon interactions, are typically evaluated perturbatively with thehelp of diagrammatic techniques [93, 103]. Even though purely classical, a nonperturba-tive possibility of measuring normal mode correlation functions, enabling phonon-resolvedinsight into heat transport at the RTA level, is given by equilibrium MD (EMD) simula-tions [95, 105–109].To illustrate these concepts in the case of CNTs and to provide examples of some finite-temperature properties of phonons by the use of EMD simulations, we consider normalmode coordinates inferred from atomic velocities asPk,ν(t) =√mCN∑l,pe−ik lae˜∗p,k,ν · u˙l,p(t), (3.2.1)which is the inverse transformation of Equation (2.1.9), describing noninteracting travelingwaves under the harmonic lattice approximation. In the absence of lattice anharmonic-ity, and disregarding quantum effects, we have for the normal mode correlation function〈Pk,ν(t)P−k,ν(0)〉 = kBT cos(ωk,νt). Hence, the phonon spectrum of CNTs may be foundthrough a Fourier transform in time Fω asS(k, ω) = Fω[3N∑ν=1〈Pk,ν(t)P−k,ν(0)〉], (3.2.2)where the function S(k, ω), in the case of harmonic interactions, yields a sum of delta peaks35along the phonon dispersion branches ων(k). The Wiener-Khinchin theorem [85] states thatFω [〈Pσ(t)P ∗σ (0)〉] = |Fω [Pσ(t)] |2. Further, using the fact that phonon eigenvectors e˜p,k,νform a complete basis of C3N , it can be checked that Equation (3.2.2) is equivalent toS(k, ω) = mCN∑α∑p∣∣∣∣∣∑le−ik laFω[u˙αl,p(t)]∣∣∣∣∣2. (3.2.3)Without relying on the harmonic lattice approximation, this form allows one to producea finite-temperature dispersion relation of phonons by taking as input a time series ofatomic velocities u˙l,p(t) obtained from EMD simulations. As example, we perform EMDsimulations of a (10,10) CNT at varying temperatures to showcase the time-space Fouriertransform of atomic velocities as per Equation (3.2.3) in Figure 3.1. Here, it can be seenthat phonon frequency bands derived from a harmonic lattice calculation, as initially shownin Figure 2.1, provide an overall good approximation for phonon frequencies at finite tem-peratures. With increasing temperature, however, both a downward shift and a broadeningof phonon frequencies can be noticed, especially so in the case of high frequency bands lyingaround 50 THz.Thomas et al. [95] proposed to use EMD-derived phonon dispersion relations as perEquation (3.2.3) directly to gain phonon-resolved insight into CNT heat transport. Morespecifically, it was suggested to consider constant wave number slices of S(k, ω) and tofit the resulting frequency peaks by a sum of Lorentzian profiles, with the extracted full-width at half-maximum of individual peaks taken as inverse phonon lifetimes τMDσ . Theassumption of Lorentzian line profiles is generally justified in cases where the anharmonicphonon self-energy correction Σσ(ω) is small, |Σσ(ω)| < ωσ, and approximately constantaround unperturbed (harmonic) frequencies ωσ so that Σσ(ω) ≈ Σσ(ωσ) [93, 103, 104]. Fol-lowing Ref. [95], one would then proceed to calculate the sum∑σ cσv2στMDσ to obtain anRTA-level estimate of the CNT thermal conductivity. From a data processing perspec-tive, resolving S(k, ω) with wave number resolution ∆k = 2pi/Na and frequency resolution∆ω/2pi = 1/nsts requires a total EMD simulation run time t = nsts during which 3NNreal-valued coordinates have to be stored sequentially ns times. At the same time, to cap-ture the full spectrum of phonon frequencies, a sufficiently small sampling period ts hasto be chosen. We find that such storage requirements become quickly prohibitive when anaccurate sampling of low-frequency, long-wavelength phonons in CNTs is desired.In this thesis, we run EMD simulations during which we perform on-the-fly calculationsof only a subset of time-evolved normal mode coordinates. In doing so, we are able toaddress the aforementioned storage bottleneck and achieve much improved statistics onlow-frequency, long-wavelength modes, at the cost of filtering out many high-frequency36Figure 3.1: Finite-temperature phonon dispersion relations of a (10,10) CNT as ob-tained from EMD simulations at temperatures 50 K, 300 K, 550 K, and 800 K. Theintensity of S(ω, k), see Equation (3.2.3), is plotted on a logarithmic scale in arbitraryunits.37modes that we identify to contribute negligibly in the summation∑σ cσv2στσ. In particular,we use phonon eigenvectors obtained from a harmonic lattice calculation to compute Pσ(t)for selected modes according to Equation (3.2.1). An attempt is then made to identify asingle peak in the corresponding spectral density of Pσ(t) and to fit a single Lorentzianprofile,|Fω [Pσ(t)] |2 = Cσ(ω − ω˜σ)2 + (2 τMDσ )−2, (3.2.4)with fitting parameters ω˜σ, 1/τMDσ and Cσ. We will discuss numerical results based onEquation (3.2.4) in Chapter 4 in an effort to shed light on system-size effects of the CNTthermal conductivity.38Chapter 4Comparison of calculation andsimulation based approaches:Effects of system size andtemperatureIn the previous two chapters, we described the two main-stream computational approachesfor predicting phononic heat transfer rates of pristine carbon nanotubes: the PBTT and MD.In each respective framework, we introduced several types of conductivity calculations thatare subject to different approximations. In the PBTT, κbal is the ballistic conductivity whichis exclusively due to harmonic atom-atom interactions. Accounting for lattice anharmonicitythrough lowest-order perturbation theory, κRTA is the zeroth-order solution for the linearizedphonon Boltzmann equation based on which a more accurate solution κiter can be foundvia iteration. The iteration serves to distinguish between normal and umklapp scatteringprocesses. In the framework of MD, direct simulations of heat conducting CNTs leadto κNEMD, whereas individual heat carrying phonons can be analyzed by means of EMDsimulations through the lens of an RTA.In the present chapter, these types of calculations are brought together and numericalresults of phonon conductivities for a selection of armchair CNTs are presented. By carefullycomparing PBTT and MD-based approaches, we aim to clarify some of the controversiesthat exist around heat transport in individual CNTs. In particular, Section 4.1 focuses onthe system-size dependence of the phonon conductivity at room temperature. Here, specialattention is devoted to the role of acoustic phonons in CNT heat transport. Subsequently,in Section 4.2, we concentrate on the temperature dependence of phonon transport. Themain points to take away from our numerical calculations will be summarized in Section 4.3.394.1 System-size dependence of the thermal conductivityfrom the PBTT and NEMDIn Figure 4.1, we show the length-dependent CNT thermal conductivity κ(L) at room tem-perature as obtained by PBTT calculations and by NEMD simulations. A direct comparisonof the two different computational approaches is presented in Figure 4.1(a), where the left,middle and right panels refer to tube diameters of 5.43 A˚, 9.49 A˚ and 13.56 A˚, respectively.To show trends with respect to the tube diameter more clearly, the solid curves of Fig-ure 4.1(a) as derived from the PBTT and the NEMD method are replotted on log-linearscale in Figures 4.1(b) and 4.1(c), respectively.For tubes of varying diameters, both computational approaches yield similar trends atthe submicron scale, as can be seen in Figure 4.1(a). At this length scale, PBTT-derivedthermal conductivity values κiter are only weakly affected by intrinsic three-phonon scat-tering processes, which can be noted from the fact that κiter(L) approaches a linear scalingκ ∼ L as L → 0, which is, by definition, the ballistic limiting behavior of κiter(L) in theabsence of lattice anharmonicity. For short tubes, irrespectively of the tube diameter, wefurther observe that NEMD-derived transport predictions κNEMD surpass correspondingnumerical estimates from the PBTT. As detailed in Section 2.2.1, ballistic upper limitsto phonon transport imply κbal(L) = γL, where the slope γ depends on per-phonon heatcapacities cσ and group velocities vσ through γ ∝∑σ cσ|vσ|. We indicate by dashed anddashed-dotted lines in Figure 4.1(a) the ballistic upper limits κbal(L) to phonon transportthat result from the quantum and classical expressions for the phonon heat capacities cσ, re-spectively. Here, it becomes apparent that the overestimation of the true (quantum) phononheat capacity in MD favors an overprediction of heat transfer rates in the NEMD method,provided that the CNTs are sufficiently short so that transport remains predominantly bal-listic. Sensible upper bounds to thermal transport in the ballistic limit are confirmed byfact that NEMD-derived predictions κNEMD fall below the dash-dotted line in Figure 4.1(a).In the transition away from the ballistic transport regime, for all tube diameters con-sidered, a cross-over length scale can be observed in Figure 4.1(a) beyond which we findκNEMD < κiter. Thus, for L > 1µm, an overestimation of NEMD-derived classical heattransfer rates is no longer guaranteed. As tube lengths go beyond the micrometer scale,intrinsic anharmonic thermal resistance becomes increasingly important and its distincttreatment in the PBTT and NEMD method becomes evident. PBTT calculations do notincorporate beyond-lowest-order phonon-phonon scattering processes, such that absolutevalues of κiter should be regarded as upper bound to the true CNT thermal conductivity.NEMD results, on the other hand, are obtained nonperturbatively, but phonon occupationnumbers of high-frequency modes are largely overestimated in MD which brings about anoverestimation of phonon-phonon scattering rates. Therefore, at the micrometer scale and40Figure 4.1: System-size dependence of the CNT thermal conductivity κ at T =300 K as predicted by PBTT calculations and NEMD simulations. (a) Comparisonof computational approaches with CNT diameter increasing from the left to the rightplot. Straight lines represent ballistic upper limits to κ, with dashed and dashed-dotted lines accounting for quantum and classical phonon heat capacities, respectively.PBTT calculations shown in (a) are performed with a 1D wave number grid comprisingN = 400 k-points. (b) PBTT results in log-linear scale obtained by imposing varyingwave number grids, with reference data taken from Lindsay et al. [48]. (c) NEMDresults in log-linear scale with reference data from Cao and Qu [49] and from Sa¨a¨skilahtiet al. [110].41beyond, we expect the true thermal conductivity of CNTs to be sandwiched between thePBTT- and NEMD-derived curves. Such a statement, however, requires caution since ourvalues of κNEMD(L) are not guaranteed to yield a lower bound to the true conductivity ina strict sense. In the NEMD method, an overestimation of phonon-phonon scattering ratesmight be masked by an overestimation of phonon heat capacities, especially so in the caseof our classical simulations at room temperature, far below the CNT Debye temperature(≈ 2200 K). We note that some authors [108, 111, 112] employed classical equilibrium oc-cupation numbers, n0,clσ = kBT/~ωσ, in the PBTT formalism. In this way, an attempt canbe made to shed light on the role of phonon statistics in giving rise to anharmonic ther-mal resistance. While this approach might be justifiable at relatively high temperaturescompared to a material’s Debye temperature, it lacks firm theoretical footing in the PBTTformalism as pointed out by Feng et al. [90]. Given the extremely high Debye temperatureof CNTs, we refrain here from such an approach.Moving towards macroscopically long CNTs, both the PBTT and NEMD method facedifferent numerical challenges. As indicated by the dashed and dashed-dotted lines inFigure 4.1(b), we find that for tube lengths L > 10µm PBTT-derived values κiter becomeincreasingly sensitive to the employed wave number grid resolution. In fact, as we willdiscuss in Section 4.1.2, our numerical calculations of κiter in the limit L → ∞ cannot beconverged to an unambiguous limiting value. As such, in the framework of the PBTT,evidence is found for the absence of a fully diffusive heat transport regime in the case ofpristine CNTs. In NEMD simulations, the maximally accessible system size is restrictedby the atom count and by an increase of simulation time required to drive longer CNTsinto a nonequilibrium steady state. As can be seen in Figures 4.1(a) and 4.1(c), in the caseof pristine CNTs, the steady increase of NEMD-derived conductivity values κNEMD in therange up to L = 10µm precludes definite statements as to the existence of a fully diffusiveheat transport regime. Looking at Figure 4.1(c), a power-law divergence κNEMD ∼ Lα with0 < α < 1 in the case of the (4,4) CNT and a weaker, possibly logarithmic, divergence inthe case of the (7,7) and (10,10) tubes might be deduced. Equally likely, NEMD resultscould be described by a convergent expression of the form κ(L) =∑σ κ∞σ (1 + Λσ/L)−1,assuming a phonon-resolved description of the thermal conductivity where a few selectedmodes exhibit finite but extremely long intrinsic mean free paths Λσ.Interestingly, regarding the tube diameter dependence of the thermal conductivity, op-posing trends can be noted from Figures 4.1(b) and 4.1(c). A decreasing trend of κiterwith increasing tube diameter D in the range 5.43 A˚ ≤ D ≤ 13.56 A˚ is found by meansof PBTT calculations. While the diameter-dependent trend of κiter appears to be ro-bust with respect to the tube length in Figure 4.1(b), it can be speculated, based onthe varying slopes of κNEMD(L) in the micrometer regime, that the opposing diameter-42dependent trend in Figure 4.1(c) eventually reverses in larger NEMD simulations. In fact,in the limit of short tubes, focusing on ballistic upper bounds, κbal(L) = γL, as describedabove, we find that classical phonon heat capacities give rise to slightly larger values ofγ in larger diameter tubes. For example, considering the (4,4) and (10,10) tubes, we findγ(4,4)/γ(10,10) ≈ 0.91 if the classical phonon heat capacities are assumed, whereas the ratioincreases to γ(4,4)/γ(10,10) ≈ 0.99 if the true (quantum) heat capacities are adopted. Tosome extent, these numbers might help to explain the contradicting diameter trends forthe finite-length CNTs as obtained by the quantum and classical computational approachesshown in Figure 4.1.4.1.1 Comparison with earlier modeling efforts of κ(L,D)PBTT calculationsPreviously, similar PBTT calculations of κiter under the inclusion of the full spectrum ofanharmonic three-phonon scattering processes were performed by Lindsay et al. [48, 56] and,more recently, by Yue et al. [47]. The former authors investigated the length, κiter(L) [56],and diameter dependence, κiter(D) [48], of the CNT thermal conductivity based on forceconstants derived from different parametrizations of the Tersoff potential. Adopting theoptimized parameter set [72], diameter-dependent values of κiter for L = 3µm given inRef. [48] are reasonably well reproduced by our calculations, as shown by the inset ofFigure 4.1(b). Considering phonon transport at macroscopic length scales, it was notedin Ref. [56] that a PBTT calculation approach in conjunction with lowest-order phonon-phonon scattering yields a tube-length divergence of κiter(L→∞), which can be attributedto low-frequency FA modes of CNTs. As we will discuss below in Section 4.1.2, at the levelof individual phonons, this observation falls into line with our numerical findings, but FAmodes are in fact found to be not the only contributors to system-size convergence issues.Yue et al. [47] studied the diameter dependence of CNT thermal conductivity in theinfinite-length limit, using force constants from ab initio calculations. Here, the authorsreported a decreasing trend of κiter(D) with increasing D which conforms to our results inthe range 5.43 A˚ ≤ D ≤ 13.56 A˚. However, it was also noted in Ref. [47] that values ofκiter(L → ∞) can be converged with respect to the employed wave number grid withoutproblems, which is in disagreement with both Lindsay et al. [56] and our findings, castingsome doubts on the numerical calculations performed in Ref. [47].NEMD simulationsUsing the optimized Tersoff potential [72], absolute values of tube-length-dependent CNTthermal conductivities were derived earlier from NEMD simulations by Cao and Qu [49]43and by Sa¨a¨skilahti et al. [110]. The former authors considered tubes of varying diametersin the range L ≤ 2.4µm, while the latter focused on a (10,10) CNT with L ≤ 4µm. Asindicated by the inset of Figure 4.1(c), our NEMD results obtained for a (10,10) CNTcompare satisfactorily with the data reported in these studies. Moreover, our NEMD-derived diameter dependence agrees with the simulation data of Cao and Qu [49] to theextent that an increase of κNEMD with increasing diameter is observed. These authors alsoarrived at the conclusion that κNEMD values diverge with a power-law dependence away fromthe ballistic regime, κNEMD ∼ Lα, where they note a tube-diameter-independent exponentα ≈ 0.21. A universal, diameter-independent scaling behavior of κNEMD(L), however, isnot confirmed by our NEMD simulations in the micrometer range. Instead, we can see inFigure 4.1(c) that the smaller diameter (4,4) CNT exhibits a more pronounced increase ofκNEMD with L for L > 1µm as compared to the larger diameter tubes of chirality (7,7) and(10,10).Given the prohibitive computational demands of simulating CNTs beyond the microm-eter scale, Sa¨a¨skilahti et al. [110] extrapolated the tube-length dependence of κNEMD basedon a spectrally resolved phonon mean-free path Λ(ω) obtained from NEMD simulations.Here, the authors identified the importance of low-frequency phonons whose dissipationrates in the limit ω → 0 were found to determine whether κNEMD saturates or diverges withL. On a similar note, in Sections 4.1.3 and 4.1.4, we study the attenuation of low-frequencyacoustic phonons in the frameworks of the PBTT and EMD simulations, respectively.4.1.2 PBTT calculations of κ(L→∞)In principle, within PBTT calculations, it is straightforward to address the infinite-sizelimit of CNT thermal transport, setting (τbsσ )−1 = 2|vσ|/L ≡ 0 for all modes. Thermalconductivity predictions resulting from such calculations are presented in Figure 4.2. Here,for varying 1D wave number grids with N distinct k-points, we show the branch-by-branchcontributions to the thermal conductivity following from the iterative solution approach,κiterν ∝ N−1∑k ck,νv2k,ντeffk,ν , and from the RTA, κRTAν ∝ N−1∑k ck,νv2k,ντ3phk,ν .Based on Figure 4.2, two important observations can be made regarding κiter(L→∞):(i) For all tube diameters considered, total thermal conductivity predictions depend heavilyon the chosen wave number grid which is mostly due to FA and TW phonon branches, whoseconductivity contributions increase drastically with finer k-grids. Even more critical, somewave number grids can be found that do not give rise to a numerically stable iteration at all.(ii) A notable variation of κiter with respect to the CNT diameter is mostly due to differingFA mode behavior. For finer and finer k-grids, the smaller the tube diameter, the fasterthe rate of divergence of FA mode contributions. These numerical results are indicative ofa tube-length divergence of the CNT thermal conductivity at macroscopic length scales.44Figure 4.2: PBTT calculation results for CNTs at T = 300 K obtained by setting(τbsσ )−1 = 2|vσ|/L ≡ 0, assuming infinitely long tubes. The horizontal axis signifies howmany equidistant wave number grid points are included in the calculation. The totalsum of per-phonon conductivities (black squares) is decomposed into contributions fromacoustic longitudinal (LA), twisting (TW), flexural (FA), and optical phonon branches.In the top and bottom row, PBTT calculation results are shown that stem from aniterative solution approach to the phonon transport equation and from a relaxationtime approximation (RTA), respectively. As indicated by dashed vertical lines, theiterative method does not converge for certain wave number grids.In comparing iterative with RTA-based calculations in Figure 4.2, an overall increase ofconductivity predictions is found by going beyond the RTA. This behaves according to thecommonly adopted notion in which the RTA is assumed to underpredict the thermal con-ductivity, given its lack of distinction between normal and umklapp scattering processes [50].Importantly, it can also be noted from the bottom panels of Figure 4.2 that an apparentgrid size convergence of κRTA(L → ∞) may in fact result from a spurious nonmonotonicbehavior of FA modes. This appears to be especially so in the case of the (10,10) CNT (seethe bottom right panel of Figure 4.2), where the TW mode contributions continue to growfor finer and finer k-grids, despite an apparent saturation of the total thermal conductivityκRTA.45Our mode-resolved PBTT calculations make it clear that acoustic phonons play a gov-erning role in CNT heat transport. These modes are predicted to clearly dominate phonontransport in the long-tube limit. To develop a deeper understanding of size effects in CNTphonon transport, in the following, we discuss in more detail the conductivity contributionsof individual acoustic phonon branches.4.1.3 PBTT calculations of acoustic phonon lifetimesIf we convert the discrete summation over N distinct wave numbers to an integration,N−1∑k →∫dk, then the PBTT-derived thermal conductivity of CNTs, both at the RTAlevel and beyond, can be written asκ =3N∑ν=1∫ pi/a−pi/adk2piκν(k), (4.1.1)with wave number-dependent per-phonon branch contributions readingκν(k) =1Accν(k)vν(k)2τν(k). (4.1.2)Here, the functions cν(k) and vν(k) are always well-behaved and upper-bounded throughoutthe whole phonon spectrum. In contrast, if we assume that the dependence of κ on the tubelength L enters through an upper limit imposed on phonon lifetimes as (τσ)−1 = (τ ισ)−1 +2|vσ|/L,1 where ι denotes “3ph” or “eff”, reflecting an RTA or beyond-RTA treatmentof phonon-phonon scatterings in infinitely extended tubes, respectively, then it becomesclear that a source of divergence in Equation (4.1.1) arises from the small k behavior of(τ ισ)−1. Namely, by applying a crude approximation of phonon-phonon scattering rates,noting (τ ισ)−1 . ωσ, a tube-length divergence of κ(L → ∞) can be linked to the gaplessacoustic phonon branches, ν ∈ {LA,TW,FA}, for which ων(|k| → 0) = 0. Previous theo-retical studies [53, 113, 114] of CNT heat transport invoked the Klemens approximation oflow-frequency anharmonic three-phonon scattering rates, (τ3ph)−1 ∼ ω2 [115, 116], to arguethat the thermal conductivity should be divergent with tube length. The Klemens approx-imation, however, fails to capture the quadratic phonon dispersion of FA modes, ωFA ∼ k2,and is therefore qualitatively wrong in the case of CNTs, as was later pointed out by Lind-say et al. [56]. Generally, if we assume a power-law dependence of acoustic phonon-phonon1The summation of inverse lifetimes is known as Matthiessen’s rule [50]. At the RTA level, a formal justi-fication for (τσ)−1 = (τ3phσ )−1 + 2|vσ|/L was presented in Section 2.2.3. Going beyond the RTA, it has to benoted that, for finite L, iteratively determined values of κiter(L) as per Equations (2.2.26) and (2.2.27) aregenerally not equal to Equations (4.1.1) and (4.1.2) when the simple substitution (τσ)−1 = (τ effσ )−1 + 2|vσ|/Lis performed. Still, the Matthiessen approach, by definition, converges toward κiter(L) as per Equa-tion (2.2.26) as L becomes large.46scattering rates, (τ ισ)−1 ∼ ka for |k| → 0, we find for those phonon branches that disperselinearly, namely, the LA and TW branches, that conductivity contributions in the limitL→∞ scale asκLA/TW(L) ∼∫0τ(k) dk ∼∫0dkka + 1/L∼{L1−1/a if a 6= 1,lnL if a = 1.(4.1.3)Similarly, given the quadratic dispersion of FA branches, corresponding branch contributionsin the long-tube limit behave asκFA(L) ∼∫0v(k)2τ(k) dk ∼∫0k2 dkka + k/L∼{L1−2/(a−1) if a 6= 3,lnL if a = 3.(4.1.4)This suggests to perform PBTT calculations for infinitely extended tubes, setting (τbsσ )−1 =2|vσ|/L ≡ 0, so that the small k behavior of acoustic phonon lifetimes, resulting exclusivelyfrom three-phonon scattering processes, can be examined.In Figure 4.3, we show numerical results of τ3phν (k) and τ effν (k) for varying tube di-ameters and ν ∈ {LA,TW,FA}. Here, we indicate by dashed lines the critical power-lawdependencies ∼ k−1 and ∼ k−3 inducing a tube-length divergence of κν(L) as per Equa-tions (4.1.3) and (4.1.4), respectively. As can be seen from the left plots in Figure 4.3,at the RTA level, no indication for a length divergence of κRTAν (L) can be deduced fromthe scaling behavior of τ3phν (k) for ν ∈ {LA,TW,FA}. Predictions of τ3phLA (k) suggest in-creasingly strong anharmonic dissipation rates of LA modes in the long-wavelength limitk → 0+. Hence, the numerical integration κRTALA ∼∫0 τLA(k) dk can be straightforwardlycarried out. Although not necessarily indicative of a tube length divergence of κRTATW (L),the scaling behavior of τ3phTW (k) suggests a significant conductivity contribution stemmingfrom long-wavelength TW modes. Such scaling behavior at the very least aggravates wavenumber grid convergence issues, when accurate estimates of κRTATW (L→∞) are desired. FAmode lifetimes τ3phFA (k) seem to plateau in the long-wavelength limit k → 0+, but also dis-play some spurious discontinuities in the small k region. Since FA mode frequencies vanishquadratically for k → 0+, corresponding energy scales become so small that the numericalsearch for energy conserving three-phonon scattering partners may fail, which in turn maycause artifacts in the predictions of τ3phFA (k → 0+). These artifacts hint at the source ofnonmonotonic behavior of κRTAFA (L→∞) as previously observed in Figure 4.2.Adopting a more accurate, beyond-RTA description of phonon transport, effective life-times τ effν (k) shown in the right-hand side of Figure 4.3 point very clearly towards a lengthdivergence of κiter(L) in the case of all tube diameters considered. Going from τ3phν (k)to τ effν (k), the scaling behavior of FA modes in particular undergoes a qualitative change,47Figure 4.3: Lifetimes of acoustic phonons of CNTs at T = 300 K as obtained inthe PBTT framework. The left and right columns show lifetimes that enter into thecalculation of κRTA and κiter, respectively. Dashed lines represent critical power lawdependencies τ ∼ k−a, which go along with a tube-length divergence of κ.explaining a divergence of κiterFA (L → ∞). This observation is in agreement with the ear-lier findings of Lindsay et al. [56], according to which nonequilibrium occupation numbersof low-frequency FA modes resulting from three-phonon interactions are responsible for adivergence of κiter(L → ∞). FA modes, however, are not the sole reason for the lack ofconvergence in the long-tube limit, as found in our PBTT calculations. Effective lifetimesof TW modes τ effTW(k) clearly approach a critical scaling behavior as well (here ∼ k−1), sothat an additional source of divergence can be attributed to κiterTW(L → ∞). As a conse-quence, when numerical grid integrations are performed, one faces ever-increasing valuesof κiter(L → ∞) as the resolution of the wave number grid increases. To the best of ourknowledge, the prominent role of low-frequency TW modes in CNT heat transport has notyet been explicitly pointed out in the literature. Unlike τ effFA (k) and τeffTW(k), effective life-48times of LA modes τ effLA(k) appear to be upper bounded as k → 0+, which suggests that theconductivity contributions stemming from LA modes remain finite in the long-tube limit,κiterLA (L→∞) <∞.4.1.4 EMD simulations of acoustic phonon lifetimesAs explained in Chapter 3 EMD simulations enable nonperturbative calculations of theclassical limit of phonon lifetimes. Given the presence of higher-order phonon-phonon in-teractions, it is interesting to consider EMD-derived trends of acoustic phonon lifetimesτMDν (k → 0+) and to test whether they agree with perturbative predictions of τ3phν (k) aspreviously shown in Figure 4.3.EMD simulations of finite-length CNTs cannot resolve arbitrarily small wave numbersk. Under periodic boundary conditions along the tube axis, only a discrete set of wavenumbers is contained in a simulation cell of N translational unit cells, cf. Equation (2.1.5).To determine acoustic phonon lifetimes by means of EMD in a regime of long wavelengths,ideally as long as possible, we simulate CNTs of chirality (3,3), (4,4), (7,7) as well as(10,10), and choose a simulation cell that comprises N = 400 unit cells. This implies amaximally resolvable wavelength of λmax = 400 a ≈ 100 nm. We further restrict our analysisto the eight smallest wave number modes so that k is an integer multiple of 2pi/λmax inthe range 2pi/λmax ≤ k ≤ 16pi/λmax. In this way, by performing on-the-fly calculationsof instantaneous normal mode coordinates Pk,ν(t) as detailed in Section 3.2, we achievean acceptable tradeoff between system size and simulation runtime in our calculations ofτMDν (k).In Figure 4.4, we show a selection of spectral densities of acoustic normal mode co-ordinates |Fω[Pk,ν(t)]|2 for ν ∈ {LA,TW,FA} as obtained from our EMD simulations ofa (10,10) CNT, that are to be used to infer τMDν (k → 0+). For tubes of varying diame-ters, the complete set of spectral densities extracted from our EMD simulations is givenin Appendix B. Under the assumption that the effects of anharmonic phonon-phonon in-teractions can be captured by just two mode-specific parameters, a frequency shift ω˜σ anda line broadening factor ∝ 1/τMDσ , we expect each mode to exhibit a single frequencypeak, approximately described by a Lorentzian profile with its center near the harmonicfrequency prediction ωσ. As Figure 4.4 clearly demonstrates, however, a single peak withLorentzian type broadening may not necessarily provide an accurate description of phononsin the low-frequency regime. In particular in the case of low-frequency FA modes, we ob-serve significant deviations from single-peak behavior in a regime where the correspondingharmonic mode frequencies of FA modes fall below 0.3 THz. At these very low FA modefrequencies, the approximations of lowest-order anharmonic perturbation theory, appliedwithin the PBTT, may become questionable and a more refined theory may be needed to49Figure 4.4: Spectral densities of acoustic phonons at T = 300 K as obtained fromEMD simulations of a (10,10) CNT comprising N = 400 unit cells. Columns arearranged in order of increasing wave number (from left to right). Frequency axesare rescaled to show deviations from harmonic frequency predictions ωσ. To estimatephonon relaxation times τMDσ , a Lorentzian fit (dashed-dotted lines) is applied in thecase of LA (top row) and TW modes (middle row). Qualitatively similar spectraldensities are found for the two degenerate FA branches.account for double peaks and other complexities in our data of |Fω[Pk,FA(t)]|2. In light ofthis, we proceed by excluding FA modes from our EMD phonon lifetime calculations. Al-though some deviations from single-peak behavior manifest also in the case of low-frequencyLA and TW modes, see for example the TW mode at k = 16pi/λmax in Figure 4.4, thesedeviations are typically far less pronounced. Hence, a Lorentzian fit can still be appliedto our data of |Fω[Pk,LA(t)]|2 and |Fω[Pk,TW(t)]|2, enabling an estimation of τMDLA (k) andτMDTW (k), respectively. Acoustic phonon lifetimes of LA and TW modes as obtained by fittingEMD-derived spectral densities are summarized in Figure 4.5.It can be observed in Figure 4.5 that LA modes attenuate faster in tubes of chirality (3,3)and (4,4) as compared to the larger diameter (7,7) and (10,10) tubes. Such a qualitative50Figure 4.5: Lifetimes of longitudinal (LA) and twisting (TW) modes of CNTs atT = 300 K as obtained from EMD simulations. Plots are arranged in order of in-creasing diameter (from left to right). Underlying MD simulations comprise N = 400translational unit cells. The eight smallest resolvable wave numbers are shown thatare commensurate with the simulation cell. Insets display corresponding phonon fre-quencies, both from a harmonic lattice calculation (drawn by lines) and from fittingthe spectral density of normal mode coordinates (drawn by plus markers).trend with respect to the tube diameter can also be noticed in PBTT-derived predictionsof τ3phLA in Figure 4.3 where LA lifetimes are found to reach larger values in larger diametertubes. Fit lines τMDLA ∼ k−a with significantly smaller exponent a in smaller diameter tubesmight be ascribed to an onset of a saturation in the limit k → 0+ which, according toour perturbatively determined values τ3phLA , sets in for larger values of k in smaller diametertubes. In the case of the TW branch, for all considered tube diameters, a monotonic increaseof phonon lifetimes τMDTW down to the smallest resolvable wave number is detected. Fit linesτMDTW ∼ k−a indicate that exponents lie in the range from 0.80 to 1.16, which conforms withan earlier EMD result, predicting τMDTW ∼ k−1.1, based on simulations of shorter (10,10)CNTs [117]. Overall, both the lack of a distinctive diameter dependence of τMDTW and ascaling exponent a ≈ 1 seem to support our perturbative predictions of τ3phTW . As such,EMD simulations clearly corroborate the predominant role of long-wavelength TW modesin CNT heat transport. Within the scope of an RTA, but going beyond the limitations oflowest-order anharmonic perturbation theory, some evidence is found that TW modes maygive rise to a tube-length divergence of the CNT thermal conductivity.51Figure 4.1: Top panel: Thermal conductivity κ vs temperature T for a (7,7) CNT ofvarying tube length L. Empty and filled markers display PBTT calculation (κiter) andNEMD simulation results (κNEMD), respectively. PBTT calculations are performedwith a 1D wave number grid comprising N = 400 k-points. Bottom panel: Same dataas above, but with κ rescaled by T . Rescaling κ serves to show possible deviationsfrom κ ∼ T−1 in the high temperature regime.4.2 Temperature dependence of the thermal conductivityfrom the PBTT and NEMDThe temperature dependence of the phonon conductivity κ(T ) of CNTs as obtained fromPBTT calculations and NEMD simulations is shown in Figure 4.1. As can be noted fromboth quantum mechanical (PBTT) and classical (NEMD) approaches, the observed tem-perature dependence of CNT heat transport strongly depends on the tube length L. Inthe framework of the PBTT, temperature dependent ballistic upper limits to heat trans-fer rates become apparent for T < 100 K. In this temperature regime, conductivity dataκiter(T ) for varying tube lengths can be collapsed onto a single curve G(T ) ∝ κ(T )/L, whichis physically equivalent to a tube length-independent thermal conductance G of CNTs in apredominantly ballistic transport regime at low temperatures. For micrometer-long CNTs52with 1µm ≤ L ≤ 10µm, κiter(T ) peaks near T ≈ 150 K, which agrees well with the PBTTcalculations of Mingo and Broido [53] and Lindsay and Broido [60]. Surpassing the peaktemperature, as the relative importance of anharmonic phonon-phonon scattering processesincreases, the scaling behavior of κiter(T ) converges toward κ ∼ T−1, which manifests in thebottom panel of Figure 4.1 as an asymptotic flattening of κiter(T )T in the high temperaturelimit. This is the expected limiting behavior of κiter(T ) given the limitations of lowest-orderanharmonic perturbation theory, which can be noticed from the fact that three-phononscattering rates, see Equation (2.2.21), depend linearly on equilibrium phonon occupationnumbers which scale as n0(T ) ≈ kBT/~ω for ~ω kBT . As can also be seen in the bottompanel of Figure 4.1, a regime exists above the peak temperature of the conductivity whereCNTs with L & 5µm exhibit a steeper than ∼ T−1 dependence. Such a steep temperaturedependence reflects a freeze-out of umklapp scattering processes. Given the energy andmomentum selection rules of three-phonon scattering processes, see Equation (2.2.15), itis easily checked that a scattering event of umklapp type has to involve at least one high-frequency phonon. Since quantum occupation numbers of high-frequency phonons vanishexponentially at small temperatures, n0(T ) ≈ exp(−~ω/kBT ) for kBT ~ω, one gener-ally assumes that the corresponding anharmonic thermal resistance vanishes in the samemanner [50, 51]. According to our PBTT calculations, such manifestations of umklappscattering, leading to a steeper than ∼ T−1 scaling, critically depend on the CNT lengthL, to the point of being completely suppressed in submicrometer long tubes.Our classical NEMD results in Figure 4.1 clearly demonstrate a tube length-dependentonset of anharmonic thermal resistance. As can be appreciated from the flattening ofκNEMD(T )T in the bottom panel of Figure 4.1, relatively expensive simulations of micro-meter-long CNTs with more than ≈ 105 atoms are required to enter a regime where latticeanharmonicity becomes the limiting factor of thermal transport. Despite the inclusion ofhigher-order phonon-phonon interactions in the NEMD method, for micrometer-long CNTsand T ≤ 800 K, no convincing evidence for a rapid fall-off ∼ T−2 at high temperature can befound, which would indicate a break-down of lowest-order anharmonic perturbation theoryin the PBTT.In light of the rather different approximation that enter into the quantum and clas-sical approaches in Figure 4.1, absolute values of κiter(T ) and κNEMD(T ) are found tocompare surprisingly well in a high-temperature regime where the effects of anharmonicphonon-phonon scattering set in. Unfortunately, although classical values of κNEMD(T ) areexpected to become more and more accurate at high temperature, we are not able to ex-tend our NEMD simulations of CNTs beyond 800 K. Running MD simulations at such hightemperatures, we find that our chosen empirical interaction potential produces spurioussignatures of bond dissociation.534.3 ConclusionsFocusing on a selection of armchair tubes and adopting an empirical atomic interactionpotential, we have pursued state-of-the-art modeling approaches of phonon transport tostudy the length, diameter, and temperature dependence of the CNT thermal conductiv-ity. With regard to tube length effects, both quantum mechanical PBTT calculations andclassical NEMD simulations, independently of another, confirm a relatively weak influenceof intrinsic anharmonic resistance in submicrometer long CNTs. In this regime, at sub-micrometer length scales, an overestimation of heat transfer in the NEMD method can bewell understood by the notion of ballistic phonon transport in the framework of the PBTT.Drawing definite conclusions about phonon-mediated heat transport in the long-tube limitproves to be challenging by means of direct NEMD simulations, which are limited to atomcounts on the order of 106. NEMD-derived conductivity values for finite-length CNTs maybe extrapolated beyond the micrometer length scale, whereby both length-divergent andconvergent fit functions are conceivable. In the PBTT, the long-tube behavior of phonontransport is governed by the dissipation of low-frequency acoustic phonons. By numericallyevaluating the full spectrum of three-phonon scattering processes, a length divergence ofthe phonon conductivity is encountered which can be traced to nonintegrable FA and TWmode lifetimes in the long-wavelength limit. These findings are partly supported by EMDsimulations, which corroborate critically low dissipation rates of TW modes in particular.Qualitatively different diameter trends as obtained from PBTT calculations [48] andNEMD simulations [49], that can be found in the literature for the same force field parametriza-tions of CNTs, are reaffirmed in our study and remain an open theoretical issue. Using thenotion of ballistic phonon transport, some evidence is found that classical phonon statisticsin MD favor larger conductivity values in larger diameter tubes, which is less so when thecorrect quantum statistics in the PBTT are adopted.In individual heat conducting CNTs, the relative importance of anharmonic thermalresistance can be tuned by either varying the tube length or the lattice temperature, whichtranspires from both our PBTT calculation results and our NEMD simulations. By extend-ing our NEMD simulations to micrometer length scales and to temperatures of 800 K, aheat transport regime is entered where anharmonic phonon-phonon interactions become aspronounced as possible, factoring in the computational limits of our all-atom simulations.In this regime, no convincing evidence for the importance of higher-order phonon-phononinteractions can be found. This gives credence to our PBTT calculations which incor-porate only three-phonon interactions. It has to be conceded, however, that conclusionsabout higher-order phonon-phonon interactions may at least in part depend on the chosenempirical interaction potential.54Chapter 5Low-frequency three-phononscattering rates and the role oftensile lattice strain5.1 IntroductionA peculiarity of both one and two-dimensional materials is the existence of low-frequencyflexural bending (FA) modes. One special feature of these out-of-line or out-of-plane vi-brations is that associated mode frequencies are hardened in the presence of tensile lat-tice strain. By taking into account the interatomic forces of prestretched low-dimensionalcrystals, a strain hardening of FA modes can be directly demonstrated by means of lattice-dynamical phonon calculations, see, for example, Refs. [118, 119] for some related results onCNTs and graphene. In this context, a notable observation concerning phonon-mediatedheat transport in graphene was made by Bonini et al. [82] who focused on the role of FAmodes in the framework of the PBTT. Namely, it was found that three-phonon scatter-ing processes with the involvement of FA modes become systematically weakened underincreasing tensile lattice strain. Specifically, by examining analytically the implications ofanharmonic perturbation theory in the limit of low phonon frequencies and by applying theRTA, it was shown in Ref. [82] that FA mode conductivity contributions in the presence oftensile strain grow without bound in the limit of large crystalline domain sizes. Althoughthese findings with regard to thermal transport in graphene were later called into questionand attributed to an inadequacy of the RTA [120, 121], the results in Ref. [82], nevertheless,give an unequivocal demonstration of the importance of lattice strain as far as low-frequencythree-phonon scattering rates in two-dimensional materials are concerned.In the case of 1D CNTs, some analytical considerations of low-frequency three-phonon55scattering rates within the PBTT formalism were first presented by Mingo and Broido [53].By simplifying the exact scattering formulas inherent to lowest-order anharmonic pertur-bation theory, they derived approximate long-wavelength scaling relations capturing thethree-phonon scattering strength of acoustic phonons and thereby pointed out the issue ofheat carrying modes with vanishing scattering rates. Later, more refined calculations ofthree-phonon scattering rates and of heat transport in free-standing CNTs were performedby us, see Chapter 4, and others [47, 48, 56, 59, 60]. With the help of numerical rou-tines, these studies retained the full complexity associated with the many-body problem ofinteracting phonons but did not present comprehensive analytical considerations of exactlow-frequency scaling relations of three-phonon scattering rates, that can be inferred fromlowest-order anharmonic perturbation theory [82]. More so, in the case of CNTs, the effectsof strain on low-frequency three-phonon scattering rates involving out-of-line FA modeshave so far not been rigorously addressed in the PBTT formalism.To fill these gaps in our knowledge of CNT heat transport, hereinafter, we aim toelucidate generic trends of acoustic mode interactions in the light of one-dimensional phononconfinement and strain hardening of low-frequency FA modes. To this end, without loss ofgenerality, we perform PBTT calculations on (4,4) CNTs at room temperature and imposea strain-dependent lattice parameter a = (1 + )a0, where  ≥ 0 denotes a variable tensilestrain amplitude. Tensile lattice strain  then enters naturally into our PBTT calculations byinfluencing second and third order atomic force constants. Importantly, while our previousefforts in Chapter 4, as well as earlier related studies on stress-free CNTs [47, 48, 56, 59,60], were concerned with reporting numerical lattice-dynamical results, here we complementour PBTT calculations with analytical considerations of strain-dependent three-phononscattering rates in the low-frequency limit. In particular, our considerations hark backto earlier CNT continuum modeling results [122–125] which help us, in part, to put ournumerical lattice-dynamical PBTT results on firm theoretical footing. In this way, at thelevel of the RTA, we derive unambiguous qualitative trends of the thermal conductivitywith respect to the tube length, whereby we are able to highlight the critical role of tensilelattice strain in ultra-long heat conducting CNTs.The present chapter is organized as follows. In regard to the low-frequency phonondynamics in CNTs, we first recall some relevant continuum modeling results in Section 5.2.Concentrating on the dissipation of long-wavelength acoustic phonons in stress-free andstretched tube configurations, we present lattice-dynamical calculations of three-phononscattering rates in Section 5.3. Here, we show numerical data at the level of individual three-phonon scattering channels and specify exact long-wavelength scaling relations, wherebydetailed analytical derivations are deferred to Appendix C. In Section 5.4, we use thesefindings to revisit RTA-level predictions of CNT phonon transport. We finally discuss56further research directions in Section 5.5 and conclude in Section 5.6.5.2 Continuum models and anharmonic scattering channelsTreated as isotropic, free-standing, linear elastic continua, both long tubes and solid rods ofradius R are known to permit three distinctive types of low-frequency waves: twisting (TW),longitudinal (LA) and two degenerate flexural bending type (FA) modes, exhibiting linear,ωLA/TW ∼ |k|, and quadratic, ωFA ∼ Rk2, axial wave dispersion in the long-wavelength limit|k|R 1 [126]. For CNTs, these acoustic waves can be reproduced by a simple continuummodel of a cylindrical surface [122, 123]. Specifically, by assuming that CNTs inherit the in-plane elastic isotropy of graphene sheets [122, 123], the acoustic wave dispersion for |k| → 0becomesωLA(k) = vLA|k|, vLA =√Yρm,ωTW(k) = vTW|k|, vTW =√µρm,ωFA(k) =vLAR√2k2,(5.2.1)where the LA and TW wave velocities are determined by the 2D Young’s modulus Y , the 2Dshear modulus µ, and the area mass density of the cylinder surface ρm. More sophisticatedcontinuum models have to be employed to account for tube chirality effects and elasticanisotropy. Considering both chiral and achiral CNTs, recently, Gupta and Kumar [125]applied a Cosserat rod theory in order to derive analytical expression for the acoustic wavedispersion under external loads. Following Ref. [125], to the lowest order in the tensilestrain amplitude , frequency hardening of FA modes manifests itself in CNTs asω2FA(k, ) = ω2FA(k, 0) + T1 k2 − T2 k4, (5.2.2)where ωFA(k, 0) denotes the FA mode dispersion in the stress-free state and the Ti’s stand forsome positive constants which relate to the tube’s stretching and bending stiffness. For ouratomistic model of a (4,4) CNT, acoustic mode frequencies obtained from lattice-dynamicalcalculations are plotted in Figure 5.1. From LA and TW mode frequencies in unstrainedtubes, we infer vLA = 21.37 km/s and vTW = 13.57 km/s, respectively. At  = 2.0 %,LA mode velocities are found to be smaller by 2.9 %, whereas TW mode velocities areslightly increased by 0.2 %. Most importantly, in the presence of tensile lattice strain, long-wavelength FA modes transition from a quadratic to a linear wave dispersion as is faithfullydescribed by Equation (5.2.2).57Figure 5.1: Log-log plot of the acoustic phonon dispersion for a (4,4) CNT withvarying tensile strain amplitude . Empty plot markers correspond to lattice-dynamicalcalculations that take only strain-dependent second order force constants as input.Solid lines are obtained by fitting the continuum predictions as per Equations (5.2.1)and (5.2.2) to phonon data in the region k ≤ 0.1pi/a. Here, ωLA = 42.56 k, ωTW =27.04 k, and ω2FA = (10.61 − 73.70 ) × 103 k4 + 1.75 × 103k2 for LA, TW, and FAmodes, respectively, with ων in units of THz and k in units of pi/a, where a is thestrain-dependent lattice constant.Continuum models might be further invoked to provide insight into acoustic phonon-phonon interaction processes. Within lowest-order anharmonic perturbation theory, fora given acoustic phonon of low frequency and long wavelength, σ ≡ (ν, k), labeled bydispersion branch ν and wave number k, two types of interacting phonon triplets (σ, σ2, σ3)are permitted by the law of energy conservation:type (i): ωσ ≈ ωσ2 ≈ ωσ3 ωOPT, (5.2.3a)type (ii): ωσ ωσ2 ≈ ωσ3 . (5.2.3b)Type (i) triplets comprise three acoustic phonons whose energies fall below the lowest op-tical phonon energy ωOPT, whereas type (ii) triplets entail two high-energy phonons ofnear identical energy. A continuum approach lends itself to systematically analyze an-harmonic scattering channels of type (i). Treating stress-free CNTs as isotropic hollowcylinders and applying tools of nonlinear elasticity theory [127], it was shown by De Mar-tino et al. [124] that possible interacting acoustic mode triplets reduce to (LA,LA,LA),58(LA,TW,TW), (LA,FA,FA), and (TW,FA,FA) for which three-phonon coupling coefficientsVσ,σ2,σ3 , see Equation (2.2.16), can be derived in terms of second and third order elas-tic constants. In particular, a noteworthy result of Ref. [124] is the quartic wave numberdependence VTWFAFA ∼ kk2k3(k2−k3)R, which is at variance with the standard cubic long-wavelength approximation Vσ,σ2,σ3 ∼ kk2k3 for type (i) triplets commonly found in classictextbooks [50, 128]. Within a lattice-dynamical treatment of three-phonon scattering chan-nels, this behavior in CNTs was previously overlooked [56]. As we show below, the weakcoupling strength between TW and FA modes gives rise to an intricate scattering behaviorof low-frequency TW modes in stress-free and stretched tubes.5.3 Channel-by-channel decomposition for stress-free andstretched configurationsTaking into account the full spectrum of three-phonon interactions including type (i) and(ii) phonon triplets, in the following, we carry out lattice-dynamical calculations of low-frequency three-phonon scattering rates for our model of a (4,4) CNT. Specifically, forvarying strain amplitudes , we analyze phonon-phonon interactions of acoustic modes at thelevel of individual three-phonon decay and absorption channels, considering a decompositionof scattering rates (see Equation (2.2.24)) of the form1τ3phν (k)= Γν(k) =∑ν2,ν312Γ−ν,ν2,ν3(k) + Γ+ν,ν2,ν3(k) (5.3.1)for small wave numbers k and ν ∈ {LA, TW, FA}. This enables us to identify type (i)and (ii) three-phonon channels and to single out the dominant contribution of anharmonicphonon-phonon dissipation in the long-wavelength limit. In turn, based on analytical con-siderations of these channels, the exact low-frequency scaling relations of three-phononscattering rates Γν ∼ rks as prescribed by lowest-order anharmonic perturbation theorycan be established.For brevity in notation, we will signify mode polarizations appearing in triple subscriptsin Equation (5.3.1) by single letters. For example, a scattering channel of type (ii) involvingone low-frequency LA mode and two high-frequency optical modes will be denoted by Γ+LOO.5.3.1 Longitudinal modesGiven the large group velocity of LA modes, selection rules imposed by energy and quasi-momentum conservation imply a severely reduced three-phonon scattering phase space avail-able to long-wavelength LA modes. In Figure 5.1, we show total LA scattering rates undervarying tensile strain amplitudes at different wavelength scales. Decomposing our data for59Figure 5.1: Anharmonic scattering rates (left) of LA modes at T = 300 K anddecomposition (right) into individual three-phonon scattering channels in the regionk ≤ 0.04pi/a. Top and bottom panels show stress-free,  = 0, and stretched,  > 0,configurations of a (4,4) CNT, respectively. Lattice-dynamical predictions as per Equa-tion (5.3.1) are represented by empty plot markers. Solid lines are obtained by sub-stituting long-wavelength approximations into Equation (5.3.1). The strain-inducedfrequency hardening of FA modes brings with it a long-wavelength scaling transitionfrom ∼ k−1/2 in stress-free to ∼ −5/2k2 in strained tubes.small k, we detect three different types of scattering channels that contribute to the total LArate in the long-wavelength limit: LA→ FA+FA, LA→ TW+TW, and LA+OPT→ OPTamong which the first largely dominates in both unstrained and stretched CNTs.In the absence of lattice strain, one finds Γ=0LA ≈ Γ−LFF ∼ k−1/2 which can be traced backto the long-wavelength scaling relations ωFA ∼ k2 and VLFF ∼ kk2k3 (see Appendix C.1).As can be noted from the top left-hand panel of Figure 5.1, the square root scaling breaksdown at smaller wavelength scales k > 0.1pi/a. This threshold conforms with the phonondispersion data in Figure 5.1, where it can be seen that FA mode frequencies deviate fromthe quadratic trend once k > 0.1pi/a, indicating a regime where atomistic details enteringinto the exact phonon dispersion become important.In the presence of tensile lattice strain,  > 0, lattice-dynamical predictions in Figure 5.1indicate a reduction of long-wavelength rates with increasing strain amplitude and a transi-tion away from the square root scaling in unstrained tubes. To explain the strain dependenceof ΓLA, we resolve Γ−LFF as per Equation (5.3.1) in the continuum limit,∑k →∫dk, whichrequires knowing the functional wave number dependencies of ωLA, ωFA and VLFF. Here, weadopt the dispersion fit lines of Figure 5.1 and, under the assumption that strain has no effecton the cubic wave number dependence of three-phonon couplings, we set VLFF = C3 kk2k3.60Representing the three-phonon coupling strength, the only unknown parameter C3 is thenfitted to lattice-dynamical data of Γ−LFF. This procedure yields the solid lines in Figure 5.1for tubes in unstrained and strained configurations. In fitting our lattice-dynamical pre-dictions, we do not observe a systematic variation of C3 under varying strain amplitudes,C3 , which leaves ωFA as the only strain dependent quantity entering into Γ−LFF. Thus,the strain induced downward shift of LA rates can be solely attributed to the frequencyhardening of FA modes.The implications of a linearized FA mode dispersion in the limit k → 0+ as per Equa-tion (5.2.2) are best understood in the phonon frequency domain, where the density of FAstates DFA within the low-frequency interval ω + dω goes from DFA ∼ ω−1/2 in unstrainedto DFA ∼ ω0 in strained systems. Loosely speaking, under increasing tensile lattice strain,there are less and less preferential FA states that a given long-wavelength LA mode candecay into. Assuming a strictly linear FA mode dispersion with strain dependent prefactor,ωFA ∼ 1/2|k|, one finds Γ>0LA ≈ Γ−LFF ∼ −5/2k2 (see Appendix C.1).The scattering amplitudes Γ−LTT and Γ+LOO, both of which do not involve long-wavelengthFA modes, remain relatively insensitive to strain. A scaling behavior Γ−LTT ∼ k2 is de-ducible from a linear dispersion ωTW ∼ |k| and a coupling coefficient VLTT ∼ kk2k3 (seeAppendix C.3). To produce the solid lines in Figure 5.1 for Γ−LTT, as before, the unknownprefactor of VLTT is taken as a fit parameter to our lattice-dynamical predictions. In thecase of optical absorption processes, numerical noise inherent to our evaluation of Equa-tion (5.3.1) becomes apparent to the extent that no continuous trend lines for Γ+LOO arediscernable. Notwithstanding the uncertainty in our prediction of Γ+LOO, data points in Fig-ure 5.1 suggest that optical absorption processes as well as TW mode interactions contributenegligibly to the dissipation of long-wavelength LA modes.5.3.2 Flexural modesGenerally speaking, the three-phonon scattering phase space for small-wave number FAmodes is reduced in that no lower-lying phonon branch providing pathways for kinemati-cally allowed decay processes exists. Our lattice-dynamical predictions of strain dependentFA scattering rates are summarized in Figure 5.2. In the long-wavelength limit k → 0+,traceable scattering contributions arise from two types of acoustic absorptions that dependcritically on tensile strain, namely FA + FA → LA as well as FA + FA → TW, and fromoptical interbranch absorptions FA + OPT→ OPT which appear unaffected by strain.As is evident in Figure 5.2, a reduction in the total scattering rate of long-wavelength FAmodes under increasing tensile lattice strain derives from the dominant process FA + FA→LA. The same three-phonon coupling coefficient VLFF enters into decay Γ−LFF and absorptionamplitudes Γ+FFL. Hence, according to our earlier observation VLFF  in the case of61Figure 5.2: Same as Figure 5.1 but for FA modes. For stress-free CNTs, total FAscattering rates converge towards a constant, whereas tensile lattice strain gives rise toan asymptotic scaling ∼ −1k2.LA modes, a reduction in the total rate should be solely ascribable to the strain-inducedhardening of FA modes. Indeed, by setting VLFF = C3kk2k3 and by fitting our lattice-dynamical predictions of Γ+FFL for different tensile strain amplitudes , we find no systematicvariation of C3 with increasing tensile strain, leaving only ωFA as a strain-dependent factor.Based on analytical considerations of Γ+FFL, adopting either a quadratic ωFA ∼ k2, or astrictly linear FA mode dispersion ωFA ∼ 1/2|k|, stress-free and stretched tubes give rise toΓ=0FA ≈ Γ+FFL ∼ k0 and Γ>0FA ≈ Γ+FFL ∼ −1k2, respectively (see Appendix C.2).In the case of the absorption channel FA + FA→ TW, the scaling relations Γ+FFT ∼ k2in unstrained and Γ+FFT ∼ −1k4 in strained tubes stem from a quartic dependence VTFF ∼kk2k3(k2 − k3) (see Appendix C.2). In our lattice-dynamical calculations of VTFF, contri-butions that are cubic in the wave number k cancel out, in agreement with elasticity theorypredictions [124]. Similar cancellation effects manifest in the case of optical absorption pro-cesses FA+OPT→ OPT. Following a standard argument [128], for a given long-wavelengthacoustic phonon (ν, k), one expects a linear scaling VνOO ∼ k to lowest order in k. In orderto explain the scaling Γ+FOO ∼ k4 of optical interbranch transitions in Figure 5.2, however,we have to assume VFOO ∼ k2 (see Appendix C.4).5.3.3 Twisting modesFor long-wavelength TW modes, we have ωFA < ωTW < ωLA and, as compared to LA and FAmodes, a kinematically less restricted three-phonon scattering phase space is anticipated. Asshown in Figure 5.3, phonon-phonon transitions that contribute to the dissipation of long-62Figure 5.3: Same as Figure 5.1 but for TW modes. A third row of panels is added todistinguish the strain regimes 0 <  < 4 % (middle) and  ≥ 4 % (bottom). Processesinvolving high-frequency optical modes are further subdivided into contributions result-ing from intrabranch, TW + OPTa → OPTa and interbranch, TW + OPTa → OPTb,transitions whose long-wavelength scaling is ∼ k and ∼ k2, respectively. For stress-freetubes, the total scattering rate scales as ∼ k1/2 in the limit k → 0+. Under finitepositive strain,  > 0, the total rate develops a dependence ∼ k. For  ≥ 4 %, opticalintrabranch transitions are kinematically forbidden and the total TW scattering rateinherits a dependence ∼ k2.wavelength TW modes can be decomposed into two types of low-frequency acoustic and twotypes of high-frequency optical scattering channels. The former channels are TW→ FA+FAand TW + TW → LA, while the latter are given by intra- and interbranch absorptions,TW + OPTa → OPTa and TW + OPTa → OPTb, respectively. The relative importance ofthese channels varies with tensile strain, causing a nontrivial strain dependence in the totalrate ΓTW.In the case of unstrained tubes,  = 0, the dominant scattering contributions stemfrom the decay into two flexural modes, Γ−TFF, and from optical intrabranch absorptions,Γ+TOaOa , both of which are roughly of the same order of magnitude within the resolved rangeof wave numbers. For the decay TW → FA + FA, we have Γ−TFF ∼ k1/2 which is derivablegiven a quartic scaling of three-phonon couplings VTFF ∼ kk2k3(k2− k3) and provided thatωFA ∼ k2 under stress-free conditions (see Appendix C.1). Considering optical absorption63processes of acoustic modes, the scaling behavior of related scattering amplitudes dependson whether momentum is transferred between two crossing optical branches or within thesame branch. If we adopt the long-wavelength approximation VTOO ∼ k [128], we find thatvariations in the scattering phase space lead to Γ+TOaOa ∼ k and Γ+TOaOb ∼ k2, explainingthe differing scalings of intra- and interbranch transitions in Figure 5.3, respectively (seeAppendix C.4).In the presence of tensile strain,  > 0, as before in the case of Γ−LFF, Γ+FFL, and Γ+FFT,frequency hardening of long-wavelength FA modes strongly suppresses the scattering am-plitude of the decay TW → FA + FA. Applying a fit to our lattice-dynamical predictionsof Γ−TFF, we again notice that strain has no influence on the coupling coefficient VTFF,whose quartic wave number dependence together with ωFA ∼ 1/2|k| yields Γ−TFF ∼ −7/2k4(see Appendix C.1). For moderately strained tubes, total TW scattering rates in the long-wavelength limit remain lower bounded by optical intrabranch transitions and the limitingbehavior transitions from ∼ k1/2 to ∼ k.Interestingly, at larger strain amplitudes, long-wavelength scattering contributions aris-ing from optical intrabranch absorption processes may become susceptible to lattice strainas well. For our model of a (4,4) CNT, we find that optical intrabranch transitions canoccur only in the regime  < 4 % beyond which the total TW rate abruptly gains a k2dependence resulting from optical interbranch transitions. While no structural changes canbe observed in the carbon lattice at  = 4 %, an explanation for the discontinuity in thetotal TW rate is found by considering the strain dependence of optical phonon frequencies,see Figure 5.4. A necessary condition for the realization of intrabranch transitions is theexistence of a supersonic optical branch ωOPT∗ . As illustrated in Figure 5.4(a), one suchbranch whose slope ∂ωOPT∗/∂k surpasses the slope of the TW branch can be singled outin the spectrum of high-frequency optical modes. Under tensile loading, the branch ωOPT∗softens, see Figure 5.4(b), and its slope decreases. Figure 5.4(c) demonstrates that for strainamplitudes  ≥ 4 % the necessary condition ∂ωOPT∗/∂k > vTW for optical intrabranch tran-sitions is no longer fulfilled. A somewhat pathological prediction reveals itself as the criticalstrain amplitude  = 4 % is approached from below. Here, scattering rates are predicted toincrease, which can be traced to the vanishing curvature, ∂2ωOPT∗/∂k2 → 0, in the wavenumber regions depicted by the insets of Figure 5.4(a), cf. Equations (C.4) and (C.4.2).In passing, we can compare TW and LA mode velocities in Figure 5.4(c) to conclude thatlong-wavelength LA modes remain protected from optical intrabranch transitions.We point out that a recent numerical study on unstrained nanowires by Rashid etal. [129] found low-frequency power-law dependencies of three-phonon scattering rates forLA, FA, and TW modes that conform to our results for unstrained tubes. From the view-point of acoustic phonons, nanowires and nanotubes are very similar. It seems therefore64Figure 5.4: Discontinuity in the strain dependence of scattering rates of long-wavelength TW modes due to optical intrabranch transitions. (a) The phonon disper-sion of a (4,4) CNT at  = 0 comprises one optical branch (OPT∗) whose slope matchesthe slope of the TW branch, enabling intrabranch transitions TW + OPTa → OPTa.(b) Frequency softening of the OPT∗ branch under tensile lattice strain. (c) The slopeof the OPT∗ branch decreases with increasing strain amplitude. At  = 4 %, the veloc-ity vTW exceeds ∂ωOPT∗/∂k and the scattering of long-wavelength TW modes due tointrabranch transitions is kinematically prohibited.natural to conjecture that considerations in the present study apply to a broader class ofone-dimensional materials, implying, for example, that a strain-induced reduction in thescattering strength of long-wavelength acoustic phonons can be achieved in nanowires aswell.5.4 Thermal conductivity under the relaxation timeapproximationWithin the PBTT formalism, the power-law dependencies of three-phonon scattering rates1/τ3phν ∼ rks, as found in Section 5.3, shed new light on size effects in CNT heat trans-port. At the RTA level, the numerical challenge of converging PBTT-derived conductivityresults with respect to the employed wave number grid resolution, encountered previouslyin Chapter 4, can be fully clarified. Even more interesting, with the long-wavelength scaling65relations 1/τ3phν ∼ rks at hand, unambiguous statements concerning the long-tube limitof phonon transport, κRTA(L→∞), can be made. In particular, for stress-free tubes, wehave the asymptotic limits 1/τ3phLA ∼ |k|−1/2, 1/τ3phTW ∼ |k|1/2, and 1/τ3phFA ∼ k0, which im-plies an intrinsically upper bounded thermal conductivity value for macroscopically longCNTs, κRTA ∼ ∫0 dk κν(k) < ∞ for L → ∞, see Equations (4.1.3) and (4.1.4). In thepresence of tensile lattice strain, on the other hand, scattering rates vanish more rapidly,1/τ3phν ∼ r|k|s with r ≤ 0 and s ≥ 1. That is to say, at the RTA level, the integration ofper-branch conductivity contributions becomes divergent and stretched CNTs are predictedto exhibit an ever-increasing thermal conductivity value in the long-tube limit.Usually, in lattice-dynamical PBTT calculations, per-phonon thermal conductivities aresampled and summed over a finite 1D wave number grid, κ = N−1∑ν∑k κν(k), with gridspacing ∆k = 2pi/Na, which requires convergence tests with respect to N . This approach,however, becomes problematic when one has to account for a divergence of anharmoniclifetimes, τ3phν (k → 0+) → ∞, as is the case for LA and FA modes in strained and forTW modes in both unstrained and strained tubes. The imposed tube length L, whichenters into calculations via 1/τ0σ = 1/τ3phσ + 2|vσ|/L, effectively filters out weakly dampedacoustic modes, but, with L increasing and the term 2|vσ|/L becoming insignificant, suf-ficiently dense wave number integration grids to converge the thermal conductivity mightbe out of reach by means of lattice-dynamical calculations. Here, to tackle these finite-grid-resolution issues for our model of a (4,4) CNT, which we assume to be of large butfinite length L, we use the asymptotic fits to dispersion laws and three-phonon scatteringrates as depicted in Figures 5.1 and 5.1–5.3 and extrapolate the wave number dependenciesof acoustic per-branch conductivities, κν(k → 0+). In this way, by numerically integratingthese extrapolated trends below a certain cut off wave number kc, κLW =∫ kc0 dk κν(k → 0+)for ν ∈ {LA,FA,TW}, we are able to fully account for long-wavelength conductivity con-tributions to κRTA, that would otherwise be truncated in standard PBTT calculations.In Figure 5.1, we combine per-phonon conductivities sampled and summed over a gridincluding N = 1000 k-points, κΣ, with integrated acoustic mode contributions κLW stem-ming from k ∈ [−2pi/1000a; 2pi/1000a], a small region around the Gamma point k = 0 notcaptured by the imposed grid, to show RTA-level conductivity predictions in the limit ofsmall strain amplitudes  and large tube length scales L. The relative importance of summedκΣ vs integrated long-wavelength contributions κLW is displayed in Figure 5.2. As can beseen in Figure 5.1, relatively long tubes L > 100µm are required to observe an enhancementof thermal transport due to a reduction in the scattering rates of long-wavelength acousticmodes in the presence of strain. For a tube with fixed but very large L, low-frequencyconductivity contributions are predicted to be significantly enhanced with increasing strainamplitude. Shorter tubes, however, at the micrometer scale and below, do not support66Figure 5.1: Lattice thermal conductivity under the relaxation time approximationfor a (4,4) CNT at T = 300 K as obtained from strain-dependent three-phonon scat-tering rates. In stress-free tubes,  = 0, the conductivity converges with tube lengthL. Under finite tensile strain,  > 0, and in the limit of long tubes, the acousticphonon contributions to the lattice thermal conductivity scale as κRTALA/FA ∼ 5/4L1/2and as κRTATW ∼ lnL (∼ L1/2) for  < 4 % (≥ 4 %).67such low phonon frequencies, and therefore a strain-induced enhancement of phonon trans-port cannot manifest. In fact, given that the majority of vibrational frequencies in solidstypically softens under tensile strain, as is captured by a positive mode-specific Gru¨neisenparameter [118, 130], a strain-induced reduction of heat transfer rates in sub-micrometerlong CNTs may be anticipated instead.Focusing on the length-divergent trends of the thermal conductivity in Figure 5.1, exactpower-law dependencies κ ∼ rLα for small nonzero strain amplitudes  > 0 can be derived.In the case of LA modes, it follows from 1/τ3phLA ∼ 5/2k2 that conductivity contributions ofLA modes scale as κRTALA ∼ 5/4L1/2. The same scaling relation is obtained in the case of FAbranch contributions, κRTAFA ∼ 5/4L1/2, which, in this case, is attributable to 1/τ3phFA ∼ −1k2as well as ωFA ∼ 1/2|k|. For TW modes, the long-wavelength scaling relation 1/τ3phTW ∼ |k|leads to a logarithmically divergent conductivity contribution in tubes subject to strainamplitudes 0 <  < 4 %. The discontinuity in the strain dependence of long-wavelengthTW scattering rates at  = 4 %, as detailed in Figure 5.4, is filtered out in short tubes butbecomes apparent if L > 10µm. For strain amplitudes  ≥ 4 %, a divergence κRTATW ∼ L1/2 isimplied by 1/τ3phTW ∼ k2. In the bottom panel of Figure 5.1, we show total thermal conduc-tivity predictions including optical branch contributions which, however, do not increasesignificantly with tube length and are therefore of minor importance in the large-L regime.As can be inferred from Figure 5.2, per-branch acoustic conductivity data sampledover a grid with N = 1000 discrete wave numbers gives a poor representation of the totalthermal conductivity as the tube length L becomes large. In unstrained tubes,  = 0,it is the square root divergence of intrinsic TW lifetimes τ3phTW (k → 0+) ∼ k−1/2 whichrequires grid calculations with N 1000 to converge conductivity predictions in the long-tube limit. Convergence issues with respect to the wave number grid resolution N aregenerally exacerbated in the presence of tensile strain,  > 0, where all three acoustic modepolarizations exhibit a singularity of the form τν(k → 0+) ∼ k−r with r ≥ 1.5.5 DiscussionTwo approximations have been employed in the present study that deserve further scrutinyin the future. For one, phonon frequency shifts at finite temperature induced by latticeanharmonicity have been assumed to be negligible. In the case of two-dimensional crystals,by taking into account the real part of the phonon self-energy Re[Σσ(ω)], it was noted [131]that out-of-plane vibrations are renormalized at finite temperature to the extent that thequadratic dispersion in stress-free samples under the harmonic lattice approximation be-comes linear in the long-wavelength limit. For unstrained tubes, the obtained scaling lawsof three-phonon scattering presented in Section 5.3 critically depend on the condition thatFA mode frequencies vanish quadratically in the long-wavelength limit, ωFA ∼ k2. Taking68Figure 5.2: Relative importance of summed κΣ vs integrated κLW thermal conductiv-ity contributions of acoustic modes as a function of tube length L and strain amplitude. κΣ corresponds to a summation over a wave number grid with N = 1000 k-points.Long-wavelength conductivity contributions are captured by κLW but not included inκΣ. A value of −1 on the vertical axis implies a convergence of thermal conductivitywith respect to the wave number grid resolution N , whereas larger values justify anextrapolation approach.a temperature-dependent, possibly linearized, FA mode dispersion as input to a perturba-tive calculation of three-phonon scattering rates, low-frequency dissipation rates of acousticmodes could be drastically decreased, as we demonstrate here by examining tubes undertensile load. In a similar context, while studying phonon-phonon interactions in a semiflex-ible monoatomic chain, Santhosh and Kumar [132] pointed out the necessity of anharmonicfrequency renormalization, since the standard perturbative treatment of harmonic phononmodes violates the notion of well-defined quasi-particles, requiring ων(k)τ3phν (k) > 1 ask → 0+. In particular, according to Ref. [132], longitudinal, ωx ∼ k, and transverse,ωy ∼ k2, modes of a semiflexible chain under the harmonic lattice approximation give riseto 1/τ3phx ∼ |k|−1/2 and 1/τ3phy ∼ k0, in agreement with our findings for LA and FA modesin unstrained tubes, cf. Sections 5.3.1 and 5.3.2. Effective long-wavelength dissipation ratesunder phonon renormalization, as predicted in Ref. [132], are instead significantly reduced.Another point worthy of future research relates to the approximation of the thermal con-ductivity. Previous studies emphasized the inadequacy of the RTA in the case of one- [48,6956, 59, 133, 134] and two-dimensional [120, 121, 135, 136] materials, noting the impor-tance of obtaining the thermal conductivity from beyond-RTA methods in the frameworkof the PBTT. It is often assumed that the RTA generally underestimates the true thermalconductivity [50, 137]. Although such a statement cannot be proven rigorously [138], anunderestimation of the thermal conductivity by the RTA when compared to beyond-RTAresults appears to be confirmed by our iterative calculations in Chapter 4, see Figure 4.2.Therefore, given the RTA-derived thermal conductivity results presented in Section 5.4, oneis tempted to conclude that CNTs under tensile strain act as superdiffusive heat conductorswhich give rise to κ ∼ Lα with α ≥ 1/2 in the long-tube limit. In the case of graphene,a similar reasoning based on the RTA was applied by Bonini et al. [82] to conclude thatthe thermal conductivity of isotropically strained sheets should be logarithmically divergentwith domain size. In pursuing numerical calculations beyond the RTA, however, some ofthe same authors later arrived at the conclusion that the conductivity of pristine graphenesheets should remain upper bounded, even in the presence of strain [120]. By taking into ac-count the power-law dependencies of three-phonon scattering rates obtained in Section 5.3,it would be desirable to go beyond the RTA, aiming for more rigorous analytical results inregard to the length-scaling of the CNT thermal conductivity.5.6 ConclusionsTaking as example a (4,4) CNT, we have performed lattice-dynamical calculations in theframework of the PBTT to address the role of acoustic phonons in low-dimensional nan-otube heat transport in light of tensile lattice strain. By supporting our calculations withanalytical considerations and by making recourse to continuum theories, we derived thegeneral long-wavelength scaling relations of acoustic phonon-phonon scattering rates thatfollow from lowest-order anharmonic perturbation theory. Based on our model, the onsetof long-wavelength behavior in CNTs of radius R is predicted at around |k|R < 0.3.As was shown earlier in the case of graphene [82], the three-phonon scattering strength oflong-wavelength acoustic modes can be reduced by tensile lattice strain whenever acoustic-flexural mode type processes are dominant. The scattering phase space for these processesis greatly reduced by a strain-induced hardening of flexural mode frequencies. In goingfrom a two-dimensional phonon scattering phase space in graphene to a one-dimensionalone, this effect becomes even more pronounced in CNTs. As compared to longitudinalmodes, long-wavelength twisting modes in CNTs couple only weakly to flexural modes.Twisting modes have divergent anharmonic lifetimes in the long-wavelength limit in bothunstrained and strained tube configurations. In standard calculation approaches, ever-increasing acoustic mode lifetimes in CNTs pose a severe numerical challenge when accuratethermal conductivity estimates for long CNTs are desired.70For CNTs under tensile lattice strain, RTA-level predictions provide strong evidencefor the nonexistence of an upper bound to thermal conductivity, suggesting the possibilityof superdiffusive heat transport in macroscopically long nanotubes. Whether beyond-RTAapproaches, or the inclusion of additional physics, such as finite-temperature effects or thesystematic inclusion of higher-order phonon-phonon scattering processes, help to renormal-ize a formally divergent thermal conductivity awaits further investigation.71Chapter 6Revisiting phonon transportbeyond the relaxation timeapproximation6.1 IntroductionIn parallel to our research efforts presented in Chapters 4 and 5, a breakthrough in thecontext of atomistic modeling of CNT heat transport was recently reported by Barbali-nardo et al. [59]. These authors presented PBTT and MD-derived data for the room tem-perature thermal conductivity of a free-standing (10,0) CNT, arriving at two remarkableconclusions: (i) going from nanometer-sized to macroscopically long CNTs, tube-length-dependent predictions of heat transfer rates can be brought into quantitative agreementacross the two rather different computational approaches, and (ii), by applying the PBTTformalism and by going beyond the RTA, a length-converged value of about 9960 Wm−1K−1can be deduced in the long-tube limit, yielding a ballpark estimate for macroscopically longCNTs in experiments. In the PBTT calculations of Ref. [59], a new beyond-RTA solutionapproach to the linearized phonon transport equation was pursued, previously not yet ex-plored in the case of CNTs. Specifically, a matrix inversion approach was deemed necessaryto settle the long-standing debate about converging vs diverging trends of a tube-length-dependent CNT thermal conductivity.In the present chapter, we comment on the PBTT calculations of Barbalinardo et al. [59].As a first step, in Sections 6.2 and 6.3, we expand on the PBTT formalism introduced inSection 2.2, to give a background for the matrix inversion approach. From the linearizedphonon transport equation, we derive the so-called collision matrix. The so-introducednotation in turn allows us to clarify the role of umklapp scattering and, furthermore, to72re-express the beyond-RTA phonon conductivity of CNTs in terms of a matrix inverse. InSection 6.4, we put this new approach of calculating heat transfer rates via matrix inversionto use, focusing explicitly on possible upper bounds of conductivity values in the long-tube limit. Here, strong numerical evidence refuting (i) and (ii) in Ref. [59] is presented.In particular, with the help of analytical insights from Chapter 5, we show that spuriouslength-convergent PBTT results in Ref. [59] emanate from low-precision force constants.6.2 The collision matrix and the necessity of umklappscatteringIn the PBTT formalism, the inseparable connection between umklapp scattering processesand finite thermal resistance due to lattice anharmonicity can be demonstrated by consid-ering the linearized phonon transport equation of Section 2.2 for a CNT in a scenario wherethere is no temperature gradient across the system, ∂T/∂z = 0. In this case, the drift termin Equation (2.2.1) vanishes and, accounting for three-phonon interactions, the linearizedtransport equation, see Equations (2.2.6) and (2.2.18), can be written as0 = − 1n0σ(n0σ + 1)n˙σ∣∣∣scatt=∑σ2,σ312[Ψσ −Ψσ2 −Ψσ3 ] Γ−σ,σ2,σ3+ [Ψσ + Ψσ2 −Ψσ3 ] Γ+σ,σ2,σ3 .(6.2.1)The functions Ψσ measure mode-specific deviations from the equilibrium phonon occupationn0σ, and, as per Equation (2.2.3), any nonvanishing heat flux is linked to a nontrivial solutionΨσ 6= 0 of Equation (6.2.1). Therefore, to model the correct physics in the absence of atemperature gradient, we might expect that Equation (6.2.1) admits only the trivial solutionΨσ = 0. The corresponding situation can be conveniently re-expressed in matrix-vectornotation as0 = C∞Ψ with C∞σσ′ =( ∑σ2,σ312Γ−σ,σ2,σ3 + Γ+σ,σ2,σ3)δσ,σ′+∑σ3Γ+σ,σ′,σ3 − Γ+σ,σ3,σ′ − Γ−σ,σ′,σ3 ,(6.2.2)where the matrix C∞ ∈ R3NN×3NN will be referred to as the collision matrix in the infinitesize limit (see further explanation below), and the vector Ψ ∈ R3NN subsumes all phonondeviation functions Ψσ for a given wave number grid with N k-points. For ∂T/∂z = 0,intuitively, we anticipate no other solution than the trivial one Ψ = 0, which is tantamountto saying that the collision matrix C∞ is nonsingular. Owing to the properties of three-phonon transition amplitudes Γ±, however, such reasoning is inaccurate and two solutionsexist that deserve special attention.73By virtue of the energy selection rule, three-phonon transition amplitudes Γ−σ,σ2,σ3 andΓ+σ,σ2,σ3 are nonzero if and only if ωσ − ωσ2 − ωσ3 = 0 and ωσ + ωσ2 − ωσ3 = 0, hold true,respectively. Hence, it is clear from Equation (6.2.1) that a nontrivial solution to 0 = C∞Ψis given by Ψ˜σ = ωσ. Fortunately, consistent with expectations, the so-obtained heat fluxvanishes, since Ψ˜σ = ων(k) is an even function in k, see Equation (2.2.3) and also Ref. [83].Underscoring the importance of umklapp scattering events, a physically more relevantsituation is described by Ψσ = k. Namely, if three-phonon scattering amplitudes in Equa-tion (6.2.1) were to be found without the inclusion of momentum-nonconserving umklappprocesses, we would find another nontrivial solution to 0 = C∞Ψ to be given by Ψσ = kwhich, in this case, is an odd function in k, leading to a nonzero heat flux in the absence ofa temperature gradient. Put differently, without umklapp and with normal scattering pro-cesses alone, there would be no mechanism by which a nonzero heat flux could be degradeddue to anharmonic phonon-phonon interactions, which would imply an infinite thermal con-ductivity. Obviously, if we consider the RTA, where the infinite-size conductivity is givenby κRTA ∝∑σ cσv2στ3phσ with lifetimes 1/τ3phσ = ∑σ2,σ3 Γ−σ,σ2,σ3/2+Γ+σ,σ2,σ3 , a finite conduc-tivity value is implied even for a situation where umklapp scattering processes are frozenout. For this reason, it is commonly assumed that the RTA generally underestimates theconductivity when compared to beyond-RTA calculation results of phonon transport [50,137]. For general three-phonon scattering spectra, comprising both normal and umklappevents, however, such an assumption may break down [138].6.3 Link to the iterative solutionExpanding on the matrix-vector formulation of phonon transport in Equation (6.2.2), PBTTconductivity calculations beyond the RTA essentially amount to solving a large system oflinear equations. To see this, we now turn on a temperature gradient and consider heattransport across CNTs of tube length L as per Equation (2.2.20), which is equivalent towritingξ = C(L)F with Cσσ′(L) = C∞σ,σ′ +2|vσ|Lδσ,σ′and ξσ = vσ~ωσkBT 2.(6.3.1)Here, an extra scattering term is added to the diagonal part of the infinite-size collisionmatrix of Equation (6.2.2) which accounts for ballistic upper limits of phonon transport.Further, as before in Section 2.2, we work with Ψ = (∆T/L)F . Now, if Equation (6.3.1)admits a unique solution for F , meaning that the matrix inverse of C(L) exists, then the74beyond-RTA phonon conductivity of CNTs follows asκinv =~2NaAckBT 2∑σ,σ′ωσωσ′vσvσ′n0σ(n0σ + 1) C−1σ,σ′(L). (6.3.2)Compared to iterative calculations introduced in Section 2.2.4, direct matrix inversion asper Equation (6.3.2) constitutes an alternative way of determining the phonon conductivitybeyond the RTA. In fact, to show the equivalence of κinv and κiter, we can use the matrix-vector notation at hand and reformulate the iterative procedure of Section 2.2.4 in terms ofdiagonal and off-diagonal terms of the collision matrix, C(L) = Cdiag(L)+Coff . In particular,it follows from Equation (2.2.26) thatF (i+1) =1Cdiag(L)ξ −1Cdiag(L)CoffF (i). (6.3.3)The so-obtained recursive expression can now be rewritten explicitly asF =∞∑p=0[− 1Cdiag(L)Coff]p 1Cdiag(L)ξ, (6.3.4)where it becomes clear that the iteration is nothing else but an approximation of the matrixinverse of C(L). Importantly, as can also be noted from Equation (6.3.4), the iteration maydiverge for L → ∞ in scenarios where the collision matrix C(L → ∞) = C∞ is no longerdiagonal dominant [139]. In this case, computationally more expensive numerical routinesfor solving ξ = C(L)F may be in demand.As was already shown in Chapter 4, iterative PBTT calculations of the CNT thermalconductivity in the macroscopic limit, κiter(L → ∞), diverge, which indicates the possi-bility of superdiffusive phonon-mediated heat transport in 1D materials. In light of Equa-tion (6.3.4), it remains to be checked whether such a prediction persists if an alternative,numerically more robust, beyond-RTA solution approach is pursued.6.4 Calculations of κ(L→∞) via direct matrix inversionFocusing on free-standing pristine CNTs at room temperature, it was claimed by Barbali-nardo et al. [59] that, through direct numerical inversion of the infinite-size collision matrixC∞, an intrinsically upper bounded phonon conductivity can be restored, κinv(L → ∞) <∞. Starting from the optimized Tersoff atomic interaction potential [72], these authorsindicated κinv(L→∞) = 9960 Wm−1K−1 for a (10,0) CNT. In the following, we pursuedirect numerical inversions of C∞ along the lines of Ref. [59] to solve ξ = C∞F as perEquation (6.3.1). Before we present our numerical results contradicting Ref. [59], we com-75ment, from a theoretical point of view, on the existence of (C∞)−1 and on the importanceof accurate force constant calculations in the framework of the PBTT. In particular, by ad-dressing the latter point, we clarify the cause for an observed spurious convergent behaviorof κinv(L→∞) in Ref. [59].6.4.1 The null space due to energy conservationAs noted before, a null space of the collision matrix C∞ is spanned by the vector Ψ˜ withcomponents Ψ˜σ = ωσ, which, in principle, implies the nonexistence of the matrix inverse(C∞)−1. Conceptually, despite noninvertibility, the system of linear equations ξ = C∞Fmay still be solved for F , given the orthogonality of the vectors ξ and Ψ˜ whose componentsare odd and even functions in k, respectively [83]. Therefore, strictly speaking, by pursuinga matrix inversion approach within the PBTT formalism, one should resort to the Moore-Penrose inverse [140].In actual numerical PBTT calculations, however, phonon triplets are sampled from thesame regular k-grid which inevitably introduces small violations of the condition of energyconservation. Hence, for a numerically constructed collision matrix C∞, the null spacespanned by Ψ˜σ = ωσ will have a tendency to be masked. Here below, following the approachof Ref. [59], we apply standard linear algebra routines to compute the actual inverse of C∞,which we then use for calculations of κinv(L→∞) according to Equation (6.3.2).6.4.2 Accuracy of third order force constantsIn Chapter 5, we analytically derived the asymptotic power-law dependencies of the kine-matically allowed three-phonon scattering channels in the low-frequency limit. Specifically,in the case of free-standing CNTs, we found the long-wavelength asymptotics of anharmonicacoustic phonon dissipation to be governed by the scattering amplitudes Γ−LFF ∼ |k|−1/2,Γ+FFL ∼ k0, and Γ−TFF ∼ |k|1/2. In numerical PBTT calculations, to correctly capture theseasymptotic wave number dependencies, we have to evaluate an expression of the form (seeAppendix C)Γ±ν,ν2,ν3(k) = ζν,ν2,ν3(k) |Vν,ν2,ν3(k)|2. (6.4.1)Taking as input only phonon frequencies, the function ζν,ν2,ν3(k) can be considered as achannel-specific “weighted phase space” [141, 142]. For the dominant scattering amplitudeslisted above, its asymptotic scaling behavior is strongly singular, ζν,ν2,ν3(k → 0+) ∼ k−awith a ≥ 9/2, see Equations (C.1.4) and (C.2.3). As such, the three-phonon couplingcoefficients Vν,ν2,ν3 are responsible for balancing the strong singularity of ζν,ν2,ν3 . Low-frequency, long-wavelengths scaling behaviors of Vν,ν2,ν3(k → 0+) are discussed in detail in76Figure 6.1: Lattice-dynamical predictions of anharmonic acoustic phonon lifetimesat T = 300 K in a (10,0) CNT, obtained using as input third order force constants(3rdFCs) at different levels of accuracy. Numerical and analytical 3rdFCs are obtainedvia finite-differencing and symbolic differentiation, respectively.Ref. [128]. In general, a crucially important point is thatVν,ν2,ν3(k → 0+) = 0 if and only if∑l1,p1Φα1,α2,α3l1,p1,l2,p2,l3,p3 = 0. (6.4.2)Here on the right hand side, the sum over third order force constants should be zero for anypair of summation indices (l1, p1), (l2, p2) or (l3, p3), which encodes the invariance of thepotential energy against rigid lattice translations. In the case of CNTs, for low-frequencyTW modes, an additional sum rule has to be obeyed by third order force constants thatderives from an invariance against rigid rotations around the tube axis [122]. Accordingly,from Equations (6.4.1) and (6.4.2), we anticipate that spurious, rapidly decreasing trendsof anharmonic acoustic phonon lifetimes τ3phν ∼ 1/Γ±ν,ν2,ν3 become apparent for k → 0+whenever such sum rules of third order anharmonic force constants are not strictly enforced.For the low-frequency acoustic modes in CNTs, we demonstrate in Figure 6.1 the highsensitivity of numerically calculated dissipation rates to the accuracy of third order forceconstants. Here, we take two different sets of third order force constants as input, both ofwhich are derived from the optimized Tersoff potential [72]. One set is generated by system-atically displacing pairs of atoms a small but finite distance h away from the equilibriumconfiguration of a 1× 1× 3 CNT supercell, so that the resulting atomic forces allow for an77approximation according to whichΦα1,α2,α3l1,p1,l2,p2,l3,p3 ≈14h2[zα1l1,p1(uα2l2,p2 = huα3l3,p3 = −h)−zα1l1,p1(uα2l2,p2 = huα3l3,p3 = h)+ zα1l1,p1(uα2l2,p2 = −huα3l3,p3 = h)−zα1l1,p1(uα2l2,p2 = −huα3l3,p3 = −h)],(6.4.3)where individual force components zαl,p are evaluated by means of analytical formulas. Indeducing the collision matrix C∞, such a finite-differencing approach was also applied byBarbalinardo et al. [59]. Approximating third order force constants as per Equation (6.4.3)is common practice in ab initio PBTT calculations and, in this connection, the importanceof enforcing sum rules a posteriori is often emphasized [2, 68, 69]. Here, however, we aimto reproduce the calculation results of Barbalinardo et al. whose computational approach,as detailed in Ref. [59], did not involve any correction procedures in conjunction with thirdorder force constant calculations as per Equation (6.4.3). With the optimized Tersoff po-tential at hand, whose functional form is continuously differentiable, the other set of thirdorder force constants is produced by explicit symbolic differentiations. In previous chap-ters, this latter approach was pursued throughout since, within the limitations of machineprecision, it ensures optimal compliance with sum rules.As can be noted from Figure 6.1, third order force constants obtained via finite differenc-ing cause spurious predictions of anharmonic acoustic phonon lifetimes. A small violationof sum rules as per Equation (6.4.2) brings with it a nonvanishing coupling coefficientin the long wavelength-limit, Vν,ν2,ν3(k → 0+) 6= 0, and in turn long-wavelength scat-tering rates as per Equation (6.4.1) blow up because contributions of ζν,ν2,ν3 and Vν,ν2,ν3no longer balance appropriately. At the RTA level, one is confronted with grid integra-tions, κRTA ∝ N−1∑σ cσv2στ3phσ , and, according to Figure 6.1, convergence issues linkedto acoustic modes with large lifetimes will be masked whenever calculations are carriedout with third order force constants of low numerical accuracy. In such cases, it is easyto miss the possibly divergent or at least slowly converging conductivity contributions oflong-wavelength acoustic modes, which were ascertained at the RTA level in Chapter 5. InSection 6.4.3, going beyond the RTA, we investigate to what extent these inaccuracies ofthird force constants propagate to predictions of κinv(L→∞).6.4.3 Numerical resultsIn Figure 6.2, we present our κinv(L → ∞) data obtained by inverting the infinite-sizecollision matrix C∞ for varying wave number grids. Results are shown for a (10,0) CNTwhich implies C∞ ∈ R120N×120N , with N (horizontal axis) denoting the discretization of78Figure 6.2: Room temperature thermal conductivity predictions κinv(L → ∞) for a(10,0) CNT as obtained by inverting C∞, cf. Equation (6.3.2), vs wave number gridresolution. Third order anharmonic force constants (3rdFCs) are used at different levelsof numerical accuracy, see Section 6.4.2 for explanations. The left (right) plot showsnumerical data based on quantum (classical) phonon statistics. Calculation results arecompared with limiting values reported in Ref. [59].the 1D wave number grid 1× 1×N . To expand on the importance of numerical accuracyas per Section 6.4.2, we calculate third order force constants entering into the computationof C∞ in two different ways: using a finite differencing approach; and more accurately,via symbolic differentiation. Furthermore, quantum and “classical” calculations are carriedout, meaning we embark on the strategy of replacing actual (quantum) phonon occupationsnumbers in the PBTT formalism with the classical high-temperature limits n0,clσ = kBT/~ωσ,irrespectively of phonon frequency. According to Ref. [59], by means of such a strategy, it ispossible to achieve convincing quantitative agreement between MD and PBTT-derived roomtemperature conductivity data, regardless of whether microscopically or macroscopicallylong free-standing CNTs are considered. Here, for a (10,0) CNT, a classical limiting valueκinv(L→∞) = 3190 Wm−1K−1 was indicated.Based on Figure 6.2, no conclusive answer in regard to the existence of an intrinsicallyupper bounded thermal transport coefficient κinv(L→∞) can be given, provided that cal-culations are performed with high-quality third order force constants. In general, for finerand finer wave number grids, the inverse of an increasingly large collision matrix C∞ hasto be computed, which implies strict upper bounds in terms of feasible calculations dueto limited computational memory. Looking at a (10,0) CNT, we find this upper boundto be given by a wave number grid with N = 400 discrete k-grid points. Therefore, itremains unclear whether or not the blue curve in the left panel of Figure 6.2 eventuallysaturates. Nonetheless, Figure 6.2 clearly demonstrates that a spurious convergent behav-ior of κinv(L → ∞) with respect to the k-grid becomes apparent, erroneously, wheneverthird order force constants are of insufficient numerical accuracy. By resorting to a finite-79differencing approach, as similarly pursued in Ref. [59], we are able to closely reproducelimiting values κinv(L → ∞) of CNT heat transport reported in Ref. [59]. These values,however, transpire to be misleading in both quantum and classical calculations. At the levelof the RTA, a systematic underprediction of CNT heat transfer rates in the long-tube limitdue to a violation of translational and rotational sum rules can be perfectly understoodbased on analytical considerations of low-frequency phonon dissipation rates. By going to ahigher level of theory, a picture is drawn by Figure 6.2 about wave number grid integrationissues that is qualitatively consistent with RTA-level predictions.6.5 ConclusionsAs an alternative to performing iterative PBTT calculations, we have considered a matrixinversion approach for computing the phonon conductivity of CNTs. Based on an accuratedescription of low-frequency anharmonic phonon scattering rates, this new machinery hasthen been applied to re-examine the possibility of tube-length-converged conductivity co-efficients in the case of free-standing CNTs. Despite some conceptual advantages in termsof numerical stability, we find that the debate about the infinite-size CNT thermal conduc-tivity cannot be settled by a brute force matrix inversion approach. At a level of theoryhigher than the RTA, both the reconciliation of PBTT and MD-derived heat transfer ratesof CNTs as well as reliable quantum mechanical estimates of room temperature heat trans-fer rates in macroscopically long CNTs remain open theoretical issues. This contradictsthe recent work of Barbalinardo et al. [59] who, unfortunately, based their conclusions onsystematic calculation errors.One important insight from our work is the high sensitivity of numerical results to theanharmonic force constant accuracy. For lattice-dynamical heat transport calculations onCNTs, we can state: The larger the tube length, the more important it becomes to enforcetranslational and rotational sum rules of anharmonic force constants. In a similar form,such a statement is expected to hold true for 1D materials in general. For CNTs, theanalytically derived scaling laws of low-frequency anharmonic scattering rates prove to bea useful tool by means of which the validity of PBTT calculation results can be checked.80Chapter 7Conclusions and future work7.1 ConclusionsBased on an empirical, realistic force model, we have pursued the two main-stream atomisticmodeling approaches, the PBTT and MD, to study intrinsic heat transport properties ofindividual CNTs. Thermal transport in these materials represents a particular challenge toatomistic modeling approaches. In the framework of the PBTT, issues arise in the light ofa 1D phonon momentum space. In MD simulations of heat transport, it is mainly the weaklattice anharmonicity of CNTs that promotes large time and length scales over which thedynamics of atoms have to be sampled. CNTs have thus long been a matter of controversyin the field of nanoscale thermal transport. The present work gives new insights into thesechallenges, marking a significant contribution to our understanding of low-dimensional heattransport in a technologically important class of materials.In Chapter 4, we compared PBTT and NEMD-derived thermal conductivity predic-tions as a function of system size and temperature. This comparison allowed us to criticallyevaluate the opposing approximations in the different approaches, namely the perturba-tive treatment of phonon-phonon interactions in the PBTT and the use of classical phononstatistics in MD. Here, we observed no signatures of beyond-lowest-order phonon-phononinteractions for the tube length scales L . 5µm and the temperature range T . 800 K ac-cessible to NEMD simulations. This naturally raises the question: Should one rely fully onquantum mechanical PBTT calculations of CNT heat transport, which in fact give accessto larger length scales and more extreme temperature ranges? Certainly, in view of possiblehigher order phonon-phonon interaction, more work is needed to establish the applicabilityof perturbative calculations in the case of CNTs at extremely high temperatures T > 1000 K(see discussion below). Perhaps more interestingly, as was further shown in Chapter 4, thereliability of the PBTT is challenged by large tube length scales as well. Pursuing an iter-81ative scheme to solve the linearized Peierls-Boltzmann transport equation, we encounteredthe problem of weakly damped acoustic phonons in CNTs. This problem directly links tothe long-standing debate over thermal conductivity convergence vs divergence in CNTs.From our numerics in the framework of the PBTT, we concluded the following: (i) Aniterative solution approach to the Peierls-Boltzmann transport equation with three-phononinteractions yields a tube-length divergence of the CNT thermal conductivity, suggestinganomalous (i.e. superdiffusive) heat transport in CNTs. (ii) The longer the tube length,the harder it becomes to numerically pinpoint the conductivity contributions of acousticphonons in CNTs. As a direct consequence of (ii), for reliable estimates of the CNT thermalconductivity, we recommend standard iterative PBTT calculations only up to tube lengthscales of roughly L . 100µm.In Chapter 5, by applying the RTA in the framework of the PBTT, we clarified thelarge-L behavior of the CNT thermal conductivity and related numerical issues stemmingfrom weakly damped acoustic phonons. Most notably, by taking input from continuumtheories and by the help of analytical considerations, we rigorously proved that, at macro-scopic tube length scales, one cannot converge RTA-level conductivities of CNTs by meansof standard wave number grid sampling approaches, with maximal k-grid resolutions (typ-ically ∆k a ∼ 10−3) being limited by computational memory. Furthermore, as far as thelarge-L behavior of thermal transport is concerned, we demonstrated the decisive role oftensile lattice strain. In going from a zero to a small, but finite, strain amplitude , we foundthat the anharmonic three-phonon scattering strength of acoustic modes changes discontin-uously in the long-wavelength limit. Attributed to a reduced phonon momentum space di-mensionality, this somewhat pathological strain behavior carries over to the predicted CNTthermal conductivity to the extent that, in the large-L limit, we have κRTA( = 0) <∞ butκRTA( > 0) → ∞. By arguing that strain cannot be avoided in experimental setups, onemay interpret these findings as proof for superdiffusive heat transport in ultra-long CNTs.However, as we also demonstrated in Chapter 5, additional discontinuities in the strain de-pendence of κRTA are encountered as L becomes large. Such discontinuities in the large-Llimit, which are rooted in the 1D phonon momentum space of CNTs, cast some conceptualdoubts on the applicability of standard PBTT calculations when applied to macroscopicallylong CNTs. This further substantiates our conclusion from Chapter 4 that the PBTT, in-cluding three-phonon interactions, can be considered as a reliable tool to quantify the CNTthermal conductivity only for microscopic rather than macroscopic tube length scales.Although our results in Chapter 5 were confined to the RTA, we subsequently showed inChapter 6 that the debate over convergence vs divergence of the CNT thermal conductivityis not easily resolved if one goes to a higher level of transport theory, as was recently claimedin Ref. [59]. In particular, by pursuing a matrix inversion approach as an alternative to an82iterative scheme for solving the linearized Peierls-Boltzmann transport equation, we againarrived at the conclusion that reliable quantum mechanical estimates of the thermal con-ductivity cannot be given for macroscopically long CNTs, that is - they are beyond thepredictive capabilities of modern-day atomistic modeling approaches.7.2 Future workThere are some aspects of phonon-mediated heat transport in CNTs that still remain poorlyunderstood. In the following, we point out a few open questions and suggest future researchdirections.• How do results from the Green-Kubo relation fit in?Among atomistic modeling approaches and independent of the PBTT as well as the NEMDmethod, the Green-Kubo (GK) relation of the thermal conductivity [143, 144] when com-bined with EMD simulations provides an alternative route toward lattice heat transport,not pursued in the present work. Most recent studies in the area [57, 59] seem to indi-cate that, by employing a GK-EMD approach, a finite but classical value of the thermalconductivity for macroscopically long CNTs can be deduced. If this is the case, then oneshould expect that quantum upper bounds exist as well and that a more refined quantumtheory for lattice transport in CNTs may be in demand. Indeed, in the case of quasi-1Dmaterials, low-frequency phonon-phonon interactions may require a nonperturbative treat-ment [132], as is also indicated by our somewhat ill-behaved EMD-derived normal modefrequency spectra in Section 4.1.4. Moreover, it would be interesting to use the GK-EMDmethod to examine the diameter dependence of the CNT thermal conductivity. As we sawin Section 4.1, the qualitative disagreement between PBTT and NEMD-derived diameterdependencies constitutes an unsolved problem.• When do higher-order phonon-phonon interactions become significant?To fully answer this question one has to resort to an ab initio modeling approach because thesignificance of higher-order phonon-phonon interactions is ultimately determined by the trueshape of the Born-Oppenheimer energy landscape. With sufficient computational resources,an ab initio MD approach may give the best possible description of the high-temperatureCNT thermal conductivity, indicating possible deviations from the T−1 behavior. Alter-natively, one may try to go beyond the three-phonon limitation in the framework of thePBTT. In fact, computer codes to tackle four-phonon interactions in the PBTT are becom-ing more prevalent [145, 146]. For different temperatures, such calculations would enable83direct comparisons of three-phonon and four-phonon scattering rates in CNTs. Moreover,some previous PBTT studies [53, 56] used four-phonon interactions in a rather ad hoc man-ner to argue that the PBTT-derived CNT thermal conductivity should remain finite in thelimit of large tube lengths. In a way similar to Chapter 5 and perhaps with the help of anempirical force model, it would be interesting to scrutinize such claims in a more systematicmanner.• How sensitive is the CNT thermal conductivity to external perturbations?In the case of truly 1D monoatomic chain models, it is often observed that superdiffusiveheat transport becomes diffusive whenever the translational invariance of the chain is bro-ken [147, 148], which is achieved through some kind of on-site pinning potential. In thecase of real-world CNTs, such an external pinning potential may be represented by a sub-strate onto which the tube is deposited. Hence, by means of NEMD simulations, it willbe intriguing to check whether the length dependence of the CNT thermal conductivity in,say, micrometer-long tubes changes if one attaches a spring to initial atomic positions.• Are electron-phonon interactions important?In a nonequilibrium heat transport scenario, phonons may also transfer momentum andenergy to electrons and vice versa. These effects can be handled by a Boltzmann transportequation that couples electronic and phononic degrees of freedom [149]. The question howsignificant these mutual interactions in heat conducting CNTs are awaits future investiga-tions.84Bibliography[1] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S.Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J.Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, “LAMMPS - a flexible simula-tion tool for particle-based materials modeling at the atomic, meso, and continuumscales”, Computer Physics Communications 271, 108171 (2022) (pg. v).[2] W. Li, J. Carrete, N. A. Katcho, and N. Mingo, “ShengBTE: a solver of the Boltz-mann transport equation for phonons”, Computer Physics Communications 185,1747–1758 (2014) (pg. v, 29, 78).[3] G. Barbalinardo, Z. Chen, N. W. Lundgren, and D. Donadio, “Efficient anhar-monic lattice dynamics calculations of thermal transport in crystalline and disorderedsolids”, Journal of Applied Physics 128, 135104 (2020) (pg. v).[4] Q. Zheng, M. Hao, R. Miao, J. Schaadt, and C. Dames, “Advances in thermal con-ductivity for energy applications: a review”, Prog. Energy 3, 012002 (2021) (pg. 1).[5] A. A. Balandin, “Nanophononics: phonon engineering in nanostructures and nanode-vices”, Journal of Nanoscience and Nanotechnology 5, 1015–1022 (2005) (pg. 1).[6] D. G. Cahill, P. V. Braun, G. Chen, D. R. Clarke, S. Fan, K. E. Goodson, P. Keblinski,W. P. King, G. D. Mahan, A. Majumdar, H. J. Maris, S. R. Phillpot, E. Pop, andL. Shi, “Nanoscale thermal transport. II. 2003–2012”, Applied Physics Reviews 1,011305 (2014) (pg. 1).[7] A. M. Marconnet, M. A. Panzer, and K. E. Goodson, “Thermal conduction phenom-ena in carbon nanotubes and related nanostructured materials”, Rev. Mod. Phys.85, 1295–1326 (2013) (pg. 1).[8] B. Kumanek and D. Janas, “Thermal conductivity of carbon nanotube networks: areview”, J. Mater. Sci. 54, 7397–7427 (2019) (pg. 1).[9] X. Gui, A. Cao, J. Wei, H. Li, Y. Jia, Z. Li, L. Fan, K. Wang, H. Zhu, and D. Wu,“Soft, highly conductive nanotube sponges and composites with controlled compress-ibility”, ACS Nano 4, 2320–2326 (2010) (pg. 1).85[10] R. L. Whitby, T. Fukuda, T. Maekawa, S. L. James, and S. V. Mikhalovsky, “Ge-ometric control and tuneable pore size distribution of buckypaper and buckydiscs”,Carbon 46, 949–956 (2008) (pg. 1).[11] N. Behabtu, C. C. Young, D. E. Tsentalovich, O. Kleinerman, X. Wang, A. W. K. Ma,E. A. Bengio, R. F. ter Waarbeek, J. J. de Jong, R. E. Hoogerwerf, S. B. Fairchild,J. B. Ferguson, B. Maruyama, J. Kono, Y. Talmon, Y. Cohen, M. J. Otto, and M.Pasquali, “Strong, light, multifunctional fibers of carbon nanotubes with ultrahighconductivity”, Science 339, 182–186 (2013) (pg. 1).[12] M. Xu, D. N. Futaba, M. Yumura, and K. Hata, “Alignment control of carbon nan-otube forest from random to nearly perfectly aligned by utilizing the crowding effect”,ACS Nano 6, 5837–5844 (2012) (pg. 1).[13] R. S. Prasher, X. J. Hu, Y. Chalopin, N. Mingo, K. Lofgreen, S. Volz, F. Cleri,and P. Keblinski, “Turning carbon nanotubes from exceptional heat conductors intoinsulators”, Phys. Rev. Lett. 102, 105901 (2009) (pg. 1).[14] L. Zhang, G. Zhang, C. Liu, and S. Fan, “High-density carbon nanotube buckypaperswith superior transport and mechanical properties”, Nano Letters 12, 4848–4852(2012) (pg. 1).[15] J. L. Blackburn, A. J. Ferguson, C. Cho, and J. C. Grunlan, “Carbon-nanotube-basedthermoelectric materials and devices”, Adv. Mater. 30, 1704386 (2018) (pg. 1).[16] W. Yu, C. Liu, and S. Fan, “Advances of CNT-based systems in thermal manage-ment”, Nano Res. 14, 2471–2490 (2021) (pg. 1).[17] E. Pop, V. Varshney, and A. K. Roy, “Thermal properties of graphene: fundamentalsand applications”, MRS Bulletin 37, 1273–1281 (2012) (pg. 2).[18] A. A. Balandin, “Phononics of graphene and related materials”, ACS Nano 14, 5170–5178 (2020) (pg. 2).[19] S. Berber, Y.-K. Kwon, and D. Toma´nek, “Unusually high thermal conductivity ofcarbon nanotubes”, Phys. Rev. Lett. 84, 4613–4616 (2000) (pg. 2, 7).[20] S Lepri, R Livi, and A Politi, “On the anomalous thermal conductivity of one-dimensional lattices”, Europhysics Letters (EPL) 43, 271–276 (1998) (pg. 2, 11).[21] O. Narayan and S. Ramaswamy, “Anomalous heat conduction in one-dimensionalmomentum-conserving systems”, Phys. Rev. Lett. 89, 200601 (2002) (pg. 2, 11).[22] A. Pereverzev, “Fermi-Pasta-Ulam β lattice: Peierls equation and anomalous heatconductivity”, Phys. Rev. E 68, 056124 (2003) (pg. 2, 11).86[23] J. Lukkarinen and H. Spohn, “Anomalous energy transport in the FPU-β chain”,Communications on Pure and Applied Mathematics 61, 1753–1786 (2008) (pg. 2,11).[24] T. Mai, A. Dhar, and O. Narayan, “Equilibration and universal heat conduction inFermi-Pasta-Ulam chains”, Phys. Rev. Lett. 98, 184301 (2007) (pg. 2, 11).[25] K. Saito and A. Dhar, “Heat conduction in a three dimensional anharmonic crystal”,Phys. Rev. Lett. 104, 040601 (2010) (pg. 2, 11).[26] R. Livi and S. Lepri, “Heat in one dimension”, Nature 421, 327–327 (2003) (pg. 2).[27] A. Dhar, “Heat transport in low-dimensional systems”, Advances in Physics 57, 457–537 (2008) (pg. 2).[28] A. Dhar and K. Saito, “Heat transport in harmonic systems”, in Thermal transportin low dimensions: from statistical physics to nanoscale heat transfer , edited by S.Lepri (Springer International Publishing, Cham, 2016), pp. 39–105 (pg. 3).[29] G. Fugallo and L. Colombo, “Calculating lattice thermal conductivity: a synopsis”,Phys. Scr. 93, 043002 (2018) (pg. 3).[30] R. Saito, G. Dresselhaus, and M. Dresselhaus, Physical properties of carbon nanotubes(Imperial College Press, London, 1998) (pg. 5).[31] S. Reich, C. Thomsen, and J. Maultzsch, Carbon nanotubes: basic concepts and phys-ical properties (Wiley-VCH, Weinheim, 2004) (pg. 6, 22).[32] J. Hone, M. Whitney, C. Piskoti, and A. Zettl, “Thermal conductivity of single-walledcarbon nanotubes”, Phys. Rev. B 59, R2514–R2516 (1999) (pg. 8, 11).[33] J. Hone, M. C. Llaguno, N. M. Nemes, A. T. Johnson, J. E. Fischer, D. A. Wal-ters, M. J. Casavant, J. Schmidt, and R. E. Smalley, “Electrical and thermal trans-port properties of magnetically aligned single wall carbon nanotube films”, AppliedPhysics Letters 77, 666–668 (2000) (pg. 8).[34] F. Lian, J. P. Llinas, Z. Li, D. Estrada, and E. Pop, “Thermal conductivity ofchirality-sorted carbon nanotube networks”, Applied Physics Letters 108, 103101(2016) (pg. 8).[35] T. Yamamoto, S. Watanabe, and K. Watanabe, “Universal features of quantized ther-mal conductance of carbon nanotubes”, Phys. Rev. Lett. 92, 075502 (2004) (pg. 8).[36] N. Mingo and D. A. Broido, “Carbon nanotube ballistic thermal conductance andits limits”, Phys. Rev. Lett. 95, 096105 (2005) (pg. 8, 12, 24, 25).87[37] C. Yu, L. Shi, Z. Yao, D. Li, and A. Majumdar, “Thermal conductance and ther-mopower of an individual single-wall carbon nanotube”, Nano Letters 5, 1842–1846(2005) (pg. 8, 9).[38] E. Pop, D. Mann, Q. Wang, K. Goodson, and H. Dai, “Thermal conductance of anindividual single-wall carbon nanotube above room temperature”, Nano Letters 6,96–100 (2006) (pg. 8, 9, 11, 13).[39] J. Liu, T. Li, Y. Hu, and X. Zhang, “Benchmark study of the length dependent ther-mal conductivity of individual suspended, pristine SWCNTs”, Nanoscale 9, 1496–1501 (2017) (pg. 8, 9).[40] V. Lee, C.-H. Wu, Z.-X. Lou, W.-L. Lee, and C.-W. Chang, “Divergent and ultrahighthermal conductivity in millimeter-long nanotubes”, Phys. Rev. Lett. 118, 135901(2017) (pg. 9).[41] J. Tersoff, “Empirical interatomic potential for carbon, with applications to amor-phous carbon”, Phys. Rev. Lett. 61, 2879–2882 (1988) (pg. 9).[42] D. W. Brenner, “Empirical potential for hydrocarbons for use in simulating thechemical vapor deposition of diamond films”, Phys. Rev. B 42, 9458–9471 (1990)(pg. 9).[43] D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, and S. B. Sin-nott, “A second-generation reactive empirical bond order (REBO) potential energyexpression for hydrocarbons”, Journal of Physics: Condensed Matter 14, 783–802(2002) (pg. 9).[44] S. J. Stuart, A. B. Tutein, and J. A. Harrison, “A reactive potential for hydrocarbonswith intermolecular interactions”, The Journal of Chemical Physics 112, 6472–6486(2000) (pg. 9).[45] R. N. Salaway and L. V. Zhigilei, “Molecular dynamics simulations of thermal con-ductivity of carbon nanotubes: resolving the effects of computational parameters”,International Journal of Heat and Mass Transfer 70, 954–964 (2014) (pg. 9).[46] J. A. Thomas, R. M. Iutzi, and A. J. H. McGaughey, “Thermal conductivity andphonon transport in empty and water-filled carbon nanotubes”, Phys. Rev. B 81,045413 (2010) (pg. 9–12).[47] S.-Y. Yue, T. Ouyang, and M. Hu, “Diameter dependence of lattice thermal conduc-tivity of single-walled carbon nanotubes: study from ab initio”, Sci. Rep. 5, 15440(2015) (pg. 9–12, 43, 56).88[48] L. Lindsay, D. A. Broido, and N. Mingo, “Diameter dependence of carbon nanotubethermal conductivity and extension to the graphene limit”, Phys. Rev. B 82, 161402(2010) (pg. 9, 10, 41, 43, 54, 56, 69).[49] A. Cao and J. Qu, “Size dependent thermal conductivity of single-walled carbonnanotubes”, Journal of Applied Physics 112, 013503 (2012) (pg. 9–12, 41, 43, 44,54).[50] J. Ziman, Electrons and Phonons (Oxford University Press, Oxford, 1960) (pg. 10,11, 17, 23, 24, 26, 27, 45, 46, 53, 59, 70, 74, 104).[51] N. W. Ashcroft and N. D. Mermin, Solid state physics (Saunders College Publishing,Fort Worth, 1976) (pg. 10, 11, 13, 53).[52] S. Maruyama, “A molecular dynamics simulation of heat conduction of a finite lengthsingle-walled carbon nanotube”, Microscale Thermophysical Engineering 7, 41–50(2003) (pg. 11, 12).[53] N. Mingo and D. A. Broido, “Length dependence of carbon nanotube thermal conduc-tivity and the “problem of long waves””, Nano Letters 5, 1221–1225 (2005) (pg. 11–13, 17, 24, 25, 46, 53, 56, 84).[54] G. Zhang and B. Li, “Thermal conductivity of nanotubes revisited: effects of chirality,isotope impurity, tube length, and temperature”, The Journal of Chemical Physics123, 114714 (2005) (pg. 11, 12).[55] D. Donadio and G. Galli, “Thermal conductivity of isolated and interacting carbonnanotubes: comparing results from molecular dynamics and the Boltzmann transportequation”, Phys. Rev. Lett. 99, 255502 (2007) (pg. 11, 12).[56] L. Lindsay, D. A. Broido, and N. Mingo, “Lattice thermal conductivity of single-walled carbon nanotubes: beyond the relaxation time approximation and phonon-phonon scattering selection rules”, Phys. Rev. B 80, 125407 (2009) (pg. 11, 12, 43,46, 48, 56, 59, 69, 84).[57] Z. Fan, H. Dong, A. Harju, and T. Ala-Nissila, “Homogeneous nonequilibrium molec-ular dynamics method for heat transport and spectral decomposition with many-body potentials”, Phys. Rev. B 99, 064308 (2019) (pg. 11, 12, 83).[58] U. Ray and D. T. Limmer, “Heat current fluctuations and anomalous transport inlow-dimensional carbon lattices”, Phys. Rev. B 100, 241409 (2019) (pg. 11, 12).[59] G. Barbalinardo, Z. Chen, H. Dong, Z. Fan, and D. Donadio, “Ultrahigh convergentthermal conductivity of carbon nanotubes from comprehensive atomistic modeling”,Phys. Rev. Lett. 127, 025902 (2021) (pg. 11, 12, 16, 56, 70, 72, 73, 75, 76, 78–80, 82,83).89[60] L. Lindsay and D. A. Broido, “Theory of thermal transport in multilayer hexagonalboron nitride and nanotubes”, Phys. Rev. B 85, 035436 (2012) (pg. 13, 53, 56).[61] J. R. Lukes and H. Zhong, “Thermal conductivity of individual single-wall carbonnanotubes”, Journal of Heat Transfer 129, 705–716 (2006) (pg. 13).[62] J. E. Turney, A. J. H. McGaughey, and C. H. Amon, “Assessing the applicabilityof quantum corrections to classical thermal conductivity predictions”, Phys. Rev. B79, 224305 (2009) (pg. 14).[63] L. Lindsay, A. Katre, A. Cepellotti, and N. Mingo, “Perspective on ab initio phononthermal transport”, Journal of Applied Physics 126, 050902 (2019) (pg. 14).[64] S. Stackhouse, L. Stixrude, and B. B. Karki, “Thermal conductivity of periclase(MgO) from first principles”, Phys. Rev. Lett. 104, 208501 (2010) (pg. 14).[65] A. Marcolongo, P. Umari, and S. Baroni, “Microscopic theory and quantum simula-tion of atomic heat transport”, Nature Physics 12, 80–84 (2016) (pg. 14).[66] C. Carbogno, R. Ramprasad, and M. Scheffler, “Ab initio Green-Kubo approach forthe thermal conductivity of solids”, Phys. Rev. Lett. 118, 175901 (2017) (pg. 14).[67] M. Puligheddu, F. Gygi, and G. Galli, “First-principles simulations of heat trans-port”, Phys. Rev. Materials 1, 060802 (2017) (pg. 14).[68] A. J. H. McGaughey, A. Jain, H.-Y. Kim, and B. Fu, “Phonon properties and thermalconductivity from first principles, lattice dynamics, and the Boltzmann transportequation”, Journal of Applied Physics 125, 011101 (2019) (pg. 14, 78).[69] W. Li, L. Lindsay, D. A. Broido, D. A. Stewart, and N. Mingo, “Thermal conductivityof bulk and nanowire Mg2SixSn1−x alloys from first principles”, Phys. Rev. B 86,174307 (2012) (pg. 14, 15, 78).[70] J. Tersoff, “New empirical approach for the structure and energy of covalent systems”,Phys. Rev. B 37, 6991–7000 (1988) (pg. 14, 97).[71] J. Tersoff, “Empirical interatomic potential for carbon, with applications to amor-phous carbon”, Phys. Rev. Lett. 61, 2879–2882 (1988) (pg. 14).[72] L. Lindsay and D. A. Broido, “Optimized Tersoff and Brenner empirical potentialparameters for lattice dynamics and phonon thermal transport in carbon nanotubesand graphene”, Phys. Rev. B 81, 205441 (2010) (pg. 14, 15, 43, 75, 77, 98).[73] M. Mohr, J. Maultzsch, E. Dobardzic, S. Reich, I. Milosevic, M. Damnjanovic, A.Bosak, M. Krisch, and C. Thomsen, “Phonon dispersion of graphite by inelastic x-rayscattering”, Phys. Rev. B 76, 035439 (2007) (pg. 14).90[74] T. D. Nguyen, “GPU-accelerated Tersoff potentials for massively parallel molecu-lar dynamics simulations”, Computer Physics Communications 212, 113–122 (2017)(pg. 15).[75] M. Born and R. Oppenheimer, “Zur Quantentheorie der Molekeln”, Annalen derPhysik 389, 457–484 (1927) (pg. 15).[76] G. D. Mahan and F. Claro, “Nonlocal theory of thermal conductivity”, Phys. Rev.B 38, 1963–1969 (1988) (pg. 15).[77] P. B. Allen, “Size effects in thermal conduction by phonons”, Phys. Rev. B 90, 054301(2014) (pg. 15).[78] P. B. Allen, “Analysis of nonlocal phonon thermal conductivity simulations showingthe ballistic to diffusive crossover”, Phys. Rev. B 97, 134307 (2018) (pg. 15).[79] A. Cepellotti and N. Marzari, “Boltzmann transport in nanostructures as a frictioneffect”, Nano Letters 17, 4675–4682 (2017) (pg. 15, 25).[80] C. Hua and L. Lindsay, “Space-time dependent thermal conductivity in nonlocalthermal transport”, Phys. Rev. B 102, 104310 (2020) (pg. 15, 25).[81] V. Chiloyan, S. Huberman, Z. Ding, J. Mendoza, A. A. Maznev, K. A. Nelson, and G.Chen, “Green’s functions of the Boltzmann transport equation with the full scatteringmatrix for phonon nanoscale transport beyond the relaxation-time approximation”,Phys. Rev. B 104, 245424 (2021) (pg. 15, 25).[82] N. Bonini, J. Garg, and N. Marzari, “Acoustic phonon lifetimes and thermal trans-port in free-standing and strained graphene”, Nano Letters 12, 2673–2678 (2012)(pg. 16, 55, 56, 70).[83] R. Peierls, Quantum theory of solids (Clarendon Press, Oxford, 1955) (pg. 17, 23,24, 27, 74, 76).[84] G. Srivastava, The physics of phonons (Taylor & Francis, New York, 1990) (pg. 17,24, 27, 29).[85] M. T. Dove, Introduction to lattice dynamics (Cambridge University Press, Cam-bridge, 1993) (pg. 22, 36).[86] R. Juneja, S. The´baud, T. Pandey, C. Polanco, D. Moseley, M. Manley, Y. Cheng,B. Winn, D. Abernathy, R. Hermann, and L. Lindsay, “Quasiparticle twist dynamicsin non-symmorphic materials”, Materials Today Physics 21, 100548 (2021) (pg. 22).[87] R. Peierls, “Zur kinetischen Theorie der Wa¨rmeleitung in Kristallen”, Annalen derPhysik 395, 1055–1101 (1929) (pg. 23, 31).91[88] L. G. C. Rego and G. Kirczenow, “Quantized thermal conductance of dielectricquantum wires”, Phys. Rev. Lett. 81, 232–235 (1998) (pg. 26).[89] J.-S. Wang, J. Wang, and J. T. Lu¨, “Quantum thermal transport in nanostructures”,The European Physical Journal B 62, 381–404 (2008) (pg. 26).[90] T. Feng and X. Ruan, “Quantum mechanical prediction of four-phonon scatteringrates and reduced thermal conductivity of solids”, Phys. Rev. B 93, 045202 (2016)(pg. 27, 42).[91] M. Omini and A. Sparavigna, “An iterative approach to the phonon Boltzmannequation in the theory of thermal conductivity”, Physica B: Condensed Matter 212,101–112 (1995) (pg. 29, 30).[92] D. A. Broido, M. Malorny, G. Birner, N. Mingo, and D. A. Stewart, “Intrinsic lat-tice thermal conductivity of semiconductors from first principles”, Applied PhysicsLetters 91, 231922 (2007) (pg. 29).[93] A. A. Maradudin and A. E. Fein, “Scattering of neutrons by an anharmonic crystal”,Phys. Rev. 128, 2589–2608 (1962) (pg. 30, 35, 36).[94] Z. Li, S. Xiong, C. Sievers, Y. Hu, Z. Fan, N. Wei, H. Bao, S. Chen, D. Donadio, andT. Ala-Nissila, “Influence of thermostatting on nonequilibrium molecular dynamicssimulations of heat conduction in solids”, The Journal of Chemical Physics 151,234105 (2019) (pg. 32, 33).[95] J. A. Thomas, J. E. Turney, R. M. Iutzi, C. H. Amon, and A. J. H. McGaughey, “Pre-dicting phonon dispersion relations and lifetimes from the spectral energy density”,Phys. Rev. B 81, 081411 (2010) (pg. 32, 35, 36).[96] Z. Fan, L. F. C. Pereira, H.-Q. Wang, J.-C. Zheng, D. Donadio, and A. Harju,“Force and heat current formulas for many-body potentials in molecular dynamicssimulations with applications to thermal conductivity calculations”, Phys. Rev. B92, 094301 (2015) (pg. 33).[97] D. Surblys, H. Matsubara, G. Kikugawa, and T. Ohara, “Application of atomicstress to compute heat flux via molecular dynamics for systems with many-bodyinteractions”, Phys. Rev. E 99, 051301 (2019) (pg. 33).[98] P. Boone, H. Babaei, and C. E. Wilmer, “Heat flux for many-body interactions:corrections to LAMMPS”, Journal of Chemical Theory and Computation 15, 5579–5587 (2019) (pg. 33).[99] T. Ikeshoji and B. Hafskjold, “Non-equilibrium molecular dynamics calculation ofheat conduction in liquid and through liquid-gas interface”, Molecular Physics 81,251–261 (1994) (pg. 33).92[100] F. Mu¨ller-Plathe, “A simple nonequilibrium molecular dynamics method for calcu-lating the thermal conductivity”, The Journal of Chemical Physics 106, 6082–6085(1997) (pg. 33, 34).[101] P. Jund and R. Jullien, “Molecular-dynamics calculation of the thermal conductivityof vitreous silica”, Phys. Rev. B 59, 13707–13711 (1999) (pg. 33).[102] P. Wirnsberger, D. Frenkel, and C. Dellago, “An enhanced version of the heat ex-change algorithm with excellent energy conservation properties”, The Journal ofChemical Physics 143, 124104 (2015) (pg. 33).[103] T. H. K. Barron and M. L. Klein, “Perturbation theory of anharmonic crystals”,in Dynamical properties of solids, edited by G. K. Horton and A. A. Maradudin(North-Holland Publishing Company, Amsterdam, 1974), pp. 391–449 (pg. 35, 36).[104] T. Sun, X. Shen, and P. B. Allen, “Phonon quasiparticles and anharmonic pertur-bation theory tested by molecular dynamics on a model system”, Phys. Rev. B 82,224304 (2010) (pg. 35, 36).[105] A. J. C. Ladd, B. Moran, and W. G. Hoover, “Lattice thermal conductivity: a com-parison of molecular dynamics and anharmonic lattice dynamics”, Phys. Rev. B 34,5058–5064 (1986) (pg. 35).[106] A. J. H. McGaughey and M. Kaviany, “Quantitative validation of the Boltzmanntransport equation phonon thermal conductivity model under the single-mode relax-ation time approximation”, Phys. Rev. B 69, 094303 (2004) (pg. 35).[107] N. de Koker, “Thermal conductivity of MgO periclase from equilibrium first princi-ples molecular dynamics”, Phys. Rev. Lett. 103, 125902 (2009) (pg. 35).[108] J. E. Turney, E. S. Landry, A. J. H. McGaughey, and C. H. Amon, “Predictingphonon properties and thermal conductivity from anharmonic lattice dynamics cal-culations and molecular dynamics simulations”, Phys. Rev. B 79, 064301 (2009)(pg. 35, 42).[109] T. Feng, B. Qiu, and X. Ruan, “Anharmonicity and necessity of phonon eigenvectorsin the phonon normal mode analysis”, Journal of Applied Physics 117, 195102 (2015)(pg. 35).[110] K. Sa¨a¨skilahti, J. Oksanen, S. Volz, and J. Tulkki, “Frequency-dependent phononmean free path in carbon nanotubes from nonequilibrium molecular dynamics”, Phys.Rev. B 91, 115426 (2015) (pg. 41, 44).[111] Y. He, I. Savic´, D. Donadio, and G. Galli, “Lattice thermal conductivity of semi-conducting bulk materials: atomistic simulations”, Phys. Chem. Chem. Phys. 14,16209–16222 (2012) (pg. 42).93[112] M. Puligheddu, Y. Xia, M. Chan, and G. Galli, “Computational prediction of latticethermal conductivity: a comparison of molecular dynamics and Boltzmann transportapproaches”, Phys. Rev. Materials 3, 085401 (2019) (pg. 42).[113] P. Chantrenne and J.-L. Barrat, “Analytical model for the thermal conductivity ofnanostructures”, Superlattices and Microstructures 35, 173–186 (2004) (pg. 46).[114] J. Wang and J.-S. Wang, “Carbon nanotube thermal transport: ballistic to diffusive”,Applied Physics Letters 88, 111909 (2006) (pg. 46).[115] P. Klemens, “Theory of the thermal conductivity of solids”, in Thermal conductivity,edited by R. Tye (Academic Press, London, 1969), pp. 1–68 (pg. 46).[116] P. Klemens and D. Pedraza, “Thermal conductivity of graphite in the basal plane”,Carbon 32, 735–741 (1994) (pg. 46).[117] Z.-Y. Ong, E. Pop, and J. Shiomi, “Reduction of phonon lifetimes and thermal con-ductivity of a carbon nanotube on amorphous silica”, Phys. Rev. B 84, 165418 (2011)(pg. 51).[118] N. Mounet and N. Marzari, “First-principles determination of the structural, vibra-tional and thermodynamic properties of diamond, graphite, and derivatives”, Phys.Rev. B 71, 205214 (2005) (pg. 55, 68).[119] A. Aghaei, K. Dayal, and R. S. Elliott, “Anomalous phonon behavior of carbonnanotubes: first-order influence of external load”, Journal of Applied Physics 113,023503 (2013) (pg. 55).[120] G. Fugallo, A. Cepellotti, L. Paulatto, M. Lazzeri, N. Marzari, and F. Mauri, “Ther-mal conductivity of graphene and graphite: collective excitations and mean freepaths”, Nano Letters 14, 6109–6114 (2014) (pg. 55, 70).[121] A. Cepellotti and N. Marzari, “Thermal transport in crystals as a kinetic theory ofrelaxons”, Phys. Rev. X 6, 041013 (2016) (pg. 55, 70).[122] V. N. Popov, V. E. Van Doren, and M. Balkanski, “Elastic properties of single-walledcarbon nanotubes”, Phys. Rev. B 61, 3078–3084 (2000) (pg. 56, 57, 77).[123] H. Suzuura and T. Ando, “Phonons and electron-phonon scattering in carbon nan-otubes”, Phys. Rev. B 65, 235412 (2002) (pg. 56, 57).[124] A. De Martino, R. Egger, and A. O. Gogolin, “Phonon-phonon interactions andphonon damping in carbon nanotubes”, Phys. Rev. B 79, 205408 (2009) (pg. 56, 58,59, 62, 104).94[125] P. Gupta and A. Kumar, “Phonons in chiral nanorods and nanotubes: a Cosserat-rod-based continuum approach”, Mathematics and Mechanics of Solids 24, 3897–3919 (2019) (pg. 56, 57).[126] E. Viola and A. Marzani, “Exact analysis of wave motions in rods and hollow cylin-ders”, in Mechanical vibration: where do we stand? , edited by I. Elishakoff (Springer,Vienna, 2007), pp. 83–104 (pg. 57).[127] L. Landau and E. Lifshitz, Theory of elasticity, Vol. 7, Course of Theoretical Physics(Pergamon, Oxford, 1970) Chap. III, pp. 101–122 (pg. 58).[128] E. Lifshitz and L. Pitaevski, Physical kinetics, Vol. 10, Course of Theoretical Physics(Pergamon, Oxford, 1981) Chap. VII, pp. 284–311 (pg. 59, 62, 64, 77, 104, 106).[129] Z. Rashid, L. Zhu, and W. Li, “Effect of confinement on anharmonic phonon scatter-ing and thermal conductivity in pristine silicon nanowires”, Phys. Rev. B 97, 075441(2018) (pg. 64).[130] C. H. Xu, C. Z. Wang, C. T. Chan, and K. M. Ho, “Theory of the thermal expansionof Si and diamond”, Phys. Rev. B 43, 5024–5027 (1991) (pg. 68).[131] K. H. Michel, S. Costamagna, and F. M. Peeters, “Theory of anharmonic phononsin two-dimensional crystals”, Phys. Rev. B 91, 134302 (2015) (pg. 68).[132] Santhosh G. and D. Kumar, “Anomalous transport and phonon renormalization in achain with transverse and longitudinal vibrations”, Phys. Rev. E 82, 011130 (2010)(pg. 69, 83).[133] X. Wang, M. Kaviany, and B. Huang, “Phonon coupling and transport in individualpolyethylene chains: a comparison study with the bulk crystal”, Nanoscale 9, 18022–18031 (2017) (pg. 70).[134] T. Pandey, C. A. Polanco, V. R. Cooper, D. S. Parker, and L. Lindsay, “Symmetry-driven phonon chirality and transport in one-dimensional and bulk Ba3N-derivedmaterials”, Phys. Rev. B 98, 241405 (2018) (pg. 70).[135] L. Lindsay, D. A. Broido, and N. Mingo, “Flexural phonons and thermal transportin graphene”, Phys. Rev. B 82, 115427 (2010) (pg. 70).[136] A. Cepellotti, G. Fugallo, L. Paulatto, M. Lazzeri, F. Mauri, and N. Marzari, “Phononhydrodynamics in two-dimensional materials”, Nature Communications 6, 6400 (2015)(pg. 70).[137] P. B. Allen, “Improved Callaway model for lattice thermal conductivity”, Phys. Rev.B 88, 144302 (2013) (pg. 70, 74).95[138] A. Chernatynskiy and S. R. Phillpot, “Evaluation of computational techniques forsolving the Boltzmann transport equation for lattice thermal conductivity calcula-tions”, Phys. Rev. B 82, 134301 (2010) (pg. 70, 74).[139] G. Fugallo, M. Lazzeri, L. Paulatto, and F. Mauri, “Ab initio variational approach forevaluating lattice thermal conductivity”, Phys. Rev. B 88, 045430 (2013) (pg. 75).[140] L. Chaput, “Direct solution to the linearized phonon Boltzmann equation”, Phys.Rev. Lett. 110, 265506 (2013) (pg. 76).[141] W. Li and N. Mingo, “Thermal conductivity of fully filled skutterudites: role of thefiller”, Phys. Rev. B 89, 184304 (2014) (pg. 76).[142] W. Li and N. Mingo, “Ultralow lattice thermal conductivity of the fully filled skut-terudite YbFe4Sb12 due to the flat avoided-crossing filler modes”, Phys. Rev. B 91,144304 (2015) (pg. 76).[143] M. S. Green, “Markoff random processes and the statistical mechanics of time-dependent phenomena. II. Irreversible processes in fluids”, The Journal of ChemicalPhysics 22, 398–413 (1954) (pg. 83).[144] R. Kubo, M. Yokota, and S. Nakajima, “Statistical-mechanical theory of irreversibleprocesses. II. Response to thermal disturbance”, Journal of the Physical Society ofJapan 12, 1203–1211 (1957) (pg. 83).[145] B. Zhang, Z. Fan, C. Y. Zhao, and X. Gu, “GPU PBTE: an efficient solver for threeand four phonon scattering rates on graphics processing units”, Journal of Physics:Condensed Matter 33, 495901 (2021) (pg. 83).[146] Z. Han, X. Yang, W. Li, T. Feng, and X. Ruan, “FourPhonon: an extension moduleto ShengBTE for computing four-phonon scattering rates and thermal conductivity”,Computer Physics Communications 270, 108179 (2022) (pg. 83).[147] S. Lepri, R. Livi, and A. Politi, “Heat transport in low dimensions: introduction andphenomenology”, in Thermal transport in low dimensions: from statistical physicsto nanoscale heat transfer , edited by S. Lepri (Springer International Publishing,Cham, 2016), pp. 1–37 (pg. 84).[148] J. Lukkarinen, “Kinetic theory of phonons in weakly anharmonic particle chains”, inThermal transport in low dimensions: from statistical physics to nanoscale heat trans-fer , edited by S. Lepri (Springer International Publishing, Cham, 2016), pp. 159–214(pg. 84).[149] N. H. Protik, C. Li, M. Pruneda, D. Broido, and P. Ordejo´n, “The elphbolt abinitio solver for the coupled electron-phonon Boltzmann transport equations”, npjComputational Materials 8, 28 (2022) (pg. 84).96Appendix AFunctional form of the TersoffpotentialFor a monoatomic system of N atoms with coordinates {r1, r2, . . . , rN}, the bond-orderpotential developed by Tersoff [70] takes the following functional formU emp(r1, r2, . . . , rN ) =12∑i,ji 6=jfc(rij) [fr(rij) + bijfa(rij)] . (A.1)Here, rij is the distance between atoms i and j, fc is a cutoff term limiting the range ofatomic interactions, fr and fa are repulsive and attractive terms, respectively, while bij is amany-body term that factors in the relative orientation of neighboring bonds, reflecting thetendency of individual covalent bonds to become weaker with increasing coordination [70].In particular, these terms have the formfc(r) =1 if r < R−∆R,12− 12sin(pi2(r −R)∆R)if |r −R| ≤ ∆R,0 otherwise,(A.2)fr(r) = A exp(−λ1r), (A.3)fa(r) = −B exp(−λ2r), (A.4)bij = (1 + βmζmij )−1/2m, (A.5)ζij =∑kk 6=ik 6=jfc(rik)g(θijk) exp(λ33(rij − rik)3), (A.6)97θijk = arccos(rTij · rikrijrik), (A.7)g(θ) = 1 +c2d2− c2d2 + (h− cos(θ))2 . (A.8)To describe CNTs, we follow Ref. [72] and adopt the parametrization of Equations (A.2)–(A.6) and (A.8) as per Table A.1. Having fc(rij) = 0 for rij > 2.1 A˚, the short-rangedTersoff potential enables a significant reduction of computational cost. By referring to asingle-layer of graphene as sketched in Figure 1.1, it can be easily checked that, under thischoice of force cutoff, there are at most 10× 32 = 90 and 52× 33 = 1404 nonzero secondand third order force constants for any given atom, respectively.Table A.1: Parametrization of the Tersoff potential for graphitic carbon, taken fromRef. [72].A = 1393.6 eV λ3 = 0.0000 A˚−1 d = 4.3484B = 430.0 eV β = 1.5724× 10−7 h = −0.930λ1 = 3.4879 A˚−1 m = 0.72751 R = 1.95 A˚λ2 = 2.2119 A˚−1 c = 38049.0 ∆R = 0.15 A˚98Appendix BSpectral densities of acousticphonons from EMDThe complete set of spectral densities of EMD-derived acoustic normal mode coordinatesconsidered in the course of Section 4.1.4 is presented below. Specifically, for varying CNTdiameters, we show |Fω[Pk,ν(t)]|2 as in Figure 4.4 but now with eight columns, each of whichcorresponding to a wave number k = (2pi/λmax)w where w = 1, 2, . . . , 8. On the basis of thedata in Figures B.1–B.4, we perform phonon lifetime calculations for LA and TW modesaccording to Equation (3.2.4). The extracted lifetimes are plotted in Figure 4.5 of the maintext.99Figure B.1: Same as Figure 4.4 but for a (3,3) CNT.Figure B.2: Same as Figure 4.4 but for a (4,4) CNT.100Figure B.3: Same as Figure 4.4 but for a (7,7) CNT.Figure B.4: Same as Figure 4.4 but showing the full range of wave numbers 2pi/λmax ≤k ≤ 16pi/λmax.101Appendix CAsymptotic analysis ofthree-phonon scattering channelsFocusing on low-frequency acoustic modes in stress-free and stretched CNTs, here we ana-lytically derive the long-wavelength scaling relations of individual three-phonon scatteringchannels Γ±ν,ν2,ν3(k → 0+), as mentioned in Section 5.3.Considering Equations (2.2.15), (2.2.17) and (2.2.21), three-phonon scattering contribu-tions stemming from phonon triplets as per Equation (5.2.3) can be written as:Γ±ν,ν2,ν3(k) =~a8m3C∫ pi/a−pi/adk2[n0ν2(k2) +12 ± 12]n0ν3(k ± k2)n0ν(k)ων(k)ων2(k2)ων3(k ± k2)× |Vν,ν2,ν3(k,±k2,−k ∓ k2)|2 δ(Ω±ν,ν2,ν3(k, k2)),(C.1)where we use the fact that, as a result of energy conservation, phonon occupation numbersmust obey the following conditions:(1n0σ+ 1)(1n0σ2+ 1)=(1n0σ3+ 1), (C.2a)(1n0σ+ 1)=(1n0σ2+ 1)(1n0σ3+ 1)(C.2b)for absorption and decay processes, respectively. In general, given the conservation ofphonon quasi-momentum, the wave number of the third phonon is chosen such that k3 =k ± k2 − 2pi j/a falls within the zone [−pi/a;pi/a]. For the case at hand, |k| pi/a andwe can restrict our attention to scattering processes of normal type, meaning j = 0. Thecondition of energy conservation in Equation (C.1) is captured by the functionΩ±ν,ν2,ν3(k, k2) = ων(k)± ων2(k2)− ων3(k ± k2). (C.3)102For acoustic phonons with branch index ν and k → 0+, we can assume n0ν(k) ≈ kBT/~ων(k),and it follows that the small k behavior of Equation (C.1) is governed byΓ±ν,ν2,ν3(k) ∼[n0ν2(k∗) + 12 ± 12]n0ν3(k ± k∗)ων2(k∗) ων3(k ± k∗)|Vν,ν2,ν3(k,±k∗,−k ∓ k∗)|2∣∣∂Ω±ν,ν2,ν3(k, k∗)/∂k2∣∣ , (C.4)with the wave number k∗ = k∗ν,ν2,ν3(k) following from Ω±ν,ν2,ν3(k, k∗) = 0. Progress onΓ±ν,ν2,ν3(k → 0+) as per Equation (C.4) can now be made by taking into account the limitingbehavior of phonon dispersions ων/ν2/ν3(|k| → 0) and three-phonon coupling coefficientsVν,ν2,ν3(k → 0+). In the following, we analyze the three-phonon channels detected by ournumerical calculations in Section 5.3.C.1 Decay into two flexural modesWe consider transitions where a long-wavelength acoustic phonon from a linear branchdecays into two flexural acoustic modes, LA/TW→ FA+FA. According to Equation (C.3),k∗ is the positive root of0 = vνk −√T1k∗2 + Y()k∗4 −√T1 (k − k∗)2 + Y() (k − k∗)4, (C.1.1)where ν ∈ {LA,TW}, and the function Y() = v2LAR2/2−T2 is introduced to simplify thenotation. The condition of energy conservation can be solved analytically in two limitingcases by either assuming a strictly quadratic dispersion of flexural modes under stress-freeconditions, ωFA ∼ k2, or by adopting a strictly linear dispersion in the presence of tensilestrain, ωFA ∼ 1/2k. In the former scenario, we havek∗ =k2+√vν2Y(0)k −k24, (C.1.2)and ∣∣∣∣∂Ω−νFF(k, k∗)∂k2∣∣∣∣ = 2√2√Y(0)vνk − Y(0)k2, (C.1.3)where both expressions in Equations (C.1.2) and (C.1.3) go as k1/2. In the low-frequencyregime, n0FA(k) ≈ kBT/~ωFA(k), and by combining Equations (C.4), (C.1.2) and (C.1.3) wefindΓ−νFF(k) ∼ k−9/2 |VνFF|2 ∼{k−1/2 if ν = LA,k1/2 if ν = TW,(C.1.4)103where the standard long-wavelength approximation [50, 128] VLFF ∼ kk2k3 applies to LA→FA + FA but cancellation effects in the three-phonon coupling coefficient lead to VTFF ∼kk2k3(k2 − k3) [124] in the case of TW→ FA + FA.If, on the other hand, a strictly linear flexural mode dispersion is assumed, ωFA ∼ 1/2k,we havek∗ =(1 +vν√T1)k2and∣∣∣∣∂Ω−νFF(k, k∗)∂k2∣∣∣∣ = 2√T1, (C.1.5)which, when plugged into Equation (C.4), leads toΓ−νFF(k) ∼ −5/2k−4|VνFF|2(1− v2ν/T1)2∼{−5/2k2 if ν = LA,−7/2k4 if ν = TW,(C.1.6)where the long-wavelength approximations of VLFF and VTFF remain unaffected by tensilestrain, as is indicated by our numerical calculations of three-phonon couplings in Section 5.3.C.2 Coalescence of two flexural modesTo determine the scaling of the absorption channel Γ+FFν3 for ν3 ∈ {LA,TW}, we again solvethe condition of energy conservation as per Equation (C.3), Ω+FFν3(k, k∗) = 0, consideringthe cases ωFA ∼ k2 and ωFA ∼ 1/2k separately. In the unstrained case,  = 0, we havek∗ =vν32√Y(0) −√v2ν34Y(0) +vν3√Y(0)k − k2, (C.2.1)and ∣∣∣∣∣∂Ω+FFν3(k, k∗)∂k2∣∣∣∣∣ =√v2ν3 + 4√Y(0)vν3 k − 4Y(0)k2. (C.2.2)Here, the expression in Equation (C.2.1) goes as k, whereas Equation (C.2.2) converges to-wards a constant. In evaluating Equation (C.4), we can assume n0FA(k) + 1 ≈ kBT/~ωFA(k)and findΓ+FFν3(k) ∼ k−6 |Vν3FF|2 ∼{k0 if ν3 = LA,k2 if ν3 = TW,(C.2.3)where it is possible to reapply the long-wavelength approximations of VLFF and VTFF sinceVν,ν2,ν3(k,±k2,−k∓k2) is symmetric under exchange of phonon indices (ν, k)↔ (ν2,±k2)↔(ν3,−k ∓ k2).104For strained tubes,  > 0, we findk∗ =√T1− vν3√T1+ vν3k and∣∣∣∣∣∂Ω+FFν3(k, k∗)∂k2∣∣∣∣∣ = √T1+ vν3 , (C.2.4)and scattering channel contributions according to Equation (C.4) scale asΓ+FFν3(k) ∼ −2k−4 |Vν3FF|2 ∼{−1k2 if ν3 = LA,−1k4 if ν3 = TW.(C.2.5)C.3 One longitudinal and two twisting modesIn a similar manner, the strain-independent scaling relationsΓ−LTT(k) ∼ k2 and Γ+TTL(k) ∼ k2 (C.3.1)can be inferred from Equation (C.4) by noting that VLTT ∼ kk2k3.C.4 Absorption by optical phononsIf a low-frequency acoustic mode, labeled by (ν, k), is absorbed by two high-frequency opticalphonons, the condition of energy conservation reads0 = Ω+νOO(k, k∗) = ων(k) + ωOPTa(k∗)− ωOPTb(k∗ + k). (C.4.1)Such transitions manifest at some finite wave number k∗ ∼ const. either in the form ofintra-, a = b, or interbranch transitions, a 6= b, where two optical branches cross. For|k| pi/a, we can expand ωOPTb around k∗ to show that∣∣∣∣∂Ω+νOO(k, k∗)∂k2∣∣∣∣ = ∣∣∣∣∂ωOPTa(k∗)∂k2 − ∂ωOPTb(k∗)∂k2− ∂2ωOPTb(k∗)∂k22k − . . .∣∣∣∣∼{k if a = b,const. if a 6= b,(C.4.2)and it follows from Equation (C.4) that transition rates scale asΓ+νOaOb(k) ∼{k−1 |VνOO|2 if a = b,|VνOO|2 if a 6= b.(C.4.3)105Three-phonon coupling coefficients for one long-wavelength acoustic and two optical phononsVνOO are of at least linear order in k [128]. Our data in Figure 5.2 for optical interbranchtransitions of FA modes suggests VFOO ∼ k2, whereas the scaling behavior of TW modes asshown in Figure 5.3 suggests VTOO ∼ k for both intra- and interbranch transitions.106
Thesis/Dissertation
eng
Physics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Rottler, Joerg
Nojeh, Alireza
Theory and molecular simulation of lattice vibrational heat transport in carbon nanotubes
Text