Model Approach to Thermal Conductivity in Hybrid Graphene–Polymer Nanocomposites

The thermal conductivity of epoxy nanocomposites filled with self-assembled hybrid nanoparticles composed of multilayered graphene nanoplatelets and anatase nanoparticles was described using an analytical model based on the effective medium approximation with a reasonable amount of input data. The proposed effective thickness approach allowed for the simplification of the thermal conductivity simulations in hybrid graphene@anatase TiO2 nanosheets by including the phenomenological thermal boundary resistance. The sensitivity of the modeled thermal conductivity to the geometrical and material parameters of filling particles and the host polymer matrix, filler’s mass concentration, self-assembling degree, and Kapitza thermal boundary resistances at emerging interfaces was numerically evaluated. A fair agreement of the calculated and measured room-temperature thermal conductivity was obtained.

The state-of-the-art investigation of polymer nanocomposites has advanced over the past decade and generated several comprehensive review articles [34][35][36][37][38][39][40][41].The key questions raised are how the filler content impacts the overall properties of the composite and to what extent the polymer-graphene interphase interaction variations, both with the morphology and surface activity of the filler, can be accounted for.
Among the other oxides employed in hybrid graphene/metal oxide and graphene oxide/metal oxide fillers, a unique place pertains to two well-known allotropies of titanium dioxide (TiO 2 ), namely rutile and anatase.Today, graphene@TiO 2 -and graphene oxide@TiO 2 -based polymer nanocomposites are promising materials that are not only used in energy storage applications but also in photocatalysis [44,45], gas sensing [46], solar cells [47], water splitting [48], and proton exchange membranes [49].Great attention has been paid to graphene@TiO 2 and graphene oxide@TiO 2 hybrid nanoparticles because titanium dioxide nanoparticles demonstrate nontoxicity and good compatibility with organic solvents and polymers, and they also display good hydrophilicity and mechanical and thermal stabilities [50][51][52].Also, the particular surface defect structure of titanium dioxide [53,54] endows it with highly reactive surface sites, which influence the interfacial charge and phonon transfer kinetics in composites with embedded TiO 2 [55,56].Following the comments of Wu et al. [42], the functions and synergistic effects of graphene/metal oxide hybrid nanoparticles embedded into a polymer matrix, which are related to the composite's transport properties, can be briefly summarized as follows: (1) Graphene nanoplatelets play the roles of 2D supports to uniformly anchor or disperse metal oxides with well-defined sizes, shapes, and crystallinity; (2) Graphene nanoplatelets act as 2D conductive templates to build a 3D conductive porous network to improve the poor electrical properties and charge transfer pathways of pure oxides; (3) Graphene nanoplatelets suppress the volume change and agglomeration of metal oxide nanoparticles; (4) Metal oxide nanoparticles suppress the re-stacking of individual graphene nanoplatelets; (5) Metal oxide nanoparticles can promote both charge and phonon out-of-plane transfer in multilayered graphene (MLG) nanoplatelets by binding dangling edge carbon bonds of neighboring monoatomic layers with the active sites located on the metal oxide nanoparticles.
Typically, the synthesis of MLG-TiO 2 nanocomposites includes in situ growth and covalent and non-covalent grafting [57].One should bear in mind that the rational design of hybrid nanoparticle architecture and the control of the morphology and phase composition of metal oxides on graphene can ensure reproducibility and a better understanding of the structure-property relationships in the nanocomposites.Also, a clear perception of the surface chemistry on graphene and metal oxides is crucial for the enhancement of the interfacial synergistic effects that are favorable for improved phonon transport.
The thermal conductivity of the composites is restricted by the high thermal interfacial resistance originating from strong interfacial phonon scattering.To obtain a thermal interface material with nanocomposites containing high thermal conductivity fillers, their concentration would thus be rather high.However, as the material is expected to be electrically insulating to prevent the malfunction of electronic devices at high working voltages, this requires low filler concentrations.Most frequently, reasonable thermal conductivities in the nanocomposites will inevitably result in enhanced electrical conductivities due to a low percolation threshold.Meanwhile, an increased thermal conductivity with filler loading and simultaneously decreased electrical conductivity were reported to arise from the attachment of functionalized molecules or nanoparticles to the filler surface [58,59].
Assuming that the temperature field obeys a diffusive equation, the kinetic theory predicts the Fourier law to be q = −κ∆T with the heat flux q, thermal conductivity κ, and temperature gradient ∆T.In solids, heat is carried by acoustic phonons and electrons so that the thermal conductivity can be represented as a two-term sum, κ = κ ph + κ e , with phonon (κ ph ) and electron (κ e ) contributions.In metals, κ e dominates due to a large concentration of free carriers ( κ ∼380-400 W/m K in copper).In insulating crystals, heat is transported via lattice vibrations, and the Boltzmann-Peierls approach [60] became one of the cornerstones in the theory of lattice thermal conductivity κ ph .
It is, however, generally agreed that heat conduction in carbon materials is usually dominated by phonons, even in graphite, which has metal-like properties [61].This is readily explained by the strong covalent sp 2 bonding and resulting efficient heat transfer via lattice vibrations, an assertion for which there is prior evidence.However, it is likely that the contribution of κ e will be significant for the doped materials [62].
In the relaxation time (τ) approximation, various scattering mechanisms, which limit the phonon mean free path (Λ = τv, with v being the phonon group velocity, which is approximated by the sound velocity when the acoustic dispersion can be ignored), are additive for the i-th scattering processes, so that τ −1 = ∑ i τ −1  i .The acoustic phonons are scattered by other phonons, lattice defects, impurities, free electrons, interfaces, etc.
The most elementary picture of heat conductivity is based on the analogy with the kinetic theory of gases when κ ph = 1  3 C p vΛ, where C p is the specific heat capacity.In a lattice, the oscillating particles are replaced with normal vibration modes, which possess well-defined wavelengths and frequencies.It is also necessary to take into account that the different modes obey different dispersion relations and thus have different group velocities v j depending on their wavenumber (or frequency).One can thus equate κ ph = ∑ j C j (ω)v 2 j (ω)τ j (ω)dω, where C j is the heat capacity, and ω is the phonon frequency, and the summation extends over all j phonon polarization branches (two transverse and one longitudinal acoustic branch).
Thermal conductivity is called intrinsic when it is limited by the crystal lattice anharmonicity.The intrinsic κ limit is reached when the crystal lattice is perfect, and phonons can only be scattered by other phonons.The anharmonic phonon interactions, which lead to finite κ in three dimensions, can be described by the Umklapp processes.The degree of the anharmonic interatomic forces is characterized by the Grüneisen parameter γ, which enters the expressions for the Umklapp scattering rates.Thermal conductivity is termed extrinsic when it is mostly limited by extrinsic effects, such as phonon scattering on rough boundaries or lattice defects [62].In epoxy-carbon composites, the electronic thermal conductivity is less than 2% of the thermal conductivity of the material and, therefore, may be considered negligible [63].
In nanostructures, κ is reduced by scattering from boundaries.Specular scattering with momentum conservation does not add to thermal resistance.Only diffuse phonon scattering from rough interfaces without momentum conservation ( p → 0 ) limits the free path Λ.
If the boundary scattering dominates, Λ, phonon dispersion is modified due to confinement, resulting in a changing υ value and complicating the size dependence of κ ph [62].C p is determined by the phonon densities of states, which are effectively modified by reducing the dimensionality.Consequently, κ ph ∼ T 3 is obtained at low temperatures in bulk materials and κ ph ∼ T 2 is seen in 2D systems.
It is unlikely, however, that the above equations can be easily generalized to the case of heterogeneous composites, such as polymer-graphene composites.Despite the long-standing interest in the development of a suitable model describing the heat transfer through the matrix-filler interphase, models of the thermal conductivity in such systems remains highly controversial [59].Our understanding of the heat transport phenomena is even more rudimentary for functionalized fillers such as graphene/metal oxide fillers.
Here, we prepared epoxy nanocomposites with MLG and embedded nanoparticles of TiO 2 and studied a temperature dependence of the thermal conductivity of the material in the range from about 40 to 300 K.The thermal conductivity of the composite was numerically computed as a function of the two-layer filler loading.The proposed softsphere approach and related concept of the effective thickness allowed us to simplify the thermal conductivity simulations by replacing a dotty interface with a continuous one and by describing interface thermal transport in the hybrid graphene@anatase TiO 2 nanosheets via the phenomenological thermal boundary resistance.

Temperature Dependence of the Thermal Conductivity
Figure 1a shows the measured temperature dependencies of the thermal conductivity of the neat epoxy (open circles) and its composites filled with hybrid MLG@TiO 2 nanoparticles with different filler contents, C f 1 :C f 2 = 0.01:0.01(closed circles) and 0.01:0.05(squares).The experimental data show aggregate behaviors that are typical of the thermal conductivity in disordered materials.Indeed, we previously reviewed extensive studies on amorphous oxides [64].
by describing interface thermal transport in the hybrid graphene@anatase TiO2 nanosheets via the phenomenological thermal boundary resistance.

