Inside and out: Surface thermodynamics from positive to negative curvature

On first glance, a nanoparticle solvated by a fluid and the same fluid confined in a nanopore would seem to be very different physical systems; however, they both can be described as a fluid at a curved surface—a convex surface in the case of the nanoparticle and concave one for the nanopore. This raises the question as to whether these two situations can be understood within a unified thermodynamic description. The central quantity governing the thermodynamics of surfaces is the surface (or interfacial) free energy, γ, which plays a significant role in a number of phenomena of scientific and technological import, including wetting, dendritic growth and morphology, nucleation, and nanoparticle stability.^1–3 Although most experimental and computational studies of solid–liquid interfaces focus on planar interfaces, the role of curvature (or shape) on interface thermodynamics has long been of interest since the pioneering work of Gibbs;⁴ however, besides recognition of the fact that the local curvature of the interface plays a role in determining the work required to deform the interface, Gibbs restricted further study to planar interfaces.

In the years since Gibbs, significant effort has been focused on the development of quantitative descriptions of the curvature dependence of γ. Perhaps the best known of these comes from the work of Tolman⁵ for a fluid at a spherical surface of radius R, in which

where γ₀ is the surface free energy of a planar interface and δ is the so-called Tolman length. A number of approaches have been since put forward to generalize the ideas of Tolman to more complex surfaces. For a 3d surface, Helfrich postulated that γ can be expanded in a power series in the average mean $(H̄)$ and average Gaussian $(K̄)$ curvatures of the surface,⁶ 

where the sum is usually truncated at low order. For a spherical wall of radius R, where $H̄=1/R$ and $K̄=1/R2$⁠, this would give

For a 2d surface, the Helfrich expression would be a Taylor series expansion for γ in powers of the average curvature of the surface $C̄$⁠, which for a circular wall (disk) would be equal to 1/R.

A more recent development is the so-called Morphometric Thermodynamics (MT),^7–9 based on a theorem from integral geometry due to Hadwiger¹⁰ that states that a valuation defined on a convex set $S$ in d dimensions can be completely described as a linear combination of d + 1 geometric measures (volume, surface area and d − 1 measures of curvature) as long as the valuation is continuous, motion invariant, and additive. In MT, the valuation is the solvation free energy, Ω, of a nanoparticle with a surface (represented by a convex set $S$⁠) embedded in a fluid medium.⁸ 

For a 3d system, the MT solvation free energy of the nanoparticle is then given by

where V, A, H, and K are the surface volume, area, integrated mean curvature, and integrated Gaussian curvature, respectively. Defining the surface free energy γ as γ ≡ (Ω − pV)/A gives

where $H̄S$ and $K̄S$ are the average mean and Gaussian curvatures of $S$⁠, respectively. For a spherical wall particle of radius R$(C̄=1/R)$⁠, the MT expression is

Here, γ₁ is equal to −2δ, where δ is the Tolman coefficient. For the 2d hard-disk fluid at a circular wall system, there is only one curvature, $C̄$⁠, giving

For a 3d fluid of hard spheres of diameter σ at a spherical hard wall, molecular-dynamics (MD) simulations^11,12 have shown that the MT description [Eq. (6)] works extremely well at low to intermediate packing fractions with deviations (beyond quadratic in 1/R) outside of the statistical error appearing only at high packing fractions (η ≡ πσ³ρ/6 > 0.42) approaching the melting density. These simulation results are in agreement with previous calculations using classical Density Functional Theory (cDFT)¹³ on the same system,¹⁴ which showed no deviations from MT below η = 0.42, with no solutions to the cDFT equations above this value. Similar agreement was seen in cDFT calculations for the hard-sphere fluid at a variety of polyhedral walls.³ For the same fluid at a cylindrical hard wall (3d), deviations from MT at low packing fractions were observed in the simulations, but they required very high-precision calculations to resolve.¹² In the case of a 2d fluid of hard disks with diameter σ at a circular hard wall of radius R, previous Monte Carlo simulations by the authors showed that the MT expression [Eq. (7)] with $C̄=1/R$ is also very accurate over a wide range of fluid packing fractions and wall curvatures.¹⁵ (Note that, in what follows, we assume, without loss of generality, that all lengths are measured in units of σ and all energies are in units of k_BT, where k_B is the Boltzmann constant.)

