Resonant phonon modes in fullerene functionalized graphene lead to large tunability of thermal conductivity without impacting the mechanical properties

Graphene is a two-dimensional material of superlatives that is made from the hexagonal arrangement of $sp2$ bonded carbon atoms.¹ Owing to its exceptional electrical performance featured by the high carrier mobilities,² high optical transparency of $∼$97.7%,³ excellent in-plane thermal conductivity of $∼2000−5000Wm−1K−1$⁠,^4–8 and ultrahigh Young’s modulus of $∼$1 TPa,⁹ graphene has been linked with a cornucopia of potential applications in devices such as supercapacitors,¹⁰ organic solar cells,^11,12 efficient chemical sensors,^13,14 field-effect transistors,¹⁵ and high-frequency nanoresonators.¹⁶ Unfortunately, the unique properties of single graphene sheets are hard to be realized as graphene commonly forms graphite crystals where the individual sheets are linked by weak van der Waals interactions.¹⁷ However, the exfoliation of graphene sheets have been possible through surface functionalization with polymers,¹⁸ inorganic materials such as ZnO,¹⁹ aromatic molecules,²⁰ and covalently bonded fullerene molecules,^21–24 which inevitably change the physical, chemical, and mechanical properties of graphene.^25–30 Therefore, along with the understanding of the physical properties of graphene, the comprehensive knowledge of the effect of chemical functionalization on these properties is imperative to fully realize the potential and the impact that these novel materials can have on the current technology.

Among the various functional groups, fullerene-graphene hybrids have been envisioned to significantly enhance the performance of polymer solar cells²⁴ and create a new class of photovoltaic cells.³¹ Moreover, the all-carbon solar cells should have higher environmental stability due to their mechanical and chemical robustness.²³ Considering these realizable applications, it is instructive to systematically study the influence of fullerene functionalization on the thermal properties of graphene sheets.

In the context of understanding their thermal properties, various experimental and computational efforts have sought to understand the modification of heat flow across graphene nanoribbons achieved through strain engineering,^32,33 inclusion of grain boundaries,^34,35 vacancies and defects,^36,37 isotopic impurities,³⁸ foldings,³⁹ and edge roughness.^40,41 Kim et al.^26,27 and Chien et al.²⁵ used computational approaches to show that the thermal conductivity of graphene nanoribbons can be drastically reduced via the chemisorption of hydrocarbons on the surface of the nanoribbons. Pettes et al.⁴² experimentally investigated the role of polymeric residues on the thermal transport properties of suspended bilayer graphene and showed a large decrease in thermal conductance due to the scattering of heat carrying phonons in graphene by the residual polymer layer. However, to our knowledge, studies focusing on the effect of fullerene functionalization on the thermal transport and mechanical properties of graphene nanoribbons have been nonexistent in prior literature, and a systematic investigation into these topics will provide a benchmark for the comprehensive understanding of the physical properties of fullerene-graphene hybrid materials.

Therefore, in this work, we study the influence of covalently bonded fullerene molecules on the surface of graphene sheets on their thermal transport and mechanical properties through a series of molecular dynamics (MD) and lattice dynamics (LD) calculations. The thermal conductivities of the fullerene functionalized graphene sheets at different periodicities and surface coverages are predicted via equilibrium molecular dynamics (EMD) simulations for a range of temperature. Our atomistic simulations reveal that the thermal conductivity of pristine graphene can be lowered by more than an order of magnitude at room temperature via the introduction of covalently bonded fullerenes on the surface of the graphene sheets. This large tunability in the thermal conductivity is attributed to a combination of resonant phonon localization effects leading to band anticrossings and vibrational scattering at the $sp3$ bonded carbon atoms. The torsional force exerted by the fullerene molecules on the graphene sheets and the number of covalent bonds formed between the two carbon allotropes are shown to significantly affect the heat flow across the hybrid structures, whereas the size of the defects (or the number of carbon atoms in the fullerenes) has a negligible influence on their thermal properties. Furthermore, we show that even for a large surface coverage, the mechanical properties of these novel materials are uncompromised due to the covalent functionalization of graphene nanoribbons with fullerene molecules.

