Thermal conductivity at solid/fluid interfaces: from adsorption to phonon scattering through the rattle effect

Abstract

Of relevance to adsorption, separation and catalysis, fluids in nanoporous materials display intriguing phenomena arising from their underlying structural, dynamical, and thermodynamic behavior. While many effects at the solid/fluid interface in nanoporous materials are well-documented, the thermal behavior of such solids subjected to fluid adsorption remains to be deciphered. Among exotic mechanisms, the so-called rattle effect, which corresponds to the decrease in the solid thermal conductivity induced by phonon scattering at the solid/fluid interface, has received only little attention. This phenomenon, which challenges existing mixing rules and effective medium approaches, has been identified in nanoporous materials (e.g. zeolite). Considering that this nanoscale effect is necessarily restricted to a small region near the solid/fluid interface, its impact for less finely divided materials (i.e. with thicker solid/fluid domains) remains to be established. Here, we address this question by employing a molecular simulation strategy to investigate thermal transport in a prototypical nanoporous silica material filled with a simple fluid. As a result, while the conventional behavior predicting an overall thermal conductivity increase upon fluid addition is qualitatively recovered in many cases, we also observe some situations where the overall thermal conductivity κ in such hybrid systems is mostly unchanged as conductivity through fluid is counterbalanced by the rattle effect at the solid surface. Understanding such intrinsic relations and the rattle effect that govern the thermal conductivity in solid/fluid systems paves the way for the design of novel materials to harness thermal processes in practical applications.

---

1 Introduction

With pressing challenges in environmental engineering (e.g. agriculture, water treatment) and energy applications (e.g. storage, conversion), the development of novel separation and catalysis processes involving nanoporous materials is urgently needed.^1,2 In this context, the large spectrum of available nanoporous solids is expected to be a game-changer by contributing for instance to the “seven key chemical separations to change the world”.³ In such processes, the large surface area and strong confinement within the porosity of the host nanoporous material are used to harness the broad range of phenomena occurring when the solid is set in contact with a fluid.^4–7 From a more fundamental viewpoint, while the transport of fluids in nanoporous materials has been reasonably clarified in the last decades, many microscopic mechanisms that rule mass and energy transfers at the nanoscale in such systems remain to be explored. In particular, most of these transport processes are often assumed to obey effective medium theories that do not properly account for the microscopic behavior of the nanoconfined fluid and nanoporous material.^8–10 Specifically, with these strategies, the properties of the nanoconfined fluid are assumed to remain identical to their bulk counterpart despite well-documented reports on interfacial effects that can affect thermal and mass transport mechanisms (phonon scattering,¹¹ heat diffusion restriction,¹² pressure-like effect,¹³ etc.).

Among the challenges associated with fluid confinement in nanoporous materials, modeling thermal transport introduces complex features such as the rattle effect – a reduction in the system's thermal conductivity caused by the presence of fluid atoms at the solid surface. This phenomenon has been observed in a wide range of materials such as skutterudites,¹⁴ carbon nanotubes,¹⁵ methane hydrates,¹⁶ metal organic frameworks (MOFs)¹⁷ and zeolites.¹¹ Formally, the rattle effect, which can amount in some cases to a 40% decrease in the solid's thermal conductivity,¹¹ arises from phonon scattering at the solid/fluid interface as fluid atoms collide with the atomic structure of the host nanoporous material.^11,17 The effect is also observed in a experimental study with HKUST-1, which further show that liquid adsorbates can reduce thermal conductivity by up to 80%, an effect attributed to enhanced phonon scattering and not captured by effective medium approximations.¹⁸ This phenomenon is necessarily omitted from mixing rules and, in general, from effective medium theories as both the fluid and solid contributions are usually assumed to be equivalent to their bulk/fluid-free counterpart (i.e. with properties as assessed in the absence of the other phase).

The above question is even more puzzling as the available literature on the topic provides somewhat contrasting pictures regarding the rattle effect.^19,20 While phonon scattering consistently appears to reduce thermal transport at the molecular level, the resulting transport properties of the fluid + solid system depend on specific features of the solid phase – such as its surface area and wall thickness – but also the solid/fluid interaction strength. This complexity makes predictions of the overall thermal conductivity through effective medium theories particularly challenging. In this context, unraveling the thermal coupling between the host solid and the confined fluid constitutes a key challenge to design novel materials and processes with tailored thermal management properties. In particular, the occurrence of rattle effects opens perspectives to modulate thermal transport at the nanoscale by considering fluids adsorbed in nanoporous materials.

Here, to tackle some of the questions reported above, we investigate thermal transport in a nanoporous system using a simple yet realistic model consisting of a silica nanopore filled with a fluid. Overall, our molecular simulations show that the rattle effect remains present for all solid/fluid interaction strengths ε due to phonon scattering at the solid/fluid interface. Depending on ε and the relative contributions of the fluid and solid phases, the overall thermal conductivity κ in such hybrid (solid + fluid) systems can be increased or unchanged with respect to that obtained for the empty solid. While the first situation (thermal conductivity increase) corresponds to the conventional result where adding a fluid into the solid porosity enhances the overall thermal conductivity, the other situation cannot be predicted using classical mixing rules. By disentangling the thermal conductivity contributions from the solid phase, the fluid phase and the solid/fluid interface, we identify the role played by the solid/fluid interaction strength ε and the adsorbed amount per unit of surface area n_a. We also provide a simple model which allows rationalizing the decrease in the solid's thermal conductivity by considering the average phonon lifetime resulting from intrinsic phonon scattering within the solid and from solid/fluid collisions at their interface.

2 Theoretical background and formalism

2.1 Thermal conductivity and microscopic energy flux

Through the prism of statistical mechanics, the thermal conductivity κ can be evaluated using the Green–Kubo formalism by computing the time autocorrelation function of the total heat current j(t) within molecular dynamics simulations. Due to the anisotropic nature of the slit pore geometry, the Green–Kubo formulation is applied to compute the thermal conductivity κ^‖ in the parallel direction to the surface using the following formula:where V is the volume of the system, k_B the Boltzmann constant and T the temperature. The autocorrelation function 〈j^∥(t)·j^∥(0)〉 is obtained by considering the heat current j(t) in the x, y plane, which is parallel to the surface.²¹ The heat current j(t) at a given time t within the system is computed as the time derivative of the sum of each atom position r_i multiplied by the (potential and kinetic) energy e_i carried by the atom . By explicitly deriving the heat current with respect to time, we arrive at the following relation:

The first term, which corresponds to the product of the atom's energy by its velocity, is defined as the convective contribution to the total current. The second term is referred to as the virial contribution as it arises from the variation in time of the total energy of the atom (it is therefore related to the transfer of potential energy between atoms in the system).

2.2 Heat current definition

