# Mechanical and thermodynamic routes to the liquid-liquid interfacial tension and mixing free energy by molecular dynamics

Rei Ogawa, Hiroki Kusudo, Takeshi Omori, Edward Smith, Laurent Joly, Samy Merabia, Yasutaka Yamaguchi

► To cite this version:

Rei Ogawa, Hiroki Kusudo, Takeshi Omori, Edward Smith, Laurent Joly, et al.. Mechanical and thermodynamic routes to the liquid-liquid interfacial tension and mixing free energy by molecular dynamics. The Journal of Chemical Physics, 2024, 161 (22), pp.224708. 10.1063/5.0238862. hal-04875797

HAL Id: hal-04875797

https://hal.science/hal-04875797v1

Submitted on 9 Jan 2025

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
# Mechanical and thermodynamic routes to the liquid-liquid interfacial tension and mixing free energy by molecular dynamics

Rei Ogawa ⚫; Hiroki Kusudo ⚫; Takeshi Omori ⚫; Edward R. Smith ⚫; Laurent Joly ⚫; Samy Merabia ⚫; Yasutaka Yamaguchi 🕥

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

Check for updates

J. Chem. Phys. 161, 224708 (2024)

https://doi.org/10.1063/5.0238862

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

CHORUS

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

View
Online

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

Export
Citation

## Articles You May Be Interested In

Measuring line tension: Thermodynamic integration during detachment of a molecular dynamics droplet

J. Chem. Phys. (June 2024)

Molecular anatomy of the pressure anisotropy in the interface of one and two component fluids: Local thermodynamic description of the interfacial tension

J. Chem. Phys. (November 2024)

Interpretation of Young's equation for a liquid droplet on a flat and smooth solid surface: Mechanical and thermodynamic routes with a simple Lennard-Jones liquid

J. Chem. Phys. (January 2019)

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

AIP
Publishing

The Journal of Chemical Physics

Special Topics Open
for Submissions

Learn More

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

AIP
Publishing
# Mechanical and thermodynamic routes to the liquid-liquid interfacial tension and mixing free energy by molecular dynamics

Cite as: J. Chem. Phys. 161, 224708 (2024); doi: 10.1063/5.0238862

Submitted: 15 September 2024 • Accepted: 26 November 2024 •

Published Online: 12 December 2024

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

View Online

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

Export Citation

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

CrossMark

Rei Ogawa,<fnref n="1" /> ID Hiroki Kusudo,<fnref n="2" /> ID Takeshi Omori,<fnref n="3" /> ID Edward R. Smith,<fnref n="4" /> ID Laurent Joly,<fnref n="5" /> ID
Samy Merabia,<fnref n="5" /> ID and Yasutaka Yamaguchi<fnref n="1" />,<fnref n="6" />,α) ID

## AFFILIATIONS

## ABSTRACT

In this study, we carried out equilibrium molecular dynamics (EMD) simulations of the liquid-liquid (LL) interface between two different Lennard-Jones components with varying miscibility, where we examined the relation between the interfacial tension and the free energy to completely isolate the two liquids using both a mechanical and thermodynamic approach. Using the mechanical approach, we obtained a stress distribution around a quasi-one-dimensional EMD system with a flat LL interface. From the stress distribution, we calculated the LL interfacial tension based on Bakker's equation, which uses the stress anisotropy around the interface, and measured how it varied with miscibility. The second approach uses thermodynamic integration by enforcing quasi-static isolation of the two liquids to calculate the free energy. This uses the same EMD systems as the mechanical approach, with both extended dry-surface and phantom-wall (PW) schemes applied. When the two components were immiscible, the mechanical interfacial tension and isolation free energy were in good agreement. When the components were miscible, the values were significantly different. From the result of PW for the case of completely mixed liquids, the difference was attributed to the additional free energy required to separate the binary mixture into single components against the osmotic pressure prior to the complete detachment of the two components. This provides a new route to obtain the free energy of mixing.

© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0238862

## I. INTRODUCTION

When two immiscible liquids, such as water and oil, coexist, an interface is typically formed between them, referred to as a liquid-liquid (LL) interface, where a LL interfacial tension arises. Such LL interface can be found in emulsions, which are ubiquitous in our daily lives, e.g., food, drink, and cosmetics, so measuring the LL interfacial free energy and/or tension is important to control the properties of the mixture.

From a microscopic point of view, Kirkwood and Buff<fnref n="1" />,<fnref n="2" /> formulated expressions for the chemical potentials of the components of gas mixtures and liquid solutions based on statistical mechanics, and they also provided a general theoretical framework to describe surface tension.<fnref n="3" /> Related to this point, with respect to a liquid-vapor (LV) or liquid-gas (LG) interface, Bakker's equation,<fnref n="4" /> based on macroscopic thermodynamics, was known before the formulation from the statistical mechanics of Kirkwood and Buff. This equation relates the macroscopic LV or LG interfacial tension to the

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

<footnote n="2">²Department of Mechanical Systems Engineering, Tohoku University, 6-6-01 Aramaki Aoba-ku, Sendai 980-8579, Japan</footnote>

<footnote n="3">³Department of Mechanical Engineering, Osaka Metropolitan University, 1-1 Gakuen-cho, Nakaku, Sakai, Osaka 599-8531, Japan</footnote>

<footnote n="4">⁴Department of Mechanical and Aerospace Engineering, Brunel University London, Uxbridge UB8 3PH, United Kingdom</footnote>

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

<footnote n="6">⁶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) Author to whom correspondence should be addressed: yamaguchi@mech.eng.osaka-u.ac.jp</footnote>
anisotropy of the microscopic local stress near the interface. For a flat LV interface normal to the x-axis, Bakker's equation is written as

$$ \gamma_{LV} = \int_{x_V^{\text{blk}}}^{x_L^{\text{blk}}} [\tau_{yy}(x) - \tau_{xx}] dx = \int_{x_V^{\text{blk}}}^{x_L^{\text{blk}}} \tau_{yy}(x) dx - \int_{x_V^{\text{blk}}}^{x_L^{\text{blk}}} \tau_{xx} dx, \quad (1) $$

where $\gamma_{LV}$ is the LV interfacial tension, and $\tau_{yy} (= \tau_{zz})$ and $\tau_{xx}$ are the diagonal components of the stress tensor tangential and normal to the interface, respectively. From the local-force balance in the direction normal to the interface, at equilibrium, $\tau_{xx}$ is constant over the entire region and equal to the isotropic pressure in the bulk with its sign inverted. By integrating the stress difference $\tau_{yy}(x) - \tau_{xx}$ existing only around the interface, i.e., by taking $x_L^{\text{blk}}$ and $x_V^{\text{blk}}$ in Eq. (1) as the bulk positions of the liquid and gas phases, respectively, the LV interfacial tension $\gamma_{LV}$ is obtained.

At present, molecular dynamics (MD) simulation is a powerful tool to investigate LV or LG interfaces as well as LL mixtures composed of various kinds of molecular pairs in silico. A quasi-one-dimensional (1D) MD system can be easily simulated, e.g., one with LV coexistence can be constructed by confining a single molecular component in a constant-volume simulation cell with periodic boundary conditions (PBCs) at a temperature between the triple point and the critical point. In such an MD system, the two integrals in the right-hand-side of Bakker's Eq. (1) can be obtained,⁵ and thus, Eq. (1) is widely used to calculate the surface tension as a standard approach in MD. This is partly because only the integral of each principal stress component in the whole system is used, which can be easily calculated in such a quasi-1D system.

Regarding the LL interface, in the early stage of MD development, Hayoun et al.⁶ simulated a quasi-1D system with an interface between two Lennard-Jones (LJ) liquids and showed the density and pressure distribution, and Benjamin⁷ wrote a review article of MD studies regarding the structure and dynamics at the LL interface. At present, the LL interfacial tension is also calculated using the integrated form of Bakker's equation similar to Eq. (1) as a definition.⁸ Furthermore, the authors have successfully extended Bakker's Eq. (1) to flat solid-liquid (SL) and solid-vapor (SV) interfaces. In this framework, through a careful choice of the SL and SV interface positions based on a mechanical force balance, the microscopic interpretation of Young's equation was clarified as the force balance exerted on the fluid particles in a finite region surrounding the contact line. It should be noted that the microscopic stress calculated by the method of planes (MoPs)⁹ was evaluated only by including the interaction forces between fluid particles, while the force from the solid particles on the fluid particles was treated as an external force. Such a treatment of interfacial tensions, including surface tension, based on a mechanical force balance is called the mechanical route.¹⁰⁻¹⁴

In addition to the mechanical route, the SL and SV interfacial tensions have been calculated as the free energy per unit area of the interface based on thermodynamic integration methods.¹³,¹⁵⁻¹⁷ In this thermodynamic route, the solid-fluid (SF) work of adhesion was calculated by either quasi-statically pushing the liquid away from the solid using a virtual wall that interacted only with the fluid (called the phantom-wall method) or by gradually reducing the SF interaction strength until the solid and fluid no longer interacted (called the dry-surface method). More concretely, the quasi-static SL isolation process was carried out while keeping constant the number of

particles N, pressure p, and temperature T, i.e., under constant $NpT$ conditions, and the SL work of adhesion $W_{SL}$ was obtained from the change in Gibbs free energy G of the system given by