A central question is whether the MT formalism can be extended to fluids confined inside a container with curved boundary—for example, within a spherical (3d) or circular (2d) container. Relative to the solvated particle case, where the fluid is outside the curved boundary, this corresponds to a wall curvature that is negative. (Note that, Hadwiger’s theorem does not apply for such concave systems, as it is limited to compact convex sets.) For example, a hard-disk fluid (2d) contained within a circle of radius R, the average curvature $C̄=−1/R$⁠. In the limit that R → ∞ (zero curvature), both the solvated particle and the confined system have the same limit, namely the fluid in contact with a flat wall and we expect that the surface free energy will be continuous in the zero curvature limit. Sitta et al.¹⁶ have examined this problem using a hard-rectangle fluid both at a circular wall and in confinement within a circular container, showing continuity in the slope of γ vs $C̄$ through the transition $(C̄=0)$ between the confined and the open systems, although only a very limited range of curvatures near zero was studied.

Here, we extend our earlier work to a hard-disk fluid confined within a hard circular container of radius R. The MT equations above will be identical, except the surface free energy is defined slightly differently as $γS≡(ΩS+pVS)/AS$⁠. The motivation for examining this in this 2d system is that the MT expression for γ is simpler (linear in 1/R) and relatively large-size systems can be simulated with modest computational cost. This system has been studied in the past to examine the density variation of the fluid as a function of confining radius,¹⁷ but this is, to our knowledge, the first systematic study of the thermodynamics of hard disks in circular confinement (negative curvature) over a wide range of confining radii. Combined with our previous study for the hard-disk fluid at a circular wall (positive curvature), this study allows us to examine the effects of curvature on the surface thermodynamics in a unified way for both bounding wall and confinement geometries.

In our previous work on the hard-sphere and hard-disk fluids at curved convex walls, we calculated the surface free energy using the Gibbs–Cahn integration method, in which we use the adsorption relation

where ρ_b is the single particle density of the bulk fluid and v_ex is the excess volume, defined as the difference (per unit area) between the volume (V) of a region completely encompassing the interface and the volume of the bulk fluid containing the same number of particles (N) as the interfacial region. The value of γ at a particular pressure (or density) is then found by determining v_ex as a function of pressure (or density) and integrating Eq. (8). For the hard-sphere system, this is simplified because one can assume that γ(P = 0) = 0.

Because of the confinement, the pressure of a hard-disk fluid encapsulated within a hard circle of radius R will differ from that of a bulk fluid with the same chemical potential as a consequence of the Young–Laplace equation and the small system size, where the curvature dependence of the surface free energy can be significant.¹⁸ Because of this, the Gibbs–Cahn method must be modified. Consider a fluid inside a container in equilibrium with a bulk fluid at chemical potential μ, we can define the surface free energy of the interface by

where Ω and Ω_b = −p_bV are the grand potential of the fluid inside the container and that of a bulk fluid with the same volume, chemical potential, and temperature, respectively. Ω is related to the internal energy U of the fluid by the Legendre transform

with differential

The differential for U for this system (including the surface work term) is given by

where p is the pressure inside the pore, which is not the same as the pressure of the bulk fluid at the same μ and T. From Eq. (9), we have

Substituting Eqs. (11) and (12) into Eq. (13) yields

The differentials in Eq. (14) are not all independent, but they are related to one another through the Gibbs–Duhem equation for the bulk fluid,

where the subscript b indicates the value of the variable in the bulk fluid reservoir. The two simultaneous equations [Eqs. (14) and (15)] can be solved using Cramer’s rule by choosing one of the extensive variables X = N, V, or S, giving

where Δp = p − p_b and [Z/X] is defined as

with z_ex defined as the interfacial excess of quantity Z. The choice of X (S, N, or V) eliminates the differential of the intensive variable conjugate to X. In this work, we choose X = N, eliminating the dμ term. We also keep T and V constant, yielding

where ρ_p = N/V and ρ_b = N_b/V_b are the average densities of the fluid in the pore and bulk reservoir, respectively.

Equation (18) can be integrated over p_b to determine the surface free energy as a function of the bulk fluid pressure, or given an equation of state, the packing fraction η_b = πρ_b/4 of the bulk fluid,