Care must be taken with multi-component systems to ensure consistency with macroscopic transport theory. Formally, the derivation of the entropy production rate leads to a coupling between the different fluxes and their underlying driving forces. The driving force for mass flux appears as the gradient of the chemical potential divided by temperature, ∇(μ/T).²² To recover the thermal conductivity as the coefficient relating the heat current j solely to the temperature gradient ∇T, it is necessary to remove the temperature dependence from the diffusive contribution. This is achieved by redefining the heat current j to exclude the energy carried by mass transport: where δ corresponds to the different species within the system. In this expression, j′ is the total energy current given by eqn (2), h^δ the partial molar enthalpy of species δ, j^δ_M the mass current of component δ, and M^δ the molar mass of component δ.^22,23 The latter expression for the microscopic heat current can be expressed as:where i ∈ δ indexes atom of species δ, v_δ is the center of mass velocity of δ species and σ_i is the atomic stress of the i^th atom which will be defined below.

We recognize in eqn (3) the convective and virial contributions to the heat current (j^con and j^vir) for a multi-component system. The above expression can be recast as a per-atom representation. The convective term can be given as the product of this relative velocity and the total energy of each atom while the virial term can be written as a sum of atomic stress σ_i distribution based on their relative positions to the geometric center of the N_ϕ interacting atoms within a given potential energy function (which may be a bond, angle, dihedral, etc.):²⁴

j^con_i,δ = e_i(v_i − v_δ), (4a)

where ϕ denotes a potential energy interaction involving N-atoms (e.g. non-bonded Lennard–Jones interactions, molecular bonds, angles, dihedrals, etc.), r^ϕ_i0 = |r_i − r^ϕ₀| is the relative position of atom i to the geometric center r^ϕ₀ of a potential energy term ϕ, F^ϕ_i the force on atom i due to the potential energy term ϕ. With these notations, the previously introduced atomic stress tensor is defined as .

3 Computational methods

3.1 Nanoporous silica and fluid models

The slit nanopore considered in this work is built from a realistic model of amorphous silica material. The procedure employed to prepare the silica walls was described in detail elsewhere.^25–27 In short, a block of crystalline silica (cristobalite) is initially melted at T = 4000 K using a molecular dynamics simulation run. The obtained system is then rapidly quenched to T = 300 K to obtain an amorphous silica block. The resulting simulation box, initially elongated along the z direction, was then cut to generate two free surfaces. The slab thickness was intentionally chosen to be L_z = 35.8 Å to match the in-plane dimensions L_x = L_y = 35.8 Å. Periodic boundary conditions were applied in the x and y directions. To ensure that the surface is realistic from a chemical viewpoint, all uncoordinated Si atoms and any O atoms not bonded to at least one Si atom are first removed and the remaining dangling bonds of the O atoms at the silica surface are saturated by adding H atoms to form silanol groups. After surface saturation, the final silica slab contained 2601 solid atoms, comprising 779 silicon, 1646 oxygen, and 176 hydrogen atoms.

The number of OH bonds is adjusted to obtain a final silanol density of 6.9 OH per nm² (this surface density is estimated by taking the geometrical surface area of the silica surface) – a value which is consistent with typical experimental conditions.²⁸ The dimensions of the surface are L_x = 35.8 Å × L_y = 35.8 Å. The directions parallel to this surface are treated as infinite by applying periodic boundary conditions in the x and y directions. As the realistic silica surface used in this work does not exhibit a strictly identical molecular structure in the directions parallel to the surface, we looked at possible anisotropy along the x and y directions. In practice, no significant differences in thermal properties were observed when considering these two directions. Therefore, the thermal conductivity was only evaluated in the directions parallel to the silica surface. In this direction, finite-size effects associated with phonon/phonon interactions are negligible. This is supported by the fact that both the thermal conductivity and vibrational density of states for bulk amorphous silica and the silica slab from which it is extracted match (see Results section discussion).

Silica is chosen here as a prototypical surface as a large body of available materials with large specific surface area are made up of oxides. The silica slab thickness of e = 35.8 Å is chosen large enough to ensure that fluid atoms on opposite sides of the silica surface do not interact. The slit nanopore is cut so that the pore width (which is defined between the two opposite surfaces) is about H = 10 nm. The interactions between atoms of the amorphous silica were described using the force field by Carré et al., which was parameterized through a fitting scheme based on Car-Parrinello molecular dynamics simulations.²⁹ The united atom model TraPPE was used to describe the molecular interaction between the fluid atoms (with parameters chosen to describe the whole methane molecule).³⁰ A Lennard–Jones (LJ) potential was also used to describe the interactions between methane molecules and the Si and O atoms of the silica surface. The LJ parameters for Si and O are derived from the CLAYFF force field³¹ and combined using the Lorentz–Berthelot combining rules.^32,33 All parameters and details are provided in the SI.

3.2 Grand canonical Monte Carlo and molecular dynamics simulations

In Grand Canonical Monte Carlo (GCMC) simulations, the system has a constant volume V and is set in equilibrium with an infinite reservoir of methane molecules imposing their chemical potential μ and temperature T.^34,35 Here, the adsorption isotherm is simulated at T = 300 K by varying the chemical potential μ of the reservoir. All simulations are performed with LAMMPS which combines Monte Carlo (MC) and Molecular Dynamics (MD) steps, therefore sampling the phase space in the grand canonical ensemble.³⁶ This hybrid scheme ensures that the system evolves through classical MD at fixed temperature, while Monte Carlo insertion and deletion moves of fluid molecules are attempted at regular intervals to enforce the matching with the prescribed chemical potential. In practice, the hybrid GCMC–MD simulations were performed with an initial equilibration of 1 × 10⁵ steps, followed by 5 × 10⁵ production steps, during which insertion and deletion attempts were carried out at every 10 MD steps.

The absolute adsorbed amount N_a(P) is estimated as the ensemble average of the number of adsorbed molecules as a function of the gas reservoir pressure P. To establish the correspondence between pressure and chemical potential, a set of reference pressures is first selected and used to define the chemical potential of an ideal gas at the target temperature. These chemical potentials are then employed as inputs for two independent GCMC simulations: the first performed within the silica nanopore system and the second carried out for a corresponding bulk methane system at the same temperature and fixed volume. After equilibration, both simulations result in an averaged number of fluid molecules. The resulting average number of particles 〈N〉 in the bulk system is then used in a NVT simulation, for which the resulting real pressures were then evaluated by assessing the average system pressure.

Starting with the configurations obtained from GCMC simulations, equilibrium molecular dynamics (EMD) simulations are conducted with LAMMPS.³⁶ An in-house C⁺⁺ code coupled to LAMMPS is used to compute the heat flux vector within the system. These molecular simulations are performed with a cutoff of 14.0 Å and using the Particle–Particle Particle–Mesh (PPPM) solver for the long-range Coulomb and van der Waals forces. All the simulations were carried with an integration time step of 0.5 fs. The molecular simulations begin with a 0.5 ns equilibration run, during which the Nosé–Hoover thermostat is applied to the silica framework (while the fluid atoms are treated in the NVE ensemble). The proper thermalization of the fluid atoms through collisions with the silica surface is verified by assessing their average kinetic energy. Once thermal equilibrium is achieved, a production run is carried out in the NVE ensemble for an additional 0.5 ns (i.e., 10⁶ timesteps). Data for the energy flux calculations are collected at evenly spaced intervals of 1 fs during the production run. To accumulate enough statistically relevant data, 200 independent simulations are performed for each configuration considered in the present study (here, configuration refers to a pore loading with a given fluid adsorbed amount N_a). These molecular simulations differ in their initialization as each simulation run starts with a unique random velocity distribution for the fluid atoms at the start of the equilibration run.

