Thermal conductivity of one-, two- and three-dimensional carbon

Carbon atoms can form structures in one, two, and three dimensions due to its unique chemical versatility. In terms of thermal conductivity, carbon polymorphs cover a wide range from very low values with amorphous carbon to very high values with diamond, carbon nanotubes and graphene. Schwarzites are a class of three-dimensional fully covalent sp$^2$-bonded carbon polymorphs, with the same local chemical environment as graphene and carbon nanotubes, but negative Gaussian curvature. We calculate the thermal conductivity of a (10,0) carbon nanotube, graphene, and two schwarzites with different curvature, by molecular dynamics simulations based on the Tersoff empirical potential. We find that schwarzites present a thermal conductivity two orders of magnitude smaller than nanotubes and graphene. The reason for such large difference is explained by anharmonic lattice dynamics calculations, which show that phonon group velocities and mean free paths are much smaller in schwarzites than in nanotubes and graphene. Their reduced thermal conductivity, in addition to tunable electronic properties, indicate that schwarzites could pave the way towards all-carbon thermoelectric technology with high conversion efficiency.


Introduction
Carbon is arguably the most versatile element of all. Besides its importance in organic matter and life itself, carbon also assembles in structures covering the whole range of dimensionalities accessible in the real world. Up to a few decades ago only three carbon allotropes were know: graphite, diamond and amorphous carbon. The discovery of fullerenes expanded the list to include buckyballs [1] and carbon nanotubes [2]. More recently, the isolation of individual graphite sheets provided the family of carbon allotropes with yet another component named graphene [3]. While buckyballs are cagelike carbon molecules, i.e. 0D structures, and carbon nanotubes (CNTs) extend in one dimension (1D), atomically thin graphene takes the form of 2D membranes. Both nanotubes and graphene can be grown to reach sizes beyond the micrometer scale. Three-dimensional (3D) fully covalently bonded sp 2 carbon forms were predicted more than two decades ago [4,5,6,7,8], and were eventually grown, although in a nonperiodic form, by supersonic cluster beam deposition [9,10]. The periodic structures proposed in the earlier theoretical works correspond to a tessellation with hexagons and larger polygons, usually heptagons and octagons, of the negatively curved minimal surface proposed by the mathematician K. H. A. Schwarz in 1890.
Physical properties of carbon-based materials, in particular nanotubes and graphene, have been under intense investigation for the last two decades. Light atomic weight and strong covalent bonding provide carbon materials with special mechanical and thermal properties, which make them appealing for technological applications [11]. For example, diamond has the highest thermal conductivity (κ) among bulk insulators. Graphite has a high in-plane κ and a very small cross-plane one. In the case of graphite this anisotropy is due to the presence of strong covalent bonds between atoms in the same plane, but only a weak Van der Waals interaction between individual planes. CNTs and graphene have come to rival diamond presenting even higher thermal conductivities. These nanostructures display very similar chemical bonding features, where carbon atoms are connected by strong sp 2 bonds. The same type of strong sp 2 bonds characterize carbon schwarzites, which in principle could also have high thermal conductivities. However, schwarzites are also porous materials with very low density, bearing a resemblance with nano-phononic crystals, which then would favour a low κ. Indeed, nano-porosity and low density have been shown to lower the heat transport capabilities of materials [12,13]. Therefore, the actual κ of either crystalline or random schwazites is difficult to predict on the basis of general considerations and remains unknown.
Because of their high thermal conductivites carbon nanostructures have been regarded as potential building blocks for thermal management and passive cooling devices [14]. In contrast, routes to reduce κ of carbon nanostructures and improve their thermoelectric figure of merit have also been explored, laying the basis for carbonbased thermoelectric technology. Theoretical studies suggest that a significant reduction of κ with respect to pristine CNTs and graphene can be achieved in anti-dotted graphene lattice [15], graphene nanoribbons with disordered edges [16] or a specific edge design [17], and chemically functionalized graphene [18]. Although promising, these approaches rely on the use of atomically thin nanostructures, which would limit the range of processing and application. It is then worth investigating whether carbon schwarzites may respond to the need for low-κ 3D carbon. , and of two f cc schwarzites (3D), namely (C 28 ) 2 (c) and (C 64 ) 2 (d), both replicated in space and seen from the (110) direction. The different Gaussian curvature between the two schwarzitic structures appears evident from radii of the connections between the tetrahedral junctions.
In this work we investigate thermal transport in 1D, 2D, and 3D sp 2 carbon materials, computing the thermal conductivity of a single wall carbon nanotube (SWCNT), suspended graphene sheets and two crystalline f cc schwarzites (figure 1) by extensive molecular dynamics (MD) simulations. Whereas the thermal conductivity of graphene and SWCNT is very high, on the order of 10 3 Wm −1 K −1 , we find that schwarzites have a thermal conductivity on the order of 10 Wm −1 K −1 , which is surprisingly low for a crystalline carbon polymorph. We use anharmonic lattice dynamics (LD) to unravel the reason behind such large differences in thermal conductivity among carbon nanostructures with different dimensionality.
The remaining is organized as follows: in the next section we review the equilibrium MD method to compute thermal conductivity, the LD calculations, and the empirical interatomic potential used to describe the covalent carbon-carbon bonds. In section 3 we report the results of our thermal conductivity calculations by MD, including detailed size convergence studies, which are critical in this type of calculations. An interpretation of the MD results is provided in section 4, in which the phonon properties of the systems are comparatively analysed via anharmonic LD. The main results are then summarized in section 5.

