# Equilibrium and Non-Equilibrium Molecular Dynamics Simulation of Thermo-Osmosis: Enhanced Effects on Polarized Graphene Surfaces

Mehdi Ouadfel, Samy Merabia, Yasutaka Yamaguchi, Laurent Joly

► To cite this version:

Mehdi Ouadfel, Samy Merabia, Yasutaka Yamaguchi, Laurent Joly. Equilibrium and Non-Equilibrium Molecular Dynamics Simulation of Thermo-Osmosis: Enhanced Effects on Polarized Graphene Surfaces. 2024. [hal-04525549](https://hal.science/hal-04525549)

HAL Id: hal-04525549

https://hal.science/hal-04525549

Preprint submitted on 29 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Equilibrium and Non-Equilibrium Molecular Dynamics Simulation of Thermo-Osmosis: Enhanced Effects on Polarized Graphene Surfaces

Mehdi Ouadfel,<fnref n="1" />, a) Samy Merabia,<fnref n="1" />, b) Yasutaka Yamaguchi,<fnref n="2" />,<fnref n="3" />, c) and Laurent Joly<fnref n="1" />, d)

(Dated: 11 July 2024)

Thermo-osmotic flows, generated by applying a thermal gradient along a liquid-solid interface, could be harnessed to convert waste heat into electricity. While this phenomenon has been known for almost a century, there is a crucial need to gain a better understanding of the molecular origins of thermo-osmosis. In this paper, we start by detailing the multiple contributions to thermo-osmosis. We then showcase three approaches to compute the thermo-osmotic coefficient using molecular dynamics; a first method based on the computation of the interfacial enthalpy excess and Derjaguin's theoretical framework, a second approach based on the computation of the interfacial entropy excess using the so-called dry-surface method, and a novel non-equilibrium method to compute the thermo-osmotic coefficient in a periodic channel. We show that the three methods align with each other, in particular for smooth surfaces. In addition, for a polarized graphene-water interface, we observe large variations of thermo-osmotic responses, and multiple changes in flow direction with increasing surface charge. Overall, this study showcases the versatility of osmotic flows and calls for experimental investigation of thermo-osmotic behavior in the vicinity of charged surfaces.

## I. INTRODUCTION

Nanofluidic systems offer great promises for energy harvesting<fnref n="1" />⁻<fnref n="3" />. At the nanoscale, osmotic flows, which are flows generated by a thermodynamic gradient along solid-liquid interfaces, can be used to generate electricity. Indeed, if one considers an aqueous electrolyte near a charged surface, the ions in the liquid will move under the influence of the electric field generated by a surface to form the electrical double layer (EDL), a layer of electrically charged liquid⁴⁻⁶. The advection of the EDL by osmotic flows creates an electric current. For example, diffusio-osmosis, a phenomenon wherein the flow is generated by a salinity gradient, creates a so-called diffusiosmotic current⁷⁻⁹.

One can also generate an osmotic flow by applying a thermal gradient along an interface, which is referred to as thermo-osmosis¹⁰⁻¹². Derjaguin and colleagues developed a theoretical framework for thermo-osmosis, which relates the thermo-osmotic flow to the interfacial enthalpy excess¹³⁻¹⁵. This approach, based on linear irreversible thermodynamics and continuum hydrodynamics, faces several issues. The enthalpy excess, originating from the interfacial layer (interacting with the wall), typically spans a range on the order of 1 nm. In this context, the validity of the equation of continuum hydrodynamics becomes questionable¹⁶. Moreover, Derjaguin's approach assumes a constant viscosity of the fluid near the surface. While this assumption is reasonable for hydrophobic

surfaces¹⁷⁻¹⁹, it is generally not true¹⁸⁻²¹. Moreover, the local enthalpy excess is calculated using the pressure tensor, which is not uniquely defined when considering heterogeneous liquids at interfaces²².

Recently, numerous theoretical works brought together the macroscopic description of thermo-osmosis with a molecular description of the interfacial region²⁰,²³⁻²⁵. Using Onsager's reciprocity relation, which links the flow generated by a thermal gradient to the heat flux due to a pressure gradient, these studies have been able to study thermo-osmosis *via* the 'mechano-caloric' route²⁵⁻²⁸. A good agreement was found between this approach and Derjaguin's approach for model systems²⁶,²⁹. However, whether Derjaguin's framework can be applied at the molecular scale to predict thermo-osmotic transport remains under debate.

In this paper, we investigate the thermo-osmotic coefficient of realistic solid-liquid interfaces, i.e., water on graphitic surfaces, using molecular dynamics (MD) simulations. We begin by detailing the different contributions to thermo-osmosis. Subsequently, we calculate the thermo-osmotic coefficient using three different methods. First, we employ a previously discussed approach based on calculating the enthalpy excess at equilibrium and Derjaguin's theoretical framework³⁰. Second, we introduce an alternative approach that involves computing the solid-liquid entropy excess at equilibrium using the so-called dry-surface method (DSM). This method enables the calculation of the solid-liquid work of adhesion, which is linked to the entropy excess, through thermodynamic integration along the solid-liquid interaction parameter.³¹⁻³³. Finally, we have developed a non-equilibrium MD-based method to compute the thermo-osmotic response in a periodic channel, in which we directly apply a thermal gradient and compute the thermo-osmotic response.

<footnote n="1">¹) Université Claude Bernard Lyon 1, CNRS, Institut Lumière Matière, UMR5306, F69100 Villeurbanne, France</footnote>

<footnote n="2">²) Department of Mechanical Engineering, Osaka University, 2-1 Yamadaoka, Suita 565-0871, Japan</footnote>

<footnote n="3">³) Water Frontier Research Center (WaTUS), Research Institute for Science & Technology, Tokyo University of Science, 1-3 Kagurazaka, Shinjuku-ku, Tokyo 162-8601, Japan</footnote>

<footnote>a) Electronic mail: mehdi.ouadfel@univ-lyon1.fr</footnote>

<footnote>b) Electronic mail: samy.merabia@univ-lyon1.fr</footnote>

<footnote>c) Electronic mail: yamaguchi@mech.eng.osaka-u.ac.jp</footnote>

<footnote>d) Electronic mail: laurent.joly@univ-lyon1.fr</footnote>
# Equilibrium and Non-Equilibrium Molecular Dynamics Simulation of Thermo-Osmosis: Enhanced Effects on Polarized Graphene Surfaces

Mehdi Ouadfel,<fnref n="1" />, a) Samy Merabia,<fnref n="1" />, b) Yasutaka Yamaguchi,<fnref n="2" />,<fnref n="3" />, c) and Laurent Joly<fnref n="1" />, d)

(Dated: 29 March 2024)

Thermo-osmotic flows, generated by applying a thermal gradient along a liquid-solid interface, could be harnessed to convert waste heat into electricity. While this phenomenon has been known for almost a century, there is a crucial need to gain a better understanding of the molecular origins of thermo-osmosis. In this paper, we start by detailing the multiple contributions to thermo-osmosis. We then showcase three approaches to compute the thermo-osmotic coefficient using molecular dynamics; a first method based on the computation of the interfacial enthalpy excess and Derjaguin's theoretical framework, a second approach based on the computation of the interfacial entropy excess using the so-called dry-surface method, and a novel non-equilibrium method to compute the thermo-osmotic coefficient in a periodic channel. We show that the three methods align with each other, in particular for smooth surfaces. In addition, for a polarized graphene-water interface, we observe large variations of thermo-osmotic responses, and multiple changes in flow direction with increasing surface charge. Overall, this study showcases the versatility of osmotic flows and calls for experimental investigation of thermo-osmotic behavior in the vicinity of charged surfaces.

## I. INTRODUCTION

Nanofluidic systems offer great promises for energy harvesting<fnref n="1" />⁻<fnref n="3" />. At the nanoscale, osmotic flows, which are flows generated by a thermodynamic gradient along solid-liquid interfaces, can be used to generate electricity. Indeed, if one considers an aqueous electrolyte near a charged surface, the ions in the liquid will move under the influence of the electric field generated by a surface to form the electrical double layer (EDL), a layer of electrically charged liquid⁴⁻⁶. The advection of the EDL by osmotic flows creates an electric current. For example, diffusio-osmosis, a phenomenon wherein the flow is generated by a salinity gradient, creates a so-called diffusio-osmotic current⁷⁻⁹.

One can also generate an osmotic flow by applying a thermal gradient along an interface, which is referred to as thermo-osmosis¹⁰⁻¹². Derjaguin and colleagues developed a theoretical framework for thermo-osmosis, which relates the thermo-osmotic flow to the interfacial enthalpy excess¹³⁻¹⁵. This approach, based on linear irreversible thermodynamics and continuum hydrodynamics, faces several issues. The enthalpy excess, originating from the interfacial layer (interacting with the wall), typically spans a range on the order of 1 nm. In this context, the validity of the equation of continuum hydrodynamics becomes questionable¹⁶. Moreover, Derjaguin's approach assumes a constant viscosity of the fluid near the surface. While this assumption is reasonable for hydrophobic

surfaces¹⁷⁻¹⁹, it is generally not true¹⁸⁻²¹. Moreover, the local enthalpy excess is calculated using the pressure tensor, which is not uniquely defined when considering heterogeneous liquids at interfaces²².

Recently, numerous theoretical works brought together the macroscopic description of thermo-osmosis with a molecular description of the interfacial region²⁰,²³⁻²⁵. Using Onsager's reciprocity relation, which links the flow generated by a thermal gradient to the heat flux due to a pressure gradient, these studies have been able to study thermo-osmosis *via* the 'mechano-caloric' route²⁵⁻²⁸. A good agreement was found between this approach and Derjaguin's approach for model systems²⁶,²⁹. However, whether Derjaguin's framework can be applied at the molecular scale to predict thermo-osmotic transport remains under debate.

In this paper, we investigate the thermo-osmotic coefficient of realistic solid-liquid interfaces, i.e., water on graphitic surfaces, using molecular dynamics (MD) simulations. We begin by detailing the different contributions to thermo-osmosis. Subsequently, we calculate the thermo-osmotic coefficient using three different methods. First, we employ a previously discussed approach based on calculating the enthalpy excess at equilibrium and Derjaguin's theoretical framework³⁰. Second, we introduce an alternative approach that involves computing the solid-liquid entropy excess at equilibrium using the so-called dry-surface method (DSM). This method enables the calculation of the solid-liquid work of adhesion, which is linked to the entropy excess, through thermodynamic integration along the solid-liquid interaction parameters.³¹⁻³³. Finally, we have developed a non-equilibrium MD-based method to compute the thermo-osmotic response in a periodic channel, in which we directly apply a thermal gradient and compute the thermo-osmotic response.

<footnote n="1">¹) Universite Claude Bernard Lyon 1, CNRS, Institut Lumière Matière, UMR5306, F69100 Villeurbanne, France</footnote>

<footnote n="2">²) Department of Mechanical Engineering, Osaka University, 2-1 Yamadaoka, Suita 565-0871, Japan</footnote>

<footnote n="3">³) Water Frontier Research Center (WaTUS), Research Institute for Science & Technology, Tokyo University of Science, 1-3 Kagurazaka, Shinjuku-ku, Tokyo 162-8601, Japan</footnote>

<footnote>a) Electronic mail: mehdi.ouadfel@univ-lyon1.fr</footnote>

<footnote>b) Electronic mail: samy.merabia@univ-lyon1.fr</footnote>

<footnote>c) Electronic mail: yamaguchi@mech.eng.osaka-u.ac.jp</footnote>

<footnote>d) Electronic mail: laurent.joly@univ-lyon1.fr</footnote>
osmotic velocity. We find a remarkable agreement between the different approaches, confirming that it is possible to apply Derjaguin theory at the molecular scale. In addition, we observe large variations of thermo-osmotic responses on graphitic surfaces, large responses in polarized graphene, and multiple changes in flow direction with increasing surface charge.

## II. THEORY

In this section, we will showcase the different contributions to thermo-osmosis. Following the work of Derjaguin and collaborators¹³⁻¹⁵, taking into account solid-liquid slip, one can solve the Stokes equation and obtain the velocity of the fluid far from the surface, i.e. the osmotic velocity¹²,²⁵,²⁶:

$$v_{\text{osm}} = \frac{1}{\eta} \int_{0}^{\infty} (z + b) f_{x}(z) dz, \quad (1)$$

in which $z = 0$ is the position of the surface, $z = \infty$ corresponds to the middle of the channel, $\eta$ is the viscosity, considered homogeneous, $b$ is the slip length³⁴, which quantifies solid-liquid slip, and $f_x$ is the force density exerted on the fluid in the $x$ direction. When a thermal gradient is applied along the interface, this force density is written¹¹

$$f_x(z) = -\delta h(z) \frac{\nabla T}{T} - \frac{E(z)^2}{2} \nabla \varepsilon. \quad (2)$$

The first term can be derived from local thermodynamic equilibrium²⁰,³⁵, it is the force density generated by a temperature gradient, wherein $\delta h(z)$ is the enthalpy excess density and $\nabla T$ is the temperature gradient. The second term, which we could refer to as temperature-dependent permittivity-osmosis, arises from a gradient of dielectric permittivity $\varepsilon$. Note that this term, at constant pressure, depends solely on the variation in permittivity induced by temperature³⁶. Here, $E(z)$ is the electric field, which is approximately equal to its transverse component $E_z(z)$ close to an electrically charged surface. Considering only the ion contribution to the enthalpy excess, using the Poisson-Boltzmann (PB) framework (here for a symmetric Z:Z salt), the force density becomes:

$$f_x(z) = - \left[ -\varepsilon V(z) \frac{d^2 V}{dz^2} + \frac{\varepsilon}{2} \left( \frac{dV}{dz} \right)^2 \right] \frac{\nabla T}{T} - \frac{1}{2} \left( \frac{dV}{dz} \right)^2 \nabla \varepsilon, \quad (3)$$

where $V(z)$, the electrostatic potential, is given by the PB theory²⁵,³⁷,³⁸. One can then integrate Eq. (1) to obtain the thermo-osmotic velocity

$$v_{\text{to}}^{\text{PB}} = v_{\text{to},\nabla T}^{\text{PB}} + v_{\text{to},\nabla \varepsilon}^{\text{PB}}, \quad (4)$$

with

$$v_{\text{to},\nabla T}^{\text{PB}} = -\frac{\nabla T/T}{2\pi \ell_{\text{B}} \eta \beta} \left\{ \begin{aligned} & -3 \ln(1 - \gamma^2) - \operatorname{asinh}^2(x) \\ & + \frac{b}{\lambda_{\text{D}}} \left[ 3x|\gamma| - 2x \operatorname{asinh}(x) \right] \end{aligned} \right\}, \quad (5)$$

and

$$v_{\text{to},\nabla \varepsilon}^{\text{PB}} = -\frac{\nabla \varepsilon / \varepsilon}{2\pi \ell_{\text{B}} \eta \beta} \left[ -\ln(1 - \gamma^2) + \frac{b}{\lambda_{\text{D}}} x |\gamma| \right], \quad (6)$$

where $\ell_{\text{B}} = \beta q^2/(4\pi\varepsilon)$ is the Bjerrum length, with $\beta = 1/(k_{\text{B}}T)$ and $q = Z\varepsilon$ the absolute ionic charge, $\lambda_{\text{D}} = 1/(\sqrt{8\pi\ell_{\text{B}}n_0})$ is the Debye length, in which $n_0$ is the bulk ion concentration, and $x = \lambda_{\text{D}}/\ell_{\text{GC}}$, where $\ell_{\text{GC}} = q/(2\pi\ell_{\text{B}}|\Sigma|)$ is the Gouy-Chapman length, with $\Sigma$ the surface charge density. Finally

$$\gamma = \frac{\operatorname{sgn}(\Sigma)}{x} \left[ -1 + \sqrt{1 + x^2} \right]. \quad (7)$$

The permittivity gradient can be expressed as a function of the thermal gradient¹¹:

$$\frac{\nabla \varepsilon}{\varepsilon} = -\tau \frac{\nabla T}{T}, \quad (8)$$

with $\tau = 1.4$ for water at room temperature. Using this relation, one can compute the two contributions to thermo-osmosis, which we quantify by using the thermosmotic coefficient:

$$M_{\text{to}} = -\frac{v_{\text{to}}}{\nabla T/T}. \quad (9)$$

Consequently, $M_{\text{to},\nabla T}^{\text{PB}} = -v_{\text{to},\nabla T}^{\text{PB}}/(\nabla T/T)$ and $M_{\text{to},\nabla \varepsilon}^{\text{PB}} = -v_{\text{to},\nabla \varepsilon}^{\text{PB}}/(\nabla T/T)$. Figure 1 shows that $M_{\text{to},\nabla T}^{\text{PB}}$ and $M_{\text{to},\nabla \varepsilon}^{\text{PB}}$ are comparable in the PB theory. However, it has been shown that the PB contribution to enthalpy excess can be negligible compared to that of water near some surfaces²⁵,³⁰. In this case, the effect of the permittivity gradient is also negligible compared to that of the enthalpy excess term of water.

In the general case, for a mixture of particles, the enthalpy excess density is defined as²⁰:

$$\delta h(z) = \sum_{i} n_{i}(z) [h_{i}(z) - h_{i}^{\text{bulk}}], \quad (10)$$

with $n_i$ and $h_i$ respectively the number density and the enthalpy per particle of the species $i$, and the superscript bulk refers to a quantity far from the surface, where it is homogeneous. The enthalpy excess is given by:

$$\Delta H = \int_{0}^{\infty} \delta h(z) dz. \quad (11)$$

Note that this definition differs from the classical definition of an excess quantity in the surface thermodynamics
osmotic velocity. We find a remarkable agreement between the different approaches, confirming that it is possible to apply Derjaguin theory at the molecular scale. In addition, we observe large variations of thermo-osmotic responses on graphitic surfaces, large responses in polarized graphene, and multiple changes in flow direction with increasing surface charge.

## II. THEORY

In this section, we will showcase the different contributions to thermo-osmosis. Taking into account solid-liquid slip, one can solve Stokes equation and obtain the velocity of the fluid far from the surface, i.e. the osmotic velocity<sup>12,25,26</sup>:

$$v_{\text{osm}} = \frac{1}{\eta} \int_{0}^{\infty} (z + b) f_{x}(z) dz, \quad (1)$$

in which $z = 0$ is the position of the surface, $z = \infty$ corresponds to the middle of the channel, $\eta$ is the viscosity, considered homogeneous, $b$ is the slip length<sup>34</sup>, which quantifies solid-liquid slip, and $f_x$ is the force density exerted on the fluid in the $x$ direction. When a thermal gradient is applied along the interface, this force density is written<sup>11</sup>

$$f_x(z) = -\delta h(z) \frac{\nabla T}{T} - \frac{E(z)^2}{2} \nabla \varepsilon. \quad (2)$$

The first term can be derived from local thermodynamic equilibrium<sup>20,35</sup>, it is the force density generated by a temperature gradient, wherein $\delta h(z)$ is the enthalpy excess density and $\nabla T$ is the temperature gradient. The second term, which we could refer to as temperature-dependent permittivity-osmosis, arises from a gradient of dielectric permittivity $\varepsilon$. Note that this term, at constant pressure, depends solely on the variation in permittivity induced by temperature<sup>36</sup>. Here, $E(z)$ is the electric field, which is approximately equal to its transverse component $E_z(z)$ close to an electrically charged surface. Considening only the ion contribution to the enthalpy excess, using the Poisson-Boltzmann (PB) framework (here for a symmetric Z:Z salt), the force density becomes:

$$f_x(z) = - \left[ -\varepsilon V(z) \frac{d^2 V}{dz^2} + \frac{\varepsilon}{2} \left( \frac{dV}{dz} \right)^2 \right] \frac{\nabla T}{T} - \frac{1}{2} \left( \frac{dV}{dz} \right)^2 \nabla \varepsilon, \quad (3)$$

where $V(z)$, the electrostatic potential, is given by the PB theory<sup>25,37,38</sup>. One can then integrate Eq. (1) to obtain the thermo-osmotic velocity

$$v_{\text{to}}^{\text{PB}} = v_{\text{to},\nabla T}^{\text{PB}} + v_{\text{to},\nabla \varepsilon}^{\text{PB}}, \quad (4)$$

with

$$v_{\text{to},\nabla T}^{\text{PB}} = -\frac{\nabla T/T}{2\pi \ell_{\text{B}} \eta \beta} \left\{ \begin{aligned} & -3 \ln(1 - \gamma^2) - \operatorname{asinh}^2(x) \\ & + \frac{b}{\lambda_{\text{D}}} \left[ 3x|\gamma| - 2x \operatorname{asinh}(x) \right] \end{aligned} \right\}, \quad (5)$$

and

$$v_{\text{to},\nabla \varepsilon}^{\text{PB}} = -\frac{\nabla \varepsilon / \varepsilon}{2\pi \ell_{\text{B}} \eta \beta} \left[ -\ln(1 - \gamma^2) + \frac{b}{\lambda_{\text{D}}} x |\gamma| \right], \quad (6)$$

where $\ell_{\text{B}} = \beta q^2/(4\pi\varepsilon)$ is the Bjerrum length, with $\beta = 1/(k_{\text{B}}T)$ and $q = Ze$ the absolute ionic charge, $\lambda_{\text{D}} = 1/(\sqrt{8\pi\ell_{\text{B}}n_0})$ is the Debye length, in which $n_0$ is the bulk ion concentration, and $x = \lambda_{\text{D}}/\ell_{\text{GC}}$, where $\ell_{\text{GC}} = q/(2\pi\ell_{\text{B}}|\Sigma|)$ is the Gouy-Chapman length, with $\Sigma$ the surface charge density. Finally

$$\gamma = \frac{\operatorname{sgn}(\Sigma)}{x} \left[ -1 + \sqrt{1 + x^2} \right]. \quad (7)$$

The permittivity gradient can be expressed as a function of the thermal gradient<sup>11</sup>:

$$\frac{\nabla \varepsilon}{\varepsilon} = -\tau \frac{\nabla T}{T}, \quad (8)$$

with $\tau = 1.4$ for water at room temperature. Using this relation, one can compute the two contributions to thermo-osmosis, which we quantify by using the thermosmotic coefficient:

$$M_{\text{to}} = -\frac{v_{\text{to}}}{\nabla T/T}. \quad (9)$$

Consequently, $M_{\text{to},\nabla T}^{\text{PB}} = -v_{\text{to},\nabla T}^{\text{PB}}/(\nabla T/T)$ and $M_{\text{to},\nabla \varepsilon}^{\text{PB}} = -v_{\text{to},\nabla \varepsilon}^{\text{PB}}/(\nabla T/T)$. Figure 1 shows that $M_{\text{to},\nabla T}^{\text{PB}}$ and $M_{\text{to},\nabla \varepsilon}^{\text{PB}}$ are comparable in the PB theory. However, it has been shown that the PB contribution to enthalpy excess can be negligible compared to that of water near some surfaces<sup>25,30</sup>. In this case, the effect of the permittivity gradient is also negligible compared to that of the enthalpy excess term of water.

In the general case, for a mixture of particles, the enthalpy excess density is defined as<sup>20</sup>:

$$\delta h(z) = \sum_{i} n_{i}(z) [h_{i}(z) - h_{i}^{\text{bulk}}], \quad (10)$$

with $n_i$ and $h_i$ respectively the number density and the enthalpy per particle of the species $i$, and the superscript bulk refers to a quantity far from the surface, where it is homogeneous. The enthalpy excess is given by:

$$\Delta H = \int_{0}^{\infty} \delta h(z) dz. \quad (11)$$

Note that this definition differs from the classical definition of an excess quantity in the surface thermodynamics
framework$^{39}$. One can then rewrite $M_{to,\nabla T}$ as:

$$M_{to,\nabla T} = -\frac{v_{to,\nabla T}}{\nabla T/T} \qquad (12)$$

$$= \frac{1}{\eta} \int_{0}^{\infty} (z+b)\delta h(z) dz \qquad (13)$$

$$= \frac{\Delta H}{\eta}(\lambda_h + b), \qquad (14)$$

where

$$\lambda_h = \frac{1}{\Delta H} \int_{0}^{\infty} z \delta h(z) dz \qquad (15)$$

denotes the extent of the layer where the liquid interacts with the wall, and is typically around 7 Å$^{30}$. Thus, $M_{to}$ can be estimated by computing $\Delta H$ at equilibrium, $\eta$, and $b$.

Another mean of computing the thermo-osmotic coefficient is to calculate the entropy excess generated by the interactions between the liquid and the surface. Indeed one can relate the enthalpy excess to the entropy excess. Introducing $\mu_i$, the chemical potential of the species $i$, and $s_i$, the entropy per particle of the species $i$, one can use the relation $\mu_i = h_i - Ts_i$ and the fact that $\mu_i$ and $T$ are homogeneous along $z$ to show that$^{20}$:

$$\delta h(z) = T\delta s(z), \qquad (16)$$

with the entropy excess density defined as:

$$\delta s(z) = \sum_{i} n_{i}(z)[s_{i}(z) - s_{i}^{\text{bulk}}]. \qquad (17)$$

The entropy excess thus writes:

$$\Delta S = \int_{0}^{\infty} \delta s(z) dz, \qquad (18)$$

and it follows that

$$\Delta H = T\Delta S. \qquad (19)$$

Therefore, one can write:

$$M_{to} = \frac{T\Delta S}{\eta}(b + \lambda_h). \qquad (20)$$

One can relate this result to the flow generated by a Marangoni stress. Indeed, the Marangoni stress can be written as a function of the enthalpy excess (and thus the entropy excess)$^{23}$:

$$\nabla\gamma = \left(\frac{\partial\gamma}{\partial T}\right)_P \nabla T = -\Delta H \frac{\nabla T}{T} = -\Delta S \nabla T, \qquad (21)$$

and the osmotic velocity becomes

$$v_{to} = \frac{\nabla\gamma}{\eta}(\lambda_h + b). \qquad (22)$$

This result illustrates an equivalence between Marangoni and thermo-osmotic flows. On low friction surfaces where $\lambda_h \ll b$, the velocity scales as $\sim \nabla\gamma b/\eta$, which is reminiscent of the scaling law of the Marangoni flow for superhydrophobic surfaces$^{40}$.

<figure><img src="image_1.png" /><figcaption>FIG. 1. Comparison of the contributions of the thermal and permittivity gradients to thermo-osmosis in the Poisson-Boltzmann theory for a polarized graphene surface, calculated for a salt concentration of 0.3 M and a slip length of 50 nm, representative of the graphene systems considered in the MD simulations.</figcaption></figure>

### III. METHODS

#### A. Systems and force field

