Functionalization mediates heat transport in graphene nanoflakes

The high thermal conductivity of graphene and few-layer graphene undergoes severe degradations through contact with the substrate. Here we show experimentally that the thermal management of a micro heater is substantially improved by introducing alternative heat-escaping channels into a graphene-based film bonded to functionalized graphene oxide through amino-silane molecules. Using a resistance temperature probe for in situ monitoring we demonstrate that the hotspot temperature was lowered by ∼28 °C for a chip operating at 1,300 W cm−2. Thermal resistance probed by pulsed photothermal reflectance measurements demonstrated an improved thermal coupling due to functionalization on the graphene–graphene oxide interface. Three functionalization molecules manifest distinct interfacial thermal transport behaviour, corroborating our atomistic calculations in unveiling the role of molecular chain length and functional groups. Molecular dynamics simulations reveal that the functionalization constrains the cross-plane phonon scattering, which in turn enhances in-plane heat conduction of the bonded graphene film by recovering the long flexural phonon lifetime.

In Fig. 1, we show a typical example of the measurement on the five samples of the control sample set (without functionalization).
The uncertainty of the heat flux measurement is due to the slight fluctuations in the direct current power supply (Agilent E3633A), especially at a high input power, but it is maintained at a low level of 1%. Standard calibration was carried out using a resistance   Figure 3 show, the heat produced by the wiring in the test chip means the experimentally measured temperature at the edge of the chip is higher than in the simulation.
The peak experimental hotspot temperature corresponds to a thermal conductivity of 80 Wm -1 K -1 for the Si test chip in the model. This is lower than the accepted thermal conductivity value of Si, but can be readily explained by the mesh convergence. Since solving the model required almost all available computational power, it was not possible to refine the mesh as much as wanted and ideally more layers of elements in the vertical direction would have been used if available. However, since the objective of this validation was to obtain a FE test chip with a similar temperature profile to the experimentally measured chip in order to accurately model the effect of the GBF, the validation was sufficient.
The GBF and Cu heat sink were then added to the simulation. The heat sink and GBF both have identical geometries as the experimental work. A schematic of the structure is contained in Supplementary Figure 4. The graphene oxide interlayer has been treated as one large thermal resistance as it only has a layer thickness of 5 nm, which is too thin to be reliably simulated. This thermal resistance was determined by the PPR measurements.
The PPR thermal resistance measurements were used to determine the effective resistance of both the GO and FGO interlayers using the following formula, where R 1 is the resistance between the GBF and GO interlayer, and R 2 is the resistance between the GO interlayer and SiO 2 . R tot denotes the total resistance of the GO/FGO interlayer and the two thermal resistances on each side. Since both the GO and FGO interlayers have a thickness of about 5 nm, the d/K term becomes negligible. Therefore, the total thermal resistances of the GO and FGO layers and their interfaces are 1.13E-7 and 3.48E-8 m 2 KW −1 respectively. Finally the thermal contact resistance between the heat spreader and the test chip with the GBF has been estimated at 1E-5 m 2 KW −1 as it just sits on the rest of the sample and this gave a result close to experiment for the bare chip.
These values were used in Figure 1c of the manuscript, where the model closely follows experimental trends. The simulation also predicts that the functionalization process offer a small additional decrease in hotspot temperature when compared to a GO interlayer without functionalization, agreeing with experimental findings.
The in-plane thermal conductivity κ in and cross-plane thermal conductivity κ cross , densities ρ and specific heat C p of the materials used in the model are contained in Supplementary   Table I.
The chip emits a constant power Q as the heat source. A constant temperature of T 0 0 = 273 K was applied to the bottom boundary of the silicon wafer. Thermal insulation condition −n · (−κ∇T ) = 0 was applied to the other outter boundaries of the setup. The heat dissipation process was simulated in the stationary regime by solving the well-known heat equation, where ρ is the mass density, C p is the specific heat, u is the velocity vector and Q is the thermal power generated in the chip.