Green-Kubo Thermal Conductivity
Linear transport coefficients can be calculated from equilibrium fluctuations of the respective conjugate flux via Green-Kubo relations [19,20]. In the case of heat transport, κ is calculated from fluctuations of the heat flux via the heat flux autocorrelation function (HFACF). For a system in equilibrium (in the absence of a temperature gradient) the heat flux averages to zero over time, but the HFACF decays with time and its decay time is proportional to the lifetime of the heat carriers (phonons). In the present work we calculate phononic thermal conductivity via atomistic equilibrium MD simulations, where the heat flux is calculated from atomic quantities such as velocity, energy and virial atomic stress. The Green-Kubo expression for each individual component of the thermal conductivity tensor can be written as where J is the heat flux, k B is Boltzmann's constant, T is the system temperature and V its volume. The HFACF should decay to zero for t longer than the correlation time of the system, and the respective integral saturates at a constant value. In practice, at long times the HFACF becomes noisy and the integral diverges because of this statistical noise. Therefore, κ αβ is taken as the stationary value of (1) before it drifts due to the accumulated noise. Finite size effects can have a large influence on the fluctuations of physical quantities during MD simulations, therefore it is necessary to ensure convergence of κ with system size as represented by the limit to large volume in (1). In general, κ is a second-order tensor with six independent components. However, symmetry can be applied to reduce the number of independent components for specific systems. In the case of 1D materials such as a SWCNT, the system can be oriented such that there is only one non-vanishing component, along its axis, κ zz . In the case of graphene, because of its hexagonal 2D symmetry, there are only two independent but equivalent non-vanishing components κ xx = κ yy . For 3D materials with cubic symmetry, such as the fcc schwarzites here considered, there are three equivalent nonzero components κ xx = κ yy = κ zz .
Finally, applying (1) to low-dimensional systems such as nanotubes and graphene requires a suitable definition of volume. Here we assign a nominal thickness of h = 0.335 nm to a sheet of graphene, such that its volume is given by the product of this thickness with the area of the sheet, V graphene = L x L y h. Similarly, the volume of a SWCNT is approximated to that of a cylindrical shell with an inner radius r − = r − h/2 and outer radius r + = r + h/2, where the radius of a (m, 0) SWCNT is defined as r = a CC √ 3m/2π and a CC = 0.144 nm is the carbon-carbon distance. In this way we can write V SW CN T = π(r 2 + − r 2 − )L z .

