Thermal conductance of the water–gold interface: The impact of the treatment of surface polarization in non-equilibrium molecular simulations

Interfacial thermal conductance (ITC), G = J_q/ΔT, or its inverse, Kapitza resistance,^1,2 quantifies the abrupt temperature drop, ΔT, associated with a heat flux, J_q, normal to an interface between dissimilar materials or fluids. The first study of ITC focused on a metal–helium interface at low temperatures.¹ Nowadays, ITC is of interest in transport theory, technological applications concerned with condensation–evaporation processes, nanomaterials, and energy transfer in biomolecules.^3–7 The nanometer length scales characteristic of electronic components, nanomaterials, and their interfaces give rise to large heat energy densities and significant heat fluxes, with characteristic dissipation times across interfaces determined by the ITC. Therefore, an important objective is to understand the microscopic mechanism determining the magnitude and dependence of the ITC with intermolecular interactions and experimental conditions (e.g., temperature).

Traditional models for heat transport, such as diffuse mismatch and the acoustic mismatch model, rely on the knowledge of bulk thermal transport properties of materials and fluids and assume elastic phonon scattering at the interface, neglecting the local interfacial structure.² At the nanoscale, interfaces often feature local spatial heterogeneity spanning nanometer length scales.⁸ The theoretical modeling of interfacial heat transport is more complex, making it necessary to incorporate interfacial modes,⁹ interfacial structural details,¹⁰ and multi-phonon scattering events.¹¹ Interestingly, recent experiments of ZnO/GaN heteroepitaxial interfaces showed that the ITC could depend on vibrational modes characteristic of the materials.¹² This concept goes beyond the vibrational mismatch between two bulk materials assumed in the diffuse mismatch model.

The importance of wettability on the ITC was analyzed recently, showing that the conductance increases with the hydrophilicity (wettability) of surfaces in contact with water.¹³ Studies of flat¹⁴ and curved-nanometer size^15,16 surfaces provide a correlation between surface wettability, surface structure,¹⁷ interfacial width,¹⁸ and thermal conductance. Recent works¹⁹ demonstrated a correlation between interfacial fluid adsorption and thermal conductance of liquid–vapor interfaces. These studies provide further impetus to develop theoretical approaches to quantify interfacial heat transport by incorporating interfacial degrees of freedom explicitly.

In this work, we investigate the interfacial thermal conductance of the gold–water interface using Non-Equilibrium Molecular Dynamics (NEMD). This interface has been studied recently using polarizable and non-polarizable models.²⁰ The authors used the so-called density readjusting embedded atom model for the polarizable model (DR-EAM). They reported a slight increase (5%–10%) in the ITC when the calculations included the polarization. There is considerable interest in understanding the role of polarizability on transport processes, and different approaches have been considered in the literature to tackle this problem. In addition to the DR-EAM approach, which resolves the atomic valence density dynamically, the Drude–Lorentz method is widely used to investigate polarization effects in molecules²¹ and solids.²² 

In this work, we quantify the interfacial thermal conductance of polarizable and non-polarizable gold–water interfaces. To address the impact of surface wettability, which is known to influence the ITC, we have chosen two models that reproduce the experimental surface energy of gold, the gold–water monolayer hydration energy and the gold–water interfacial tension. We will show that the predicted interfacial thermal conductances obtained with the polarizable and non-polarizable models are significantly different. We explain this observation in terms of a coupling of internal degrees of freedom of the polarizable force fields and the density of states of liquid water, using for the latter both rigid and fully flexible water models. Our work provides the proof of principle on the lack of a direct correlation between wettability and thermal conductance due to the role of internal degrees of freedom of the gold surface on one side of the interface. We further explore approaches to resolve this problem and predict consistent ITCs by preserving polarization effects in interfacial heat transport computer simulations.

