Coarse-grained models have increasingly been used in large-scale particle-based simulations. However, due to their lack of degrees of freedom, it is a priori unlikely that they straightforwardly represent thermal properties with the same accuracy as their atomistic counterparts. We take a first step in addressing the impact of liquid coarse-graining on interfacial heat conduction by showing that an atomistic and a coarse-grained model of water may yield similar values of the Kapitza conductance on few-layer graphene with interactions ranging from hydrophobic to mildly hydrophilic. By design the water models employed yield similar liquid layer structures on the graphene surfaces. Moreover, they share common vibration properties close to the surfaces and thus couple with the vibrations of graphene in a similar way. These common properties explain why they yield similar Kapitza conductance values despite their bulk thermal conductivity differing by more than a factor of two.

The Kapitza conductance G (also called thermal boundary conductance or interfacial thermal conductance) and the reciprocal quantity called Kapitza resistance are important to quantify heat transfer through interfaces at the nanoscale. Knowledge of the underlying physical principles that determine G is of great relevance in many practical situations where interfaces are present, like in nanofluids or nanocomposites. A typical technological problem that requires better understanding of interfacial heat transfer is the difficulty to dissipate heat in modern microelectronic devices. As Balandin put it, “Heat is one of the worst enemies of electronics.”1 Numerous experimental and theoretical studies as well as molecular dynamics (MD) simulations have been carried out to tackle such problems. In this context, the interfacial thermal behavior of water on solid substrates has received special attention. Such systems can be seen as prototypes that can be addressed through realistic models in particle-based simulations to better understand the mechanisms of heat transfer between hard and soft materials. We briefly discuss below the related findings on which our work builds. Interested readers may refer to Refs. 2 and 3 or 4 and references therein for a broader overview.

It was observed on surfaces decorated with self-assembled monolayers with different head groups that G generally increases with the water-solid interaction strength approximately through a linear relationship between G and the work of adhesion WSL.5–8 This quantity is connected to the contact angle θ, which is a visible macroscopic measure of the affinity of a given liquid for a given surface. For nonpolar surfaces, WSL and θ are connected through WSL = γ(1 + cos θ), where γ is the liquid surface tension. In MD simulations, such a linear relationship was found on solids like graphite or silicon whose interaction strength with water was tuned through force-field parameters.2,3 However, it was also found that G strongly depends on the thickness of the substrate2 and on its crystallographic orientation.9 

Recent studies on nonpolar surfaces have demonstrated that G is strongly influenced by the water molecules in the immediate vicinity of the substrate.2,9 In this region, water adopts a layer structure with oscillations in the number density perpendicular to the surface. However, the layer structure tends to disappear when WSL decreases.10 Some of the water properties within the layers like orientation dynamics11 and local compressibility12,13 are altered compared to their bulk counterparts. In fact, local compressibility in the layers decreases when WSL increases. Thus, although these water layers are liquid, changes induced by the interaction with the surface indicate rigidification compared with bulk water. In a contribution inspired by a phonon theory of liquids,14,15 Mérabia et al. noted that “interfacial energy transport at solid/water interfaces is to a very large extent mediated by high-frequency vibrations, for which liquid water behaves as a solid.”4 Note that these authors developed their approach based on bulk water vibration properties. However, the special role played by the high-frequency vibrations they mention suggests that the increased rigid character of the water structure within the adsorbed layers associated with the rise of WSL may promote interfacial heat conduction.

In a recent work,16 we showed that single-site coarse-grained (CG) models may also adopt the layer structure of water next to nonpolar surfaces like few-layer graphene. More precisely, we found that the atomistic model SPC/E and the CG model mW described below yield very similar layer structures and values of WSL on graphene surfaces upon varying the strength of the water-carbon interaction. In light of the reasoning above that connects the water layer structure on surfaces with WSL and G, we hypothesize that mW might yield variation of G against WSL similar to SPC/E despite the missing atomistic details and the corresponding reduced number of degrees of freedom. The validity of this assumption is demonstrated and rationalized below.