Anharmonic lattice dynamics and the Boltzmann-Peierls equation
In the harmonic approximation of LD the dynamical matrix D is defined as: where V is the potential energy of the system at equilibrium, m indicates atomic masses, r the distances between pairs of atoms, and the indexes i and j group the atomic and the Cartesian indexes. In a system of N particles D is then a 3N × 3N matrix. Solving the eigenvalue problem: provides the frequencies ω λ (q) and the normalized displacement vectors e λ (q) of the normal modes of the system. To compute the phonon lifetimes, one has to consider phonon-phonon scattering processes (normal and Umklapp). Generally the main scattering contribution comes from 3-phonon scattering processes, in which either two phonons (ω λ ,ω µ ) annihilate into a third phonon (ω ν ), or one phonon (ω λ ) gives rise to two phonons (ω µ ,ω ν ). Energy and momentum conservation determine the following selection rules for the two processes: where Q is a reciprocal lattice vector. The lifetime τ λ (q) of phonon λ at q-point q is given by [21,22]: where V λµν are the coupling constants between phonons λ, µ and ν, and n are the equilibrium phonon populations, which can be chosen according to either quantum (Bose-Einstein, BE) or classical statistics. Although the correct statistics for phonons is BE, it is possible to use classical occupancies to estimate the entity of quantum effects and to compare to classical molecular dynamics simulations, which necessarily obey classical statistics [23]. V λµν are the third-order anharmonic coefficients of the interatomic potential V projected on the normal mode coordinates: where i, j, k express Cartesian indexes and r i,j,k are atomic coordinates as in (2). Phonon propagation is described by the Boltzmann transport equation (BTE): The righthand side of (8) accounts for the scattering processes that create or destroy phonons. These include the already mentioned third order anharmonic scattering processes, but also higher order terms, as well as impurity and boundary scattering. The non-equilibrium occupation function can be written as n λ = n 0 λ + δn, and, assuming small temperature gradients (∇T ), one can linearize (8) by treating δn perturbatively. In stationary conditions the first term of (8) vanishes, and the linearized BTE is written in the form This linearized BTE can be solved at different levels of accuracy and complexity. The simplest approach is the relaxation time approximation (RTA). By this approach one calculates the relaxation time of each phonon, assuming that all the other modes have the equilibrium occupation (n 0 λ ). The resulting expression for κ is the sum over the contribution of all independent phonon modes where C λ is the specific heat per unit volume of each vibrational state. By making use of anharmonic LD as outlined above, one can determine the anharmonic phonon lifetimes, usually limited to three phonon processes. The relaxation times from three-phonon scattering processes are given by (6).

Tersoff Interatomic Potential
In both our LD calculations and MD simulations the interaction between carbon atoms is described by the bond order empirical potential proposed by Tersoff [24,25]. The Tersoff potential has been designed to model covalently bonded systems and has been extensively applied to Si, C and Ge systems [26,27]. The original parameter set provided by Tersoff was developed to reproduce structural properties of bulk semiconductors such as silicon, amorphous carbon and diamond, and has been recently re-optimized to accurately reproduce the vibrational properties of sp 2 carbon nanostructures [28]. In particular, this set of parameters provides a very good description of acoustic phonons, which are the major heat carriers in a material.
In the Tersoff formulation, the potential energy between atoms i and j can be written in the form where f R (r ij ) is a repulsive term, f A (r ij ) is an attractive term and b(θ ijk ) is a three body term which controls the strength of the attractive term depending on the number of neighbours k of atoms i and j, and on the angle between these atoms, θ ijk . The term f C (r ij ) is a cutoff function which limits the range of atomic interactions in the system. Typically, the cutoff radius is such that only first neighbour interactions are considered. Due to its bond order character, the Tersoff potential can accurately describe both sp 2 and sp 3 carbon systems, as well as different carbon structures with a single set of parameters. Therefore, we apply the Tersoff potential in all simulations presented in this work, with the parameter set given in [28]. ‡

