Reducing Kapitza resistance of graphene–paraffin interfaces by alkyl functionalisation

was previously shown possible in the liquid state. For application in heat storage, the thermal properties in solid state are equally important. In the present study, the most optimal reduction of ITR in the solid state was shown for octyl functional groups; in both molten and solid samples, a striking reduction by 95% as compared to pristine graphene was achieved.


H I G H L I G H T S
• Pulse method was found better suitable for ITR measurements in solid state.• In solid state, there is no difference in ITR between octyl and hexadecyl.• Possible reason for above is orientation of the functional groups.• Reduction of ITR up to 96% was achieved.

A B S T R A C T
Paraffin waxes are a promising material for heat storage.Improvement of the thermal conductivity can be obtained by adding graphene nanofillers, but interfacial thermal resistance (ITR) forms a problem.Reduction of ITR by functionalisation of the graphene was previously shown possible in the liquid state.For application in heat storage, the thermal properties in solid state are equally important.In the present study, the most optimal reduction of ITR in the solid state was shown for octyl functional groups; in both molten and solid samples, a striking reduction by 95% as compared to pristine graphene was achieved.

Introduction
In the global search for more environmentally friendly ways of living, a great challenge lies in heating of spaces we work and live in.A promising method of heating, that can contribute to pollutionfree generation of warmth, is the application of heat buffers in heat pumps [1].In this way, heat pumps can be active at times of plenty of renewable energy, and the produced heat can be withdrawn from the buffer whenever desirable.
Paraffin waxes, which are abundantly available at low cost, have high potential for heat storage.Their large latent heat [2] and a melting temperature around human body temperature make them applicable for use in heat buffers [3].By choice of the (average) molecular weight, the best matching melting temperature for the targeted application can be obtained.This makes paraffin suitable for buffers used for heating as well as buffers used for cooling.
However, heat needs to flow into or out of heat buffers at high power in order to facilitate temperature control of domestic spaces.The low thermal conductivity (TC) of paraffin waxes of roughly 0.2 W/(m K) [4] makes charging and discharging of heat buffers containing pure paraffin inapplicable.Addition of highly thermally conductive graphene nanofillers could solve this problem.Multiple studies have shown significantly larger TC of wax-graphene composites as compared to the TC of pure paraffin [5,6] and the improvement of TC has grown over the recent years.
The solution to the low TC by introducing graphene nanofillers to the paraffin comes with a problem of its own.Where heat can flow through graphene easily, thereby aiding the transport of heat to/from large amounts of wax, the interfacial thermal resistance (ITR) acts as a bottleneck [7].Previous studies into reduction of ITR have shown that this can be achieved by functionalisation of the graphene with alkyl functional groups giving very good results in liquid paraffin [8,9].
To the best of our knowledge, no research was done so far on ITR reduction in solid paraffin.This is equally important as the reduction of ITR of molten wax for successful use of paraffin in heat buffers, in which melting and crystallisation take place during charging and discharging, respectively.In the present study, the effect of alkyl functionalisation of graphene on ITR with paraffin wax was investigated by molecular dynamics simulation.Firstly, results from two measuring techniques are compared for validation, showing good correspondence between the two in liquid phase and higher suitability of one for use in the solid phase.Secondly, the influence of alkyl length on ITR is studied in both molten and solid paraffin.Thirdly, the variation of area density of octyl functionalisation of the graphene was looked into.Lastly, a possible reason for the differences of ITR between molten and solid wax was examined.The highest reduction of ITR achieved in both liquid and solid state was a striking 95%, with octyl functionalisation at 0.23 nm −2 area density.ITR of paraffin-graphene composites with pristine graphene was found to be (8.1 ± 0.4) ⋅ 10 −9 m 2 K/W and (13.3 ± 0.5) ⋅ 10 −9 m 2 K/W in liquid and solid state respectively.
The structure of the paper is as follows.In Section 2, the details of simulations and measurements are revealed.In Section 3 the main results are presented and discussed.Finally, the conclusions and outlook are formulated in Section 4.