Temperature Dependence of the Thermal Conductivity
Figure 1a shows the measured temperature dependencies of the thermal conductivity of the neat epoxy (open circles) and its composites filled with hybrid MLG@TiO2 nanoparticles with different filler contents,  :  = 0.01:0.01(closed circles) and 0.01:0.05(squares).The experimental data show aggregate behaviors that are typical of the thermal conductivity in disordered materials.Indeed, we previously reviewed extensive studies on amorphous oxides [64].
As an example, the computed () behavior in amorphous SiO2 (a-SiO2) is depicted in Figure 2. The () dependence can be described by calculating the total mean free path Λ as a combination of various phonon scattering mechanisms using the Matthiessen rule, where the sum runs over all scattering mechanisms, e.g., phonon-phonon, phononboundary, phonon-point defect scattering, etc.
(a) (b)  As an example, the computed κ(T) behavior in amorphous SiO 2 (a-SiO 2 ) is depicted in Figure 2. The κ(T) dependence can be described by calculating the total mean free path Λ as a combination of various phonon scattering mechanisms using the Matthiessen rule, where the sum runs over all scattering mechanisms, e.g., phonon-phonon, phonon-boundary, phonon-point defect scattering, etc.In disordered materials, the glass-like κ(T) dependence is characterized by the initial κ∼ T 2 rise at very low temperatures, followed by a characteristic plateau around 10 K (marked by arrow S in Figure 2) and a further slow increase, which is weaker than κ∼ T 2 , with the increasing temperature [66].Such κ(T) dependence is exemplified in Figure 2.
The physical mechanisms that may be responsible for such behavior have been extensively debated.The T 2 dependence at the lowest T has been shown to be well described by the scattering on statistically distributed localized tunneling entities [67,68].The plateau and gradually increasing κ at higher temperatures are well described by incorporating the soft-potential model [69,70].In disordered materials, the glass-like () dependence is characterized by the initial ~ rise at very low temperatures, followed by a characteristic plateau around 10 K (marked by arrow S in Figure 2) and a further slow increase, which is weaker than ~ , with the increasing temperature [66].Such () dependence is exemplified in Figure 2.
The physical mechanisms that may be responsible for such behavior have been extensively debated.The  dependence at the lowest  has been shown to be well described by the scattering on statistically distributed localized tunneling entities [67,68].The plateau and gradually increasing  at higher temperatures are well described by incorporating the soft-potential model [69,70].
The gradual increase in the thermal conductivity shown by the circles in Figure 1 is hence in qualitative agreement with the above prediction.The () dependence given by the squares is observed to saturate above ≈ 150 K.A small plateau can also be seen between 150 and 200 K for the net epoxy (open circles).The inclusion of MLG@TiO2 nanoparticles increases the thermal conductivity (closed circles in Figure 1).Increasing the concentration of TiO2 shows a further enhancement in  (squares compared with closed circles in Figure 1).It is not immediately obvious how the above models can be combined for the interpretation of the experimental results.

Development of the Model Formalism
In constructing a suitable model to describe the experimental data, we begin by picking the thermal conductivities of epoxy, graphene, and anatase TiO2 along with the thermal boundary resistances, which are given in Table 1.The gradual increase in the thermal conductivity shown by the circles in Figure 1 is hence in qualitative agreement with the above prediction.The κ(T) dependence given by the squares is observed to saturate above ≈ 150 K.A small plateau can also be seen between 150 and 200 K for the net epoxy (open circles).The inclusion of MLG@TiO 2 nanoparticles increases the thermal conductivity (closed circles in Figure 1).Increasing the concentration of TiO 2 shows a further enhancement in κ (squares compared with closed circles in Figure 1).It is not immediately obvious how the above models can be combined for the interpretation of the experimental results.

Development of the Model Formalism
In constructing a suitable model to describe the experimental data, we begin by picking the thermal conductivities of epoxy, graphene, and anatase TiO 2 along with the thermal boundary resistances, which are given in Table 1.It should be noted that the theoretical calculations based on Maxwell Garnett's effective media approximation [75] showed that, at low loadings (ϕ f1 ≤ 2%), the thermal conductivity does not depend on the graphene-epoxy thermal boundary resistance (r 1 ) when it varies in the range from 0.35 4b in ref. [76]).
It has been widely observed that the introduction of the Kapitza thermal boundary resistance between the MLG nanoplatelets and a host polymer matrix into the model can provide readable accounts of the thermal conductivity experiments and theoretical approaches.In hindsight, the heat transfer enhancement in nanocomposites is largely determined by the thermal boundary resistance, but the conventional heat conduction models for the interfaces fail to have the correct value of the thermal boundary resistance.
We begin with the development of our model by introducing the Kapitza resistance [79] R K between the thermally conductive filler nanoparticle with the thermal conductivity κ f and the 3D polymer matrix (κ m ), which arises from the large difference in the phonon density of states of the low-dimensional and bulk media.The Kapitza resistance is in series with the particle thermal resistance f /κ f , where f is the filler particle size [77].
For simplicity, let us consider the equivalent thermal resistance as f κ f e f f = f κ f + R K and assume that the effective thermal conductivity of the particle is [77] where η f is the Kapitza length.This can be introduced to describe the length scale, which is defined by the equality of the interface thermal resistance to the heat flux and the thermal resistance originating from the bulk [80].The existence of the thermal boundary resistance leads to a temperature drop ∆θ across the interface and [80] where Q i is the heat flux across the interface, whereas the subscripts m and f indicate the matrix and filler, respectively.The most widely accepted case for removing R K from the thermal conductivity expression is extending the problem to a large-scale system, when the interface thermal resistance can usually be neglected.By increasing the interparticle separation distance and filler size f a smaller portion of f κ f e f f is related to the interface.
Here, we approximate hybrid MLG@TiO 2 nanoparticles using the two bounding surfaces denoted by the subscripts 1 (MLG) and 2 (TiO 2 ) so that where R K1 and R K2 stand for the MLG-TiO 2 and TiO 2 -epoxy interfaces, respectively.The effective thermal conductivity of the hybrid nanoparticle (κ h,e f f ) is then readily obtained as However, the above approach implies the existence of a distinct boundary between the constituent phases.Moreover, there are no contributions in Equations ( 2)-( 5) to account for the geometrical shape of the filling particles.Therefore, to account for the interphase layers surrounding the fillers of different shapes, the thermal conductivity of the interphase regions (κ int ) was determined using the Kapitza thermal boundary resistance R K as [12] where h int is the interphase layer thickness.It should be noted that, since the interphase layers are composed of segments of macromolecular chains, the values of κ int are supposed to be close to κ m .It is therefore possible to roughly estimate a spread of the h int values provided that the corresponding R K value is known (or vice versa).
Using Equation (6), the above κ f ,e f f can be modified (see Supplementary Note S1).A graphical data analysis of both κ f ,e f f and κ h,e f f at varying interphase layer parameters is performed in Section 2.5.The numerical estimates of η f and the normalized effective thermal conductivities κ f ,e f f /κ f are given in Table 2 for the anatase, the in-plane (κ g11 = κ g22 ) and cross-plane (κ g33 ) thermal conductivities of the MLGs, and the hybrid particles of the sandwich-like MLG@TiO 2 and of the graphene-wrapped anatase architectures.In the evaluation of κ int , the value of 10 nm was taken for h int (see comments in Section 2.5).

Discussing the Temperature Dependence of the Thermal Conductivity in Nanocomposites
We fitted the measured κ(T) data to a power-law dependence κ(T) ∝ T n .The applied analysis yields several temperature regions of such dependence with the marked exponent values of n, as given by the lines in Figure 1b.
For the epoxy resins, a linear temperature dependence κ(T) ∝ T was observed in the region from 10 to 300 K [81,82].In our neat epoxy, n varies from about 0.4 to 0.5 upon increasing T (open circles in Figure 1b).This discrepancy can arise from a much lower cross-linking degree in our epoxy compared to that reported in the earlier studies.
Forming the epoxy/MLG@TiO 2 composite does not affect the n values remarkably (closed circles and squares in Figure 1b).Meanwhile, increasing the TiO 2 content exhibits the saturation of κ, which shows up prominently above T ≈ 150 K (squares).
These tendencies may be compared with theoretical and experimental results from the literature.In particular, single-layer graphene nanoribbons follow the κ(T) ∝ T n behavior below 100 K with the n value varying from 1.5 to 2.0 for different phonon modes [83,84].Graphite behaves as κ(T) ∝ T 2.5 [85].The theoretical calculations for the anatase TiO 2 single crystals show that κ(T) ∝ T −1 is within the temperature range of 100-500 K [86] with the dominating Umklapp phonon scattering.For the nanosized anatase particles, about 100 nm in diameter, the measured temperature dependence of κ exhibits a sinusoidal-like growth below 120 K, followed by a plateau in the range of 120-180 K with κ ≈ 0.8 W m −1 K −1 , and then it shows a gradual decrease [87].
It can therefore be assumed from the above analysis that the behavior of κ in anatase TiO 2 may have an impact on the saturation observed in the curve with the squares in Figure 1.As the κ C (T) dependence of our epoxy/MLG@TiO 2 composite samples in the whole studied temperature range cannot be explained quantitatively by the temperature behavior of its constituents, the contributions of the filler-polymer interphase regions should then be taken into account.A substantial difference between the open and closed circle data in Figure 1 above ≈ 150 K can now be thought to be assigned to the increasing volume of the interphase region with the increasing mass concentrations of the MLG and TiO 2 fillers.
More specifically, to calculate the temperature dependences of κ f ,e f f (T) (see Section 2.2) and to extract the partial contributions of MLG and TiO 2 into Equation ( 2), the temperature dependence of the Kapitza resistance R K (T) for all interfaces in the composite is therefore needed to provide basic information for the understanding of the measured κ C (T).
There are two widely used models for evaluating R K , the acoustic mismatch model (AMM) [88,89] and the diffuse mismatch model (DMM) [85,88,89].Both models give R K (T) ∝ T 3 below 50 K if the phonon dispersion can be neglected [88].However, the tem- perature dependence of the acoustic wave velocities incorporated into the models cannot be neglected at an increasing T, even if the velocity dispersion is neglected.Moreover, in the case of graphene, the phonon dispersion of the TA modes, which provide a dominant contribution to R K , cannot be neglected either.Therefore, it may be expected that the R K (T) dependence may diverge from the T 3 trend in the temperature range of interest.It has been observed that R K (T) at the solid-solid interface decreases with an increasing T from 100 to 300 K, with a smaller decreasing rate above 200 K [90].A decreasing R K value increases the thermal conductivity of the interphase region.
We may thus assume that a similar temperature behavior of R K is observed at the graphene-anatase (R K ≡ r 2,1 ), graphene-epoxy (R K ≡ r 1 ), and anatase-epoxy (R K ≡ r 2 ) interfaces.Therefore, the increasing κ C with the increasing T, followed by their saturations at higher loadings, is due not only to the increasing κ 1,2 value but also to the decreasing (stabilizing) r 1,2 .It may be suggested that this methodology can be used to describe r 1 (T) at the graphene-epoxy interface and then to incorporate r 1 (T), r 2 (T), and r 2,1 (T) in our model.However, a thorough theoretical analysis of κ C (T) is beyond the scope of this study.