We performed simulations in which water was represented by the SPC/E17 and mW18 models. While SPC/E is a rigid atomistic model with the three atoms of water, mW is a single-site single-molecule CG model. Hydrogen bonds are explicitly represented in SPC/E through electrostatic interactions. As for mW, Molinero and Moore mentioned that “mW mimics the hydrogen-bonded structure of water through the introduction of a nonbonded angular dependent term (i.e., a Stillinger-Weber three-body potential19) that encourages tetrahedral configurations.” Graphene was modeled with atomistic details following Alexeev et al.2 Water and graphene interacted through a 12-6 Lennard-Jones potential. The wetting strength was changed by tuning the potential energy parameter εSL to obtain work of adhesion values ranging from 6 to about 100 mJ/m2 (see the supplementary material for detailed information). Note that WSL ≈ 100 mJ/m2 corresponds to θ ≈ 65°, which is a value typically obtained in experiments on highly ordered pyrolytic graphite that take surface contamination by airborne compounds into account.20 On monolayer graphene free of contamination, θ = 85° ± 5° (WSL = 78 ± 6 mJ/m2).21 Note that the carbon surfaces with the lowest values of WSL in our study should be regarded as model systems employed to enhance our understanding.

The Kapitza conductance between water and few-layer graphene was computed at an average temperature of 300 K and a pressure of 1 atm using the MD simulation package LAMMPS.22 The non-equilibrium MD (NEMD) setup employed was similar to the implementation used in Ref. 2 with 12 carbon layers on each side of the simulation cell, as shown in Fig. S1 of the supplementary material. Further information on the simulation procedure is given in the supplementary material. The thermal conductivity k of water in the sandwiched films was evaluated as the ratio of the heat flux Q and the temperature gradient in the film. For both SPC/E and mW at all WSL, the values obtained were in quantitative agreement with the corresponding liquid bulk thermal conductivity simulations (see below). Finally, it should be mentioned that the thermal conductance measured by the setup above yields the combined conductance of the graphene-graphene and graphene-water interfaces as explained in Ref. 23. Further, the absolute values of thermal properties may depend on the force field and the method employed.24 However, these facts should not alter the conclusions we draw as we compare two models with the exact same setup. The way the number of graphene layers influences G should be identical for both the studied water models, and what we refer to as G in fact accounts for the contribution of the solid-liquid and the inter-layer graphene conductances.

We briefly discuss some bulk and interfacial structure properties of SPC/E and mW that are relevant for the question we address. Both models yield very similar densities (990 kg/m3 for SPC/E and 998 kg/m3 for mW). They also lead to the typical short-range tetrahedral order of water that arises from hydrogen bonds (not shown). Owing to these similarities, it can be seen in Fig. 1 that SPC/E and mW lead to comparable distributions of the number density of water on graphene layers with the typical liquid layer structure perpendicular to the interface. When the water-surface coupling decreases, the water layer structure gradually disappears for both models and vanishes for WSL ≈ 15–20 mJ/m2. At WSL ≈ 6 mJ/m2 in Fig. 1, no layer structure is visible for SPC/E, while the maximum at about 0.45 nm for mW is not representative of the solid-liquid interaction but of the liquid-vapor interface. When water interacts with nonpolar substrates like graphene, tetrahedral order in the first layer of adsorbed molecules is perturbed due to the surface. However, for both water models, we found that molecules in the second adsorbed layer and beyond adopt short-range tetrahedral order like in the liquid bulk, independently of the solid-liquid interaction strength. Thus, despite the difference in the number of degrees of freedom between SPC/E and mW, these models adapt their structure in response to the presence of a nonpolar surface in a similar way. Furthermore, both models yield very similar structures of water at the interface. In contrast, the lack of degrees of freedom in mW strongly affects the properties related to thermal transport. We found through equilibrium simulations (see the supplementary material) that the isobaric specific heat capacity of SPC/E and mW was 84.9 J/mol K and 33.2 J/mol K, respectively. Both models yield different values of the thermal conductivity k. It was found through NEMD simulations that k = 0.83 ± 0.04 W/Km for SPC/E and k = 0.34 ± 0.01 W/Km for mW (see the supplementary material).

FIG. 1.

Number density distribution of water on graphite (located at z ≤ 0) calculated at two different wetting strengths for both the models. The coordinate z is the distance from the surface.

FIG. 1.

Number density distribution of water on graphite (located at z ≤ 0) calculated at two different wetting strengths for both the models. The coordinate z is the distance from the surface.

Close modal

