Ultrahigh yet convergent thermal conductivity of carbon nanotubes from comprehensive atomistic modeling

We consider a $(10,0)$ semiconducting SWCNT, for which we model the interatomic interactions using the empirical bond-order Tersoff potential [36] using the parameters optimized to reproduce the phonon dispersion relations of graphene [37]. We first compute the thermal conductivity of an infinitely long SWCNT by both EMD and ALD simulations with periodic boundary conditions. In EMD simulations, $\kappa$ is computed in the Green-Kubo formalism as an infinite time integral of the heat flux auto-correlation function obtained from a sufficiently large ensemble of EMD runs in the microcanonical ensemble [38, 39, 40, 41]. We performed EMD simulations using the GPUMD package [42] for periodic $(10,0)$ SWCNT cells of lengths $52$–$207$ nm at 300 K. While the Green-Kubo integral would not converge for 1D systems with anomalous heat transport [2], our EMD simulations converge with respect to both correlation time and cell length and provide a value of $\kappa=2408\pm 83$  Wm^-1K^-1 (see Supplementary Material, Table ST1, Figures S1, S2).

Whereas the convergence of $\kappa$ in EMD simulations is supported by plenty of evidence, ALD-BTE calculations have so far struggled to provide an estimate for the thermal conductivity of SWCNT in the infinite length limit. ALD-BTE allows one to compute $\kappa$ from the nonequilibrium phonon population $\delta n_{\mu}=n_{\mu}-\bar{n}_{\mu}$ obtained by solving the linearized BTE [43], which for a system extended in one dimension reads:

where $\Gamma_{\mu\mu^{\prime}}$ is the scattering tensor that accounts for all the phonon scattering processes [44]. Here we consider only anharmonic scattering approximated to the first order of anharmonicity, i.e. considering only three-phonon scattering processes. The resulting thermal conductivity is

where $N_{q}$ is the number of $q$ points used to integrate the first Brillouin zone, $c_{\mu}$ and $v_{\mu}$ are the heat capacity and group velocity of mode $\mu$ ²²2The index $\mu$ contracts both the polarization and q-point index.. Due to prominent hydrodynamic effects [46, 47], the phonon mean free path (MFP) $\lambda^{\textrm{pp}}_{\mu}$ cannot be approximated by the relaxation time approximation, and needs to be computed by inverting the scattering tensor as $\lambda^{\textrm{pp}}_{\mu}=\sum_{\mu^{\prime}}(\Gamma^{-1})_{\mu\mu^{\prime}}v_{\mu^{\prime}}$. Matrix inversion is often circumvented by a self-consistent approach [48, 49] that was applied in the previous ALD-BTE calculations for SWCNT [33, 34]. However, the self-consistent series converges only if the off-diagonal elements of the scattering tensor are smaller than the diagonal ones in modulus $|\Gamma_{\mu\mu^{\prime}}/\Gamma_{\mu\mu}|<1$ with $\mu\neq\mu^{\prime}$. [50] This is not the case for the low-frequency modes of the (10,0) SWCNT (Figure S5), hence we directly invert the scattering tensor [51], with well-converged $q$-point sampling (Figure S4) ³³3$\kappa$ is computed using the $\kappa$ALDo software and a system of 5 replicas of the 40 atoms CNT fundamental unit cell. The second and third order interatomic force constants are calculated by finite displacement using the USER-PHONON package in LAMMPS [74]. The Brillouin zone is sampled using a 200 $q$-points reciprocal space mesh, which has been tested to be converged. The second and third-order derivatives of the interatomic force constants are obtained by finite difference displacement, as explained in [51].. ALD-BTE presents the advantage that it can be carried out both with the correct QM phonon statistics and in the classical limit, which is taken by solving Eq. (1) with $\hbar\rightarrow 0$. Both classical and QM calculations converge but to very different values of $\kappa$. The classical limit at 300 K gives $\kappa_{cl}=3190$  Wm^-1K^-1 which compares well to the value computed by EMD, especially considering that in ALD-BTE anharmonic effects are truncated, leading to overestimated $\kappa$ even at low temperature [53, 54, 55] The QM calculation gives $\kappa_{QM}=9,960$  Wm^-1K^-1 in the same ballpark as the experiments on mm-long SWCNTs [17]. The difference between $\kappa_{QM}$ and $\kappa_{cl}$ is among the largest for any material at RT: it stems from both the high Debye temperature $\Theta_{D}$ and the low dimensionality of CNTs. Such striking difference indicates that classical simulations for carbon-based nanomaterials can be used to find qualitative trends, but cannot provide quantitative predictions. Quantum effects become less important at high temperature or for materials with lower $\Theta_{D}$, for which MD provides reliable results..