Theoretical Simulation and Numerical Calculations
A great number of theoretical models for evaluating the thermal conductivity of composite materials have been developed [12].To construct an appropriate model aimed to quantitatively analyze the temperature-loading dependences of κ in our nanocomposites, three assumptions were taken into account: (1) The polymer matrix contains three solid phases, i.e., MLGs (phase 1), anatase TiO 2 nanoparticles (phase 2), and self-assembled hybrid MLG@TiO 2 nanoparticles (phase 3); (2) Every particle embedded into the matrix is surrounded by an interphase layer, which has a molecular structure that is determined by the particle-matrix interaction and differs from the polymer structure beyond the layer; (3) All phases are dispersed randomly and uniformly within the matrix.
We start from the model developed by Su et al. [12], where the authors obtained an algebraic equation to calculate the loading dependence of the thermal conductivity for graphene-filled polymer nanocomposites, taking the thermal boundary resistance into account.To adapt the original model of Su et al. to describe a nanocomposite filled with two-layer ("hybrid") nanoparticles, a certain geometrical configuration of the hybrid particles should be assumed.For the case of graphene-based hybrid nanosheets, a large number of possible configurations have been proposed in the literature (for a brief survey, see, for example, [42]).Here, we use the so-called sandwich-like configuration (SLG) for the MLG@TiO 2 hybrid nanosheets proposed earlier in ref. [91].Different geometrical versions of the SLG, including parallelepiped-shaped or spheroid-shaped graphene nanoplatelets assembled with either hard or soft spherical anatase nanoparticles, are given in Figure 3 (see also Supplementary Note S1).
algebraic equation to calculate the loading dependence of the thermal conductivity for graphene-filled polymer nanocomposites, taking the thermal boundary resistance into account.To adapt the original model of Su et al. to describe a nanocomposite filled with two-layer ("hybrid") nanoparticles, a certain geometrical configuration of the hybrid particles should be assumed.For the case of graphene-based hybrid nanosheets, a large number of possible configurations have been proposed in the literature (for a brief survey, see, for example, [42]).Here, we use the so-called sandwich-like configuration (SLG) for the MLG@TiO2 hybrid nanosheets proposed earlier in ref. [91].Different geometrical versions of the SLG, including parallelepiped-shaped or spheroid-shaped graphene nanoplatelets assembled with either hard or soft spherical anatase nanoparticles, are given in Figure 3 (see also Supplementary Note S1).In our case of a three-phase composite material, the basic equation of the model (see Equation (4) in ref. [12]) takes the following form:  ( ,  ,  ,  ,  ) •  (  ,  ,  ,  ,  ,  ) = 0 (7) In our case of a three-phase composite material, the basic equation of the model (see Equation (4) in ref. [12]) takes the following form: where ϕ n represents the volume concentrations of the constituent phases (the subscripts n = 0, 1, 2, and 3 correspond to the polymer matrix, free MLG nanoplatelets, free anatase nanoparticles, and hybrid MLG@anatase nanoparticles, respectively) and C 1 and C 2 are the mass concentration of graphene nanoplatelets and anatase nanoparticles, respectively, p 1 and p 2 are the relative mass portions of the free (unassembled) MLG and anatase particles, respectively, (0 ≤ p 1,2 ≤ 1), ρ n denotes the mass densities, {L 1,2 } is the set of geometrical sizes, r 1,2 and r 21 are the Kapitza thermal boundary resistances, κ n denotes the thermal conductivities (κ 1 includes both in-plane (κ g11 ) and out-of-plane (κ g33 ) components), and κ c is the thermal conductivity of the composite.The parameters p 1 and p 2 have been introduced to describe the self-assembling effect itself.A case of p 1 = p 2 = 1 means that the effect is absent (i.e., the case of two filling phases "1" and "2" in a host matrix), whereas the opposite extreme case of p 1 = p 2 = 0 implies complete self-assembling when there is a single (hybrid) constituent phase "1@2" in a matrix.The intermediate situation of 0 < p 1,2 < 1 corresponds to a three-phase composite comprising fillers "1" and "2" and their hybrid "1@2".It should be noted, however, that there is no direct evidence for the self-assembling of graphene nanoplatelets and anatase nanoparticles.This will require further experimental evidence to secure the validity of utilizing the self-assembling properties partly because of the fact that the non-oxidized graphene can only exist in a suspension.
The explicit expressions for the functions Φ n and ϕ n are given in Supplementary Notes S1 and S2, respectively.The functions have no obvious physical meaning and are used only to reduce Equation (7) to the compact form.In general, the model incorporates a whole number of both geometrical and material parameters related to the matrix and each constituent phase.The parameters are taken as independent variables to calculate the thermal conductivity κ C of the composite.
By using the model in the "soft-spheres" approximation, first of all, we complete a fitting procedure to match the calculated values of κ C to the measured values κ C,E , supposing that an extreme case of complete self-assembling (i.e., p 1 = p 2 = 0) takes place.In doing so, we input κ m , κ 2 , κ g33 , and r 1 from Table 1, and then assign fixed values for the interface layers' thicknesses h 1 = h 2 = 24 nm, the nanoplatelet's aspect ratio α = 0.01 κ g11 , r 1 and r 2 , and then try to satisfy the target condition in the form κ C = κ C,E by varying the graphene-anatase thermal boundary resistance r 21 .The parameters used in the fitting procedure are summarized in Table 3.
Table 3. Numerical evaluation of graphene@anatase thermal boundary resistance r 21 .As expected from Table 3, the resulting values of r 21 do not depend on r 1 .Since the numerical values for r 2 are lacking in the literature, they were varied by up to one order of magnitude (the values in the brackets and without the brackets in Table 3).Decreasing r 2 increases r 21 from 11% to 12.6% for C 1 :C 2 = 0.01:0.01 or from 46% to 52% for C 1 :C 2 = 0.01:0.05,where the smaller percentages correspond to lower κ g11 values.One can also see that the certain increment in r 21 is achieved while κ g11 is doubled from 600 to 1200 W•m −1 •K −1 .The resulting magnitude of r 21 (from about 0.25 × 10 −7 to 0.62 × 10 −7 m 2 •K −1 •W −1 ) seems to be appropriate.We can, for example, compare this value with the calculated room-temperature thermal boundary resistance of the Cu/Si interface, which is 5 Experimentally, it is rather difficult to distinguish between the thermal resistance of a thin layer and the thermal resistance of the substrate.Zheng et al. [92] treated the graphite/Ni interfaces as one effective layer and fit the time domain thermoreflectance (TDTR) data to two free parameters.A total interface thermal resistance of ≈2.3 × 10 −8 m 2 •K −1 •W −1 was obtained.The effective cross-plane thermal conductivity of ≈3.3 W m −1 K −1 was reported, whereas the in-plane thermal conductivity varied between ≈650 and 1000 W m −1 K −1 .
Multilayered graphene nanoplatelets typically contain structural defects, including point and edge defects, surface wrinkles, mesovoids, etc.The effects of defects on the interface's thermal conductivity in graphene/epoxy composites were observed by employing molecular dynamics simulations [93].Among the various defect types, Stone-Wales defects are the most effective in improving the interface's thermal conductivity in the composite.This is due to the enhanced vibration intensity of the out-of-plane low-frequency phonons, which could enhance the phonon coupling between the epoxy and graphene, thus increasing the conductivity.In turn, increasing the number of graphene layers in the MLG enhances the phonon scattering, thus decreasing the thermal conductivity.When the layer number exceeds four, the MLG properties become close to that of graphite, and the interface thermal conductivity stabilizes at the value close to that of graphite.
The variation in the interfacial thermal resistance with respect to the functionalization, isotope, and acetylenic linkages in graphene is also systematically examined using molecular dynamics [94].The graphene-epoxy interfacial thermal resistance can be remarkably reduced by using the covalent and noncovalent functionalization.Among different functional groups, covalent butyl is the most effective substance in reducing the resistance.
Strictly speaking, the thermal boundary resistance for the in-plane graphene-epoxy interface (r 1X = r 1Y ) differs from the cross-plane resistance (r 1Z ) (see Supplementary Note S1).These should be completely different thermal contacts since the graphene phonon density of states is shifted to a higher energy due to the in-plane vibration compared to the crossplane case.The sound velocity is also significantly different in the in-plane and cross-plane configurations.Meanwhile, an average value r 1av of r 1 can only be measured because the partial contributions of r 1X and r 1Z cannot be separated from each other.To obtain a quantitative estimate of the error incurred on the approximation of r 1 ≈ r 1X = r 1Z , the value of r 1av for parallelepiped-shaped graphene nanoplatelets can be expressed as Inputting r 1X = r 1Y , L Z = 50 nm, L X = L Y = 5 × 10 3 nm (average linear dimensions of our nanoplatelets) yields r 1 ≈ 1.0144•r 1Z (see Supplementary Note S1).Thus, approximating r 1 ≈ r 1X = r 1Z would result in a negligible error in the calculation of the thermal conductivity κ C of the graphene-epoxy nanocomposites.One may expect that the error might be even less in the hybrid graphene@anatase TiO 2 composite due to the thermal screening of graphene nanoplatelets by titanium dioxide fillers.Therefore, exploiting the case r 1X = r 1Y may seem, at first glance, to not bring any new qualitative conclusions compared to those given above.
We note that the properties of composites appear to depend not only on the size but also on the aspect ratio of the fillers α.It turns out that the interphase layer thicknesses, h 1 and h 2 , are sensitive to the polymer molecular structure and surface reactivity of the fillers.
One can therefore consider the dependence of r 21 on the aspect ratio α of the MLGs at varying h 1 = h 2 , as shown in Figure 4.It is seen that the value of r 21 depends rather sensitively on the magnitude of α, especially at a high α > 200, when the nanoplatelets can be considered as nanoribbons.In contrast, r 21 depends only slightly on h 1 and h 2 in the range of 10.0 nm ≤ h 1 , h 2 ≤ 30.0 nm.sensitively on the magnitude of α, especially at a high  > 200, when the nanoplatelets can be considered as nanoribbons.In contrast,  depends only slightly on ℎ and ℎ in the range of 10.0 nm ≤ ℎ , ℎ ≤ 30.0 nm.
The results in Figure 4 may help to increase the accuracy in the experimental determination of  .To carry this out properly, one first needs to obtain the distribution function () of the MLGs.Averaging the basic Equation ( 7) over an interval of () will then give the value of  , which, however, is beyond the scope of this paper.It should be borne in mind that, to make numerical estimates of the impacts of the interphase layers on the structural characteristics of a three-phase polymer nanocomposite using the relations given in Supplementary Note S1, it is necessary, in addition to the empirical parameters  and  , to predict two independent quantitative parameters of the interphase layers, namely the volume fractions  , and mass densities Different colors correspond to various values of interphase layer thicknesses: h 1 = h 2 = 10 nm (black curves), 15 nm (blue curves), 20 nm (magenta curves), and 24 nm (red curves).Solid lines correspond to r 2 = 1.0 × 10 −9 m 2 •K −1 •W −1 and dashed lines correspond to The results in Figure 4 may help to increase the accuracy in the experimental determination of r 21 .To carry this out properly, one first needs to obtain the distribution function F(α) of the MLGs.Averaging the basic Equation ( 7) over an interval of F(α) will then give the value of r 21 , which, however, is beyond the scope of this paper.
It should be borne in mind that, to make numerical estimates of the impacts of the interphase layers on the structural characteristics of a three-phase polymer nanocomposite using the relations given in Supplementary Note S1, it is necessary, in addition to the empirical parameters p 1 and p 2 , to predict two independent quantitative parameters of the interphase layers, namely the volume fractions ϕ i1,2 and mass densities ρ i1,2 .Conducting an experimental evaluation of p 1,2 is a difficult task.They depend on the technique used to prepare the graphene-anatase solution and the adhesion properties of the interacting surfaces.The actual process of self-assembling is stochastic, and different nanoplatelets have different covering degrees and spotted structures.It is evident that the probability of formation of each subsequent anatase layer decreases.Therefore, we may expect an increasing p 1,2 value when the occupation of platelets increases with the enhancing C 2 and input either a linear or exponential approximation for the p 1,2 (C 2 ) dependence.
We therefore suppose that the assumption of the constant values of p 1,2 at low loadings will give us essential physical insight.There is, of course, an uncertainty in choosing p 1 and p 2 .However, we may avoid the ambiguity in the numerical estimates by inputting p 1 = p 2 when the geometrical characteristics of the hybrid nanoparticle become independent of both p 1 and p 2 .
For simplicity, we postulate that the layer's shape is similar to that of the nanoparticle itself, and the layer's spatial structure is homogeneous along the surface of the nanoparticle and varies only along the surface's normal direction.With such an assumption, ϕ i1,2 can be determined in terms of the particle dimensions and layer thicknesses h i1,2 (see Supplementary Note S1).The molecular dynamic simulations testify that the h i values are limited to about 2R G , with R G being the radius of gyration of the polymer matrix [95,96].This restriction on h i is a logical consequence of the assumption that the formation of the interphase layer occurs only as a result of the effects of the molecular adhesive forces on the polymer-nanoparticle interface.For instance, for the DGEBA epoxy resins cured with tri-ethylenetetramine (TETA) hardener (C 6 H 18 N 4 ), the R G value estimate is 10-14 nm [97].The radius R G depends on cyclic deformation characteristics such as the operation frequency and loading direction used in the molecular dynamics simulations.Therefore, without the loss of generality, we assume that h 1 = h 2 and input h 1 = h 2 = 24 nm, supposing an average value of 12 nm for R G .The numerical calculations show that an increasing h i value will result in an increase in the nanocomposite's thermal conductivity κ C .
As for the values of ρ i1,2 , they obviously cannot diverse considerably from the epoxy density ρ 0 .The earlier molecular dynamic simulations showed that a range of polymer density fluctuations depends on the nature of the molecular interactions between a filling particle and the polymer macromolecules and it is mainly localized in the close vicinity of the particle of about a few angstroms [98].On the other hand, using simple considerations based on the free volume concept of the polymer structure allowed us to evaluate the upper limit for the interface layer density for the DGEBA epoxies as 1.064ρ 0 [4].Therefore, we assume 0.95ρ 0 ≤ ρ i1,2 ≤ 1.05ρ 0 to evaluate ϕ i1,2 .This spread in the value of ρ appears to cause a negligible impact on the calculation results.
Table 4 gives the values of the filler's (ϕ 1,2,3 ) and layer's (ϕ i1,2 ) volume concentrations, which were calculated by using Equations (S38)-(S44) in Supplementary Note S2.The difference creases with the increasing p 1 = p 2 and reaches its maximal value of 1.67•10 −5 for the case of p 1 = p 2 = 1.0, n = 1, C 1 = 0.01, C 2 = 0.05.One can see that the interphase areas occupy even greater volume in the composite compared to that in the fillers, only if p 1 = p 2 = 1.0 when ϕ 3 = 0, although the above interphase layer thicknesses of h 1 = h 2 = 10 nm are not so large.Table 4. Volume concentrations of the filler (ϕ 1,2,3 ) and interphase layer (ϕ i1,2 ).Top and bottom values correspond to ρ i1 = ρ i2 = 0.95ρ 0 and ρ i1 = ρ i2 = 1.05ρ 0 , respectively.It can be concluded from these calculations that the chosen sizes of the filler's loading allow one to avoid the interphase layer overlap and percolation effects.This will be discussed further below.
The overlapping effect can be avoided at low values of C 1,2 .Indeed, the average distance L P between adjacent nanoparticles (with linear dimension R P whose volume fraction ϕ p is distributed homogeneously within the polymer matrix) can be approximated by [99] where the values of ϕ p (C 1 , C 2 ) can be evaluated by using the relations given in Supple- mentary Note S2.By taking R P = 25 nm for the anatase nanoparticles and inputting , we find a minimal value of L 2 ≈ 55 nm.An extreme case of L 2,min = 2h 2 = 20.0 nm occurs at C 2,max = 0.2407.For graphene nanoplatelets with R P = L 2 2 = 25 nm, we obtain a similar result, L 1,min = 2h 1 = 20.0 nm at C 1,max = 0.2556.
In contrast, when the "soft-spheres" model is employed, the anatase covering thickness h = h(C 1 , C 2 ) is a continuous function of both C 1 and C 2 , whereas α = α(h) = (L Z + 2h)/(L X + 2h) is a continuous function of h.Then, ϕ 3PT can be derived from where h is a solution of with v 1 = v 1,1 and v 1 = v 1,2 given by Equations (S23) and (S24) in Supplementary Note S1 for parallelepiped-and spheroid-shaped nanoplatelets, respectively.The numerical solution of Equation ( 11) for p 1 = p 2 = 0 (complete assembling) gives the set of values of the percolation thresholds presented in Table 5 for the concentrations that are reachable in our experimental samples.Interestingly, varying the mass density ρ i1,2 by 5% around the value of ρ m at C 1 = 0.01 and C 2 = 0.05 caused the percolation threshold to change to not more than 0.01% and 0.023% for the parallelepiped-and spheroid-shaped configurations, respectively.It is apparent, therefore, that the effects of the interphase layer overlaps, and percolation cannot be deduced from the experimental κ data.

