# Legendre-Fenchel transforms capture layering transitions in porous media

Olav Galteland  PoreLab, Department of Chemistry, Norwegian University of
Science and Technology  Eivind Bering  PoreLab, Department of Physics,
Norwegian University of Science and Technology  Kim Kristiansen  PoreLab,
Department of Chemistry, Norwegian University of Science and Technology  Dick
Bedeaux  PoreLab, Department of Chemistry, Norwegian University of Science and
Technology  Signe Kjelstrup  PoreLab, Department of Chemistry, Norwegian
University of Science and Technology

###### Abstract

We have investigated the state of a nanoconfined fluid in a slit pore in the
canonical and isobaric ensembles. The systems were simulated with molecular
dynamics simulations. The fluid has a transition to a close-packed structure
when the height of the slit approaches the particle diameter. The Helmholtz
energy is a non-convex function of the slit height if the number of particles
does not exceed that of one monolayer. As a consequence, the Legendre
transform cannot be applied to obtain the Gibbs energy. The Gibbs energy of a
non-deformable slit pore can be transformed into the Helmholtz energy of a
deformable slit pore using the Legendre-Fenchel transform. The Legendre-
Fenchel transform corresponds to the Maxwell construction of equal areas.

##  1 Introduction

Over the last years, there has been an increasing number of observations of
phase transitions in confined fluids. Fluids can for instance change their
critical temperature by several tens of degrees [1, 2], a two-dimensional
layer at an interface may develop more than one structure [3], and adsorption
to droplets may depend largely on droplet size [4, 5]. Classical Gibbs
thermodynamics ceases to exist on the nanoscale. The need for the inclusion of
shape and size has been met in several ways. Gibbs and followers included
curvature as a variable [6] to deal with droplet size dependence. Not only
size and shape will matter for the outcome of the analysis of simulations; the
thermodynamic properties will also depend on the small system’s environment or
the set of variables that control the system (the ensemble) according to Hill
[7, 8]. Dong [9] argued that thermodynamic variables, like the surface
tension, change as we shrink the small system, and proposed to add as variable
the integral surface tension, to complement the normal (differential) surface
tension. The pressure of fluids in porous media is of special interest as its
gradient is the main driving force for mass transport [10, 11].

A systematic way to address these problems was given by Hill [7, 8] already 50
years ago. We have argued that the problems are best addressed by his method
[12, 13, 14, 3], because the method provides a general description of small
systems. We have for instance been able to write scaling laws for small system
variables [15], and a new equilibrium criterion for pressure was developed for
two-phase equilibria in slit pores [15]. Hill’s method is therefore our first
choice when the aim is to learn more about structural transitions in confined
fluids or how variables change. We shall find here that a transforming
procedure exists in terms of the Legendre-Fenchel transform. This is a more
general transform than the Legendre transform and can be used for large as
well as small systems. It will enable us to compute the Gibbs energy from the
Helmholtz energy and vice versa.

We have recently reported the changes in free energy during polymer
stretching. The free energy depends on the conditions used, whether the
polymer is stretched at controlled length or force [16, 17]. For sufficiently
short polymers, the Helmholtz energy is a non-convex function of the
controlled length of the polymer. A convex function has a non-negative second
derivative everywhere. To transform from the Helmholtz to the Gibbs energy the
Legendre-Fenchel transform could be applied. Similarly, we observed that the
grand potential of a fluid in a slit pore is a non-convex function of the
distance between the parallel plates at constant chemical potential and
temperature [3]. Also, the Helmholtz energy of solid colloids in solution is
non-convex as a function of the controlled distance between the colloids,
keeping the temperature, volume, and the number of particles constant [18].
Both cases can be explained by a disjoining pressure (also known as the
solvation pressure) [19], which is the excess normal pressure relative to the
bulk pressure due to the packing of fluid particles between the solids. The
observations mentioned all stem from size effects. To obtain the free energy
of the corresponding constant pressure ensemble the Legendre-Fenchel transform
must be applied, and not the Legendre transform.

The three pillars that science progresses from are theory, experiments, and
simulations. The theory part is lacking in nanotechnology. Energy converting
devices are abundant, but there is little available general theory of energy
conversion for the nanoscale. The laws of energy conversion are the laws of
thermodynamics, and the question we are asking is which form these take.
Nanothermodynamics has been constructed, mostly by adding terms to Gibbs’s
classical formulation for large-scale systems. While this mending procedure
may serve the purpose in some special cases, it does not present us with a
systematic procedure to be used as a general tool. Here we argue that the
theory of Hill presented more than 50 years ago, presents an underused
opportunity for a systematic procedure. To show the advantage of this
approach, we study a transition between two structural regimes in a molecular
fluid in a porous media model and document the applicability of an important
tool, namely the Legendre-Fenchel transform.

Several problems arise for systems dependent on size and shape. To which
extent can we still use the thermodynamic tools on nanoscale systems that
apply to macroscale systems? When additional independent variables are needed
in the Gibbs equation, from which pool do we draw them and how?

The free energies of a molecular fluid in a slit pore will be investigated at
two conditions; at constant volume and constant normal pressure. In the first
case, the system is in the canonical ensemble, and in the second case, the
system is in the isobaric ensemble. In the thermodynamic limit, where the free
energy is convex, we can transform the free energy from one ensemble to
another using the Legendre transform. For small systems, this is not always
the case. However, in the case of polymer stretching, we have recently found
that the Legendre-Fenchel transform can be used [17]. This experience has led
us to wonder whether Legendre-Fenchel transforms can be used also for fluids
confined to the slit pore, thereby motivating this paper. In the isochoric
ensemble, the volume of the slit pore is controlled, and we will consider this
as a simple model of a non-deformable porous medium. In the isobaric ensemble,
the pore normal pressure is controlled, and the pore volume can fluctuate. We
will consider this as a simple model of a deformable porous medium.