SUREMENTS
Pulsed Photothermal Reflectance (PPR) method is used to measure the thermal interface resistances. Supplementary Figure 5 shows the setup of the pulsed photothermal reflection measurement used in this study. To enhance the heat absorption, a gold layer is evaporated on the surface of the GO and FGO layers after dip coating. The sample is first excited by a Nd:YAG laser pulse. This causes a fast rise in the surface temperature and then followed by a relaxation. The change of surface temperature is monitored by a probe laser which reflects off from the samples' surface. Since the relaxation time is governed by the thermal properties of the underlying layers and interfacial thermal resistance between the layers, by obtaining the temperature profile one can extract the thermal properties of the layers and thermal interface resistance between the layers through using a heat conduction model.  Table I of the manuscript and the error for thermal interface resistance measurements are about 20%.

TRANSMISSION USING ATOMISTIC GREEN'S FUNCTION
In the scattering theory, the system contains three coupled subsystems: two semi-infinite leads connected through the scattering region. The heat flux flowing in along the system axis writes where ℏω k is the energy quantum of the phonon mode k, v g,kz is the phonon group velocity of the phonon mode k z , n L,R is the phonon number on the left and right reservoir following the We note that d 3 k = dk x dk y dk z and v g,kz dk z = ∂ω/∂k z dk z = dω. The Eq. (4) reduces to Hence we identify the spectral phonon transmission function Ξ(ω) = t ω g(ω) where g(ω) = dk x dk y is the projected phonon density of states in the non-periodic directions of the system.
We probe the spectral phonon transmission function Ξ(ω) by atomistic Greens Function (AGF) [1,2] and the thermal conductance can be obtained by following the Landauer formula: where ω and ω max are the energy and the Debye frequencies. T refers to the mean temperature of the system, k B and ℏ represent the Boltzmann and the reduced Plancks constants, respectively. The transmission Ξ(ω) is obtained from a nonequilibrium Greens function approach as Tr[Γ L G s Γ R G + s ]. The advanced and retarded Green functions G + s and G s can be deduced from where ∆ is an infinitesimal imaginary part that maintains the causality of the Greens function and Σ L = K ab g L K + ab , Σ R = K ab g R K + ab are the self-energies of the left and right leads, the "+" exponent indicating the Hermitian conjugation. The leads refer to the infinite stacking of graphene layers. Finally, g L and g R refer to the surface Greens functions of the left and the right leads, while K s and K ab are the force constant matrices derived from the potential, for the molecule and between neighboring graphene sheet, respectively. The expression of the transmission also includes Γ .

SUPPLEMENTARY NOTE 5. ELECTRON TRANSPORT AND CONTRIBUTION TO THE CROSS-PLANE THERMAL CONDUCTANCE
Ab initio calculations were carried out using the quantum chemistry DFT code SIESTA Heat transport in graphene oxide has been a hot topic recently [5][6][7]. Present studies on the thermal transport of graphene oxide [1][2][3] have all focused on the influence of oxidation on the in-plane thermal conductivity of the graphene oxide. However, none of them studied the cross-plane heat conduction and/or the impact of oxidation on the heat flow in the adjacent graphene sheets of a graphene oxide -graphene hybrid system. Hence, we indeed agree with our Reviewer to conduct in-depth calculations on this so-far undiscovered subject.
We clarify the impact of oxidation of the graphene substrate on the in-plane and crossplane thermal transport in the supported graphene film. As is shown in Supplementary   Figure 9(a), the functionalized graphene layer at the bottom is oxidized on both sides of the functionalization. The procedure of our molecular dynamics simulation is the same as in the Ref. [1]. In Supplementary Figure 9(b)., we plot thermal conductivity κ of the graphene film and its thermal resistance R with the functionalized graphene oxide support versus the oxidation coverage.
For a very weak oxidation, the thermal resistance R of the graphene film and its oxidized substrate and its thermal conductivity κ are very close to that for a non-oxidized graphene support. As the oxygen coverage increases on the graphene support, the thermal resistance start to decrease. Such a reduction is due to the van-der-Waals (vdW) interaction between the oxygen atoms and the above graphene film, in addition to the direct graphene-graphene vdW interaction. Consequently, the basal-plane thermal conductivity of the graphene film decreases along with the increase of the oxidation coverage. Such a decrease is due to better thermal coupling between the graphene oxide support and the graphene film. However, as can be observed from Supplementary Figure 9(b), the reduction in both film-support thermal resistance and the basal-plane thermal conductivity remains weak, due to the weak contribution of the vdW forces between the oxidation and the graphene film above.