We computed the thermo-osmotic coefficient of an aqueous electrolyte confined between graphitic walls using MD with the LAMMPS package$^{41}$. We used sodium chloride (NaCl) as the salt, with a bulk concentration $n_0 \sim 0.3$ M, corresponding to a Debye length $\lambda_D \sim 5.7$ Å. We simulated the aqueous electrolyte using the Madrid-2019 force field$^{42}$, based on TIP4P/2005$^{43}$, a rigid non-polarizable water model, and scaled charged for the ions, $|q_{Na}| = |q_{Cl}| = 0.85$ e. The walls were kept frozen. Unless specified, for water-carbon interactions, we took the parameters from Ref. 44, $\epsilon_{CO} = 0.114$ kcal/mol and $\sigma_{CO} = 3.28$ Å, where C and O correspond to the carbon and oxygen atoms, respectively; this corresponds to a contact angle of around 85°, see Appendix B. To assess our theoretical predictions, we also considered graphene-like surfaces with different wetting properties, which was achieved by modifying $\epsilon_{CO}$, see Appendix B for details. We applied the Lorentz-Berthelot mixing rules to set the carbon-ions interaction parameters. We truncated the LJ interactions as well as the Coulombic interactions at 10 Å. For the Coulombic interactions, we used the particle-particle particle-mesh (PPPM) solver for long-range corrections, with a relative error in forces of $10^{-4}$. We used periodic boundary conditions along the $x$ and $y$ directions.

We replicated the experimental structure of graphene, with an inter-atomic distance of 1.42 Å$^{45}$. To study the impact of surface defects on thermo-osmotic flows, we considered two types of graphene surfaces (Fig. 2): pristine graphene (PG) and graphene oxide (GO). We modeled GO by randomly distributing hydroxyl groups on the
framework$^{39}$. One can then rewrite $M_{to,\nabla T}$ as:

$$M_{to,\nabla T} = -\frac{v_{to,\nabla T}}{\nabla T/T} \qquad (12)$$

$$= \frac{1}{\eta} \int_{0}^{\infty} (z+b)\delta h(z) dz \qquad (13)$$

$$= \frac{\Delta H}{\eta}(\lambda_h + b), \qquad (14)$$

where

$$\lambda_h = \frac{1}{\Delta H} \int_{0}^{\infty} z \delta h(z) dz \qquad (15)$$

denotes the extent of the layer where the liquid interacts with the wall, and is around $7\AA^{30}$. Thus, $M_{to}$ can be estimated by computing $\Delta H$ at equilibrium, $\eta$, and $b$.

Another mean of computing the thermo-osmotic coefficient is to calculate the entropy excess generated by the interactions between the liquid and the surface. Indeed one can relate the enthalpy excess to the entropy excess. Introducing $\mu_i$, the chemical potential of the species $i$, and $s_i$, the entropy per particle of the species $i$, one can use the relation $\mu_i = h_i - T s_i$ and the fact that $\mu_i$ and $T$ are homogeneous along $z$ to show that$^{20}$:

$$\delta h(z) = T \delta s(z), \qquad (16)$$

with the entropy excess density defined as:

$$\delta s(z) = \sum_{i} n_{i}(z)[s_{i}(z) - s_{i}^{\text{bulk}}]. \qquad (17)$$

The entropy excess thus writes:

$$\Delta S = \int_{0}^{\infty} \delta s(z) dz, \qquad (18)$$

and it follows that

$$\Delta H = T \Delta S. \qquad (19)$$

Therefore, one can write:

$$M_{to} = \frac{T \Delta S}{\eta}(b + \lambda_h). \qquad (20)$$

One can relate this result to the flow generated by a Marangoni stress. Indeed, the Marangoni stress can be written as a function of the enthalpy excess (and thus the entropy excess)$^{23}$:

$$\nabla \gamma = \left(\frac{\partial \gamma}{\partial T}\right)_P \nabla T = -\Delta H \frac{\nabla T}{T} = -\Delta S \nabla T, \qquad (21)$$

and the osmotic velocity becomes

$$v_{to} = \frac{\nabla \gamma}{\eta}(\lambda_h + b). \qquad (22)$$

This result illustrates an equivalence between Marangoni and thermo-osmotic flows. On low friction surfaces where $\lambda_h \ll b$, the velocity scales as $\sim \nabla \gamma b/\eta$, which is reminiscent of the scaling law of the Marangoni flow for superhydrophobic surfaces$^{40}$.

<figure><img src="image_2.png" /><figcaption>FIG. 1. Comparison of the contributions of the thermal and permittivity gradients to thermo-osmosis in the Poisson-Boltzmann theory for a polarized graphene surface.</figcaption></figure>

### III. METHODS

#### A. Systems and force field

We computed the thermo-osmotic coefficient of an aqueous electrolyte confined between graphitic walls using MD with the LAMMPS package$^{41}$. We used sodium chloride (NaCl) as the salt, with a bulk concentration $n_0 \sim 0.3$ M, corresponding to a Debye length $\lambda_D \sim 5.7$ Å. We simulated the aqueous electrolyte using the Madrid-2019 force field$^{42}$, based on TIP4P/2005$^{43}$, a rigid nonpolarizable water model, and scaled charged for the ions, $|q_{Na}| = |q_{Cl}| = 0.85$ e. The walls were kept frozen. For water-graphene interactions, we took the parameters from Ref. 44, $\epsilon_{CO} = 0.114$ kcal/mol and $\sigma_{CO} = 3.28$ Å, where C and O correspond to the carbon and oxygen atoms, respectively. To assess our theoretical predictions, we also considered graphene-like surfaces with different wetting properties, which was achieved by modifying $\epsilon_{CO}$. We applied Lorentz-Berthelot mixing rules to set the carbon-ions interaction parameters. We truncated the LJ interactions as well as the Coulombic interactions at $10$ Å. For the Coulombic interactions, we used the particle-particle particle-mesh (PPPM) solver for long-range corrections, with a relative error in forces of $10^{-4}$. We used periodic boundary conditions along the $x$ and $y$ directions.

We replicated the experimental structure of graphene, with an inter-atomic distance of $1.42$ Å$^{45}$. To study the impact of surface defects on thermo-osmotic flows, we considered two types of graphene surfaces (Fig. 2): pristine graphene (PG) and graphene oxide (GO). We modeled GO by randomly distributing hydroxyl groups on the surface, with various oxidation rates $\xi = N_O/N_C$, where $N_O$ and $N_C$ are the number of hydroxyl and carbon atoms on the surface, respectively. This approach aligns with established methodologies in previous MD papers$^{46-48}$.
<figure><img src="image_2.png" /><figcaption>FIG. 2. Different types of graphitic surfaces considered: a neutral pristine graphene (PG) surface (a), a polarized PG surface (b), a neutral graphene oxide (GO) surface (c), and an electrically charged GO surface, for which the hydrogen from the hydroxyl groups are removed to create charged groups (d). The snapshots were produced with the software Ovito$^{51}$.</figcaption></figure>

surface, with various oxidation rates $\xi = N_O/N_C$, where $N_O$ and $N_C$ are the number of hydroxyl and carbon atoms on the surface, respectively. This approach aligns with established methodologies in previous MD papers$^{46–48}$. The interaction parameters for these hydroxyl groups were based on the phenol parameters within the all-atom optimized potentials for liquid simulation (OPLS-AA) force field$^{49}$.

The graphene sheets were either kept neutral or charged. To charge the PG surface uniformly, we assigned the same charge to all carbon atoms $q = \Sigma S/N_{wall}$, where S is the surface area of the sheet, and $N_{wall}$ is the number of carbon atoms, resulting in a surface charge density $\Sigma$. To charge the GO, we remove the hydrogen atoms from the hydroxyl groups. In practice, this was achieved by setting the charge of the carbon, oxygen and hydrogen of the hydroxyl groups to $q_{C\_phenol} = 0$ e, $q_{O\_phenol} = -0.85$ e and $q_{H\_phenol} = 0$ e, respectively, to match the rescaled charge of the ions. Counter-ions were added to the system to keep it electrically neutral.

In all simulations, we froze the bottom surface. We used the top surface as a piston during an equilibration phase that lasted 1 ns, before fixing it at its equilibrium position to set the pressure to 1 atm, for all simulations except the DSM (see details of DSM simulations in Section III C).

For the systems with a temperature gradient, we used a NVE time integration, and fixed the fluid temperature at 280 and 360 K for the cold and hot region, respectively, via a canonical sampling thermostat that uses global velocity rescaling with Hamiltonian dynamics$^{50}$, applied in the y and z directions, following a previous thermo-osmosis study$^{24}$, with a damping time of 200 fs. The mean temperature was 320 K. For other systems, we used a Nosé-Hoover thermostat to set the temperature to 320 K with a damping time of 200 fs. For all simulations, the timestep was 2 fs.

## B. Enthalpy excess route

To assess the predictions of Derjaguin’s theory, we calculated the thermo-osmotic coefficient using Eq. (14). To this end, we computed the enthalpy excess by following

the procedure described in Ref. 30, with a system composed of 1500 water molecules and 9 pairs of ions, confined between graphene walls of dimensions $L_x \sim 34$ Å and $L_y \sim 26$ Å, for a separation between the walls of approximately 50 Å. The enthalpy excess density is computed with Eq. (10), in which the enthalpy per particle is defined as$^{30}$:

$$h_i(z) = u_i(z) + \frac{p^{\parallel}(z)}{n_{\text{tot}}(z)}, \quad (23)$$

with $u_i$ the internal energy of species i, $n_{tot}$ the total number density, and $p^{\parallel} = p^{xx} = p^{yy}$ the components of the virial pressure tensor parallel to the surface. Indeed, the pressure is anisotropic near the surface, and one must consider the parallel components of the pressure tensor to compute the enthalpy density$^{20,23,24}$; see in particular a detailed discussion in the supplementary material of Ref. 23. To estimate the error on the enthalpy excess, we computed the statistical error based on independent measurements.

We also calculated the slip length for these systems. The slip length is defined as$^{52}$:

$$b = \frac{\eta}{\lambda}, \quad (24)$$

where $\eta$ is the viscosity of the bulk liquid, and $\lambda$ is the solid-liquid friction coefficient. One way to compute these quantities is to use Green-Kubo formulas, at equilibrium. We computed the viscosity in independent simulations of a bulk liquid, using the Green-Kubo relation$^{53}$:

$$\eta = \frac{V}{k_B T} \int_0^\infty \langle p_{\alpha\beta}(t) p_{\alpha\beta}(0) \rangle dt, \quad (25)$$

with V the volume of the simulation box, and $p_{\alpha\beta}$ the independent components of the traceless pressure tensor. Similarly, the friction coefficient is calculated at equilibrium using$^{52}$:

$$\lambda = \frac{1}{A k_B T} \int_0^\infty \langle F_\alpha(t) F_\alpha(0) \rangle dt, \quad (26)$$

where A is the wall surface and F is the force acting on the wall along the $\alpha$ direction, $\alpha = x, y$.

One can also compute the slip length using NEMD simulations and the Navier boundary condition$^{54,55}$:

$$b = \frac{v_s}{\dot{\gamma}}, \quad (27)$$

where $\dot{\gamma}$ is the bulk shear rate and $v_s$ is the slip velocity, defined as the difference between the wall velocity and the velocity of the fluid at the position of the hydrodynamic wall, given by $\dot{\gamma}h/2$, where h, the hydrodynamic height of the liquid, is given by$^{56}$:

$$h = \frac{M}{\rho_{\text{bulk}} A}, \quad (28)$$

with M the total mass of the fluid and $\rho_{bulk}$ the bulk mass density. To compute those quantities, we moved the
<figure><img src="image_3.png" /><figcaption>FIG. 2. Different types of graphitic surfaces considered: a neutral pristine graphene (PG) surface (a), a polarized PG surface (b), a neutral graphene oxide (GO) surface (c), and an electrically charged GO surface, for which the hydrogen from the hydroxyl groups are removed to create charged groups (d). The snapshots were produced with the software Ovito$^{51}$.</figcaption></figure>

The interaction parameters for these hydroxyl groups were based on the phenol parameters within the all-atom optimized potentials for liquid simulation (OPLS-AA) force field$^{49}$.

The graphene sheets were either kept neutral or charged. To charge the PG surface uniformly, we assigned the same charge to all carbon atoms $q = \Sigma S/N_{wall}$, where S is the surface area of the sheet, and $N_{wall}$ is the number of carbon atoms, resulting in a surface charge density $\Sigma$. To charge the GO, we remove the hydrogen atoms from the hydroxyl groups. In practice, this was achieved by setting the charge of the carbon, oxygen and hydrogen of the hydroxyl groups to $q_{C_{phenol}} = 0$ e, $q_{O_{phenol}} = -0.85$ e and $q_{H_{phenol}} = 0$ e, respectively, to match the rescaled charge of the ions. Counter-ions were added to the system to keep it electrically neutral.

In all simulations, we froze the bottom surface. We used the top surface as a piston during an equilibration phase that lasted 1 ns, before fixing it at its equilibrium position to set the pressure to 1 atm, for all simulations except the DSM (see details of DSM simulations in Section III C).

For the systems with a temperature gradient, we used a NVE time integration, and fixed the fluid temperature at 280 and 360 K for the cold and hot region, respectively, via a canonical sampling thermostat that uses global velocity rescaling with Hamiltonian dynamics$^{50}$, applied in the y and z directions, following a previous thermo-osmosis study$^{24}$, with a damping time of 200 fs. The mean temperature was 320 K. For other systems, we used a Nosé-Hoover thermostat to set the temperature to 320 K with a damping time of 200 fs. For all simulations, the timestep was 2 fs.

## B. Enthalpy excess route

To assess the predictions of Derjaguin’s theory, we calculated the thermo-osmotic coefficient using Eq. (14). To this end, we computed the enthalpy excess by following the procedure described in Ref. 30, with a system composed of 1500 water molecules and 9 pairs of ions, confined between graphene walls of dimensions $L_x \sim 34$ Å and $L_y \sim 26$ Å, for a separation between the walls of

approximately 50 Å. The enthalpy excess density is computed with Eq. (10), in which the enthalpy per particle is defined as$^{30}$:

$$h_i(z) = u_i(z) + \frac{p^{\parallel}(z)}{n_{\text{tot}}(z)}, \quad (23)$$

with $u_i$ the internal energy of species i, $n_{tot}$ the total number density, and $p^{\parallel} = p^{xx} = p^{yy}$ the components of the virial pressure tensor parallel to the surface. Indeed, the pressure is anisotropic near the surface, and one must consider the parallel components of the pressure tensor to compute the enthalpy density$^{20,24}$.

We also calculated the slip length for these systems. The slip length is defined as$^{52}$:

$$b = \frac{\eta}{\lambda}, \quad (24)$$

where $\eta$ is the viscosity of the bulk liquid, and $\lambda$ is the solid-liquid friction coefficient. One way to compute these quantities is to use Green-Kubo formulas, at equilibrium. The Green-Kubo relation for the viscosity is$^{53}$:

$$\eta = \frac{V}{k_B T} \int_0^{\infty} \langle p_{\alpha\beta}(t) p_{\alpha\beta}(0) \rangle dt, \quad (25)$$

with V the volume of the simulation box, and $p_{\alpha\beta}$ the independent components of the traceless pressure tensor. Similarly, the friction coefficient is calculated at equilibrium using$^{52}$:

$$\lambda = \frac{1}{A k_B T} \int_0^{\infty} \langle F_{\alpha}(t) F_{\alpha}(0) \rangle dt, \quad (26)$$

where A is the wall surface and F is the force acting on the wall along the $\alpha$ direction, $\alpha = x, y$.

One can also compute the slip length using NEMD simulations and the Navier boundary condition$^{54,55}$:

$$b = \frac{v_s}{\dot{\gamma}}, \quad (27)$$

where $\dot{\gamma}$ is the bulk shear rate and $v_s$ is the slip velocity, defined as the difference between the wall velocity and the velocity of the fluid at the position of the hydrodynamic wall, given by $\dot{\gamma}h/2$, where h, the hydrodynamic height of the liquid, is given by$^{56}$:

$$h = \frac{M}{\rho_{\text{bulk}} A}, \quad (28)$$

with M the total mass of the fluid and $\rho_{bulk}$ the bulk mass density. To compute those quantities, we moved the walls in opposite parallel directions, at a constant speed $V_x \in [10, 50]$ m/s, generating a linear velocity profile far from the wall. We verified that the resulting quantities remained in the linear response regime.
walls in opposite parallel directions, at a constant speed $V_x \in [10, 50]$ m/s, generating a linear velocity profile far from the wall. We verified that the resulting quantities remained in the linear response regime. Green-Kubo measurements were preferred for systems with high slip lengths, while non-equilibrium simulations were used for smaller slip lengths.

### C. Entropy excess route (DSM)

The entropy excess can be derived from the solid-liquid entropy per unit area $s_{sl}$, which can be computed using the DSM. It is based on the thermodynamic integration of the solid-liquid work of adhesion $W_{sl}$ along a coupling parameter $\kappa$, that quasi-statically weakens the solid-liquid LJ interactions to a point $\kappa_0$ where the liquid is no longer influenced by the solid. Here, the system is only composed of one graphitic surface and an aqueous electrolyte. The entropy difference is computed using the following relation³¹:

$$W_{sl}(\kappa) = \delta u_{sl}(\kappa) - T\delta s_{sl}(\kappa). \quad (29)$$

Here, $\delta A(\kappa) = A(\kappa) - A(\kappa_0)$ is the difference of the quantity $A$ between its value at $\kappa$ and a reference state, defined by $\kappa_0$.

The work of adhesion can be expressed in terms of Helmholtz free energy:

$$W_{sl}(\kappa) = \left(\frac{\Delta F}{A}\right)_{N,V,T} \quad (30)$$

$$= -\frac{1}{A}[F(\kappa) - F(\kappa_0)] \quad (31)$$

$$= -\frac{1}{A} \int_{\kappa_0}^{\kappa} \left\langle \frac{\partial \phi(\kappa')}{\partial \kappa'} \right\rangle_{N,V,T} d\kappa', \quad (32)$$

with $\phi$ the solid-liquid potential. The surface is not charged here, thus we only have LJ solid-liquid interactions:

$$\phi_{sl(LJ)}(\kappa) = \kappa \sum_{i \in \text{liquid}} \sum_{j \in \text{wall}} 4\epsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{6} \right], \quad (33)$$

its derivative is given by

$$\frac{\partial \phi_{sl(LJ)}(\kappa)}{\partial \kappa} = \phi_{sl(LJ)}(\kappa = 1). \quad (34)$$

In practice we assume ergodicity and we substitute ensemble averages by time averages:

$$\left\langle \frac{\partial \phi(\kappa)}{\partial \kappa} \right\rangle = \frac{\overline{\partial \phi(\kappa)}}{\partial \kappa} = \frac{1}{t_{\text{sim}}} \int_{0}^{t_{\text{sim}}} \frac{\partial \phi(\kappa)}{\partial \kappa} dt. \quad (35)$$

The integration was performed on 30 points, $\kappa \in [0, 1.5]$, to cover all the wettings used in NEMD. For each point, the equilibration phase lasted 0.6 ns and we saved the

<figure><img src="image_3.png" /><figcaption>FIG. 3. System used for the NEMD computation of the thermo-osmotic coefficient. A temperature gradient generates a thermo-osmotic flow in region 1, it is converted to a Poiseuille flow in region 2, because the surface in this region is thermo-osmosis-insensitive, i.e. the enthalpy excess is null for this surface.</figcaption></figure>

trajectory of the simulation for 1 ns, with a period of 1 ps. We then applied Eq. (32) on those trajectory files, by setting $\kappa$ to 1, with the “rerun” command of LAMMPS.

We also computed $\delta u_{sl}(\kappa)$ with the rerun command, by turning off all interactions parameters except for the solid-liquid ones, and computing the total energy per unit area. We thus computed $\delta s_{sl}$ using Eq. (29). Finally,

$$s_{sl}(\kappa) = \Delta S = \delta s_{sl}(\kappa) + s_{sl}(\kappa_0), \quad (36)$$

in which we have considered the solid-liquid entropy to be equal to the entropy excess. We discuss this hypothesis in Appendix A. The solid-liquid entropy of the reference state $s_{sl}(\kappa_0)$ is unknown. We choose the value of $s_{sl}(\kappa_0)$ so that the solid-liquid entropy $s_{sl}(\kappa)$ vanishes for $\epsilon_{CO} = 0.105$ kcal/mol, to have a vanishing entropy excess when $\Delta H = 0$. To estimate the error on the entropy excess, we computed the statistical error based on block averaging of the measured values during a single simulation.

### D. Non-equilibrium computation of the thermo-osmotic coefficient

Another mean of computing the thermo-osmotic coefficient is to use Eq. (9), i.e. generating a thermo-osmotic flow by applying a temperature gradient along the channel in the $x$ direction. Dealing with periodic boundary conditions in the $x$ direction, one can only apply the desired temperature gradient locally, not in the entire system. Here we divide the system in 2 regions: in region 1, we apply a temperature gradient $\nabla T_1$, and in region 2, we apply $\nabla T_2$, such that $L_1 \nabla T_1 = -L_2 \nabla T_2 = \Delta T$ (Fig. 3). Let us first consider an electrically neutral surface: the average velocity is given by the sum of a thermo-osmotic and a Poiseuille flow. It is uniform along the channel:

$$\overline{v}_{\text{tot}} = \overline{v}_{\text{to},1} + \overline{v}_{\text{P},1} = \overline{v}_{\text{to},2} + \overline{v}_{\text{P},2} \quad (37)$$
### C. Entropy excess route (DSM)

The entropy excess can be derived from the solid-liquid entropy per unit area $s_{sl}$, which can be computed using the DSM. It is based on the thermodynamic integration of the solid-liquid work of adhesion $W_{sl}$ along a coupling parameter $\kappa$, that slowly weakens the solid-liquid LJ interactions to a point $\kappa_0$ where the liquid is no longer influenced by the solid. Here, the system is only composed of one graphitic surface and an aqueous electrolyte. The entropy difference is computed using the following relation³¹:

$$W_{sl}(\kappa) = \delta u_{sl}(\kappa) - T \delta s_{sl}(\kappa). \quad (29)$$

Here, $\delta A(\kappa) = A(\kappa) - A(\kappa_0)$ is the difference of the quantity $A$ between its value at $\kappa$ and a reference state, defined by $\kappa_0$.

The work of adhesion can be expressed in terms of Helmholtz free energy:

$$W_{sl}(\kappa) = \left(\frac{\Delta F}{A}\right)_{N,V,T} \quad (30)$$

$$= -\frac{1}{A}[F(\kappa) - F(\kappa_0)] \quad (31)$$

$$= -\frac{1}{A} \int_{\kappa_0}^{\kappa} \left\langle \frac{\partial \phi(\kappa')}{\partial \kappa'} \right\rangle_{N,V,T} d\kappa', \quad (32)$$

with $\phi$ the solid-liquid potential. The surface is not charged here, thus we only have LJ solid-liquid interactions:

$$\phi_{sl(LJ)}(\kappa) = \kappa \sum_{i \in \text{liquid}} \sum_{j \in \text{wall}} 4\epsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{6} \right], \quad (33)$$

its derivative is given by

$$\frac{\partial \phi_{sl(LJ)}(\kappa)}{\partial \kappa} = \phi_{sl(LJ)}(\kappa = 1). \quad (34)$$

In practice we assume ergodicity and we substitute ensemble averages by time averages:

$$\left\langle \frac{\partial \phi(\kappa)}{\partial \kappa} \right\rangle = \frac{\overline{\partial \phi(\kappa)}}{\partial \kappa} = \frac{1}{t_{\text{sim}}} \int_{0}^{t_{\text{sim}}} \frac{\partial \phi(\kappa)}{\partial \kappa} dt. \quad (35)$$

The integration was performed on 30 points, $\kappa \in [0, 1.5]$, to cover all the wettings used in NEMD. For each point, the equilibration phase lasted 0.6 ns and we saved the trajectory of the simulation for 1 ns, with a period of 1 ps. We then applied Eq. (32) on those trajectory files, by setting $\kappa$ to 1, with the “rerun” command of LAMMPS.

We also computed $\delta u_{sl}(\kappa)$ with the rerun command, by turning off all interactions parameters except for the solid-liquid ones, and computing the total energy per unit area. We thus computed $\delta s_{sl}$ using Eq. (29). Finally,

$$s_{sl}(\kappa) = \Delta S = \delta s_{sl}(\kappa) + s_{sl}(\kappa_0), \quad (36)$$

<figure><img src="image_4.png" /><figcaption>FIG. 3. System used for the NEMD computation of the thermo-osmotic coefficient. A temperature gradient generates a thermo-osmotic flow in region 1, it is converted to a Poiseuille flow in region 2, because the surface in this region is thermo-osmosis-insensitive, i.e. the enthalpy excess is null for this surface.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4. Enthalpy excess as a function of the carbon-oxygen interaction parameters $\epsilon_{CO}$ for a PG surface (blue line). The curve is fitted linearly (red curve) to determine the point at which $\Delta H = 0$, which corresponds to $\epsilon_{CO} = 0.105$ kcal/mol.</figcaption></figure>

in which we have considered the solid-liquid entropy to be equal to the entropy excess. We discuss this hypothesis in Appendix A. The solid-liquid entropy of the reference state $s_{sl}(\kappa_0)$ is unknown. We choose the value of $s_{sl}(\kappa_0)$ so that the solid-liquid entropy $s_{sl}(\kappa)$ vanishes for $\epsilon_{CO} = 0.105$ kcal/mol, to have a vanishing entropy excess when $\Delta H = 0$.

### D. Non-equilibrium computation of the thermo-osmotic coefficient

Another mean of computing the thermo-osmotic coefficient is to use Eq. (9), i.e. generating a thermo-osmotic flow by applying a temperature gradient along the channel in the $x$ direction. Dealing with periodic boundary conditions in the $x$ direction, one can only apply the de-
<figure><img src="image_4.png" /><figcaption>FIG. 4. Enthalpy excess as a function of the carbon-oxygen interaction parameters $\varepsilon_{CO}$ for a PG surface (blue line). The curve is fitted linearly (red curve) to determine the point at which $\Delta H = 0$, which corresponds to $\varepsilon_{CO} = 0.105$ kcal/mol.</figcaption></figure>

where $\overline{v_{to,i}}$ is the average thermo-osmotic velocity over the channel thickness, and $\overline{v_{P,i}}$ is the average Poiseuille flow, in region $i$. The average Poiseuille velocity is given by¹²:

$$ \overline{v_{P,i}} = -\frac{\nabla p_i h^2}{12\eta} \left(1 + 6\frac{b_i}{h}\right), \quad (38) $$

where $h$ is the height of the channel, $b_i$ is the slip length, and $\nabla p_i$ is the pressure gradient in region $i$. $\nabla p_1 = \Delta p/L_1$ and $\nabla p_2 = -\Delta p/L_2$, with $\Delta p$ the pressure difference. For an electrically neutral surface, the thermo-osmotic flow rate is given by:

$$ \overline{v_{to,i}} = \frac{1}{h} \int_{0}^{h} v_{to,i}(z) dz. \quad (39) $$

When the channel is large enough, i.e. $h \gg \lambda_h$, $v_{to}(z) \approx v_{to}(\infty)$, the average thermo-osmotic velocity becomes:

$$ \overline{v_{to,i}} = v_{to,\nabla T}^{i} = -M_{to,i} \frac{\nabla T_{i}}{T}, \quad (40) $$