Discussing the Loading Dependencies
Figure 5 shows the loading dependencies of the thermal conductivities of the MLG@TiO 2epoxy nanocomposites calculated for different cases of interphase interactions.Figure 5a corresponds to the case when the filler-epoxy macromolecule interactions are neglected so that interphase layers cannot be formed.If the fillers remain unassembled, curve 1 is reproduced.In the opposite case of the fillers that are fully assembled in hybrid particles, curves 2 and 3 are obtained.A summary of the data in Figure 5a shows the following trends: Figure 6 shows the dependencies of  in the composite with fully assembled graphene and anatase particles on the thermal boundary resistance at the grapheneanatase interface ( ) computed at the fixed values of the thermal boundary resistance at the anatase-epoxy interface ( ).As expected, a decreasing  value results in a rapid increase in  .Instead, it is useful to note that decreasing the  value (or increasing either  or  ) increases the ratio  ( )/ ( ), e.g., 1.43 for  = 1.0 × 10 −7 m 2 K W −1 (curves 3 and 9 in Figure 6) against 6.51 for  = 1.0 × 10 −9 m 2 K W −1 (curve 10). Figure 7 shows similar dependencies upon  at fixed values of  .The computed results show a similar trend, as observed in Figure 6.Thus, an increased  value is observed upon the decreasing  value, while an increased  ( )/ ( ) value is realized with the decreased  , increased  , or increased  , from 1.55 for  = 1.0 × 10 −7 m 2 K W −1 (curves 3 and 9 in Figure 7) to 7.76 for  = 1.0 × 10 −9 m 2 K W −1 (curve 10).
Remarkably, a sensitivity of  to varying  increases with the increasing  .A general expression for the effective thermal conductivity  , of a particle embedded into a polymer matrix has been obtained in [101] in the effective media (1) κ C increases with the increasing C 2 regardless of whether or not the assembling effect occurs; (2) The particle assembling decreases κ C , and the greater the assembling degree p 1 + p 2 , the smaller the observed κ C value; (3) κ C decreases with the increasing graphene-anatase thermal boundary resistance r 21 , i.e., with the weakening of the MLG-TiO 2 interphase interaction.
It is also interesting to note that κ C increases with the increasing κ g11 , κ g33 or κ 2 values, as expected.
Figure 5b represents the case when the fillers interact with the epoxy macromolecular chains, and thus, interphase layers are formed.Curve 4, similarly to curve 1, corresponds to unassembled fillers, while curves 5 and 6 correspond to fully assembled fillers.One can see that including an interphase layer results in a decreasing κ C value.The other trends deduced from Figure 5a are also seen in Figure 5b.
Figure 6 shows the dependencies of κ C in the composite with fully assembled graphene and anatase particles on the thermal boundary resistance at the graphene-anatase interface (r 21 ) computed at the fixed values of the thermal boundary resistance at the anatase-epoxy interface (r 2 ).As expected, a decreasing r 21 value results in a rapid increase in κ C .Instead, it is useful to note that decreasing the r 2 value (or increasing either C 2 or κ g11 ) increases the ratio κ C (r 21min )/κ C (r 21max ), e.g., 1.43 for r 2 = 1.0 × 10 −7 m 2 K W −1 (curves 3 and 9 in Figure 6) against 6.51 for r 2 = 1.0 × 10 −9 m 2 K W −1 (curve 10).Evidently, the value of the thermal boundary resistance at the filler-polymer interface ( and  in our case) can vary via the chemical modification of the surfaces of the filling particles.This surface functionalization involves the conjugation of suitable functional molecular groups with the active surface sites of the particles.Generally speaking, the functionalization of fillers with the requisite moieties depends on the chemistry of the fillers.Recently, efficient techniques for the functionalization of carbon nanomaterials were developed (see refs.[102,103] and references therein), the most important of which may be classified as either covalent or non-covalent functionalization.
Covalent functionalization alters the state of bond connectivity between the fillers Figure 7 shows similar dependencies upon r 2 at fixed values of r 21 .The computed results show a similar trend, as observed in Figure 6.Thus, an increased κ C value is observed upon the decreasing r 2 value, while an increased κ C (r 2min )/κ C (r 2max ) value is realized with the decreased r 21 , increased C 2 , or increased κ g11 , from 1.55 for r 21 = 1.0 × 10 −7 m 2 K W −1 (curves 3 and 9 in Figure 7) to 7.76 for r 21 = 1.0 × 10 −9 m 2 K W −1 (curve 10).Remarkably, a sensitivity of κ C to varying κ g11 increases with the increasing C 2 .
polymer's macromolecules and the dangling carbon bonds located on the skeletons of graphene nanoplatelets (mainly on its edges and in structural defects).For example, the covalent interactions of oxidized graphene's basal and lateral surfaces with epoxy are dominated by hydroxyl, epoxy, and carboxylic molecular groups [104].As a result of functionalization, the translational symmetry of the outer graphene nanosheets of MLGnanoplatelets is disrupted by changing the  carbon atoms to  carbon atoms, and both the electronic and phonon transport properties of the sheets are influenced.Moreover, an enhanced conjugation between fillers and a host polymer matrix can happen through noncovalent interactions by attaching some specific secondary ligands on the filler's surface.Various approaches of filler-to-ligand noncovalent binding include electrostatic interaction, π-π stacking, and H-bonding [103].For MLG fillers, noncovalent functionalization, in contrast to covalent functionalization, has no effect on the π-π stacking of the multilayered nanoplatelet [105].
The grafted individual moieties or ligand-moiety chains in a functionalized filler can act as heat transfer channels between the filler and a host polymer matrix, decreasing  A general expression for the effective thermal conductivity κ f ,e f f of a particle embedded into a polymer matrix has been obtained in [101] in the effective media approximation.We have used this expression to write κ f ,e f f for free anatase particles (κ 2,e f f , Equation (S8) in Supplementary Note S1), free graphene platelets (κ g11,e f f and κ g33,e f f , Equation (S11)), and hybrid graphene@anatase particles (κ h,e f f , Equation (S14)).
Figure 8 shows the dependencies of normalized effective thermal conductivities of filling particles embedded into epoxy on correspondent thermal boundary resistances.
Thus, as expected, the effective thermal conductivities of the filling particles increase with the decreasing r value because this decrease is accompanied by an increased conductivity κ int of the interphase layers given by Equation (6).
for the complete self-assembling (at  1 =  2 = 0) of graphene and anatase in an epoxy polymer matrix.In this case, the variation in the TiO2 diameter does not affect the calculation results.It is also important to note that the effect of incomplete coverage of the graphene surface with anatase (observed in Figure 4) on the calculated thermal conductivity is very small at the small numbers of filling particles  and  that were used in the experiments (last three paragraphs in Supplementary Note S1).