4 Results

4.1 Fluid adsorption in nanoporous silica

With insulating materials such as amorphous silica, thermal transfer occurs through atomic lattice vibrations (phonons) since the electronic contribution to heat conductivity is negligible. To assess the impact of fluid adsorption on energy transfer in nanoporous materials, we employ here a realistic model of nanoporous silica. As discussed in the Methods section, the pore surface is built from a numerical model of amorphous silica as obtained using a cook-and-quench procedure. Fig. 1(a) provides a magnified view of the final silica structure obtained using this procedure. Fig. 1(b) also shows a SiO₄ tetrahedron taken from the core of the silica slab and a SiO₃ (OH) silanol group located at its external surface. This silica slab is used to prepare a simulation box corresponding to a nanoporous slit having a pore size H = 10 nm. Fig. 1(c) illustrates the silica nanopore filled with methane at a pressure P = 200 atm and a temperature T = 300 K. Each methane molecule is represented as a spherical Lennard–Jones particle with parameters taken from the TraPPE force field.³⁰ While this model adopts a simplified description of the fluid by neglecting internal degrees of freedom, previous molecular simulations have shown excellent agreement between the predicted thermal conductivity and experimental measurements (with only small deviations observed at low densities).¹¹

Fig. 1 Fluid adsorption and confinement in nanoporous silica. (a) Molecular configuration of the amorphous silica surface. The silicon, oxygen, and hydrogen atoms are shown as yellow, red, and white spheres, respectively. (b) Magnified views of the SiO₄ tetrahedra forming the amorphous silica material and the SiO₃ (OH) silanol groups at its surface. (c) Fluid atoms (cyan spheres) confined within the silica nanopore of a pore size H = 10 nm at a pressure of approximately P = 200 atm and a temperature T = 300 K. This silica model and pore geometry were used for all production simulations.

The number of fluid atoms confined within the silica nanopore is obtained using Monte Carlo simulations in the Grand Canonical ensemble (GCMC) as discussed in the Methods section. This approach ensures that the chemical potential and temperature of the atoms inside the nanoporosity match those of a bulk phase under specific conditions. Consequently, considering a pure fluid here, the gas pressure is also uniquely defined for a given set of temperature and chemical potential. The adsorption isotherm N_a(P) for methane confined at T = 300 K within the silica nanopore is presented in Fig. 2(a). The absolute adsorbed amount N_a is determined as the ensemble average of the number of adsorbed atoms as a function of the gas reservoir pressure P. This adsorption isotherm exhibits a Type I (Langmuir) behavior – with a steep initial uptake at low pressures followed by slower and slower increase until adsorption saturation. Despite the large pore size H, no capillary condensation occurs as methane is above its bulk critical temperature T_c = 190.56 K. For three pressures taken along the adsorption isotherm, we show in Fig. 2(b) the fluid density profile ρ_a(z) along the direction perpendicular to the pore surfaces. These density profiles exhibit a pronounced peak near the pore surface, which is indicative of the adsorbed layer formed due to the interactions between the fluid atoms and the solid atoms. Beyond this adsorbed layer, the fluid density becomes more uniform to recover the bulk fluid density in the pore center.

Fig. 2 Fluid adsorption, density and structure. (a) Methane adsorption isotherm N_a(P) at room temperature T within the silica nanopore as obtained using Grand Canonical Monte Carlo simulations. Note the use of a log scale for the pressure axis to better visualize the overall adsorption behavior. The color points represent some specific adsorbed amounts for which we show the density profiles in panel (b). The data shown in this adsorption isotherm are obtained for the nominal solid/fluid interaction strength ε₀. The purple dashed line corresponds to a Langmuir isotherm adjusted to the adsorption data, given by N_a(P) = N_satbP/(1 + bP), where N_sat is the saturation adsorbed amount and b is the Langmuir affinity constant. The inset shows the excess adsorbed amount N^ex_a as a function of the absolute adsorbed amount N_a for the three solid/fluid interaction strengths ε: 0.5ε₀, ε₀, and 2ε₀. Details about the calculation of N^ex_a(P) are given in the text. Error bars are smaller than the marker size. (b) Typical density profiles ρ_a(z) along the direction z perpendicular to the pore surface at different pressures P. The color code corresponds to the data points that are shown in the adsorption isotherm in panel (a). For the adsorbed amount N_a corresponding to P ∼ 200 atm, the inset shows the density profiles ρ_a(z) obtained for the three solid/fluid interaction strengths: 0.5ε₀, ε₀, and 2ε₀.

To explore the role of the solid/fluid affinity on adsorption and thermal transport, we varied the solid/fluid interaction strength ε. The value of the nominal interaction strength ε₀ is given in the SI. In more detail, while maintaining the absorbed density N_a at a constant value, the solid/fluid interaction strength was set to ε = ε₀, 0.5ε₀, and 2ε₀ where ε₀ denotes the nominal interaction strength (see Methods section). As this change affects the fluid energy landscape near the solid surface, the density profile ρ_a(z) is found to evolve upon changing ε. The chemical potential μ of the confined fluid evolves here as we change the solid/fluid interaction strength ε while keeping the total number of atoms constant (this explains why the fluid density in the pore center changes upon changing ε). The inset in Fig. 2(b) depicts the density profiles corresponding to each interaction strength ε as determined for the total adsorbed amount N_a = 6.6 molec./nm³. The stronger solid/fluid interactions lead to a stronger adsorption peak as molecules from the pore central region migrate towards the surface. To quantify this effect, we introduce the excess density N^ex_a which measures the excess of fluid molecules with respect to the bulk density under the same thermodynamic conditions ρ_a(0) (taken as the density at the position z = 0 which corresponds to the pore center):

where H corresponds to the pore width. The inset in Fig. 2(a) shows the excess adsorbed amount N^ex_a as a function of the adsorbed amount N_a for the three different solid/fluid interaction strengths ε. The behavior of N^ex_a as a function of N_a varies significantly depending on the solid/fluid interaction strength. For ε₀ and 2ε₀, the excess density N^ex_a initially increases with N_a to reach a maximum around N_a ∼ 3 molec./nm³. Upon further increasing N_a, the excess density N^ex_a decreases as the bulk density becomes larger and the adsorbed density near the surface saturates. In contrast, for ε = 0.5ε₀, the excess density is close to zero or even negative for all N_a as the adsorbed density near the surface remains small.

4.2 Vibrational spectra and solid/fluid dynamical coupling