The paper is outlined as follows. We give the theoretical background in
section 2, including the Legendre-Fenchel transform in section 2.1. We proceed
to present the simulation technique in section 3 and show in section 4, that
the Helmholtz energy of an isochoric slit pore can be transformed into the
Gibbs energy of an isobaric slit pore using the Legendre-Fenchel transform.
The findings are discussed and perspectives are pointed out. In short, we
shall see that systems that are small in Hill’s sense have additional
transitions than the bulk systems have. The system is more restricted when we
control the fluid height than when we control the normal pressure.

###  1.1 System description.

We investigate a single-phase fluid in a slit pore. This can be seen as a
simple model suited to bring out the features described above. The system
consists of a fluid placed between two parallel solid walls, see Fig. 1. The
system has periodic boundary conditions in the $y$\- and $z$-directions. In
the canonical ensemble, the walls do not move, while in the isobaric ensemble,
the top wall can move in the $x$-direction and will act as a piston with
controlled normal pressure on the fluid.

Figure 1: The fluid is placed between two parallel solid walls. The walls are
separated by a distance $h$ (height) and the side lengths of the walls are
$L$. In the canonical ensemble the height $h$ is controlled and the normal
pressure $P$ fluctuates, while in the isobaric ensemble the normal pressure
$P$ is controlled and the height $h$ fluctuates.

In the canonical ensemble the volume of the system is controlled, and the
normal pressure $P$ fluctuates. In the isobaric ensemble, the normal pressure
$P$ is controlled, and the volume $V=L^{2}h=Ah$ fluctuates, where $A=L^{2}$ is
the fluid-solid surface area of one of the walls. The side lengths in the
$y$\- and $z$-directions are fixed equal to $L$, and it is only the distance
between the walls $h$ (height) that fluctuates. The side lengths $L$ are much
larger than the height such that the system may be considered to be
independent of the surface area $A$.

##  2 Theory

In the canonical ensemble, the Helmholtz energy describes the maximum
obtainable work of the system. The total differential of the Helmholtz energy
is

$\text{d}F(N,h,T)=-S\text{d}T-PA\text{d}h+\mu\text{d}N,$ (1)

where $N$ is the number of fluid particles, $T$ is the temperature, $S$ is the
entropy, $P$ is the normal pressure, $A$ is the fluid-solid surface area, and
$\mu$ is the chemical potential. In the isobaric ensemble, it is the Gibbs
energy that describes the maximum obtainable work of the system. The total
differential of the Gibbs energy is

$\text{d}G(N,P,T)=-S\text{d}T+Ah\text{d}P+\mu\text{d}N.$ (2)

The normal pressure and height are defined in terms of the free energies as

$P\equiv-\frac{1}{A}\left(\frac{\partial F}{\partial
h}\right)_{T,N}\quad\text{and}\quad h\equiv\frac{1}{A}\left(\frac{\partial
G}{\partial P}\right)_{T,N}.$ (3)

The difference in the specific Helmholtz energy $\Delta f=F/M$, where $M$ is
the total mass of the fluid, is calculated by integrating the mean normal
pressure as a function of the volume in the canonical ensemble

$\Delta f=f(N,h,T)-f(N,h_{0},T)=-\frac{A}{M}\int_{h_{0}}^{h}\langle
P\rangle\text{d}h^{\prime},$ (4)

where $\langle P\rangle$ is the mean normal pressure. Similarly, the
difference in the specific Gibbs energy $\Delta g=G/M$ is calculated by
integrating the mean height as a function of the normal pressure in the
isobaric ensemble

$\Delta g=g(N,P,T)-g(N,P_{0},T)=\frac{A}{M}\int_{P}^{P_{0}}\langle
h\rangle\text{d}P^{\prime},$ (5)

where $\langle h\rangle$ is the mean height, and $f_{0}=f(N,h_{0},T)$ and
$g_{0}=g(N,P_{0},T)$ are reference states at a large height $h_{0}$ and low
pressure $P_{0}$. The difference in the specific entropy of the system in the
isochoric ensemble is

$\Delta s=\frac{1}{T}(\Delta u-\Delta f),$ (6)

where $\Delta u$ is the change in the specific internal energy of the system.
The changes in specific entropy and internal energy are relative to a system
at height $h_{0}$.

###  2.1 The Legendre-Fenchel transform.

In the isobaric ensemble, the control variables are temperature $T$, normal
pressure $P$, and the number of fluid particles $N$. The normal pressure is
equal to the absolute force acting on the walls divided by the surface area
$A$. This is equal to the normal component of the mechanical pressure tensor
of the fluid [3]. The isobaric partition function can be obtained by a Laplace
transform of the canonical partition function

$\mathcal{Z}(N,P,T)=\beta PA\int_{0}^{\infty}Z(N,h,T)\exp(-\beta
PAh)\text{d}h,$ (7)

where $Z(N,h,T)$ is the canonical partition function and
$\beta=(k_{\text{B}}T)^{-1}$ where $k_{\text{B}}$ is the Boltzmann constant.
The Helmholtz and Gibbs energies are given by the corresponding partition
functions

$F(N,h,T)=-k_{\text{B}}T\ln Z(N,h,T)$ (8)

and

$G(N,P,T)=-k_{\text{B}}T\ln\mathcal{Z}(N,P,T),$ (9)

respectively. From the three equations above it follows that the Gibbs energy
can be obtained from the Helmholtz energy,

$\begin{split}&\exp[-\beta G(N,P,T)]=\\\ &\beta
PA\int_{0}^{\infty}\exp[-\beta(F(N,h,T)+PAh)]\text{d}h.\end{split}$ (10)

For sufficiently high surface number densities $\Gamma=N/A$, the system is
large and the Helmholtz energy $F$ is a differentiable and convex function of
the height, and the above expression reduces to the Legendre transform of the
Helmholtz energy to the Gibbs energy,