Materials and Methods
MLG nanoplatelets for our experiments were prepared from the flakes of thermally expanded graphite by using the electrochemical technique described by Xia et al. [113].The structural properties of the nanoplatelets were tested by using X-ray diffraction (XRD) spectroscopy as described in our previous report [14].A distinct difference between the XRD patterns of the initial graphite flakes and graphene sheets was observed [14].Moreover, the Raman spectra of single-layer graphene particles and similarly fabricated samples of multilayer graphene particles were compared by Gorelov et al. [10], illustrating that our multilayer graphene nanoplatelets are composed of loosely bound graphene sheets.One can therefore suggest that the fabricated particles are multilayered graphene nanoplatelets rather than exfoliated graphite platelets.
To prevent the resulting MLG material from oxidizing, it was kept as the suspension.The particles were about 5 × 5 µm 2 in in-plane dimensions and 790 m 2 /g in specific surface area.The value of the specific surface area was determined by measuring the amount of physically adsorbed nitrogen gas from adsorption-desorption isotherms according to the standard Brunauer, Emmett, and Teller (BET) method [114] by using an Autosorb Station 3 apparatus of Quantachrome Instruments (Boynton Beach, FL, USA).
The TiO2-anatase nanoparticles were deposited on the MLG by adding them into the initial ethanol-based suspension of the MLG before its ultrasonic treatment.The specific surface area of the anatase particles was ~1500 m 2 /g.
It should be mentioned that the curves presented in Figures 6 and S2 predict the bounds for κ C , where they can be tailored by varying the thermal boundary resistances, r 2 and/or r 21 (in the case of complete self-assembling), if both C 1 and C 2 remain unchanged.
Evidently, the value of the thermal boundary resistance at the filler-polymer interface (r 1 and r 2 in our case) can vary via the chemical modification of the surfaces of the filling particles.This surface functionalization involves the conjugation of suitable functional molecular groups with the active surface sites of the particles.Generally speaking, the functionalization of fillers with the requisite moieties depends on the chemistry of the fillers.Recently, efficient techniques for the functionalization of carbon nanomaterials were developed (see refs.[102,103] and references therein), the most important of which may be classified as either covalent or non-covalent functionalization.
Covalent functionalization alters the state of bond connectivity between the fillers and a host polymer matrix.The functional groups form the covalent linkage between the polymer's macromolecules and the dangling carbon bonds located on the skeletons of graphene nanoplatelets (mainly on its edges and in structural defects).For example, the covalent interactions of oxidized graphene's basal and lateral surfaces with epoxy are dominated by hydroxyl, epoxy, and carboxylic molecular groups [104].As a result of functionalization, the translational symmetry of the outer graphene nanosheets of MLGnanoplatelets is disrupted by changing the sp 2 carbon atoms to sp 3 carbon atoms, and both the electronic and phonon transport properties of the sheets are influenced.
Moreover, an enhanced conjugation between fillers and a host polymer matrix can happen through noncovalent interactions by attaching some specific secondary ligands on the filler's surface.Various approaches of filler-to-ligand noncovalent binding include electrostatic interaction, π-π stacking, and H-bonding [103].For MLG fillers, noncovalent functionalization, in contrast to covalent functionalization, has no effect on the π-π stacking of the multilayered nanoplatelet [105].
The grafted individual moieties or ligand-moiety chains in a functionalized filler can act as heat transfer channels between the filler and a host polymer matrix, decreasing phonon scattering at this interface and thus reducing the thermal boundary resistance r or increasing the heat exchange between the filler and the matrix [106,107].In other words, grafted species increase the overlapping phonon density of states of the functionalized filler and the polymer matrix.
On the other hand, the surface functionalization breaks the regular surface structure of the fillers and reduces the intrinsic mean free path length of the phonons, acting as additional scattering centers for the phonons.Molecular dynamic simulations showed that the in-plane thermal conductivity of infinitely long graphene nanoribbons (κ ∞ ) decreases sharply with the increasing the grafting density g d [107].Specifically, κ ∞ reduces to about 14% from that for ungrafted graphene at g d = 1.6%.Furthermore, the thermal conductivity of grafted graphene at g d = 12.5% is only 3.86% of that for ungrafted graphene.Qualitatively similar effects may be expected for graphene functionalized by anatase nanoparticles.
Therefore, the impact of the filler functionalization on the composite's thermal conductivity appears to be twofold due to the competition between these two phenomena, an enhanced thermal conductance at the interface and the reduced intrinsic thermal conductivity of the filler.A certain optimal value g d,opt may be expected for filling particles, which is required to maximize the ratio κ C g d,opt /κ C (g d = 0).
We can further extend the relevance of this model by providing the way in which these functionalization mechanisms contribute together to produce effective thermal conductivity.One possible way to account for the functionalized graphene is to consider the intrinsic thermal conductivities of graphene (κ g11 , κ g22 , κ g33 ) and anatase (κ 2 ) and the corresponding thermal boundary resistances (r 1X , r 1Y , r 1Z , r 2 , r 21 ) as variables that are dependent on the grafting densities g d1 and g d2 .Thus, the resulting set of functions including κ g11 (g d1 ), κ g22 (g d1 ), κ g33 (g d1 ), κ 2 (g d2 ), r 1X (g d1 ), r 1Y (g d1 ), r 1Z (g d1 ), r 2 (g d2 ), and r 21 (g d1 , g d2 ) should be prescribed.Very little information related to these functions is available.For example, κ g11 and r 1 were empirically related to g d1 for graphene nanoplatelets embedded in the polyamide-6,6 matrix [107], although the dependences look somewhat intricate.
It should also be noted that the synergetic effects of functionalization or assembling fillers can be more pronounced when the filler surfaces are subjected to preliminary purification.It is well known that different atomic and molecular impurities conjugated with filler surfaces have detrimental influences on the intrinsic and interfacial thermal properties of the filler [108][109][110][111][112]. In particular, the trace metallic and metalloid impurities [108], acidic residues [109], hydrocarbones [110], and even silicon [111] are ubiquitous on graphene surfaces.
In this study, no effort was undertaken to intentionally purify or functionalize the fillers, neither graphene nor anatase, because the main aim of this study is to reveal the thermal synergetic effects of assembling graphene and anatase.Therefore, we have no need to cater for laborious and time-consuming techniques for the purification surface of the fillers and the subsequent quantitative estimation of the type and amount of residual functional groups.
Finally, we compare the proposed model predictions and experimental results shown in Figure 1a.The open and closed triangles give the distributions of the calculated values of the thermal conductivity in our composites at different key parameters of the model.Varying the aspect ratio α from 0.005 to 0.02, the varying thermal boundary resistances r 2 and r 21 from 10 −5 to 10 −10 m 2 •W −1 •K −1 , and varying the thicknesses of the interphase layers h 1 = h 2 from 10 to 30 nm yield a fair agreement of the calculated and experimental results.The bottom mesh of points (open triangles) characterizes the hypothetical situation of extremely small contact thermal resistances between the matrix and the fillers (minimum phonon exchange between the phases of the composite), whereas the upper mesh (closed triangles) corresponds to a strong interphase interaction with a significant overlap of the phonon density spectra of the fillers and the matrix.These calculations were performed in the soft spheres approximation (complete coverage of the surfaces of MLGs with TiO 2 ) for the complete self-assembling (at p 1 = p 2 = 0) of graphene and anatase in an epoxy polymer matrix.In this case, the variation in the TiO 2 diameter does not affect the calculation results.It is also important to note that the effect of incomplete coverage of the graphene surface with anatase (observed in Figure 4) on the calculated thermal conductivity is very small at the small numbers of filling particles C 1 and C 2 that were used in the experiments (last three paragraphs in Supplementary Note S1).