NAMICS
Adaptive intermolecular reactive empirical bond order (AIREBO) potential [8] was used to simulate the graphene's C-C interactions. AIREBO potential describes the bonded C-C interactions (r < 2Å, where r is the distance between two C atoms) with the reactive empirical bond order (REBO) potential of Brenner and uses Lennard-Jones (LJ) potential for the non-bonded Van der Waals C-C interactions (2Å < r < cutoff radius). The longrange electromagnetic and the short-range repulsive-attractive interactions in the molecule is taken into account through the ReaxFF potential, which uses distance-dependent bondorder functions to represent the contributions of chemical bonding to the potential energy [9]. The initial inter-plane distance of the graphene sheets is set to be 3.35Å. After the insertion of the molecule, the distance increases to 7.5Å.
By using EMD, the trajectories in the phase space, i.e., the time-dependent positions and velocities of the carbon atoms, are computed. The temperature T can be readily calculated from the velocities of the atoms in the simulation domain by using the time dependent kinetic energies. The thermal resistance R between two adjacent graphene sheets with temperature difference ∆T could be calculated by the following equation [10]: where k B is the Boltzmann constant, N 1 and N 2 refer to the number of degrees of freedom of the graphene subsystems sheets in interaction.
The in-plane thermal conductivity κ is calculated by the Green-Kubo formula: where V is the system volume, T refer to the equilibrium temperature of the system and Therefore the error in the estimation of the thermal conductivity of graphene introduced by using classical phonon distribution should be within a tolerable range.

MOLECULE
The total thermal conductance G t of the junction comprises of the conductance G m through a single molecule and the conductance G g through direct interaction of graphene sheets. Their relation writes simply, G g is proportional to the graphene contact area A, i.e., G g (A) ∝ A. We calculate with molecular dynamics the total thermal conductances for graphene systems of different contact area, as shown in Supplementary Figure 11. From the figure, we fit the data with a linear function and the y-intercept is the conductance of the molecule, i.e. G m = G t (A = 0). Thus we get G m = 82 pW K −1 .