$G_{\text{L}}(N,P,T)=G(N,h,T)+PhA.$ (11)

For low surface number densities $\Gamma$, the system is small, and the
Helmholtz energy is non-convex and the integral in equation 10 does not reduce
to the Legendre transform. If the distribution of the normal pressure is
sharply peaked it can however be calculated with a saddlepoint approximation
[17, 20, 21],

$G_{\text{LF}}(N,P,T)=\min_{h}(F(N,h,T)+PhA).$ (12)

This is the Legendre-Fenchel (LF) transform of the Helmholtz energy $F$ to the
Gibbs energy $G_{\text{LF}}$. The LF transform returns only convex functions.
If we apply it again,

$F^{**}=\max_{p}(G_{\text{LF}}+PhA)=G_{\text{LF}}+PhA,$ (13)

we obtain the convex envelope of the Helmholtz energy $F^{**}$. The convex
envelope is the largest function satisfying $F^{**}\leq F$, which is only
equal to the original Helmholtz energy $F$ if it is a convex function. In
other words, the LF transform is not self-inverse if the function is non-
convex [20, 21]. The LF transform can be defined as either the maximum or
minimum. Since $G_{\text{LF}}$ must be convex, we can also obtain $F^{**}$
from a Legendre transform of $G_{\text{LF}}$.

The Maxwell construction of equal areas for liquid-vapor coexistence is
equivalent to the convex envelope of the Helmholtz energy $F^{**}$. The equal
area rule states that for a liquid-vapor coexistence the system follows a
constant pressure $P_{\text{eq}}$ from volume $V_{\text{l}}$ to $V_{\text{g}}$
when the system evaporates, and conversely for condensation. The two volumes
$V_{\text{l}}$ and $V_{\text{g}}$ at the pressure $P_{\text{eq}}$ are the
binodal points of the pressure-volume curve. The equal area rule states

$\int_{V_{\text{l}}}^{V_{\text{g}}}P(V)\text{d}V=P_{\text{eq}}(V_{\text{l}}-V_{\text{g}})$
(14)

where $P(V)$ is a cubic equation of state, for example, the van der Waals
equation, below the critical point. This corresponds to finding the double
tangent line of the non-convex Helmholtz energy. The double tangent line is a
tangent of the function at two different points. A double tangent line does
not exist for convex functions. The double tangent line is exactly the convex
envelope, as it is the largest function that satisfies $F^{**}\leq F$.

##  3 Simulation details

A fluid between two parallel solid walls in the canonical and the isobaric
ensemble was investigated with molecular dynamics simulations using LAMMPS
[22, 23]. The temperature of the fluid was controlled using the Nosé-Hoover
thermostat to $T=2.26T_{c}=2\varepsilon/k_{\text{B}}$, where $T_{c}$ is the
critical temperature of the bulk fluid [24].

The walls were made up of solid particles in a face-centered cubic lattice
with a number density $n_{s}^{*}=1/\sigma^{3}$ corresponding to a lattice
constant $a=2^{2/3}\sigma$, where $\sigma$ is the fluid particle diameter.
Each wall had $N^{p}=5\times 10^{4}$ solid particles and side lengths
$L=100a\approx 159\sigma$. Each solid particle had a mass equal to the fluid
particle mass, implying that the mass of the top solid wall (piston) was
$5\times 10^{4}$ times greater than a fluid particle. The mass of each
particle were equal to $m$, which in reduced units is equal to one. The total
mass of the fluid particles were $M=Nm$. In the canonical ensemble, the solid
particles were fixed in space and could not move. In the isobaric ensemble,
the solid particles in the lower wall were fixed in space, while the solid
particles in the piston were free to move as a single rigid body in the
$x$-direction. The piston could not rotate or move in the $y$\- or
$z$-directions. The fluid particles were placed between the walls. The system
is visualized in Fig. 2 using OVITO [25]. The fluid particles are drawn in red
and the solid particles are drawn in blue.

Figure 2: A rendering of simulations in the isobaric ensemble with surface
number density $\Gamma=0.81$. The top solid wall is the free to move in the
normal direction and acts as a piston on the fluid with a controlled normal
pressure $P$. The height $h$ is shown in (a) and is defined as the minimum
distance between the center of a solid particle in the bottom solid wall and
the center of a solid particle in the piston. The normal pressure and mean
height were equal to (a) $P=0.164\varepsilon/\sigma^{3}$ and $\langle
h\rangle=(10.2\pm 0.09)\sigma$, (b) $P=2.88\varepsilon/\sigma^{3}$ and
$\langle h\rangle=(2.530\pm 0.009)\sigma$ and, (c)
$P=2.89\varepsilon/\sigma^{3}$ and $\langle h\rangle=(1.865\pm 0.003)\sigma$.

The fluid-fluid and fluid-solid particles interacted with the Lennard-
Jones/spline potential [24],

$u^{\text{LJ/s}}(r)=\begin{cases}4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]&\text{if
}r<r_{s},\\\ a(r-r_{c})^{2}+b(r-r_{c})^{3}&\text{if }r_{s}<r<r_{c},\\\
0&\text{else},\end{cases}$ (15)