Materials and Methods
MLG nanoplatelets for our experiments were prepared from the flakes of thermally expanded graphite by using the electrochemical technique described by Xia et al. [113].The structural properties of the nanoplatelets were tested by using X-ray diffraction (XRD) spectroscopy as described in our previous report [14].A distinct difference between the XRD patterns of the initial graphite flakes and graphene sheets was observed [14].Moreover, the Raman spectra of single-layer graphene particles and similarly fabricated samples of multilayer graphene particles were compared by Gorelov et al. [10], illustrating that our multilayer graphene nanoplatelets are composed of loosely bound graphene sheets.One can therefore suggest that the fabricated particles are multilayered graphene nanoplatelets rather than exfoliated graphite platelets.
To prevent the resulting MLG material from oxidizing, it was kept as the suspension.The particles were about 5 × 5 µm 2 in in-plane dimensions and 790 m 2 /g in specific surface area.The value of the specific surface area was determined by measuring the amount of physically adsorbed nitrogen gas from adsorption-desorption isotherms according to the standard Brunauer, Emmett, and Teller (BET) method [114] by using an Autosorb Station 3 apparatus of Quantachrome Instruments (Boynton Beach, FL, USA).
The TiO 2 -anatase nanoparticles were deposited on the MLG by adding them into the initial ethanol-based suspension of the MLG before its ultrasonic treatment.The specific surface area of the anatase particles was ~1500 m 2 /g.
The commercially available CHS-EPOXY 520 (SpolChemie a.s., Usti nad Labem, Czech Republic) DGEBA-epoxy resin based on Bisphenol A diglycidyl ether (DGEBA), of epoxy group content (E-Index) 5.21-5.50mol/kg, with an EEW (Epoxy Equivalent Weight) of 182-192 g/mol was used as the neat resin.Diethylenetriamine (abbreviated as Dien or DETA) was used as a curing agent.DETA is a nitrogen-containing organic compound with the formula HN(CH 2 CH 2 NH 2 ) 2 [115].The epoxide to the hardener mass ratio was kept constant at 7:1.Details of the curing process can be found elsewhere [116].
The MLG mass loading values (C 1 ) for the nanocomposites were 0.5%, 1%, 2%, and 5%, whereas the TiO 2 mass loading values (C 2 ) were taken to be 0.5%, 1%, and 5%.To prepare the nanocomposites of the prescribed content (C 1 , C 2 ), the following procedure was used: The epoxide olygomer to the hardener mass ratio (R OH ) was kept constant at 7:1, i.e., R OH = m 0 m h = 7.We hence started by taking the epoxide olygomer portion with the mass m 0 = 105 g poured out into a quartz glass.Then, we added the fillers in the following two steps: first, we added the 10% graphene-ethanol suspension portion with the graphene mass and second, we added the dry anatase powder with the mass m a given by where m e = m 0 + m h = m 0 + 1 7 m 0 = 120 g.The as-prepared liquid composites were ultrasonically mixed until homogeneous suspensions were obtained.Then, they were stored in a vacuum chamber to remove ethanol and air bubbles.Finally, a hardener portion of m h = 15 g was added to the treated compound, and the resultant mixture was mixed.After adding the hardener to the suspensions, their further polymerization occurred at room temperature for 72 h.Once again, the suspensions were stored in a vacuum for the first two hours at a reduced pressure of about 10 −2 mbar to remove all the emerging gaseous constituents and thus prevent the formation of pores within the samples.
The morphology features of the nanocomposite were investigated with scanning electron microscopy (SEM) using a JEOL JSM-6490 microscope and with transmission electron microscopy (TEM) using a JEOL JEM-2100 microscope.Figure 9 shows the SEM image of the TiO 2 particles decorating the graphene nanosheets (a), and TEM images of a hybrid graphene@anatase TiO 2 nanocomposite are shown in (b) and (c).It is seen in Figure 9a that the TiO 2 particles are well dispersed over the MLG surfaces.The highresolution TEM image shown in Figure 9c displays a lattice fringe of about 0.33 nm, which can be attributed to the anatase (101) plane [117].
Molecules 2023, 28, x FOR PEER REVIEW 22 of 32 where  is the distance between the sensors, and  and  are the temperatures measured by sensors A and B, respectively.The samples were cylindrically shaped with a diameter of 12.5 mm and a height of 6 mm.Two holes with diameters of 1 mm were drilled in the samples from the top and bottom surfaces to the center, in which the temperature sensors were placed (A and B in Figure 10).After that, the holes were filled with epoxy resin for better thermal contact between the sample and the temperature sensors.BAP64-02 NXP diodes were used as temperature sensors.The diodes had rather small dimensions of 1.2 × 0.8 × 0.6 mm 3 , which made it possible to place them inside the sample without difficulty.It was revealed by the microscopy that the TiO 2 particles have a diameter of ~30-70 nm, while the thickness of the MLG nanoplatelets ranges from about 30 to 150 nm.One needs to keep in mind that these dimensions are randomly distributed through irregular MLG@TiO 2 platelets.So, one needs to choose the values of the MLG and TiO 2 thicknesses (both are taken to be 50 nm above), and then obtain the numerical continuation solution branches in κ upon varying parameters.Some results describing varying κ C upon varying MLG thickness are given in Supplementary Note S3.
Here, we used the two-probe measurement method to measure thermal conductivity [118].Our automatic measuring system is schematically sketched in Figure 10.The sample was placed between two copper discs (heater and thermostat in Figure 10).An electric resistance served as a heater, and the power P released on it was determined with the electric current and voltage applied to the resistance.Additionally, the heat flux through the sample is Q = P/S, where S is the top and bottom surface area of the sample.A thermostat was used to keep the temperature of the top surface of the sample constant.By maintaining the stationary heat flux, one can write the following relation to evaluate the thermal conductivity [119] where L is the distance between the sensors, and T A and T B are the temperatures measured by sensors A and B, respectively.The samples were cylindrically shaped with a diameter of 12.5 mm and a height of 6 mm.Two holes with diameters of 1 mm were drilled in the samples from the top and bottom surfaces to the center, in which the temperature sensors were placed (A and B in Figure 10).After that, the holes were filled with epoxy resin for better thermal contact between the sample and the temperature sensors.BAP64-02 NXP diodes were used as temperature sensors.The diodes had rather small dimensions of 1.2 × 0.8 × 0.6 mm 3 , which made it possible to place them inside the sample without difficulty.
We believe the diode temperature sensor was small enough to give rise to correct measurements since the ratio of the cross section of the sensor to that of the sample is less than 1%.There is little work reported on the usage of diode sensors to obtain thermal conductivity measurements of epoxy systems [120].It may therefore be expected that the diode cannot lead itself to significantly interrupt the temperature field of the measurement.
Furthermore, the main channel of heat delivery to the pn junction is not its plastic cover, but rather the metal terminals to which the electrical lead wires are soldered.The thermal resistance to energy incoming from the junction to the soldering point is 85 K/W according to the diode performance specifications.Hence, though the diode is too big to be ignored and inevitably affects the characteristics of the temperature field, it is immediately obvious that the temperature sensor here is formed by the metal leads, which are completely immersed in the epoxy matrix.The surface area of the lead wires is comparable to the contact area of a thermocouple that is more commonly used to minimize the sensor influence.It is also obvious that the thermal resistance between the metal leads and epoxy resin will be more or less the same in both cases.
Moreover, temperature measurements are significantly affected by the thermal time constant.The thermal time constant of the system is the time required to come to within 1/e (≈37%) of the fixed value after a step thermal disturbance of the system [119].The thermal time constant of the temperature sensor must be significantly less than the thermal time constant of the sample.Thus, the thermal resistance between the sensor and the sample will not significantly affect the measurement of the sample temperature.The thermal time constant of the diode is where R t is the resistance of the regions next to the junction, and C t = m d C d is the lumped thermal capacitance of the diode with the mass m d and specific heat C d .Taking R t = 85 K/W given above and C d = 800-1200 J•kg −1 •K −1 for the diode molding compound yields τ t on the order of hundreds of ms in the worst case.Essentially, our measurements were carried out in a quasi-static regime such that, after setting the required temperature on the controller, the measuring system waited for dozens of minutes to achieve temperature stabilization.This waiting time depended on the temperature and was set to provide the rate of heating of the sample of less than 2 K/hour.This waiting time exceeded the thermal time constant of the diode.The sensors therefore did not appear to have an appreciable effect on the temperature measurements.
We also implemented a 2D finite element method (FEM) to model the temperature distribution using the COMSOL 5.3a and typical thermal conductivity parameters for diode materials, which are κ = 400 W/mK for copper and κ = 0.58-0.67W/mK for the diode forming mass [121].The simulation results are given in Figures 11 and 12.It is seen in Figure 11 that the temperature field around the diodes is not significantly distorted.Other factors to be taken into account when assessing the accuracy of the measurements are as follows: Copper wires with a thickness of 50 µm and a length of 100 mm were used to feed the diodes.Their thermal resistance is much higher than the thermal resistance of the sample, allowing for a heat leakage to be quenched.The epoxy resin without hardener was used to give a better thermal contact between the sample, heater, and cooler.A plastic screw defining a point contact constriction with a heat-