We now discuss the variation of G for SPC/E and mW depending on WSL. It can be seen in Fig. 2 that G increases with WSL for both SPC/E and mW. The variation is linear for WSL > 40 mJ/m2. These observations are in line with several earlier studies on solid-liquid systems.2,8,9,25 Moreover, the numerical magnitude of G is in close proximity to the values of Alexeev et al.2 It can also be observed that below WSL = 40 mJ/m2 the rate with which G increases is slightly larger than above, in agreement with the finding of Xue et al. on the effect of the wetting regime on 1/G.25 The most striking result is that SPC/E and mW lead to variations of G against WSL that are quantitatively similar. Note that G seems to diverge between the two models at the greatest WSL values, although the effect is weak. In fact, additional simulations in the fully wetting regime (θ = 0°) at WSL ≈ 230 mJ/m2 showed that even for such a strong solid-liquid interaction, the difference in G between SPC/E and mW is similar to the difference at 100 mJ/m2 (note that G ≈ 45 MW/m2 K at that WSL value). The result in Fig. 2 is therefore in strong contrast with k for the two models and is in line with the idea that interfacial properties like G cannot be trivially obtained by combining individual bulk properties.26 

FIG. 2.

Interfacial thermal conductance depending on the solid-liquid work of adhesion for SPC/E (red) and mW (green).

FIG. 2.

Interfacial thermal conductance depending on the solid-liquid work of adhesion for SPC/E (red) and mW (green).

Close modal

To understand the behavior of G with respect to WSL, a spectral analysis was performed. We first discuss the equilibrium bulk and interfacial vibration properties of the water models through their vibration spectrum obtained from

Pω=imiṽiω2,
(1)

where the summation is performed over all liquid particles, m is their respective mass, and ṽiω is the Fourier-transform of the velocity time-series. Further information on the calculation of P(ω) is given in the supplementary material. It can be seen in Fig. 3 that SPC/E has a bulk spectrum that spans a broader frequency range than mW. The atomistic model in the bulk yields a relatively narrow band centered on 1.5 THz and a broader one centered on 15 THz. These bands correspond to translational and rotational librations, respectively.27 When SPC/E interacts with the surface, the maximum at 1.5 THz is shifted to about 2.5 THz while the band for the rotational librations is slightly red-shifted. Due to the presence of the surface, molecular packing in the first adsorbed layer of water is tighter and rotational motions are hampered. Consequently, the magnitude of P(ω) for the rotational librations is underrepresented compared to the translational ones. In addition to a band at 1 THz for translational librations, mW in the bulk yields another narrow band centered on 5 THz, which presumably arises from the three-body interaction potential. Like SPC/E, the presence of the surface markedly influences the spectrum of mW. The bulk band at 1 THz is blue-shifted to 2 THz, while the intensity of the band at 5 THz slightly decreases.

FIG. 3.

Vibration spectrum for bulk water for SPC/E (blue) and mW (green) along with the vibration spectrum of the water molecules in the first adsorbed layer for SPC/E (red) and mW (black) at WSL ≈ 100 mJ/m2. All spectra were normalized to ease the comparison.

FIG. 3.

Vibration spectrum for bulk water for SPC/E (blue) and mW (green) along with the vibration spectrum of the water molecules in the first adsorbed layer for SPC/E (red) and mW (black) at WSL ≈ 100 mJ/m2. All spectra were normalized to ease the comparison.

Close modal

The water vibrations on the surface in the region between 2 and 2.5 THz turn out to be of crucial importance as far as G is concerned, as revealed by the frequency-dependent heat flux q(ω) at the interfaces. Note that the frequency integral of q(ω) yields the heat flux Q used to calculate G.28 That quantity was characterized as follows:29,28

qω=2AδtiF̃iωṽi*ω,
(2)

where A is the cross-sectional area, δt is the time interval (2 fs here) between two consecutive snapshots, and the quantities in the bracket refer to the Fourier transforms of forces and velocities. The bracket refers to a non-equilibrium average. When the summation is performed on the solid atoms (index i), the force in Eq. (2) is the total force on the i-th solid atom due to the liquid particles. Conversely, the summation may be performed on liquid particles. In this case, the force on the i-th particle arises from interactions with the solid atoms. We show in Fig. 4(a) the variation of q(ω) for mW and SPC/E with the summation performed on the carbon atoms at WSL ≈ 100 mJ/m2. The coupling between both models and the carbon atoms around 2–2.5 THz can clearly be seen. Further coupling between water models and the graphene vibrations between 20 and 40 THz can be observed. The phonon density of states of the graphene layers is reported in Fig. S4 (supplementary material), where it can be seen that substantial density of states comprising the acoustic out-of-plane modes exists in this spectral region. Strikingly, both water models’ vibrations couple in a similar way with graphene. This observation explains why G takes similar values for mW and SPC/E at WSL ≈ 100 mJ/m2. The strong effect of the water vibrations around 2–2.5 THz on G is even more clearly visible in Fig. 4(b) (red and green curves) where the flux is calculated by performing the summation on the liquid particles in Eq. (2). The curves for q(ω) for both models are similar, with a slight frequency increase for SPC/E, as discussed above regarding P(ω). When WSL decreases to ≈6 mJ/m2, the maximum in q(ω) appears at lower frequencies for both models [Figs. 4(b) and 4(c)]. In this weak wetting regime, the water layer structure has vanished (Fig. 1); the interface then closely resembles a liquid-vapor-like interface. The collisions of the water molecules with the surface occur less frequently than for larger WSL values. The frequency integral of the curves for q(ω) in Figs. 4(a) and 4(c) reveals that the energy exchange at low frequency vibrations increases as WSL decreases (see Fig. S5 in the supplementary material).