Models and methods
In this study, the focus is on the interactions between paraffin and graphene at the interface, specifically the Interfacial Thermal Resistance (ITR), also known as Kapitza resistance [10,11].As shown by Nazarychev et al. [12], thermal transport in paraffin waxes is modelled best by the fully atomistic General Amber Force Field (GAFF) [13].The optimised Tersoff force field [14], which gives the most realistic results for thermal conductivity within graphene, can only be used within a layer of graphene and not for the interaction between graphene and paraffin.Given that it is purely the heat transport through the interface that we are interested in, GAFF was also used for the graphene atoms as well as for the alkyl groups.A full overview of force-field parameters is provided in the Supplemental Information.
The paraffin wax in the present study consisted purely of monodisperse eicosane, C 20 H 42 .The modelled pristine graphene sheets measuring 9.9 × 9.8 nm 2 contained 3680 atoms.To prevent folding of the graphene, the x-and y-dimensions of the simulation boxes were kept constant and only the length in z-direction was allowed to change.All boundary conditions of the simulation box were periodic.To the graphene, alkyl chains of various lengths and in various amounts were attached.At a constant area density of 0.26 nm −2 , alkyl groups of one type chosen from ethyl C 2 H 5 , butyl C 4 H 9 , octyl C 8 H 17 or hexadecyl C 16 H 33 , were attached.Furthermore, solely with octyl functionalisation, the area density was varied from 9.3 ⋅ 10 −2 nm −2 to 2.3 nm −2 , the latter of which is larger than common simulation and experimental values [5,8,15].
Eicosane chains as well as (functionalised) graphene sheets were prepared using AmberTools [16] after which they were converted to GROMACS [17] format using AnteChamber PYthon Parser interfacE (ACPYPE) [18].A single molecule of eicosane was replicated to acquire the desired number, leading to a very ordered structure as described in our previous study [19].Two graphene sheets were combined with roughly 1800 eicosane molecules, half of which between and the other half outside the graphene sheets to give a sample symmetric in zdirection.The initial order of the eicosane was disposed of by heating the samples up to 800 K for 1 ns in the NVT ensemble using the Berendsen thermostat.This was followed by equilibration at 450 K, which is well above the melting temperature of eicosane [4], and at 1 atm pressure in the NPT ensemble for 1 ns with 1 fs timesteps.Samples with the wax in semi-crystalline state were produced by slow cooling of the composites at a rate of 1 K/ns, which was shown to be sufficiently slow for crystallisation to occur, in the recent study of Nazarychev et al. [20].To speed up cooling simulations, these were executed with timesteps of 2 fs and with use of the LINCS [21] constraint algorithm for the C-H bonds.Whereas GROMACS was used for the equilibration and cooling of the samples because of its high speed, LAMMPS [22] was used for measurements due to its integrated functions.Conversion of the simulational files between these software packages was done by use of InterMol [23].Prior to measurements, the samples were re-equilibrated in the NVT with for 100 ps using the Nosé-Hoover thermostat.This equilibration as well as the measurements were run with 0.5 fs timesteps and without any constraints.
Multiple methods are available to measure Kapitza resistance, e.g. the relatively new one developed by Alosious et al. [24][25][26].Two wellestablished techniques, extensively described by e.g.Song et al. [27] were used in the present study.Firstly, the Reverse Non-Equilibrium Molecular Dynamics (RNEMD) developed by Müller-Plathe [28], normally used for measurement of thermal conductivity, was applied to create a temperature difference between the centre and bottom of the simulation box as illustrated in Fig. 1.In this method, the temperature difference is caused by swapping of velocities of atoms.Therefore, once steady state is reached, the heat flux is exactly known by with  [W/m 2 ] the heat flux, v ℎ and v  [m/s] the velocities of the hottest and coldest particles with mass  [kg] at each transfer,  [s] the time of measurement in steady state and  the area [m 2 ] through which the heat flux flows, in this case the interface is parallel to the xy-plane.Subsequently, the temperature drop at the interfaces is measured, which is defined as half the temperature drop from the eicosane directly above to that directly below the graphene.The Kapitza resistance can then be calculated using the relation In the RNEMD method, the domain is divided into an even number of slabs, parallel to the graphene sheets (dark green).By swapping of velocities, the bottom slab is cooled down (blue) while the first slab above the centre is heated (red).Because of periodic boundary conditions, a periodic copy of the cold bottom slab is present directly above the simulation box displayed and equal amounts of heat flow to the top and bottom.In the heat pulse method, no division into slabs is implemented.Heating occurs in the graphene sheets including functional groups, i.e. all of the dark green parts of the image.
with   [m 2 K/W] the Kapitza resistance,  the temperature drop over the interface and  the heat flux.As known from literature, the RNEMD method is sensitive to density fluctuations [29].Combined with the high temperature difference needed for a clear temperature profile to develop, this sensitivity made the RNEMD method less suitable for measurements of the crystallised samples at 250 K.It was therefore only used, in the present study, for molten samples at 450 K as a reference to compare the results from the second method to.Secondly, a heat pulse method was used [30,31].In this method, the graphene sheets including functional groups (as shown in Fig. 1) are abruptly heated in the NVT ensemble by temperature rescaling for 1 ps while the rest of the simulation box is kept at 450 K or 250 K in liquid and solid state simulations respectively.Subsequently, in the NVE ensemble the decay of the temperature difference between the graphene and the closest paraffin is measured.The Kapitza resistance is then determined from the time constant of this exponential decay: with  the measured time constant and  the heat capacity of graphene.The latter was determined from fluctuations of internal energy during 200 ps runs in NVT ensemble with 0.5 fs timesteps: with   [J/(kg K)] the heat capacity at constant volume and  [J] the internal energy.For pristine graphene the   value was found to be 1.04 ⋅ 10 3 J/(kg K), which is consistent with literature [32].Data of the functionalised graphene sheets is provided in the Supplemental Information.
The uncertainty of ITR data was determined from the standard deviation of separate ITR measurements of the two graphene sheets.
During measurements, special care should be taken regarding the influence of rigid bonds.As discussed by e.g.Olarte-Plata and Bresme [33] and Sanderson et al. [34], the reduction in degrees of freedom (DoF) by constraints can lead to erroneous temperature values, especially at interfaces.In the present study, no constraints are imposed during measurement simulations.The local temperature values, obtained from the local atomic kinetic energy with the compute ke/atom and fix ave/chunk commands, and used for the calculation of ITR are therefore reliable.In the SI, an example is included of how constraints can lead to large errors of obtained temperatures if the reduction in DoF is not taken into consideration.
Measurements of how far the alkyl chains of functionalised graphene reach into the paraffin were done by post-processing in MAT-LAB.Since the graphene sheets cross the periodic boundary conditions and are (nearly) flat, the distance between these and the free ends of the alkyl chains can easily be measured.The thickness of the paraffin layers was measured in the semi-crystalline samples at 250 K and then applied to the measurements both at 250 K and 450 K.
Screenshots of the simulation boxes were made with VMD [35].