$$ W_{SL} \equiv \frac{\Delta G}{A} = \gamma_{S0} + \gamma_{L0} - \gamma_{SL}, \quad (2) $$

where A is the interfacial area and $\gamma$ is the interfacial free energy per unit area with corresponding interfaces denoted by the subscripts (S: solid, L: liquid, and 0: depletion layer). Note that $\Delta G$ excluded the work exerted on the environment.¹³ The solid-vapor work of adhesion $W_{SV}$ was also evaluated by similar schemes. As a result, it was shown that the mechanically obtained SL and SV interfacial tensions and thermodynamically obtained works of adhesion agreed well for a LJ fluid on a simple crystal surface.

In this study, we investigate the possible extension of the mechanical and thermodynamic routes to the LL interface with a focus on the following points: (1) whether Bakker's equation as a mechanical route can be extended to the LL interface, and (2) whether the mechanical route corresponds to the thermodynamic route for the LL interface. Related to the second point, we tested two methods for the thermodynamic route to validate the results. In addition, we show that one of the thermodynamic methods can provide a clear insight into the isolation free energy, i.e., the free energy needed to separate the mixture into single components—which includes the free energy of mixing and the change in interfacial free energy. This also gives access to the osmotic pressure of the liquid mixtures.

## II. METHOD

### A. Simulation setup

All the simulations were carried out using an in-house code, which was used in our previous work and thoroughly validated for cases of LJ fluids with interfaces, including those with a contact line.¹³,¹⁷⁻²³ Figure 1 shows the quasi-1D equilibrium MD (EMD) system studied. Two kinds of LJ particles, denoted by $\alpha$ and $\beta$, having the same mass of $m = 6.634 \times 10^{-26}$ kg, equal to the mass of an argon atom of 40 amu, were confined between two walls on the left and right ends, both parallel to the yz-plane shown in light brown. The numbers of particles $\alpha$ and $\beta$ were both set equal to 2000. Although these LJ particles can also form a vapor phase, we write them as “liquid particles” for clarity in the following as we are only considering liquids. The interaction potentials between the same components $\phi_{\alpha\alpha}$ and $\phi_{\beta\beta}$ were expressed as

$$ \phi_{\alpha\alpha}(r_{ij}) = \phi_{\beta\beta}(r_{ij}) = 4\epsilon \left[ \left( \frac{\sigma}{r_{ij}} \right)^{12} - \left( \frac{\sigma}{r_{ij}} \right)^{6} \right], \quad (3) $$

where $r_{ij}$ denotes the distance between particles i and j of the same component $\alpha$ or $\beta$, and $\sigma$ and $\epsilon$ were the LJ length and energy potential parameters, which we set as $\sigma = 0.34$ nm and $\epsilon = 1.67 \times 10^{-21}$ J. The interaction potential between different components $\phi_{\alpha\beta}(r_{ij})$ was expressed by multiplying the interaction potential in Eq. (3) by $\eta$ as follows:

$$ \phi_{\alpha\beta}(r_{ij}) = \eta\phi_{\alpha\alpha}(r_{ij}), \quad (4) $$

where $\eta$ was set between 0.01 and 1 as a variable parameter.
<figure><img src="image_11.png" /><figcaption>FIG. 1. (a) Setup of the equilibrium simulation system with a liquid-liquid interface controlled by a miscibility parameter $\eta$. A pressure control piston (PC) and a semipermeable phantom wall (PW) were located on each end of the system in the x-direction (left and right). (b) Side snapshots of the equilibrium systems for the free-energy calculation by the thermodynamic integration using the (i) extended dry surface (DS) and (ii) extended PW methods. For the extended DS method (i), the miscibility parameter $\eta$ was varied while keeping the PWs at $x_1^{PW}$ and $-x_1^{PW}$ far from the liquid; thus, the PWs as well as the PCs are shown only on the top panel of (i). For the PW method (ii), the PW positions $x_1^{PW}$ and $-x_1^{PW}$ were varied while keeping $\eta$ unchanged.</figcaption></figure>

Periodic boundary conditions were adopted in y- and z-directions, and the cell size in these directions, $L_y$ and $L_z$, were both 4 nm. In the following, the system yz cross sectional area is denoted $A = L_yL_z$. In addition, we located a pressure control (PC) wall and a semipermeable phantom wall (PW) on each end of the system in the x-direction (left and right), all of which were parallel to the yz-plane and interacted with the liquid particles as a unique function of the distance given by