Simulation Details
All structures are generated with a crystalline atomic arrangement corresponding to the optimized configuration with periodic boundary conditions. Each system is thermalized for 100 ps with a Nosé-Hoover thermostat [29] set at 300 K while its volume is kept constant. The system is them coupled to a barostat set at zero pressure and allowed to expand or contract the volume in order to achieve a relaxed configuration. This step depends on the system size and can last from 1 ns to 10 ns. The relaxed structure is used as the starting configuration for every simulation performed afterwards, where the velocities are set to 300 K and the system is again coupled to a Nosé-Hoover thermostat for 1 ns in order to decorrelate the initial configurations. Finally, the thermostat is decoupled from the system and the heat flux calculations are performed under microcanonical conditions, being recorded every 5 MD steps. During MD simulations the equations of motion are integrated with a velocity Verlet integrator with a timestep of 1 fs.
In order to increase the accuracy of κ at least 10 independent configurations are averaged, and the uncertainty is taken as the standard deviation. In the case of graphene and schwarzite, the final value of κ is averaged over the equivalent diagonal components of the conductivity tensor. The total simulation time of each configuration determines the accuracy of the HFACF. For the SWCNT a total time of 120 ns was required for each configuration. In the case of graphene, each simulation was at least 60 ns long, and for the schwarzite structures at least 40 ns.
Size convergence of κ is a very important aspect of MD-based heat transport simulations, not only because of the size effects on fluctuations, but also because ‡ In the particular case of fcc-(C 28 ) 2 , the distance between second neighbours can be very close to the cutoff radius R = 0.21 nm of the chosen parameter set. This issue can cause large peaks in the atomic forces during an MD simulation, which then requires a very small integration timestep. In order to keep the same timestep for all simulations, we shortened the cutoff radius to R = 0.19 nm exclusively for the fcc-(C 28 ) 2 systems. This change in cutoff radius does not affect any relevant physical property of the material. the phonon mean-free-path is limited by the size of the simulation cell. Therefore, we simulate supercells of increasing size with periodic boundary conditions. This is particularly important for low-dimensional systems with long wavelength phonons where the thermal conductivity might diverge [30,31,32,33]. In effect, the convergence trend of κ with system size provides an indirect estimate regarding the role of certain phonons in the heat transport process [34,35].

Carbon Nanotubes
Carbon nanotubes are excellent heat conductors. Regarding their electronic properties, SWCNTs can behave as metals or semiconductors, depending on their chirality. In semiconducting nanotubes, such as the (10,0) SWCNT, heat is mostly carried by phonons. The high value of κ found in CNTs is due to a combination of low dimensionality and a crystalline structure which favours long wavelength phonons [36].
Whereas 1D momentum conserving systems present a power-law dependence of the thermal conductivity on the length [30,31,32,33], κ ∼ L α , with α ∼ 0.33 -0.4, CNTs are not purely 1D and thus may have a large but finite κ. Although the size convergence of κ for CNTs is still debated [37,38], anharmonic LD calculations and equilibrium MD simulations show that the conductivity of SWCNTs seems to converge with length in the limit of infinitely extended systems [39,34]. Divergent κ and breakdown of Fourier's law for CNTs may be observed in the presence of large thermal fluxes [40].
Using (1) we calculate κ for SWCNT with supercells of increasing length along the nanotube axis as shown in figure 2. For very short supercell length, there is poor sampling of low-frequency flexural modes and the thermal conductivity is underestimated. Another contributing factor to this underestimation is the excessive scattering of long wavelength phonons due to the finite size of the supercell. Both limitations can be mitigated by increasing the length of the supercell in order to account for all relevant acoustic phonon modes appropriately. The results converge to κ = 1700 ± 200 Wm −1 K −1 for L > 5 nm. It is worth noting that this numerical value is in agreement with a previous work which employed the original parameter set provided by Tersoff, κ = 1750 ± 230 Wm −1 K −1 [34,41].

Graphene
This one atom thick 2D material presents remarkable physical properties which hold much promise to future technological applications. Its very high thermal conductivity and long phonon mean free path rival those of nanotubes. Similarly to CNTs, graphene is a low dimensional material and as such its thermal conductivity should, in principle, diverge with the sheet size. In the particular case of 2D materials, the conductivity shows a logarithmic divergence, such that κ ∼ log L [30, 31, 32, 33].   Nonetheless, by performing equilibrium MD simulations and calculating the thermal conductivity by the Green-Kubo formula in (1), we observe a convergence of κ with the size of the graphene supercell, as shown in figure 3. The reported conductivity values are calculated as the average over two components of the conductivity tensor, which are equivalent due to hexagonal lattice symmetry, such that κ = (κ xx + κ yy )/2. Each graphene supercell is approximately square (L = L x ≈ L y ), with the armchair edge oriented along the cartesian x-axis. The smallest supercell has ≈ 200 atoms and the largest ≈ 3 · 10 5 . We obtain a converged value for the conductivity κ = 1015 ± 120 Wm −1 K −1 , for a supercell with ≈ 2 · 10 4 carbon atoms [35]. Similarly to the case of SWCNT, the size convergence of κ for graphene has also been verified by anharmonic LD calculations [42].
An interesting aspect of the convergence trend shown in figure 3 is that, differently from the SWCNT case, κ decreases as the supercell is increased. Once again, as the supercell is increased we can probe longer wavelength phonons, which in graphene are typically acoustic flexural modes. In direct contrast to SWCNT, the role of these flexural modes in graphene is to limit κ, avoiding the log divergence expected for 2D systems. We have also recently reported that in the absence of such flexural modes the conductivity of graphene diverges as expected [35].