NEMD allows one to study thermal transport in an analogous way to the experiments. Two local thermostats apply a temperature difference $\Delta T$ at opposite sides of a system with finite length $L$ to generate a nonequilibrium heat flux $J$ [56, 57, 58]. Evaluating ${J}$ at a stationary state and probing the temperature difference $\Delta T$ between the two reservoirs, with temperature controlled through Langevin dynamics, provides the thermal conductance per unit area $G=J/\Delta T$ and the effective length-dependent thermal conductivity $\kappa(L)=GL$ [59, 60]. Here we performed NEMD simulations to compute $\kappa(L)$ for the $(10,0)$ SWCNT with $L$ from 10 nm to 10 $\mu$m using GPUMD [42] exploiting a simulation protocol carefully evaluated for carbon nanostructures [59] (see SM, Fig. S3, Tab. S1, including Refs. [61, 62]). The NEMD results displayed in Fig. 1a show that $\kappa(L)$ apparently diverges with length, at least up to 10 $\mu$m, which is at the edge of the length accessible to NEMD simulations. Also these results are in agreement with former works, but they do not contradict the former EMD predictions, as $\kappa(L)$ remains below the infinite limit $\kappa(\infty)$ computed by EMD.

To compute the length-dependent thermal conductivity $\kappa(L)$ of the SWCNT at QM level, we use the ALD-BTE method. The customary way of introducing finite length in ALD-BTE is to correct the anharmonic phonon MFPs with a boundary scattering term through Matthiessen’s rule [63]: $1/\lambda^{L}_{\mu}=1/\lambda^{\textrm{diff}}_{\mu}+1/{L}$, which correctly describes both the ballistic ($L\rightarrow 0$) and the diffusive ($L\rightarrow\infty$) limits. Figure 1a shows a large discrepancy between ALD-BTE results obtained using the Matthiessen rule and the reference NEMD results: for any $L>10$ nm the ALD-BTE calculations systematically overestimate $\kappa(L)$. The reason for the failure of the Matthiessen rule is that it does not account correctly for the boundary conditions for Eq. (1). Specifically, in Eq. (1) the first equality is justified under the assumption of local thermal equilibrium $\frac{\partial n_{\mu}}{\partial x}=\frac{\partial n_{\mu}}{\partial T}\nabla T=\frac{\partial\bar{n}_{\mu}}{\partial T}\nabla T$, which is not the case in systems with finite length, as suggested by the nonlinear temperature profiles of our NEMD simulations (Figure S3). Using the McKelvey-Shockley (MS) flux method, Maassen and Lundström proposed a solution of the BTE for finite systems with consistent boundary conditions [64]. This method gives up the assumption of local thermal equilibrium and treats separately forward and backward fluxes along the transport direction. The finite-length thermal conductivity resulting from this approach for a 1D systems is:

in which the boundary scattering term $L/(2|v_{\mu^{\prime}}|)$ is added to the diagonal elements of $\Gamma_{\mu\mu^{\prime}}$. $\kappa(L)$ computed using the MS method in the classical limit is in excellent agreement with the reference NEMD calculations (Figure 1a). The agreement between carefully implemented NEMD simulations and well-converged ALD-BTE calculations with the correct boundary conditions gives solid ground to the finite-size MS approach [60] and resolves a long-standing puzzle about nanoscale heat transport.

We can now exploit ALD-BTE to calculate $\kappa(L)$ with quantum statistics to obtain predictions directly comparable to experiments for any length (Figure 1b). Quantum effects influence both the modal heat capacity and the MFPs (see Equations 2 and  3). Due to equipartition in classical statistical mechanics $c_{\mu}$(classical) is always larger than $c_{\mu}$(QM). Conversely, QM MFPs are longer than the classical ones (Figure 2), because the number of phonon scattering processes is lower in QM as high frequency modes are less populated. For this reason, the RT $\kappa_{QM}$ is much larger than $\kappa_{cl}$. At finite lengths, however, we observe a crossover at $L\sim$1 $\mu$m, below which $\kappa(L)_{cl}>\kappa(L)_{QM}$ as the effect on modal heat capacity dominates over that on phonon MFPs, which are limited by the length of the system. Furthermore, Figure 1b shows that the (10,0) SWCNT must be longer than the maximum thermal phonon MFP, i.e. few millimeters, for $\kappa(L)$ to converge to the infinite size limit. These results reconcile equilibrium and nonequilibrium simulations, ALD-BTE calculations and experimental measurements, all painting the same coherent picture of thermal transport in SWCNTs: the infinite-length limit of $\kappa$ is finite, but such limit is reached for mm-long SWCNTs, which are difficult to probe experimentally.