where $\nabla T_i$ is the temperature gradient in region $i$; $\nabla T_1 = \Delta T/L_1$ and $\nabla T_2 = -\Delta T/L_2$, with $\Delta T$ the temperature difference. When the surfaces in region 1 and 2 are identical, conservation of flow rate predicts that Poiseuille and thermo-osmotic flows perfectly cancel each other, so that total flow is zero. One solution to this problem is to define the surface in region 2 such that $\Delta H = 0$. We will call such a surface thermo-osmosis-insensitive (TOI).

To obtain a TOI surface for the neutral PG wall, we start from a graphene surface and we slightly modify its solid-liquid interaction parameters, $\varepsilon_{CO} \in [0.08, 0.15]$ kcal/mol in our case, which allows us to realize a surface with a null enthalpy excess for $\varepsilon_{CO} = 0.105$ kcal/mol (Fig. 4). Note that changing the $\varepsilon_{CO}$ amounts to changing the contact angle (see Appendix B for details on the

calculation of the contact angle). We applied the same procedure for the polarized PG, neutral and charged GO surfaces: for each value of the surface charge or oxidation rate, we varied $\varepsilon_{CO}$ to find the value canceling $\Delta H$.

Considering a TOI surface, the equation of the average flow velocity becomes:

$$ \overline{v_{tot}} = -M_{to} \frac{\nabla T_1}{T} + \overline{v_{P,1}} = \overline{v_{P,2}}, \quad (41) $$

where we set $M_{to,1} = M_{to}$. Using the expression of $\overline{v_{P,i}}$, Eq. (38), where $\nabla p_1 = \Delta p/L_1$ and $\nabla p_2 = -\Delta p/L_2$, one can obtain the expression of $\Delta p$:

$$ \Delta p = \overline{v_{tot}} \times \frac{12\eta L_2}{h^2 (1 + 6 b_2/h)}, \quad (42) $$

and then relate the thermo-osmotic coefficient to the total flow velocity:

$$ M_{to} = -\frac{\overline{v_{tot}}}{\nabla T_1/T} \left[ 1 + \frac{L_2 (1 + 6 b_1/h)}{L_1 (1 + 6 b_2/h)} \right]. \quad (43) $$

When considering an electrically charged surface, one must take into account the permittivity-osmotic flow generated in region 2. The total flow rate becomes:

$$ \overline{v_{tot}} = -M_{to} \frac{\nabla T_1}{T} + \overline{v_{P,1}} = \overline{v_{to,\nabla\varepsilon,2}} + \overline{v_{P,2}}, \quad (44) $$

with

$$ \overline{v_{to,\nabla\varepsilon,2}} = v_{to,\nabla\varepsilon,2}^{\text{PB}}, \quad (45) $$

if we consider the channel large enough compared to the EDL, $h \gg \min(\lambda_D, \ell_{GC})$. The thermo-osmotic coefficient is then given by:

$$ M_{to} = -\frac{1}{\nabla T_1/T} \left\{ \left( \overline{v_{tot}} - v_{to,\nabla\varepsilon,2}^{\text{PB}} \right) \left[ 1 + \frac{L_2 (1 + 6 b_1/h)}{L_1 (1 + 6 b_2/h)} \right] + v_{to,\nabla\varepsilon,2}^{\text{PB}} \right\}. \quad (46) $$

We fixed $L_1 \sim 100$, $L_2 \sim 20$ Å, and $L_y \sim 26$ Å. The size of the thermostated regions was set to 10 Å. We verified that the chosen temperature difference, $\Delta T = 80$ K, lied within the linear response regime. For NEMD simulations, we performed only one measure per data point; the uncertainty originates in the uncertainty on the enthalpy excess of the thermo-osmosis-insensitive surface.

# IV. RESULTS AND DISCUSSION

## A. Comparison of the different approaches

We begin by comparing the different methods used to compute the thermo-osmotic coefficient and we assess their agreement. Figure 5 presents the thermo-osmotic
sired temperature gradient locally, not in the entire system. Here we divide the system in 2 regions: in region 1, we apply a temperature gradient $\nabla T_1$, and in region 2, we apply $\nabla T_2$, such that $L_1\nabla T_1 = -L_2\nabla T_2 = \Delta T$ (Fig. 3). Let us first consider an electrically neutral surface: the average velocity is given by the sum of a thermo-osmotic and a Poiseuille flow. It is uniform along the channel:

$$ \overline{v}_{\text{tot}} = \overline{v}_{\text{to,1}} + \overline{v}_{\text{P,1}} = \overline{v}_{\text{to,2}} + \overline{v}_{\text{P,2}} \quad (37) $$

where $\overline{v}_{\text{to,i}}$ is the average thermo-osmotic velocity over the channel thickness, and $\overline{v}_{\text{P,i}}$ is the average Poiseuille flow, in region i. The average Poiseuille velocity is given by¹²:

$$ \overline{v}_{\text{P,i}} = -\frac{\nabla p_i h^2}{12\eta} \left(1 + 6\frac{b_i}{h}\right), \quad (38) $$

where h is the height of the channel, $b_i$ is the slip length, and $\nabla p_i$ is the pressure gradient in region i. $\nabla p_1 = \Delta p/L_1$ and $\nabla p_2 = -\Delta p/L_2$, with $\Delta p$ the pressure difference. For an electrically neutral surface, the thermo-osmotic flow rate is given by:

$$ \overline{v}_{\text{to,i}} = \frac{1}{h} \int_{0}^{h} v_{\text{to,i}}(z) dz. \quad (39) $$

When the channel is large enough, i.e. $h \gg \lambda_h$, $v_{\text{to}}(z) \approx v_{\text{to}}(\infty)$, the average thermo-osmotic velocity becomes:

$$ \overline{v}_{\text{to,i}} = v_{\text{to},\nabla T}^{i} = -M_{\text{to},i} \frac{\nabla T_{i}}{T}, \quad (40) $$

where $\nabla T_i$ is the temperature gradient in region i; $\nabla T_1 = \Delta T/L_1$ and $\nabla T_2 = -\Delta T/L_2$, with $\Delta T$ the temperature difference. When the surfaces in region 1 and 2 are identical, conservation of flow rate predicts that Poiseuille and thermo-osmotic flows perfectly cancel each other, so that total flow is zero. One solution to this problem is to define the surface in region 2 such that $\Delta H = 0$. We will call such a surface thermo-osmosis-insensitive (TOI).

To obtain a TOI surface for the neutral PG wall, we start from a graphene surface and we slightly modify its solid-liquid interaction parameters, $\epsilon_{CO} \in [0.08, 0.15]$ kcal/mol in our case, which allows us to realize a surface with a null enthalpy excess for $\epsilon_{CO} = 0.105$ kcal/mol (Fig. 4). Note that changing the $\epsilon_{CO}$ amounts to changing the contact angle (see Appendix B for details on the calculation of the contact angle). We applied the same procedure for the polarized PG, neutral and charged GO surfaces: for each value of the surface charge or oxidation rate, we varied $\epsilon_{CO}$ to find the value canceling $\Delta H$.

Considering a TOI surface, the equation of the average flow velocity becomes:

$$ \overline{v}_{\text{tot}} = -M_{\text{to}} \frac{\nabla T_1}{T} + \overline{v}_{\text{P,1}} = \overline{v}_{\text{P,2}}, \quad (41) $$

where we set $M_{\text{to,1}} = M_{\text{to}}$. The thermo-osmotic coefficient is computed by correcting the total flow with the Poiseuille flow:

$$ M_{\text{to}} = -\frac{\overline{v}_{\text{tot}}}{\nabla T_1/T} \left[ 1 + \frac{L_2 (1 + 6 b_1 / h)}{L_1 (1 + 6 b_2 / h)} \right]. \quad (42) $$

When considering an electrically charged surface, one must take into account the permittivity-osmotic flow generated in region 2. The total flow rate becomes:

$$ \overline{v}_{\text{tot}} = -M_{\text{to}} \frac{\nabla T_1}{T} + \overline{v}_{\text{P,1}} = \overline{v}_{\text{to},\nabla \epsilon,2} + \overline{v}_{\text{P,2}}, \quad (43) $$

with

$$ \overline{v}_{\text{to},\nabla \epsilon,2} = v_{\text{to},\nabla \epsilon,2}^{\text{PB}}, \quad (44) $$

if we consider the channel large enough compared to the EDL, $h \gg \min(\lambda_D, \ell_{GC})$. The thermo-osmotic coefficient is then given by:

$$ M_{\text{to}} = -\frac{1}{\nabla T_1/T} \left\{ \left( \overline{v}_{\text{tot}} - v_{\text{to},\nabla\epsilon,2}^{\text{PB}} \right) \left[ 1 + \frac{L_2 (1 + 6 b_1 / h)}{L_1 (1 + 6 b_2 / h)} \right] + v_{\text{to},\nabla\epsilon,2}^{\text{PB}} \right\}. \quad (45) $$

We fixed $L_1 \sim 100$, $L_2 \sim 20$ Å, and $L_y \sim 26$ Å. The size of the thermostated regions was set to 10 Å. We verified that the chosen temperature difference, $\Delta T = 80$ K, lied within the linear response regime.

# IV. RESULTS AND DISCUSSION

## A. Comparison of the different approaches

We begin by comparing the different methods used to compute the thermo-osmotic coefficient and we assess their agreement. Figure 5 presents the thermo-osmotic coefficient $M_{\text{to}}$, on a graphene-like surface with tunable wetting properties, as a function of the contact angle. We found a remarkable agreement between Derjaguin's approach (both through the enthalpy excess route and the entropy excess route) and the direct computation of the thermo-osmotic coefficient using NEMD. The calculations based on the enthalpy or the entropy excess give additional information on the contributions to $M_{\text{to}}$, i.e. the driving force $\nabla\gamma = -\Delta H\nabla T/T = -\Delta S\nabla T$, and hydrodynamic properties (the slip length and viscosity). The agreement between the calculation of entropy excess with the DSM and the computation of enthalpy excess, see Fig. 5.c, tends to validate the computation of the enthalpy excess using the pressure tensor. The DSM is much more accurate than the direct enthalpy excess calculation, while being faster in terms of computational time. A first limitation of the DSM to compute the entropy excess is the presence of a parameter in Eq. (36), the solid-liquid entropy of the reference state $\kappa_0$. We determined this value by finding the point at which the entropy excess is zero, for which $M_{\text{to}}$ vanishes (based on the excess enthalpy route results). Another limitation of this method is that it cannot be used with electrically charged surface. Indeed, it requires turning off the solid-liquid Coulombic interactions, which in practice means
<figure><img src="image_5.png" /></figure>

FIG. 5. Thermo-osmotic coefficient of a pristine graphene-like surface with varying wetting properties, as a function of the contact angle. In figure a), the blue squares represent the values of the NEMD simulations, wherein a flow is generated by applying a thermal gradient along the interface (Eq. 43). The orange circles represent the curve generated by the computation of the enthalpy excess and the slip length, following Derjaguin's approach (Eq. 14). The green triangles represent the results of the computation of the entropy excess through the dry-surface method (DSM) (Eq. 20). Panel b) represents the slip length. Panel c) represents the enthalpy excess (orange circles) and the entropy excess (green triangles). Some error bars do not appear because they are smaller than the symbols.

coefficient $M_{\text{to}}$, on a graphene-like surface with tunable wetting properties, as a function of the contact angle. We found a remarkable agreement between Derjaguin's approach (both through the enthalpy excess route and the entropy excess route) and the direct computation of the thermo-osmotic coefficient using NEMD. The calculations based on the enthalpy or the entropy excess give additional information on the contributions to $M_{\text{to}}$, i.e. the driving force $\nabla\gamma = -\Delta H\nabla T/T = -\Delta S\nabla T$, and hydrodynamic properties (the slip length and viscosity). The agreement between the calculation of entropy excess with the DSM and the computation of enthalpy excess, see Fig. 5.c, tends to validate the computation of the enthalpy excess using the pressure tensor. The DSM is much more accurate than the direct enthalpy excess calculation, while being faster in terms of computational

<figure><img src="image_6.png" /></figure>

FIG. 6. Thermo-osmotic coefficient of a polarized graphene surface as a function of the surface charge density. The non-equilibrium computation (blue squares) agrees well with Derjaguin's formula based on the enthalpy excess (orange circles). The Poisson-Boltzmann (PB) prediction to thermo-osmosis is also represented (green line). As discussed in the text, it is negligible compared to the contribution of water enthalpy excess. Panels b) and c) represent respectively the slip length and the enthalpy excess. Some error bars do not appear because they are smaller than the symbols.

time. A first limitation of the DSM to compute the entropy excess is the presence of a parameter in Eq. (36), the solid-liquid entropy of the reference state $\kappa_0$. We determined this value by finding the point at which the entropy excess is zero, for which $M_{\text{to}}$ vanishes (based on the excess enthalpy route results). Another limitation of this method is that it cannot be used with electrically charged surface. Indeed, it requires turning off the solid-liquid Coulombic interactions, which in practice means removing the charge of the wall, leading to non-neutral systems. In contrast, both NEMD and $\Delta H$ route can be used on charged surfaces, and provide consistent results (Fig. 6).

Returning to the NEMD computation, the agreement between this method and the methods based on the computation of the enthalpy and the entropy excess tends to confirm the validity of Derjaguin's approach to compute the thermo-osmotic coefficient for smooth surfaces. In particular, Eq. (1) assumes that the viscosity is homogeneous
<figure><img src="image_6.png" /></figure>