fcc-(C 28 ) 2 and fcc-(C 64 ) 2 Schwarzites
Starting from a single graphene sheet, it is possible to generate "closed" (positively curved) sp 2 structures, such as buckyballs, by replacing hexagonal rings by pentagons. Analogously, replacing hexagons by heptagons and octagons it is possible to generate open (3-dimensional negatively curved) sp 2 carbon crystalline structures [4]. Periodic negatively curved crystalline schwarzites have been predicted to be more stable than the fullerenes with equal absolute curvature, however it has so far not been possible to synthesize them, because 7-and 8-membered rings tend to get stabilized by adjacency. The formation of abutting large rings promotes the growth of purely graphitic spongy carbon, which retains a highly ordered short-range structure but forms self-affine amorphous minimal surface [9,43,10].
Crystalline schwarzites can be of either octahedral P-type or tetrahedral D-type, made of units with six or four branches, respectively. P-type units form simple cubic periodic structures whereas D-type units arrange themselves in diamond structures. Dtype schwarzites containing only 6 and 7-membered rings are particularly interesting, because it is possible to establish a one to one correspondence to fullerenes by replacing the 5-membered rings with 7-membered rings. The smallest schwarzite in this class is made of units of 28 carbon atoms forming solely 7-membered rings.
To address heat transport in this class of materials we consider the smallest D-type schwarzite, (C 28 ) 2 , and a relatively large one, (C 64 ) 2 made of 6-and 8-membered rings (Fig 1c and d). Stability, elastic properties and electronic structure of these systems were formerly characterized by semi-empirical and first-principles simulations [7,44].
These previous works showed that the smallest schwarzites, like nanotubes, can be either metallic or semiconducting and that there is a close relation between topology and electronic properties. Even though (C 28 ) 2 and (C 64 ) 2 are metallic, (C 28 ) 2 displays very localized states, and both have a low density of electronic states at the Fermi level, which leads to poor electrical conductivity and electronic thermal conductivity. In addition the electronic structure of larger systems with smaller Gaussian curvature would resemble that of rippled or distorted graphene, with an electronic energy gap at the Dirac point that directly depends on the curvature, thus suggesting the possibility for gap engineering. [45]   As we did for SWCNT and graphene we have computed the phononic thermal conductivity of the D-type schwarzites (C 28 ) 2 and (C 64 ) 2 by equilibrium MD simulations of cubic supercells, in order to verify its size convergence. Once again, κ is calculated as the average over the three equivalent components of the conductivity tensor. Our simulations show that κ converges for supercells with linear dimension of about 8 nm containing about 3·10 4 atoms. Surprisingly the bulk thermal conductivity of both schwarzites is 21 ± 2 Wm −1 K −1 , about 50 and 100 times lower than that of graphene and (10,0) SWCNT, respectively.
The density of (C 28 ) 2 and (C 64 ) 2 schwarzites are 1.3 and 1.2 g/cm 3 , i.e. 0.38 and 0.34 the density of diamond and about one half the density of graphene with a nominal thickness of 0.335 nm. Since κ is linearly proportional to the atomic density, these density ratios alone do not justify such a large reduction of the conductivity in schwarzites relative to other crystalline carbon structures.