Results and discussion
The results of the present study are presented in three parts.Firstly, reduction of Kapitza resistance in molten samples is discussed.Secondly, the novel subject of Interfacial Thermal Resistance (ITR) in solid samples is investigated.Lastly, the effect of area density of functionalisation on ITR is studied.
In our previous study [36], it was shown that next to a graphene sheet, an aligned layer of paraffin forms in otherwise molten wax at 450 K.Because of the empty space between the graphene and this first layer, as shown in Fig. 2(a), it is expected that thermal transport perpendicularly through the wax-graphene interface will be inhibited.From literature, e.g. the studies by Wang et al. [9] and Liu et al. [8], it is known that alkyl functionalisation of graphene is an effective method to lower the ITR, also known as Kapitza resistance.In Fig. 2(b), the density profile of a composite with a sheet of graphene functionalised with octyl at a very high area density of 2.3 nm −2 is shown.
Comparing the density profiles of the composites with pristine and with functionalised graphene, we see a clear difference.The empty space which is present at the interface with pristine graphene, is filled at the interface with functionalised graphene.In analogy to a bridge enabling traffic to cross a valley, it is expected that the functional groups will increase heat transport between the graphene and nearby paraffin.To study this more closely, two effects were systematically looked into.Firstly, the Kapitza resistance of samples with functional groups of various lengths at a fixed area density was measured.Secondly, while keeping the alkyl chain length constant, the effect of area density on ITR was quantified.The Kapitza resistance of the molten samples was measured with two techniques, as described in Section 2. As shown in Fig. 3, the results of these two methods agree well.Additionally, the Kapitza resistance of pristine graphene obtained from these measurements, being (8.1 ± 0.4) ⋅ 10 −9 and (8.9 ± 0.5) ⋅ 10 −9 m 2 K/W for the pulse and RNEMD techniques respectively, are in agreement with those in literature [9].With increasing length of the alkyl chains, Kapitza resistance decreases.This is in accordance with previous studies, such as by [5,8].Whereas the longest alkyl group included in their studies was butyl, we include octyl and hexadecyl in the present investigation.At 450 K, with the paraffin wax in molten state, the ITR is shown to decrease further with lengthening of the alkyl groups, reduced by half for ethyl and by almost 70% for hexadecyl.
For successful optimisation of graphene-paraffin composites for use in heat buffers, knowledge of ITR in the solid state is essential given its important role in heat transport during charging and discharging of the buffer.As discussed in Section 2, the RNEMD technique was used mainly as a means to verify the data obtained by the pulse technique.In the following measurements only the pulse method is applied.As far as we know, no research was done so far into the ITR of graphene-wax composites in semi-crystalline state.In Fig. 4, the ITR of the samples with varying alkyl length at 250 K are given.The sample with pristine graphene had ITR (13.3 ± 0.5) ⋅ 10 −9 m 2 K/W, which is significantly larger than at 450 K. Contrary to the previous measurements at high temperature, these do not show an ongoing decrease of ITR with increased alkyl length.Instead, a plateau In the sample with pristine graphene (blue solid line), there is a clear absence of material directly next to the graphene sheet.For graphene functionalised with 2.3 nm −2 octyl chains (red dashed line), the extrema are much less pronounced, especially the minimum next to the graphene sheet.It is therefore expected that heat will flow more easily through the interfaces of functionalised graphene.Fig. 3. Relative Kapitza resistance of graphene-eicosane nanocomposites at 450 K in which graphene was functionalised with ethyl, butyl, octyl or hexadecyl, all at area density 0.26 nm −2 .In blue dots the results obtained with the pulse method and in red diamonds obtained with the RNEMD technique are shown.All data are normalised with regards to   of pristine graphene measured by the pulse method, (8.1 ± 0.4) ⋅ 10 −9 m 2 K/W.Results from both methods agree well qualitatively, with a steady decrease of   with increasing alkyl chain-length.For further measurements, it is therefore justified to use only one technique.Fig. 4. Relative Kapitza resistance of semi-crystalline graphene-eicosane nanocomposites at 250 K in which the graphene was functionalised with ethyl, butyl, octyl or hexadecyl, all at area density 0.26 nm −2 .  = (13.3± 0.5) ⋅ 10 −9 m 2 K/W.While at shorter alkyl chain lengths a longer length corresponds to a clearly lower ITR, the results for octyl and hexadecyl are nearly identical.In solid wax, hexadecyl therefore does not have the advantage over octyl that it has in the molten state.
is and there is (almost) no difference between the samples with octyl and hexadecyl functionalisation.
A possible explanation for this unexpected result can be found in the penetration depth of the alkyl chains into the wax.This is illustrated in Fig. 5(a) and quantified in Fig. 5(b).It can be observed that while ethyl chains are only just long enough to reach the layer of paraffin closest to the graphene sheet and butyl reaches well into this layer, octyl extends into the second layer of wax.The hexadecyl chains, however, have a different orientation with many of the chains being ''collapsed'' and parallel to the graphene sheet instead of penetrating more distant paraffin layers.While the hexadecyl chains are longer and therefore have more contact area with surrounding paraffin molecules, the better penetration of octyl into the second layer of wax has the advantage that more heat can flow into this second layer, thereby contributing to heating of the wax layer closest to the graphene from two sides.
Possible reasons for the better penetration of octyl into paraffin as compared to hexadecyl are molecular weight segregation and chain flexibility.As discussed in our previous study [36], molecular weight segregation of paraffin takes place in the vicinity of graphene fillers.It was shown that during slow cooling of a graphene-paraffin composite with bidisperse wax, triacontane (C 30 H 62 ) chains crystallise against the graphene and simultaneously strong segregation takes place, by which the shorter decane (C 10 H 22 ) molecules are moved further away from the graphene.In the samples of the present study, a similar phenomenon is observed.The hexadecyl chains (C 16 H 33 ), which are very close in weight to the eicosane (C 20 H 42 ) chains, solidify at roughly the same temperature as the eicosane wax and therefore partially crystallise as well.The octyl chains however have a much lower weight than the eicosane molecules, meaning they have a much lower melting temperature.The orientation of the octyl perpendicular to the eicosane layers can be considered as molecular weight segregation of these chains which are anchored to the graphene sheet.Another explanation of the different orientation of octyl and hexadecyl can be sought in the flexibility of the chains.Although the freedom of movement of neighbouring carbon atoms in relation to each other in these alkyl groups is equal, the larger length of the hexadecyl chains gives more room for the chain end to bend towards the graphene surface and align itself parallel to it.
The equal reduction of ITR in the semi-crystalline samples with octyl and hexadecyl functionalisation leads to the conclusion that octyl is the preferable functional group since the same effect can be achieved with half the weight percentage of alkyl chains attached to the graphene.Therefore, the effect of the area density of functional groups on ITR will be studied for octyl.In Fig. 6, the Kapitza resistance of samples with various area densities of octyl chains is shown.At area densities within common range [5,8,15], Kapitza resistance decreases steadily with increasing number of octyl chains per unit area graphene surface.The extremely high value of 2.3 nm −2 therefore is included to see if this decrease persists.The strikingly low value of less than 5% of ITR compared to that of pristine graphene shows that this is indeed the case.In the semi-crystalline samples, the same trend is observed; ITR decreases with increasing area density of octyl chains.For all densities, the decrease of ITR as compared to that of the non-functionalised graphene is stronger at 250 K than at 450 K.A possible reason for this lies in the mobility of paraffin molecules in the layers directly adjacent to the graphene.Measurements of mean square displacement show a significantly lower value at 250 K than at 450 K.In molten state there is therefore still some thermal transport by convection while at 250 K there is none.The functional groups thus have a larger relative contribution to heat flow in solid state than in molten state.