SUPPLEMENTARY NOTE 9. CRITICAL LAYER NUMBER OF GRAPHENE FOR THERMAL PROPERTY SWITCH
We identify here the factors that determine the critical layer number for the thermal property change. These factors are I) number density of the functionalization molecules, II) functional groups of the molecule, and III) substrate of the graphene film.
i) For a diluted molecule density, as shown in Supplementary Figure 12a, the ripples are formed not only on the graphene layer which is directly covalently bonded to the molecule, but also on the top layer which has no direct covalent bonds with the underlying layer. Such geometrical defects are phonon scatterers and in this case the two layers of graphene with functionalization has lower basal-plane thermal conductivity compared to that without. In this density range, the ripples are well separated from each other and increasing the molecule density will decrease linearly the phonon MFP.
The critical number l c for the trend switch will be higher than 2 since more graphene sheets are required to recover the flat conformation, i.e. l c > 2. The precise layer number l c will depend on the molecule density.
ii) For a high molecule density, as shown in Supplementary Figure 12b, the ripples no longer remain separated and start to merge with their neighbours, recovering the flat surface in an extended area. Such merging hence recovers the interrupted phonon MFP and alleviates the detrimental influence of the separate ripples on the basal-plane thermal conductivity of graphene. In this case the two layers of graphene with functionalization have higher basal-plane thermal conductivity compared to that without functionalization, i.e. l c = 2.
• II. Functional group of the molecule.
To examine the effect of other functional groups, we use another functionalization molecule named HLK-5 which has been experimentally applied to reduce the thermal resistance between carbon nanotube (CNT) fillers and the polymer matrix [11]. The difference between the HLK-5 molecule and the amino-sliane used in the present study is that its amino group is doubly bonded to an aromatic ring of graphene/CNT. Such a double bond has a stronger bond strength than the single bond of the amino-silane molecule, thus resulting in a more severe static graphene distortion. Therefore in this case the two layers of graphene with functionalization have lower basal-plane thermal conductivity compared to that without.
• III. Substrate of the graphene film.
To investigate the effect of a different substrate on the critical number of layer, we use here a crystalline silica substrate to support the graphene film. The silica substrate could be easily functionalized with amino-silane molecules. The molecular dynamiscs simulation parameters are the same as in [12]. As can be seen from Supplementary The normal mode coordinate q(k, t) writes as a spatial Fourier Transform of the atomic displacements u j (t) [13,14] q(k, t) = and its time derivativeq(k, t) where j is the atom index, m j , r 0 j and e j (k) refer to the mass, the equilibrium position and the corresponding eigenvector of the atom j.
The potential and kinetic energies of the normal mode are and The total energy of the normal mode is One can demonstrate that where Γ(k) is the linewidth, equal to 1/[2τ (k)]. Thus, the lifetime can be approximated as To calculate the phonon dispersion curve, we first compute the spectral energy density of the mode k which correponds to a double Fourier Transform of the atomic velocities in Eq. (13) to the reciprocal space of frequency and wave vector In a harmonic system, phonons do not scatter. Therefore their E SED (k, ω) and exhibit delta-function type peaks in vibrational spectrum. In real material systems, phonons scatter and thus have finite lifetime due to the anharmonic interactions, leading to the linewidth broadening and thus giving Lorentzian peak shapes, as shown in Supplementary Figure 15.

HANCEMENT DUE TO INCREASED RELAXATION TIME
With the presence of the molecule, the ZA phonon branch has a finite frequency value at the Brillouin zone center, i.e. ω(k = Γ) > 0, as compared to the zero frequency for graphene without the molecule. Such massive ZA branch reduces the phonon group velocity v g in the frequency region near the zone center, which is susceptible to reduce the in-plane thermal conductivity if the phonon relaxation time τ is not changed, since the thermal conductivity C v is the specific heat per volume. However, as it can be seen in Fig. (6) of the manuscript, that the phonon relaxation time is substantially increased by the presence of the molecule over the entire phonon spectrum.
In order to clarify the competing effect of the increase phonon relaxation time and the decreased group velocity, we calculate the spectral in-plane thermal conductivity κ(ω) of the ZA branch by following the single-mode-relaxation-time approximation based on the kinetic theory, where C v (ω), v g (ω) and τ (ω) refer to the specific heat per volume, phonon group velocity and relaxation time extracted from Fig. (6) of the manuscript. The spectral thermal conductivity κ(ω) of the ZA branch is shown in Supplementary Figure 16. At low frequencies (ω/2π ≤ 2.2 THz), the substantially increased relaxation time compensates the effect of the reduction in the phonon group velocity v g by the massive ZA branch due to the presence of the molecules. At larger frequencies, the in-plane thermal conductivity is enhanced mostly due to the increase relaxation time. This remains valid for the rest of the phonon branches.

