Gas Sorption Characterization of Porous Materials Employing a Statistical Theory for Bethe Lattices

In the present work, a recently developed statistical theory for adsorption and desorption processes in mesoporous solids, modeled by random Bethe lattices, has been applied to obtain pore size distributions and interpore connectivity from sorption isotherms in real random porous materials, employing a robust and validated methodology. Using the experimental adsorption–desorption N2 isotherms at 77.4 K on Vycor glass, a porous material with random pore structure, we demonstrate the solution of the inverse problem resulting in extracted pore size distribution and interpore connectivity, notably different from the predictions of earlier theories. The results presented are corroborated by the analysis of 3D digital images of reconstructed Vycor porous glass, showing excellent agreement between the predictions of geometric analysis and the new statistical theory.


Algorithmic details on extracting structural properties from sorption isotherms
To seek a solution of the volume Pore Size Distribution (PSD), φ(x), in the adsorptiondesorption integral equations, we must first discretize it setting, () ≅ , using a certain number of discrete pore sizes, called nodes.Hence, we must first set the limits of integration in pore size, x, from [0, ∞] to a finite region of [  ,   ], with xmin and xmax, being the respective lower and upper pore sizes of the PSD that we want to extract.Evidently, these values are not known a priori, and a trial-and-error procedure is needed to improve the accuracy and resolution of the extracted PSD.

S2
The above equations can be written for a set of Ma pressure points for adsorption, and Md pressure points for desorption, with a total number of pressure point, M=Ma+Md.N is the number of discrete cylinder sizes for which we have input kernels.Thus, we employ the following general set of equations, for i=1,M and j=1,N (N < M ): Or equivalently, we need to minimize: with, and, Parameters   (𝑎𝑑𝑠) ,   (𝑎𝑑𝑠) ,   (𝑎𝑑𝑠) ,   (𝑑𝑒𝑠) ,   (𝑑𝑒𝑠) Similarly, function  ,, =  , (, , , ), is calculated from: ] From the expressions in eqs.(S4) through (S8) it is evident that the matrices K (ads) and K (des) have a complex dependence on relative pressure, p/p0, as well as on the unknown variables, z, f, and φ (through  , ,  , ).Hence, eq. ( S2) is non-linear in terms on the unknown variables and must be solved employing a triple iterating scheme to determine f, z, and φ as follows: In the first iteration loop, in f, at the k th iteration, f=fk, while in the second iteration loop, in z, at the l th iteration,  =   .Then we solve eq.(S2) for f=fk,,  =   , to get the volume PSD, φ that minimizes, the residual,   ().Finally, we come with a set of optimal solutions in PSD, φ, among which we select the one that gives the minimum, value of   () for a variety of values of parameters,  ∈ [0, 1] , and  ∈ [2, ∞).The above procedure is illustrated schematically in Figure S1: Note that the solution eq.(S2) for a fixed set of values, f=fk,,  =   , to get the volume PSD, φ is not straightforward, because this equation is non-linear in φ.Hence, we choose to solve eq.(S2) iteratively, as follows: 1. We start by assigning an initial guess for φ employing a random number generator from a uniform number distribution in (0,1).
2. Accordingly, we compute K (ads) , K (des) from eqs. (S3-S8) and solve the simple linear east square problem: 3. Then we update to the next iteration using a simple relaxation scheme to avoid the solution from diverging:

S6
Note that in eq.(S9), we use w=0.5, which is conservative but ensures convergence down to practically machine accuracy.
We repeat steps 2-3, until we achieve converge in φ based on the following criterion: A routine solving linearly constrained linear least-squares problems based on a twophase (primal) quadratic programming method has been implemented.The routine has been successfully employed in extracting PSDs from sorption isotherms in SDCM structures 1 .Issues regarding convergence, sensitivity, and uniqueness of solution, have been considered and checked in a similar fashion with our previous work on SDCM structures 1 .Stringent convergence tolerances have been imposed (accepted error tolerance Ltol~1x10 -12 ) and up to 500 iterations have been required in some cases.In all cases the procedure converged to the same final solution for the extracted PSD starting from different initial guesses using the procedure described in step 1,

