# Statistical modeling of equilibrium phase transition in confined fluids

Gunjan Auti<EMAIL_ADDRESS>Department of Mechanical
Engineering, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan  Soumyadeep Paul  Department of
Mechanical Engineering, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan  Department of Mechanical
Engineering, Stanford University, Building 530, 440 Escondido Mall, Stanford,
California 94305, USA  Wei-Lun Hsu  Department of Mechanical Engineering, The
University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan  Shohei Chiashi  Department of
Mechanical Engineering, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan  Shigeo Maruyama  Department of
Mechanical Engineering, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan  Hirofumi Daiguji
<EMAIL_ADDRESS>Department of Mechanical Engineering, The
University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan

###### Abstract

The phase transition of confined fluids in mesoporous materials deviates from
that of bulk fluids due to the interactions with the surrounding heterogeneous
structure. For example, adsorbed fluids in metal-organic-frameworks (MOFs)
have atypical phase characteristics such as capillary condensation and higher-
order phase transitions due to a strong heterogeneous field. Considering a
many-body problem in the presence of a nonuniform external field, we model the
host-guest and guest-guest interactions in MOFs. To solve the three-
dimensional Ising model, we use the mean-field theory to approximate the
guest-guest interactions and Mayer’s $f$-functions to describe the host-guest
interactions in a unit cell. Later, using Hill’s theory of nanothermodynamics,
we define differential thermodynamic functions to understand the distribution
of intensive properties and integral thermodynamic functions to explain the
phase transition in confined fluids. The investigation reveals a distinct
behavior where fluids confined in larger pores undergo a discontinuous (first-
order) phase transition, whereas those confined in smaller pores experience a
continuous (higher-order) phase transition. Furthermore, the results indicate
that the free-energy barrier for phase transitions is lower in confined fluids
than in bulk fluids giving rise to a lower condensation pressure relative to
the bulk saturation pressure. Finally, the integral thermodynamic functions
are succinctly presented in the form of a phase diagram, marking an initial
step toward a more practical approach for understanding the phase behavior of
confined fluids.

††preprint: APS/123-QED

##  I Introduction

Confined fluids can be gas (vapor) trapped in nanobubbles [1], adsorbed gas in
porous structures [2, 3, 4], natural gas trapped in shale and tight rock
formations [5], or biomolecules trapped in cells [6]. Fluids such as these
enclosed within restricted spaces have distinct physical and thermodynamic
characteristics that differ from those of bulk fluids. These unique
characteristics include phenomena such as an atypical phase transition and
packing polymorphism [7, 8, 9, 10], a shift in the freezing and melting points
[11], an anomalously low dielectric constant [12], and very high hydrodynamic
slippage [13]. Although these phenomena have been widely observed and
reported, how heterogeneity and their multiscale nature affect the fluid
characteristics still lacks a comprehensive thermodynamic understanding.

The atypical thermodynamic properties are due to the heterogeneous
interactions and steric hindrance of confined fluids. Generally, these
interactions take the form of van der Waals forces and the inverse radial
dependence of these cohesive interactions translates into a layered
distribution of density near the surface, which creates anisotropy [14, 15,
16, 17, 18]. Gibbs laid the foundation for modeling how heterogeneous
interactions affect fluid properties by formulating surface thermodynamics,
where he introduced the concept of the Gibbs surface excess [19]. Later, Hill
put forth a thermodynamic approach to model small systems [20]. Since then,
numerous models describing how heterogeneity affects the fluid properties have
been proposed. For example, various self-consistent field models have been
presented based on the analogous Bogoliubov–Born–Green–Kirkwood–Yvon hierarchy
describing a system containing a large number of interacting particles [21,
22, 23, 24, 25, 26, 27, 28, 29]. These models can be solved for special cases
but become much more complex for a different sets of conditions.

To make the model more realizable under general conditions, Sircar and Mayers
[2, 30] proposed a semi-empirical formulation for low-concentration
adsorption. They coined the term “ideal adsorbed phase” to indicate a behavior
that differs from that of bulk fluids. Nicholson [31, 32] later presented an
elaborate molecular theory focusing on the adsorption of lattice gas. This
model shows qualitatively how the thermodynamic properties vary but does not
focus on the phase transition of the adsorbed fluid. Martinez et al. [33]
predicted adsorption isotherms using a two-dimensional statistical associating
fluid theory for a square well potential on a flat surface.

In parallel, with the advancement of molecular simulations, Evans [34, 3]
undertook extensive Monte Carlo simulations to elucidate phase transitions in
mesoporous slits. Schmidt and Löwen [35] introduced a hard-sphere model to
investigate the freezing transition between parallel plates, presenting a
phase diagram through Monte Carlo simulations. Subsequently, numerous studies
used brute force molecular simulations to produce both qualitative and
quantitative phase diagrams. For instance, Kimura and Maruyama [36] explored
boiling phase transitions and cluster formation using molecular dynamics
simulations, while Radhakrishnan et al. [37] used biased potentials and
umbrella sampling in grand canonical Monte Carlo (GCMC) simulations to obtain
thermodynamic properties. Takaiwa et al. [38] presented a phase diagram of
water in carbon nanotubes. Zhou et al. [39] conducted ab initio simulations
using density functional theory and grand canonical algorithms to illustrate a
surface phase diagram.

More recently, these molecular simulation techniques have been improved by
incorporating various machine-learning algorithms [40, 41, 42, 43]. Such
methodologies aim to extract specific properties of adsorbed fluids, reducing
computational demands and facilitating the prediction of certain thermodynamic
properties. However, even the most current deep-learning models are black
boxes [44, 45], meaning that they perform the majority of their calculations
internally, so that significant thermodynamic information is overlooked. This
hinders a comprehensive understanding of the physics underlying adsorption and
confinement. Moreover, deep-learning algorithms are essentially the methods of
statistical mechanics, as explained by Lin et al. [46]. Consequently,
understanding the thermodynamics of confined fluids requires statistical
methods and analytical models based on appropriate approximations.

Statistical models have been used to understand various multiscale processes.
For example, Košmrlj and Nelson [47] gave a model for the thin shells and
argued that large spherical shells are unstable due to thermally generated
pressure using statistical mechanics. Goodrich et al. [48] formulated a
statistical model for nanocluster formation in the crystallization process,
and Molina et al. [49] experimentally described the many-body interactions
that occur in confined space for self-organizing of droplets. In a similar
way, to understand the multiscale process of phase transition in confined
fluids, we have developed a semi-analytical statistical model.

Figure 1: Statistical model for adsorbed fluid. (a) Schematic representation
of the system under consideration. The argon molecules (red) in bulk with
chemical potential $\mu^{\text{bulk}}$ are in equilibrium with the system of
adsorbed molecules with chemical potential $\mu^{\text{ads}}$ in a metal-
organic framework (MOF) with chemical potential $\mu^{\text{frame}}$. The
relative density distribution obtained from the proposed model for an MOF with
a 12 Å unit cell is shown at (b) low relative pressure ($p/p_{0}=0.04$) and
(c) after saturation ($p/p_{0}=0.20$). Similar distribution functions are
plotted for a 24 Å unit cell MOF at (d) low relative pressure ($p/p_{0}=0.10$)
and (e) after saturation ($p/p_{0}=0.70$). The deep blue shade represents the
heterogeneity (metal or ligand) and all distributions are plotted for a cross
section of the unit cell at $z=a/2$.

The present paper introduces a three-dimensional (3D) Ising model for argon
confined in a cubic metal-organic framework [MOF, see Fig. 1(a)]. This
approach considers the nonuniformity of the external field by decoupling
homogeneous and heterogeneous interactions. The homogeneous interactions are
considered through mean-field theory and the heterogeneous interactions are
approximated using Mayer’s $f$-functions. The nonuniform interactions lead to
a nonuniform density distribution in the pores, as depicted in Figs.
1(b)–1(e). To account for the thermodynamic properties of such a distribution,
differential (local) and integral (global) intensive thermodynamic functions
can be defined. For example, Figs. 1(b) and 1(c) show the expected relative
density distribution (a differential property) of argon at relative pressures,
$p/p_{0}=0.04$ and $p/p_{0}=0.20$, respectively, for a small pore ($a=12$ Å).
Even for extremely low relative pressure, the pores are close to saturation.
Conversely, Figs. 1(d) and 1(e) show the expected relative density
distribution of argon at relative pressures $p/p_{0}=0.10$ and $p/p_{0}=0.70$,
respectively, for a large pore ($a=24$ Å). Here, the layered adsorption occurs
near the heterogeneity at lower relative pressure and a uniform density
distribution develops at higher relative pressure. Later in this paper, we
derive the integral thermodynamic functions using Hill’s thermodynamics for
small systems. Based on these integral properties, we discuss the phase
transition of confined fluids and compare it with that of bulk fluids.

This paper discusses the statistical modeling of the confined fluid while
highlighting the phase transition during capillary condensation. We also
showcase the phase diagram of the adsorbed fluid and discuss the similarities
and key differences vis-á-vis the bulk fluid. The remainder of the paper is
organized as follows: Section II.2 details a mathematical derivation of the
Ising model in a grand canonical ensemble for a confined fluid. Section II.3
derives the relevant thermodynamic properties of the adsorbed fluids. Section
III.1 explains the benchmarking of the proposed model with GCMC simulations.
Section III.2 discusses the phase transition as a function of pore size based
on a double-well potential. Sec. III.3 introduces the phase diagram for the
adsorbed fluid. Section IV.1 then discusses how this model may be used in
different contexts (e.g., where the confinement effect is significant).
Finally, Sec. IV.2 presents the conclusions of this paper.