Here, the dependence of the bulk pressure on packing fraction can be determined to sufficient accuracy from the equation of state (EOS) of Kolafa and Rottner.¹⁹ 

To determine v_ex [and, thus, γ through Eq. (19)] for the hard-disk fluid confined within a hard circular wall of radius R, we perform a series of Grand Canonical Monte Carlo (GCMC) simulations (constant μ, V, and T) using a modified version of the HOOMD-blue simulation software.^20,21 The software was modified to allow the grand canonical integration scheme to be paired with external potentials. The GCMC simulations require as input the chemical potential (or fugacity) of the bulk fluid, from which the bulk packing fraction can be obtained either through simulation of the bulk fluid or through the application of a high accuracy EOS. To perform simulations of the confined fluid with the same bulk packing fraction as our previous work on positive curvature walls, we apply the Kolafa and Rottner EOS to determine the fugacity by numerical integration

where the compression factor Z = βP/ρ.

The excess volume for the hard-disk fluid/hard circular wall system was calculated directly from average density in the pore using Eq. (18) using GCMC for a range of R from 19/20 = 0.95 to 30 (corresponding to curvatures ranging from −20/19 to −1/30) and bulk packing fractions ranging from 0 to 0.6049. The specific packing fractions were chosen to coincide with those studied for our earlier simulations on the hard-disk fluid/hard circular convex wall system¹⁵ for ease of comparison. (For reference, our highest densities are below the first-order fluid/hexatic phase transition for the hard-disk fluid at about η = 0.7.^22–24) Data for R = ∞ (planar wall) are taken from Ref. 25. The number of GCMC cycles for each simulation was chosen to give similar levels of statistical precision and ranged from 50 × 10⁹ for the smallest cavities at low density to 1.2 × 10⁹ for the largest cavities at high density. The surface free energy γ as a function of η was then obtained through numerical integration of Eq. (19) using the Kolafa–Rottner EOS for the pressure. Simpson’s rule was applied for the numerical integration for all of the even numbered bins in η. For the odd numbered bins (2n + 1), the trapezoid rule was applied to the final bin. For all γ data points, the estimated numerical integration error was smaller than the statistical error.

The simulation results are useful for the validation of theoretical approaches to surface thermodynamics. In this work, we examine two approaches. The first uses a formulation of classical Density Functional Theory (cDFT)¹³ for 2d hard particles by Roth, Mecke, and Oettel (RMO).²⁶ Within cDFT, the grand potential functional Ω[ρ(r)] is a functional of the single particle density ρ(r) with two important properties: (i) The equilibrium density profile ρ₀(r) minimizes the functional; and (ii) at its minimum, the functional reduces to the grand potential of system Ω = Ω[ρ₀(r)]. By minimizing the functional via the variational principle

one obtains the equilibrium structure and thermodynamics at the same time. Hence, one can directly calculate the surface free energy γ and the excess volume v_ex. The shape of the cavity enters the calculation via the external potential and can be arbitrary. This method makes no assumptions as to the functional form of the curvature dependence of v_ex and γ, and, therefore, it can predict potential deviations from the MT form.

Here, we employ a fine 2d Cartesian grid, which allows us to make use of efficient convolutions with the help of fast Fourier transforms (FFT)²⁷ on central processing units (CPUs) and graphics processing units (GPUs)²⁸ for the computation of weighted densities. This also makes it possible to treat the circular cavity and the planar geometry on equal footing. Clearly, the Cartesian coordinates introduce small numerical artifacts in the case of circular cavities. In order to minimize artifacts, we employ a square grid of 2048 × 2048 points, and we have verified that the Gibbs adsorption theorem is satisfied.²⁶ For the circular cavities, considered here, it would be possible to make use of the polar symmetry of the problem and reduce the cDFT problem to an effective 1d one, however, at the cost of losing the possibility to employ standard FFT methods for convolutions. In the effective 1d case, one ends up with a Hankel transform, which can be evaluated using FFT methods when a logarithmic grid is employed.^29,30

The generalized scaled particle theory (gSPT) method used here is similar to that used in our previous paper on the hard-disk fluid at a planar wall,²⁵ modified to apply to a curved confining wall. While strictly speaking, this approach applies directly to a hard-disk fluid outside a hard circular wall (R > 0), we will be comparing it to both the positive and negative R cases. As in the previous paper, the theory begins with the original scaled particle theory (SPT) expression for the excess free energy density Φ, which reads