We performed simulations of the gold(111)–water interface and calculated the ITC by setting an explicit heat flux in the direction perpendicular to the interface plane. We employed the SPC/E water model²³ and a flexible water model based on the SPC model (mSPCFw).²⁴ For gold, we employed non-polarizable (NP)²⁵ and polarizable (P) force fields.²² The polarization calculations rely on the Lorentz–Drude oscillator model, which is widely used to include polarization effects in the simulation of condensed phases. The polarizability of the gold atoms was modeled by adding an electron “shell” to each gold “core,” with shell and core sites bonded through a harmonic bond with force constant, k_e. This model reproduces accurately the image potential induced by a charge next to a metallic surface.²² The P and NP models predict the density and surface free energy of gold, and the gold–water interfacial tension 1.23 J/m² (P) and 1.24 J/m² (NP), in good agreement with experiments (see Refs. 22 and 25). The surface free energies of Au(111) [1.55 mJ/m² (P) and 1.54 mJ/m² (NP)] are also in excellent agreement with experiments (1.54 mJ/m²²⁶). In the Appendix, we include additional details on the force fields employed and simulation protocols employed in this work.

We show in Fig. 1 the simulation set-up employed in this work, consisting of a gold slab in contact with water. The thermal gradient was set up by thermostatting the water molecules (cold) at the edges of the simulation box and the entirety of the gold slab (hot) (see the Appendix for details). All the simulations were performed at constant volume using a simulation cell with dimensions, 113.06 × 28.84 × 24.97 ų, with the volume adjusted to ensure that the density of water far from the gold surface was ρ_m = 1.0 g/cm³. Following pre-equilibration using a Nosé–Hoover thermostat (with damping time τ_D = 0.1 ps) at the average temperature of the hot and cold regions, we activated the thermostats. The non-equilibrium system reached the stationary state after 10⁵ steps (∼100 ps). We show in Fig. 1 representative temperature profiles for the P and NP models. Despite having the same interfacial free energy, the models predict very different temperature jumps, ΔT, and very different thermal conductances (see Fig. 2). The thermal conductances were calculated using equation G = J_q/ΔT, with the heat flux J_q obtained from the continuity equation (see the Appendix for additional details).

The temperature profile of water molecules can be obtained using the velocity of the center of mass and the rotational velocities. In most simulation codes, the trajectories contain atomic positions and velocities. Analyzing these atomic velocities to compute temperature profiles is advantageous at interfaces, where the water molecules might adopt preferred orientations. There are examples in the literature of oscillatory temperature profiles (see, e.g., Ref. 28). Hence, we calculate below the number of degrees of freedom needed to extract temperatures from the atomic velocities.

Since the SPC/E water model used in the present work is rigid, care must be taken in order to account for the reduced degrees of freedom for the water molecules. The total number of degrees of freedom for the water molecule is six, three translational and three vibrational. To compute the temperature for each atom type, these six degrees of freedom are usually distributed in the three atoms, resulting in two degrees of freedom per atom. A more accurate description can be achieved by considering the molecular geometry. In the oxygen frame of reference, the coordinates of the hydrogen atoms can be written as

where r₀ and θ₀ are the equilibrium bond length and angle. The center of mass of the molecule (in the oxygen frame of reference) is

where m_H and m_O are the masses of hydrogen and oxygen, respectively. The moments of inertia of the molecule are given by

Assuming equipartition of kinetic energy, the angular velocities of the water molecule can be written as

where ɛ = k_BT/2 is the energy per degree of freedom. The linear velocities of the molecule are

and the kinetic energy for the hydrogen and oxygen atoms

The number of degrees of freedom needed to calculate the temperature are given by $Ndof=2KkBT$⁠. Solving the system of equations for the SPC/E geometry (r₀ = 1.0 Å and θ₀ = 109.47°) gives the following effective degrees of freedom for the atoms in the water molecule:

Figure 2 shows the temperature profiles for the non-polarizable case by considering the number of degrees of freedom per atom given in Eqs. (8) and (9). The red symbols, corresponding to two degrees of freedom per atom, show significant deviations with respect to the temperature profiles using the effective degrees of freedom per atom, possibly due to the orientation of the water molecules in the first few interfacial layers.

The temperature jump was calculated by fitting the temperature profiles from each replica (for the gold and water molecules) to the heat diffusion equation

where J_q is the heat flux, κ is the thermal conductivity, and x₀ and T₀ represent a reference coordinate and temperature, respectively. The fittings were restricted to a narrow interfacial region of $∼7$ Å (three sampling points for each phase with bin width $a0/3$⁠; see Fig. 1) where κ can be assumed to be approximately constant.. The temperature jump was found by evaluating the fittings at the interface, which is located at $xi=±a0N3/6$⁠, where a₀ = 4.078 Å is the gold lattice parameter and N = 15 is the number of gold layers in the slab (see Fig. 1). This relationship follows from the distance between (111) planes, $a0/3$ and N Au planes of atoms centered at x = 0. The location of the interface at x_i corresponds to the plane between the outermost gold layer and the water layer in direct contact with gold. The heat flux, J_q, was calculated by averaging the energy exchange rate at the hot and cold thermostats, and at stationary conditions, $Q̇stat$⁠, as $Jq=Q̇stat/2A$⁠, where A is the cross-sectional area of the simulation box.

We note that previous studies extrapolated the temperature profiles from the bulk to the interface to computer temperature jumps.²⁸ In that work, the interfacial temperature was not well defined, and it showed an oscillatory behavior. However, our temperature profiles using the procedure detailed in Sec. II B are well defined and do not show oscillations. The calculation of the temperature jump using the interfacial profiles takes into account the changes in density and thermal transport of the solvation layers as well as the changes in the dynamics of the interfacial molecules with respect to the bulk liquid, which is manifested in a different vibrational density of the states (VDoS) (see Fig. 7 in Ref. 15 for a comparison of interfacial and bulk VDoS in a simple liquid).

The phonon dispersion relations were obtained using equilibrium molecular dynamics simulations and the method and auxiliary code phana described in Ref. 29. The dynamical matrix was constructed from the atomic displacements of one system containing 1728 gold atoms (12 × 12 × 12 supercell) during 1 ns, at a temperature of 300 K. For the P model, we used only the core (gold) atoms.

The interfacial VDOS was computed from the velocity autocorrelation function (VACF) of atoms initially belonging to the interfacial region (defined as the first peak in the density profiles for water, and the first atomic layer for gold), during a time window of 4 ps. After this time, the atoms considered for the VACF were re-defined. The expression for the DOS is given by

where $vα,ξ(0)vα,ξ(t)$ is the VACF of atom type α, with perpendicular (x) or parallel (y, z) component ξ to the interface. The intensity of the VDOS for the hydrogen atoms shown in the results section includes a factor of 2.

The P model predicts thermal conductances $∼2$ to 3 times higher than the non-polarizable model in a wide range of temperatures (see Fig. 3). The thermal conductances obtained with the NP model are similar to those reported in simulations of gold nanoparticles [150–170 MW/(K m²)].¹⁶ Similar values have also been reported in simulations of water boiling on flat gold surfaces, 120 MW/(K m²).³⁰ In a recent article, the thermal conductance of polarizable surfaces was computed using the density readjusting embedded atom method,²⁰ reporting 279 ± 16 MW/(K m²), close to the ITC reported also in that work in the absence of polarizability, 254 ± 17 MW/(K m²). These values are much lower than the ITC obtained with the polarizable model employed here, $∼540$ MW/(K m²). The inclusion of flexibility in the water model results in an even larger ITC, $∼875$ MW/(K m²). We explain below the physical origin of the different ITCs.