##  II Model

###  II.1 General assumptions

We focus on the classical regime, where quantum effects may be disregarded.
Consequently, the van der Waals potentials generated by different sources are
treated as additive, following the established principles of intermolecular
forces [50]. Furthermore, the analysis assumes equilibrium conditions
throughout. This assumption is based on the premise that external conditions,
such as temperature and pressure, are time-independent. Given such conditions,
macroscopic quantities can be expressed in terms of microscopic average
values, distribution functions, or probabilities.

Furthermore, mean-field theory is used to solve the many-body problem of
adsorbents in a potential well of nonuniform depth. This assumption is valid
under conditions of low adsorbent concentration because interactions between
adsorbed molecules are negligible at such concentrations. Moreover, the mean-
field theory remains valid at and above the saturation point for capillary
condensation because the distribution of fluid molecules within pores becomes
uniform, leading to a mean-field effect. However, the mean-field theory may
not be accurate in a high-density gas-like regime. In such cases where the
adsorbents are relatively concentrated, a density distribution around the
heterogeneity creates an anisotropy. However, this effect is not considered in
the current model.

Given these assumptions, the current work provides a framework for analyzing
and understanding the behavior of fluids in confined spaces, particularly in
the context of adsorption in MOFs. These assumptions allow for simplified
models and calculations, enabling insights into the thermodynamic properties
and phase transitions of the adsorbed fluids.

To clarify the formulation of the model, we briefly revisit the fundamental
concepts of Hill’s nanothermodynamics [51] in the context of the current
problem (a detailed derivation and discussion are available in Ref. [52]). In
our case, fluid confined in the pore of a modeled MOF is in equilibrium with
the surrounding bulk fluid. The argon molecules confined within the MOF are
the “system” in this investigation. To understand the thermodynamic
characteristics of this system, we subdivide it into an ensemble of $\eta$
small, equivalent, distinguishable, independent systems, as shown in Fig.
1(a). Therefore, assuming the total volume is constant, the system described
here at equilibrium gives

$\displaystyle dE_{t}=TdS_{t}+\mu dN_{t}+\xi d\eta,$ (1)

where $E_{t}$ is the total energy of the system, $T$ is the temperature,
$S_{t}$ is the total entropy of the system, $\mu$ is the chemical potential,
$N_{t}$ are the total number of molecules in the system, $\xi$ is the
subdivision potential, and $\eta$ is the number of subdivisions.

Equation (1) resembles Gibb’s equilibrium equation for a two-component system
with $N_{t}$ being the number of molecules. Since $\eta$ is the number of
subdivisions then the total volume $V_{t}=\eta V$, where $V$ is the volume of
each subdivision. Therefore, we consider the work associated with varying
$\eta$ at pressure $p$ by adding the work of expansion, $-p\eta dV$, in Eq.
(1) to obtain

$\displaystyle dE_{t}=TdS_{t}-p\eta dV+\mu dN_{t}+\xi d\eta,$ (2)

where $-p\eta\equiv\partial E_{t}/\partial V$. We now use Hill’s definition of
subdivision potential, $\xi\equiv-\hat{p}V$, where $\hat{p}$ is the integral
pressure [52]. Integrating Eq. (2) gives the equilibrium equation for the
total system:

$\displaystyle E_{t}=TS_{t}+\mu N_{t}-\hat{p}V\eta.$ (3)

From here, it is straightforward to show for a grand canonical ensemble that

$\displaystyle\hat{p}V=k_{\text{B}}T\ln\Xi,$ (4)

where $k_{\text{B}}$ is the Boltzmann constant and $\Xi$ is the grand
partition function for the small chosen system of volume $V$.

Hereinafter, we use Hill’s notation [52] where the hat ($\hat{~{}}$) denotes
the integral intensive thermodynamic function and for any extensive
thermodynamic function $\alpha$, and

$\displaystyle\bar{\alpha}\equiv\frac{1}{V}\int_{V}\alpha dV$ (5)

denotes the integral extensive thermodynamic function. In contrast, symbols
without hat or bar represent differential thermodynamic functions.

###  II.2 Ising framework

To understand the phase transition, the Ising model plays a crucial role. The
inherent complexity in the 3D Ising model coupled with the external non-
uniform field presents significant challenges. However, we have made certain
assumptions, mentioned in the prior section to obtain an approximate solution.
The formulation for the confined fluids is as follows:

Let $\mathbf{p}=(p_{1},p_{2},\dots,p_{N})$,
$\mathbf{q}=(q_{1},q_{2},\dots,q_{N})$ be the momentum and position
coordinates, respectively, in phase space for a system of $N$ molecules
confined in a framework creating a potential $U_{ma}(\mathbf{q})$. The
Hamiltonian $\mathcal{H}$ is then

$\displaystyle\mathcal{H}(\mathbf{p},\mathbf{q})$
$\displaystyle=\mathcal{K}(\mathbf{p},\mathbf{q})+\mathcal{U}(\mathbf{p},\mathbf{q})$
(6a)
$\displaystyle=\sum_{i}^{N}\frac{p^{2}_{i}}{2m}+\sum_{i}^{N-1}\sum_{j>i}^{N}\Phi(\mathbf{r}_{i}-\mathbf{r}_{j})+\sum_{i}^{N}U_{ma}(q_{i}),$
(6b)

where $\mathcal{K}$ and $\mathcal{U}$ are kinetic and potential energy
contributions to the Hamiltonian, respectively. $m$ is mass of the particle,
$\Phi(\mathbf{r}_{i}-\mathbf{r}_{j})$ is the potential between molecules $i$
and $j$ as a function of the spatial coordinate $\mathbf{r}$. This
intermolecular potential can be approximated as a field defined in terms of
the phase-space coordinate $U_{aa}(\mathbf{q})$:

$\displaystyle\sum_{i}^{N}U_{aa}(q_{i})\equiv\sum_{i}^{N-1}\sum_{j>i}^{N}\Phi(\mathbf{r}_{i}-\mathbf{r}_{j}).$
(7)

####  II.2.1 Canonical ensemble

Equation (6b) shows that the kinetic-energy term is independent of the
position coordinate $\mathbf{q}$ and the potential-energy term is independent
of the momentum coordinate $\mathbf{p}$. Assuming that the confining framework
is stationary and only fluid particles contribute to the kinetic energy, we
separate the variables and define the canonical partition function
$\mathbb{Z}$ as the product of the kinetic contribution $Z_{k-aa}$ and the
configurational contribution $Z_{q}$ [refer Eq. 8]. Additionally, the kinetic
energy in turn depends only on the temperature of the reservoir. Therefore, we
focus on the solution of the configurational partition.

$\displaystyle\mathbb{Z}=Z_{k-aa}Z_{q},$ (8)

where the configurational contribution is

$\displaystyle
Z_{q}=\frac{1}{V^{N}}\int_{V}e^{-\mathcal{U}(q_{1},q_{2},\dots,q_{N})/\tau}dq_{1}\dots
dq_{N},$ (9)

with $\tau=k_{\text{B}}T$. The potential energy is a combination of potentials
created by adsorbate-adsorbate interaction and the MOF-adsorbate interaction
as described in Eq. 10

$\displaystyle\mathcal{U}(\mathbf{q})=U_{aa}(\mathbf{q})+U_{ma}(\mathbf{q}).$
(10)

$U_{ma}$ is the potential energy due to the framework at any given position,
which implies

$\displaystyle U_{ma}(\mathbf{q})$
$\displaystyle=U_{ma}(q_{1})+U_{ma}(q_{2})+\dots+U_{ma}(q_{N}).$ (11)

The configurational partition function [Eq. (9)] is

$\displaystyle Z_{q}=\frac{1}{V^{N}}\int_{V}$
$\displaystyle\exp\\{-[U_{ma}(q_{1})+\cdots+U_{ma}(q_{N})$ $\displaystyle+$
$\displaystyle U_{aa}(q_{1})+\cdots+U_{aa}(q_{N})]/\tau\\}dq_{1}\cdots
dq_{N}.$ (12)

For the extreme case in the absence of any external field, where
$U_{ma}(q_{1},\dots,q_{N})=0$, Eq. (12) resembles the configurational
partition function for a bulk fluid.

The complexity due to the inclusion of an external field may be treated in
several ways. For example, Travalloni et al. [53] used the extension of the
generalized van der Waals theory to model the confined fluids. Simon et al.
[54, 55] took as a lattice model the adsorbed gas and the different
orientations of the flexible framework and defined a transfer matrix for a
rather complex problem. Poluektov [28] analyzed the self-consistent field
model for classical systems using a one-dimensional perturbation theory. Singh
et al. [56] proposed decoupling the two types of interactions and
approximating the solution. Recently, Dong et al. [57] used Gibbs-surface
thermodynamics to define the problem in appropriate independent variables and
obtained an analytical solution for the special case of confinement between
parallel sheets.

We follow the approach of Singh et al. [56] of decoupling the interactions and
approximating the effect of the external potential using Mayer’s $f$-functions
[58, 59]. Let $f_{i}$ be defined as follows:

$\displaystyle f_{i}\equiv e^{-U_{ma}(q_{i})/\tau}-1.$ (13)

Equation (12) can then be written as