FIG. 5. Thermo-osmotic coefficient of a pristine graphene-like surface with varying wetting properties, as a function of the contact angle. In figure a), the blue squares represent the values of the NEMD simulations, wherein a flow is generated by applying a thermal gradient along the interface (Eq. 42). The orange circles represent the curve generated by the computation of the enthalpy excess and the slip length, following Derjaguin's approach (Eq. 14). The green triangles represent the results of the computation of the entropy excess through the dry-surface method (DSM) (Eq. 20). Panel b) represents the slip length. Panel c) represents the enthalpy excess (orange circles) and the entropy excess (green triangles).

removing the charge of the wall, leading to non-neutral systems. In contrast, both NEMD and $\Delta H$ route can be used on charged surfaces, and provide consistent results (Fig. 6).

Returning to the NEMD computation, the agreement between this method and the methods based on the computation of the enthalpy and the entropy excess tends to confirm the validity of Derjaguin's approach to compute the thermo-osmotic coefficient for smooth surfaces. NEMD is comparable to enthalpy excess computation in terms of computational time while being much more precise. On the other hand, it requires finding a TOI surface, i.e., computing the enthalpy excess for a range of interaction parameters and finding the surface that nullifies the enthalpy excess. Another limitation of the method is that it provides limited microscopic information, as it only returns one quantity, $M_{to}$.

<figure><img src="image_7.png" /></figure>

FIG. 6. Thermo-osmotic coefficient of a polarized graphene surface as a function of the surface charge density. The non-equilibrium computation (blue squares) agrees well with Derjaguin's formula based on the enthalpy excess (orange circles). The Poisson-Boltzmann (PB) prediction to thermo-osmosis is also represented (green line). As discussed in the text, it is negligible compared to the contribution of water enthalpy excess. Panels b) and c) represent respectively the slip length and the enthalpy excess.

The agreement between the methods remains consistent for the neutral GO surface (Fig. 7). However, for the charged GO surface, there is a clear departure from Derjaguin's theory (Fig. 8). Two reasons could explain this difference; the first reason is the assumption that viscosity is constant along the channel. While this assumption is reasonable for smooth surfaces, it may not be true for defective ones. In Ref. 17, the authors computed the velocity profile of water near diamond interfaces by generating a Couette flow in the channel. They considered two types of surfaces: a smooth diamond surface and a defective surface with hydroxyl groups as defects. They found that the presence of hydroxyl groups could significantly increase the viscosity near the surface. Consequently, the thermo-osmotic velocity predicted by Derjaguin's formula could be overestimated$^{30}$.

Another possible reason involves the distribution of enthalpy excess. In Fig. 9, we computed the map of enthalpy excess for both the pristine and defective surfaces. While the enthalpy excess of the smooth surface is ho-
<figure><img src="image_7.png" /></figure>

FIG. 7. Thermo-osmotic coefficient of a neutral graphene oxide surface as a function of the oxidation rate. We used two methods to compute the thermo-osmotic coefficient, non-equilibrium simulations (blue squares) and the enthalpy excess route (orange circles). The inset shows the slip length as a function of the oxidation rate (blue circles), the red curve represents the prediction of Eq. (48) with $\sigma_h = 1.2$ Å. The thickness of the interaction layer $\lambda_h$ is shown with a dotted line for comparison. Some error bars do not appear because they are smaller than the symbols.

<figure><img src="image_8.png" /></figure>

FIG. 8. Thermo-osmotic coefficient as a function of the surface charge density for a charged graphene oxide surface. The blue squares represent the non-equilibrium calculation and the orange circles show the enthalpy excess route. Inset: slip length as a function of the surface charge density in blue circles, theoretical prediction in red (using Eq. (47) with $\sigma_h = 8.8$ Å). The thickness of the interaction layer $\lambda_h$ is shown with a dotted line for comparison. Some error bars do not appear because they are smaller than the symbols.

neous, and deviations from this assumption should have a strong impact on the osmotic response, see Appendix C. Therefore, our results tend to show the viscosity does not vary significantly close to the walls considered here. This is compatible with previous work showing a constant viscosity close to hydrophobic walls¹⁸. NEMD is

<figure><img src="image_9.png" /></figure>

FIG. 9. Enthalpy excess map for a pristine graphene surface (a) and for a graphene oxide surface (b), the black points represent the position of the hydroxyl groups. Units of the colorbar are in J/m².

comparable to enthalpy excess computation in terms of computational time while being much more precise. On the other hand, it requires finding a TOI surface, i.e., computing the enthalpy excess for a range of interaction parameters and finding the surface that nullifies the enthalpy excess. Another limitation of the method is that it provides limited microscopic information, as it only returns one quantity, $M_{to}$.

The agreement between the methods remains consistent for the neutral GO surface (Fig. 7). However, for the charged GO surface, there is a clear departure from Derjaguin's theory (Fig. 8). Two reasons could explain this difference; the first reason is the assumption that viscosity is constant along the channel. While this assumption is reasonable for smooth surfaces, it may not be true for defective ones. In Ref. 17, the authors computed the velocity profile of water near diamond interfaces by generating a Couette flow in the channel. They considered two types of surfaces: a smooth diamond surface and a defective surface with hydroxyl groups as defects. They found that the presence of hydroxyl groups could significantly increase the viscosity near the surface. Consequently, the thermo-osmotic velocity predicted by Derjaguin's formula could be overestimated, see Appendix C and Ref. 30.

Another possible reason involves the distribution of enthalpy excess. In Fig. 9, we computed the map of enthalpy excess for both the pristine and defective surfaces using a NEMD simulation with a flow along the *x* direction. While the enthalpy excess of the smooth surface is homogeneous on the *xy* plane, it becomes heterogeneous with the GO surface. Specifically, the enthalpy excess is not constant along the *y* direction. The non-uniform enthalpy excess near the surface could explain the deviation from Derjaguin's theory, which assumes a uniform enthalpy excess.
<figure><img src="image_8.png" /><figcaption>FIG. 7. Thermo-osmotic coefficient of a neutral graphene oxide surface as a function of the oxidation rate. We used two methods to compute the thermo-osmotic coefficient, non-equilibrium simulations (blue squares) and the enthalpy excess route (orange circles). The inset shows the slip length as a function of the oxidation rate (blue circles), the red curve represents the prediction of Eq. (47) with $\sigma_h = 1.2$ Å. The thickness of the interaction layer $\lambda_h$ is shown with a dotted line for comparison.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8. Thermo-osmotic coefficient as a function of the surface charge density for a charged graphene oxide surface. The blue squares represent the non-equilibrium calculation and the orange circles show the enthalpy excess route. Inset: slip length as a function of the surface charge density in blue circles, theoretical prediction in red (using Eq. (46) with $\sigma_h = 8.8$ Å). The thickness of the interaction layer $\lambda_h$ is shown with a dotted line for comparison.</figcaption></figure>

mogeneous on the $xy$ plane, it becomes heterogeneous with the GO surface. Specifically, the enthalpy excess is not constant along the $y$ direction. The non-uniform enthalpy excess near the surface could explain the deviation from Derjaguin's theory, which assumes a uniform enthalpy excess.

<figure><img src="image_10.png" /><figcaption>FIG. 9. Enthalpy excess map for a pristine graphene surface (a) and for a graphene oxide surface (b), the black points represent the position of the hydroxyl groups. Units of the colorbar are in J/m$^2$.</figcaption></figure>

## B. Thermo-osmosis on homogeneous and heterogeneous surfaces

Now, let us examine the physical interpretation of the results on thermo-osmosis on graphene-like surfaces. In Fig. 5, one can observe that the thermo-osmotic coefficient varies considerably with the different wettability conditions considered. For hydrophilic surfaces, the flow is thermophilic, i.e. the liquid flows toward the hot side. However, for hydrophobic surfaces, the flow is thermophobic. The change of sign in the thermo-osmotic coefficient when modifying the wetting has been observed experimentally$^{35}$ and predicted with MD simulations$^{26,30}$ in other systems. The change of direction of the thermo-osmotic flow can be understood by the change of sign of the enthalpy excess (Fig. 5.c).

The thermo-osmotic coefficient computed here is one order of magnitude greater than that of previous MD works: it was found to be around $10^{-8}$ m$^2$/s for model LJ systems$^{26}$ and water-silica nanochannel$^{28}$. Although the enthalpy excess values are similar across the systems, the exceptionally high slip length of water on graphene surfaces leads to a massive amplification of the thermo-osmotic flow.

In Fig. 6, one can see a significant impact of the surface charge density on the thermo-osmotic coefficient, with a fluid flow enhancement by a factor of $\sim 3$ for the highest value of $M_{to}$ compared to the electrically neutral graphene surface. Once again, we observe a change of sign when changing the surface charge density, and this behavior can be fully explained by the variation of the enthalpy excess with $\Sigma$. The variation of the enthalpy excess has already been characterized in Ref. 30 for a model system, it is parabolic-like. For the pristine graphene considered here, the enthalpy excess is negative for $\Sigma = 0$ mC/m$^2$, which gives two points where $\Delta H = 0$, and thus two points for which the thermo-osmotic flow changes its direction. Finally, as stated in Sec. II, the
## B. Thermo-osmosis on homogeneous and heterogeneous surfaces

Now, let us examine the physical interpretation of the results on thermo-osmosis on graphene-like surfaces. In Fig. 5, one can observe that the thermo-osmotic coefficient varies considerably with the different wettability conditions considered. For hydrophilic surfaces, the flow is thermophilic, i.e. the liquid flows toward the hot side. However, for hydrophobic surfaces, the flow is thermophilic. The change of sign in the thermo-osmotic coefficient when modifying the wetting has been observed experimentally³⁵ and predicted with MD simulations²⁶,³⁰ in other systems. The change of direction of the thermo-osmotic flow can be understood by the change of sign of the enthalpy excess (Fig. 5.c).

The thermo-osmotic coefficient computed here is one order of magnitude greater than that of previous MD works: it was found to be around $10^{-8}$ m²/s for model LJ systems²⁶ and water-silica nanochannel²⁸. Although the enthalpy excess values are similar across the systems, the exceptionally high slip length of water on graphene surfaces leads to a massive amplification of the thermo-osmotic flow.

In Fig. 6, one can see a significant impact of the surface charge density on the thermo-osmotic coefficient, with a fluid flow enhancement by a factor of ~3 for the highest value of $M_{to}$ compared to the electrically neutral graphene surface. Once again, we observe a change of sign when changing the surface charge density, and this behavior can be fully explained by the variation of the enthalpy excess with $\Sigma$. The variation of the enthalpy excess has already been characterized in Ref. 30 for a model system, it is parabolic-like. For the pristine graphene considered here, the enthalpy excess is negative for $\Sigma = 0$ mC/m², which gives two points where $\Delta H = 0$, and thus two points for which the thermo-osmotic flow changes its direction. Finally, as stated in Sec. II, the Poisson-Boltzmann (PB) contribution to the enthalpy excess, and thus the thermo-osmotic coefficient, is negligible compared to the contribution of water molecules.

Figures 7 and 8 represent the thermo-osmotic coefficient for a neutral and a charged GO surface, respectively, along with their slip length. When considering defective surfaces, the presence of defects greatly reduces the solid-liquid slip⁵⁷, thereby diminishing the osmotic flow generated. The neutral GO surface still exhibits a relatively high thermo-osmotic coefficient for very small oxidation rates. The attenuation is more pronounced for charged defects, as the electrostatic interactions increase solid-liquid friction. In this case, the fluid flow quickly drops to values that are too small to be measured with our method.

The variation of the slip length with the surface density is given in Ref. 55 for a heterogeneous surface:

$$b = \frac{b_0}{1 + 3\pi\sigma_h b_0 |\Sigma|/q}, \quad (47)$$

with $b_0$ the slip length of the neutral surface, $\sigma_h$ the hydrodynamic radius of the defect plus counter-ion assembly, treated here as a fitting parameter. We used $\sigma_h = 8.8$ Å for our purpose. Similarly, for a neutral defective surface, one can write the slip length as:

$$b = \frac{b_0}{1 + 3\pi\sigma_h b_0 n_C \xi}, \quad (48)$$

with $n_C$ the number density of the carbon atoms, and $\sigma_h = 1.2$ Å the hydrodynamic radius of the hydroxyl groups. In both cases, the slip length quickly decreases even at low oxidation degree. Therefore, one can expect the experimental measured thermo-osmotic coefficient of a graphene nanochannel to be lower than what we have found numerically, due to the presence of defects.

## V. CONCLUSION

In this paper, we have presented three methods to compute the thermo-osmotic coefficient using molecular simulations. The first two are quite similar in that they are based on Derjaguin's theoretical framework. The first one requires to compute the enthalpy excess at equilibrium. The second one uses the dry-surface method to compute the solid-liquid entropy, which is linked to the enthalpy excess. The last method is based on NEMD, the thermo-osmotic velocity is computed directly by applying a thermal gradient in the liquid along the solid-liquid interface. We found a good agreement between the three methods, validating the relevance of Derjaguin's approach to study thermo-osmosis, for homogeneous surfaces. Even though the methods agree one with each other, they require different simulation times and are not equally accurate. For instance, the computation of the solid-liquid entropy using dry-surface method is much more accurate than the calculation of the enthalpy excess, requires less computational time, but is limited to electrically neutral surfaces. The NEMD calculation is very precise but requires to find a surface with a vanishing enthalpy excess, increasing the computational cost of the method.

We applied this methods to study the thermo-osmotic slip of charged graphitic surface. Although the neutral pristine graphene surface displays a relatively large thermo-osmotic coefficient due to its high slip length, we found that polarizing the graphene surface could significantly enhance the strength of the thermo-osmotic flow, by a factor 3 to 4 for the highest values. We found that this amplification is due to the variation of the enthalpy excess with the surface charge density. Our findings suggest that by controlling the surface charge density, one can precisely make the liquid flow in either directions in the channel, or not flow at all. Overall, these results call for experimental exploration of thermo-osmotic flow in the vicinity of electrically charged graphene surfaces.
Poisson-Boltzmann (PB) contribution to the enthalpy excess, and thus the thermo-osmotic coefficient, is negligible compared to the contribution of water molecules.