The results presented in Fig. 3 are surprising since the polarizable and non-polarizable models predict the same gold–water interfacial tension and gold surface energy. Additional simulations using the flexible water model and the non-polarizable surface show that the ITC does not depend on the internal degrees of freedom of the water molecule. Hence, to understand the origin of the differences between the P and NP models, we analyzed the interfacial structure of water, the phonon dispersion curves of gold, and the vibrational density of states of gold and water.

The density profiles of P and NP models do not show noticeable differences (see Fig. 3, bottom). Both profiles show significant layering, reflecting the strength of the water–gold interaction. The strong layering is consistent with experimental observations,³¹ indicating a significant deviation of the interfacial structure with respect to the bulk one. This highlights the importance of quantifying the temperature jump using the temperature gradients of the interfacial layers (see Fig. 1). Within the water layers in direct contact with gold, the hydrogens feature a stronger preference to approach the gold surface leading to a preferred orientation of the water molecules in that solvation layer (see Fig. 10 in the Appendix). Overall, these results do not reveal a correlation between the thermal conductance and the adsorption of water molecules on the gold surface since the latter is the same for polarizable and non-polarizable surfaces. The independence of the adsorption of water with the polarizability of the gold surface is consistent with recent works of ionic liquids confined between metallic surfaces.³² These studies showed that polarization does not modify the interfacial structure of the ionic liquid–metal interface significantly. Recent studies of ionic liquids nano-confined between polarizable surfaces did also show a minor impact of polarization on friction forces.³³ 

The ITC has been described before using the density of states of the two bulk phases in contact.² To rationalize our results, we analyzed the phonon dispersion curves of “bulk” gold using the method discussed in Ref. 29. Figure 4 shows the phonon dispersion relations for the polarizable and non-polarizable gold models. There is no significant difference between the phonon dispersion for the core particles of the polarizable model and the non-polarizable gold. Hence, the differences reported in Fig. 2 cannot be rationalized on the basis of differences in the phonon dispersion of the core atoms.

We extended our analysis of dynamic properties by computing the vibrational density of states (VDOS) of interfacial gold and water using the velocity autocorrelation functions. The low-frequency regions are fairly similar for both P and NP models (see Fig. 5). However, at high frequencies, the P model features a vibrational band, which emerges from the harmonically coupled core–shell charges. These high-frequency bands overlap with the hydrogen VDOS of water. We argue that this overlap is responsible for the enhanced thermal conductance reported in Fig. 3-top. To test this hypothesis, we performed additional simulations and shifted the high-frequency modes associated with the core–shell harmonic bond. We achieved a systematic blue-shift by reducing the mass of the shell-site, m_s (see Fig. 6-top and bottom). The timestep of these simulations was adjusted accordingly in order to generate stable trajectories (see  Appendix A for details). While the structural properties remain indistinguishable from the original model, we found that the ITC varied strongly with the frequency shift, featuring a non-monotonic dependence on the mass of the shell particle (see Fig. 7). The coupling between water and gold is particularly strong at m_s ∼ 4–6, and the ITC reaches a maximum at about 2250 MW/(m² K), thus maximizing the interfacial heat transfer. At m_s < 1, the ITC decreases, and it converges to the ITC calculated using the NP model. We find similar behavior in the simulations performed with the flexible water model. The ITC reaches significant values at m_s ∼ 4, and for small Drude masses m_s, the ITC converges to values much higher than those obtained with the rigid water model. This overestimation of the ITC with the flexible model indicates the existence of significant coupling between the gold and water internal vibrations. We note that while the mSPCFw model has been fitted to reproduce the vibration modes of water, the fact that is a classical model results in a too high occupancy of the high frequency modes, leading to the over-prediction of the ITC. This over-prediction highlights the inadequacy of the flexible model (in a classical simulation) to describe the thermal transport across the polarizable gold–water interface.

To rationalize the non-monotonicity of the ITC with m_s we computed the overlap of the hydrogen and high-frequency bands of the P-gold by evaluating the following “overlap” integral:

where D_Au,s and D_H represent the VDOS of the shell atoms in gold and the hydrogen atoms in the water molecules, respectively. Figure 7 (middle panel) and Fig. 8 (middle panel) show that the overlap integral features a maximum at approximately the same m_s values where ITC is maximized. Furthermore, the overlap decreases significantly for m_s < 4 highlighting the correlation between the ITC and the degree of overlap of the vibrational modes in gold and water. We present in Fig. 7 (bottom panel) and Fig. 8 (bottom panel) a more concrete correlation between the ITC and the overlap of vibrational modes. The ITC clearly increases with the degree of overlap. The turning point at $S/Sms=1.0$ is due the maximum in G and the overlapping integral. We note that in both branches, G increases with the overlap integral as expected.

We have shown that internal degrees of freedom in polarizable solids simulated with the Lorentz–Drude model lead to significant changes in the ITC of solid–liquid interfaces and an overestimation of the ITC over the values obtained with non-polarizable models. The different thermal conductances obtained with polarizable and non-polarizable models can be explained by considering the internal degrees of freedom of polarizable solid. The coupling of the core–shell internal degrees of freedom with the librational modes of the water molecules leads to a significant enhancement of the ITC. The ITC can be modified systematically by introducing a blue shift in the core–shell vibrational modes toward higher frequencies, reducing the overlap of core–shell and water vibrational modes. For minimal overlaps of the corresponding VDOS, the ITC converges toward the values predicted by non-polarizable models. Thus, the systematic modification of the mass of the shell particle in the Drude oscillator provides a route to eliminate the over-prediction of the ITC, preserving the polarization response of the surface. Since the latter depends on the polarizability of the gold atoms through α = q²/(2k_r), it is, therefore, independent of the mass m_s. However, we obtain convergence between the polarizable and non-polarizable ITCs at small masses, m_s = 0.2, which requires very short timesteps, δt ∼ 0.1 fs, therefore increasing the computational cost of the simulations.

In addition to the SPC/E rigid model, we have presented results for the mSPCFw flexible water model, which features additional vibrational bands owing to the bending and bond stretching modes. In this case, reducing m_s does not result in the convergence of the ITC to the values obtained with the P and NP models using a rigid water model. Instead, the ITC is overestimated for all masses investigated here, highlighting the limitations of “classical” flexible water models to compute the ITC when they are combined with polarizable surfaces modeled with the Lorentz–Drude model. Our results indicate that the classical Drude oscillator models, widely used in computational studies of materials,²¹ should be employed with care, as they might lead to enhanced heat transport across interfaces.

Finally, we note that beyond the numerical calculation of the ITC examined here, the comparison of the results obtained with polarizable and non-polarizable models indicates clearly that the surface free energy cannot be used as a single variable to establish correlations for the ITC, namely, the ITC is not directly correlated with the interfacial free energy. Indeed, the two models studied here predict essentially the same interfacial free energies, but the ITC is largely different. Our results call for careful consideration of the fitting procedures to simulate ITC, which should focus on interfacial free energies, wettability, and vibrational coupling effects specific to the models employed in the computations.

We acknowledge the Leverhulme Trust for the award of Grant No. RPG-2018-384, the Imperial College London RCS High Performance Computing facility, and the UK Materials and Molecular Modeling Hub for computational resources (partially funded by EPSRC through Grant Nos. EP/P020194/1 and EP/T022213/1).

The authors have no conflicts to disclose.

The data that support the findings of this study are available within the article.

