Random Surface Statistical Associating Fluid Theory:
Adsorption of nAlkanes on Rough Surface
Abstract
Adsorption properties of chain fluids are of interest from both fundamental and industrial points of view. Density Functional Theory (DFT) based models are among the most appropriate techniques allowing to describe surface phenomena. At the same time Statistical Associating Fluid Theory (SAFT) successfully describes bulk PVT properties of chainfluids. In this publication we have developed novel version of SAFTDFT approach entitled RSSAFT which is capable to describe adsorption of short hydrocarbons on geometrically rough surface. Major advantage of our theory is application to adsorption on natural roughs surfaces with normal and lateral heterogeneity. For this reason we have proposed workflow where surface of real solid sample is analyzed using theoretical approach developed in our previous work Aslyamov and Khlyupin (2017) and experimentally by means of low temperature adsorption isotherm measurements for simple fluids. As result RSSAFT can predict adsorption properties of chain fluids taking into account geometry of the surface sample under the consideration. In order to test our workflow we have investigated hexane adsorption on carbon black with initially unknown geometry. Theoretical predictions for hexane adsorption at 303K and 293K fit corresponding experimental data well.
pacs:
Valid PACS appear hereI Introduction
Confined fluid thermodynamic properties become crucial when one deals with oilgas production from unconventional reservoirs, where the fluid often is contained in pores of nanoscopic scale. It is well known that adsorption phenomena play crucial role in capacity characterization of nanoporous materials. For this reason theoretical model which is capable to predict adsorption properties of chain fluid is highly desirable Liu et al. (2017). Also industrial needs induce high interest to theory accounting for real adsorbate surface geometry. Indeed, all natural materials are geometrically heterogeneous and nanoroughness has strong influence on surface phenomena Aslyamov and Khlyupin (2017). Thus, in this publication we have proposed theoretical approach describing chainmolecules fluid adsorption on surface with real rough geometry. As result PVTproperties of confined fluid can be described in framework of this theoretical method.
One of the most useful theoretical approaches to describe thermodynamic properties of real fluid is Statistical Associating Fluid Theory (SAFT) based EOS. Among the variety of known SAFT versions we could emphasize VRSAFT due to explicit connection with molecular structure and various successful applications to hydrocarbons modeling. Moreover, the model is defined by a set of fluid molecular parameters such as size and characteristic interaction energy. This fact and a wide application motivate the study of PVTproperties for various chain fluids by means of SAFT. In this work we have considered SAFTVR version for Mie potential, developed by Lafitte et. all Lafitte et al. (2013). This model demonstrates excellent match with experiment for wide range of hydrocarbons Lafitte et al. (2013).
Overwhelming part of publications considering chain molecules in terms of SAFT are dedicated to homogeneous bulk case. However, nonuniform spatial distribution of confined fluids can not be described using bulk EOS only and inhomogeneous extension is needed. In order to take into account impact of spatial boundaries on fluids PVT properties Density Functional Theory (DFT) is often used. Indeed, in literature one can find a variety of surface phenomena which was described in terms of DFT. Despite the fact, that description of simple fluids using DFT has a long history and now is a standard practice, the extension to the case of SAFTEOS is an actual problem Liu et al. (2017); MartínezRuiz et al. (2017); Schindler et al. (2013). Rigorous way to account for chain structure of molecules in DFT approach is developed in work Yu and Wu (2002). However this model does not consider fluidfluid interaction which is necessary for accurate description of thermodynamic properties even in bulk case. Molecular interaction can be described in terms of VRSAFT GilVillegas et al. (1997) where central role plays radial distribution function of hard sphere fluid. This bulk model was used as basis in order to develop inhomogeneous SAFT version describing interfacial behavior of chain molecules. In work Schindler et al. (2013) authors considered adsorption phenomenon using one of the most popular version of DFT and SAFTVR. Results of this work demonstrate qualitative agreement with computer simulations, however authors did not consider adsorption isotherms and corresponding experimental data. Later DFTSAFT approach was verified in comparison with experimental data in the case of monomer fluids Malheiro et al. (2014). In these publications Schindler et al. (2013); Malheiro et al. (2014) Radial Distribution Function (RDF) was calculated from result Chang and Sandler (1994) where analytical expression was obtained, however this form of RDF does not provide correct expression for density derivative which is needed in DFT calculations. Alternative way was developed in work Franco et al. (2017) where meanvalue theorem was applied in order to calculate confinement term in Helmholtz free energy. This approach contains an assumption of a squarewell potential for the fluidsolid interactions. This fact is serious restriction for further comparison with experiments for natural nonideal materials, since this simplified potential is unable to describe geometrical roughness.
In a way similar to the approaches discussed above, we have used VRSAFT Lafitte et al. (2013) as the basis. In another words, our version of SAFTDFT entitled RSSAFT (RS means Random Surface) contains VRSAFT Lafitte et al. (2013) as homogeneous bulk limit (far enough from the geometrical boundaries). At the same time RSSAFT has two major advantages. Firstly, we have developed novel analytical form for fluidfluid interaction contribution. This expression matches with very popular perturbative expansion from work Lafitte et al. (2013) in bulk and can be extended for inhomogeneous fluid. Also this explicit form is derived in terms of Lambert special functions that allows to obtain analytical expression for density derivative. Secondly, we have used our previous results Aslyamov and Khlyupin (2017) to account for real surface geometry at nanoscale. Thus, RSSAFT describes adsorption properties of fluid with accurate bulkEOS taking into account surface nanoroughness both in normal and lateral directions.
For wide applications it is necessary to link theoretical predictions with characteristics of real solid samples. Low temperature adsorption analysis is widely used for characterization of pore and surface structures Landers et al. (2013). Theoretical model can be a part of workflow containing analysis of simple gas adsorption at low temperature and further predictions for more complex fluid adsorption on the same solid sample. For example, in Mitchell et al. (2015) authors obtained pore size distribution (PSD) for activated carbon from analysis of nitrogen adsorption at 77K and then calculated nalkane fluid adsorption taking into account obtained PSD. For adsorption calculations the authors used model from Schindler et al. (2013) which was developed for inhomogeneous fluid near ideal smooth surface.
In the current publication we propose the workflow allowing to predict adsorption properties of nalkanes on natural surface. The major point of this approach is surface geometry analysis of adsorbent sample. In publication Aslyamov and Khlyupin (2017) two of us shown that geometrical model can be obtained from adsorption isotherms measurements. The most appropriate way for this procedure is low temperature simple gas adsorption. Thus, experimentally measured adsorption isotherm for simple fluid on adsorbent sample is used as input data for our workflow. In the framework of RSDFT it is possible to fit experimental data and to construct adsorbent surface model corresponding to geometrical roughness of the solid sample. Finally, presented in the current manuscript theory RSSAFT is capable to describe nalkanes adsorption properties taking into account obtained surface characteristics. Thus, RSSAFT workflow allows to predict surface phenomena for chain fluids on natural material using routine experimental data with simple fluids. Example of this workflow can be found in Fig. 1. As one can see experimental data for low temperature nitrogen adsorption is used in order to tune RSDFT and to obtain solid surface model. As result RSSAFT predicts hexane adsorption isotherms on the solid sample using obtained geometry characteristics.
In order to test our theoretical approach we have investigated hexane adsorption on carbon black with initially unknown surface geometry. In accordance with our workflow we have experimentally measured adsorption isotherm for nitrogen at 77K. We have used RSDFT to obtain geometry roughness characteristics in normal and lateral directions. Finally, we have performed experiments on hexane adsorption on the same solid sample at several temperatures. The comparison demonstrates that RSSAFT predictions fit experimental results well.
Ii Inhomogeneous version of SAFT
The general idea is modification of known bulk SAFT Helmholtz free energy using certain version of DFT as was done previously Liu et al. (2017); MartínezRuiz et al. (2017); Schindler et al. (2013). As we discussed above in our study VRSAFT plays the role of the basic model corresponding to homogeneous limit. More precisely, our bulk fluid corresponds to onecomponent system of molecules with chainstructure consisting of spherical segments with diameter . Molecular interaction is defined by monomer fluidfluid Miepotential between the segments:
(1) 
where is the distance between the segments, is the potential depth, and are attractive and repulsive exponents, respectively. In the case of bulk fluid this system was successfully described by VRSAFT in work Lafitte et al. (2013). In accordance with SAFT formalism Helmholtz free energy can be expressed in terms of ideal, monomersegment and chain contributions as
(2) 
In formulation of ideal contribution we have followed to work Yu and Wu (2002) where ideal part contains translational degree of freedom and takes into account chain connectivity via bonding potential. This expression accounts configuration of chain in three dimensional space which is defined by a set of coordinates . Ideal term accounts direct effect of connectivity while chain term corresponds to indirect interactions due to the excluded volume effects. Chain contribution is defined from Lafitte et al. (2013) using inhomogeneous density distribution.
In accordance with formulation of VRSAFT Lafitte et al. (2013) we have used thermodynamic perturbation theory for monomer contribution :
where is reference system, and are the first two perturbations terms. It is possible to map the reference fluid to hard spheres system with another effective diameter Barker and Henderson (1967a):
Thus, reference system corresponds to gas of hard spheres with new diameter . For inhomogeneous hard sphere contribution we have used result of Fundamental Measure Theory Roth (2010) which equals to corresponding term from Lafitte et al. (2013) at bulk limit.
In order to calculate perturbation terms authors of Lafitte et al. (2013) used BarkerHenderson theory, where the first two perturbation terms are defined using hard sphere radial distribution function RDF only. Thus, explicit expressions for density of Helmholtz free energy and ( is the number of molecules) have the following form:
(3)  
The second expression in (3) is different from classical one Barker and Henderson (1967b), because additional multiple is applied Lafitte et al. (2013), where is correction prefactor which depends on density of fluid (we have used expression from Lafitte et al. (2013)), is isothermal compressibility of reference system (HS fluid), and are the number segment density and the pressure of HS fluid, respectively. Parameter is function of dimensionless density and can be calculated from the Carnahan Starling compressibility Carnahan and Starling (1969)
It is possible to extend results for the first two perturbative terms to the case of inhomogeneous fluids, following the ideas proposed by Toxvaerd Toxvaerd (1981). At the same time inhomogeneous case contains serious issues related to the absence of explicit representation for inhomogeneous RDF . However a lot of studies note that there is adequate approximations allowing to consider bulk RDF at averaged density , which has the following expression:
Thus the first two perturbation terms at the case of inhomogeneous fluid can be represented as the follows:
(4)  
For further work it is important to note the two points: 1) contrary to the bulk case approximations which were performed in Lafitte et al. (2013) can not be used here due to abrupt oscillations of near the walls; 2) DFT approach for chain molecules Yu and Wu (2002) demands knowledge of functional derivatives of Helmholtz free energy. For these reasons in this work we have used novel form of RDF allowing to obtain the potentials (4) and the density derivative as analytical expressions.
RDF is obtained in PercusYevick (PY) assumptions and corresponding Verlet Weis (VW) modification have the following form:
(5) 
(6) 
detailed description of density functions ,, ,, , , and system parameter , where can be found in Appendix A. This result (5) was obtained as direct calculation of famous Wertheim’s solution of OrnsteinZernike equation in PY approximation. The final expression was modified in accordance with VerletWeiss corrections. As one can see from Fig. 2 expression (5) fits simulations results well. Key difference of this RDF form and alternative published versions is that potentials (4) and their density derivatives have analytical expressions.
In order to demonstrate appropriateness of our form of RDF for calculations of perturbations terms we started from bulk case. Integrals (3) can be calculated analytically using RDF (5), analytical results can be found in Appendix A. As one can see from Fig. 2 these perturbation terms at bulk case have good agreement with numerical results from work Lafitte et al. (2013).
It is standard to describe inhomogeneous fluids by means of grand canonical thermodynamic potential :
(7) 
where is the external potential, is the chemical potential, is chain configuration density (dependence on is crucial difference from the standard DFT approach). Our interest is focused on description of fluids near solid wall, for this reason external potential is defined by interaction between solid molecules and fluid segments. Equilibrium density distribution corresponds to the solution of the following equation:
(8) 
as was demonstrated in work Yu and Wu (2002) in the case when terms of and are depended on segment density solution of (8) can be written as function of . In assumption of surface radial symmetry external potential can be represented as function of one spatial coordinate only , where is the distance between the chainsegment and the surface Khlyupin and Aslyamov (2017).
(9) 
where and function is determined from recurrence formula:
(10) 
with and . For numerical calculations it is more convenient to consider , where is segment bulk density which satisfies to equation (9) at the limit . Let us consider (9) far enough from the walls that , that is equivalent to , then one can obtain the following equation:
(11) 
After substituting of this result into equation (9) the following equation can be written:
(12) 
This equation can be solved by method of simple iteration.
Iii Surface geometry and external potential
It already was assumed that the system has radial symmetry and density distribution depends on zcoordinate only. In terms of external potential this situation is equivalent to . Such model allows to investigate gas adsorption on surface of amorphous materials.
Usually natural materials are rough at nanoscale and surface geometry significantly influences the properties of confined fluids. One of the most modern model describing real surface geometry is correlated Gaussian random process Herminghaus (2012). In the frame of this approach random rough surface can be characterized by two natural parameters:

Parameter corresponds to variance of the fluctuating surface roughness height in the normal direction.

Parameter (the correlation length) determines the characteristic scale along the lateral direction of the correlation function decay. For example, the correlation length of white noise equals zero, which reflects the independence of heights at any two points in the lateral plane.
In work Khlyupin and Aslyamov (2017) two of us obtained fluidsolid potential for rough surface using correlated random model. Our approach is based on Free Energy Averaging Technique presented in the work Forte et al. (2014). Also advanced theory of Markovian random processes and the first passage time probability problem were applied as part of averaging procedure (see work Khlyupin and Aslyamov (2017) for more details). Proposed model may be applied for wide range of correlation functions of the random solid surfaces. As result potential accounts for surface roughness via and depends on zcoordinate only. Thus, after this modification density distribution equation (9) is still applicable. Also surface geometry restricts configurational space available for fluid molecule. This fact modifies result of integration . One can calculate this twodimensional integral taking into account only permitted domains for fluids at a certain level as
(13) 
where is the total area and is the part of which is free from solid media at level . One can find this area from the following expression:
(14) 
Adsorption isotherm in the case of rough surface with known parameters can be calculated using external potential Khlyupin and Aslyamov (2017) and expression (13):
(15) 
Iv Results
In this section we have applied the developed theory for consideration of gas adsorption on surface of carbon black using developed theory. Despite carbon black is often used in adsorption investigations its surface geometry is not clear. A lot of publications claim nonideal structure of surface at nanoscale. However none of these studies defined the roughness based on direct measurements. For this reason surface characterization of the adsorbate sample is serious challenge.
Our major aim is describe nalkane adsorption on certain carbon sample. In order to make desired theoretical predictions taking into account real surface geometry we have proposed the workflow containing the following two steps:

Considered solid sample is investigated experimentally by simple fluid (argon, nitrogen) adsorption at low temperature. RSDFT allows to obtain geometrical characteristics of rough surfaces form adsorption isotherms analysis, both hight variance (roughness in normal direction) and correlation length (lateral characteristic).

Developed theory RSSAFT can use obtained from RSDFT characteristics of rough surface geometry in further theoretical predictions of complex (chain) fluid adsorption properties.
As result the theory can provide accurate PVTproperties in the bulk phase for wide range of nalkanes and also take into account roughness of adsorbate sample. Finally in this section one can compare theoretical predictions with experimental data for nalkane adsorption.
iv.1 Analysis of surface geometry
Simple gas adsorption on ungraphitized carbon black was investigated in work Aslyamov and Khlyupin (2017) where RSDFT analysis demonstrated significant geometrical roughness of the surface. In accordance with RSDFT approach, information about surface geometry can be obtained from adsorption isotherms analysis. More precisely, smooth form (without steps) of the isotherm is related to geometrically rough structure of surface. Thus, the best fit of experimental data by RSDFT calculations allows to define roughness parameter and lateral characteristic . This information is sufficient to construct surface profile corresponding to investigated solid sample.
We have performed experimental measurements of surface geometry using ultrahighpurity nitrogen (99.999%) as the adsorbate. Experimental adsorption isotherms were measured at 77K (Fig. 3) using an ASAP2020 volumetric adsorption analyzer from Micromeritics, USA and allowed to obtain geometry characteristics of the carbon surface. Carbon black powder provided by SigmaAldrich was used as representative carbon material for both theoretical and experimental study. Prior the isotherm measurements, the sample was subjected to degassed under vacuum at elevated temperatures (473K) overnight (12 h) to remove any adsorbed species. The BrunauerEmmettTeller (BET) equation was applied to obtain specific surface area (SSA). For studied carbon black sample it was 1.6 . The saturation pressure of nitrogen at a liquid nitrogen temperature was determined every two hours during the experiment using a vapor pressure thermometer.
iv.2 nAlkane adsorption on real rough surface
We have performed experimental measurements of nalkane adsorption on the same solid sample as in the case of nitrogen at several temperatures. The highpurity hexane (99.0%) vapor adsorption isotherms were determined by ASAP2020 using vapor adsorption option. Hexane is placed in a special vapor tube and degassed prior the adsorption measurements. The set of vapor adsorption isotherms on carbon black surface was measured at 293K and 303K (Fig. 4). The recirculation bath (recirculating chiller with thermostat Julabo FP50) was used to control the temperature of the sample tube along its whole length.
Our developed version of SAFT entitled RSSAFT allows to use geometrical parameters which are obtained from nitrogen adsorption analysis. Thus, predictions of the theory can be compared with experimental data for hexane. As was noted above the bulk limit of our version of SAFT coincides with very popular SAFTVR. This fact allows us to use published parameters Garrido et al. (2017) for fluid description Table 1. Thus, our version of SAFT can describe bulk PVT properties of hydrocarbons with good accuracy.
M  , K  , Å  , K  , Å  