Figures 7 and 8 represent the thermo-osmotic coefficient for a neutral and a charged GO surface, respectively, along with their slip length. When considering defective surfaces, the presence of defects greatly reduces the solid-liquid slip⁵⁷, thereby diminishing the osmotic flow generated. The neutral GO surface still exhibits a relatively high thermo-osmotic coefficient for very small oxidation rates. The attenuation is more pronounced for charged defects, as the electrostatic interactions increase solid-liquid friction. In this case, the fluid flow quickly drops to values that are too small to be measured with our method.

The variation of the slip length with the surface density is given in Ref. 55 for a heterogeneous surface:

$$b = \frac{b_0}{1 + 3\pi\sigma_h b_0 |\Sigma| / q}, \quad (46)$$

with $b_0$ the slip length of the neutral surface, $\sigma_h$ the hydrodynamic radius of the defect plus counterion assembly, treated here as a fitting parameter. We used $\sigma_h = 8.8$ Å for our purpose. Similarly, for a neutral defective surface, one can write the slip length as:

$$b = \frac{b_0}{1 + 3\pi\sigma_h b_0 n_C \xi}, \quad (47)$$

with $n_C$ the number density of the carbon atoms, and $\sigma_h = 1.2$ Å the hydrodynamic radius of the hydroxyl groups. In both cases, the slip length quickly decreases even at low oxidation degree. Therefore, one can expect the experimental measured thermo-osmotic coefficient of a graphene nanochannel to be lower than what we have found numerically, due to the presence of defects.

## V. CONCLUSION

In this paper, we have presented three methods to compute the thermo-osmotic coefficient using molecular simulations. The first two are quite similar in that they are based on Derjaguin's theoretical framework. The first one requires to compute the enthalpy excess at equilibrium. The second one uses the dry-surface method to compute the solid-liquid entropy, which is linked to the enthalpy excess. The last method is based on NEMD, the thermo-osmotic velocity is computed directly by applying a thermal gradient in the liquid along the solid-liquid interface. We found a good agreement between the three methods, validating the relevance of Derjaguin's approach to study thermo-osmosis, for homogeneous surfaces. Even though the methods agree one with each other, they require different simulation times and are not equally accurate. For instance, the computation of the solid-liquid entropy using dry-surface method is much more accurate than the calculation of the enthalpy excess, requires less computational time, but is limited to

electrically neutral surfaces. The NEMD calculation is very precise but requires to find a surface with a vanishing enthalpy excess, increasing the computational cost of the method.

We applied this method to study the thermo-osmotic slip of charged graphitic surface. Although the neutral pristine graphene surface displays a relatively large thermo-osmotic coefficient due to its high slip length, we found that polarizing the graphene surface could significantly enhance the strength of the thermo-osmotic flow, by a factor 3 to 4 for the highest values. We found that this amplification is due to the variation of the enthalpy excess with the surface charge density. Our findings suggest that by controlling the surface charge density, one can precisely make the liquid flow in either directions in the channel, or not flow at all. Overall, these results call for experimental exploration of thermo-osmotic flow in the vicinity of electrically charged graphene surfaces.

## ACKNOWLEDGMENTS

The authors thank D. Pandey and S. Hardt for fruitful discussions. This work was supported by the ANR, Project ANR-21-CE50-0042 smoothE. This work used the HPC resources from the CNRS/IN2P3 Computing Center (Lyon - France) and from the PSMN mesocenter in Lyon. This manuscript is distributed under a Creative Commons Attribution | 4.0 International license.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

## APPENDIX A: EQUIVALENCE BETWEEN ENTROPY EXCESS AND SOLID-LIQUID ENTROPY

In the framework of surface thermodynamics³⁹, the solid-liquid entropy can be written as:

$$s_{\text{sl}} = \int_{0}^{\infty} dz \left[ \sum_{i} (n_{i}(z)s_{i}(z) - n_{i}^{\text{bulk}}s_{i}^{\text{bulk}}) \right]. \quad (48)$$

The entropy excess $\Delta S$ is then expressed as:

$$\Delta S = s_{\text{sl}} - \sum_{i} s_{i}^{\text{bulk}} \int_{0}^{\infty} dz (n_{i}(z) - n_{i}^{\text{bulk}}) \quad (49)$$

$$= s_{\text{sl}} - \sum_{i} s_{i}^{\text{bulk}} \Gamma_{i}, \quad (50)$$

where $\Gamma_i$ is the adsorption of species i. In principle, there is a difference between $\Delta S$ and $s_{sl}$, and this difference depends on the wetting properties of the surface, as the adsorption depends on the wetting. However, in practice,
## ACKNOWLEDGMENTS

The authors thank D. Pandey and S. Hardt for fruitful discussions. This work was supported by the ANR, Project ANR-21-CE50-0042 smoothE. This work used the HPC resources from the CNRS/IN2P3 Computing Center (Lyon - France) and from the PSMN mesocenter in Lyon. This manuscript is distributed under a Creative Commons Attribution | 4.0 International license.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

## APPENDIX A: EQUIVALENCE BETWEEN ENTROPY EXCESS AND SOLID-LIQUID ENTROPY

In the framework of surface thermodynamics³⁹, the solid-liquid entropy can be written as:

$$s_{sl} = \int_{0}^{\infty} dz \left[ \sum_{i} (n_{i}(z)s_{i}(z) - n_{i}^{\text{bulk}}s_{i}^{\text{bulk}}) \right]. \quad (49)$$

The entropy excess $\Delta S$ is then expressed as:

$$\Delta S = s_{sl} - \sum_{i} s_{i}^{\text{bulk}} \int_{0}^{\infty} dz (n_{i}(z) - n_{i}^{\text{bulk}}) \quad (50)$$

$$= s_{sl} - \sum_{i} s_{i}^{\text{bulk}} \Gamma_{i}, \quad (51)$$

where $\Gamma_i$ is the adsorption of species i. In principle, there is a difference between $\Delta S$ and $s_{sl}$, and this difference depends on the wetting properties of the surface, as the adsorption depends on the wetting. However, in practice, the calculation of the solid-liquid entropy matches very well with the enthalpy excess route, and with NEMD results (Fig. 5), suggesting an equivalence between $\Delta S$ and $s_{sl}$.

## APPENDIX B: CONTACT ANGLE COMPUTATION

To characterize the wetting of the solid-liquid interface, we performed simulations to compute the contact angle of liquid droplets on the graphene-like surfaces. We consider two systems of different size to quantify finite-size effects. The first system featured 1000 water molecules enclosed between graphene-like walls of size $L_x \sim L_y \sim 100$ Å. The walls were separated by a distance of 60 Å, to ensure that the droplet formed on one wall does not interact with the other wall. The second system was composed of 4000 water molecules. The size of walls was $L_x \sim L_y \sim 180$ Å. They were kept at a distance of 100 Å. We ran an equilibration phase for 0.5 ns,

<figure><img src="image_10.png" /><figcaption>FIG. 10. Contact angle as a function of the liquid-wall interaction parameter. We compare the contact angle for two size of droplets, 1000 water molecules (orange squares) and 4000 water molecules (blue circles).</figcaption></figure>

the production period lasted 4 ns. We computed the contact angle by fitting the average density profile along $z$ by⁵⁸:

$$\rho_{\text{app}}(z) = \frac{(\rho_{\text{l}} - \rho_{\text{v}})\pi r(z)^2}{L_{\text{x}}L_{\text{y}}} + \rho_{\text{v}}, \quad (52)$$

where $\rho_{\text{app}}$ is the apparent density profile, $\rho_{\text{l}}$ and $\rho_{\text{v}}$ are the real liquid and vapor densities, and $r(z)$ is the radius of the droplet section, which is equal to $\sqrt{R^2 - (z - z_0)^2}$, with $R$ the radius and $z_0$ the center of the sphere. The contact angle is then given by:

$$\theta = \arccos\left(\frac{z_{\text{s}} - z_{0}}{R}\right), \quad (53)$$

where $z_s$ is the hydrodynamic wall position⁵⁶.

The variation of the contact angle with $\epsilon_{CO}$ is presented in Fig. 10. The contact angle values given in the main text are for the 4000-molecule water droplet, which are more accurate because they are less affected by finite-size effects.

## APPENDIX C: OSMOTIC VELOCITY WITH A HETEROGENEOUS VISCOSITY

In this appendix, we provide extensions of Eq. (1), which relates the osmotic velocity to the force density exerted on the fluid in the $x$ direction, for a heterogeneous (local) viscosity. Denoting $z$ the distance to the wall, the Stokes equation writes:

$$-\eta(z) \frac{d^2 v}{dz^2} = f_x(z), \quad (54)$$

with $\eta(z)$ the viscosity, $f_x(z)$ the external force density, and $v(z)$ the velocity along the $x$ direction.
the calculation of the solid-liquid entropy matches very well with the enthalpy excess route, and with NEMD results (Fig. 5), suggesting an equivalence between $\Delta S$ and $s_{sl}$.

## APPENDIX B: CONTACT ANGLE COMPUTATION

To characterize the wetting of the solid-liquid interface, we performed simulations to compute the contact angle of liquid droplets on the graphene-like surfaces. We consider two systems of different size to quantify finite-size effects. The first system featured 1000 water molecules enclosed between graphene-like walls of size $L_x \sim L_y \sim 100$ Å. The walls were separated by a distance of 60 Å, to ensure that the droplet formed on one wall does not interact with the other wall. The second system was composed of 4000 water molecules. The size of walls was $L_x \sim L_y \sim 180$ Å. They were kept at a distance of 100 Å. We ran an equilibration phase for 0.5 ns, the production period lasted 4 ns. We computed the contact angle by fitting the average density profile along $z$ by$^{58}$:

$$ \rho_{\text{app}}(z) = \frac{(\rho_l - \rho_v)\pi r(z)^2}{L_x L_y} + \rho_v, \quad (51) $$

where $\rho_{app}$ is the apparent density profile, $\rho_l$ and $\rho_v$ are the real liquid and vapor densities, and $r(z)$ is the radius of the droplet section, which is equal to $\sqrt{R^2 - (z - z_0)^2}$, with $R$ the radius and $z_0$ the center of the sphere. The contact angle is then given by:

$$ \theta = \arccos\left(\frac{z_s - z_0}{R}\right), \quad (52) $$

where $z_s$ is the hydrodynamic wall position$^{56}$.

The variation of the contact angle with $\varepsilon_{CO}$ is presented in Fig. 10. The contact angle values given in the main text are for the 4000-molecule water droplet, which are more accurate because they are less affected by finite-size effects.

## VI. REFERENCES

¹R. B. Schoch, J. Han, and P. Renaud, Reviews of Modern Physics 80, 839 (2008).

²L. Bocquet and E. Charlaix, Chem. Soc. Rev. 39, 1073 (2010).

³N. Kavokine, R. R. Netz, and L. Bocquet, Annual Review of Fluid Mechanics 53, 377 (2021).

⁴D. C. Grahame, Chemical Reviews 41, 441 (1947).

⁵J. Israelachvili, Intermolecular and Surface Forces (Academic Press, 2011).

⁶R. Hartkamp, A.-L. Biance, L. Fu, J.-F. Dufrêche, O. Bonhomme, and L. Joly, Current Opinion in Colloid & Interface Science 37, 101 (2018).

⁷A. Siria, P. Poncharal, A.-L. Biance, R. Fulcrand, X. Blase, S. T. Purcell, and L. Bocquet, Nature 494, 455 (2013).

⁸J. C. Fair and J. F. Osterle, The Journal of Chemical Physics 54, 3307 (1971).

<figure><img src="image_11.png" /><figcaption>FIG. 10. Contact angle as a function of the liquid-wall interaction parameter. We compare the contact angle for two size of droplets, 1000 water molecules (orange squares) and 4000 water molecules (blue circles).</figcaption></figure>

⁹T. Mouterde and L. Bocquet, The European Physical Journal E 41, 148 (2018).

¹⁰V. M. Barragán and S. Kjelstrup, Journal of Non-Equilibrium Thermodynamics 42, 10.1515/jnet-2016-0088 (2017).

¹¹A. Würger, Reports on Progress in Physics 73, 126601 (2010).

¹²C. Herrero, A. Allemand, S. Merabia, A.-L. Biance, and L. Joly, arXiv preprint arXiv:2204.13522 (2022).

¹³B. V. Derjaguin, N. V. Churaev, and V. M. Muller, Surface Forces (Springer US, Boston, MA, 1987).

¹⁴B. V. Derjaguin and G. P. Sidorenkov, CR Acad. Sci. URSS 32, 622 (1941).

¹⁵J. Anderson, Annual Review of Fluid Mechanics 21, 61 (1989).

¹⁶K. P. Travis, B. D. Todd, and D. J. Evans, Physical Review E 55, 4288 (1997).

¹⁷C. Sendner, D. Horinek, L. Bocquet, and R. R. Netz, Langmuir 25, 10768 (2009).

¹⁸D. J. Bonthuis and R. R. Netz, The Journal of Physical Chemistry B 117, 11397 (2013).

¹⁹T.-D. Li, J. Gao, R. Szoszkiewicz, U. Landman, and E. Riedo, Physical Review B 75, 115415 (2007).

²⁰R. Ganti, Y. Liu, and D. Frenkel, Physical Review Letters 119, 038002 (2017).

²¹Y. Leng and P. T. Cummings, Physical Review Letters 94, 026101 (2005).

²²K. Shi, E. R. Smith, E. E. Santiso, and K. E. Gubbins, The Journal of Chemical Physics 158, 040901 (2023).