Before discussing the overall fluid + solid thermal conductivity, we first discuss the impact of fluid adsorption on the vibrational properties of the solid and fluid.^11,15,37 Throughout the remainder of this paper, α and β refer to the fluid and solid phases, respectively. As shown in Fig. 3(a), we computed the vibrational density of states (VDOS), i.e. D(ω)^β, of the silica atoms within 3 Å layers taken at different positions with respect to the silica surface. We consider the following regions: surface [S], subsurface [S*], and bulk [B]. D(ω)^β is obtained for each region by taking the Fourier transform of the velocity autocorrelation functions where the sum runs over all atoms i in the corresponding layer. The VDOS for each region are shown in Fig. 3(b) where we also compare the results for the empty and filled nanopores (the adsorbed fluid is taken at a pressure of 200 atm, i.e. N_a = 6.6 molec./nm³). The VDOS for the silica slab D^β(ω) exhibits some features that reflect its chemical and amorphous nature. Notably, the peak observed in the high-energy range (ħω ∼ 450 meV) can be assigned to the stretching of the hydroxyl groups at the surface. This is consistent with experimental measurements reporting silanol stretching frequencies around ∼3600 cm⁻¹ (equivalent to 446 meV).³⁸ Additionally, D^β(ω) displays a broad distribution of vibrational frequencies, which arise from the disordered atomic structure of amorphous silica (this contrasts with crystalline materials where phonon modes appear as sharp peaks in the VDOS due to the strong positional and orientation ordering). As expected, the VDOS for the fluid phase D^α(ω) does not display an assembly of discrete narrow peaks, but rather a continuous distribution of states in the low energy range corresponding to ħω ∼ 0–25 meV.

Fig. 3 Surface and volume vibrational density of states. (a) Molecular configuration showing the surface [S], subsurface [S*] and bulk [B] regions considered within the amorphous silica material. (b) Vibrational density of states (VDOS) for the confined fluid D^α(ω) at a pressure P = 200 atm and for the amorphous silica material D^β(ω). The latter are computed from the atom velocity autocorrelation functions for the different systems, as described in the main text. Solid lines correspond to the VDOS curves smoothed using a moving-average window of 25 points. The corresponding raw spectra are shown in the same colors with higher opacity. The black line indicates the density of states D^β(ω) for bulk silica, while the blue and red lines indicate the density of states D^β(ω) for the silica wall in the presence and absence of the confined fluid, respectively. The purple dashed line indicates the density of states D^α(ω) for the confined fluid. Each panel represents D^β(ω) calculated for a 3 Å thick layer as depicted in panel (a): (from top to bottom) Surface (S), subsurface (S*), and bulk central part (B). The densities of states D(ω) for the bulk silica and confined fluid are identical across all panels. All simulations are conducted at T = 300 K. The data D(ω) are shown the linear scale and are avspatialover the three spacial directions. Adsorption modifies the surface VDOS, consistent with interfacial phonon scattering.

In agreement with previous results on fluid within zeolite cages,¹¹ the vibrational modes of the fluid atoms at the pore surface shows a weak coupling with those of the solid atoms. However, this coupling is difficult to observe here as it only affects atoms directly located at the solid/fluid interface. Indeed, considering the large pore width H = 10 nm, the fraction of fluid atoms at the solid surface remains overall very small. In other words, even if the adsorbed density N_a is large, a great number of molecules are distributed in the large porous domain available to the fluid – therefore, making the contribution of adsorbed fluid atoms small in the overall VDOS D^α(ω). As a result, D^α(ω)for the confined fluid remains very close to that of the bulk phase. In contrast, the presence of fluid atoms in contact with the solid surface impacts the vibrational properties of the solid surface layers [S] and [S*]. This effect can be seen in the low energy range where the fluid density of states D^α(ω) is non-zero. The adsorption-induced changes in D^β(ω) indicate that the solid atoms at the solid/fluid interface display a different vibrational behavior in the presence of adsorbed fluid atoms. As will be discussed in detail below, this result provides evidence for phonon scattering by the fluid at the surface – the so-called rattle effect – and the resulting decrease in the solid's contribution to the thermal conductivity. In contrast, as expected, the fluid phase does not impact the VDOS for the solid atoms in the bulk region of the silica slab. The VDOS for bulk amorphous silica, which was used to prepare the silica nanopore, closely matches the VDOS obtained for the bulk region [B] in the silica wall. This further corroborates that the rattle mechanism stems from the coupling between the fluid and the solid at their interface.

4.3 Thermal conductivity of fluid-filled nanoporous silica

The domain sizes in the system under study lead to fluid and solid phases with quite distinct contributions to the thermal conductivities compared to their bulk counterpart value. It is therefore important to assess how each phase impacts the effective thermal transport within the fluid + solid system. As discussed in Section 2, the heat current j can be computed as a sum of individual atomic terms. For an individual atom, this involves summing over its convective j^con and virial j^vir terms according to eqn (4a) and (4b), respectively. On the one hand, the convective term j^con accounts for the heat current contribution corresponding to the kinetic energy of atoms. On the other hand, the virial term j^vir stands for the potential energy exchange through the multiple interactions experienced by an atom as described by the atomic stress tensor σ_i defined in Section 2 [eqn (4b)]. For a given fluid atom i, this atomic stress tensor σ_i is expressed as a sum of interactions with atoms in the same phase (αα) and with the atoms of the other phase (αβ). This yields a total stress which can be written as σ_i = σ_i,(αα) + σ_i,(αβ) and a virial contribution which can be decomposed as j^vir_i = j^vir_i,(αα) + j^vir_i,(αβ).

Fig. 4(a) illustrates the total heat current autocorrelation function 〈j^∥(t)· j^∥(0)〉 at T = 300 K and P = 200 atm along the direction parallel to the surface. The effective thermal conductivity κ is obtained by integrating the heat flux autocorrelation function in eqn (1) as shown in Fig. 4(b). Fig. 4(c) shows the effective thermal conductivity κ^‖ in the direction parallel to the interface, with contributions from both the solid and fluid phases, as a function of the fluid density N_a within the nanopore. Changing the fluid/solid interaction strength ε impacts both the rattle effect at the solid/fluid interface but also the adsorption layer density. As a result, when considering the effect ε, the different trends in Fig. 4(c) illustrate the competition between the changes in the thermal conductivity through the solid phase and the changes in the thermal conductivity through the adsorbed layer. This point will be further discussed in the next section. The thermal conductivity of the empty silica nanopore is – a value slightly larger but quite consistent with the experimental thermal conductivity for amorphous silica (∼1.30 W m⁻¹ K⁻¹).³⁹ Considering the empty solid for which the fluid conductivity is necessarily zero (since N_a = 0), it is interesting to note that one recovers the simple mixing rule which predicts that the effective thermal conductivity is given by where ϕ_α = 0.642 and ϕ_β = 0.358 are the volume fraction occupied by the fluid and solid phases while and are the heat conductivity of the fluid and solid phases. In more detail, this yields κ^‖ = 0.565 W m⁻¹ K⁻¹ which corresponds to the horizontal red dashed line in Fig. 4(c). As expected, regardless of the interaction strength ε, κ^‖ tends to this value as N_a → 0. In contrast, as soon as the adsorbed amount n_a becomes non-negligible, the interaction strength ε controls the affinity of the N_a fluid atoms with the solid surface and defines the density of the adsorbed layer [Fig. 2(b)]. Although the overall thermal conductivity may remain approximately unchanged for weak solid–fluid interactions, the solid contribution κ_β systematically decreases due to the rattle effect (as will be shown later in detail). However, while κ_β decreases with n_a, we note that the intrinsic contribution κ_ββ to κ_β remains unchanged for all adsorbed amounts and fluid-interaction strengths. These data, which will be presented in detail in Fig. 5, indicate that the total conductivity κ remains unaffected by n_a for the lowest interaction strength as a result of the compensation between the rattle effect and the increase due to the addition of a fluid contribution. Therefore, changes in the interaction strength ε lead to different solid-fluid interfacial behavior and its contribution to the effective thermal conductivity. In the next section, this contribution will be analyzed to explain the different behaviors observed in the heat transport.