where $r=|\boldsymbol{r}_{j}-\boldsymbol{r}_{i}|$ is the distance between
particle $i$ and $j$, $\sigma$ is the particle diameter, and $\varepsilon$ is
the well-depth of the interactions. The parameters
$r_{s}=(26/7)^{(1/6)}\sigma\approx 1.24\sigma$,
$a=-24192/3211(\varepsilon/r_{s}^{2})$,
$b=-387072/61009(\varepsilon/r_{s}^{3})$ and $r_{c}=67/48r_{s}\approx
1.74\sigma$ were set such that the potential energy and the force were
continuous at $r_{s}$ and $r_{c}$. An advantage of the Lennard-Jones/spline
potential is that the cut-off is much shorter than the regular Lennard-Jones
potential with a typical cut-off at $2.5\sigma$, which considerably decreases
the computation time. See for more details on the properties of the Lennard-
Jones/spline potential the work of Hafskjold et al [24] and Kristiansen [26].
The fluid-fluid and fluid-solid interactions were equal, while the solid-solid
interaction was zero. The size of the timesteps was set equal to $\delta
t=0.002$. All units are in reduced Lennard-Jones units, meaning that they are
reduced with the mass, particle diameter $\sigma$, the minimum of the
interaction potential $\varepsilon$, and the Boltzmann constant
$k_{\text{B}}$.

The system was initialized by creating one slab of fluid particles between two
slabs of solid particles, all in a face-centered cubic lattice. The fluid
particles were initialized with a velocity such that the temperature was equal
to $T=2\varepsilon/k_{\text{B}}$, and they were free to move for $10^{4}$
timesteps to melt the face-centered cubic lattice of the fluid. In the
canonical ensemble, the piston was moved with a constant velocity for $10^{5}$
timesteps to reach the desired height $h$, from here on defined as the minimum
distance between the center of a solid particle in the bottom solid wall and
the center of a solid particle in the piston. The controlled height was in the
range $h\in[1.7,110]\sigma$. Then the position of the piston was fixed, and
the system was run for $10^{6}$ timesteps to calculate the mean normal
pressure, $\langle P\rangle$.

The mean normal pressure was calculated as the arithmetic mean of the
instantaneous forces of the particles in the piston at time $t$ at time
intervals $\Delta t=0.2$

$\langle P\rangle=\frac{1}{An}\sum_{t=0}^{n}\sum_{i=1}^{N^{p}}f_{x,i,t},$ (16)

where $f_{x,i,t}$ is the force in the $x$-direction acting on the solid
particle $i$ in the piston at time $t$, $N^{p}$ is the number of solid
particles, and $n=10^{6}(\delta t/\Delta t)=10^{4}$ is the number of samples
where $\delta t$ is the size of the timestep. Alternatively, the normal
pressure could be calculated as the normal component of the mechanical
pressure tensor in the fluid [3]. This has been done as a consistency check.

In the isobaric ensemble, an external force was added to the piston after
melting the face-centered cubic lattice of the fluid. The simulations with
controlled normal pressures were run sequentially from low to high normal
pressure to obtain the compression curve, and from high to low normal pressure
to obtain the expansion curve. This was done to reach all available states of
the system in the isobaric ensemble. The controlled normal pressures were in
the range $P\in[0.0045,25]\varepsilon/\sigma^{3}$. Each of the simulations was
run for $10^{6}$ timesteps to calculate the mean height $\langle h\rangle$.

An external force $f_{x}$ in the $x$-direction was applied to each solid
particle in the piston,

$f_{x}=\frac{PA}{N^{p}}$ (17)

where $N^{p}$ is the number of particles in the piston. The mean height
$\langle h\rangle$ was calculated as the arithmetic mean of the instantaneous
height at time $t$ at time intervals of $\Delta t=0.2$

$\langle h\rangle=\frac{1}{n}\sum_{t}^{n}h_{t}.$ (18)

The mean of the specific internal energy was calculated as the mean potential
energy plus the mean kinetic energy of the fluid particles,

$\langle
u\rangle=\frac{1}{M}\left(\sum_{i=1}^{N}\sum_{j>i}^{N}u^{\text{LJ/s}}(r)+\frac{1}{2}\sum_{i=1}^{N}m(\boldsymbol{v}_{i}\cdot\boldsymbol{v}_{i})\right)$
(19)

where $\boldsymbol{v}_{i}$ is the velocity of particle $i$, $m$ is the mass of
each fluid particle and $M=Nm$ is the total mass of the fluid particles. The
specific internal energy was used together with the specific Helmholtz energy
to calculate the specific entropy. The reference states $f_{0}$, $u_{0}$, and
$s_{0}$ were calculated at height $h_{0}=110\sigma_{0}$, and $g_{0}$ was
calculated at pressure $P_{0}$ corresponding to a mean height $\langle
h_{0}\rangle=110\sigma_{0}$.

##  4 Results and discussion

Figure 3: Normal pressure as a function of the height in isochoric conditions
for varying surface number densities. The insert is an enlargement of the
region where structural transitions occur. Figure 4: Specific internal energy
as a function of the height in isochoric conditions for varying surface number
densities. Figure 5: Difference in specific entropy as a function of the
height in isochoric conditions for varying surface number densities. Figure 6:
Difference in specific Helmholtz energy as a function of height for isochoric
conditions for varying surface number densities.

The results are illustrated in Figs. 3 to 10. Fig. 3 shows the normal
pressure-height relationship in isochoric conditions for various surface
number densities. Figs. 4 to 6, give the specific internal energy, difference
in specific entropy, and difference in specific Helmholtz energy,
respectively, all properties as a function of a controlled height. The
difference in specific entropy and Helmholtz energy are given as the
difference relative to the reference state at $h_{0}=110\sigma$. Fig. 7 is a
visualisation of the monolayer in isochoric conditions when close to hexagonal
packing appears.

The impact of the choice of environmental control variables on the normal
pressure-height relationship, the basis of the findings reported, is
illustrated in Fig. 8. This provides a basis for examination of Legendre and
the Legendre-Fenchel transforms, see Figs. 9 and 10. The various results will
now be explained and discussed.

###  4.1 The normal pressure, specific internal energy, entropy, and Helmholtz
energy as a function of height.

The normal pressure-height relationship as a function of height and for
various number densities was presented in Fig. 3. The thermodynamic limit
behavior is seen for large heights or large surface number densities. The
large system has a differentiable convex Helmholtz energy, approximately when
the height $h>3\sigma$ and when the surface number densities are $\Gamma\geq
1.18$. The typical small system behaviour appears for heights $h<3\sigma$ and
surface densities $\Gamma<1.18$.