DUCTIVITY OF GRAPHENE ON THE MOLECULE DEN-SITY
We clarify the thermal conductivity dependence on the number density ρ of the functionalization molecules. We identify three regions where the thermal conductivity demonstrates different behavior following the molecule number density: • very diluted regime, ρ ≤ 0.01 nm −2 : For a very low molecule density, the basalplane thermal conductivity shows very slight reduction compared to the case "wo/ molecule", as shown in Supplementary Figure 17. In fact, when the number density of the molecule is further decreased (e.g. ρ = 0.0025 nm −2 corresponding to one single molecule on a 20 nm × 20 nm basal-plane area of graphene), the thermal conductivity of the "w/ molecule" case equals to that of the case "wo/ molecule" with the precision of the error bar.
• Region I (diluted regime, 0.01 nm −2 ≤ ρ ≤ 0.055 nm −2 ): When the molecule density increases, the ripples formed on the graphene surface which are caused by the point defects on the lower graphene layer due to the covalent functionalization start to scatter phonons and reduce the phonon meam-free-path (MFP). Note that in this regime, the ripples are well separated from each other and increasing the molecule density will decrease linearly the phonon MFP.
• Region II (intermediate regime, 0.055 ≤ ρ ≤ 0.11 nm −2 ): When ρ further in-creases, the ripples no longer remain separated and start to merge with their neighbours. Such merging recovers the interrupted phonon MFP and alleviates the detrimental influence of the separate ripples on the basal-plane thermal conductivity of graphene. In addition, the increase of the thermal resistance R of the supported graphene and its functionalized substrate starts to gain importance since the molecules reduces their van-der-Waals (vdW) interaction by intercalation. The heat channels through the molecules are weak compared to the direct graphene-graphene interlayer vdW interaction and therefore the overall cross-plane thermal resistance is increased. Such increased R will further favor the in-plane heat conduction thus enhancing the basalplane thermal conductivity.
• Region III (dense regime, ρ ≥ 0.11 nm −2 ): When the molecule density continues to increases, the molecules introduce numerous channels to transport heat in the cross-plane direction, whereas the thermal conductance of direct graphene-graphene interaction remains unchanged because the interlayer distance is kept stable by the molecules.
Therefore the overall substrate-graphene thermal conductance is increased and the basal-plane thermal conductivity starts to fall.
In heat-spreader applications using chemical functionalization, the molecule density is generally large. In the work of T. Meier et al., a Self-Assembled Monolayer with a number density of 4.67 nm −2 (see manuscript Ref. [38]), which is in the region III of Fig. R22, were used to study the cross-plane heat transport from a silica substrate to a gold atomic-force microscope tip. Therefore, we present in Fig. 5. of the manuscript the basal-plane thermal conductivity for ρ ≥ 0.11 nm −2 .

SUPPLEMENTARY NOTE 13. RAMAN SPECTRA OF GBF
Raman spectroscopy was performed on the GBF before and after spin coating the functionalized graphene oxide (FGO) layer between 10-20 nm thick, as shown in where Supplementary Figure 18. The D band, G Band and 2D band of the Raman spectra was at 1334 cm −1 , 1577 cm −1 and 2657 cm −1 , respectively, before spin coating. The noise level is represented by D band, with the relatively intense D band signifying the presence of defects due to chemical reduction. In contrast, after the spin coating, the intensity ratio of I D /I G was 1.4, compared to that of 0.9 before spin coating, with this attributed to the presence of FGO on the GBF surface. The multilayered structure of the GBF is evidenced by I 2D /I G < 1.

SUPPLEMENTARY NOTE 14. FTIR SPECTRA OF GBF
Fourier transform infrared spectroscopy (Spectrum Two, PerkinElmer) was conducted to study the functional groups on the GBF before and after spin coating with FGO. As shown in Supplementary Figure 19, before spin coating FGO, no significant peak is present between the range of 2600-3100 cm −1 and 600-1600 cm −1 . After spin coating, the peaks at 2884 cm −1 and 2929 cm −1 emerged, corrresponding to stretching vibrations of CH and CH −2 of the 3-Aminopropyl triethoxysilane (APTES) molecules. Moreover, peaks at 1034 cm −1 and 691 cm −1 indicate the presence of Si-O-Si and -Si-C-group, which provides more evidence for functionalization.