Conclusion
Paraffin waxes are a promising material for heat storage.For fast charging and discharging of heat buffers, a much higher TC than that Fig. 6.Relative Kapitza resistance of graphene-eicosane samples at 250 K (blue dots) and 450 K (red diamonds) functionalised with octyl at various area densities.From 0.09 to 0.51 nm −2 the resulting ITR decreases steadily, leading to the question if this decrease will persist at even higher area densities.At 2.3 nm −2 , the ITR has indeed declined to only 6% (liquid) and 4% (solid) of the ITR of pristine graphene in the respective state. of bulk paraffin is needed, which can be achieved by adding graphene nanofillers.Interfacial Thermal Resistance (ITR) forms a bottleneck in the heat transport between graphene and paraffin.Functionalisation of the graphene can reduce ITR, thereby improving the effective TC of the graphene-paraffin composite.
Functional groups of ethyl, butyl, octyl and hexadecyl were attached to graphene sheets in MD simulation boxes otherwise containing eicosane.The area density of the alkyl chains was kept constant.When measuring the resulting Kapitza resistance in molten samples at 450 K, ITR was found to decrease with increasing alkyl chain length.Because the heat transport in solid wax is equally important for the application in heat storage, ITR was also measured at 250 K.In the solid state, after an initial decrease of ITR with increased alkyl length, no significant difference was observed in ITR between the samples with octyl and hexadecyl functional groups.Octyl therefore is the preferable alkyl type given the same ITR reduction at half the weight percentage of functional groups.
The effect of area density of functionalisation with octyl was measured.A steady decrease of ITR with increased area density was observed in both the molten and solid samples.The very large reduction by 94% (liquid) and 96% (solid) was achieved with 2.3 nm −2 octyl functionalisation.Given the higher ITR of samples containing pristine graphene in solid state, (8.1 ± 0.4) ⋅ 10 −9 m 2 K/W and (13.3 ± 0.5) ⋅ 10 −9 m 2 K/W in liquid and solid state respectively, the absolute reduction in the latter is striking.
Possible directions of further research are finite graphene flakes and functionalisation of the edges, as in e.g. the study by Konatham et al. [30].Furthermore, the influence of the type of thermostat and of the use of constraints on crystallisation of the wax could be studied in more depth following the findings of Olarte-Plata and Bresme [33] and of Sanderson et al. [34].