Conclusions
A theoretical model based on the effective medium approximation was developed to calculate the loading dependence of the thermal conductivity  in hybrid graphenepolymer nanocomposites and evaluate its sensitivity to the geometrical and material parameters of filling particles and the host polymer matrix, filler's mass concentration, self-assembling degree, and Kapitza thermal boundary resistances at emerging interfaces.The hybrid graphene@anatase nanosheets were modeled as sandwich-like structures of either parallelepiped-or spheroid-like configurations, where the geometrical and material characteristics of the anatase layer were treated in a continuous media approximation for two extreme cases representing anatase nanoparticles as either soft or hard solid spheres.The differences between the model results obtained with the above two configurations increased with the decreasing the nanosheet's aspect ratio.The model was also capable of predicting the () dependencies in the nanocomposites if the thermal parameters of the constituent phases and the thermal boundary resistances were known.The soft-sphere approach and related concept of the effective thickness allowed us to simplify the simulations and describe the interface thermal transport in the hybrid graphene@anatase nanosheets via the phenomenological thermal boundary resistance.
We also measured the temperature dependence of the thermal conductivity of epoxy nanocomposites filled with self-assembled hybrid nanoparticles composed of multilayered graphene nanoplatelets and anatase-TiO2 nanoparticles within a temperature range from 40 K to 320 K.When the concentration of anatase in the composites increased, a saturation tendency in the () dependence was observed above ≈150 K, which was possibly due to a saturation of the thermal conductivity of the anatase.A fair agreement between the calculated and measured room-temperature thermal conductivity was obtained.
The numerical results clearly illustrate the thermal screening effect of a lowconductive anatase layer that covers high-conductive graphene nanoplatelets.The observed growing thermal conductivity after raising the temperature is not only related to the increasing intrinsic thermal conductivities of both the fillers and host epoxy but is also subject to temperature-induced variations in the thermal boundary resistances.
An effective tailoring of  is possible only if the anatase-to-graphene mass ratio does not exceed a certain bound, which corresponds to the sheer covering of graphene The temperature variation across the dashed line at x = 0.7 mm in Figure 11c from the bottom to the top sample surface is shown in Figure 12.The dashed line simulates the variation in the absence of the diode sensors, whereas the solid line approximates the one with the embedded diodes.The dashed line in Figure 11c passes through the centers of the diode electrodes, which are heat-sensitive elements in our case.One can therefore see a perturbation in the temperature field at around y ≈ ±1 mm just near the metal electrodes.This is due to the much higher thermal conductivity of the electrode than the epoxy resin.However, at the outer side of the electrodes (at y ≈ ±1.15 mm in Figure 12), the solidand dashed-line temperatures are practically identical.This ensures that the experimental temperature data obtained from these measurements are accurate and reliable when taking into account the thicknesses of the electrodes and when measuring the inter-sensor distance between their external edges.
Other factors to be taken into account when assessing the accuracy of the measurements are as follows: Copper wires with a thickness of 50 µm and a length of 100 mm were used to feed the diodes.Their thermal resistance is much higher than the thermal resistance of the sample, allowing for a heat leakage to be quenched.The epoxy resin without hardener was used to give a better thermal contact between the sample, heater, and cooler.A plastic screw defining a point contact constriction with a heat-insulating fastener was used to clamp the sample, which yielded a heat loss of less than 1%.With the above assumptions, the errors in the measured thermal conductivities do not exceed about 5%.
The sample was placed in a closed-cycle cryostat (CS204, Advanced Research Systems) that was used to vary its temperature.The temperature in the cryostat was controlled using a Lake Shore 331S controller through the computer's RS232 digital port.Measurements were carried out in the temperature range from 40 to 300 K.