For smaller heights or surface number densities, the normal pressure goes
through a local minimum and maximum as the height changes. This implies that
the Helmholtz energy is non-convex. The specific internal energy has also a
minimum, see Fig. 4. The entropy is monotonically increasing, except for the
case $\Gamma=1.18$ which has a local minimum, see Fig. 5. The Helmholtz energy
in Fig. 6 captures the trade-off between the internal energy and the entropy.
The Helmholtz energy as a function of height is non-convex for surface
densities $\Gamma<1.18$, which entails that the Legendre transform can not be
applied. The insert in Fig. 6 magnifies the interesting region. The emerging
structures are stabilized by the ability of the system to go to lower energy
and higher entropy. The curves indicate a smooth structural transition.

Figure 7: Top down view of fluid particles in isochoric conditions which has
formed a hexagonal monolayer at height $h=2^{2/3}\sigma\approx 1.59\sigma$.
The solid particles are not shown. The mean normal pressure is $\langle
P\rangle=(9.4\pm 0.1)\varepsilon/\sigma^{3}$, specific internal energy
$\langle u\rangle=(-4.39\pm 0.03)\varepsilon/m$, difference in the specific
entropy $\langle\Delta s\rangle=-9.7k_{\text{B}}/m$, and difference in
specific Helmholtz energy is $\langle\Delta f\rangle=11.4\varepsilon/m$ [25].

###  4.2 A small system phase transition.

Figs. 3 to 6 showed a family of curves that represent the small system in a
transition region for small values of $h$. This transition is special for a
small thermodynamic system, it disappears when surface number density
increases. How can we understand better the behavior of the particles in this
region?

For the system to change from a fluid to a close-packing structure (face-
centered cubic or hexagonal close-packed) without defects, the surface number
density $\Gamma$ must be

$\Gamma=\frac{k}{r^{2}},$ (20)

where $k$ is the number of layers (for a monolayer $k=1$) and $r$ is the
distance between the fluid particles. The total potential energy is at a
minimum at the distance $r=2^{1/6}\sigma\approx 1.12\sigma$, which is where
the Lennard-Jones/spline interaction potential is at a minimum. The
interaction potential is short-ranged, it is zero for distances larger than
the cut-off at $r_{c}\approx 1.74\sigma$. The surface number density for a
monolayer ($k=1$) is $\Gamma\approx 0.79$, and for a double layer ($k=2$) it
is $\Gamma\approx 1.59$. An odd number of layers can complete the face-
centered cubic lattice of walls without faults, while an even number of layers
must have an odd number of stacking faults.

In the extreme case of a monolayer, we assume that the fluid is packed in such
a way that each fluid particle lies on average at a distance
$r_{0}=2^{1/6}\sigma$ away from eight solid particles and four fluid
particles. In addition, there are two solid particle neighbours and four fluid
particle neighbours each at a distance $r_{a}=2^{2/3}\sigma$. All other
particles lie beyond the cutoff distance $r_{c}$ and do not contribute to the
potential energy.

The minimum specific potential energy of a monolayer is then

$e_{p,0}=10u^{\text{LJ/s}}(r_{0})+4u^{\text{LJ/s}}(r_{a})\approx-10.39\varepsilon/m,$
(21)

where $m$ is the particle mass. Taking into account thermal fluctuations, it
is

$\langle e_{p}\rangle=e_{p,0}+a\langle r^{2}\rangle$ (22)

where $a$ is a constant. Assuming no defects, the average deviation from the
average position can be set to zero. The quadratic term represents the lowest
order correction to the mean specific potential energy due to thermal
fluctuations making the fluid particle spend on average more time away from
the average position. By the equipartition theorem, we then obtain the mean
potential energy of a fluid particle

$\langle
e_{p}\rangle=e_{p,0}+\frac{3}{2}k_{\text{B}}T\approx-7.39\varepsilon/m.$ (23)

The expected minimum mean specific internal energy is the sum of the mean
specific potential and kinetic energies

$\langle u\rangle_{\text{min}}=\langle e_{p}\rangle+\langle
e_{k}\rangle\approx-4.39\varepsilon/m$ (24)

where we have used that the mean specific kinetic energy is $\langle
e_{k}\rangle=3k_{\text{B}}T/(2m)=3\varepsilon/m$ at a temperature
$T=2\varepsilon/k_{\text{B}}$.

The observed minimum of the specific internal energy for $\Gamma=0.79$ is
indeed $\langle u\rangle=(-4.39\pm 0.03)\varepsilon/m$ at
$h=2^{2/3}\sigma\approx 1.59\sigma$, see Fig. 4. The simulated structure is
visualised in Fig. 7. The figure shows that the fluid particles have formed a
hexagonal monolayer layer. The mean normal pressure is $\langle
P\rangle=(9.4\pm 0.1)\varepsilon/\sigma^{3}$, difference in the specific
entropy $\langle\Delta s\rangle=-9.7k_{\text{B}}/m$, and difference in
specific Helmholtz energy is $\langle\Delta f\rangle=11.4\varepsilon/m$. The
specific internal energy is at a minimum at this height, however, the
Helmholtz energy is not.

A transition from a fluid to a close-packed structure under stress does not
imply a first-order phase transition, a continuous transition from a fluid to
a close-packed structure packing may occur. In other words, the free energy
can be smooth and continuous during the transition.