$\displaystyle
Z_{q}=\frac{1}{V^{N}}\int_{V}e^{-U_{aa}(\mathbf{q})/\tau}\prod_{i=1}^{N}(1+f_{i})dq_{1}\cdots
dq_{N},$ (14)

where

$\displaystyle\prod_{i=1}^{N}(1+f_{i})=1+\sum_{i=1}^{N}f_{i}+\sum_{i=1}^{N-1}\sum_{j>i}^{N}f_{i}f_{j}+{O}(f_{i}^{3}).$
(15)

Including the individual molecular terms and ignoring the higher-order terms
in Eq. (15) yields,

$\displaystyle Z_{q}=$
$\displaystyle\,\frac{1}{V^{N}}\int_{V}e^{-U_{aa}(\mathbf{q})/\tau}\left(1+\sum_{i=1}^{N}f_{i}\right)dq_{1}\cdots
dq_{N}$ (16)

This expansion is analogous to the first-order fluid-fluid interactions
described by Mayer [59]. In this formulation, Mayer’s $f$-functions are
applied to simplify the heterogeneous interactions and not the fluid-fluid
interactions.

Assuming that, despite the density distribution resulting from the external
potential, fluid particles collectively establish a mean field to render the
analytical solution tractable, we can write

$\displaystyle U_{aa}(\mathbf{q})=Nu_{aa},$ (17)

where $u_{aa}$ is the mean field independent of the positional coordinate in
phase space. Thus, $u_{aa}$ can be moved outside the integral and Eq. (16)
takes the form

$\displaystyle Z_{q}=$
$\displaystyle\,\frac{1}{V^{N}}\int_{V}e^{-Nu_{aa}(\mathbf{q})/\tau}dq_{1}\cdots
dq_{N}$
$\displaystyle+\frac{e^{-Nu_{aa}(\mathbf{q})/\tau}}{V}\left(\int_{V}f_{1}dq_{1}+\cdots+\int_{V}f_{N}dq_{N}\right).$
(18)

The first term in Eq. (18) describes the configurational partition function of
the bulk phase in the absence of any external field; let it be $Z_{q-aa}$. In
addition, we can define

$\displaystyle\phi$ $\displaystyle\equiv\frac{1}{V}\int_{V}f_{i}dq_{i}$
$\displaystyle=\frac{1}{V}\int_{V}\left(e^{-U_{ma}(q_{i})/\tau}-1\right)dq_{i}.$
(19)

Using this definition,

$\displaystyle Z_{q}\approx Z_{q-aa}\left(1+N\phi\right),$ (20)

where $Z_{q-aa}$ is the configurational partition function of the mean-field
bulk fluid. Putting this back into Eq. (8) gives

$\displaystyle\mathbb{Z}$ $\displaystyle=Z_{k-aa}Z_{q-aa}\left(1+N\phi\right)$
$\displaystyle=\mathbb{Z}_{\text{bulk}}\left(1+N\phi\right).$ (21)

Given that $U_{ma}$ is a spatially varying potential inside the mesopores, we
consider an infinitesimal volume $dV=dxdydz$ at coordinate $q_{i}=(x,y,z)$. We
assume a uniform potential $U_{ma}(x,y,z)$ in the infinitesimal volume, which
implies that the canonical partition function for an infinitesimal volume is

$\displaystyle\mathbb{Z}=\mathbb{Z}_{\text{bulk}}[1+N\phi(x,y,z)].$ (22)

####  II.2.2 Grand canonical ensemble

The equilibrium assumption implies that the chemical potential
$\mu^{\text{bulk}}$ of the bulk phase equals the chemical potential
$\mu^{\text{total}}$ of the argon inside the nanospace. The chemical potential
$\mu^{\text{total}}$ of argon inside the nanospace consists of contributions
from other argon inside nanospace ($\mu^{\text{ads}}$) and from the framework
atoms ($\mu^{\text{frame}}$) that form the heterogeneity. Both contributions
are made through intermolecular forces, not intramolecular forces:

$\displaystyle\mu^{\text{bulk}}(p_{\text{ext}},T)$
$\displaystyle=\mu^{\text{total}}$ (23a)
$\displaystyle=\mu^{\text{ads}}+\mu^{\text{frame}}.$ (23b)

Similar to the excess chemical potential defined by Widom [60], the chemical
potential of the adsorbed phase combines the intramolecular chemical potential
$\mu^{\text{ads}}$ from the adsorbed phase and the excess chemical potential
$\mu^{\text{frame}}$ due to the framework. We use the mean value of the
interaction energies in the unit cell to calculate the excess chemical
potential $\mu^{\text{frame}}$:

$\mu^{\text{frame}}=-\tau\ln\left\langle\exp\left(\frac{-U_{ma}(\mathbf{q})}{\tau}\right)\right\rangle.$
(24)

In addition, the unit-cell volume $V$ is fixed and the temperature $T$ is
controlled externally. Therefore, the modeling is done in a grand canonical
ensemble ($\mu^{\text{ads}},V,T$). Moreover, if a partition function of a
system of particles can be obtained, the relevant quantities of interest, such
as density, pressure, entropy, and free energy can be derived. To this end, we
start the model by defining the grand partition function $\Xi_{\text{ads}}$
for the adsorbed fluid [52]:

$\displaystyle\Xi_{\text{ads}}$
$\displaystyle=\sum_{N=0}^{\infty}\int\frac{d^{N}\mathbf{q}d^{N}\mathbf{p}}{h^{3N}N!}e^{-[\mathcal{H}(\mathbf{p},\mathbf{q})-\mu^{\text{ads}}N]/\tau}$
(25a) $\displaystyle=\sum_{N=0}^{\infty}\mathbb{Z}e^{N\mu^{\text{ads}}/\tau},$
(25b)

where $h$ is Planck’s constant and $\mathbb{Z}$ is the canonical partition
function for an ensemble of $N$ molecules. Combining Eqs. (22) and (25b) gives

$\displaystyle\Xi_{\text{ads}}(\mu^{\text{ads}},V,T)=\Xi_{\text{bulk}}(\mu^{\text{ads}},V,T)\left(1+\phi\langle
N\rangle^{\mu^{\text{ads}}}_{\text{bulk}}\right).$ (26)

Under extreme conditions where there is no external potential (i.e., no
confinement), the partition function $\Xi_{\text{ads}}$ simplifies to the
partition function $\Xi_{\text{bulk}}$ of the bulk fluid. This ensures the
consistency and coherence of the equation, particularly in scenarios where
confinement effects are negligible.

###  II.3 Thermodynamic properties

When examining intensive thermodynamic functions for a bulk fluid, such as
pressure or chemical potential, it is common to assume that the fluid is
homogeneous, meaning that the properties of the fluid are uniform throughout
the entire volume under consideration. However, in the context of confined
fluids, the presence of heterogeneous interactions introduces nonuniformities
in the intensive thermodynamic functions. To address this distribution of
properties, Hill introduced both differential and integral thermodynamic
functions [61, 52]. Differential thermodynamic functions are defined at a
specific point in space, indicating their local nature. Conversely, integral
thermodynamic functions extend their definition across the entire volume of
the system, providing a global characterization.

Given that the phase of any substance is defined for a group of molecules [62,
63], understanding the phase of the fluid requires that the integral
properties be considered. Therefore, in the exploration of phase transitions,
emphasis is placed on global (i.e., integral) thermodynamic properties rather
than local (i.e., differential) thermodynamic properties.

Relevant thermodynamic properties such as the grand potential
$\bar{\Omega}_{\text{ads}}$, the expected number $\overline{\langle
N_{\text{ads}}\rangle}$ of molecules adsorbed, the pressure
$\hat{p}_{\text{ads}}$ of the adsorbed phase, entropy $\bar{S}_{\text{ads}}$,
enthalpy $\bar{H}_{\text{ads}}$, Helmholtz free energy $\bar{F}_{\text{ads}}$,
and Gibbs free energy $\bar{G}_{\text{ads}}$ can all be obtained from the
grand partition function $\Xi_{\text{ads}}$, as shown in the following
section.

####  II.3.1 Grand potential $\bar{\Omega}_{\text{ads}}$

The grand potential can be obtained from its definition

$\displaystyle\Omega_{\text{ads}}(x,y,z)$
$\displaystyle=-\tau\ln[\Xi_{\text{ads}}(x,y,z)]$
$\displaystyle=-\tau\ln(\Xi_{\text{bulk}})-\tau\ln[1+\phi(x,y,z)]$
$\displaystyle=\Omega_{\text{bulk}}-\tau\ln[1+\phi(x,y,z)].$ (27)

where $\Omega_{\text{bulk}}$ is the grand potential of the bulk fluid. Since
we are interested in the integral properties inside the MOF pore, we spatially
average the grand potential:

$\displaystyle\bar{\Omega}_{\text{ads}}$
$\displaystyle=\frac{1}{V}\int_{V}\Omega_{\text{ads}}(x,y,z)dxdydz$
$\displaystyle=\hat{p}V\equiv\xi.$ (28)

####  II.3.2 Pressure $\hat{p}_{\text{ads}}$ of adsorbed phase

For any macroscopic system, we can write

$\displaystyle E_{t}-TS_{t}-\mu N_{t}=-\tau\ln\Xi_{t}.$ (29)

Combining Eqs. (3), (4), and (29), we obtain the integral pressure
$\hat{p}_{\text{ads}}$ of the adsorbed fluid in terms of the grand partition
function:

$\displaystyle\xi\equiv\hat{p}_{\text{ads}}V=\tau\ln\Xi_{\text{ads}}=-\bar{\Omega}_{\text{ads}}.$
(30)

Therefore,

$\displaystyle\hat{p}_{\text{ads}}=-\frac{(\bar{\Omega}_{\text{ads}})_{\mu,T}}{V}.$
(31)

####  II.3.3 Expected number $\overline{\langle N_{\text{ads}}\rangle}$ of
molecules adsorbed

The expected number $\overline{\langle N_{\text{ads}}\rangle}$ of molecules
adsorbed can be calculated as follows:

$\displaystyle\langle N_{\text{ads}}(x,y,z)\rangle$
$\displaystyle=\frac{\sum_{N}N\mathbb{Z}_{\text{ads}}e^{\mu^{\text{ads}}N/\tau}}{\Xi_{\text{ads}}}$
(32a) $\displaystyle=\frac{\langle
N\rangle^{\mu^{\text{ads}}}_{\text{bulk}}+\langle
N^{2}\rangle^{\mu^{\text{ads}}}_{\text{bulk}}\phi(x,y,z)}{1+\langle
N\rangle^{\mu^{\text{ads}}}_{\text{bulk}}\phi(x,y,z)},$ (32b)

where $\langle N\rangle^{\mu^{\text{ads}}}_{\text{bulk}}$ is the average
number of molecules that would be present in the bulk if the chemical
potential were $\mu^{\text{ads}}$. Note that the term $\langle
N_{\text{ads}}(x,y,z)\rangle$ is not the actual number of molecules at a given
position but the expected number. Therefore, the expected total number of
molecules in the unit cell is

$\displaystyle\overline{\langle
N_{\text{ads}}\rangle}=\frac{1}{V}\int_{V}\langle N_{\text{ads}}(x,y,z)\rangle
dxdydz.$ (33)

The number of molecules confined within a system is influenced by two factors:
intermolecular interactions and heterogeneous interactions. The maximum number
of molecules is capped by the volume of the unit cell. These interactions
collectively contribute to the effective potential experienced by the
adsorbate molecules. Equation (26) uses a decoupling approach to separate
these two interactions. The intermolecular interaction among the confined
molecules is considered independently, following which the heterogeneous
interactions are incorporated as an additional potential term. The volume
constraint is accounted for by capping the summation in Eq. (32) at
$N_{\text{max}}$ such that

$\displaystyle
N_{\text{max}}=\left\lfloor\frac{V_{\text{cell}}-V_{\text{m}}-V_{\text{l}}}{b}\right\rfloor,$
(34)

where $\lfloor\cdot\rfloor$ is the floor function, $V_{\text{cell}}$ is the
unit-cell volume, $V_{\text{m}}$ and $V_{\text{l}}$ are the volume of the unit
cell occupied by the metal and ligand, respectively,
$b=\sqrt{2}\sigma_{a}^{3}$ is the volume of each molecule according to van der
Waals theory, and $\sigma_{a}$ is the Lennard–Jones size parameter for
adsorbates.

####  II.3.4 Other thermodynamic functions

The following standard thermodynamic relations still apply to the integral
values:

$\displaystyle\text{Entropy:}~{}(\bar{S}_{\text{ads}})$
$\displaystyle=\Bigl{(}\partialderivative{\bar{\Omega}_{\text{ads}}}{T}\Bigr{)}_{p,\mu^{\text{ads}}},$
(35) $\displaystyle\text{Gibbs free energy:}~{}(\bar{G}_{\text{ads}})$
$\displaystyle=\mu^{\text{ads}}\overline{\langle N_{\text{ads}}\rangle},$ (36)
$\displaystyle\text{Helmholtz free energy:}~{}(\bar{F}_{\text{ads}})$
$\displaystyle=\bar{\Omega}_{\text{ads}}+\mu^{\text{ads}}\overline{\langle
N_{\text{ads}}\rangle},$ (37)
$\displaystyle\text{Enthalpy:}~{}(\bar{H}_{\text{ads}})$
$\displaystyle=\bar{G}_{\text{ads}}+T\bar{S}_{\text{ads}}.$ (38)

##  III Results and Discussions

###  III.1 Benchmarking

Figure 2: Benchmarking of the model with GCMC simulations. (a) Model metal-
organic framework (MOF) with metal oxides and ligands; $a$ is the unit-cell
length. (b) Potential-energy distribution inside the unit cell of this MOF.
(c) Adsorption isotherm for argon obtained from the model (green line) and the
GCMC simulations (pink diamonds). Capillary condensation occurs at relative
pressure $p/p_{\text{0}}\approx 0.6$. The insets (i)–(iii) show the
distribution of argon atoms inside the unit cell obtained from the GCMC
simulation and the expected relative density distribution (normalized with the
maximum local density for better contrast) obtained from the proposed model
presented at the corresponding locations in the isotherm. The benchmarking is
done for argon adsorption in a model MOF with unit-cell length $a=24$ Å, LJ
parameters $\sigma_{m}=5$ Å, $\varepsilon_{m}=120$ K, and temperature $T=120$
K.

We assume that bulk argon behaves as a van der Waals fluid, so we can write
the canonical partition function of the bulk fluid as

$\displaystyle\mathbb{Z}_{\text{bulk}}=\frac{1}{N!}\left(\frac{V-Nb}{\lambda^{3}_{T}}\right)^{N}\exp\left(\frac{-aN^{2}}{V\tau}\right),$
(39)

where $a$ and $b$ are the van-der Waals coefficients for argon and
$\lambda_{T}$ is the thermal de Broglie wavelength [64]. Based on this
assumption various cases with different pore sizes and temperatures were
analyzed using an in-house GPU-accelerated Python code (available upon
request).

This proposed model is benchmarked using a GCMC simulation for a model MOF
with a cubic unit cell. Metal oxides occupy the vertices of the cube and
ligands are located on the edges of the cube, as shown in Fig. 2(a). To
maintain the consistency with the statistical model, we used averaged LJ
parameters for metal-oxide and aromatic rings to model this framework [65].
The GCMC simulations were performed using the RASPA [66] simulation software
package. All simulations included a 50 000-cycle equilibration period and a
100 000-cycle production run. In these simulations, the structure of all the
frameworks is considered rigid; that is, all species of the framework are held
fixed at their crystallographic positions. The argon atoms can move in three
different ways in the GCMC simulation: translation, rotation, and swap. The
interaction between MOF and argon and between argon atoms was modeled using
the LJ potential function and the Lorentz–Berthelot mixing rule. Sample
crystallographic information (.cif) files and the LJ parameters are listed in
the electronic supporting information (ESI) [67] Sec. S1 A.

Figure 2 shows the benchmarking results for the 24 Å unit cell. Figures 2(a)
and 2(b) show the unit-cell structure and the potential-energy distribution,
respectively. Figure 2(c) compares the adsorption isotherm derived from the
GCMC simulation with that produced by the proposed model, showing that the two
curves are consistent. Additionally, Fig. 2(c) shows in the three lower panels
the molecular distribution obtained through the GCMC simulation along with the
density distribution within the pore at positions (i), (ii), and (iii). The
proposed model correctly captures the previously observed trend of layered
adsorption [14, 15, 16, 17, 18]. Specifically, at lower relative pressures,
adsorption predominantly occurs near the heterogeneity, as shown in Fig.
2(c)(i). As the relative pressure increases but before capillary condensation,
a distinct layering of adsorbed molecules is evident in Fig. 2(c)(ii).
Finally, beyond capillary condensation, the pore becomes saturated with a
density distribution resembling that of the bulk liquid [Fig. 2(c)(iii)].
Additional benchmarking results are presented in Sec. S1. Moreover, Sec. S2 of
the ESI [67] provides a brief parametric investigation that supports the
hypothesis made in our previous paper [68] that the thermodynamic properties
of confined fluids are a function of the confinement parameter
$\Psi\equiv\sigma_{ma}/a$.

Figure S1 of the ESI [67] shows that the isotherms for ultrasmall pore size
(10 Å) and larger pore sizes (24 Å) are consistent with the results obtained
from GCMC simulations. This result is attributed to the assumption of a
uniform field generated by the fluid molecules adsorbed in the cavity. In
ultrasmall pores, the variability in the field resulting from argon adsorption
at different positions is effectively equivalent to a uniform distribution.
Similarly, for larger pore sizes, the distribution of molecules near
adsorption sites remains independent at lower concentrations, thereby creating
a uniform field inside the cavity, as depicted in Fig. 2(c)(i). At higher
concentrations (that is, post capillary condensation), argon molecules densely
occupy the cavity, resulting in a uniform field, as shown in Fig. 2(c)(iii).
For concentrations between these extremes, the adsorption isotherm slightly
deviates from the GCMC values, as shown in Fig. 2(c)(ii).

Likewise, for the medium pore sizes shown in Figs. S1(b)–S1(e) of the ESI
[67], the isotherm produced by the proposed model deviates from the GCMC
isotherm. This discrepancy arises when the assumption of a uniform field in
the cavity is no longer valid. This problem could be addressed through an
iterative process, where the density distribution obtained from the proposed
model serves as the initial guess. However, this problem is beyond the scope
of present paper.

###  III.2 Phase transition in confinement

####  III.2.1 Types of phase transitions