Fig. 4 Thermal conductivity κ of a fluid-filled silica nanopore. (a) Normalized time autocorrelation function of the heat current 〈j^∥(t)·j^∥(0)〉 at T = 300 K and N_a ≈ 6.6 molec./nm³ along the direction parallel to the pore surface. (b) Thermal conductivity κ^‖ as computed using the Green–Kubo formalism for P = 200 atm, T = 300 K and N_a ≈ 6.6 molec./nm³. The variable t represents the upper time boundary in the Green–Kubo integral given in eqn (1), i.e. the cumulative integral of the heat-current autocorrelation function 〈j^∥(t)·j^∥(0)〉 up to time t. (c) Thermal conductivity parallel to the pore surface κ^‖ as a function of the adsorbed amount N_a. The data are presented for the three solid/fluid interaction strengths ε: 0.5ε₀, ε₀, and 2ε₀ with the legend shown within the panel. The horizontal red dashed line represents the thermal conductivity of the empty nanoporous silica.

Fig. 5 Fluid, solid and interface contributions to thermal transport. (a)–(c) present the different contributions to the thermal conductivity parallel to the pore surface κ^‖ as a function of the number of molecules per unit of surface area n_a: solid contribution κ^‖_β in red, fluid contribution κ^‖_α in cyan, and solid/fluid interface contribution κ^‖_αβ in purple. We also show molecular configurations that illustrate the energy flux along the silica material, the energy flux through the confined fluid and the energy flux at the solid/fluid interface, respectively. (d) Reciprocal characteristic time τ_a for phonon scattering due to solid/fluid collisions normalized to the intrinsic scattering time τ₀ in the solid phase as a function of the surface density of adsorbed molecules within the pore n_a (see text). The linear dashed lines represent the linear trend for for low values of n_a. The inset shows the function f(n_a) which represents the dependence of the interfacial scattering rate 1/τ_a on the adsorbed density n_a. Across all panels, the triangles, circles, and squares represent the data obtained for the three solid/fluid interaction strengths ε, respectively: 0.5ε₀, ε₀, and 2ε₀.

Fig. S1 of the SI compares the effective thermal conductivity κ of the system as obtained in our molecular simulation and that predicted by using a classical mixing rule. While both datasets follow the same trend as a function of N_a, the thermal conductivity by MD is slightly above the one predicted by the classical mixing rule. This may be explained by the thermal conductivity increase due to virial contribution on adsorbed fluid molecules, which is not take under consideration in the classical model. Moreover, we note that the breakdown of classical mixing rules can already be seen from the results shown in Fig. 4(c). Indeed, the mixing rule approach predicts a unique dependence of the effective thermal conductivity on the fluid adsorbed amount – in stark contrast with our molecular simulation results. In other words, our results reveal that changes in the solid/fluid interaction strength leads to qualitatively different dependences of κ(n_a) on (n_a) which cannot be captured by conventional mixing rules.

5 Discussion

5.1 Solid phase, fluid phase and solid/fluid interface contributions

To shed light on the impact of adsorption on thermal transport in the fluid-filled nanoporous silica, the different contributions to the overall thermal conductivity were calculated. In what follows, we decompose the total thermal conductivity κ into distinct contributions to the heat current autocorrelation function 〈j^∥(t)· j^∥(0)〉. As discussed in Section 2, the instantaneous heat current j(t) can be expressed as the sum of contributions from different groups of atoms. Accordingly, j(t) = j_α + j_β + j_αβ where j_α, j_β, and j_αβ are the fluid, solid, and solid/fluid interface contributions, respectively. In more detail, the fluid contribution j_α to the heat current includes (i) the convective term j^con_α from all fluid atoms and the virial term j^vir_(αα) from fluid/fluid interactions. Similarly, the solid phase contribution j_β includes (i) the convective term j^con_β from all solid atoms and (ii) the virial term j^vir_(ββ) from solid/solid interactions. Finally, the solid/fluid interfacial contribution j_αβ to the heat current is equivalent to the virial term j^vir_(αβ) arising from solid/fluid interactions.

Using the heat current decomposition into three contributions, the overall thermal conductivity can be decomposed into 9 = 3 × 3 terms. Therefore, we define the following component contributions by the solid, fluid, and solid/fluid interface:

κ^‖_β= κ_β,α + κ_β,αβ + κ_β,β, (6a)

κ^‖_α= κ_α,α + κ_α,αβ + κ_α,β, (6b)

κ^‖_αβ= κ_αβ,α + κ_αβ,αβ + κ_αβ,β. (6c)

These 9 terms narrow down to 6 terms owing to the symmetrical nature of the heat flux correlations 〈j_δ(0)· j_δ′(t)〉 (with δ, δ′ = α, β, αβ): . All six independent thermal conductivity terms contributing to eqn (6) are reported in the SI. With this convention, the total in-plane thermal conductivity can be written as

κ^‖ = κ^‖_α+ κ^‖_β+ κ^‖_αβ. (7)

Fig. 5 shows the three contributions κ^‖_α, κ^‖_β and κ^‖_αβ as a function of the adsorbed amount per unit of surface area n_a. We consider n_a rather than the total adsorbed amount N_a as the latter does not capture how the fluid atoms distribute between the adsorbed layer at the pore surface and the bulk-like density in the pore center. The adsorbed amount is defined as , where z₁ ∼ 0 and z₂ define the boundaries of the adsorbed peak in the vicinity of the silica surface. This definition of the adsorbed amount at given thermodynamic conditions allows rationalizing the variations in the overall thermal conductivity κ when changing the solid/fluid interaction strength ε₀. All thermal conductivity contributions are evaluated for the different solid/fluid interaction strengths (see the caption of Fig. 5 for the labeling of each dataset). The schematic views in Fig. 5 illustrate the different contributions to the thermal conductivity as described in what follows. The top molecular configuration illustrates the phonons within the solid, which are responsible for thermal transport across amorphous silica. The middle molecular configuration illustrates the energy transfer between interacting fluid atoms within the confined phase. Finally, the bottom molecular configuration illustrates the energy transfer across the interface through the rattle effect between solid and fluid atoms.