Hereafter, we analyze the origin of the large difference between $\kappa_{cl}$ and $\kappa_{QM}$, and why heat transport in SWCNTs is not anomalous as opposed to 1D models. Figure 3a displays the comparison between the classical and QM frequency-resolved thermal conductivity $\kappa(\omega)$ for the infinite-length SWCNT. The shaded areas represent the differential contribution $\Delta\kappa/\Delta\omega$. The contribution to $\kappa_{QM}$ from frequencies lower than 1 THz amounts to about 50% of the total, and that below 10 THz to 90%. In the classical calculation the contribution from the low-frequency modes is substantially reduced, as they are involved in a larger number of anharmonic scattering processes. The difference stems from the 2-to-1 phonon scattering processes that involve one high-frequency and one low-frequency mode combining into another high-frequency phonon. Classical statistics overestimates these processes due to the excess population of high-frequency modes, which are depleted according to the correct QM statistics. Figure 3b and c show the dispersion relations near the $\Gamma$ point with the color of the dots indicating the 3-phonon scattering phase space [65], i.e. the number of scattering processes available to each phonon (b) and the mode-resolved contribution to the total thermal conductivity (c). The longitudinal and torsional acoustic (LA and RA) modes are involved in relatively few scattering processes and provide the largest contribution to $\kappa$. The two degenerate transverse acoustic modes (TA) also substantially contribute to $\kappa$, but they also provide a large number of scattering channels. Conversely, higher order flexural modes have a large scattering phase space without contributing to heat transport mostly because of their small group velocity. Figure 3a also shows the finite length effect on the relative modal contribution to $\kappa$. As the long MFP of low-frequency modes is truncated by the system boundaries, the main contribution to $\kappa$ shifts toward higher-frequencies modes with shorter intrinsic $\lambda_{pp}$.

Finally, we investigate the difference between the SWCNT, a 3D system extended in one dimension with converging $\kappa(L)$, and a system with the same structure in which the atomic coordinates perpendicular to the CNT axis are constrained, so that the dynamics is 1D. The phonon dispersion relations of this model, compared to those of the 3D SWCNT are displayed in Figure S7. ALD-BTE calculations suggest that the lifetime of the longitudinal acoustic mode diverges for $q\longrightarrow 0$ as $\tau\propto q^{-\eta}$, with $\eta=1.76$ (quantum) and 1.96 (classical), thus supporting the divergence of $\kappa(L)\propto L^{\alpha}$ with $\alpha=0.43$ (quantum) and 0.48 (classical). Trends and exponents are compatible with both classical and quantum lattice dynamics calculations of anomalous transport in the FPU model [66, 67], showing that the suppression of the dynamics in the plane orthogonal to the transport direction restores anomalous thermal transport with $\kappa\propto L^{\alpha}$ as in 1D models.

We showed that different atomistic approaches, MD and ALD-BTE, independently provide the evidence that the thermal conductivity of SWCNTs has a finite bulk limit, in contrast with anomalous transport observed in 1D models in 3D space [68]. This result does not conflict with experimental observations of increasing $\kappa(L)$ for mm-long nanotubes, as the computed convergence length at RT is beyond 1 mm. To capture the correct magnitude of $\kappa$ it is essential to treat phonons quantum mechanically, meaning that classical MD simulations can provide qualitative trends but not accurate predictions below $\Theta_{D}$. For the (10,0) SWCNT considered here, the infinite length limit of $\kappa$ at the QM level is about three larger than in the classical limit. Modal analysis shows that both the two fundamental and the optical higher order flexural modes provide the main source for anharmonic scattering channels that make the thermal conductivity of SWCNTs finite. Our work reconciles formerly divergent numerical approaches to compute thermal transport in CNTs, and highlight the fundamental differences between 1D models that showcase anomalous heat transport, and nanostructures with finite thermal conductivity in the bulk limit, thus providing a substantial contribution to solve a critical problem across condensed matter and statistical physics. We expect this fundamental result to hold not only for CNTs, but also for nanowires and other nanomaterials that exhibit strong dependence of $\kappa$ on length at the macroscopic scale [69]. Divergence may still occur for suspended polymer threads [70, 71] and aligned atomic chains [72].