Hexane  2  376.35  4.508  6  19.26  67.2  3.954 
In order to minimize the number of free parameters we have fixed solidfluid characteristic diameter in accordance with LorentzBerthelot rules Maitland et al. (1981) , where is commonly used value of carbon molecule diameter. Geometrical parameters of surface are fixed too and are defined from nitrogen adsorption RSDFT analysis. Thus, we have used only one tuning parameter which is characteristic solid fluid energy .
Comparison of experimental results and RSSAFT calculation can be found in Fig. 4. As on can see the theory fits experimental data well, especially at small pressures. Also the theoretical predictions correctly reflect temperature changes: the start of adsorption isotherm in terms of relative pressure is significantly shifted to the right with the temperature increase. The most noticeable deviations from experiments the theory demonstrates in vicinity of saturation pressure . At these conditions theoretical adsorption isotherms are steplike that is not confirmed by experiments. According to the literature steplike adsorption isotherms relay to well defined structure of adsorbed fluid. Indeed, in our theory all adsorbed chain molecules have parallel orientation to the surface level. On the other hand, pressures near correspond to significant fluid fluctuations which can broke well defined layering. This is exactly what one can see from experimental adsorption isotherms. Despite difference in form of isotherms at pressures near theoretical results have a good quantitative agreement with experiment data.
V Conclusion
Crossover of DFT approach and SAFT model for chain structured molecules is one of the most perspective way to describe chain fluids surface phenomena. In this publication we have proposed novel version of DFTSAFT approach entitled RSSAFT. The key differences of our model are concluded in alternative analytical form for intermolecular term of Helmholtz energy and random surface approach which allows take into account geometrical roughness in both normal and lateral directions.
Developed theory can be applied for prediction of short hydrocarbons adsorption properties on the surface of certain solid sample. We have formulated workflow where information about geometrical roughness is obtained from experiment and RSDFT calculations for low temperature adsorption of simple fluids. As was demonstrated in work Aslyamov and Khlyupin (2017) RSDFT allows to construct model of rough surface for certain solid sample which is defined by two geometrical parameters. These surface characteristics are used in further RSSAFT calculations for nalkane adsorption on the same solid sample.
In oder to demonstrate capability of our approach we have compared RSSAFT predictions with experiment. We have considered hexane adsorption on carbon black sample with initially unknown surface geometry. We have used the same carbon black solid sample in all experiments which conclude: nitrogen adsorption at 77K, hexane adsorption at 293K and 303K. Theoretical RSDFT analysis demonstrates rough geometry of sample. RSSAFT predictions for considered surface fit experimental data well.
Appendix A Hard Sphere RDF
The initial description of HS fluid was performed by Wertheim Wertheim (1963) and Thile Thiele (1963). They obtained the solution of OrnsteinZernike equation in PercusYevick (PY) approximation Percus and Yevick (1958). However, analytical result was obtained only for the Laplace transform of product
(16) 
Wertheim Wertheim (1963) obtained Laplace image as analytical function. According to him RDF has the following form:
(17) 
where, is point on the real coordinate of the complex plane, such that is greater than the real part of all singularities of the integrand.
Here integral (17) is calculated by direct method using residue theorem of complex analysis. For determination of singularity points it is necessary to solve following equation in respect to variable (denominator of (17) equals to zero):
(18) 
After substitution of expressions (A), it transforms into transcendental equation for variable :
(19) 
Equation (A) has infinite number of roots on complex plane. Let us start with obvious root which is pole of the third rang, also it is unique real solution of (A). The other roots are conjugated complex simple poles , where , are real and imaginary parts of complex number. Thus, in accordance to residue theorem, expression (17) can be rewritten as:
(20) 
where is derivative with respect to at the point , here the first term “1”corresponds to the residue at point , the second term is sum over all simple complex poles. Simpler expression can be obtained after summing conjugated poles:
where
Expression (A) contains only real functions which are depended on . Thus, as one can see from (A), in order to calculate RDF one needs only the distribution of roots .
Let us consider new equation which is the limit of equation (A):
(23) 
By introduction of a new variable it is possible to rewrite (23) in more simple form: . Such equation can be solved exactly in terms of Lambert functions Scott et al. (2006); Corless et al. (1996)
After simple modifications, the above equation can be written as . Thus, solution of equation (23) has the following form
(24) 
where enumerates complex branch of Lambert function.
Using exact solution (24) as the limit, the solution of (A) can be written as series of :
(25) 
where coefficients depend only on density and can be found after substitution of (25) in (A). For the aims of this work it will be enough to consider only six terms in the sum (25). Corresponding coefficients have the following explicit expressions:
(26)  
Thus zeros of (A) are defined by the following partition sum:
(27) 
The next question is summation over (index of zeros) in expression (A). It is easy to verify, that for all poles , and , then contribution of exponential nth term in sum (A) rapidly decreases, when the number increases. Thus, for accurate result summation over all terms in (A) is not required and (A) can be rewritten as analytical expression with terms (in our study we have used ):
(28) 
Despite wide applications, PY approximation has two weak points: the contact value is too low at high density; phase of oscillation at large distance differs from the exact one. In order to correct these artifacts construction of VerletWeis can be used Verlet and Weis (1972). This is achieved by introduction of a modified density and modified HS diameter in order to correct RDF oscillation. Verlet and Weis, also, proposed an additional term which improved contact value , so corrected RDF has following form:
where parameters and can be found from
The form of expression of added term in (A) coincides with analytical result (A). This fact helps the process of further calculations of corrected form (A).
References
 Aslyamov and Khlyupin (2017) T. Aslyamov and A. Khlyupin, The Journal of chemical physics 147, 154703 (2017).
 Liu et al. (2017) J. Liu, L. Wang, S. Xi, D. Asthagiri, and W. G. Chapman, Langmuir 33, 11189 (2017).
 Lafitte et al. (2013) T. Lafitte, A. Apostolakou, C. Avendaño, A. Galindo, C. S. Adjiman, E. A. Müller, and G. Jackson, The Journal of chemical physics 139, 154504 (2013).
 MartínezRuiz et al. (2017) F. J. MartínezRuiz, F. J. Blas, A. I. M.V. Bravo, J. M. Míguez, and L. G. MacDowell, Physical Chemistry Chemical Physics 19, 12296 (2017).
 Schindler et al. (2013) B. J. Schindler, L. A. Mitchell, C. McCabe, P. T. Cummings, and M. D. LeVan, The Journal of Physical Chemistry C 117, 21337 (2013).
 Yu and Wu (2002) Y.X. Yu and J. Wu, The Journal of chemical physics 117, 2368 (2002).
 GilVillegas et al. (1997) A. GilVillegas, A. Galindo, P. J. Whitehead, S. J. Mills, G. Jackson, and A. N. Burgess, The Journal of chemical physics 106, 4168 (1997).
 Malheiro et al. (2014) C. Malheiro, B. Mendiboure, F. Plantier, F. J. Blas, and C. Miqueu, The Journal of chemical physics 140, 134707 (2014).
 Chang and Sandler (1994) J. Chang and S. I. Sandler, Molecular Physics 81, 735 (1994).
 Franco et al. (2017) L. F. Franco, I. G. Economou, and M. Castier, Langmuir 33, 11291 (2017).
 Landers et al. (2013) J. Landers, G. Y. Gor, and A. V. Neimark, Colloids and Surfaces A: Physicochemical and Engineering Aspects 437, 3 (2013).
 Mitchell et al. (2015) L. A. Mitchell, B. J. Schindler, G. Das, M. C. dos Ramos, C. McCabe, P. T. Cummings, and M. D. LeVan, The Journal of Physical Chemistry C 119, 1457 (2015).
 Barker and Henderson (1967a) J. A. Barker and D. Henderson, The Journal of Chemical Physics 47, 4714 (1967a).
 Roth (2010) R. Roth, Journal of Physics: Condensed Matter 22, 063102 (2010).
 Barker and Henderson (1967b) J. A. Barker and D. Henderson, The Journal of Chemical Physics 47, 2856 (1967b).
 Carnahan and Starling (1969) N. F. Carnahan and K. E. Starling, The Journal of Chemical Physics 51, 635 (1969).
 Toxvaerd (1981) S. Toxvaerd, The Journal of Chemical Physics 74, 1998 (1981).
 Khlyupin and Aslyamov (2017) A. Khlyupin and T. Aslyamov, Journal of Statistical Physics 167, 1519 (2017).
 Herminghaus (2012) S. Herminghaus, Physical review letters 109, 236102 (2012).
 Forte et al. (2014) E. Forte, A. J. Haslam, G. Jackson, and E. A. Müller, Physical Chemistry Chemical Physics 16, 19165 (2014).
 Garrido et al. (2017) J. M. Garrido, M. Cartes, and A. Mejía, The Journal of Supercritical Fluids 129, 83 (2017).
 Maitland et al. (1981) G. Maitland, M. Rigby, E. Smith, and W. Wakeham, “Intermolecular forces,” (1981).
 Wertheim (1963) M. Wertheim, Physical Review Letters 10, 321 (1963).
 Thiele (1963) E. Thiele, The Journal of Chemical Physics 39, 474 (1963).
 Percus and Yevick (1958) J. K. Percus and G. J. Yevick, Physical Review 110, 1 (1958).
 Scott et al. (2006) T. C. Scott, R. Mann, and R. E. Martinez Ii, Applicable Algebra in Engineering, Communication and Computing 17, 41 (2006).
 Corless et al. (1996) R. M. Corless, G. H. Gonnet, D. E. Hare, D. J. Jeffrey, and D. E. Knuth, Advances in Computational mathematics 5, 329 (1996).
 Verlet and Weis (1972) L. Verlet and J.J. Weis, Physical Review A 5, 939 (1972).