5.1.1 Solid contribution κ^‖_β. As shown in Fig. 5(a), the solid's thermal conductivity κ^‖_β= κ_β,α + κ_β,αβ + κ_β,β decreases as the adsorbed amount n_a increases – therefore, indicating that the solid contribution to thermal transport is governed by the rattle effect. The gray points in the same plot correspond to the solid/solid interaction contribution κ_β,β, which remains constant and independent of both n_a and ε. This indicates that κ_β,β ∼ κ_β(0) represents the bulk solid thermal conductivity (shown as the dashed line in the Fig. 5). Therefore, the source of the rattle effect arises from the terms κ_β,α and κ_β,αβ, which assume negative values. While the solid/fluid interaction strength ε does not lead to drastic changes in κ^‖_β, a small yet noticeable influence is observed (this dependence may also be affected by how n_a is defined). This leads to larger fluctuations on κ_β when evaluated from the superposition of three datasets corresponding to different solid–fluid interaction strengths ε. Nevertheless, when considered individually, each dataset exhibits variations that remain within the associated statistical uncertainties. The rattle effect originates from momentum transfer between the fluid and solid and, thus, primarily depends on the frequency of collisions with the solid surface. It is therefore more sensitive to the local fluid density near the interface than to the adsorption energy itself. Nevertheless, since the adsorbed layer density strongly depends on ε, the rattle effect indirectly reflects variations in the solid/fluid interaction. These findings are consistent with previous results on thermal transport in zeolites subjected to fluid adsorption.¹¹ Moreover, the underlying phonon scattering at the solid/fluid interface is a function of the fluid adsorbed amount at the surface as will be rationalized in detail in the next section by considering a simple stochastic model for the phonon scattering rate.

5.1.2 Fluid contribution κ_α. Fig. 5(b) shows the fluid contribution κ^‖_α to the overall thermal conductivity κ as a function of the adsorbed amount n_a. As expected, the fluid contribution to the thermal conductivity is nearly independent of the solid/fluid interaction strength. The fluid contribution κ^‖_α scales as n_a². As show in Fig. S2 of the SI, there are three main contributions for κ^‖_α: (i) κ^‖_αα which is the dominant term that scales as n_a; (ii) κ^‖_α,αβ whose behavior strongly depends on the interaction strength ε value; (iii) the contribution from the heat current correlation between the solid and fluid, which results in a negative contribution that decreases κ^‖_α. Although the dominant term κ^‖_αα scales with n_a, the combined effect of the cross-correlation terms, particularly κ^‖_α,αβ and the negative solid–fluid contribution, introduces an additional density dependence. This leads to the total κ^‖_α to effectively scale as n_a² when all contributions are summed.

It is believed that the influence of the bulk-like region on κ^‖_α is minimal, given that the effective parallel thermal conductivity arises from surface and bulk-like contributions acting in parallel. Indeed, because κ_α is observed to be a unique function of the surface adsorbate density n_a (see Fig. 5(b)), the effective property appears to be primarily governed by the surface region. This suggests that fluid molecules outside the surface layer contribute weakly to heat transport. Moreover, the higher local density within the surface region, compared to the bulk-like region, enhances its thermal conductivity and further supports this interpretation. However, we note that, for the reasons provided above, while n_a is the controlling parameter discussed in this work, it is not possible to rule out that it also governs other confinement geometries.

5.1.3 Interface contribution κ_αβ. As illustrated in Fig. 5(c), the solid/fluid interaction appears to have a small effect on the thermal conductivity across their interface. However, the interfacial contribution, κ^‖_αβ, which acts as a thermal bridge between the solid and liquid, remains negligible compared to the contributions from the solid and fluid phases as shown in Fig. 5(a) and (b) (despite a slight increase of κ^‖_αβ with both n_a and ε). Notably, the main contribution to κ^‖_αβ arises from the cross-term κ_αβ,α while the other interfacial terms barely contribute. Physically, this suggests that direct heat transfer across the interface region constitutes only a small fraction of the total thermal conductivity. Instead, heat transport is dominated by conduction within the bulk solid and fluid phases – with the interface serving primarily as a boundary rather than an efficient heat conduction channel. This behavior can be attributed to the molecular structure and dynamics near the interface: although stronger solid/fluid affinity enhances adsorption and coupling at the interface, the interfacial region may exhibit reduced mobility that limits its thermal conductivity. Moreover, this thermal interface is the locus of the collisions between fluid and solid atoms – which leads to incoherent and, therefore, less efficient heat transport. In other words, heat is dissipated through unmatched vibration of solid and fluid atoms at interface, which lead to collisions in their trajectories and, hence, the rattle effect.¹¹

5.2 Phonon lifetime and the rattle effect

In insulating materials such as amorphous silica, the thermal conductivity κ arises from phonon transport. We use here a simple statistical physics approach based on the microscopic kinetic theory of phonons to quantify phonon scattering in heterogeneous (i.e., solid + fluid) systems. Specifically, we model thermal transport in the solid structure by considering its vibrational modes while accounting for the so-called rattle effect. We begin with the microscopic expression for the thermal conductivity, which includes all phonon modes and their associated heat capacities within the solid phase:Here, v(q, p) represents the phonon group velocity, τ(q, p) the phonon lifetime, and C_p(q, p) the heat capacity per unit of volume for a given phonon mode (a phonon mode is characterized by its wavevector q and polarization p). The term cos θ_x(q) accounts for the velocity projection along a given Cartesian direction. The second equality is simply obtained by assuming that the material is isotropic and that both the group velocity and phonon lifetime are constant –v(q, p) ∼ v and τ(q, p) ∼ τ.

To incorporate the effect of the adsorbed fluid n_a on phonon scattering via the rattle mechanism, we propose a stochastic model based on Matthiessen's rule.⁴⁰ In this model, we implicitly assume that scattering occurs through two mechanisms independent of each other: (1) intrinsic phonon scattering within the solid (which is characterized by a relaxation rate τ₀⁻¹ that is independent of n_a and ε) and (2) phonon scattering at the solid/fluid interface due to collisions between fluid and solid atoms (which is characterized by a relaxation rate τ_a⁻¹ that depends on n_a and ε). With this simple model, the number of phonons N(t + dt) not scattered at time t + dt is written as the number of phonons N(t) not scattered at time t minus the number of phonons scattered over the time interval dt: N(t + dt) = N(t)[1 − dt/τ₀ − dt/τ_a(n_a, ε)]. Solving this differential equation leads to an exponential decay for the number of phonons:^11,17

We now use this expression for the phonon lifetime τ(n_a,ε) to estimate the thermal conductivity via the second equality in eqn (8). By assuming that the solid's heat capacity and speed of sound are independent of n_a and by noting that κ^‖_β(n_a = 0) ∝ τ₀ while κ^‖_β(n_a, ε) ∝ τ(n_a, ε), we apply eqn (9) to obtain that the collision-induced phonon scattering rate (normalized to by τ₀):