Fullerene molecules are covalently bonded onto the surface of the graphene sheets with side lengths $L$ through periodic inclusions, $P$⁠, to form the decorated graphene-fullerene sheets as shown in Fig. 1(a). Along with varying the period of fullerene inclusion, we also study the effect of 2-sided fullerene functionalization [as shown in the schematics in Figs. 1(b) and 1(c) for the side and top views, respectively] to assess the role of surface coverage on the thermal properties of our decorated graphene sheets. All MD simulations are performed with the LAMMPS code,⁴³ and the interactions between the carbon atoms are described by the adaptive intermolecular reactive empirical bond order (AIREBO) potential.⁴⁴ We apply periodic boundary conditions in the $x$- and $y$-directions, and we set the cross-sectional area for our computational domains to be $∼324nm2$⁠. The hexagonal faces of the fullerenes are attached to the hexagonal rings in the graphene sheets via [6+6] cycloaddition, which forms six C–C covalent bonds. First-principle calculations have shown that this procedure forms stable carbon-based structures.⁴⁵ Along with the decorated graphene sheets with six C–C covalent bonds, we also study structures formed by the [2+2] cycloaddition (with two C–C covalent bonds, which have also been shown to form stable structures) to access the role of covalent bonds on their thermal properties.

The computational domains are equilibrated under the Nose–Hoover thermostat and barostat,⁴⁶ which is the isothermal–isobaric ensemble with the number of particles, pressure, and temperature of the system held constant for a total of 2 ns at ambient pressure. Following isothermal-isobaric ensemble, we prescribe a canonical ensemble with the number of particles (N), volume (V), and temperature (T) held constant for another 1 ns to fully equilibrate the structures. After equilibration, the thermal conductivities of our graphene and decorated graphene sheets are predicted via the Green–Kubo (GK) approach under the EMD framework. Care must be taken to ensure that the thermal conductivities predicted from our simulations are not influenced by size effects; therefore, $L$ for our simulations is chosen to produce converged values of thermal conductivities for all structures.

In the GK formalism, the thermal conductivities of our pristine graphene and the decorated graphene sheets along the $x$- and $y$-directions are calculated as

Here, $t$ is the time, $T$ and $V$ are the temperature and volume of the systems, respectively, and $⟨Sx,y(t)Sx,y(0)⟩$ is the component of the heat current autocorrelation function (HCACF) in the $x$- or $y$-directions. Depending on the thermal conductivity of the computational domain, the total correlation time period for the integration of the HCACF ranges from 30 ps to 200 ps, which ensures that the HCACF decays to zero as shown in the inset of Fig. 2(a) for a 1-sided functionalized graphene with 2.5 nm period; the HCACF for pristine graphene is computed for a total of 200 ps as shown in Fig. 2(b); however, the inclusion of fullerenes leads to a faster decay in the HCACF and a shorter correlation time is required to calculate a converged thermal conductivity. The heat current is computed every 10 time steps during the data collection period, after which, integration is carried out to calculate the thermal conductivities. The converged thermal conductivity for the 1-sided functionalized graphene with 2.5 nm period length is determined by averaging from 5 ps to 30 ps as shown in Fig. 2(a). For all of our thermal conductivity calculations, we specify the thickness of the graphene layers as 3.4 Å, which is the van der Waals thickness. Along with the choice of the potential and the size of the computational domain, this definition of graphene layer thickness can drastically alter the quantitative thermal conductivities reported in the literature.^47,48 However, we note that our calculations for pristine graphene sheets are in good agreement with MD-predicted thermal conductivity of graphene described by a similar interatomic potential.⁴⁷ We also note that since the main goal of this work is to establish a comparative investigation into the effect of fullerene functionalization on the thermal transport properties graphene sheets, we do not try to consolidate and explain the differences between the prior literature on the quantitative thermal conductivity predictions of graphene using MD simulations.

The GK approach has been extensively used to predict the lattice thermal conductivity of various types of crystalline as well as amorphous material systems.^49–57 However, there has been considerable ambiguity in efficiently calculating the thermal conductivity via Eq. (1) due to uncertainties associated with finite simulation times and domain sizes.^{50,52,55,58–62} Therefore, to efficiently calculate the thermal conductivity via GK calculations, various approaches, such as exponential fitting of the HCACF,^52,58,59 and the integration of appropriately truncated HCACF have been implemented.^55,59,63 For the exponential fitting procedure, the time constants that are typically derived from the HCACF can be attributed to average phonon scattering rates. Specifically, the time constant related to the rapid decay (⁠$1/τs$⁠) corresponds to the scattering rate of short wavelength phonons, whereas the time constant related to the slower decay (⁠$1/τl$⁠) of the HCACF is related to the scattering rate of the long wavelength phonons and their respective contributions to thermal conductivity.^58,59,64 This approach of determining the average lifetimes for phonons will be used to support some of the conclusions made in this work. However, to calculate the thermal conductivities (as mentioned above), we employ the direct integration of the HCACF in accordance with the criterion set forth in Ref. 55 to accurately determine thermal conductivities by choosing the appropriate correlation time to calculate the HCACF as well as implementing sufficiently long simulation times to reduce uncertainties. Typically, correlation times have to be long enough so that the HCACF profile crosses zero multiple times during the integration period. We find that for our decorated graphene sheets, a correlation time of $∼$30 to 100 ps and a total integration time of 20 ns are sufficiently long enough to accurately determine the thermal conductivities as shown in Fig. 2(a); readers are referred to Refs. 55 and 59 for detailed investigations on how to reduce uncertainties in EMD calculations. For our calculations, the uncertainties are determined from three independent simulations, which results in a total of six sets of converged thermal conductivities (as calculated in the $x$- and $y$-directions for each simulation).