Fig. 1 .
Fig.1.Schematic drawing of a simulation box illustrating the heating techniques used for measurements of Kapitza resistance.In the RNEMD method, the domain is divided into an even number of slabs, parallel to the graphene sheets (dark green).By swapping of velocities, the bottom slab is cooled down (blue) while the first slab above the centre is heated (red).Because of periodic boundary conditions, a periodic copy of the cold bottom slab is present directly above the simulation box displayed and equal amounts of heat flow to the top and bottom.In the heat pulse method, no division into slabs is implemented.Heating occurs in the graphene sheets including functional groups, i.e. all of the dark green parts of the image.

Fig. 2 .
Fig. 2. (a) Graphene-paraffin composite at 450 K, zoomed in on the interface.Directly left and right of the graphene, an aligned layer of paraffin has formed in the otherwise molten eicosane.At the interfaces between the two materials, there is a void.Because of symmetry, two interfaces are depicted; from left to right first the interface from wax to graphene and then the second interface from graphene to wax.(b) Density profile of carbon atoms close to the graphene, averaged over all interfaces.In the sample with pristine graphene (blue solid line), there is a clear absence of material directly next to the graphene sheet.For graphene functionalised with 2.3 nm −2 octyl chains (red dashed line), the extrema are much less pronounced, especially the minimum next to the graphene sheet.It is therefore expected that heat will flow more easily through the interfaces of functionalised graphene.

Fig. 5 .
Fig. 5. (a) Graphene functionalised with ethyl (top left), butyl (top right), octyl (bottom left) and hexadecyl (bottom right).The wax has crystallised into aligned layers parallel to the graphene.While ethyl barely touches the paraffin and butyl reaches into the layer of paraffin closest to the graphene because of their length, octyl extends well into the second layer of wax.For hexadecyl a different orientation is observed; many of the chains have collapsed and do not reach further than the closest layer of paraffin.(b) Quantification of layer reached by the alkyl chain ends.The majority of octyl chains reach into the second layer of paraffin, while of hexadecyl most of the ends are in the layer closest to the graphene.