Statistical modeling of equilibrium phase transition in confined fluids

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.

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

$$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

$$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:

$$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

$$\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

$$\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

$$\mathcal{H}(\mathbf{p},\mathbf{q}) =\mathcal{K}(\mathbf{p},\mathbf{q})+\mathcal{U}(\mathbf{p},\mathbf{q})$$ (6a)

$$=\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})$:

$$\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.

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

where the configurational contribution is

$$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

$$\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

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

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

$$Z_{q}=\frac{1}{V^{N}}\int_{V} \exp\{-[U_{ma}(q_{1})+\cdots+U_{ma}(q_{N})$$

$$+ 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:

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

Equation (12) can then be written as

$$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

$$\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,

$$Z_{q}= \,\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

$$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

$$Z_{q}= \,\frac{1}{V^{N}}\int_{V}e^{-Nu_{aa}(\mathbf{q})/\tau}dq_{1}\cdots dq_{N}$$

$$+\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

$$\phi \equiv\frac{1}{V}\int_{V}f_{i}dq_{i}$$

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

Using this definition,

$$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

$$\mathbb{Z} =Z_{k-aa}Z_{q-aa}\left(1+N\phi\right)$$

$$=\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

$$\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:

$$\mu^{\text{bulk}}(p_{\text{ext}},T) =\mu^{\text{total}}$$ (23a)

$$=\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]:

$$\Xi_{\text{ads}} =\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)

$$=\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

$$\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

$$\Omega_{\text{ads}}(x,y,z) =-\tau\ln[\Xi_{\text{ads}}(x,y,z)]$$

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

$$=\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:

$$\bar{\Omega}_{\text{ads}} =\frac{1}{V}\int_{V}\Omega_{\text{ads}}(x,y,z)dxdydz$$

$$=\hat{p}V\equiv\xi.$$ (28)

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

For any macroscopic system, we can write

$$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:

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

Therefore,

$$\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:

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

$$=\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

$$\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. (II.3.3) at $N_{\text{max}}$ such that

$$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:

$$\text{Entropy:}~{}(\bar{S}_{\text{ads}}) =\Bigl{(}\partialderivative{\bar{\Omega}_{\text{ads}}}{T}\Bigr{)}_{p,\mu^{\text{ads}}},$$ (35)

$$\text{Gibbs free energy:}~{}(\bar{G}_{\text{ads}}) =\mu^{\text{ads}}\overline{\langle N_{\text{ads}}\rangle},$$ (36)

$$\text{Helmholtz free energy:}~{}(\bar{F}_{\text{ads}}) =\bar{\Omega}_{\text{ads}}+\mu^{\text{ads}}\overline{\langle N_{\text{ads}}\rangle},$$ (37)

$$\text{Enthalpy:}~{}(\bar{H}_{\text{ads}}) =\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

$$\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)]:

$$F_{\text{ads}}(N,V,T) =-\tau\ln[\mathbb{Z}_{\text{bulk}}(1+N\phi)\bigr{]}$$

$$=-N\tau\left\{\ln\left[\frac{(V-Nb)}{\lambda^{3}_{T}N}\right]+1\right\}-\frac{N^{2}a}{V}$$

$$-\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:

$$\mu^{\text{ads}}=\left(\partialderivative{F_{\text{ads}}}{N}\right)_{V,T},$$

$$\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:

$$\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,

$$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 H${}_{2}$ 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 H${}_{2}$ 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.