The relaxation time appearing in the left-hand side of eqn (10) corresponds to the ratio of the intrinsic phonon relaxation time τ₀ to the relaxation time associated exclusively with fluid-induced scattering, τ_a(n_a,ε). The right-hand-side term arises from separating the total and fluid-induced scattering mechanisms through Matthiessen's rule.

We note that the phonon relaxation time τ introduced in eqn (8) (and further decomposed in τ₀ and τ_a, which ratio is given by eqn (10)) characterizes heat transport through the solid phase – and not through the overall hybrid (i.e. solid + fluid) system. As a result, changes in τ due to the presence of adsorbed fluid at the solid/fluid interface directly impacts the solid thermal conductivity κ_β. Fig. 5(d) shows the normalized phonon scattering rate, , as a function of the adsorbed amount n_a for the three different solid/fluid interaction strengths: 0.5ε₀, ε₀, and 2ε₀. The data shows that stronger solid-fluid interactions lead to a decrease in the phonon scattering rate. This result, which appears counterintuitive at first, can be rationalized as follows. Strong fluid/solid interactions lead to fluid atoms that are tightly physisorbed to the surface and, therefore, dynamically coupled to the solid molecular structure. As a result, such significant binding leads to frequent collisions between fluid and solid atoms which, in turn, leads to a decrease in the phonon lifetime and hence a decrease in the solid's heat conductivity. In fact, in agreement with prior observations, phonon scattering at the fluid/solid interface mostly depends on collisions between solid and fluid atoms and is therefore governed by momentum exchange at the solid surface.¹¹

These data reveal a strong dependence of on ε with a crossover in the scattering behavior upon increasing n_a. Here, the dilute regime refers to the limit in which adsorbed molecules contribute independently to phonon scattering – such that each fluid atom acts as an independent phonon scatterer, leading to a linear scaling of with n_a (even if intermolecular interactions within the adsorbed phase are non-negligible). This linear trend is confirmed by the dash-dotted line in Fig. 5(d). To quantify deviations from this ideal behavior corresponding to infinite dilution (Henry regime), we model the scattering rate as τ_a(n_a,ε)⁻¹ ∼ f(n_a,ε)τ⁰_a(ε)⁻¹, where τ⁰_a(ε) ∼ τ_a(n_a → 0, ε) is the reference scattering time in the dilute limit. For consistency, we define , with n⁰_a= 10 molec./nm² (however, we note that our results remain qualitatively independent of the exact choice as long as n⁰_a lies within the linear regime). The inset of Fig. 5(d) shows the resulting function . We expect f(n_a,ε) to grow linearly at low n_a as each fluid atom scatters phonon independently of each other. In fact, while the data do not firmly show this linear regime in Fig. 5(d) with n_a in the low adsorbed amount range, it is theoretically expected. Indeed, in the limit of nearly non-interacting fluid molecules (Henry regime of the adsorption isotherm), we expect each fluid molecule to contribute independently of each other to the rattle effect – therefore, leading to a linear regime f(n_a) ∼ n_a. In contrast, as can be observed in the inset of Fig. 5(d) upon increasing n_a, f(n_a) becomes non linear and tends to a plateau as saturation of the rattle effect occurs when the nanoconfined fluid becomes very crowded. For the interaction strength ε = 0.5ε₀, the function is close to unity over the evaluated range of n_a. This indicates that, for such weak interactions, the fluid-induced scattering rate is constant and close to that in the dilute regime (no cooperative effect between adsorbed molecules in terms of phonon scattering). Moreover, the corresponding rates exhibit larger relative statistical fluctuations and do not display the monotonic trends observed for the other interactions strengths ε₀ and 2ε₀.

Formally, f(n_a) describes the dependence on the adsorbed amount n_a of the scattering rate due to fluid molecules at the surface (the latter is normalized by the scattering rate in the limit of a vanishingly small number of fluid molecules ). Within the silica pore, f(n_a,ε) grows linearly at low n_a as each fluid atom scatters phonons independently of each other. On the other hand, as n_a increases, fluid/fluid interactions give rise to cooperative scattering which results in a increase in f(n_a,ε) followed by saturation as the interface becomes crowded/saturated at very large n_a. Interestingly, this transition from independent (linear) to cooperative (non linear) behavior of the so-called “rattle effect” is not observed for methane confined within zeolites,¹¹ where f(n_a) remains linear until saturation is reached. This qualitative difference likely stems from pore morphology/topology effects: severe confinement and very large surface areas in zeolites do not lead to collective fluid effects within the adsorbed phase as fluid/fluid interactions are negligible over fluid/solid interactions. In contrast, the open silica surface studied here leads to large fluid/fluid interactions at the surface, which result in a different dependence of the scattering rate on the adsorbed amount f(n_a).

Although differences in the phonon spectra for crystalline zeolite and amorphous silica are expected, we note that f(n_a,ε) relates directly to the scattering of phonons by fluid molecules at the solid/liquid interface regardless of the intrinsic features of each solid phase. In other words, as can be seen from our formal treatment using Matthiessen's rules, fluid-driven scattering of the phonon modes adds to the intrinsic phonon scattering independently of the type of phonon spectrum (crystalline versus amorphous material).

6 Conclusions

We conducted a quantitative assessment of the thermal conductivity in a prototypical system consisting of a nanoporous material filled with a confined fluid. For this, we simulated a system composed of a simple fluid confined in a realistic model of nanoporous silica made up of amorphous silica walls. By considering different solid/fluid interaction strengths ε and adsorbed amounts N_a, we investigate the parameters that give rise to the rattle effect. Despite the use of thick silica walls, we find that the rattle effect remains a non-negligible phenomenon that affects the overall thermal conductivity in such solid/fluid systems. We identify that the thermal conductivity and its fluid phase, solid phase, and solid/fluid interface contributions can be rationalized by taking into account the adsorbed density n_a. While the thermal conductivity of the solid phase κ^‖_β is decreased by the rattle effect, the thermal conductivity of the fluid phase κ^‖_α increases as the fluid density increases within the pore. The thermal conductivity contribution from both fluid and solid phases (κ^‖_α and κ^‖_β, respectively) follows a unique and general relation with the adsorbed amount – regardless of the interaction strength between the solid and fluid atoms. We also evaluate the role of the solid/fluid interface contribution κ^‖_αβ to the thermal conductivity. This term is a function of the adsorbed density n_a but also slightly dependent of the adsorbent/adsorbate interaction strength ε.

The connection between the thermal conductivity in solid/fluid systems and the fluid interfacial/adsorbed density suggests that the rattle mechanism is an important contribution. This microscopic phenomenon, which arises from phonon scattering induced by collisions of the fluid atoms with the solid surface, leads to a decrease in the solid's contribution to thermal conductivity. As for the fluid contribution to the overall thermal conductivity, we find that increasing n_a leads to a much larger number of collisions – which increases the virial (i.e. potential energy) contribution to the thermal conductivity. While our results suggest that the rattle effect is a general mechanism that should be taken into account in any solid/fluid systems, further evaluation is needed to assess its impact for more complex systems such as when polyatomic molecules or mixtures are considered. In practice, our findings highlight the inherent limitations of conventional mixing rules for porous materials systems. Such macroscopic approaches, which are routinely used in engineering applications, neglect the coupling between the solid and fluid phases and its consequences. As a result, describing thermal transport in such complex solid/fluid phases requires an atomistic-level approach to correctly capture the underlying mechanisms. Finally, the findings reported in the present result suggest that the solid/fluid coupling can be harnessed to engineer specific nanodevices with controlled effective thermal conductivity. Future studies may explore the impact of varying fluid types and nanopores geometries on the rattle effect as well as the effect of temperature and pressure on phonon scattering dynamics. These investigations could ultimately lead to the development of materials tailored for specific thermal management applications. This opens perspectives for the design of materials with tailored and optimized thermal properties.