Consider for comparison the familiar pressure-volume isotherms of cubic
equations of state, for example, the van der Waals equation, for temperatures
below the critical point. The Helmholtz energy is a non-convex function of
volume. The binodal curve intersects the pressure-volume isotherms at two
points $a=(p_{\text{eq}},V_{\text{l}})$ and $b=(p_{\text{eq}},V_{\text{g}})$
for $T<T_{c}$. The part of the isotherm between point $a$ and $b$ is known as
the van der Waals loop. The Helmholtz energy is a non-convex function in this
region. The line between the points $a$ and $b$ corresponds to the Maxwell
construction of equal areas. See equation 14. The double Legendre-Fenchel
transform of the specific Helmholtz energy $f$ gives its convex envelope
$f^{**}$, which corresponds to the Maxwell construction of equal areas. The
spinodal curve intersects the isotherm at the local minimum and maximum, the
spinodal region is a subset of the binodal region. The binodal region is
metastable, while the spinodal region is unstable. Experimentally it is
observed that fluids do not necessarily follow the van der Waals loop, but
rather the straight line connecting points $a$ and $b$. This is a first-order
phase transition, as the pressure is non-smooth at the points $a$ and $b$. The
system is free to decompose in the binodal and spinodal region, which is
energetically more favorable. During liquid-vapor phase decomposition the
pressure is constant and equal to $P_{\text{eq}}$. The constant pressure
corresponds to the double tangent line in the free energy. This double tangent
line is the largest convex curve that satisfies $f^{**}\leq f$, which is
exactly its convex envelope.

The states in the binodal and spinodal regions can be stable due to the
restrictions that the confinement imposes on the system. The coexistence of
fluid and close-packed structures is not possible in these simulations, as
this would imply that there would be regions with differing heights. This
could be possible if the walls were free to rotate or deform, however, the
system is restricted such that the height is everywhere the same. In isochoric
conditions, there is a smooth transition from a fluid to a close-packed
structure. In isobaric conditions, the system is less restricted, and the
system undergoes a first-order phase transition when it enters the spinodal
region, see Figs. 8 (center). This is because the spinodal region is unstable.
The set of control variables provide different stable states with their
different restrictions on the system.

Figure 8: Normal pressure as a function of the height for isochoric and for
isobaric conditions. The surface number density in each figure is
$\Gamma=0.22$, $\Gamma=0.81$, and $\Gamma=1.18$ from top to bottom. The top
and middle figures shows the isobaric expansion and compression curves.

Figure 9: Specific Gibbs energy as a function of normal pressure. The surface
number densities from top to bottom are $\Gamma=0.22$, $\Gamma=0.81$, and
$\Gamma=1.18$, respectively.

Figure 10: Specific Gibbs energy as a function of height. The surface number
density in the figures from top to bottom are $\Gamma=0.22$, $\Gamma=0.81$,
and $\Gamma=1.18$, respectively.

###  4.3 Response functions.

In the isobaric conditions, it is useful to consider the response function

$K_{N,T}=-A\left(\frac{\partial h}{\partial
P}\right)_{N,T}=-\left(\frac{\partial^{2}G}{\partial P^{2}}\right)_{N,T}\geq
0,$ (25)

where $K_{N,T}/Ah$ is the isothermal compressibility [27]. As the states with
negative $K_{N,T}$ cannot be stable under height fluctuations, $K_{N,T}$ is
restricted to be non-negative. This can be seen in Fig. 8, where the states
with negative $K_{N,T}$ in isochoric conditions are not available in isobaric
conditions. The region with negative compressibility (positive slope)
corresponds to the spinodal region, which is the region between the local
minimum and maximum. In the isochoric ensemble, the height is a control
variable, and states with negative isothermal compressibility can be realized.
In other words, the response function

$\frac{1}{K_{N,T}}=-\frac{1}{A}\left(\frac{\partial P}{\partial
h}\right)_{N,T}=\frac{1}{A^{2}}\left(\frac{\partial^{2}F}{\partial
h^{2}}\right)_{N,T}$ (26)

can be negative. This can be seen in Fig. 8.

The function $g_{\text{LF}}(N,P,T)$ presents two non-smooth points, indicating
a first-order phase transition. The isothermal compressibility in isochoric
conditions is negative, and is well defined in all available states. In
isobaric conditions, the isothermal compressibility is undefined when the free
energy is non-smooth. It is well defined and non-negative for all other
states. The reason for this is that the normal pressure is increased by
applying a directed force that seeks to compress the system. If the system was
to increase its volume in response to this compression force, the fluid center
of mass would have to move in the direction opposite to the applied force,
thus violating conservation of linear momentum.

###  4.4 The specific Gibbs energy as a function of normal pressure. Legendre-
Fenchel transforms.

Figs. 8 presents the normal pressure-height relationships for isochoric and
isobaric conditions. The states available for the system to follow in
isochoric conditions are blue, while isobaric conditions provide states given
by the orange and green curves. The orange curve is for isobaric compression
and the green curve is for isobaric expansion. Once the maximum in a curve is
reached during isobaric compression, the system will switch to a smaller
height. Alternatively, by isobaric expansion all points near the minimum
become accessible. The dotted line is the height computed from the derivative
of the Legendre-Fenchel transform $g_{\text{LF}}$ with respect to the normal
pressure $P$, see equation 3. This curve gives the Maxwell construction of
equal areas, or in other words the straight line across the van der Waals loop
or the binodal region. For isobaric conditions, the system is metastable in
the binodal region and unstable in the spinodal region.

The specific Gibbs energy is presented in Figs. 9 and 10 as a function of
normal pressure and height, respectively. For the large thermodynamic system,
i.e. the bottom panels where $\Gamma=1.18$, the specific Gibbs energy of
compression (orange curve) coincides with the Legendre transform (blue curve),
and the Legendre-Fenchel transform (black dotted curve) of the specific
Helmholtz energy. This is the expected behavior of large systems.

At lower surface number densities, the Legendre transform ceases to apply. But
we can still understand the system in terms of its thermodynamic properties.
The normal pressure-height curves form a van der Waals loop for isochoric
conditions. The specific Gibbs energy from compression differs from that of
expansion for isobaric conditions. The Legendre-Fenchel transform (dotted
curve) follows the Legendre transform (blue curve), except in the van der
Waals loop, which is cut out by the Legendre-Fenchel transform.