We also study the mechanical properties under uniaxial tension for our structures. The computational domains are deformed in the $x$-direction at a strain rate of $108s−1$⁠, while we impose zero pressure in the $y$-direction under the isothermal-isobaric ensemble. The strain is calculated every 0.1 ps as, (⁠$Lx−L$⁠)/L, where $Lx$ is the current length in the direction of the applied strain. The stress of the entire structure is also outputted every 0.1 ps to generate a stress–strain relationship for our structures.

Figure 3 shows the thermal conductivity of 1- and 2-sided functionalized graphenes as a function of period length calculated at room temperature. The decrease in period length monotonically decreases the thermal conductivity of the 1- and 2-sided functionalized graphenes. For reference, the predicted thermal conductivity of a pristine graphene by utilizing the Green–Kubo approach is $∼208Wm−1K−1$ at room temperature. Note, for our pristine graphene, a longer integration time period of 200 ps is required to ensure that the HCACF completely decays to zero. Prior literature on EMD calculations for materials with similar thermal conductivities have also demonstrated that longer correlation times are needed for accurate thermal conductivity predictions.^50,55,61 The decorated graphene sheets with short period lengths possess much lower thermal conductivities in comparison to the thermal conductivity of a pristine graphene sheet. For the shortest period (⁠$P=1.2nm$⁠) 2-sided functionalized graphene, we calculate a thermal conductivity of $∼15Wm−1K−1$ at room temperature, which is more than an order of magnitude lower than that of pristine graphene and marks a $∼$93% reduction, the reasons for which will be discussed in detail below.

Figure 4(a) shows the temperature dependent thermal conductivities in the 100 K–600 K temperature range for our 1- and 2-sided functionalized graphenes with different period thicknesses. For comparison, we also plot the temperature dependent thermal conductivity for our pristine graphene that shows a $∼T−1$ trend indicative of Umklapp scattering processes dominating thermal transport. In contrast, the thermal conductivity for the 1- and 2-sided functionalized graphene sheets with 3.5 nm period lengths show a drastically reduced temperature dependence (⁠$κ∼T−0.18$ and $T−0.14$⁠, respectively). Figure 4(a) also shows the temperature dependent thermal conductivity for a 2-sided functionalized graphene with 2 nm period, which shows a complete lack of temperature dependence for the range of temperature studied in this work (⁠$κ∼T−0.07$⁠). The increase in the lack of temperature dependencies for our decorated graphene sheets with larger surface coverages of fullerene molecules is indicative of the reduction in the influence of anharmonic Umklapp scattering processes, dictating heat transfer across these systems.

To investigate the effect of the number of covalent bonds on the thermal transport properties of our decorated graphene sheets, we plot the temperature dependent thermal conductivities of our 2-sided functionalized graphene with 3.5 nm period formed with 2 (solid triangle symbols) and 6 (solid circle symbols) C–C bonds in Fig. 4(b). For the temperature range studied in this work, the reduction in the number of covalent bonds leads to an increase in their thermal conductivity as well as their temperature dependency. The increase in the thermal conductivity, as well as its temperature dependence for the structures made via the [2+2] cycloaddition approach, arises mainly because of lower defect scattering due to a lower number of covalent bonds and also because of the fact that anharmonicity largely affects the higher-frequency phonons at higher temperatures. Thus, a larger number of covalent bonds between the fullerenes and the graphene sheet are more efficient in scattering higher-frequency phonons. Furthermore, the AIREBO potential takes into account the dependence on dihedral angles (contrary to the Tersoff potential that has been extensively used to describe thermal transport across carbon-based structures⁶⁵) and implements the torsional four-body potential for all dihedrals. The reduction in the number of C–C bonds between the fullerenes and the graphene sheets reduces the torsional force. It is important to note that the functionalization of fullerenes on the graphene sheets leads to the transformation of the in-plane $sp2$ C–C bonds into the out-of-plane $sp3$-like C–C bonds.^66,67 The torsional four-body term leads to a higher scattering of vibrations at these $sp3$ C–C bonds, which results in the further decrease in thermal conductivity of the structures formed by the [6+6] cycloaddition method in comparison to the thermal conductivity of the structures formed via the [2+2] cycloaddition method.