$$ \phi_{\text{fw}}(x_i') = \frac{4}{5} \pi \varepsilon_{\text{fw}} \rho_s x_i'^2 \left( \frac{\sigma_{\text{fw}}}{x_i'} \right)^{12}, \quad (5) $$

where $x_i'$ is the distance between particle i at $x_i$ and corresponding wall. This potential field corresponds to a mean potential field formed by a single wall layer of uniformly distributed solid particles with an area number density $\rho_s = (3.61)^2 \text{ nm}^{-2}$, which interact with the liquid particles through the LJ potential only with the repulsive term [Eq. (3)] without $(\sigma/r_{ij})^6$ in the RHS] with the energy and length parameters being $\varepsilon_{\text{fw}} = 1.29 \times 10^{-21} \text{ J}$ and $\sigma_{\text{fw}} = 0.345 \text{ nm}$, respectively. The presented results are not sensitive to the choice

of parameter values, which are set the same as in our previous studies,¹³,¹⁷ as long as the interaction is repulsive and short-ranged.

For the PCs on the left and right at $x_{\text{left}}^{\text{PC}}$ and $x_{\text{right}}^{\text{PC}}$, respectively, shown in light brown in Fig. 1(a), we used

$$ \phi_{\text{liquid-pc}}(x_i') = \phi_{\text{fw}}(x_i'), \qquad (6) \\ x_i' = x_i - x_{\text{left}}^{\text{PC}}, \quad x_i' = x_{\text{right}}^{\text{PC}} - x_i \quad (i \in \alpha, \beta). $$

By adjusting the positions $x_{\text{left}}^{\text{PC}}$ and $x_{\text{right}}^{\text{PC}}$ of the walls as pistons, the system pressure was maintained at a constant value of $p_{\text{set}} \approx 1 \text{ MPa}$.

On the other hand, for the PWs on the left and right at $x_{\alpha}^{\text{PW}}$ and $x_{\beta}^{\text{PW}}$, respectively, we applied

$$ \phi_{\text{liquid-pw}\alpha}(x_i') = \phi_{\text{fw}}(x_i'), \quad x_i' = x_{\alpha}^{\text{PW}} - x_i \quad (i \in \alpha) \quad (7) $$

and

$$ \phi_{\text{liquid-pw}\beta}(x_i') = \phi_{\text{fw}}(x_i'), \quad x_i' = x_i - x_{\beta}^{\text{PW}} \quad (i \in \beta) \quad (8) $$

with setting

$$ x_{\alpha}^{\text{PW}} = x^{\text{PW}}, \quad \text{and} \quad x_{\beta}^{\text{PW}} = -x^{\text{PW}}. \quad (9) $$
With this setting, the PWs interacted only with α or β and worked as semipermeable membranes located symmetric to the yz-plane.

The system temperature T was controlled at 85 K by applying a velocity rescaling thermostat to the liquid particles located at less than 1.5 nm from the potential walls, only for the velocity components in the y- and z-directions. These thermostat regions were sufficiently far away from the interface, and no direct thermostating was applied to the region near the interface so that this thermostat had no effects on the presented results.

For the calculation of the stress distribution as a mechanical route and for the dry surface (DS) method as a thermodynamic route shown in Fig. 1(b-i), the symmetric PW position was fixed at $x^{\text{PW}}(x_\alpha^{\text{PW}} = -x_\beta^{\text{PW}}) = 8$ nm, at which the PWs were sufficiently far away from the liquid so that they did not interact with the liquids. We denote this value of $x^{\text{PW}} = 8$ nm as $x_1^{\text{PW}}$ hereafter. With this setting, a quasi-1D system under constant $NpT$ was obtained as an equilibrium state for each $\eta$ value. The two liquids were completely mixed without interface at $\eta = 1$ because both liquids are identical. By decreasing $\eta$, the two liquids were separated at $\eta = 0.85$ and formed a flat LL interface, and at $\eta = 0.01$, the two liquids were isolated with an empty region between two liquids. This value of $\eta = 0.01$ is denoted as $\eta_0$ hereafter.

For the extended phantom wall (PW) method in Fig. 1(b-ii) as another thermodynamic route described below, the PW position $x^{\text{PW}}$ was changed while keeping $\eta$ unchanged. With the decrease of $x^{\text{PW}}$, the two liquids were separated by the semipermeable PWs, and they were isolated at $x^{\text{PW}} = 0.2$ nm in this case. This value of $x^{\text{PW}} = 0.2$ nm is denoted as $x_0^{\text{PW}}$ hereafter.

## B. Mechanical route

Here, we describe how Bakker's equation is used to calculate the LL interfacial tension from the stress distribution using a mechanical route. Equation (1) extended to the LL interface can be written as

$$ \gamma_{\alpha\beta} = \int_{x_{\alpha}^{\text{blk}}}^{x_{\beta}^{\text{blk}}} [\tau_{yy}(x) - \tau_{xx}] dx, \quad (10) $$

where $\gamma_{\alpha\beta}$ is the LL interfacial tension, and $x_{\alpha}^{\text{blk}}$ and $x_{\beta}^{\text{blk}}$ ($> x_{\alpha}^{\text{blk}}$) are the bulk positions of the liquid phases α and β, respectively, where the stress is isotropic.

We chose to compute explicitly the stress distribution rather than just the stress integral, as is usually done,⁸,²⁴ because the stress distribution provides more insight into the molecular origin of surface tension.²⁵⁻²⁷ Note that, while the stress integral is always well-defined and can be computed using standard MD packages such as LAMMPS²⁸ and GROMACS,²⁹ the local stress lacks a unique definition. Only certain approaches, such as the volume average (VA) method or the method of planes (MoPs), ensure that $\tau_{xx}$ is constant over the entire region.⁹,¹³,¹⁸,³⁰ As detailed below, we used the VA method implemented in our in-house code in this work. We emphasize that, regardless of the local stress definition, the stress integral is always the same²⁴ and provides the same interfacial tension as obtained through the stress distribution.

In general, the stress tensor $\tau$ is expressed by the sum of a kinetic term $\tau^{\text{kin}}$ and an interaction term $\tau^{\text{int}}$ as

$$ \tau = \tau^{\text{kin}} + \tau^{\text{int}}. \quad (11) $$

In this study, the liquid was composed of two monoatomic components, α and β, and we included all contributions from the two components in each term. At first, the kinetic energy term $\tau^{\text{kin}}$ can be written as

$$ \tau^{\text{kin}} = \tau^{\text{kin},\alpha} + \tau^{\text{kin},\beta}, \quad (12) $$

where the superscript “kin,α,” for instance, denotes the contribution from α particles to the kinetic term. On the other hand, the interaction term $\tau^{\text{int}}$ is written as

$$ \tau^{\text{int}} = \tau^{\text{int},\alpha\alpha} + \tau^{\text{int},\beta\beta} + \tau^{\text{int},\alpha\beta}, \quad (13) $$

where the superscript “int,αβ” denotes, for instance, the contribution from the intermolecular interaction between α and β particles to the interaction term. We verified that the stress definition with Eqs. (11)–(13) satisfied the local mechanical balance

$$ \nabla \cdot \tau = 0 \quad (14) $$

in equilibrium systems at an arbitrary point in the absence of the external field,³⁰ i.e., except near the PCs and PWs. In practice, we calculated the stress distribution by the volume average (VA) method¹⁸,³⁰ and checked if

$$ \tau_{xx} = \text{const} \quad (15) $$

was satisfied in the present quasi-1D system.

We applied the VA approach for local flat regions with a thickness $\delta x$ and a volume $V^{\text{CV}}$ of $A\delta x$, where the superscript “CV” stands for “control volume.” The two terms of the kinetic contribution in the RHS of Eq. (12) were expressed by

$$ \tau^{\text{kin},\xi} = - \frac{1}{V^{\text{CV}}} \left\langle \sum_{i \in \xi}^{N_{\xi}} m_i \boldsymbol{v}_i \boldsymbol{v}_i \theta_i \right\rangle, \quad (\xi = \alpha, \beta), \quad (16) $$

where $r_i$ and $v_i$ are the position and velocity vectors of ith $\xi$ particle, and $\theta_i$ is a function that is one if the particle is in the local volume or zero otherwise, and the summation is taken over all $N_\xi$ particles in the system. The angle brackets denote the ensemble average; in practice, this ensemble average is usually substituted by the time average in steady-state MD systems, including EMD ones,²¹,²³ and a moving time average in non-equilibrium MD.³¹ On the other hand, for a simple two-body interaction between the particles, the three terms of the interaction contribution in the RHS of Eq. (13) are separated into the following:

$$ \tau^{\text{int},\xi\xi} = \begin{cases} -\frac{1}{V^{\text{CV}}} \left\langle \sum_{i,j(>i)}^{N_{\xi}-1} \sum_{\xi}^{N_{\xi}} w_{ij}^{\text{CV}} \boldsymbol{r}_{ij} \otimes \boldsymbol{f}_{ij} \right\rangle, & (\xi = \zeta = \alpha, \beta), \\ -\frac{1}{V^{\text{CV}}} \left\langle \sum_{i \in \xi}^{N_{\xi}} \sum_{j \in \xi}^{N_{\xi}} w_{ij}^{\text{CV}} \boldsymbol{r}_{ij} \otimes \boldsymbol{f}_{ij} \right\rangle, & (\xi = \alpha, \zeta = \beta), \end{cases} \quad (17) $$

where $r_{ij} = r_j - r_i$ and $f_{ij}$ are the relative position vector and force exerted from particle i to j, whereas $w_{ij}^{\text{CV}}$ denotes the weighting function given as the length fraction of the straight line segment connecting particle i and j in the CV. A mathematically proper expression for the Cartesian coordinate system is given in Ref. 30. We obtained the stress distribution with a CV thickness $\delta x$ of 1.24 × 10⁻¹ nm and time-averaging over 200 ns.
Indeed, as naturally expected from equilibrium momentum balance, Eq. (15) was satisfied for all systems with varying $η$ values, and based on this result, the interfacial tension $γ_{αβ}$ was obtained by Eq. (10) using the stress distributions obtained by the VA approach.

### C. Thermodynamic route

As a thermodynamic approach, Leroy et al.¹⁰,¹¹,³² proposed to calculate the SL and SV interfacial tensions using the thermodynamic integration (TI) method. Generally, TI is a method to calculate the free energy difference between two different equilibrium systems by connecting them by a thermodynamically reversible quasi-static route.³³⁻³⁵ For example, one introduces a coupling parameter $λ$ into the Hamiltonian of the system $H(Γ, λ)$, with the phase variable $Γ$ composed of the positions and momenta of all constituent particles. Let $G(λ)$ be the Gibbs free energy of the system as a function of $λ$. Then, the difference between a target system at $λ = λ_1$ and a reference system at $λ = λ_0$ is obtained as follows:

$$G(λ_1) - G(λ_0) = \int_{λ_0}^{λ_1} \left\langle \frac{\partial H(λ')}{\partial λ'} \right\rangle dλ', \quad (18)$$

where $\left\langle \frac{\partial H(λ')}{\partial λ'} \right\rangle$ denotes the equilibrium ensemble average of $\frac{\partial H(Γ, λ')}{\partial λ'}$ in the phase space of $Γ$, which is substituted by the time average in the EMD systems in this study.³⁶ By embedding $λ$ into the Hamiltonian in the right-hand side (RHS) so that $H(Γ, λ)$ can analytically be differentiable by $λ$, the free energy difference in Eq. (18) can be calculated by numerically integrating $\left\langle \frac{\partial H}{\partial λ'} \right\rangle$ obtained in each equilibrium system with a discrete $λ'$ value between $λ_0$ and $λ_1$.

Two implementations of the TI are used: the phantom wall (PW) method¹⁰,³² and the dry surface (DS) method.¹¹ This provides a comparison of the two approaches as well as ensuring the TI is performed consistently. Conceptually, the PW works like a pair of nets pulled through the liquid, each catching only one particle type to separate the two liquids. Meanwhile, the DS slowly changes the interaction of the two liquids, encouraging them to separate.

In the PW method, a wall is introduced that interacts only with the liquid particles through a short-range repulsive potential function. This is applied to a quasi-1D system with a flat SL interface. The PW is set parallel to the interface, and the liquid is expelled by quasi-statically moving PW starting from the solid side and moving to the liquid side under constant $NpT$ conditions. In this method, the PW position is linked with the coupling parameter $λ$ in the system Hamiltonian, and $\frac{\partial H}{\partial λ'}$ corresponds to the force exerted by the PW on the system, i.e., the quasi-static work exerted on the system is calculated by the integral in Eq. (18). This thermodynamic minimum work corresponds to the free energy difference between a system with a target SL interface at $λ = λ_1$ and a reference system with a solid surface exposed to vacuum and a PW-liquid interface achieved at $λ = λ_0$.¹⁶,³⁷

On the other hand, for the DS method, the parameter $λ$ is embedded into the SL interaction potential energy and the free energy difference of the target system relative to the reference system with a “dry” solid surface, in which the SL interaction is very weak or almost repulsive.¹³,¹⁵,¹⁹ In this method, $\left\langle \frac{\partial H}{\partial λ'} \right\rangle$ in Eq. (18) corresponds to the total SL interaction energy of the system.

Although the two methods give similar results, the DS method allows parameterization of the SL interfacial tension as a function

of the liquid interaction parameter with a lower computational cost than the PW method.¹³,¹⁵ Indeed, with the PW method, one needs to carry out a SL stripping process (by moving the PW to quasi-statically expel the liquid from the solid as mentioned above) for each SL interaction parameter; however, the PW method is simple and, therefore, applicable to charged systems, including long-range interactions. In addition, the PW method gives a direct intuitive link with the thermodynamic minimum work.

In this study, we used the DS and PW methods, both extended for the liquid-liquid interface. In both methods, we evaluated the thermodynamic minimum work needed to change from a target system to a reference system with $α$ and $β$ completely isolated without mixing, where the contribution to the free energy difference vanished. Note that the implementations of the reference systems were different for the DS and PW methods as in Figs. 1(b-i) and 1(b-ii); however, they are assumed equivalent from a thermodynamic point of view.

The details of the two methods are described in Appendix A. The basic point is that the thermodynamic equilibrium state of the present $NpT$ constant system is determined by giving two variables of miscibility $η$ and the phantom wall position $x^{PW}$, which correspond to the positions of the symmetric semipermeable PWs at $±x^{PW}$, and we change only one as the coupling parameter for the TI while keeping the other unchanged. Note that the average piston positions $\langle x_{\text{left}}^{\text{pc}} \rangle$ and $\langle x_{\text{right}}^{\text{pc}} \rangle$ are dependent variables determined by $(η, x^{PW})$ through the control pressure $p$. In the DS method, we set the miscibility parameter $η$ in Eq. (4) as the coupling parameter for the TI while keeping $x^{PW} = x_1^{PW}$ unchanged and reproduced the reference system by setting $η → 0$ and keeping $η$ positive as in the original DS method.¹¹ In the present case, we set the system at $η = η_0$ (= 0.01) as the reference system with completely isolated liquids as shown in Fig. 1(b-i). The value of $x_1^{PW}$ was set large enough that the PWs were located far away from the liquid (behind the pressure control pistons) and did not interact with the liquid particles, i.e., the PWs had no contribution to the system Hamiltonian for the present DS procedure. Then, based on Eq. (18), the free energy difference between the target system at $(η, x^{PW}) = (η, x_1^{PW})$ and the reference system at $(η_0, x_1^{PW})$ is calculated by

$$G(η, x_1^{\text{PW}}) - G(η_0, x_1^{\text{PW}}) = \int_{η_0}^{η} \left\langle \frac{\partial H(η', x_1^{\text{PW}})}{\partial η'} \right\rangle dη', \quad (19)$$

where the integrand $\left\langle \frac{\partial H(η', x_1^{PW})}{\partial η'} \right\rangle$ in the RHS can be easily obtained in the MD system as in the original DS procedure. According to the second law of thermodynamics, this corresponds to the sum of the minimum work needed for the change from the target system to the reference system and the quasi-static work exerted on the pressure control pistons as the environment under constant $NpT$. By separating the work on the PCs per unit area given by $ΔW_{\text{pc}}[(η, x_1^{PW}) → (η_0, x_1^{PW})]$ [see Eq. (A2) for the definition], we define $W_{\text{iso}}^{\text{DS}}(η)$ as the “work of isolation” calculated by the DS method as

$$W_{\text{iso}}^{\text{DS}}(η) = -\frac{G(η, x_1^{\text{PW}}) - G(η_0, x_1^{\text{PW}})}{A} - \Delta W_{\text{pc}}[(η, x_1^{\text{PW}}) → (η_0, x_1^{\text{PW}})]. \quad (20)$$

In this study, we prepared multiple equilibrium systems with different miscibility parameters $η \in [0.01, 1]$ and calculated the time
average of the average LL potential energy over 20 ns for each equilibrium system to numerically integrate the RHS of Eq. (19) [see Eq. (A7) in Appendix A]. With this procedure, the dependence of $W_{iso}^{DS}(\eta)$ on $\eta$ was clearly observed through the numerical integration with respect to $\eta$.

On the other hand, in the PW method, we used $x^{PW}$ as the coupling parameter while keeping $\eta$ unchanged and reproduced the reference system by pushing the PWs as in Fig. 1(b-ii), i.e., changing $x^{PW}$ from $x_1^{PW}$ to $x_0^{PW}$ so that the two liquids were completely isolated. In this case, the free energy difference of $G(\eta, x^{PW})$ is given by

$$G(\eta, x_1^{PW}) - G(\eta, x_0^{PW}) = \int_{x_0^{PW}}^{x_1^{PW}} \left\langle \frac{\partial \mathcal{H}(\eta, x^{PW})}{\partial x^{PW'}} \right\rangle dx^{PW'}. \quad (21)$$

It is obvious from the PW-liquid potential functions in Eqs. (7)–(9) that the Hamiltonian derivative is reduced to the forces on the PWs as

$$\left\langle \frac{\partial \mathcal{H}(\eta, x^{PW})}{\partial x^{PW}} \right\rangle = -\left\langle F_{\text{liquid-pw}\alpha}(\eta, x^{PW}) \right\rangle + \left\langle F_{\text{liquid-pw}\beta}(\eta, x^{PW}) \right\rangle, \quad (22)$$

where $\langle F_{\text{liquid-pw}\alpha}(\eta, x^{PW}) \rangle$ and $\langle F_{\text{liquid-pw}\beta}(\eta, x^{PW}) \rangle$ denote the forces on the phantom-walls $\alpha$ and $\beta$ from the corresponding liquids, respectively. Note that

$$\langle F_{\text{liquid-pw}\alpha} \rangle \approx -\langle F_{\text{liquid-pw}\beta} \rangle \ge 0 \quad (23)$$

for the forces because of the repulsive setting. Hence, the RHS of Eq. (21) corresponds to the work exerted by the two PWs (see Appendix A for details). Similar to Eq. (20), we define $W_{iso}^{PW}(\eta)$ as the work of isolation calculated by the PW method as

$$W_{iso}^{PW}(\eta) = -\frac{G(\eta, x_1^{PW}) - G(\eta, x_0^{PW})}{A} - \Delta W_{pc}[(\eta, x_1^{PW}) \rightarrow (\eta, x_0^{PW})]. \quad (24)$$

Note that the PW position $x^{PW}$ was varied while keeping $\eta$ constant. In practice, we calculated $W_{iso}^{PW}(\eta)$ for several discrete $\eta$ values. This was in contrast to the calculation of $W_{iso}^{DS}(\eta)$, which was calculated through the numerical integration with respect to $\eta$.

We compared $W_{iso}^{PW}(\eta)$, $W_{iso}^{DS}(\eta)$, and the interfacial tension $\gamma_{\alpha\beta}(\eta)$ for various $\eta$ values. We validated that the two TI approaches were carried out quasi-statically by tracing equilibrium points for

<figure><img src="image_12.png" /><figcaption>FIG. 2. (a) Distributions of the diagonal stress components $\tau_{yy}$ (black) and $\tau_{xx}$ (green) calculated by the VA and the densities $\rho_{\alpha}$ (blue) and $\rho_{\beta}$ (red) of $\alpha$ and $\beta$ components for the systems at $\eta = 0.01, 0.4, 0.8$, and 1. Enlarged snapshots of the systems around the interface are also shown. (b) Interfacial tension $\gamma_{\alpha\beta}$ calculated from the stress distribution by Eq. (10) as the mechanical route. The value of $2\gamma_{LV}$ obtained in an independent system is also displayed.</figcaption></figure>
both reversible paths through the comparison between $W_{iso}^{PW}(\eta)$ and $W_{iso}^{DS}(\eta)$. This will be described in Sec. III B with Fig. 3.

### III. RESULTS AND DISCUSSION

#### A. Stress distribution and resulting interfacial tension

Figure 2(a) shows the distributions of the diagonal stress components $\tau_{yy}$ and $\tau_{xx}$ calculated by the VA and the density distributions of each component for the systems at $\eta = 0.01, 0.4, 0.8$, and 1. Enlarged snapshots of the systems around the interface are shown above the distributions. As observed in the density distributions and snapshots, two liquids completely isolated at $\eta = 0.01$ were gradually mixed with the increase of $\eta$, and one component dissolved into the other at $\eta = 0.8$ where, for instance, $\rho_{\alpha}$ was non-zero even at the right boundary away from the interface. The two liquids were completely mixed at $\eta = 1$, and the density was homogeneous without forming an interface because the two particles were identical. Regarding the stress distributions, by using the stress definition described in Sec. II B, the uniformity of $\tau_{xx}$ in Eq. (15), which was consistent with the mechanical equilibrium condition in the direction normal to the interface, was satisfied for all $\eta$ values even at the interfaces with a steep spatial change in the densities. In addition, the constant value $\tau_{xx}$ was equal to $\tau_{yy}$ in the bulk regions away from the interface, i.e., the stress isotropy was satisfied there. The isotropic value $\tau_{xx} = \tau_{yy}$) in bulk is negative because it is equal to the bulk pressure with its sign reverted. Based on these results, we obtained the LL interfacial tension $\gamma_{\alpha\beta}(\eta)$ by Eq. (10) as a function of $\eta$ using the present definition. This corresponds to the integral for the regions filled with light red in Fig. 2(a).

Two distinct regions with $\tau_{yy} > \tau_{xx}$ separated at the boundary of two liquids existed for $\eta = 0.01(=\eta_0)$ [Fig. 2(a-i)], where the two liquids were isolated and the sum of the densities $\rho_{\alpha} + \rho_{\beta}$ was almost zero at the boundary. The interface is equivalent to the system with two interfaces between the liquid (L) and vapor (V), each with a surface tension of $\gamma_{LV}$, i.e.,

$$ \gamma_{\alpha\beta}(\eta_0) \approx 2\gamma_{LV}. \qquad (25) $$

<figure><img src="image_13.png" /><figcaption>FIG. 3. Comparison between the relative interfacial tension $-(\gamma_{\alpha\beta} - 2\gamma_{LV})$ obtained with the mechanical route and work of isolation $W_{iso}$ for various $\eta$. Error bars for $W_{iso}^{DS}$ are not shown for better visualization here (see Fig. 5). Side snapshots in Fig. 1(b) with the corresponding $(\eta, x^{PW})$ value are appended for some systems.</figcaption></figure>

This can also be intuitively understood from the distribution of $\tau_{yy}$ with two peaks displayed in Fig. 2(a-i). Indeed, we independently performed an EMD simulation of a quasi-1D single component LV system at coexistence with two flat LV interfaces at the same temperature and calculated the value of $\gamma_{LV}^{13}$ and confirmed that the resulting $2\gamma_{LV}$ was consistent with the value of $\gamma_{\alpha\beta}(\eta_0)$ at 22.96 $\pm 0.03 \times 10^{-3}$ N/m in the present study as indicated in Eq. (25). Note that the transverse stress $\tau_{yy}$ reached very large positive values $\sim +10^7$ Pa, as estimated by Rowlinson and Widom,³⁸ in comparison to the bulk value around $-10^6$ Pa.

The two isolated peaks of $\tau_{yy}$ were merged with the increase of $\eta$, and the maximum values of $\tau_{yy}$ became smaller; thus, the resulting integral was also reduced [Figs. 2(a-ii) and 2(a-iii)]. At $\eta = 1$, the two stress components $\tau_{yy}$ and $\tau_{xx}$ were equal, and the corresponding integral in Eq. (10) was zero [Fig. 2(a-iv)]. This physically means that no interface exists and the mechanical interfacial tension satisfies

$$ \gamma_{\alpha\beta}(\eta = 1) = 0. \qquad (26) $$

In principle, the LL interfacial tension $\gamma_{\alpha\beta}$ ranges between 0 and $2\gamma_{LV}$, corresponding to completely miscible and immiscible cases, respectively. Figure 2(b) shows $\gamma_{\alpha\beta}(\eta)$ calculated by Eq. (10) for each value of $\eta$. As expected, $\gamma_{\alpha\beta}(\eta)$ monotonically decreased with the increase of $\eta$, and it reached zero at around $\eta = 0.85$. Above this threshold value, the two liquids were mixed, and no LL interface was formed, as shown in Fig. 1(b-i). Note that this threshold value of $\eta$, above which the two components are completely mixed, is similar to the critical point for a single component system. However, different from the critical point, this threshold value of $\eta$ depends on the temperature. We used $-[\gamma_{\alpha\beta}(\eta) - 2\gamma_{LV}]$ as the difference from $\gamma_{\alpha\beta}(\eta_0)$ indicated in the right vertical axis for the comparison with the works of isolation below.

#### B. Comparison of interfacial tension and works of isolation

We compared the mechanical interfacial tension $\gamma_{\alpha\beta}(\eta)$ and the works for isolation $W_{iso}^{DS}(\eta)$ and $W_{iso}^{PW}(\eta)$ obtained by the DS
and extended-PW methods, respectively, for various $\eta$. Prior to that, we set the standard basis for the comparison through the difference of $\gamma_{\alpha\beta}(\eta)$ and $W_{iso}(\eta)$, considering two representative cases of $\eta = \eta_0$ and $\eta = 1$, where the values could be evaluated from the physical meanings. As shown in Fig. 2, $\gamma_{\alpha\beta}(\eta_0) = 2\gamma_{LV}$ and $\gamma_{\alpha\beta}(1) = 0$ in Eqs. (25) and (26), respectively, hold for $\gamma_{\alpha\beta}(\eta)$. Meanwhile, for the works of isolation, it is obvious that

$$W_{iso}^{DS}(\eta_0) = 0, \quad \text{and} \quad W_{iso}^{PW}(\eta_0) = 0, \qquad (27)$$

where the second equality is because no work is given to the system for further separation of the two isolated liquids by the PWs. On the other hand, for the case of $\eta = 1$, where $\alpha$ and $\beta$ are identical, a certain positive work $W_{iso}^{PW}(\eta = 1)$ is needed to isolate $\alpha$ and $\beta$ from this mixed state, although the correspondence with the mechanical interfacial tension $\gamma_{\alpha\beta}(1)$ is not clear. Considering these points, we compared $-[\gamma_{\alpha\beta}(\eta) - 2\gamma_{LV}]$ using the values obtained by the mechanical route with $W_{iso}(\eta)$ obtained by the two thermodynamic routes. Note that the value $2\gamma_{LV}$ corresponds to the work per unit area required to divide a single component liquid into two isolated parts with a depletion region equivalent to the gas phase between them as in the reference system. In addition, note that the unit N/m for the mechanical route is equivalent to the unit J/m² for the thermodynamic route.

Figure 3 displays the comparison among relative interfacial tension $-[\gamma_{\alpha\beta}(\eta) - 2\gamma_{LV}]$ and the works of isolation $W_{iso}^{DS}(\eta)$ and $W_{iso}^{PW}(\eta)$ as a function of the miscibility parameter $\eta$. Note that the error bars for $W_{iso}^{DS}$ are not shown in this figure for better visualization; however, they are relatively small as shown in Fig. 5 in Appendix B. We start from the comparison between $W_{iso}^{DS}(\eta)$ and $W_{iso}^{PW}(\eta)$. As described in Subsection II C, $W_{iso}^{DS}(\eta)$ was obtained as a quasi-smooth function by the extended-DS method. On the other hand, $W_{iso}^{PW}(\eta)$ was obtained for discrete $\eta$ values. Nevertheless, it was shown that the two values $W_{iso}^{DS}(\eta)$ and $W_{iso}^{PW}(\eta)$ gave the same result, although the thermodynamic integration paths were toward equivalent but different reference systems at $(\eta, x_{1}^{PW}) = (\eta_0, x_{1}^{PW})$ and $(\eta, x_{0}^{PW})$, respectively. This match also indicates that the works of isolation $W_{iso}$ were correctly obtained by the two TI methods by quasi-statically tracing the equilibrium thermodynamic points along the reversible TI paths. Note that the error bars were larger for the PW method because the force on the PWs in Eq. (22) was used, which was subject to larger thermal fluctuations (see Fig. 4).

Regarding the comparison between $-[\gamma_{\alpha\beta}(\eta) - 2\gamma_{LV}]$ obtained by the mechanical route and $W_{iso}(\eta)$ by the thermodynamic routes, they matched well for small $\eta$ values; however, with the increase of $\eta$ from about $\eta = 0.7$, the difference became large and it steeply increased above $\eta = 0.8$. For instance, the difference for the completely mixing case at $\eta = 1$ was about $193 \times 10^{-3}$ N/m, which was much larger than $-(\gamma_{\alpha\beta} - 2\gamma_{LV})$ of about $23 \times 10^{-3}$ N/m. Briefly, the difference was because extra work is needed to separate the components in the bulk liquid regions in addition to the work needed to completely detach the two liquids into two parts, as schematically illustrated in the right-hand side of Fig. 3. This point will be further examined in the next Subsection III C.

<figure><img src="image_14.png" /><figcaption>(b) Schematic of the force balance among PW-force and pressures ($x^{PW} = 3.5$ nm)</figcaption></figure>

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

<figure><img src="image_16.png" /><figcaption>FIG. 4. (a) Force on the right phantom wall (PW) per unit area upon the calculation of the work of isolation $W_{iso}^{PW}(\eta)$ by the extended-PW method in the completely miscible case ($\eta = 1$). The inset corresponds to Fig. 1(a-ii). (b) Schematic of the force balance among the force on the left and right PWs, the pressure of the two single component liquids, and that of the mixed liquid in the center between the two PWs.</figcaption></figure>

### C. Decomposition of the work of isolation for the completely miscible case

To examine the difference between the mechanical interfacial tension and the thermodynamic work of isolation observed in Fig. 3, we looked into the intermediate process of the PW method in the completely miscible case ($\eta = 1$), focusing on the force on the PWs. Figure 4(a) shows the force on the right PW per unit area $\frac{\langle F_{liquid-pwa} \rangle}{A}$ upon the calculation of the work of isolation $W_{iso}^{PW}(\eta)$ by the PW method at $\eta = 1$, where the inset corresponds to Fig. 1(a-ii). Note that the $\alpha-\beta$ interface was not formed in the $\alpha-\beta$ mixture because the two liquids were identical in this system at $\eta = 1$. As the PW
entered the liquid with the increase of $x^{\text{PW}}$, almost constant force was exerted on the PW as observed at $x^{\text{PW}} = 5, 3.5$, and 2 nm, where the liquid separation by the semipermeable PWs proceeded with reducing the volume of the mixed liquid. The force $\frac{F_{\text{liquid-pwa}}}{A}$ was noticeably larger than the control pressure at these states. This force corresponds to the osmotic pressure,³⁹,⁴⁰ as will be discussed in Sub-section III D. The force steeply rose up just before the liquid was completely isolated at $x^{\text{PW}} = 0.5$ nm, and it decreased down to the control pressure at $x^{\text{PW}} = 0.2$ nm. No work was needed to further separate the liquids where the work done by the PWs on the system and that on the PCs were balanced.

As indicated with the arrows and snapshots in Fig. 3, the TI was started from $(\eta, x^{\text{PW}}) = (1, x_1^{\text{PW}})$ to equivalent different systems, i.e., to $(\eta_0, x_1^{\text{PW}})$ for the DS and to $(\eta, x_0^{\text{PW}})$ for the PW methods. Indeed, along the path of the PW method, the liquids were gradually separated until the intermediate state had a LL interface without sandwiching vapor, as shown in the snapshot at $x^{\text{PW}} = 0.5$ nm in the inset of Fig. 4. From this intermediate state, the two liquids were completely detached by the PWs to achieve the reference system having two isolated LV interfaces. Now, for the case at $\eta = 1$ where $\alpha$ and $\beta$ were identical, we assume $\gamma_{\alpha\beta}$ to be zero at the intermediate state even though two liquids are separated by the semipermeable PWs. Then, letting $x_{\text{sep}}^{\text{PW}}(\eta)$ be the PW position of the intermediate state, the minimum work $\Delta W[(\eta, x_{\text{sep}}^{\text{PW}}(\eta)) \rightarrow (\eta, x_0^{\text{PW}})]|_{\eta=1}$ needed for the change from this intermediate system at $(\eta, x^{\text{PW}}) = (1, x_{\text{sep}}^{\text{PW}}(1))$ to the reference system at $(1, x_0^{\text{PW}})$ is equal to $2\gamma_{\text{LV}}$, i.e.,

$$ \Delta W[(\eta, x_{\text{sep}}^{\text{PW}}(\eta)) \rightarrow (\eta, x_0^{\text{PW}})]|_{\eta=1} = -[\gamma_{\alpha\beta}(\eta = 1) - 2\gamma_{\text{LV}}] = 2\gamma_{\text{LV}}. \quad (28) $$

Note that $x_{\text{sep}}^{\text{PW}}(\eta)|_{\eta=1}$ was about 0.5 nm as indicated in Fig. 4, but in general, $x_{\text{sep}}^{\text{PW}}(\eta)$ is given as a function of $\eta$ because the position slightly depends on the mixing feature governed by $\eta$. We define this as the “detachment work” $W_{\text{det}}(\eta)$ as

$$ W_{\text{det}}(\eta) \equiv \Delta W[(\eta, x_{\text{sep}}^{\text{PW}}(\eta)) \rightarrow (\eta, x_0^{\text{PW}})], \quad (29) $$

which satisfies

$$ W_{\text{det}}(1) = 2\gamma_{\text{LV}} \quad (30) $$

for a specific case at $\eta = 1$. In addition, we introduce the minimum work from the target system $(\eta, x_1^{\text{PW}})$ to the intermediate state at $(\eta, x_{\text{sep}}^{\text{PW}}(\eta))$ to separate the liquids defined by

$$ W_{\text{sep}}(\eta) \equiv \Delta W[(\eta, x_1^{\text{PW}}) \rightarrow (\eta, x_{\text{sep}}^{\text{PW}}(\eta))] = W_{\text{iso}}(\eta) - W_{\text{det}}(\eta), \quad (31) $$

which we call the “separation work.”

The value of $W_{\text{sep}}(1)$ is estimated for a specific case at $\eta = 1$ from a viewpoint of configurational entropy here. Indeed, $\eta = 1$ corresponds to an ideal mixture for which the mixing free energy is purely of entropic origin. In particular, since $\alpha$ and $\beta$ were identical under the constant temperature $T$ and pressure $p$ in this case, the two states should have the same liquid structure with the same internal energy $U$ and volume $V$ even though liquids $\alpha$ and $\beta$ were separated on the left and right sides, respectively. In other words, the difference between the two states is that particles of $\alpha$ and $\beta$ can move in the whole space between the two PCs in the mix-state, whereas each kind of particle can move only within half of the space. Let $G$ and $S$ be the Gibbs free energy and entropy, respectively, and let “mix” and “sep” be the subscripts for completely mixed target and intermediate

states, respectively. Then, it follows for $W_{\text{sep}}(1)$ that

$$ \begin{aligned} W_{\text{sep}}(1) &= \frac{1}{A} (G_{\text{sep}} - G_{\text{mix}})|_{\eta=1} \\ &= \frac{1}{A} [(U_{\text{sep}} + pV_{\text{sep}} - TS_{\text{sep}}) - (U_{\text{mix}} + pV_{\text{mix}} - TS_{\text{mix}})]|_{\eta=1} \\ &= \frac{T}{A} (S_{\text{mix}} - S_{\text{sep}})|_{\eta=1} \end{aligned} \quad (32) $$

because $U_{\text{mix}} = U_{\text{sep}}$ and $V_{\text{mix}} = V_{\text{sep}}$ are assumed. The entropy difference can be estimated by the possible volume available for liquids $\alpha$ and $\beta$ composed of $N_{\alpha}$ and $N_{\beta}$ particles, respectively, as in a thought experiment of ideal gas separation by using a semipermeable membrane in ordinary thermodynamics and statistical mechanics textbooks⁴¹ as follows:

$$ \begin{aligned} (S_{\text{mix}} - S_{\text{sep}})|_{\eta=1} &\equiv (S_{\text{mix}}^{\alpha} - S_{\text{sep}}^{\alpha})|_{\eta=1} + (S_{\text{mix}}^{\beta} - S_{\text{sep}}^{\beta})|_{\eta=1} \\ &= N_{\alpha} k_{\text{B}} \ln \left. \frac{V_{\text{mix}}}{V_{\text{sep}}^{\alpha}} \right|_{\eta=1} + N_{\beta} k_{\text{B}} \ln \left. \frac{V_{\text{mix}}}{V_{\text{sep}}^{\beta}} \right|_{\eta=1} \\ &= 2N_{\alpha} k_{\text{B}} \ln \left. \frac{V_{\text{mix}}}{V_{\text{sep}}^{\alpha}} \right|_{\eta=1}, \end{aligned} \quad (33) $$

where $k_B$ is the Boltzmann constant; and $V_{\text{mix}}$, $V_{\text{sep}}^{\alpha}$, and $V_{\text{sep}}^{\beta}$ are the volumes of the liquid in the mixed state and those for $\alpha$ and $\beta$ parts in the intermediate state. Note that $N_{\alpha} = N_{\beta}$ and $V_{\alpha} = V_{\beta}$ were used for the final equality. The resulting $W_{\text{sep}}(1)$ obtained by Eq. (31) using the DS result of $W_{\text{iso}}^{\text{DS}}(1)$ and $W_{\text{det}}(1) = 2\gamma_{\text{LV}}$ was

$$ W_{\text{sep}}(1) = W_{\text{iso}}^{\text{DS}}(1) - 2\gamma_{\text{LV}} = (193.0 \pm 1.0) \times 10^{-3} \text{ J/m}^2. \quad (34) $$

On the other hand, the entropy difference estimated by Eqs. (32) and (33) was

$$ \frac{2N_{\alpha}k_{\text{B}}T}{A} \ln \left. \frac{V_{\text{mix}}}{V_{\text{sep}}^{\alpha}} \right|_{\eta=1} = (185.7 \pm 9.6) \times 10^{-3} \text{ J/m}^2, \quad (35) $$

where the PW position and average positions of the PCs for the system at $(\eta, x^{\text{PW}}) = (1, 0.5)$ nm were used to roughly estimate the volume $V_{\text{sep}}^{\alpha}$. Indeed, the two agreed well, and this indicated that the work of isolation included the interfacial tension and the mixing free energy. Note that the method can also be used for non-ideal mixtures, i.e., for the case of $\eta \neq 1$, to extract both the interfacial free energy change and the free energy of mixing, which will then contain both enthalpic and entropic contributions.

## D. Osmotic pressure

We discuss the meaning of the constant force per unit area $\frac{F_{\text{liquid-pwa}}}{A}$, which was larger than the control pressure observed during the liquid separation in Fig. 4(a). Figure 4(b) illustrates the schematic of the force balance at $x^{\text{PW}} = 3.5$ nm. The pressure was controlled at $p_{\text{set}}$ by the pressure control pistons on both ends of the system. This indicated that the two single component liquids on the left and right, both between a piston and a phantom wall (PW), had a pressure of $p_{\text{set}}$. On the other hand, the mixed liquid between the PWs in the center was subject to the pressure of the single component liquids as well as the PWs. Hence, the constant force per unit
area shown with the dotted purple line in Fig. 4(a) corresponded to the osmotic pressure, as discussed in previous work.39,40

## IV. CONCLUSION

We performed molecular dynamics simulations of a liquid-liquid interface between two different Lennard-Jones liquids with various miscibilities and evaluated the interfacial tension using both mechanical and thermodynamic routes. In the case of the mechanical route, the local stress normal to the interface, computed through the volume average method, was observed to be constant over the entire region. From the stress distributions, we calculated the liquid-liquid interfacial tension obtained by using Bakker's equation for various miscibilities and compared it with the free energy obtained by the thermodynamic integration, where the extended dry-surface and phantom-wall schemes were used to quasi-statically isolate the two liquids under constant pressure and temperature conditions. When the two components were immiscible, the mechanical and thermodynamic results were in good agreement, whereas when they were miscible, the values were significantly different. This difference was attributed to the additional free energy required to separate the binary liquid into two components, i.e., the free energy of mixing. In the phantom wall setup, it was possible to disentangle the free energy of mixing, which corresponded to the work of the osmotic pressure acting on the phantom wall prior to the complete detachment of the two components and the change in interfacial free energy occurring upon detachment. In the ideal mixture case, we showed that the free energy of mixing corresponded to the entropy difference between the mixed state and the separated state. For non-ideal mixtures, the PW method provides the full free energy of mixing—including an enthalpic and an entropic contribution, together with the osmotic pressure of the mixtures.

## ACKNOWLEDGMENTS

We thank Haruki Oga for the fruitful discussion as a former member of Y.Y.'s group at Osaka University. H.K., T.O., and Y.Y. were supported by the JSPS KAKENHI grant (Nos. JP23KJ0090, JP23H01346, and JP22H01400), Japan, respectively. Y.Y. was also supported by the JST CREST grant (No. JPMJCR18I1), Japan. We discussed this research during the NEMD Conference held on 13th and 14th June 2024, partly supported by JSPS Bilateral Joint Research Seminars (Japan–UK, No. 220249903).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Rei Ogawa:** Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Project administration (equal); Validation (lead); Visualization (lead); Writing – original draft (equal); Writing – review & editing (supporting).

**Hiroki Kusudo:** Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Project administration (supporting); Resources (supporting); Software (supporting); Supervision (supporting); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal).

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

## DATA AVAILABILITY

Based on Eq. (18), the free energy difference between the target system at $(\eta, x_1^{\text{pw}}) = (\eta_0, x_1^{\text{pw}})$ and the reference system at $(\eta_0, x_1^{\text{pw}})$ is calculated by

$$G(\eta, x_1^{\text{pw}}) - G(\eta_0, x_1^{\text{pw}}) = \int_{\eta_0}^{\eta} \left( \frac{\partial \mathcal{H}(\eta', x_1^{\text{pw}})}{\partial \eta'} \right) d\eta'. \quad (A1)$$

Let the minimum work exerted on the PCs and that needed for the change from the target system both per unit area be defined by $\Delta W_{\text{pc}}[(\eta, x_1^{\text{pw}}) \to (\eta_0, x_1^{\text{pw}})]$ and $\Delta W[(\eta, x_1^{\text{pw}}) \to (\eta_0, x_1^{\text{pw}})]$, respectively, it follows for the former that

$$\begin{aligned} \Delta W_{\text{pc}}[(\eta, x_1^{\text{pw}}) \to (\eta_0, x_1^{\text{pw}})] \\ = p_{\text{set}}[\langle x_{\text{right}}^{\text{pc}}(\eta_0, x_1^{\text{pw}}) \rangle - \langle x_{\text{right}}^{\text{pc}}(\eta, x_1^{\text{pw}}) \rangle] \\ + p_{\text{set}}[\langle x_{\text{left}}^{\text{pc}}(\eta, x_1^{\text{pw}}) \rangle - \langle x_{\text{left}}^{\text{pc}}(\eta_0, x_1^{\text{pw}}) \rangle]. \end{aligned} \quad (A2)$$

Hence, the latter can be obtained by

$$\Delta W[(\eta, x_1^{\text{pw}}) \to (\eta_0, x_1^{\text{pw}})] = - \frac{G(\eta, x_1^{\text{pw}}) - G(\eta_0, x_1^{\text{pw}})}{A} - \Delta W_{\text{pc}}[(\eta, x_1^{\text{pw}}) \to (\eta_0, x_1^{\text{pw}})]. \quad (A3)$$

We define the “work of isolation” by the DS denoted by $W_{\text{iso}}^{\text{DS}}(\eta)$ as this difference in this study, i.e.,

$$W_{\text{iso}}^{\text{DS}}(\eta) \equiv \Delta W[(\eta, x_1^{\text{pw}}) \to (\eta_0, x_1^{\text{pw}})]. \quad (A4)$$

Equation (A3) with the definition in Eq. (A4) corresponds to Eq. (20).
Regarding Eq. (A1), let $\Phi_{a\beta}(\Gamma, \eta)$ be the sum of the potential energy between the different liquids in Eq. (4),

$$ \Phi_{a\beta}(\Gamma, \eta, x_1^{\text{pw}}) = \sum_{i \in \alpha} \sum_{j \in \beta} \phi_{a\beta}(r_{ij}) = \sum_{i \in \alpha} \sum_{j \in \beta} \eta \phi_{aa}(r_{ij}), \quad (A5) $$

then the integrand $\left(\frac{\partial \mathcal{H}(\eta', x_1^{\text{pw}})}{\partial \eta'}\right)$ in Eq. (19) is given by

$$ \begin{aligned} \left\langle \frac{\partial \mathcal{H}(\eta', x_1^{\text{pw}})}{\partial \eta'} \right\rangle &= \left\langle \frac{\partial \Phi_{a\beta}(\eta', x_1^{\text{pw}})}{\partial \eta'} \right\rangle = \left\langle \frac{\partial}{\partial \eta'} \sum_{i \in \alpha} \sum_{j \in \beta} \phi_{a\beta}(r_{ij}) \right\rangle \\ &= \left\langle \sum_{i \in \alpha} \sum_{j \in \beta} \phi_{a\alpha}(r_{ij}) \right\rangle = \frac{1}{\eta'} \left\langle \sum_{i \in \alpha} \sum_{j \in \beta} \phi_{a\beta}(r_{ij}) \right\rangle \\ &= \frac{\langle \Phi_{a\beta}(\eta', x_1^{\text{pw}}) \rangle}{\eta'} \end{aligned} \quad (A6) $$

where $\langle \Phi_{a\beta}(\eta', x_1^{\text{pw}}) \rangle$ can be easily obtained in the MD simulation. By inserting Eq. (A6) into Eq. (19) and further inserting it into Eq. (A3), the work of isolation results in

$$ W_{\text{iso}}^{\text{DS}}(\eta) = \int_{\eta_0}^{\eta} \frac{1}{\eta'} \left[ -\frac{\langle \Phi_{a\beta}(\eta', x_1^{\text{pw}}) \rangle}{A} \right] d\eta' \\ - \Delta W_{\text{pc}} [(\eta, x_1^{\text{pw}}) \rightarrow (\eta_0, x_1^{\text{pw}})], \quad (A7) $$

where $\langle \Phi_{a\beta}(\eta', x_1^{\text{pw}}) \rangle_A$ is the average LL potential energy per area with its sign reverted.

In this study, we prepared multiple equilibrium systems with different values of the miscibility parameter $\eta \in [0.01, 1]$ and calculated the time average of the average LL potential energy over 20 ns for each equilibrium system to numerically integrate the first term of the RHS of Eq. (A7).

On the other hand, in the PW method, we used the PW position $x^{\text{PW}}$ as the coupling parameter and reproduced the reference system by pushing the PWs as in Fig. 1(b-ii), i.e., decreasing $x^{\text{PW}}$ from $x_1^{\text{PW}}$ down to $x_0^{\text{PW}}$ so that the two liquids were completely isolated. For this work, we define the work of isolation by the extended-PW denoted by $W_{\text{iso}}^{\text{PW}}$ given by

$$ W_{\text{iso}}^{\text{PW}}(\eta) \equiv \Delta W \left[ (\eta, x_1^{\text{PW}}) \rightarrow (\eta, x_0^{\text{PW}}) \right]. \quad (A8) $$

In this case, the difference of G is written by

$$ G(\eta, x_1^{\text{pw}}) - G(\eta, x_0^{\text{pw}}) = \int_{x_0^{\text{pw}}}^{x_1^{\text{pw}}} \left\langle \frac{\partial \mathcal{H}(\eta, x^{\text{PW}})}{\partial x^{\text{PW}}} \right\rangle dx^{\text{PW}}. \quad (A9) $$

It is obvious from the PW-liquid potential functions in Eqs. (7)–(9) that the Hamiltonian derivative is reduced to the forces on the PWs as

$$ \begin{aligned} \left\langle \frac{\partial \mathcal{H}(\eta, x^{\text{PW}})}{\partial x^{\text{PW}}} \right\rangle &= \langle -F_{\text{liquid-pw}\alpha}(\eta, x^{\text{PW}}) + F_{\text{liquid-pw}\beta}(\eta, x^{\text{PW}}) \rangle \\ &= \langle -F_{\text{liquid-pw}\alpha}(\eta, x^{\text{PW}}) \rangle + \langle F_{\text{liquid-pw}\beta}(\eta, x^{\text{PW}}) \rangle, \end{aligned} \quad (A10) $$

where the average forces from the liquid $\alpha$ on PW$\alpha$ (right) and $\beta$ on PW$\beta$ (left) are defined by

$$ \langle F_{\text{liquid-pw}\alpha}(\eta, x^{\text{PW}}) \rangle = \left\langle -\sum_{i \in \alpha} \frac{\partial \phi_{\text{liquid-pw}\alpha}(x_i')}{\partial x_{\alpha}^{\text{PW}}} \right\rangle \quad (A11) $$

and

$$ \langle F_{\text{liquid-pw}\beta}(\eta, x^{\text{PW}}) \rangle = \left\langle \sum_{i \in \beta} \frac{\partial \phi_{\text{liquid-pw}\beta}(x_i')}{\partial x_{\beta}^{\text{PW}}} \right\rangle, \quad (A12) $$

respectively. By using Eqs. (A10) and (A3), the work of isolation for the present systems is obtained by

$$ \begin{aligned} W_{\text{iso}}^{\text{PW}}(\eta) &= -\frac{G(\eta, x_1^{\text{PW}}) - G(\eta, x_0^{\text{PW}})}{A} - \Delta W_{\text{pc}} [(\eta, x_1^{\text{PW}}) \rightarrow (\eta, x_0^{\text{PW}})] \\ &= -\int_{x_1^{\text{PW}}}^{x_0^{\text{PW}}} \frac{\langle F_{\text{liquid-pw}\alpha}(\eta, x^{\text{PW}}) \rangle}{A} dx^{\text{PW}} \\ &\quad + \int_{x_1^{\text{PW}}}^{x_0^{\text{PW}}} \frac{\langle F_{\text{liquid-pw}\beta}(\eta, x^{\text{PW}}) \rangle}{A} dx^{\text{PW}} \\ &\quad - \Delta W_{\text{pc}} [(\eta, x_1^{\text{PW}}) \rightarrow (\eta, x_0^{\text{PW}})]. \end{aligned} \quad (A13) $$

Note $x_1^{\text{PW}} > x_0^{\text{PW}}$, and also note that

$$ \frac{\langle F_{\text{liquid-pw}\alpha} \rangle}{A} \approx - \frac{\langle F_{\text{liquid-pw}\beta} \rangle}{A} \geq 0 \quad (A14) $$

holds for the forces from the property of $\phi_{\text{liquid-pw}}$ and symmetry.

## APPENDIX B: WORK OF ISOLATION BY THE EXTENDED-DS METHOD

Figure 5(a) shows the work of isolation $W_{\text{iso}}^{\text{DS}}(\eta)$ obtained by the DS method in Eq. (A7) as a function of the miscibility parameter $\eta$, where the average potential energy between $\alpha$ and $\beta$ per

<figure><img src="image_17.png" /><figcaption>FIG. 5. (a) Work of isolation $W_{\text{iso}}^{\text{DS}}$ obtained by the DS method in Eq. (A7), and (b) average potential energy between different components $\alpha$ and $\beta$ for various miscibility parameters $\eta$. The inset corresponds to Fig. 1(a-i). The error bars for (b) were smaller than the size of the symbol.</figcaption></figure>
area $-\frac{(\Phi_{ap}(\eta x_1^{PW}))}{A}$ in Fig. 5(b) was used. Note that the PWs did not interact with the liquids at the PW position $x^{PW} = x_1^{PW}$, and the PWs are not shown in the inset. The potential energy $\Phi_{\alpha\beta}$ was always negative and was almost zero for small $\eta$ values because the two liquids were not mixed, as observed in the enlarged snapshot at $\eta = 0.4$ in Fig. 2. Hence, the resulting $W_{iso}$ monotonically increased with a small gradient up to about $\eta = 0.4$. For $\eta > 0.5$, $-\frac{(\Phi_{ap})}{A}$ considerably increased, and the increase became steep, especially above around $\eta = 0.7$, and the resulting $W_{iso}^{DS}(\eta)$ also showed a large increase around $0.7 < \eta < 0.85$. Finally, above around $\eta = 0.85$, the steep increase of $-\frac{(\Phi_{ap})}{A}$ calmed down. This turning point value of $\eta = 0.85$ matched the point above which $\gamma_{\alpha\beta}$ became zero in Fig. 2.

## REFERENCES

¹J. G. Kirkwood, *J. Chem. Phys.* **3**, 300 (1935).

²J. G. Kirkwood and F. P. Buff, *J. Chem. Phys.* **19**, 774 (1951).

³J. G. Kirkwood and F. P. Buff, *J. Chem. Phys.* **17**, 338 (1949).

⁴G. Bakker, *Kapillarität und Oberflächenspannung* (Wien-Harms, 1928), Vol. 6.

⁵M. P. Allen and D. J. Tildesley, *Computer Simulation of Liquids* (Oxford University Press, 1987).

⁶M. Hayoun, M. Meyer, M. Mareschal, G. Ciccotti, and P. Turq, “Molecular dynamics simulation of a liquid-liquid interface,” in *Chemical Reactivity in Liquids: Fundamental Aspects*, edited by M. Moreau and P. Turq (Springer US, Boston, MA, 1988), pp. 279–286.

⁷I. Benjamin, *Annu. Rev. Phys. Chem.* **48**, 407 (1997).

⁸E. Feria, J. Algaba, J. M. Míguez, A. Mejia, and F. J. Blas, *Phys. Chem. Chem. Phys.* **24**, 5371 (2022).

⁹B. D. Todd, D. J. Evans, and P. J. Daivis, *Phys. Rev. E* **52**, 1627 (1995).

¹⁰F. Leroy, D. J. V. A. Dos Santos, and F. Müller-Plathe, *Macromol. Rapid Commun.* **30**, 864 (2009).

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

¹²M. Kanduč, *J. Chem. Phys.* **147**, 174701 (2017).

¹³Y. Yamaguchi, H. Kusudo, D. Surblys, T. Omori, and G. Kikugawa, *J. Chem. Phys.* **150**, 044701 (2019).

¹⁴A. Russo, M. A. Durán-Olivencia, S. Kalliadasis, and R. Hartkamp, *J. Chem. Phys.* **150**, 214705 (2019).

¹⁵D. Surblys, F. Leroy, Y. Yamaguchi, and F. Müller-Plathe, *J. Chem. Phys.* **148**, 134707 (2018).

¹⁶C. Bistafa, D. Surblys, H. Kusudo, and Y. Yamaguchi, *J. Chem. Phys.* **155**, 064703 (2021).

¹⁷M. Shintaku, H. Oga, H. Kusudo, E. R. Smith, T. Omori, and Y. Yamaguchi, *J. Chem. Phys.* **160**, 224502 (2024).

¹⁸S. Nishida, D. Surblys, Y. Yamaguchi, K. Kuroda, M. Kagawa, T. Nakajima, and H. Fujimura, *J. Chem. Phys.* **140**, 074707 (2014).

¹⁹H. Kusudo, T. Omori, and Y. Yamaguchi, *J. Chem. Phys.* **151**, 154501 (2019).

²⁰Y. Imaizumi, T. Omori, H. Kusudo, C. Bistafa, and Y. Yamaguchi, *J. Chem. Phys.* **153**, 034701 (2020).

²¹H. Kusudo, T. Omori, and Y. Yamaguchi, *J. Chem. Phys.* **155**, 184103 (2021).

²²K. Watanabe, H. Kusudo, C. Bistafa, T. Omori, and Y. Yamaguchi, *J. Chem. Phys.* **156**, 054701 (2022).

²³H. Kusudo, T. Omori, L. Joly, and Y. Yamaguchi, *J. Chem. Phys.* **159**, 161102 (2023).

²⁴J. P. R. B. Walton, D. J. Tildesley, J. S. Rowlinson, and J. R. Henderson, *Mol. Phys.* **48**, 1357 (1983).

²⁵M. Sega, B. Fábián, G. Horvai, and P. Jedlovszky, *J. Phys. Chem. C* **120**, 27468 (2016).

²⁶G. Hantal, B. Fábián, M. Sega, and P. Jedlovszky, *J. Mol. Liq.* **306**, 112872 (2020).

²⁷T. Omori and Y. Yamaguchi, *J. Chem. Phys.* **161**, 204704 (2024).

²⁸S. Plimpton, *J. Comput. Phys.* **117**, 1 (1995).

²⁹B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, *J. Chem. Theory Comput.* **4**, 435 (2008).

³⁰K. Shi, E. R. Smith, E. E. Santiso, and K. E. Gubbins, *J. Chem. Phys.* **158**, 040901 (2023).

³¹B. D. Todd and P. J. Daivis, *Nonequilibrium Molecular Dynamics: Theory, Algorithms and Applications* (Cambridge University Press, 2017).

³²F. Leroy and F. Müller-Plathe, *J. Chem. Phys.* **133**, 044110 (2010).

³³D. Frenkel and B. Smit, in *Understanding Molecular Simulation: From Algorithms to Applications*, 2nd ed., *Computational Science Series Vol. 1* (Academic Press, 2001), pp. 168–172.

³⁴R. L. Davidchack and B. B. Laird, *J. Chem. Phys.* **118**, 7651 (2003).

³⁵N. di Pasquale and R. L. Davidchack, *J. Phys. Chem. A* **126**, 2134 (2022).

³⁶Y. Yamaguchi, H. Kusudo, C. Bistafa, D. Surblys, T. Omori, and G. Kikugawa, *ECS Trans.* **108**, 93 (2022).

³⁷D. Surblys, Y. Yamaguchi, K. Kuroda, M. Kagawa, T. Nakajima, and H. Fujimura, *J. Chem. Phys.* **140**, 034505 (2014).

³⁸J. S. Rowlinson and B. Widom, *Molecular Theory of Capillarity* (Dover, 1982), p. 45.

³⁹Y. Luo and B. Roux, *J. Phys. Chem. Lett.* **1**, 183 (2010).

⁴⁰M. Bley, M. Duvail, P. Guilbaud, and J.-F. Dufrêche, *J. Phys. Chem. B* **121**, 9647 (2017).

⁴¹W. Greiner, L. Neise, H. Stöcker, and D. Rischke, *Thermodynamics and Statistical Mechanics (Classical Theoretical Physics)*, 3rd ed. (Springer, 1995), pp. 43–48.