The underlying distribution of normal pressures in the isochoric conditions is
highly peaked, implying that the conditions for the saddlepoint approximation,
necessary for the Legendre-Fenchel transform, are obeyed. The small error bar
in the normal pressure, which is the standard deviation, testifies to this,
see Fig. 8. The mean relative standard deviation is $1.1$%, $0.4$%, $0.09$%
for $\Gamma=0.22$, $\Gamma=0.81$, and $\Gamma=1.18$, respectively.

The results of the slit pore simulations show that structural changes are more
restricted when in isochoric conditions than in isobaric conditions. The
findings are similar to observations of polymer stretching [16]. Also here a
transition was found between states. But then the different regimes referred
to the type of degrees of freedom of the molecule (active or frozen rotational
or stretching degrees of freedom) [16, 17]. We have seen that the Legendre-
Fenchel transforms apply to two widely different cases, so we may pose the
question: Does the transform apply in general to energy conversion in small
systems? From the mathematics point of view, this seems likely [28, 29]. More
data is needed before we may conclude, but this study brings out an
interesting perspective. If the answer is yes, we may have a new tool for
energy conversion in small systems.

When the Helmholtz energy is non-convex in the isochoric ensemble the system
exhibits hysteresis in the isobaric ensemble. Hysteresis means that the state
of the system depends on the system’s history. The specific Gibbs energy
depends on whether the system comes from a compressed or expanded state. In
the isobaric ensemble, available states can be explored by first compressing
and then expanding the system. One of the two available states at the same
controlled normal pressure in the compression and expansion curves are
metastable, the other state is stable. The Legendre-Fenchel transform gives
the stable states. Since the piston cannot deform the energy barrier to go
from a metastable state to a stable state is high, which is the reason for the
pronounced hysteresis in this system. The loop created in the specific Gibbs
energy profile in isobaric conditions is interesting. The existence of a loop
means in principle, that work can be extracted from the loop by only two
steps, namely compression, and expansion. The slit pore may serve as a very
first model for deformable porous media in this context.

##  5 Conclusions

The normal pressure varies with the height. The excess normal pressure, the
normal pressure minus the bulk pressure, is often called the disjoining or
solvation pressure [19, 3]. The disjoining pressure is typically defined in
the grand canonical ensemble (an open system). The disjoining pressure
oscillates with a period equal to the fluid particle diameter because of the
fluid packing between the walls. The mechanism of the excess normal pressure
is the same here as for an open system, but the behaviour of the normal
pressure observed in the closed system is different than for an open system.
The observation of disjoining pressure is not new, see for example the work by
Israelachvili [19]. It has been observed in simulations as well as
experimentally. However, Israelachvili does not apply thermodynamics and the
machinery that follows it. In this work, we have expanded upon the work by
Israelachvili and others by giving a thermodynamic description. As a result,
we get additional relations that can be used, for example, Maxwell relations.

In isochoric conditions, there is a smooth transition from a fluid to a close-
packed structure. The Helmholtz energy of this smooth transition is non-convex
for small surface number densities. The system is restricted in such a way
that there must be a smooth transition. It is not possible for the system to
have coexistence of a fluid and a close-packed structure because that would
imply that the system could have two different heights at the same time. A
possible future generalization could be to allow the solid walls to deform in
a less restricted manner such that it would be possible for the system to have
varying heights. In this way the system would allow for coexistence of a fluid
and a close-packed structure. We hypothesize that the hysteresis in such a
system will be less pronounced. Vapor-liquid coexistence is not restricted in
this way, the vapor and liquid can have two different volumes. Since the
height is controlled, the system is also restricted to that height. However,
the normal pressure is controlled in isobaric conditions. As a result, there
is a first-order phase transition from a fluid to a close-packed structure for
small surface number densities. There is no coexistence of fluid and close-
packed structure for the same reason as for isochoric conditions. The Legendre
transform does not apply to non-convex free energies, and we have shown the
Legendre-Fenchel transform must be applied.

Small systems in Hill’s sense are not extensive. They are characterized by
giving different responses to their ensemble of control. In the present case,
we have studied and compared two small systems, i.e. isochoric and isobaric
fluids confined to slit pores. Despite their smallness, we have found that
they are perfectly well describable by thermodynamics when the theory is
adjusted to deal with smallness. One such adjustment means to use Legendre-
Fenchel transforms, rather than Legendre transforms. Doing that, we have shown
that the specific Helmholtz energy can be transformed into the specific Gibbs
energy. The findings are general, and support the systematic approach of Hill
for descriptions of other small systems.

## Author Contributions

O.G. contributed to formal analysis, investigation, methodology, software, and
visualisation. D.B. and S.K contributed to supervision. All authors
contributed to conceptualization, writing original drafts, reviewing, and
editing. All authors have read and agreed to the published version of the
manuscript.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

The simulations were performed on resources provided by UNINETT Sigma2 - the
National Infrastructure for High Performance Computing and Data Storage in
Norway. We thank the Research Council of Norway through its Centres of
Excellence funding scheme, project number 262644, PoreLab.

