Legendre-Fenchel transforms capture layering transitions in porous media

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,³ 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⁶ 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⁹ 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,^3,12–14 because the method provides a general description of small systems. In this work, we investigate a slit pore of large surface area and small height. The system has periodic boundary conditions in the directions tangential to the slit pore wall. In this way, the system is confined, but the number of particles is large. We have for instance been able to write scaling laws for small system variables,¹⁵ and a new equilibrium criterion for pressure was developed for two-phase equilibria in slit pores.¹⁵ 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. A first-order phase transition is defined as a discontinuity in the relation between the normal pressure and the slit pore height. In confinement a phase transition is no longer sharp, the transition is smooth,⁸ this is by definition not a first-order phase transition. 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.³ 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.¹⁸ Both cases can be explained by a disjoining pressure (also known as the solvation pressure),¹⁹ 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.¹⁷ 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.

Fig. 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²h = Ah fluctuates, where A = L² 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

dF(N,h,A,T) = −SdT − PAdh + 2γdA + μdN, (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, γ is the fluid-solid surface tension and μ 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

dG(N,P,A,T) = −SdT + AhdP + 2γdA + μdN. (2)

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

The difference in the specific Helmholtz energy Δ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

where 〈P〉 is the mean normal pressure. Similarly, the difference in the specific Gibbs energy Δg = G/M is calculated by integrating the mean height as a function of the normal pressure in the isobaric ensemblewhere 〈h〉 is the mean height, and f₀ = f(N,h₀,A,T) and g₀ = g(N,P₀,A,T) are reference states at a large height h₀ and low pressure P₀. The difference in the specific entropy of the system in the isochoric ensemble iswhere Δ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₀.

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.³ The isobaric partition function can be obtained by a Laplace transform of the canonical partition functionwhere Z(N,h,A,T) is the canonical partition function and β = (k_BT)⁻¹ where k_B is the Boltzmann constant. The Helmholtz and Gibbs energies are given by the corresponding partition functions

F(N,h,A,T) = −k_BT ln Z(N,h,A,T) (8)

and

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

For sufficiently high surface number densities Γ = 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_L(N,P,A,T) = F(N,h,A,T) + PhA. (11)

For low surface number densities Γ, the system is small, and the Helmholtz energy is non-convex and the integral in eqn (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

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

we obtain the convex envelope of the Helmholtz energy F**. The convex envelope is the largest function satisfying F** ≤ 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_LF must be convex, we can also obtain F** from a Legendre transform of G_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_eq from volume V_l to V_g when the system evaporates, and conversely for condensation. The two volumes V_l and V_g at the pressure P_eq are the binodal points of the pressure-volume curve. The equal area rule states

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** ≤ 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ε/k_B, where T_c = 0.885 is the critical temperature of the bulk fluid.²⁴ The temperature was controlled above the critical point to avoid the liquid-vapor coexistence region.

The walls were made up of solid particles in a face-centered cubic lattice with a number density corresponding to a lattice constant a = 2^2/3σ, where σ is the fluid particle diameter. Each wall had N^p = 5 × 10⁴ solid particles and side lengths L = 100a ≈ 159σ. Each solid particle had a mass equal to the fluid particle mass, implying that the mass of the top solid wall (piston) was 5 × 10⁴ 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.²⁵ The fluid particles are drawn in red and the solid particles are drawn in blue. The simulation box had periodic boundary conditions in the y- and z-directions.

Fig. 2 A rendering of simulations in the isobaric ensemble with surface number density Γ = 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ε/σ³ and 〈h〉 = (10.2 ± 0.09)σ, (b) P = 2.88ε/σ³ and 〈h〉 = (2.530 ± 0.009)σ and, (c) P = 2.89ε/σ³ and 〈h〉 = (1.865 ± 0.003)σ.

The fluid–fluid and fluid–solid particles interacted with the Lennard–Jones/spline potential,²⁴

where r = |r_j − r_i| is the distance between particle i and j, σ is the particle diameter, and ε is the well-depth of the interactions. The parameter r_s = (26/7)^(1/6)σ ≈ 1.24σ is the inflection point of the Lennard-Jones potential, and the parameters a = −24 192/3211(ε/r_s2), b = −387 072/61 009(ε/r³_s) and r_c = 67/48r_s ≈ 1.74σ were set such that the potential energy and the force were continuous and smooth at r_s. The potential is suitable to model noble gases or methane. 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σ, 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.²⁴ and Kristiansen.²⁶ 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 δt = 0.002. All units are in reduced Lennard-Jones units, meaning that they are reduced with the mass, particle diameter σ, the minimum of the interaction potential ε, and the Boltzmann constant k_B. The Lennard-Jones/spline potential is suitable to model argon or methane. However, we do not.

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ε/k_B, and they were free to move for 10⁴ 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⁵ 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 ∈ [1.7, 110]σ. Then the position of the piston was fixed, and the system was run for 10⁶ timesteps to calculate the mean normal pressure, 〈P〉.

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 Δt = 0.2

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⁶(δt/Δt) = 10⁴ is the number of samples where δ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.³ 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 ∈ [0.0045, 25]ε/σ³. Each of the simulations was run for 10⁶ timesteps to calculate the mean height 〈h〉.

An external force f_x in the x-direction was applied to each solid particle in the piston,

where N^p is the number of particles in the piston. The mean height 〈h〉 was calculated as the arithmetic mean of the instantaneous height at time t at time intervals of Δt = 0.2

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

where 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₀, u₀, and s₀ were calculated at height h₀ = 110σ₀, and g₀ was calculated at pressure P₀ corresponding to a mean height 〈h₀〉 = 110σ₀.

4 Results and discussion

The results are illustrated in Fig. 3 to 10. Fig. 3 shows the normal pressure-height relationship in isochoric conditions for various surface number densities. To investigate that the results are independent of the surface area, we have simulated a system with a surface area four times larger (results not shown). We found that the results were independent of the surface area. Fig. 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₀ = 110σ. Fig. 7 visualises two cases in isochoric conditions when close to hexagonal packing appears.

Fig. 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.

Fig. 4 Specific internal energy as a function of the height in isochoric conditions for varying surface number densities.

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 Fig. 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σ and when the surface number densities are Γ ≥ 1.18. The typical small system behaviour appears for heights h < 3σ and surface densities Γ < 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 Γ = 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 Γ < 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.

Fig. 5 Difference in specific entropy as a function of the height in isochoric conditions for varying surface number densities.

Fig. 6 Difference in specific Helmholtz energy as a function of height for isochoric conditions for varying surface number densities.

4.2 A small system phase transition

Fig. 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. The top panel in Fig. 7 illustrates the fluid particles at height h = 2^2/3σ ≈ 1.59σ and surface number density Γ = 0.22. The fluid particles cluster together and form dense patches. How can we understand better the behavior of the particles in this region?

Fig. 7 Top down view of fluid particles in isochoric conditions at height h = 2^2/3σ ≈ 1.59σ and surface number densities (top) Γ = 0.22 and (bottom) Γ = 0.79. In the bottom case the particles form a monolayer without defects. The solid particles are not shown. The mean normal pressure is 〈P〉 = (9.4 ± 0.1)ε/σ³, specific internal energy 〈u〉 = (−4.39 ± 0.03)ε/m, difference in the specific entropy 〈Δs〉 = −9.7k_B/m, and difference in specific Helmholtz energy is 〈Δf〉 = 11.4ε/m.²⁵

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 Γ must be

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σ ≈ 1.12σ, 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 ≈ 1.74σ. The surface number density for a monolayer (k = 1) is Γ ≈ 0.79, and for a double layer (k = 2) it is Γ ≈ 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₀ = 2^1/6σ 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σ. 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^LJ/s(r₀) + 4u^LJ/s(r_a) ≈ −10.39ε/m, (21)

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

〈e_p〉 = e_p,0 + a〈r²〉 (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

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

〈u〉_min = 〈e_p〉 + 〈e_k〉 ≈ −4.39ε/m (24)

where we have used that the mean specific kinetic energy is 〈e_k〉 = 3k_BT/(2m) = 3ε/m at a temperature T = 2ε/k_B.

The observed minimum of the specific internal energy for Γ = 0.79 is indeed 〈u〉 = (−4.39 ± 0.03)ε/m at h = 2^2/3σ ≈ 1.59σ, see Fig. 4. The simulated structure is visualised in the bottom panel of Fig. 7. The figure shows that the fluid particles have formed a hexagonal monolayer layer. The mean normal pressure is 〈P〉 = (9.4 ± 0.1)ε/σ³, difference in the specific entropy 〈Δs〉 = −9.7k_B/m, and difference in specific Helmholtz energy is 〈Δf〉 = 11.4ε/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_eq, V_l) and b = (p_eq, V_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 eqn (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_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** ≤ 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 Fig. 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.

Fig. 8 Normal pressure as a function of the height for isochoric and for isobaric conditions. The surface number densities from top to bottom are Γ = 0.22, Γ = 0.81, and Γ = 1.18, respectively. The top and middle figures shows the isobaric expansion and compression curves.

4.3 Response functions

In the isobaric conditions, it is useful to consider the response functionwhere K_N,T/Ah is the isothermal compressibility.²⁷ 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 functioncan be negative. This can be seen in Fig. 8.

The function g_LF(N,P,A,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

Fig. 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_LF with respect to the normal pressure P, see eqn (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 Fig. 9 and 10 as a function of normal pressure and height, respectively. For the large thermodynamic system, i.e. the bottom panels where Γ = 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.

Fig. 9 Specific Gibbs energy as a function of normal pressure. The surface number densities from top to bottom are Γ = 0.22, Γ = 0.81, and Γ = 1.18, respectively.

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

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 Γ = 0.22, Γ = 0.81, and Γ = 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.¹⁶ 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 and perspectives

The normal pressure varies with the height of the slit pore. The excess normal pressure, the normal pressure minus the bulk pressure, is often called the disjoining or solvation pressure.^3,19 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.¹⁹ In the height-controlled system, there is no first-order phase transition. We observe a smooth transition from a fluid to a close-packed structure with oscillating normal pressure, in a manner similar to disjoining pressures. Disjoining pressures have been observed in simulations as well as by experiments. 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 of this work 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 might 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 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. The system is large in the surface area A, and small in the height h. 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 by 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.

Notes and references

1. J. Brennan and W. Dong, J. Chem. Phys., 2002, 116, 8948–8958.
2. J. K. Brennan and W. Dong, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2003, 67, 031503.
3. O. Galteland, D. Bedeaux and S. Kjelstrup, Nanomaterials, 2021, 11, 165.
4. B. A. Strøm, J. He, D. Bedeaux and S. Kjelstrup, Nanomaterials, 2020, 10, 1691.
5. B. A. Strøm, D. Bedeaux and S. K. Schnell, Nanomaterials, 2021, 11, 431.
6. Ø. Wilhelmsen, D. Bedeaux and D. Reguera, J. Chem. Phys., 2015, 142, 064706.
7. T. L. Hill, Thermodynamics of small systems, part 1, Benjamin, 1963.
8. T. L. Hill, Thermodynamics of small systems, part 2, Benjamin, 1964.
9. W. Dong, Proc. Natl. Acad. Sci. U. S. A., 2021, 118 DOI:10.1073/pnas.2019873118.
10. S. Kjelstrup, D. Bedeaux, A. Hansen, B. Hafskjold and O. Galteland, Front. Phys., 2018, 6, 126.
11. S. Kjelstrup, D. Bedeaux, A. Hansen, B. Hafskjold and O. Galteland, Front. Phys., 2019, 6, 150.
12. D. Bedeaux, S. Kjelstrup and S. K. Schnell, Nanothermodynamics. General theory, NTNU, Trondheim, Norway, 2020.
13. O. Galteland, D. Bedeaux, B. Hafskjold and S. Kjelstrup, Front. Phys., 2019, 7, 60.
14. M. Erdős, O. Galteland, D. Bedeaux, S. Kjelstrup, O. A. Moultos and T. J. Vlugt, Nanomaterials, 2020, 10, 293.
15. M. T. Rauter, O. Galteland, M. Erdős, O. A. Moultos, T. J. Vlugt, S. K. Schnell, D. Bedeaux and S. Kjelstrup, Nanomaterials, 2020, 10, 608.
16. E. Bering, S. Kjelstrup, D. Bedeaux, J. M. Rubi and A. S. de Wijn, J. Phys. Chem. B, 2020, 124, 8909–8917.
17. E. Bering, D. Bedeaux, S. Kjelstrup, A. S. de Wijn, I. Latella and J. M. Rubi, Nanomaterials, 2020, 10, 2355.
18. O. Galteland, F. Bresme and B. Hafskjold, Langmuir, 2020, 36, 14530–14538.
19. J. N. Israelachvili, Intermolecular and surface forces, Academic Press, 2015.
20. H. Touchette, Legendre-Fenchel transforms in a nutshell, 2005, available at: https://appliedmaths.sun.ac.za/htouchette/archive/notes/lfth2.pdf, Nov. 2021.
21. R. T. Rockafellar, Convex analysis, Princeton University Press, 2015.
22. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
23. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in't Veld, A. Kohlmeyer, S. G. Moore and T. D. Nguyen, et al. , Comput. Phys. Commun., 2021, 108171.
24. B. Hafskjold, K. P. Travis, A. B. Hass, M. Hammer, A. Aasen and Ø. Wilhelmsen, Mol. Phys., 2019, 117, 3754–3769.
25. A. Stukowski, Model. Simul. Mater. Sci. Eng., 2009, 18, 015012.
26. K. R. Kristiansen, Front. Phys., 2020, 8, 271.
27. A. Campa, L. Casetti, I. Latella, A. Pérez-Madrid and S. Ruffo, Entropy, 2018, 20, 907.
28. J. Commons, Y.-J. Yang and H. Qian, 2021, arXiv preprint arXiv:2108.08948.
29. H. Qian, 2021, arXiv preprint arXiv:2109.12806.

---