Since the inception of fullerenes and 3D carbon nanotubes with cylindrical
pores, the phase-transitions inside these structures have been discussed [69].
Multiple studies show that freezing of water confined in these structures may
occur continuously or discontinuously [70, 71]. The results in this section
show that this transition depends on the pore size. Smaller pores produce
continuous phase transitions, whereas larger pores produce discontinuous
(first-order) phase transitions.

The Helmholtz free energy can be expressed in terms of the canonical partition
function [Eq. (22)]:

$\displaystyle F_{\text{ads}}(N,V,T)$
$\displaystyle=-\tau\ln[\mathbb{Z}_{\text{bulk}}(1+N\phi)\bigr{]}$
$\displaystyle=-N\tau\left\\{\ln\left[\frac{(V-Nb)}{\lambda^{3}_{T}N}\right]+1\right\\}-\frac{N^{2}a}{V}$
$\displaystyle-\tau\ln[1+N\phi(x,y,z)\bigr{]}.$ (40)

Therefore, we define the differential chemical potential $\mu^{\text{ads}}$
and integral chemical potential $\hat{\mu}^{\text{ads}}$ as follows:

$\displaystyle\mu^{\text{ads}}=\left(\partialderivative{F_{\text{ads}}}{N}\right)_{V,T},$
$\displaystyle\hat{\mu}^{\text{ads}}=\left(\partialderivative{\bar{F}_{\text{ads}}}{N}\right)_{V,T}=\left(\frac{\bar{F}_{\text{ads}}}{N}\right)_{V,T}.$
(41)

Figure 3: Types of phase transitions. Potential wells for fluid confined in
(a) 11 Å pores and (b) 24 Å pores. Entropy variation is plotted as a function
of relative bulk pressure confined in (c) 11 Å pores and (d) 24 Å pores at
temperatures ranging from 110 K - 150 K.

Using the formulation in the canonical ensemble, we calculate the grand
potential $\omega$ as follows:

$\displaystyle\omega=\bar{F}_{\text{ads}}-N\hat{\mu}^{\text{ads}}_{\text{sat}}.$
(42)

where $\hat{\mu}^{\text{ads}}_{\text{sat}}$ is the chemical potential at which
capillary condensation occurs for a given temperature.

Figures 3(a) and 3(b) show the grand potential thus obtained plotted as a
function of density. For the small pore size (11 Å), two phenomena occur.
First, the double well vanishes at a much lower temperature than for the large
pore size, indicating a lower critical temperature for fluids confined in
small pores. Consequently, the entropy variation in Fig. 3(c) shows that a
continuous phase transition occurs beyond 130 K. Second, the energy barrier to
cross the well is less than the thermal noise ($\Delta E_{a}\approx
0.015{k}_{\text{B}}T$ for 110 K), implying that the system spontaneously jumps
between the wells. As a result, a minute step occurs in the entropy variation
at 110 K at a lower relative pressure, as shown in Fig. 3(c). This implies
that two different phases exist. However,

$\displaystyle P_{\text{A}1}/P_{\text{A}2}=\exp(\Delta
E_{a}/k_{\text{B}}T)\approx 1,$ (43)

implying that the probability of the system being in state A1 or A3
($P_{\text{A}1}$) approximately equals the probability of the system being in
state A2 ($P_{\text{A}2}$). Therefore, these two phases are practically
indistinguishable. In addition, we hypothesize that the absence of hysteresis
during the adsorption-desorption loop is because the required activation
energy is negligible. This hypothesis is consistent with published data that
show that the Type-I adsorption isotherm for H2 adsorption in IRMOF-1
(generally observed for small pores) has no hysteresis [72].

In contrast, for the large pore size (24 Å), the barrier height is
significantly larger ($\Delta E_{a}\approx 15k_{\text{B}}T$), highlighting a
clear distinction between the “gas-like adsorbed phase” (i.e., state B1) and
the “capillary condensed phase” (i.e., state B3), as shown in Figs. 3(b) and
3(d). However, at 150 K, the double well completely vanishes, and the adsorbed
fluid exhibits a single well, which is characteristic of the supercritical
bulk fluid, suggesting that the confined fluid has a critical point.
Furthermore, based on the activation energy required for capillary
condensation, we hypothesize that hysteresis occurs during the adsorption-
desorption process for large pores. As seen for the adsorption of H2 in
MOF-253 with a Type-V adsorption isotherm (generally observed for large
pores), significant hysteresis occurs in the adsorption-desorption loop [73].

####  III.2.2 Capillary condensation

Figure 4(a) shows that the phase transition for the adsorbed fluid occurs at a
lower relative pressure than for the bulk phase transition. To elucidate this,
we compare the nucleation of droplets in a bulk fluid with that in a confined
fluid.

As with macroscale condensation, capillary condensation starts with the
heterogeneous nucleation of drop clusters that subsequently grow to form
droplets [74, 75, 76]. Consistent with classical nucleation theory [77], the
nucleation rate depends on the nucleation barrier $\Delta G^{*}$, which is the
difference between the interface free energy and the fluid free energy. A
lower nucleation barrier corresponds to a higher nucleation rate.

Figure 4: Energy barrier for capillary condensation. (a) Isotherm of relative
number density vs relative pressure for adsorbed argon and the bulk argon.
Adsorbed argon condenses at a lower pressure than bulk argon
($p/p_{\text{0}}\approx 0.6$). $N_{\text{max}}$ is the maximum number of
molecules that can be adsorbed. (b) The energy barrier for the phase
transition of bulk argon is higher than that for confined argon ($\Delta
E_{b}>\Delta E_{a}$).

Figure 4(b) compares the energy barrier $\Delta E_{b}$ for the bulk fluid with
that for the adsorbed fluid ($\Delta E_{a}$). The result shows that the free
energy for the adsorbed fluid is notably lower than that for the bulk fluid,
implying that the free-energy barrier $\Delta G^{*}$ for the confined drop
nucleation is less than that for nonconfined drop nucleation. Consequently, at
a given temperature, the condensation pressure for the adsorbed fluid is less
than that for the bulk fluid, as shown in Fig. 4(a).

###  III.3 Phase diagram of adsorbed fluid

Figure 5: Phase diagram for argon adsorbed in a metal-organic framework (MOF).
The integral properties are plotted in the form of a phase diagram for the
adsorbed fluid. Different perspectives of the (a) 3D $N\text{-}p\text{-}h$
diagram are plotted in the insets: (b) $p\text{-}h$, (c) $N\text{-}p$, and (d)
$N\text{-}h$. The red highlighted areas are the phase-coexistence region
constructed using the boundaries of discontinuities on either end of each
isotherm. Using this coexistence region, a critical point for capillary
condensation is depicted with a red dot. The isotherms are plotted every 6 K
and one isotherm in each inset is highlighted with green, to show the shape of
the isotherm. Each isotherm follows the direction
A$\rightarrow$B$\rightarrow$C$\rightarrow$D$\rightarrow$E$\rightarrow$F$\rightarrow$G.

From an application perspective, the phase diagram is an essential tool for
designing engineering processes. For example, Lilley and Prasher [78]
presented a qualitative phase diagram for crystallization of salts in an
ionocaloric refrigeration cycle. In a previous paper [79], we introduced the
concept of a 3D phase diagram as a valuable tool for the design of hybrid
compression-adsorption heat-pump systems. The statistical model discussed
earlier provides a framework for constructing such a phase diagram. Previous
research has also aimed to construct phase diagrams for confined fluids. For
example, numerous attempts have been made to construct such phase diagrams by
using Monte Carlo simulations [80, 81, 37, 82]. However, such simulations are
computationally demanding, limiting the number of adsorption isotherms that
can be generated. Consequently, acquiring the requisite thermodynamic
properties to create a phase diagram for adsorbed fluids is a significant
challenge.

Radhakrishnan et al. [37] solved this problem by applying umbrella sampling
and bias potentials to compute the system’s free energy. However, such an
approach necessitates a priori knowledge of the process, and the convergence
of their method depends heavily on the selected collective variables. Lum and
Chandler [83] addressed this issue within the framework of statistical
mechanics, albeit for a specific scenario, by deriving a phase diagram for
vapor confined within a cylindrical pore. Unfortunately, this approach
overlooks the nonuniform characteristics inherent in the external field at
this length scale, thus lacking generalizability. In contrast, Travalloni et
al. [53] constructed a phase diagram employing square-well potentials to
account for external heterogeneity. Nevertheless, this particular study does
not consider the excess chemical potential induced by external interactions,
leading to unrealistically high pressures during capillary condensation.

Figure 5(a) illustrates the 3D phase diagram of adsorbed argon within a 24 Å
model of an MOF. The three axes correspond to distinct state variables: the
external pressure is denoted $p_{\text{bulk}}$, the enthalpy per unit mass of
adsorbed argon is denoted $h_{\text{ads}}$, and the number of argon molecules
adsorbed per unit cell of the model MOF is denoted $N_{\text{ads}}$. Figures
5(b)–5(d) portray the projections of this phase diagram from three distinct
orientations, generating the $p\text{-}h$, $N\text{-}p$, and $N\text{-}h$
diagrams. Notably, the isotherms displayed contain discontinuities up to a
certain temperature, indicating that two distinct phases coexist. The
connection of these discontinuous points yields saturation lines, which are
accentuated in red. The highlighted area signifies the coexistence region
corresponding to the capillary condensed phase and the gas-like adsorbed
phase. In this phase diagram, A$\rightarrow$B adsorption occurs near the
heterogeneity-forming layered structure inside the pore, whereas
B$\rightarrow$C$\rightarrow$D$\rightarrow$E depicts the coexistence region.
E$\rightarrow$F shows that the capillary condensed liquid density increases
with pressure. Finally, beyond the bulk saturation line F$\rightarrow$G, the
pore is completely filled.