## References

  * [1] JK Brennan and W Dong.  Phase transitions of one-component fluids adsorbed in random porous media: Monte Carlo simulations.  The Journal of chemical physics, 116(20):8948–8958, 2002.
  * [2] John K Brennan and Wei Dong.  Molecular simulation of the vapor-liquid phase behavior of Lennard-Jones mixtures in porous solids.  Physical review E, 67(3):031503, 2003.
  * [3] Olav Galteland, Dick Bedeaux, and Signe Kjelstrup.  Nanothermodynamic description and molecular simulation of a single-phase fluid in a slit pore.  Nanomaterials, 11(1):165, 2021.
  * [4] Bjørn A Strøm, Jianying He, Dick Bedeaux, and Signe Kjelstrup.  When thermodynamic properties of adsorbed films depend on size: Fundamental theory and case study.  Nanomaterials, 10(9):1691, 2020.
  * [5] Bjørn A Strøm, Dick Bedeaux, and Sondre K Schnell.  Adsorption of an Ideal Gas on a Small Spherical Adsorbent.  Nanomaterials, 11(2):431, 2021.
  * [6] Øivind Wilhelmsen, Dick Bedeaux, and David Reguera.  Tolman length and rigidity constants of the Lennard-Jones fluid.  The Journal of chemical physics, 142(6):064706, 2015.
  * [7] Terrell L Hill.  Thermodynamics of small systems, part 1.  Benjamin, 1963.
  * [8] Terrell L Hill.  Thermodynamics of small systems, part 2.  Benjamin, 1964.
  * [9] W. Dong.  Thermodynamics of interfaces extended to nanoscales by introducing integral and differential surface tensions.  Proceedings of the National Academy of Sciences, 118(3), 2021.
  * [10] Signe Kjelstrup, Dick Bedeaux, Alex Hansen, Bjørn Hafskjold, and Olav Galteland.  Non-isothermal transport of multi-phase fluids in porous media. the entropy production.  Frontiers in Physics, 6:126, 2018.
  * [11] Signe Kjelstrup, Dick Bedeaux, Alex Hansen, Bjørn Hafskjold, and Olav Galteland.  Non-isothermal transport of multi-phase fluids in porous media. Constitutive equations.  Frontiers in Physics, 6:150, 2019.
  * [12] Dick Bedeaux, Signe Kjelstrup, and Sondre K Schnell.  Nanothermodynamics. General theory.  NTNU, Trondheim, Norway, 2020.
  * [13] Olav Galteland, Dick Bedeaux, Bjørn Hafskjold, and Signe Kjelstrup.  Pressures inside a nano-porous medium. The case of a single phase fluid.  Frontiers in Physics, 7:60, 2019.
  * [14] Máté Erdős, Olav Galteland, Dick Bedeaux, Signe Kjelstrup, Othonas A Moultos, and Thijs JH Vlugt.  Gibbs ensemble Monte Carlo simulation of fluids in confinement: Relation between the differential and integral pressures.  Nanomaterials, 10(2):293, 2020.
  * [15] Michael T Rauter, Olav Galteland, Máté Erdős, Othonas A Moultos, Thijs JH Vlugt, Sondre K Schnell, Dick Bedeaux, and Signe Kjelstrup.  Two-phase equilibrium conditions in nanopores.  Nanomaterials, 10(4):608, 2020.
  * [16] Eivind Bering, Signe Kjelstrup, Dick Bedeaux, J Miguel Rubi, and Astrid S de Wijn.  Entropy production beyond the thermodynamic limit from single-molecule stretching simulations.  The Journal of Physical Chemistry B, 124(40):8909–8917, 2020.
  * [17] Eivind Bering, Dick Bedeaux, Signe Kjelstrup, Astrid S de Wijn, Ivan Latella, and J Miguel Rubi.  A Legendre–Fenchel Transform for Molecular Stretching Energies.  Nanomaterials, 10(12):2355, 2020.
  * [18] Olav Galteland, Fernando Bresme, and Bjørn Hafskjold.  Solvent-Mediated Forces between Ellipsoidal Nanoparticles Adsorbed at Liquid–Vapor Interfaces.  Langmuir, 36(48):14530–14538, 2020.
  * [19] Jacob N Israelachvili.  Intermolecular and surface forces.  Academic press, 2015.
  * [20] Hugo Touchette.  Legendre-Fenchel transforms in a nutshell, 2005.  available at: https://appliedmaths.sun.ac.za/h̃touchette/archive/notes/lfth2.pdf (Nov. 2021).
  * [21] Ralph Tyrell Rockafellar.  Convex analysis.  Princeton university press, 2015.
  * [22] Steve Plimpton.  Fast parallel algorithms for short-range molecular dynamics.  Journal of computational physics, 117(1):1–19, 1995.
  * [23] Aidan P Thompson, H Metin Aktulga, Richard Berger, Dan S Bolintineanu, W Michael Brown, Paul S Crozier, Pieter J in’t Veld, Axel Kohlmeyer, Stan G Moore, Trung Dac Nguyen, et al.  LAMMPS-A flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.  Computer Physics Communications, page 108171, 2021.
  * [24] Bjørn Hafskjold, Karl Patrick Travis, Amanda Bailey Hass, Morten Hammer, Ailo Aasen, and Øivind Wilhelmsen.  Thermodynamic properties of the 3D Lennard-Jones/spline model.  Molecular Physics, 117(23-24):3754–3769, 2019.
  * [25] Alexander Stukowski.  Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool.  Modelling and Simulation in Materials Science and Engineering, 18(1):015012, 2009.
  * [26] Kim R. Kristiansen.  Transport Properties of the Simple Lennard-Jones/Spline Fluid I: Binary Scattering and High-Accuracy Low-Density Transport Coefficients.  Frontiers in Physics, 8:271, 2020.
  * [27] Alessandro Campa, Lapo Casetti, Ivan Latella, Agustín Pérez-Madrid, and Stefano Ruffo.  Concavity, response functions and replica energy.  Entropy, 20(12):907, 2018.
  * [28] Jeffrey Commons, Ying-Jen Yang, and Hong Qian.  Duality Symmetry, Two Entropy Functions, and an Eigenvalue Problem in Gibbs’ Theory.  arXiv preprint arXiv:2108.08948, 2021.
  * [29] Hong Qian.  Thermodynamic Behavior of Statistical Event Counting in Time: Independent and Correlated Measurements.  arXiv preprint arXiv:2109.12806, 2021.