Discussion
The schwarzites considered here are crystalline and periodic, there is neither mass disorder nor any other source of phonon scattering, and the type of chemical bonding is analogous to the other sp 2 forms of carbon considered. Therefore, the only possible reason for such a dramatic difference in the heat transport coefficient lies in substantial differences in the phonon properties at either harmonic or anharmonic level, which arise from the imposition of a negative curvature to graphene in order to form a 3D structure. We have characterized the vibrational properties of (10,0) SWCNT, graphene and schwarzite computing phonon dispersion curves and group velocities by harmonic LD. We have also calculated the phonon lifetimes by anharmonic LD as detailed in section 2. Phonon dispersion curves are shown in figure 5. While graphene displays three acoustic and three optical branches, SWCNT and schwarzites have a large number of optical branches, stemming from the relatively large number of degrees of freedom in the unit cell. Nonetheless, the modes with frequency up to ∼ 45 THz in SWCNT are very dispersive high order flexural modes [36], whereas for schwarzites we observe that the optical modes have relatively flat dispersions, and limit the frequency range of the acoustic phonons between 0 and ∼7 THz. There is also a remarkable difference with respect to the acoustic modes of graphene, where the longitudinal acoustic branch extends beyond 40 THz. These differences among phonon band structures imply large Figure 6. Phonons group velocities for (10,0) SWCNT, graphene, and schwarzites fcc-(C 28 ) 2 and fcc-(C 64 ) 2 , as a function of frequency calculated from phonon dispersions within the harmonic approximation. Phonons in schwarzite present much lower group velocities which accounts, at least in part, to the reduction in κ as shown in (10). differences in the phonon group velocities of carbon materials of different dimensionality, as shown in figure 6. Group velocities are much lower in schwarzites, when compared to SWCNT and graphene. Following (10), phonon group velocities contribute as a quadratic term to κ, therefore the observed suppression of v g (ω, q) justifies, at least qualitatively, the reduction of κ in schwarzites.
To get a quantitative picture of the specific phonon contributions to the thermal conductivity of schwarzites, we have computed the phonon lifetimes of (C 28 ) 2 and (C 64 ) 2 by evaluating (6) for the conventional cubic cells using 6x6x6 and 4x4x4 shifted q-point meshes [46,47]. The low frequency phonon lifetimes of the schwarzites are larger than those of the flexural modes of graphene, but decrease faster (τ ∝ ω −1.8 ), so that there is a cross-over at relatively low frequency ∼ 5 THz, after which the phonon modes of graphene not only have much larger group velocities, but also larger relaxation times. It is worth noting that the lifetimes of both the longitudinal and the transverse acoustic modes of graphene are nearly frequency independent [42], so that they also contribute to amplify the difference in κ between graphene and schwarzites.
The total thermal conductivity from anharmonic LD is given by (10). The possibility to perform anharmonic LD calculations using either classical or quantum phonon distribution functions permit the evaluation of the errors occurring in MD simulations from neglecting nuclear quantum effects. While we get a very good agreement between MD and anharmonic LD with classical statistics, using the correct quantum distribution functions we get about twice as large κ for both schwarzites. For (C 28 ) 2 κ cl = 28 Wm −1 K −1 and κ QM = 57 Wm −1 K −1 . For (C 64 ) 2 κ cl = 26 Wm −1 K −1 and κ QM = 64 Wm −1 K −1 . The comparison between quantum and classical normalized accumulation functions (figure 7), i.e. the cumulative contribution to κ as a function of the mean free path, shows that quantum effects largely enhance the contribution of phonons with micrometer long mean free path. These are typically low-frequency phonons with wavelength of several tens of nanometers, therefore their contribution may be attenuated by further disorder at the same length-scale.

Conclusions
In conclusion, we have shown that the thermal conductivity of two D-type schwarzites converge with supercell size to a value κ = 21 ± 2 Wm −1 K −1 . In spite of its carbon sp 2 bonded crystalline structure, the schwarzites have a thermal conductivity about 50 and 80 times lower than SWCNT and graphene, calculated by classical molecular dynamics at equilibrium with the same empirical potential. This remarkable difference between the thermal conductivities of similarly sp 2 -bonded carbon materials cannot be explained solely by the difference in atomic density among these materials. The reduction of κ stems from the change in dimensionality from low-D to 3D, which drastically changes the nature of the phonon modes in the system. The absence of flexural modes and the large number of atoms in the unit cells of schwarzites limits the frequency range of the acoustic modes and greatly reduces their group velocities. Nevertheless, the accumulation functions show that phonons with long mean free path larger than 1 µm provide the major contribution to heat transport. A 3D all-carbon material such as schwarzite, with such low thermal conductivity and its corresponding electronic properties might provide the route for high efficiency carbon-based thermoelectrics. In fact the observed reduction of κ from graphene to schwarzites is comparable to that that can be achieved by chemical functionalisation of graphene [18], but schwarzites present the advantage of being stable three-dimensional structures. In addition the thermal conductivity of schwarzites may be further reduced by intercalation, or naturally occurring mesoscale disorder, which may contribute to scatter longer mean free path phonons.