Conclusions
A theoretical model based on the effective medium approximation was developed to calculate the loading dependence of the thermal conductivity κ in hybrid graphenepolymer nanocomposites and evaluate its sensitivity to the geometrical and material parameters of filling particles and the host polymer matrix, filler's mass concentration, selfassembling degree, and Kapitza thermal boundary resistances at emerging interfaces.The hybrid graphene@anatase nanosheets were modeled as sandwich-like structures of either parallelepiped-or spheroid-like configurations, where the geometrical and material characteristics of the anatase layer were treated in a continuous media approximation for two extreme cases representing anatase nanoparticles as either soft or hard solid spheres.The differences between the model results obtained with the above two configurations increased with the decreasing the nanosheet's aspect ratio.The model was also capable of predicting the κ(T) dependencies in the nanocomposites if the thermal parameters of the constituent phases and the thermal boundary resistances were known.The soft-sphere approach and related concept of the effective thickness allowed us to simplify the simulations and describe the interface thermal transport in the hybrid graphene@anatase nanosheets via the phenomenological thermal boundary resistance.
We also measured the temperature dependence of the thermal conductivity of epoxy nanocomposites filled with self-assembled hybrid nanoparticles composed of multilayered graphene nanoplatelets and anatase-TiO 2 nanoparticles within a temperature range from 40 K to 320 K.When the concentration of anatase in the composites increased, a saturation tendency in the κ(T) dependence was observed above ≈150 K, which was possibly due to a saturation of the thermal conductivity of the anatase.A fair agreement between the calculated and measured room-temperature thermal conductivity was obtained.
The numerical results clearly illustrate the thermal screening effect of a low-conductive anatase layer that covers high-conductive graphene nanoplatelets.The observed growing thermal conductivity after raising the temperature is not only related to the increasing intrinsic thermal conductivities of both the fillers and host epoxy but is also subject to temperature-induced variations in the thermal boundary resistances.
An effective tailoring of κ is possible only if the anatase-to-graphene mass ratio does not exceed a certain bound, which corresponds to the sheer covering of graphene nanoplatelets by anatase nanoparticles.It is furthermore expected that κ will increase further when decreasing the thermal boundary resistances at the anatase-graphene interface.The latter can be realized by using more sophisticated techniques to chemically graft anatase to the graphene surface.However, when the concentration of hybrid nanosheets approaches its percolation thresholds with an increasing anatase-to-graphene concentration ratio, the temperature dependence of the composite's thermal conductivity weakens.It is therefore suggested that varying the filler concentration that covers the graphene nanoplatelets offers interesting opportunities of applied relevance in the area of self-assembling engineering of graphene surfaces.

Nomenclature of the Model Parameters
L X , L Y , L Z linear dimensions of graphene nanoplatelets, nm f filler particle size, nm α = L Z /L X thickness-to-length aspect ratio of graphene nanoplatelets R 2 average radius of anatase nanoparticles, nm h int interphase layer thickness, nm h 1 thickness of the graphene-epoxy interphase layer, nm h 2 thickness of the anatase-epoxy interphase layer, nm

Figure 1 .
Figure 1.Measured temperature dependencies of  for the neat epoxy (open circles) and epoxy-MLG@TiO2 nanocomposites with different mass loading ratios for MLG (  ) and TiO2 ( ):  :  = 0.01:0.01(closed circles) and 0.01:0.05(squares) on lin-lin (a) and log-log (b) scales.Lines in (b) are fitted to a power law () ∝  with the exponent values  .Triangles in (a) are the calculated results at 300 K (for a detailed description, see Section 2.5).They are arbitrarily shifted along the temperature axis for better comparison.

Figure 1 .
Figure 1.Measured temperature dependencies of κ for the neat epoxy (open circles) and epoxy-MLG@TiO 2 nanocomposites with different mass loading ratios for MLG (C f 1 ) and TiO 2 (C f 2 ): C f 1 : C f 2 = 0.01:0.01(closed circles) and 0.01:0.05(squares) on lin-lin (a) and log-log (b) scales.Lines in (b) are fitted to a power law κ(T) ∝ T n with the exponent values n.Triangles in (a) are the calculated results at 300 K (for a detailed description, see Section 2.5).They are arbitrarily shifted along the temperature axis for better comparison.

Figure 2 .
Figure 2. A log-lin plot of the thermal conductivity of amorphous SiO2, calculated by using the data obtained by Goodson et al. [65].

Figure 2 .
Figure 2. A log-lin plot of the thermal conductivity of amorphous SiO 2 , calculated by using the data obtained by Goodson et al. [65].

Figure 3 .
Figure 3. Regular geometrical architectures of hybrid MLG@TiO2 nanosheets configured in a sandwich-like structure: (a) parallelepiped-shaped graphene nanoplatelet covered entirely with hard spherical anatase nanoparticles; (b) parallelepiped-shaped graphene nanoplatelet covered partially with hard spherical anatase nanoparticles; (c) geometrical effect of tight packing of hard spheres on a flat surface; (d) soft-sphere approach to model anatase wrapping of graphene nanoplatelet as a solid layer with the effective thickness of ℎ .

Figure 3 .
Figure 3. Regular geometrical architectures of hybrid MLG@TiO 2 nanosheets configured in a sandwich-like structure: (a) parallelepiped-shaped graphene nanoplatelet covered entirely with hard spherical anatase nanoparticles; (b) parallelepiped-shaped graphene nanoplatelet covered partially with hard spherical anatase nanoparticles; (c) geometrical effect of tight packing of hard spheres on a flat surface; (d) soft-sphere approach to model anatase wrapping of graphene nanoplatelet as a solid layer with the effective thickness of h 21 .

*
The values of r 21 without the brackets correspond to r 2 without the brackets.** The values in brackets represent appropriate r 2 and r 21 values.

Figure 9 .
Figure 9. (a) SEM image of anatase TiO2 particles decorating graphene nanosheets, (b) TEM image of the graphene flake with TiO2 particles on it marked by dotted circles, and (c) high-resolution TEM image of a hybrid graphene flake@anatase TiO2 nanocomposite.The scale bars are 500 nm in (a) and 50 nm in (b).

Figure 9 .
Figure 9. (a) SEM image of anatase TiO 2 particles decorating graphene nanosheets, (b) TEM image of the graphene flake with TiO 2 particles on it marked by dotted circles, and (c) high-resolution TEM image of a hybrid graphene flake@anatase TiO 2 nanocomposite.The scale bars are 500 nm in (a) and 50 nm in (b).

Figure 9 .
Figure 9. (a) SEM image of anatase TiO2 particles decorating graphene nanosheets, (b) TEM image of the graphene flake with TiO2 particles on it marked by dotted circles, and (c) high-resolution TEM image of a hybrid graphene flake@anatase TiO2 nanocomposite.The scale bars are 500 nm in (a) and 50 nm in (b).

Figure 10 .
Figure 10.Schematics of the setup used for the thermal conductivity measurements.A and B are the temperature sensors embedded into the sample.L is the distance between the sensors.Thermostat is formed by mounting the top sample edge onto the cold head of the cryostat used (T-controlled heat sink).

Molecules 2023 ,Figure 11 .
Figure 11.Simulated temperature distribution in the sample plane across the diode sensors in full color (left-hand images) and contour (right-hand images) plots.(a,b)-without the diodes; (c,d)with the temperature sensors embedded into the sample (see Figure10).The upper sample surface is heated.The dashed line at  = 0.7 mm in (c) shows the temperature distribution given in Figure12.

Figure 11 .
Figure 11.Simulated temperature distribution in the sample plane across the diode sensors in full color (left-hand images) and contour (right-hand images) plots.(a,b)-without the diodes;(c,d)-with the temperature sensors embedded into the sample (see Figure10).The upper sample surface is heated.The dashed line at x = 0.7 mm in (c) shows the temperature distribution given in Figure12.

Figure 12 .
Figure 12.The temperature variation across the dashed line at  = 0.7 mm shown in Figure 11c with (solid line) and without (dashed line) the diode sensors.

Figure 12 .
Figure 12.The temperature variation across the dashed line at x = 0.7 mm shown in Figure 11c with (solid line) and without (dashed line) the diode sensors.

Table 1 .
Physical parameters of the constituent phases.
2 •K η m, f = R K κ m, f Kapitza length (subscript m− matrix, f − filler), nm C 1 mass concentration of graphene nanoplatelets in the composite C 2 mass concentration of anatase nanoparticles in the composite ϕ n volume concentration of constituent phases (subscripts n = 0, 1, 2, and 3 correspond to the polymer matrix, free MLG nanoplatelets, free anatase nanoparticles, and hybrid MLG@anatase nanoparticles, respectively); see Equations (S39)-(S45) (Supplementary Note S2) ϕ i1 volume portion of the graphene-epoxy interphase layer ϕ i2 volume portion of the anatase-epoxy interphase layer ϕ i3 volume portion of the interphase layer surrounding graphene@anatase hybrid nanoparticles ϕ i30 volume portion of an anatase in graphene@anatase hybrid nanoparticles Φ n functions defined by Equations (S1)-(S4) (Supplementary Note S1) F 1,2,3 volume factors, which determine the ratio of the interphase layer's volume to the particle's volume p 1,2 mass portions of free (unassembled) fillers, graphene (subscript 1), and anatase (subscript 2) ρ n mass densities of the polymer matrix (n = 0), graphene nanoplatelets (n = 1), and anatase nanoparticles (n = 2), kg•m −3 ρ i1,2 mass densities of the graphene-epoxy (subscript 1) and anatase-epoxy (subscript 2) interphase layers m n , v n , s n masses (in kg), volumes (in m 3 ), and surface areas (in m 2 ) of graphene nanoplatelets (n = 1) and anatase nanoparticles (n = 2) S ij (α) two-rank Eshelby tensor N 21 total number of anatase nanoparticles per single graphene nanoplatelet N 0 maximum number of anatase nanoparticles that can be packed tightly over the surface of a single graphene nanoplatelet NL 2 hypothetical number of layers of anatase particles covering graphene nanoplatelets H 21,e f f effective thickness of an anatase covering graphene nanoplatelets in the "hard-spheres" approximation, nm R G radius of gyration of the polymer matrix g d grafting density