Furthermore, a detailed examination of the $p\text{-}h$ diagram in Fig. 5(b)
shows that, at low pressure, the enthalpy of the adsorbed fluid contrasts with
that of the bulk fluid. As the pressure decreases, the magnitude of the
enthalpy of the adsorbed fluid increases. This phenomenon can be understood by
considering the occurrence of layered adsorption near the metallic
heterogeneous site. The heightened cohesive interaction with the adsorption
site liberates additional energy, increasing the enthalpy at lower pressures.
ESI Fig. S8 shows the bulk argon $p\text{-}h$ diagram.

An important observation from the phase diagram is the absence of
discontinuities in the isotherms beyond a specific temperature, resembling the
behavior of bulk fluids. This temperature is denoted the critical point for
capillary condensation. Beyond this critical point, capillary condensation for
the adsorbed fluid ceases, resulting in a lack of stepwise behavior in the
adsorption isotherm. Notably, the critical point for capillary condensation is
positioned at a lower temperature compared with the bulk critical point of
argon (151 K, 48.5 bar). This difference is attributed to the excess chemical
potential of the adsorbed fluid vis-á-vis the bulk fluid, a consequence of
heterogeneous interactions. Similar results of reduction in critical pressure
of the liquid-liquid phase transition have been observed for water in a salt
solution, where the salt ions act as the heterogeneity [84, 85].

This phase diagram provides a basis to understand the phase transition of
confined fluids. Note that a critical order parameter analysis for the “gas-
like” adsorbed phase to “capillary-condesed” liquid phase and the analytical
construction of the co-existence region still remains to be addressed.

##  IV Outlook and Conclusions

###  IV.1 Outlook

This paper proposes a method to predict the integral pressure inside MOF pores
(Sec. II). Figure S3(d) of the ESI [67] shows the calculated integral
pressure. $\hat{p}_{\text{ads}}$ is significantly greater than the bulk
pressure $p$ and jumps discontinuously upon crossing the pore boundary.
Therefore, we define a disjoining pressure $\Pi_{d}$ [86, 50] such that

$\Pi_{d}\equiv\hat{p}_{\text{ads}}-p_{\text{bulk}}.$ (44)

In a similar way, when considering a constant-pressure ensemble ($NPT$), we
establish the integral chemical potential such that
$\xi\equiv\hat{\mu}^{\text{ads}}N$. Subsequently, a disjoining chemical
potential $M_{d}$ [87] is defined as

$M_{d}\equiv\hat{\mu}^{\text{ads}}-\mu^{\text{bulk}}.$ (45)

The proposed formulation can then be extended in other contexts where wall
effects and nonuniform external fields are significant, enabling us to
calculate these disjoining quantities to explain phenomena other than
adsorption.

For example, numerous attempts have been made to explain the stability of
surface nanobubbles using the Young–Laplace equation and considering external
conditions such as surface charge [88, 89] and contact-line pinning due to
hydrophilic and hydrophobic heterogeneities [90, 89]. However, a consensus
regarding the exact factors contributing to surface-bubble stability remains
elusive. The effect of surface heterogeneity clearly explains the excessive
pressure inside these bubbles in terms of a disjoining pressure, offering a
plausible explanation for their stability [91].

Similarly, experimental observations reveal that the nucleation temperature
for heterogeneous boiling is consistently lower than predicted by classical
nucleation theory [92, 93, 94]. Numerous models [95, 96, 97] and complex
simulations [98, 36, 99] have been proposed to explain this disparity between
experimental and theoretical results. Many such models consider liquid-vapor
surface tension through the Young–Laplace equation and liquid-solid surface
tension through contact-angle variation. However, to the best of our
knowledge, none of these models account for the wall effect, leading to a
disjoining pressure between vapor clusters and the liquid, thereby lowering
the free-energy barrier for phase transition, as illustrated in Fig. 4(b). The
lower free-energy barrier in turn explain the lower the nucleation temperature
of heterogeneous boiling [100].

Additionally, in the context of water transport through nanoscale osmotic
membranes, a recent study [101] argued that no chemical potential gradient
exists for transport within the membrane. However, the nanometer-scale
membrane porosity allows surface effects to trigger a disjoining chemical
potential, thereby creating a chemical-potential gradient that drives the
transport process. Considering how surface effects have a heightened impact on
the nanoscale will allow a better understanding of water transport in osmotic
membranes.

###  IV.2 Conclusions

In conclusion, our investigation offers a statistical approach to address the
3D Ising model of phase transitions in confined fluids, producing reasonably
accurate results. By applying the proposed model and integrating Hill’s theory
of nanothermodynamics, we derive both the differential (local) and integral
(global) thermodynamic properties of the adsorbed fluid. The proposed model
has practical utility for predicting the behavior of adsorbed fluids within
porous structures, facilitating the design of materials tailored to specific
requirements. The key insights derived from this model are outlined as
follows:

First, the nature of the phase transition in the confined fluid is determined
by the extent of confinement, specifically the pore size. In small pores, the
activation-energy barrier for phase transition is approximately 0.01
$k_{\text{B}}T$, significantly less than the thermal noise. Consequently, the
phase transition occurs spontaneously. Additionally, the unstable and
metastable states, while theoretically existent, are practically
indistinguishable from the stable state. Conversely, in the case of large
pores, the activation-energy barrier for phase transition is of the order of
10 $k_{\text{B}}T$, clearly distinguishing between the gas-like adsorbed phase
and the capillary condensed phase.

Second, owing to additional interactions with the surface, the free-energy
barrier for phase transitions in confined fluids is lower than in bulk fluids.
This reduced energy barrier implies that condensation inside MOF pores occurs
at a lower pressure for a given temperature, explaining the lower capillary
condensation pressure.

Finally, the model proposed in this paper consolidates the integral
thermodynamic properties in the form of a phase diagram for confined fluids.
The phase diagram resembles the bulk fluid phase diagram except for the higher
enthalpy released at lower pressure and the lower critical temperature and
pressure due heterogenous interactions.

###### Acknowledgements.

We acknowledge funding by JST, CREST (Grant No. JPMJCR17I3).