where ξ₂ = πσ²ρ/4 = η, ξ₁ = πσρ, and ξ₀ = ρ. Here, the scaled particle variable ξ₂ is equivalent to the packing fraction η. This free energy is then generalized as follows:

where f and g are functions of the dimensionless variable ξ₂.

The functions f and g can be specified by requiring agreement with a imposed equation of state via the compressibility factor

which amounts to requiring that

This can be then combined with the usual scaled particle relation

One can readily solve Eqs. (25) and (26) for the functions f and g, leading to

and

To extend this to curved confining walls, we use a binary mixture approach; the excess chemical potential of a single hard disk of radius R added to the system can be obtained as

from which the interfacial free energy for the interface with curvature 1/R is obtained as

Using the functions f and g, this results in

Compressibility factors for the hard-disk fluid include the result from scaled particle theory

which corresponds to f(η) = g(η) = 1 and, hence,

The more accurate Solana EOS³¹ reads

and leads to

with

At very large densities, the SMC2 EOS³¹ with $ZSMC2=ZSol−50000η40(1−η)3$ is even more accurate. The resulting functions f and g, and hence γ, can be calculated analytically, but the expressions are lengthy.

The results for the excess volume v_ex and the surface free energy γ for the confined system as functions of η are plotted in Fig. 1 over the full range of R. To better show the dependence on curvature, we combine the data from the current calculations with those from our previous work for the fluid outside a circular wall¹⁵ to show in Fig. 2 the excess surface volume v_ex a function of the curvature $Ĉ$ over the full range of positive and negative curvatures for selected packing fractions. (Data for all packing fractions and wall curvatures studied can be found in the supplementary material—for clarity, not all packing fractions studied are included in the plots.) Morphometric thermodynamics predicts that v_ex is linear in the curvature. This is seen in Fig. 2 to hold except for high negative curvature at the highest packing fractions.

At this point, it is useful to note that, in this work, we define R to represent the actual radius of the circular wall. An alternative definition that is frequently seen uses a radius R′, defined by the point of closest contact of the hard particle to the circular wall.³² For the 2d system here, this gives R′ = R ± 1/2, with the “−” sign for the confined system and the “+” for solvated disk. The data from this work can be converted to the alternate reference using the definition of v_ex. If we denote the excess volume in the alternate reference as $vex′$⁠, we have

Because the grand potential Ω is a well-defined property of the system, independent of the definition of the radius, we obtain

where the symbol “′” denotes surface thermodynamic quantities measured using the R′ radius definition. The first non-Hadwiger coefficient γ₂, however, is invariant to the choice of radius definition. We have chosen the radius convention used here and in our earlier work because it gives a value for γ₀ that is positive and is easier to extend to hard-disk fluid mixtures of different-sized disks.

To see the deviations from MT at high η more clearly, Fig. 3 shows the excess volume vs curvature for the highest packing fractions studied. At the very highest values of η, the slope of the v_ex vs $C̄$ curve changes sign for negative curvatures and v_ex exhibits a minimum at $C̄=0$⁠. These significant deviations from the MT ansatz at high η are due to the presence of packing constraints in the confined systems that lead to a density maximum far smaller than close packing in the bulk hard-disk fluid. In other words, at high densities the ratio ρ_p/ρ_b in Eq. (18) changes from an increasing function of η to a decreasing function of η due to a plateau in η_p as maximum packing in the pore is approached.

Figures 2 and 3 also show the predictions from cDFT in comparison to the simulation data. The agreement at low to intermediate packing fractions is excellent. At the higher packing fractions (Fig. 3), the Density Functional Theory (DFT) results for v_vex show some deviation from the simulation results for the positive curvature systems; however, the agreement for negative curvatures is remarkably good, with DFT successfully reproducing the nonlinear deviations from the MT predictions.