We also plot the temperature dependent thermal conductivity calculated without including the four-body torsional term in the AIREBO potential for our 2-sided functionalized graphene with 3.5 nm period (hollow triangle symbols) in Fig. 4(b). Similar to the effect of reducing the number of bonds, excluding the torsional term leads to a larger temperature dependence and an increase in the thermal conductivity, further solidifying our conclusion that covalent $sp3$-like C–C bonds severely scatter vibrations in our decorated graphene sheets via torsional forces, consequently leading to a massive reduction in thermal conductivity. This is supported by previous computational works that have shown that the rigid body rotational motion of the fullerene superatoms can severely impact the thermal transport properties of these types of carbon based materials.^67–70 

We also investigate the effect of increasing the number of carbon atoms in the fullerene molecules by repeating our simulations for different allotropes of fullerenes for the case of 1-sided functionalized fullerene with 2.5 nm period, the results of which are shown in Fig. 4(c). Within statistical uncertainties, the thermal conductivities are similar for all fullerene allotropes, suggesting that changing the number of carbon atoms in the fullerene molecules has a negligible influence on the thermal conductivity of the decorated fullerenes. We have confirmed this by running additional simulations on our decorated graphene sheets with 1.6 nm period and different fullerene allotropes to comprehensively conclude that the size of the fullerene molecules does not significantly impact their thermal conductivities. This implies that the $sp3$-like C–C bonds between the graphene and the fullerenes, and not the size of the defects, are responsible in dictating their thermal properties.

Next, to investigate the origin of the large thermal conductivity reductions of our decorated graphene sheets in more detail, we use harmonic LD calculations to obtain the dispersion relations in the high symmetry directions for decorated graphene sheets with varying period lengths. Figures 5(a), 5(b), and 5(c) show the dispersion relations calculated from supercells of 2.5 nm lengths for decorated graphene sheets with varying period lengths and varying number of C–C bonds between the fullerenes and the graphene. Note, since the number of phonon modes increases with increasing size of the computational domain, we have used similar lengths for the supercells in the LD calculations. The functionalization of the graphene sheets with fullerenes results in a series of resonant bands across the entire Brillouin zone. The hybridization of the acoustic bands with these resonant bands is characterized by the splitting of frequencies into up and down branches leading to mini bandgaps. Similar hybridization effects have been observed in our previous work for fullerene inclusion on the side walls of carbon nanotubes,⁶⁷ experimentally demonstrated for acoustic and nanophononic metamaterials, and also has been theorized for coupled oscillators by solving the equations of motion.^71–78 These works have shown that the nanoresonators, which in our case are the fullerene molecules, introduce standing waves that hybridize with the modes in the underlying membrane (the graphene sheet) that ultimately leads to avoided level crossings and flat bands that characterize localized modes with greatly reduced group velocities.⁶⁷ This localized resonant method has been used to manipulate waves with characteristic wavelengths much larger than the size of the resonators, thus offering a unique way to lower the thermal conductivity without affecting the intrinsic composition of the underlying material. The effect of resonant modes on the dispersion relations due to strongly and weakly adsorbed atoms on the surface has been previously studied theoretically in Ref. 79.

For our fullerene-graphene systems, the hybridization and localization of eigenmodes for our decorated graphene sheets is a direct consequence of the strong covalent nature of the bonds formed between the graphene and the fullerene molecules. The nature of the C–C bonds between the fullerenes and the graphene sheets is shown to impact the resonant modes as observed by comparing Figs. 5(a) and 5(b) for decorated graphene sheets with 2.5 nm periods but with six and two C–C bonds, respectively. Reducing the number of bonds leads to a reduction in the hybridization effect for the low frequency vibrations as highlighted in Fig. 5(a) with a comparatively flatter dispersion relation, which ultimately influences the thermal conductivity of the decorated graphene sheets [as shown in Fig. 4(b), where the sheets made by the [2+2] cycloaddition method showed an enhanced thermal conductivity as compared to those made with the [6+6] cycloaddition method]. We also show that the increase in the number of carbon bonds from the increase in the periodicity and surface coverage of fullerenes leads to stronger hybridization effects as observed by comparing Fig. 5(a) for the 2.5 nm period and Fig. 5(c) for the 1.2 nm period lengths [as highlighted in Fig. 5(c)]. It should also be noted that the increase in the number of atoms in the fullerene molecules and, therefore, the size of the fullerenes did not affect the hybridization of the phonon dispersion, further solidifying our conclusion that the C–C bonding environment and not the size of the “defects” (fullerene molecules) are responsible in the modification of the thermal properties of graphene nanoribbons.