Solution stabilization and smoothing using cubic B-splines
A major drawback of the above methodology is that the number of nodes in pores sizes, N, must be much smaller than the number of points in pressure, M, to get a stable PSD solution 2-3 , a constraint that limits the resolution of the PSD and thus the accuracy of the solution.The above issue is resolved by considering the solution () to be represented by piecewise continuous functions, such as cubic B-splines [2][3] , where   are the unknown coefficients to be determined,   are cubic B-splines with equally spaced knots in the range [  ,   ], and N is now the number of cubic Bsplines in the linear combination in Eq.(S-11).Note that each cubic B-spline is a cubic polynomial in four subsequent subintervals between the knots and it is set equal to zero outside this region.This function is non-negative, and its first two derivatives are also continuous 3 .Then, Eq. (S-2a) is written as follows: , = ∑    (  ;   ; )   =1 (S-12a) where the transformed Kernels using cubic B-splines are written as: With Nx being the number of points used to evaluate the integral in Eq.(S-12b), in full accord with the methodology proposed in Ref. 3 .Thereafter, Eqs.(S-12a) and (S-12b) replace Eq. (S-2a) in the optimization problem used to extract PSD, , and  from the adsorption-desorption isotherms.The advantage of employing the above procedure is that one can use a very fine resolution in terms of the number of pore sizes, Nx, ensuring a high accuracy in the solution of the various integrals in Eqs.(Error!Reference source not found.)and (Error!Reference source not found.),while solving the problems for a small number of B-spline nodes [2][3] , N, so that N<<M (but Nx>>M).

Production of Sorption Kernels of N2 at 77.4 K on Cylindrical Pores of Silica
We follow the method developed by Bonnet and Wolf 4 , and later by Morishige [5][6] , who have modified the Broekoff and De Boer (BDB) method [7][8] , to include solid-fluid interaction at the atomic level, and thermal activation-evaporation to account for activated condensation and cavitation.According to these studies [4][5][6] , and along the framework of a continuum thermodynamic approach, the grand potential of the cylindrical bubble of radius R, per unit pore length, with respect to the filled pore is: Where Δρ is the difference between bulk liquid and saturation density, Δμ is the difference in the chemical potential between the external vapor in contact with a solid and saturated vapor, γ is the liquid-vapor surface tension (assumed constant and equal to its bulk value), and Uc is the solid-fluid interaction potential for a cylindrical pore.
Substituting Δμ, by its relation,  =    ln ( ) , where kB is the Boltzmann constant, and f/f0 is the relative fugacity, which is practically equal to relative pressure, p/p0, for N2 at T=77.4 K, we get: Then we can determine the various characteristic sizes with respect to relative pressure for spontaneous condensation (nucleation), growth (equilibrium), and thermally activated condensation and evaporation (cavitation), applying the various conditions for () according to Bonnet-Wolf (BW) theory 4 .To apply thermally activated condensation and evaporation mechanisms in BW theory we need to impose an energy barrier, Ec.Following Morishige 5 , we set   = 60  , since this value gives the best fit of critical pore sizes with condensation and cavitation relative pressures for N2 on KIT-5, an ordered mesoporous silica for which there are molecular simulation results 9 .
Details on the implementation of the BW theory can be found elsewhere 4 .
From (S13b) we can also determine the statistical film thickness t, by setting   = 0 , and finding the root, Rm , of this equation 6 : Then t is determined from the definition: From eqs. (S14a), (S14b) it is evident that the statistical thickness, t, depends on relative pressure p/p0 and pore radius, R0, in accord with BDB theory.Consequently, t is employed to determine the various sorption kernels, in terms of the amount adsorbed,   , as follows [10][11] : where is the respective transition pressure,   is the pore volume of the material, and   is the adsorbate (liquid) density.
The basic difference between the work of Bonnet and Wolf with that of Morishige is the solid-fluid (SiO2-N2) interaction potential employed.Bonnet and Wolf have considered the attractive part of the cylindrical Lenard Jones 9-3 (CLJ 9-3) potential, in accord with Saam and Cole's theory 12 .According to this potential, the fluid interacts with a continuum solid of density ρs, as follows 13 : ,9−3 (,  0 ) = 2    3   [ 6 (,  0 ,   ) −  3 (,  0 ,   )] with In equations (S16a, S16b), r, is the fluid molecule's distance from the center of the cylindrical pore, R0, is the pore radius, Γ, is the Gamma function, F, is the hypergeometric function, while σsf, and εsf, are the Lenard Jones length and energy parameters for the interaction between the fluid molecules and the solid.Note that n in eq. ( S16b) must be either an integer or a half-integer greater than 1/2 14 .
This interaction potential considers pore walls in the form a continuum solid of infinite thickness and predicts a shallow minimum located too close to the cylindrical surface.
Evidently, some part of the atomic structure of the solid material must be included to Note that all parameters in eq.(S17), including n in eq.(S17b), are defined as previously, except from R0, which is defined as the distance between the centers of oxygen atoms in the first layer of the pore wall and the pore center 6,13 .This radius is

Figure S1 .
Figure S1.Schematic logical block diagram of the iterative procedure to extract the optimal f, z, φ, from adsorption-desorption isotherms.

Figure S2 .
Figure S2.Logical block diagram of the iterative procedure to obtain optimal PSD, φ, for a pair of finite size parameter fk, and pore connectivity, zl, values.