FIG. 4.

The normalized frequency-dependent heat flux calculated from the carbon atoms at WSL ≈ 100 mJ/m2 and WSL ≈ 6 mJ/m2 is presented in (a) and (c), respectively, for SPC/E (red) and mW (green). (b) Normalized frequency-dependent heat flux calculated from water WSL ≈ 100 mJ/m2 (red for SPC/E and green for mW) and WSL ≈ 6 mJ/m2 (blue for SPC/E and black for mW).

FIG. 4.

The normalized frequency-dependent heat flux calculated from the carbon atoms at WSL ≈ 100 mJ/m2 and WSL ≈ 6 mJ/m2 is presented in (a) and (c), respectively, for SPC/E (red) and mW (green). (b) Normalized frequency-dependent heat flux calculated from water WSL ≈ 100 mJ/m2 (red for SPC/E and green for mW) and WSL ≈ 6 mJ/m2 (blue for SPC/E and black for mW).

Close modal

To conclude, the general mechanism that leads G to increase with WSL (Fig. 2) may be summarized as follows. When the interaction between water and graphene increases, the packing of water molecules in the vicinity of the surface is enhanced. From the vibration point of view, tighter packing results in the increase of the frequency of the water translational librations, yielding enhanced coupling with graphene. In turn, the enhanced vibrational coupling between water and graphene yields an increase in G. Even though mW and SPC/E differ in their degrees of freedom, they lead to very similar interfacial water structures on graphene at each value of WSL that covers a broad range from hydrophobic to mildly hydrophilic. This common behavior concerning structure presumably results from the representation of short-range tetrahedral order arising from hydrogen bonds by both models although through different interaction potentials. Water adsorption on the surfaces influences the vibration properties of the models in a comparable way. Owing to the fact that both SPC/E and mW vibration properties behave similarly next to the substrate at all interaction couplings, G increases with WSL in a similar way for both models. Our work demonstrates the close connection between the interfacial liquid structure and Kapitza conductance. Furthermore, it takes a first step in addressing the effect of coarse-graining on interfacial heat transfer by showing that liquid CG models may in fact be able to reproduce the behavior of atomistic models under specific circumstances regarding their structure behavior. Such CG models may prove useful to address the interfacial thermal properties of large-scale systems that are out of reach of atomistic simulations. For the systems considered here, the CG simulations were roughly twice as fast as the atomistic ones with the computation time being strongly influenced by the carbon-carbon interactions.

See supplementary material for information on simulation protocols, system sizes (Sec. I), interaction parameters (Sec. II), details on the evaluation of temperature discontinuities (Secs. III and IV), details on bulk thermal conductivity and heat capacity evaluation (Sec. V), details on power spectrum evaluation (Secs. VI and VII), and cumulative heat flux corresponding to Figs. 4(a) and 4(b) (Sec. VIII).

V.R.A. thanks Professor Dr. Florian Müller-Plathe and Dr. Samy Mérabia for fruitful discussions and Dr. Haoxue Han for useful discussions. This research was funded by the German Research Foundation (DFG) through the SFB TRR 146 Multiscale Simulation Methods for Soft Matter Systems. We acknowledge the HPC Center at TU Darmstadt for computer time on the Lichtenberg Cluster and the John von Neumann Institute for Computing (NIC) at Jülich Supercomputing Center for computer time on JURECA.