To visualize the hybridization effect discussed above, we plot the eigenvectors for the hybridized modes at $∼$0.42 THz for the two supercells of 1.2 nm and 2.5 nm period lengths in Figs. 5(d) and 5(e), respectively. The arrows shown in the figure describe the direction and magnitude of the particular eigenmode calculated from our LD calculations, which shows a hybridized mode resulting from the coupling between the phonon modes in the graphene and the fullerene resonators. This can be identified near to the resonant frequency [star symbols in Figs. 5(a) and 5(c)], where the mismatch in the vibrational amplitudes in the graphene sheet and in the fullerenes results in the flat bands. The mismatch in the vibrational amplitudes in the fullerene resonators and the graphene sheet causes the atoms in the fullerene molecules to vibrate in the opposite direction to those of the atoms in the graphene sheet as depicted by the arrows in Figs. 5(d) and 5(e). As mentioned above, this localization effect for the modes is characterized by the flat bands in the dispersion relation at both subwavelength and superwavelength frequencies that have highly reduced group velocities, which ultimately leads to the reduction in the thermal conductivity of the functionalized graphene sheets as shown in Fig. 2(b).

To further support our conclusion that the resonant modes effectively scatter the long wavelength as well as the shorter wavelength phonons in these materials, we compare the time constants derived from the exponential fitting of the HCACF as shown in Fig. 2(b) (and discussed in Sec. II). For our pristine graphene, we determine a $τs∼1.4ps$ and $τl∼43ps$⁠. The inclusion of fullerenes at 3.5 nm period length is shown to reduce the lifetime of long wavelength phonons drastically (⁠$τl∼9.1ps$⁠). Further decrease in the period length to 1.2 nm results in $τl∼0.4ps$⁠, which suggests that the increase in surface coverage of fullerene molecules leads to a drastic reduction in the average lifetime of low frequency modes that carry a significant amount of heat in graphene. Along with $τl$⁠, the faster decay associated with short wavelengths also decreases with fullerene inclusion, which further supports the claim that both subwavelength and superwavelength frequencies are scattered by the fullerene resonators.

Finally, to investigate the effect of fullerene functionalization on the mechanical properties of our graphene sheets, we plot the stress–strain relationship in Fig. 6(a), which is used to calculate Young’s modulus of our structures. The slope of the linear region as shown by the dashed-line in Fig. 6 gives the Young’s modulus, which for our pristine graphene is $∼$0.9 TPa consistent with the results of earlier experimental and computational studies.^9,80 The similarity of the stress–strain response for all of our structures (including the calculated Young’s modulus values), regardless of the surface coverage of fullerenes suggests that the functionalization does not compromise the mechanical strength of the graphene. Furthermore, we also investigate the influence of the size of the fullerenes on the mechanical properties for our decorated graphene sheets in Fig. 6(b), which shows Young’s modulus calculated for our structures with 1.6 nm period lengths with different fullerene allotropes. As is clear, increasing the number of carbon atoms has a negligible influence on Young’s modulus; note, the number of C–C bonds between the fullerene and the graphene also does not influence their Young’s modulus. This suggests that the functionalization of fullerenes on graphene sheets has a negligible influence on their mechanical properties, whereas their thermal conductivity can be tuned over a large range with different surface coverages. This is in contrast to the results of absorption of hydrogen and methyl functional groups on graphene sheets that drastically deteriorate the Young’s modulus.^29,81

In conclusion, through systematic atomistic simulations, we have shown that covalently bonded fullerene molecules on the surface of graphene sheets can provide a unique pathway to tune the thermal transport properties by controlling the surface coverage. Our results show that the thermal conductivity of fullerene functionalized graphene with short period lengths can be lowered by as much as $∼$93% compared to the thermal conductivity of a pristine graphene, which is attributed to a combination of resonant phonon localization effects, leading to band anticrossings and vibrational scattering at the $sp3$ bonded carbon atoms. The main advantage of the approach, presented in this work, to tune the thermal conductivity of graphene lies in the unique opportunity where the electronic and the mechanical properties of the graphene are left unhindered.⁸² This could provide an avenue to engineer high thermoelectric efficiencies by massively reducing the vibrational thermal conductivity while enhancing the electronic conductivity in functionalized graphenes.

We appreciate support from the Office of Naval Research, Grant No. N00014-15-1-2769.