The results for the surface free energy, γ, as a function of curvature for various packing fractions, are shown in Fig. 4, combined with the data from our previous work for positive curvature.¹⁵ For clarity, the same quantities for the highest packing fractions studied are shown in Fig. 5. Moreover, the solid lines in Figs. 4 and 5 show the results for γ from cDFT, which agree extremely well with the simulation results, even up to high negative curvature and high packing fractions. The better agreement between DFT and simulation for γ vs that for v_ex can be understood by noting that γ is obtained by integration from the excess volume, which decreases the relative differences at high η. The surface free energy follows the MT prediction at low η to intermediate η; however, at the higher packing fractions studies, there is significant deviation from linearity, especially at high negative curvature, due to the effects of confinement. For R = −3 (i.e. $C̄=−1/3$⁠), the pore is six fluid particle diameters across and can accommodate fewer than 32 particles at maximum packing. The additivity condition of Hadwiger’s theorem implies that when correlations in the fluid are on the same length scale as the wall, the geometric foundation of MT no longer holds. However, Fig. 5 shows that even for a highly confined fluid at a high density, where correlations with the wall extended throughout the cavity, MT has strong predictive power.

At high packing fractions, the DFT results show an unusual maximum in between the highest two negative curvatures studied by the simulation: R = 0.95 and $1.07735;C̄≈1.0526$ and 0.9282, respectively. These two curvatures correspond to system sizes in which only one and two hard disks, respectively, can fit within the circle at maximum packing. For these systems, the grand partition function can be calculated exactly,³³ and, therefore, we have exact results for γ for the two largest negative curvatures,

where f = e^βμ is the fugacity and R′ = R − 1/2. For these two highest negative curvatures, the simulation results agree with the exact expressions above within the error bars (see the comparison table included in the supplementary material). In addition, the agreement with the DFT results, while not perfect, is quite good and shows that the maximum predicted by the DFT is real and not an artifact (Fig. 6).

To examine the deviations from MT, we have performed a least-squares fit of the values of γ at each packing fraction to a quadratic in $C̄=1/R$ separately for positive (convex) and negative (concave) curvatures. For $C̄≥0$⁠, we use the entire range of curvatures in the fit; however, because of the strong deviations at high curvature for the concave systems, we use only those curvatures in the interval [−1/3, 0] (five data points). The results for γ₀ are presented in Fig. 7 along with the predictions from both SPT and gSPT [Eqs. (33) and (36), respectively]. The results for γ₀ for the convex and concave systems agree within the statistical error of the calculation, so only the data for the convex system are shown. The SPT expression underestimates the value of γ₀, especially at high packing fractions; however, the agreement for gSPT with the simulations is excellent over the entire range. The results for the linear and quadratic coefficients—γ₁ and γ₂, respectively—are shown in Fig. 8. Here, we see that the predictions of MT, i.e., γ₂ = 0 and γ₁(R < 0) = γ₁(R > 0) hold up to packing fractions of about 0.45, showing a surprisingly large range of validity of MT even for moderate negative curvatures and moderate packing fractions.

The work presented here presents the most comprehensive study of the surface free energy γ of curved surfaces through the zero-total-curvature limit to the knowledge of the authors. Using Monte Carlo simulation, classical Density Functional Theory (cDFT), and a generalized scaled particle theory (gSPT), we find that in the hard-disk fluid at circular hard walls, the predictions of morphometric thermodynamics are highly accurate over a large range of fluid conditions and wall sizes, extending from highly confined fluids to solvated walls on the same size scale as the fluid particles. This indicates that even under conditions where the concavity and additivity requirements of Hadwiger’s theorem are clearly violated, the MT framework has the potential to provide strong predictive power in assessing surface free energies. At high densities and high negative curvatures, both the excess volume (v_ex) and γ exhibit significant deviations from the predictions of MT due to strong correlations in these dense highly confined systems. The results for γ calculated using cDFT show excellent agreement with the simulation results at all curvatures and packing fractions and both the simulation and cDFT show agreement with the exact values of γ for highly confined systems that can contain at most one or two particles. In addition, gSPT is successful at accurately predicting γ for packing fractions and curvatures where MT works well.

Included in the supplementary material is a table of the raw data for the excess volume v_ex and the surface free energy γ from the GCMC simulations for all pore sizes R and packing fractions studied. In addition, for the two smallest pores containing a maximum of one or two disks, respectively, a table is included that directly compares the simulation and cDFT values of the surface free energy γ with the results for γ from the exact partition functions.

B.B.L. and S.C.M. acknowledge support from the National Science Foundation for partial funding of this work under Grant No. CHE-1465226.

The authors have no conflicts to disclose.

The data that support the findings of this study are available within the article and its supplementary material.