## References

  * Taverna  _et al._ [2008] D. Taverna, M. Kociak, O. Stéphan, A. Fabre, E. Finot, B. Décamps, and C. Colliex, Probing physical properties of confined fluids within individual nanobubbles, Phys. Rev. Lett. 100, 035301 (2008).
  * Sircar and Myers [1970] S. E. Sircar and A. L. Myers, Statistical thermodynamics of adsorption from liquid mixtures on solids. i. ideal adsorbed phase, J. Phys. Chem. 74, 2828 (1970).
  * Evans [1990] R. Evans, Fluids adsorbed in narrow pores: phase equilibria and structure, J. Phys. Condens. Matter. 2, 8989 (1990).
  * Morris and Wheatley [2008] R. E. Morris and P. S. Wheatley, Gas storage in nanoporous materials, Angew. Chem. Int. Ed. 47, 4966 (2008).
  * Jew  _et al._ [2022] A. D. Jew, J. L. Druhan, M. Ihme, A. R. Kovscek, I. Battiato, J. P. Kaszuba, J. R. Bargar, and G. E. Brown Jr, Chemical and reactive transport processes associated with hydraulic fracturing of unconventional oil/gas shales, Chem. Rev. 122, 9198 (2022).
  * Krishnamurthy and Cornell [2012] V. Krishnamurthy and B. Cornell, Engineering aspects of biological ion channels—from biosensors to computational models for permeation, Protoplasma 249, 3 (2012).
  * Alabarse  _et al._ [2012] F. Alabarse, J. Haines, O. Cambon, C. Levelut, D. Bourgogne, A. Haidoux, D. Granier, and B. Coasne, Freezing of water confined at the nanoscale, Phys. Rev. Lett. 109, 035701 (2012).
  * Algara-Siller  _et al._ [2015] G. Algara-Siller, O. Lehtinen, F. Wang, R. R. Nair, U. Kaiser, H. Wu, A. K. Geim, and I. V. Grigorieva, Square ice in graphene nanocapillaries, Nature 519, 443 (2015).
  * Gimondi and Salvalaglio [2018] I. Gimondi and M. Salvalaglio, $CO_{2}$ packing polymorphism under confinement in cylindrical nanopores, Mol. Syst. Des. Eng. 3, 243 (2018).
  * Kapil  _et al._ [2022] V. Kapil, C. Schran, A. Zen, J. Chen, C. J. Pickard, and A. Michaelides, The first-principles phase diagram of monolayer nanoconfined water, Nature 609, 512 (2022).
  * Hamada  _et al._ [2007] Y. Hamada, K. Koga, and H. Tanaka, Phase equilibria and interfacial tension of fluids confined in narrow pores, J. Chem. Phys. 127, 084908 (2007).
  * Fumagalli  _et al._ [2018] L. Fumagalli, A. Esfandiar, R. Fabregas, S. Hu, P. Ares, A. Janardanan, Q. Yang, B. Radha, T. Taniguchi, K. Watanabe, _et al._ , Anomalously low dielectric constant of confined water, Science 360, 1339 (2018).
  * Kavokine  _et al._ [2022] N. Kavokine, M.-L. Bocquet, and L. Bocquet, Fluctuation-induced quantum friction in nanoscale water flows, Nature 602, 84 (2022).
  * De Gennes [1981] P. G. De Gennes, Polymer solutions near an interface. adsorption and depletion layers, Macromolecules 14, 1637 (1981).
  * Fukui  _et al._ [1999] K. Fukui, H. Onishi, and Y. Iwasawa, Imaging of atomic-scale structure of oxide surfaces and adsorbed molecules by noncontact atomic force microscopy, App. Surf. Sci. 140, 259 (1999).
  * Fukuma  _et al._ [2010] T. Fukuma, Y. Ueda, S. Yoshioka, and H. Asakawa, Atomic-scale distribution of water molecules at the mica-water interface visualized by three-dimensional scanning force microscopy, Phys. Rev. Lett. 104, 016101 (2010).
  * Page  _et al._ [2014] A. J. Page, A. Elbourne, R. Stefanovic, M. A. Addicoat, G. G. Warr, K. Voïtchovsky, and R. Atkin, 3-dimensional atomic scale structure of the ionic liquid–graphite interface elucidated by AM-AFM and quantum chemical simulations, Nanoscale 6, 8100 (2014).
  * Wang and Hadjiconstantinou [2018] G. J. Wang and N. G. Hadjiconstantinou, Layered fluid structure and anomalous diffusion under nanoconfinement, Langmuir 34, 6976 (2018).
  * Gibbs [1928] J. W. Gibbs,  _The collected works of JW Gibbs_ (Longmans, Green, 1928).
  * Hill [1962] T. L. Hill, Thermodynamics of small systems, J. Chem. Phys. 36, 3182 (1962).
  * Ono and Kondo [1960] S. Ono and S. Kondo, Molecular theory of surface tension in liquids, in  _Structure of Liquids / Struktur der Flüssigkeiten_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 1960) pp. 134–280.
  * Steele and Ross [1960] W. A. Steele and M. Ross, Distribution Functions of a Fluid in an External Potential Field: Application to Physical Adsorption, J. Chem. Phys. 33, 464 (1960).
  * Wertheim [1984a] M. S. Wertheim, Fluids with highly directional attractive forces. I. Statistical thermodynamics, J. Stat. Phys. 35, 19 (1984a).
  * Wertheim [1984b] M. S. Wertheim, Fluids with highly directional attractive forces. II. Thermodynamic perturbation theory and integral equations, J. Stat. Phys. 35, 35 (1984b).
  * Wertheim [1986a] M. Wertheim, Fluids with highly directional attractive forces. III. Multiple attraction sites, J. Stat. Phys. 42, 459 (1986a).
  * Wertheim [1986b] M. Wertheim, Fluids with highly directional attractive forces. IV. Equilibrium polymerization, J. Stat. Phys. 42, 477 (1986b).
  * Biesheuvel  _et al._ [2006] P. M. Biesheuvel, F. A. Leermakers, and M. A. C. Stuart, Self-consistent field theory of protein adsorption in a non-gaussian polyelectrolyte brush, Phys. Rev. E 73, 011802 (2006).
  * Poluektov [2015] Y. M. Poluektov, Thermodynamic perturbation theory for classical systems based on self-consistent field model, Ukr. J. Phys. 60, 553 (2015).
  * Laktionov  _et al._ [2022] M. Y. Laktionov, O. V. Shavykin, F. A. Leermakers, E. B. Zhulina, and O. V. Borisov, Colloidal particles interacting with a polymer brush: A self-consistent field theory, Phys. Chem. Chem. Phys. 24, 8463 (2022).
  * Sircar [1985] S. Sircar, Excess properties and thermodynamics of multicomponent gas adsorption, J. Chem. Soc., Faraday Trans. 1 81, 1527 (1985).
  * Nicholson [1975] D. Nicholson, Molecular theory of adsorption in pore spaces. Part 1.—Isotherms for simple lattice models, J. Chem. Soc., Faraday Trans. 1 71, 238 (1975).
  * Nicholson [1976] D. Nicholson, Molecular theory of adsorption in pore spaces. part 2.—thermodynamic and molecular lattice model descriptions of capillary condensation, J. Chem. Soc., Faraday Trans. 1 72, 29 (1976).
  * Martinez  _et al._ [2007] A. Martinez, M. Castro, C. McCabe, and A. Gil-Villegas, Predicting adsorption isotherms using a two-dimensional statistical associating fluid theory, J. Chem. Phys. 126 (2007).
  * Evans  _et al._ [1986] R. Evans, U. M. B. Marconi, and P. Tarazona, Fluids in narrow pores: Adsorption, capillary condensation, and critical points, J. Chem. Phys. 84, 2376 (1986).
  * Schmidt and Löwen [1997] M. Schmidt and H. Löwen, Phase diagram of hard spheres confined between two parallel plates, Phys. Rev. E 55, 7228 (1997).
  * Kimura and Maruyama [2002] T. Kimura and S. Maruyama, Molecular dynamics simulation of heterogeneous nucleation of a liquid droplet on a solid surface, Micr. Thermo. Eng. 6, 3 (2002).
  * Radhakrishnan  _et al._ [2002] R. Radhakrishnan, K. E. Gubbins, and M. Sliwinska-Bartkowiak, Global phase diagrams for freezing in porous media, J. Chem. Phys. 116, 1147 (2002).
  * Takaiwa  _et al._ [2008] D. Takaiwa, I. Hatano, K. Koga, and H. Tanaka, Phase diagram of water in carbon nanotubes, Proc. Natl. Acad. Sci. 105, 39 (2008).
  * Zhou  _et al._ [2019] Y. Zhou, M. Scheffler, and L. M. Ghiringhelli, Determining surface phase diagrams including anharmonic effects, Phys. Rev. B 100, 174106 (2019).
  * Zhang  _et al._ [2019] Z. Zhang, J. A. Schott, M. Liu, H. Chen, X. Lu, B. G. Sumpter, J. Fu, and S. Dai, Prediction of carbon dioxide adsorption via deep learning, Angew. Chem. Int. Ed. 131, 265 (2019).
  * Gurnani  _et al._ [2021] R. Gurnani, Z. Yu, C. Kim, D. S. Sholl, and R. Ramprasad, Interpretable machine learning-based predictions of methane uptake isotherms in metal–organic frameworks, Chem. Mater. 33, 3543 (2021).
  * Fung  _et al._ [2021] V. Fung, G. Hu, P. Ganesh, and B. G. Sumpter, Machine learned features from density of states for accurate adsorption energy prediction, Nat. Comm. 12, 88 (2021).
  * Cui  _et al._ [2023] J. Cui, F. Wu, W. Zhang, L. Yang, J. Hu, Y. Fang, P. Ye, Q. Zhang, X. Suo, Y. Mo, _et al._ , Direct prediction of gas adsorption via spatial atom interaction learning, Nat. Comm. 14, 7043 (2023).
  * Castelvecchi [2016] D. Castelvecchi, Can we open the black box of AI?, Nat. News 538, 20 (2016).
  * Rudin [2019] C. Rudin, Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead, Nat. Mach. Intell. 1, 206 (2019).
  * Lin  _et al._ [2017] H. W. Lin, M. Tegmark, and D. Rolnick, Why does deep and cheap learning work so well?, J. Stat. Phys. 168, 1223 (2017).
  * Košmrlj and Nelson [2017] A. Košmrlj and D. R. Nelson, Statistical mechanics of thin spherical shells, Phys. Rev. X 7, 011002 (2017).
  * Goodrich  _et al._ [2021] C. P. Goodrich, E. M. King, S. S. Schoenholz, E. D. Cubuk, and M. P. Brenner, Designing self-assembling kinetics with differentiable statistical physics models, Proc. Natl. Acad. Sci. 118, e2024083118 (2021).
  * Molina  _et al._ [2021] A. Molina, S. Kumar, S. Karpitschka, and M. Prakash, Droplet tilings for rapid exploration of spatially constrained many-body systems, Proc. Natl. Acad. Sci. 118, e2020014118 (2021).
  * Israelachvili [2011] J. N. Israelachvili,  _Intermolecular and surface forces_ (Academic press, 2011).
  * Hill [2001] T. L. Hill, A different approach to nanothermodynamics, Nano Lett. 1, 273 (2001).
  * Hill [1994] T. L. Hill,  _Thermodynamics of small systems_ (Courier Corporation, 1994).
  * Travalloni  _et al._ [2010] L. Travalloni, M. Castier, F. W. Tavares, and S. I. Sandler, Critical behavior of pure confined fluids from an extension of the van der waals equation of state, J. Supercrit. Fluids 55, 455 (2010).
  * Simon  _et al._ [2014] C. M. Simon, J. Kim, L.-C. Lin, R. L. Martin, M. Haranczyk, and B. Smit, Optimizing nanoporous materials for gas storage, Phys. Chem. Chem. Phys. 16, 5499 (2014).
  * Simon  _et al._ [2017] C. M. Simon, E. Braun, C. Carraro, and B. Smit, Statistical mechanical model of gas adsorption in porous crystals with dynamic moieties, Proc. Natl. Acad. Sci. 114, E287 (2017).
  * Singh  _et al._ [2022] N. Singh, F. Simeski, and M. Ihme, Computing thermodynamic properties of fluids augmented by nanoconfinement: Application to pressurized methane, J. Phys. Chem. B 126, 8623 (2022).
  * Dong  _et al._ [2023] W. Dong, T. Franosch, and R. Schilling, Thermodynamics, statistical mechanics and the vanishing pore width limit of confined fluids, Comm. Phys. 6, 161 (2023).
  * Mayer [1937] J. E. Mayer, The statistical mechanics of condensing systems. i, J. Chem. Phys. 5, 67 (1937).
  * Mayer and Montroll [1941] J. E. Mayer and E. Montroll, Molecular distribution, J. Chem. Phys. 9, 2 (1941).
  * Widom [1963] B. Widom, Some topics in the theory of fluids, J. Chem. Phys. 39, 2808 (1963).
  * Hill and Greenschlag [1961] T. L. Hill and S. Greenschlag, Statistical mechanics of monatomic systems in an external periodic potential field. I. Introduction, virial expansion, and classical second virial coefficient, J. Chem. Phys. 34, 1538 (1961).
  * Stanley [1971] H. E. Stanley,  _Phase transitions and critical phenomena_ , Vol. 7 (Clarendon Press, Oxford, 1971).
  * Nicolis and Prigogine [1989] G. Nicolis and I. Prigogine,  _Exploring complexity: An introduction_ (1989).
  * Kittel and Kroemer [1980] C. Kittel and H. Kroemer, Thermal physics freeman, San Francisco (1980).
  * You  _et al._ [2022] X. You, Y. Li, H. Mo, and Y. Gui, Theoretical studies on lennard-jones parameters of benzene and polycyclic aromatic hydrocarbons, Faraday Discuss. 238, 103 (2022).
  * Dubbeldam  _et al._ [2016] D. Dubbeldam, S. Calero, D. E. Ellis, and R. Q. Snurr, Raspa: molecular simulation software for adsorption and diffusion in flexible nanoporous materials, Mol. Simul. 42, 81 (2016).
  * [67] Supplementary information, see Supplemental Material at [URL will be inserted by publisher] for further information.
  * Auti  _et al._ [2023] G. Auti, Y. Kametani, H. Kimura, S. Paul, W.-L. Hsu, S. Kusaka, R. Matsuda, T. Uemura, S. Chiashi, and H. Daiguji, Effect of pore size on heat release from CO2 adsorption in MIL-101, MOF-177, and UiO-66, J. Mat. Chem. A 11, 20043 (2023).
  * Ball [1993] P. Ball, New horizons in inner space, Nature 361, 297 (1993).
  * Koga  _et al._ [2001] K. Koga, G. Gao, H. Tanaka, and X. C. Zeng, Formation of ordered ice nanotubes inside carbon nanotubes, Nature 412, 802 (2001).
  * Han  _et al._ [2010] S. Han, M. Choi, P. Kumar, and H. E. Stanley, Phase transitions in confined water nanofilms, Nat. Phys. 6, 685 (2010).
  * Rowsell  _et al._ [2004] J. L. Rowsell, A. R. Millward, K. S. Park, and O. M. Yaghi, Hydrogen sorption in functionalized metal- organic frameworks, J. Am. Chem. Soc. 126, 5666 (2004).
  * Yaghi [2012] O. M. Yaghi,  _Hydrogen storage in metal-organic frameworks_ , Tech. Rep. (University of California-Los Angeles, Los Angeles, CA, 2012).
  * Bocquet  _et al._ [1998] L. Bocquet, E. Charlaix, S. Ciliberto, and J. Crassous, Moisture-induced ageing in granular media and the kinetics of capillary condensation, Nature 396, 735 (1998).
  * Restagno  _et al._ [2000] F. Restagno, L. Bocquet, and T. Biben, Metastability and nucleation in capillary condensation, Phys. Rev. Lett. 84, 2433 (2000).
  * Maeda and Israelachvili [2002] N. Maeda and J. N. Israelachvili, Nanoscale mechanisms of evaporation, condensation and nucleation in confined geometries, J. Phys. Chem. B 106, 3534 (2002).
  * Debenedetti [1996] P. G. Debenedetti,  _Metastable liquids: concepts and principles_ (Princeton university press, 1996).
  * Lilley and Prasher [2022] D. Lilley and R. Prasher, Ionocaloric refrigeration cycle, Science 378, 1344 (2022).
  * Shamim  _et al._ [2022] J. A. Shamim, G. Auti, H. Kimura, S. Fei, W.-L. Hsu, H. Daiguji, and A. Majumdar, Concept of a hybrid compression-adsorption heat pump cycle, Cell Rep. Phys. Sci. , 101131 (2022).
  * Page and Monson [1996] K. Page and P. Monson, Monte carlo calculations of phase diagrams for a fluid confined in a disordered porous material, Phys. Rev. E 54, 6557 (1996).
  * Radhakrishnan  _et al._ [2000] R. Radhakrishnan, K. E. Gubbins, and M. Sliwinska-Bartkowiak, Effect of the fluid-wall interaction on freezing of confined fluids: Toward the development of a global phase diagram, J. Chem. Phys. 112, 11048 (2000).
  * De Grandis  _et al._ [2007] V. De Grandis, P. Gallo, and M. Rovere, The phase diagram of confined fluids, J. Mol. Liq. 134, 90 (2007).
  * Lum and Chandler [1998] K. Lum and D. Chandler, Phase diagram and free energies of vapor films and tubes for a confined fluid, Int. J. Thermophys. 19, 845 (1998).
  * Gallo  _et al._ [2016] P. Gallo, K. Amann-Winkel, C. A. Angell, M. A. Anisimov, F. Caupin, C. Chakravarty, E. Lascaris, T. Loerting, A. Z. Panagiotopoulos, J. Russo, _et al._ , Water: A tale of two liquids, Chem. Rev. 116, 7463 (2016).
  * Perin and Gallo [2023] L. Perin and P. Gallo, Phase diagram of aqueous solutions of LiCl: A study of concentration effects on the anomalies of water, J. Phys. Chem. B 127, 4613 (2023).
  * Derjaguin and Kusakov [1936] B. Derjaguin and M. Kusakov, The properties of thin layers of liquids, Proc. Acad. Sci. USSR Chem. Series 5, 741 (1936).
  * Dong [2023] W. Dong, Nanoscale thermodynamics needs the concept of a disjoining chemical potential, Nat. Comm. 14, 1824 (2023).
  * Shi  _et al._ [2021] X. Shi, S. Xue, T. Marhaba, and W. Zhang, Probing internal pressures and long-term stability of nanobubbles in water, Langmuir 37, 2514 (2021).
  * Tan  _et al._ [2021] B. H. Tan, H. An, and C.-D. Ohl, Stability of surface and bulk nanobubbles, Curr. Opin. Colloid Interface Sci 53, 101428 (2021).
  * Maheshwari  _et al._ [2016] S. Maheshwari, M. van der Hoef, X. Zhang, and D. Lohse, Stability of surface nanobubbles: A molecular dynamics study, Langmuir 32, 11116 (2016).
  * Svetovoy  _et al._ [2016] V. B. Svetovoy, I. Devic, J. H. Snoeijer, and D. Lohse, Effect of disjoining pressure on surface nanobubbles, Langmuir 32, 11188 (2016).
  * Witharana  _et al._ [2012] S. Witharana, B. Phillips, S. Strobel, H. Kim, T. McKrell, J.-B. Chang, J. Buongiorno, K. Berggren, L. Chen, and Y. Ding, Bubble nucleation on nano-to micro-size cavities and posts: An experimental validation of classical theory, J. App. Phys. 112 (2012).
  * Chen  _et al._ [2019] J. Chen, K. Zhou, Y. Wang, J. Gao, T. Yuan, J. Pang, S. Tang, H.-Y. Chen, and W. Wang, Measuring the activation energy barrier for the nucleation of single nanosized vapor bubbles, Proc. Natl. Acad. Sci. 116, 12678 (2019).
  * Paul  _et al._ [2020] S. Paul, W.-L. Hsu, M. Magnini, L. Mason, Y. Ho, O. K. Matar, and H. Daiguji, Single-bubble dynamics in nanopores: Transition between homogeneous and heterogeneous nucleation, Phys. Rev. Res. 2, 043400 (2020).
  * Lutsko [2018] J. F. Lutsko, Systematically extending classical nucleation theory, New J. Phys. 20, 103015 (2018).
  * Abyzov and Schmelzer [2014] A. S. Abyzov and J. W. Schmelzer, Heterogeneous nucleation in solutions: Generalized gibbs’ approach, J. Chem. Phys. 140 (2014).
  * Magaletti  _et al._ [2020] F. Magaletti, A. Georgoulas, and M. Marengo, Unraveling low nucleation temperatures in pool boiling through fluctuating hydrodynamics simulations, Int. J. Multiph. Flow 130, 103356 (2020).
  * Shen and Debenedetti [2001] V. K. Shen and P. G. Debenedetti, Density-functional study of homogeneous bubble nucleation in the stretched Lennard-Jones fluid, J. Chem. Phys. 114, 4149 (2001).
  * Lutsko [2011] J. F. Lutsko, Density functional theory of inhomogeneous liquids. iv. squared-gradient approximation and classical nucleation theory, J. Chem. Phys. 134 (2011).
  * Majumdar and Mezic [1999] A. Majumdar and I. Mezic, Instability of ultra-thin water films and the mechanism of droplet formation on hydrophilic surfaces, J. Heat Transfer 121, 964 (1999).
  * Song  _et al._ [2021] L. Song, M. Heiranian, and M. Elimelech, True driving force and characteristics of water transport in osmotic membranes, Desalination 520, 115360 (2021).