1.
2.
D.
Alexeev
,
J.
Chen
,
J. H.
Walther
,
K. P.
Giapis
,
P.
Angelikopoulos
, and
P.
Koumoutsakos
,
Nano Lett.
15
,
5744
(
2015
).
3.
B.
Ramos-Alvarado
,
S.
Kumar
, and
G. P.
Peterson
,
J. Phys. Chem. Lett.
7
,
3497
(
2016
).
4.
S.
Merabia
,
J.
Lombard
, and
A.
Alkurdi
,
Int. J. Heat Mass Transfer
100
,
287
(
2016
).
5.
N.
Shenogina
,
R.
Godawat
,
P.
Keblinski
, and
S.
Garde
,
Phys. Rev. Lett.
102
,
156101
(
2009
).
6.
Z.
Ge
,
D. G.
Cahill
, and
P. V.
Braun
,
Phys. Rev. Lett.
96
,
186101
(
2006
).
7.
H.
Acharya
,
N. J.
Mozdzierz
,
P.
Keblinski
, and
S.
Garde
,
Ind. Eng. Chem. Res.
51
,
1767
(
2012
).
8.
H.
Harikrishna
,
W. A.
Ducker
, and
S. T.
Huxtable
,
Appl. Phys. Lett.
102
,
251606
(
2013
).
9.
B.
Ramos-Alvarado
and
S.
Kumar
,
J. Phys. Chem. C
121
,
11380
(
2017
).
10.
F.
Sedlmeier
,
J.
Janecek
,
C.
Sendner
,
L.
Bocquet
,
R. R.
Netz
, and
D.
Horinek
,
Biointerphases
3
,
FC23
(
2008
).
11.
G.
Stirnemann
,
P. J.
Rossky
,
J. T.
Hynes
, and
D.
Laage
,
Faraday Discuss.
146
,
263
(
2010
).
12.
R.
Evans
and
M. C.
Stewart
,
J. Phys.: Condens. Matter
27
,
194111
(
2015
).
13.
S. N.
Jamadagni
,
R.
Godawat
, and
S.
Garde
,
Annu. Rev. Chem. Biomol. Eng.
2
,
147
(
2011
).
14.
J.
Frenkel
, in
Kinetic Theory of Liquids
, edited by
R. H.
Fowler
,
P.
Kapitza
, and
N. F.
Mott
(
Oxford University Press
,
1947
).
15.
D.
Bolmatov
,
V. V.
Brazhkin
, and
K.
Trachenko
,
Sci. Rep.
2
,
1
(
2012
).
16.
V. R.
Ardham
and
F.
Leroy
,
J. Chem. Phys.
147
,
074702
(
2017
).
17.
H. J. C.
Berendsen
,
J. R.
Grigera
, and
T. P.
Straatsma
,
J. Phys. Chem.
91
,
6269
(
1987
).
18.
V.
Molinero
and
E. B.
Moore
,
J. Phys. Chem. B
113
,
4008
(
2009
).
19.
F. H.
Stillinger
and
T. A.
Weber
,
Phys. Rev. B
31
,
5262
(
1985
).
20.
A.
Kozbial
,
C.
Trouba
,
H.
Liu
, and
L.
Li
,
Langmuir
33
,
959
(
2017
).
21.
T.
Ondarçuhu
,
V.
Thomas
,
M.
Nuñez
,
E.
Dujardin
,
A.
Rahman
,
C. T.
Black
, and
A.
Checco
,
Sci. Rep.
6
,
24237
(
2016
).
22.
S.
Plimpton
,
J. Comput. Phys.
117
,
1
(
1995
).
23.
L.
Hu
,
T.
Desai
, and
P.
Keblinski
,
Phys. Rev. B
83
,
195423
(
2011
).
24.
M.
Zhang
,
E.
Lussetti
,
L. E. S.
de Souza
, and
F.
Müller-Plathe
,
J. Phys. Chem. B
109
,
15060
(
2005
).
25.
L.
Xue
,
P.
Keblinski
,
S. R.
Phillpot
,
S.
Choi
, and
J. A.
Eastman
,
J. Chem. Phys.
118
,
337
(
2003
).
26.
K.
Gordiz
and
A.
Henry
,
New J. Phys.
17
,
103002
(
2015
).
27.
S. T.
Lin
,
P. K.
Maiti
, and
W. A.
Goddard
 III
,
J. Phys. Chem. B
114
,
8191
(
2010
).
28.
K.
Sääskilahti
,
J.
Oksanen
,
J.
Tulkki
, and
S.
Volz
,
Phys. Rev. B
90
,
134312
(
2014
).
29.
A.
Giri
,
J. L.
Braun
, and
P. E.
Hopkins
,
J. Phys. Chem. C
120
,
24847
24856
(
2016
).

Supplementary Material