²³R. Ganti, Y. Liu, and D. Frenkel, Physical Review Letters 121, 068002 (2018).

²⁴P. Anzini, Z. Filiberti, and A. Parola, Physical Review E 106, 024116 (2022).

²⁵C. Herrero, M. De San Féliciano, S. Merabia, and L. Joly, Nanoscale 14, 626 (2022).

²⁶L. Fu, S. Merabia, and L. Joly, Physical Review Letters 119, 214501 (2017).

²⁷W. Q. Chen, M. Sedighi, and A. P. Jivkov, Nanoscale 13, 1696 (2021).

²⁸W. Q. Chen, A. P. Jivkov, and M. Sedighi, ACS Applied Materials & Interfaces, acsami.3c02559 (2023).

²⁹K. Proesmans and D. Frenkel, The Journal of Chemical Physics 151, 124109 (2019).

³⁰M. Ouadfel, M. De San Féliciano, C. Herrero, S. Merabia, and L. Joly, Physical Chemistry Chemical Physics 25, 24321 (2023).

³¹D. Surblys, F. Leroy, Y. Yamaguchi, and F. Müller-Plathe, The Journal of Chemical Physics 148, 134707 (2018).
Denoting $b_0$ the local slip length at the wall, $b_0 = \eta(0)/\lambda$, with $\lambda$ the interfacial friction coefficient, the slip boundary condition at the wall writes:

$$v(0) = b_0 \left. \frac{dv}{dz} \right|_{z=0} . \qquad (55)$$

Denoting $\eta_{\text{bulk}}$ the viscosity far from the wall, one can solve Eq. (54) with the slip boundary condition and a vanishing shear rate far from the wall to obtain the osmotic velocity, i.e. the plateau of velocity far from the wall:

$$v_{\text{osm}} = \frac{1}{\eta_{\text{bulk}}} \int_{0}^{\infty} (z + b_0) f_x(z) \frac{\eta_{\text{bulk}}}{\eta(z)} dz. \quad (56)$$

Similarly, for a no-slip surface with a stagnant layer of thickness $z_s$, the osmotic velocity writes:

$$v_{\text{osm}} = \frac{1}{\eta_{\text{bulk}}} \int_{z_s}^{\infty} (z - z_s) f_x(z) \frac{\eta_{\text{bulk}}}{\eta(z)} dz. \quad (57)$$

Therefore, any local increase in the viscosity amounts to effectively reducing the driving force density by a factor $\eta(z)/\eta_{\text{bulk}}$ as compared to Eq. (1) assuming a homogeneous viscosity.

## VI. REFERENCES

¹R. B. Schoch, J. Han, and P. Renaud, Reviews of Modern Physics 80, 839 (2008).

²L. Bocquet and E. Charlaix, Chem. Soc. Rev. 39, 1073 (2010).

³N. Kavokine, R. R. Netz, and L. Bocquet, Annual Review of Fluid Mechanics 53, 377 (2021).

⁴D. C. Grahame, Chemical Reviews 41, 441 (1947).

⁵J. Israelachvili, *Intermolecular and Surface Forces* (Academic Press, 2011).

⁶R. Hartkamp, A.-L. Biance, L. Fu, J.-F. Dufrêche, O. Bonhomme, and L. Joly, Current Opinion in Colloid & Interface Science 37, 101 (2018).

⁷A. Siria, P. Poncharal, A.-L. Biance, R. Fulcrand, X. Blase, S. T. Purcell, and L. Bocquet, Nature 494, 455 (2013).

⁸J. C. Fair and J. F. Osterle, The Journal of Chemical Physics 54, 3307 (1971).

⁹T. Mouterde and L. Bocquet, The European Physical Journal E 41, 148 (2018).

¹⁰V. M. Barragán and S. Kjelstrup, Journal of Non-Equilibrium Thermodynamics 42, 10.1515/jnet-2016-0088 (2017).

¹¹A. Würger, Reports on Progress in Physics 73, 126601 (2010).

¹²C. Herrero, A. Allemand, S. Merabia, A.-L. Biance, and L. Joly, arXiv preprint arXiv:2204.13522 (2022).

¹³B. V. Derjaguin, N. V. Churaev, and V. M. Muller, *Surface Forces* (Springer US, Boston, MA, 1987).

¹⁴B. V. Derjaguin and G. P. Sidorenkov, CR Acad. Sci. URSS 32, 622 (1941).

¹⁵J. Anderson, Annual Review of Fluid Mechanics 21, 61 (1989).

¹⁶K. P. Travis, B. D. Todd, and D. J. Evans, Physical Review E 55, 4288 (1997).

¹⁷C. Sendner, D. Horinek, L. Bocquet, and R. R. Netz, Langmuir 25, 10768 (2009).

¹⁸D. J. Bonthuis and R. R. Netz, The Journal of Physical Chemistry B 117, 11397 (2013).

¹⁹T.-D. Li, J. Gao, R. Szoszkiewicz, U. Landman, and E. Riedo, Physical Review B 75, 115415 (2007).

²⁰R. Ganti, Y. Liu, and D. Frenkel, Physical Review Letters 119, 038002 (2017).

²¹Y. Leng and P. T. Cummings, Physical Review Letters 94, 026101 (2005).

²²K. Shi, E. R. Smith, E. E. Santiso, and K. E. Gubbins, The Journal of Chemical Physics 158, 040901 (2023).

²³R. Ganti, Y. Liu, and D. Frenkel, Physical Review Letters 121, 068002 (2018).

²⁴P. Anzini, Z. Filiberti, and A. Parola, Physical Review E 106, 024116 (2022).

²⁵C. Herrero, M. De San Féliciano, S. Merabia, and L. Joly, Nanoscale 14, 626 (2022).

²⁶L. Fu, S. Merabia, and L. Joly, Physical Review Letters 119, 214501 (2017).

²⁷W. Q. Chen, M. Sedighi, and A. P. Jivkov, Nanoscale 13, 1696 (2021).

²⁸W. Q. Chen, A. P. Jivkov, and M. Sedighi, ACS Applied Materials & Interfaces, acsami.3c02559 (2023).

²⁹K. Proesmans and D. Frenkel, The Journal of Chemical Physics 151, 124109 (2019).

³⁰M. Ouadfel, M. De San Féliciano, C. Herrero, S. Merabia, and L. Joly, Physical Chemistry Chemical Physics 25, 24321 (2023).

³¹D. Surblys, F. Leroy, Y. Yamaguchi, and F. Müller-Plathe, The Journal of Chemical Physics 148, 134707 (2018).

³²F. Leroy and F. Müller-Plathe, Langmuir 31, 8335 (2015).

³³V. R. Ardham, G. Deichmann, N. F. A. Van Der Vegt, and F. Leroy, The Journal of Chemical Physics 143, 243135 (2015).

³⁴L. Bocquet and J.-L. Barrat, Soft Matter 3, 685 (2007).

³⁵A. P. Bregulla, A. Würger, K. Günther, M. Mertig, and F. Cichos, Physical Review Letters 116, 188303 (2016).

³⁶A. Catenaccio, Y. Daruich, and C. Magallanes, Chemical Physics Letters 367, 669 (2003).

³⁷C. Herrero and L. Joly, Poisson-Boltzmann formulary: Second edition (2022), arXiv:2105.00720 [cond-mat].

³⁸R. Blossey, *The Poisson-Boltzmann Equation: An Introduction*, SpringerBriefs in Physics (Springer International Publishing, Cham, 2023).

³⁹J. S. Rowlinson and B. Widom, *Molecular theory of capillarity*, dover ed ed. (Dover Publications, Mineola, N.Y, 2002).

⁴⁰A. Gao, H.-J. Butt, W. Steffen, and C. Schönecker, Langmuir 37, 8677 (2011).

⁴¹A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Computer Physics Communications 271, 108171 (2022).

⁴²I. M. Zeron, J. L. F. Abascal, and C. Vega, The Journal of Chemical Physics 151, 134504 (2019).

⁴³J. L. F. Abascal and C. Vega, The Journal of Chemical Physics 123, 234505 (2005).

⁴⁴K. Falk, F. Sedlmeier, L. Joly, R. R. Netz, and L. Bocquet, Nano Letters 10, 4067 (2010).

⁴⁵G. Yang, L. Li, W. B. Lee, and M. C. Ng, Science and Technology of Advanced Materials 19, 613 (2018).

⁴⁶D. Bahamon and L. F. Vega, Molecular Physics 117, 3703 (2019).

⁴⁷N. Wei, X. Peng, and Z. Xu, Physical Review E 89, 012113 (2014).

⁴⁸A. K. Giri, F. Teixeira, and M. N. D. Cordeiro, Desalination 460, 1 (2019).

⁴⁹W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, Journal of the American Chemical Society 118, 11225 (1996).

⁵⁰G. Bussi, D. Donadio, and M. Parrinello, The Journal of Chemical Physics 126, 014101 (2007).

⁵¹A. Stukowski, Modelling and Simulation in Materials Science and Engineering 18, 015012 (2010).

⁵²L. Bocquet and J.-L. Barrat, The Journal of Chemical Physics 139, 044704 (2013).

⁵³D. Alfé and M. J. Gillan, Physical Review Letters 81, 5161 (1998).

⁵⁴B. Cross, C. Barraud, C. Picard, L. Léger, F. Restagno, and
³²F. Leroy and F. Müller-Plathe, *Langmuir* **31**, 8335 (2015).

³³V. R. Ardham, G. Deichmann, N. F. A. Van Der Vegt, and F. Leroy, *The Journal of Chemical Physics* **143**, 243135 (2015).

³⁴L. Bocquet and J.-L. Barrat, *Soft Matter* **3**, 685 (2007).

³⁵A. P. Bregulla, A. Würger, K. Günther, M. Mertig, and F. Cichos, *Physical Review Letters* **116**, 188303 (2016).

³⁶A. Catenaccio, Y. Daruich, and C. Magallanes, *Chemical Physics Letters* **367**, 669 (2003).

³⁷C. Herrero and L. Joly, *Poisson-Boltzmann formulary: Second edition* (2022), arXiv:2105.00720 [cond-mat].

³⁸R. Blossey, *The Poisson-Boltzmann Equation: An Introduction*, SpringerBriefs in Physics (Springer International Publishing, Cham, 2023).

³⁹J. S. Rowlinson and B. Widom, *Molecular theory of capillarity*, dover ed ed. (Dover Publications, Mineola, N.Y., 2002).

⁴⁰A. Gao, H.-J. Butt, W. Steffen, and C. Schönecker, *Langmuir* **37**, 8677 (2021).

⁴¹A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, *Computer Physics Communications* **271**, 108171 (2022).

⁴²I. M. Zeron, J. L. F. Abascal, and C. Vega, *The Journal of Chemical Physics* **151**, 134504 (2019).

⁴³J. L. F. Abascal and C. Vega, *The Journal of Chemical Physics* **123**, 234505 (2005).

⁴⁴K. Falk, F. Sedlmeier, L. Joly, R. R. Netz, and L. Bocquet, *Nano Letters* **10**, 4067 (2010).

⁴⁵G. Yang, L. Li, W. B. Lee, and M. C. Ng, *Science and Technology of Advanced Materials* **19**, 613 (2018).

⁴⁶D. Bahamon and L. F. Vega, *Molecular Physics* **117**, 3703 (2019).

⁴⁷N. Wei, X. Peng, and Z. Xu, *Physical Review E* **89**, 012113 (2014).

⁴⁸A. K. Giri, F. Teixeira, and M. N. D. Cordeiro, *Desalination* **460**, 1 (2019).

⁴⁹W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, *Journal of the American Chemical Society* **118**, 11225 (1996).

⁵⁰G. Bussi, D. Donadio, and M. Parrinello, *The Journal of Chemical Physics* **126**, 014101 (2007).

⁵¹A. Stukowski, *Modelling and Simulation in Materials Science and Engineering* **18**, 015012 (2010).

⁵²L. Bocquet and J.-L. Barrat, *The Journal of Chemical Physics* **139**, 044704 (2013).

⁵³D. Alfè and M. J. Gillan, *Physical Review Letters* **81**, 5161 (1998).

⁵⁴B. Cross, C. Barraud, C. Picard, L. Léger, F. Restagno, and E. Charlaix, *Physical Review Fluids* **3**, 062001 (2018).

⁵⁵Y. Xie, L. Fu, T. Niehaus, and L. Joly, *Physical Review Letters* **125**, 014501 (2020).

⁵⁶C. Herrero, T. Omori, Y. Yamaguchi, and L. Joly, *The Journal of Chemical Physics* **151**, 041103 (2019).

⁵⁷E. Mangaud, M.-L. Bocquet, L. Bocquet, and B. Rotenberg, *The Journal of Chemical Physics* **156**, 044703 (2022).

⁵⁸L. Fu, L. Joly, and S. Merabia, *Physical Review Letters* **123**, 138001 (2019).
E. Charlaix, Physical Review Fluids 3, 062001 (2018).

<footnote n="55">⁵⁵Y. Xie, L. Fu, T. Niehaus, and L. Joly, Physical Review Letters **125**, 014501 (2020).</footnote>

<footnote n="56">⁵⁶C. Herrero, T. Omori, Y. Yamaguchi, and L. Joly, The Journal of Chemical Physics **151**, 041103 (2019).</footnote>

<footnote n="57">⁵⁷E. Mangaud, M.-L. Bocquet, L. Bocquet, and B. Rotenberg, The Journal of Chemical Physics **156**, 044703 (2022).</footnote>

<footnote n="58">⁵⁸L. Fu, L. Joly, and S. Merabia, Physical Review Letters **123**, 138001 (2019).</footnote>