The non-equilibrium simulations were performed using hot and cold thermostatting regions for gold and water, respectively. The gold slab had a lattice parameter a = 4.078 Å and consisted of 15, 10, and 10 atomic layers in the x, y, and z directions, respectively. The dimensions of the slab correspond to 35.3 × 28.84 × 24.97 ų. The gold slab was placed in contact with water at an average density ρ_m = 1.01 g/cm³, which results in a density of ρ_m = 1.0 g/cm³ in the bulk region. The thermostats are set as follows: all the atoms in the gold slab are thermostatted using the Langevin thermostat with a time constant of 2.5 ps (see Ref. 34 for further information on the impact of the thermostat on the ITC). For water, we defined a thermostatting region at a distance of $∼22$ Å from the interface. The temperature difference between the hot and cold thermostats was set to 50 K. This temperature difference generates a heat flux, with a magnitude determined by the thermal conductivities of gold, water, and interfacial thermal conductance. The total volume of the simulation box is ∼113.06 × 28.84 × 24.97 ų. Figure 9 shows the internal energy exchanged by the thermostats for non-polarizable (blue) and polarizable (red) models. The symmetry of the energy exchanged in the hot and cold thermostats indicates good conservation of energy in the system.

The simulations were carried out using a Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)³⁵ and the following procedure. First, we equilibrated the system to the target volume using a Nosé–Hoover thermostat at a temperature T = (T_cold + T_hot)/2 and a thermostat damping parameter of 0.1 ps. After 10 ps, the hot and cold thermostats are introduced in their respective regions using a canonical velocity rescaling³⁶ for the water molecules and a Langevin thermostat for the gold atoms, using a damping parameter of 2.5 ps. The damping parameter employed here for the Langevin thermostat does not modify appreciably the vibrational density of states of the solid and enables a good thermalization of both core and shell particles in the Drude oscillator (see Ref. 34). The temperature profile is allowed to reach the stationary state for 100 ps. Following the equilibration period, we computed the temperature profiles and heat flux over 400 ps. The timestep was set according to the mass of the shell particles: δt = 0.1 (m_s = 0.2–0.3 u), δt = 0.2 fs (m_s = 0.4 u), δt = 0.25 fs (m_s = 0.5–0.7 u), δt = 0.4 fs (m_s = 0.8–0.9 u), and δt = 0.5 fs (m_s ≥ 1.0 u). For the non-polarizable model, we set δt = 1.0 fs for the rigid and δt = 0.5 fs for the flexible water models, respectively. We used 10–20 independent replicas to compute the data reported.

The non-polarizable model is defined in terms of a Lennard-Jones 12-6 potential, parameterized by Heinz et al.²⁵ based on the density as well as the surface energy of the (111) surface. In conjunction with the SPC/E water model and geometric combination rules, the model provides a good estimate²⁵ of the gold–water interfacial tension (1.24 J/m² vs an experimental upper limit of 1.47 J/m²).³⁷ 

To include polarization effects, Geada et al.²² extended the Lennard-Jones model of Ref. 25 by adding a core–shell Drude oscillator pair for each gold atom. The core–shell pair interacts through a harmonic spring only. This model reproduces the classical image potential, lattice parameters, surface energy, and gold–water interfacial tension (1.24 J/m²). The polarizable model consists of core particles (with charge q = +1.0e) attached with a harmonic spring to dummy “electrons” with charge q = −1.0e. We report in Table I the bonded and non-bonded interactions.

Water was modeled using the rigid SPC/E model,²³ and the flexible water model (mSPCFw) parameterized by Smirnov²⁴ to reproduce the vibrational density of states of a water molecule in the gas phase. The flexible model, based on the SPCfw model,³⁸ uses a harmonic potential for the H–O–H angle, and a Morse potential for the O–H bond with parameters summarized in Table II,

Lennard-Jones cross interactions between different atoms were obtained using geometric combination rules for σ and ɛ.

The number density profiles of oxygen and hydrogen atoms for both the polarizable and non-polarizable gold models and the rigid SPC/E water model are shown in Fig. 10. We note that, for the polarizable gold model, there is an enhancement of the number density of hydrogen atoms at the gold surface. This is consistent with a reorientation of the water molecule with the hydrogen atoms pointing toward the gold, as discussed in the main text.