Author contributions

B. C. and L. F. M. F conceived the research. N. F. S. carried out the molecular simulations. N. F. S., T. R. G., C. H. and B. C. analyzed the data and developed the model. All authors contributed to the analysis of the results. N. F. S and B. C. wrote the manuscript with inputs from all authors.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article, including the molecular configurations generated and/or analyzed in the framework of this study are available at Science Data Bank at https://doi.org/10.57760/sciencedb.30816.

Supplementary information (SI): additional data, corresponding to the force field parameters, comparison of simulation data with the classical mixing rules and details about the thermal conductivity contributions defined in the text. See DOI: https://doi.org/10.1039/d5nr04591f.

Acknowledgements

This work was partially supported by São Paulo Research Foundation (FAPESP) grant no. 2018/02713-8 and 2021/11105-4, CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico), and CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior). This work was also carried within the framework of a joint project HTD-POM (ANR-23-CE08-0038-01) funded by the French Agence Nationale pour la Recherche (ANR) and the Brazilian agency FAPESP (grant no. 2023/03170-6). All (or most of) the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities agency. N. F. S. and B. C. are grateful to J. L. Barrat (Univ. Grenoble Alpes) for stimulating discussions.

References

1. G. Sneddon, A. Greenaway and H. H. Yiu, Adv. Energy Mater., 2014, 4, 1301873.
2. E. Kianfar and H. Sayadi, Carbon Lett., 2022, 32, 1645–1669.
3. D. S. Sholl and R. P. Lively, Nature, 2016, 532, 435–437.
4. L. Bocquet, Nat. Mater., 2020, 19, 254–256.
5. N. Kavokine, R. R. Netz and L. Bocquet, Annu. Rev. Fluid Mech., 2021, 53, 377–410.
6. A. Schlaich, J.-L. Barrat and B. Coasne, Chem. Rev., 2025, 125, 2561–2624.
7. G. Gonella, E. H. Backus, Y. Nagata, D. J. Bonthuis, P. Loche, A. Schlaich, R. R. Netz, A. Kühnle, I. T. McCrum and M. T. Koper, et al., Nat. Rev. Chem., 2021, 5, 466–485.
8. Y. F. Han, S. Li, H.-D. Liu and W. Cui, Therm. Sci., 2020, 24, 3749–3756.
9. D. Raabe, Modell. Simul. Mater. Sci. Eng., 2004, 12, R13.
10. Y.-L. He and W.-Q. Tao, J. Heat Transfer, 2015, 137, 090801.
11. N. F. Souza, C. Picard, L. F. M. Franco and B. Coasne, J. Phys. Chem. B, 2024, 128, 2516–2527.
12. R. Alipour, M. Ghoranneviss, M. Mirzaee and A. Jafari, Heat Mass Transfer, 2014, 50, 1727–1735.
13. X. Liang, Chin. Sci. Bull., 2007, 52, 2457–2472.
14. G. Nolas, G. Slack, D. Morelli, T. Tritt and A. Ehrlich, J. Appl. Phys., 1996, 79, 4002–4008.
15. J. A. Thomas, R. M. Iutzi and A. J. McGaughey, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 045413.
16. E. J. Rosenbaum, N. J. English, J. K. Johnson, D. W. Shaw and R. P. Warzinski, J. Phys. Chem. B, 2007, 111, 13194–13205.
17. H. Babaei and C. E. Wilmer, Phys. Rev. Lett., 2016, 116, 025902.
18. H. Babaei, M. E. DeCoster, M. Jeong, Z. M. Hassan, T. Islamoglu, H. Baumgart, A. J. H. McGaughey, E. Redel, O. K. Farha, P. E. Hopkins, J. A. Malen and C. E. Wilmer, Nat. Commun., 2020, 11, 4010.
19. V. V. Murashov, J. Condens. Matter Phys., 1999, 11, 1261.
20. C.-Y. Chen and D. I. Kopelevich, Mol. Simul., 2008, 34, 155–167.
21. H. Mizuno, S. Mossa and J.-L. Barrat, Sci. Rep., 2015, 5, 14116.
22. S. R. De Groot and P. Mazur, in Non-equilibrium thermodynamics, Courier Corporation, 2013.
23. G. Paolini and G. Ciccotti, Phys. Rev. A, 1987, 35, 5156.
24. D. Surblys, H. Matsubara, G. Kikugawa and T. Ohara, Phys. Rev. E, 2019, 99, 051301.
25. B. Coasne, L. Viau and A. Vioux, J. Phys. Chem. Lett., 2011, 2, 1150–1154.
26. L. N. Ho, S. Clauzier, Y. Schuurman, D. Farrusseng and B. Coasne, J. Phys. Chem. Lett., 2013, 4, 2274–2278.
27. S. Hocine, R. Hartkamp, B. Siboulet, M. Duvail, B. Coasne, P. Turq and J.-F. Dufrêche, J. Phys. Chem. C, 2016, 120, 963–973.
28. L. Zhuravlev, Langmuir, 1987, 3, 316–318.
29. A. Carré, J. Horbach, S. Ispas and W. Kob, Europhys. Lett., 2008, 82, 17001.
30. B. L. Eggimann, A. J. Sunnarborg, H. D. Stern, A. P. Bliss and J. I. Siepmann, Mol. Simul., 2014, 40, 101–105.
31. R. T. Cygan, J.-J. Liang and A. G. Kalinichev, J. Phys. Chem. B, 2004, 108, 1255–1266.
32. H. A. Lorentz, Ann. Phys., 1881, 248, 127–136.
33. D. Berthelot, Compt. Rendus, 1898, 126, 15.
34. L. N. Ho, Y. Schuurman, D. Farrusseng and B. Coasne, J. Phys. Chem. C, 2015, 119, 21547–21554.
35. B. Coasne and P. Ugliengo, Langmuir, 2012, 28, 11131–11141.
36. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
37. N. J. English and S. T. John, Phys. Rev. Lett., 2009, 103, 015901.
38. L. Dalstein, E. Potapova and E. Tyrode, Phys. Chem. Chem. Phys., 2017, 19, 10343–10349.
39. S. Andersson and L. Dzhavadov, J. Condens. Matter Phys., 1992, 4, 6209.
40. N. W. Ashcroft and N. Mermin, Solid state, Holt Rinehart and Winston, 1976, vol. 1.

---