Free energy landscape within the hysteresis regime for fluids confined in disordered mesoporous solids

Thermodynamics and density distributions of fluids and solids confined in structurally disordered mesoporous materials become extremely complex.^1–4 Phase transitions in them typically exhibit hysteretic behavior, signaling the occurrence of out-of-equilibrium phenomena.⁵ One of the intriguing questions emerging in this respect is the extent to which the experimental observations made in these systems can be described within the framework of equilibrium thermodynamics. This and related questions were the focus of a series of studies^6–8 in which mean field theory (MFT) of the lattice gas was used to explore fluid behavior under the action of a random field. The latter was created by including into the lattice a matrix of randomly distributed particles exerting a solid–fluid interaction potential with the fluid sites. With this statistical mechanical model, capturing the most essential properties of disordered pore networks, several important conclusions have been drawn. First of all, the structural disorder has been shown to result in a very complex grand free energy landscape with a multitude of the local metastable thermodynamic states in which the system becomes captured. It was demonstrated that the emerging collective effects, completely absent in collections of independent pores [referred to here as the independent pore model (IPM)] and brought about by the pore interconnectivity, may lead to out-of-equilibrium phase transformations. The novel insight into the fundamentals of these processes gained had direct implications for the characterization sciences. In particular, it was argued that because the collective phenomena are associated with simultaneous transitions (avalanches) over many pores with different pore sizes, rather than in single pores with identical pore sizes, the structure determination in the frame of classical approaches might be meaningless.

Even though the model provided novel qualitative insight, several practical questions remained to be explored. Among them, an intriguing problem is the relationship between the variation of the chemical potential and the spectrum of the grand free energies the system passes through along with this variation. In other words, in how many different states can the system be prepared experimentally at a given chemical potential, and what is the spectrum of the free energies associated with these states? What are the rules determining a particular trajectory among the large spectrum of possible states? A distinct question concerns the equilibrium transition line, i.e., the one connecting the loci of the lowest free energy minima for different chemical potentials. Its existence has been addressed in earlier work,^6,7 but how it is related to the structure of the pore network remains largely unknown. It was argued convincingly that these global grand free energy minima are unattainable in experiments. Nevertheless, it is certainly of interest to explore just how closely they might be approached.

Another important question is related to the internal dynamics of the system following a change in the chemical potential of the system as controlled by the thermodynamic bath surrounding the system. It has been shown experimentally that, in structurally disordered pore systems, there is a relaxation process toward the global grand free energy minimum due to activated hopping between free energy minima.⁹ The relaxation rates are found to slow down with time resulting in a significant time separation between the relaxation dynamics and dynamics in response to the driving force.⁸ The hopping process between the local energy minima is accompanied by the changes in the phase composition. On considering, for example, gas sorption, the system equilibration is accompanied by mass transfer to (release) or from (uptake) the external bath.^10,11 An important question in this regard is: if a certain state within the hysteresis region is attained via some experiment, will the system relaxation be accompanied by uptake or release? In other words, in which direction will the phase composition change due to internal dynamics?

To investigate the questions outlined above, computer simulations of fluids and solids confined in model pore networks can be used, but they are typically extremely demanding computationally.^12,13 Recently, a rigorous statistical mechanical theory for phase equilibria for the simplest pore network model, i.e., for one-dimensional chains of pores with the statistical disorder, in what follows referred to as statistically disordered chain model (SDCM), was developed.¹⁴ The model considers chain-like arrangements of individual pores. In these pores, the occurrences of the first-order phase transitions depend on thermodynamic boundary conditions as determined by the phase states in the two adjacent pores. In this way, due to the inter-pore connections, phase transformation in one pore may eliminate the nucleation barrier in the adjacent pore and, thus, change the transition mechanism there from nucleation-controlled to the one determined by the phase growth from the adjacent pore. Because thermodynamic conditions for nucleation and invasion are different, including also their pore-size dependencies, this gives rise to a very complex interplay between the two transition mechanisms in the entire chain. Consequently, this breaks down unambiguity between the phase state and pore size. Notably, the unequivocal connection between the phase state and pore size is the most essential property of IPM, thus, SDCM turns out to be crucially different from IPM and to mimic more accurately properties of real materials.

Before giving some details about SDCM, it is worth noting that it captures qualitatively the most essential features relevant for phase equilibria in real disordered materials. It should not be confused that the chain-like pore configuration of SDCM represents one-dimensional systems, while real materials contain three-dimensional networks. What really distinguishes SDCM from them is only the number of the inter-pore connections, two in SDCM, and two to several in real materials (with the exception of IPM-like systems in which each pore is connected to all other pores in the sample via the external bulk phase). Because, as it has been noted earlier, the pore interconnections are responsible for the alterations of the boundary conditions for the adjacent pores, different interconnection numbers will result in the different weighting of the nucleation and phase growth mechanisms. Hence, the effects of different interconnection numbers are only quantitative, resulting in different percolation properties. This is in contrast to genuine dimensionality effects in phase transitions, like in the Ising model, where no phase transitions possible for spin chains at finite temperatures.^15,16

Because the most essential physical properties of fluids and solids in pore networks are captured by SDCM, the results obtained for SDCM allowed to explain many features observed in the experiments self-consistently and without need for making some additional assumptions.⁴ Here, we briefly list some of these findings demonstrating the validity and potentials of SDCM. Perhaps one of the most important feature of SDCM is the fact that it correctly reproduces the scanning behavior (chemical potential reversal upon incomplete transformation),⁵ while all IPM-based models fail completely. The number of inter-pore connections in SDCM is fixed at two, and the chains are assumed to be statistically disordered. Hence, the details of phase equilibria are determined by the pore size distribution (PSD), largely deciding on the phase growth probabilities, as well as the chain length, strongly affecting the nucleation probabilities for short chains. The chain-length dependency is a unique attribute of SDCM, explaining the experimentally observed grain-size dependency of the phase transitions in granular porous materials.^4,17 Importantly, in the limit of the chain length equal to one, the SDCM results accurately reproduce that obtained within IPM. SDCM reveals that PSD plays an essential role in distinguishing between weak and strong disorder cases in phase transitions. In this way, SDCM helps rationalizing why phase transition behaviors in two random porous glasses, CPGs, and Vycor porous glass, differ substantially.^18,19 Irrespective of similar pore space topology in them, the results obtained within SDCM reveal that these two materials belong to different disorder classes. This results in avalanches dominating phase transitions in CPGs so that each grain in the sample behaves as a single pore rendering the behavior IPM-like. In contrast, in the strong disorder case for Vycor porous glass, phase coexistence in the monolith is observed over broad pressure and temperature ranges.

An important advantage of the model is that solutions describing different details of phase equilibria, such as distributions of the coexisting phases along the pore space and attained along any arbitrary variation of the chemical potential, can be obtained analytically. With these distributions known, any physical quantity of interest can be computed. In particular, the free energies associated with the thermodynamic states can be calculated within some thermodynamic models. In this work, we address the problems stated above by studying the free energy landscapes in structurally disordered pore chains for solid–liquid and gas–liquid equilibria. By using analytically solvable capillary approximation ( Appendix A) and its analog for freezing and melting ( Appendix B), we obtain important conclusions about the variations of the free energy along both boundary and scanning transitions. We show that the analysis of the results obtained allow us to establish more general properties, such as the location of the global free energy minima, the validity of which extends beyond the simple chain model and applies for more complex pore networks as well. To demonstrate that our results are not just restricted to the capillary approximation, we show selected results obtained with the use of mean-field theory for lattice gas. These results are in complete agreement with the analytical results obtained within the capillary approximation.

The essence of the SDCM, also known in the literature as the corrugated or serially connected pore model, is shown in Fig. 1.^4,20–28 It is composed of cylindrical pore sections with different pore diameters x and identical lengths l serially connected to each other. The pore diameters x_i, where i is the index of a section in the chain with the total number L of the pore sections, are distributed statistically and are drawn from a volumetric pore size distribution (PSD) function φ(x). Thus, the model may be regarded as consisting of a large number of the individual chains, each of the length L, and with random statistics in each chain. In this way, the results obtained can be considered as averaged over many realizations of structural disorder.

The theory we use is based on the single segment (pore section) thermodynamics, as fully described by the set of respective phase transition kernels, and by incorporating the effects of the thermodynamic states of neighboring segments. For two phases A (liquid) and B (gas or solid) in the SDCM, two different phase change scenarios in the individual pore sections need to be considered.²⁶ First of all, in a pore section containing the phase A and having at both pore openings contact to the pore sections filled with the same phase A, the new phase B can only be nucleated. The relevant examples are cavitation^29–31 or homogeneous nucleation of ice^32,33 or liquid bridging on the contrary.³⁴ We refer to all of these mechanisms as nucleation. When nucleation has occurred, the newly nucleated phase fills the entire pore section. Note that the transition points³⁴ for nucleation, either in temperature or in pressure, depend on the pore size x as encoded in the respective transition kernels. Second, if a pore section contains the phase A and has a contact at the pore opening to the phase B, then the latter can trigger phase transformation by the invasion of the phase B. Some literature examples are gas or ice invasion and advanced adsorption or advanced melting.^17,35 This mechanism is, in what follows, referred to as phase growth. The temperature and pressure at which it takes place also depend on the pore size x as determined by the equilibrium transition kernel.

This short description of the phase change scenarios indicates that the phase configuration in a chain (or in collections of chains) cannot be predicted without knowing an exact sequence of the pore sections with different pore sizes.³⁶ This problem can, however, be circumvented by considering large ensembles of random chains, such as modeled by SDCM. In this way, one operates with the probabilities of having certain sequences and computes statistical averages over large ensembles of the pore chains with random disorder statistics.³⁷ Exactly this approach has been used in Ref. 14, and the problem was solved rigorously by using the mean-field-like statistical methods based on the combinatorial analysis.^38–40 The respective solutions can be found elsewhere,^14,26,41 and here we only summarize the most essential points.

In order to understand all steps involved in the computations of the phase states, let us consider the boundary adsorption transition as a representative example. At relatively low pressures, all pore sections in the chains contain only the adsorbed liquid film coexisting with the gaseous phase in the pore core part. This is the initial state for establishing the phase equilibria along the adsorption branch. Let us now consider the system at an arbitrary pressure P and assume that all chains have a length L, which is the number of the pore sections in a chain. It can be stated definitively that, if P is higher than the critical pressure P_c,n(x) to nucleate the liquid bridge in a pore section with the pore size x, then this pore section will be filled with the capillary condensed phase following the nucleation event. Thus, the number N of the pore sections in a chain containing the capillary condensed phase due to nucleation can be found by analyzing PSDs, in line with the original work by Mason.³⁸ Because of statistical disorder, these N sections are randomly distributed along the chain. Each such pore section changes the transition mechanism in the two adjacent pore sections, which can be filled by the capillary-condensed phase if P is higher than the critical pressure for phase growth (invasion) P_c,g(x). Indeed, because P_c,g(x) is lower than P_c,n(x), in addition to those pore sections already filled due to nucleation, some pore sections will be filled due to phase growth (advanced adsorption). The probabilities for such events are determined, once again, by PSD deciding on the probability that the pore size is sufficiently large so that P > P_c,g(x). In this way, each pore section filled with the capillary-condensed liquid due to nucleation will give rise to the formation of a continuous domain composed of pore sections filled with the capillary-condensed phase with the length λ. In our original work, Ref. 14, the accurate derivations for obtaining the distributions of the domain length f(λ), and the parts of the entire pore system filled with the capillary-condensed phase are outlined in details. The desorption branch behavior is treated similarly with the only exception that now the initial state is that where all pore sections are filled with the capillary-condensed phase. The distinct difference to the adsorption case is, however, the fact that two pore sections at the chain ends have contact to the external gas phase. Somewhat similarly, the scanning behaviors are analyzed by explicitly taking account of the initial phase configuration before the scanning curve is computed.

Thus, the phase coexistence properties, either between gaseous and liquid domains or between solid and liquid domains, are quantified in terms of (i) the relative volumetric PSDs associated with the pore sections containing either the phases A or B and (ii) the length distributions f(λ) of the continuous domains with a single phase spanning over many pore sections. The two former distributions are denoted by ψ(x) (by default describing the pores containing the liquid phase) and ψ′(x) (by default describing the pores containing the gaseous/solid phase). By definition ψ(x) + ψ′(x) = φ(x), where φ(x) is the total PSD, i.e., our analysis within SDCM yields which part of PSD contains capillary-condensed phase and in which part the gaseous or solid phases will be found. Notably, ψ(x) and ψ′(x) resulting for the SDCM differ fundamentally from the similar distribution functions obtained for the IPM. Indeed, in the IPM, there is unambiguous correspondence between the pore size x and the phase state in this pore. That means that PSD is simply subdivided into two parts below and above some critical pore size. In contrast, in SDCM, two different pore sections having identical pore sizes may be found in two different thermodynamic states, i.e., containing either phase A or B. The latter is the result of cooperative effects completely absent in IPM.

With known distributions functions, ψ(x), ψ′(x), and f(λ), resulting at different temperatures and pressures for different histories of how a current thermodynamic state was attained, any quantity of interest can be computed. In particular, the average phase composition, θ, describing the (relative) fraction of the liquid phase in the pore system, typically measured in the conventional gas sorption or thermoporometry experiments, is obtained as

where ζ indicates the external control parameter (relative pressure p/p₀ for gas sorption and relative temperature T/T₀ for freezing/melting experiments). In Eq. (1), the kernels K_liq(x, ζ) and K_g,s(x, ζ) describe the relative amounts of liquid in the pores containing liquid or gaseous/solid phases, respectively. In the latter case, K_g,s(x, ζ) reflects the physical adsorption of gas molecules on the pore walls or the occurrence of the non-frozen layers between solid core and pore walls. The kernels used for the calculations in the present work are as follows: (i) for melting and freezing, the kernels for water as described in Ref. 41 were used, (ii) for capillary condensation and evaporation, the kernels were constructed by combining the Kelvin and the Halsey-Hill equations with the thermodynamic parameters for nitrogen. Note that the choice of kernels has only quantitative, but not qualitative impact.

The main focus of the present work is to explore the evolutions of the system free energies with the variations of the chemical potential. Let us consider first the gas–liquid equilibria. As mentioned earlier, for establishing the qualitative behavior, we will use the capillary approximation. It can be shown readily (see for the derivation details  Appendix A) that, in the capillary approximation, the grand potential Ω of the fluid in a cylindrical pore with diameter x in contact with the gas reservoir containing gas at a relative pressure ζ = p/p₀ is equal to

where A is the interfacial area between coexisting gaseous phase in the pore core and liquid layer physisorbed on the pore walls, γ_gl is the surface free energy of this interface per unit area, R_g is the universal gas constant, V_gas is the volume occupied by the gaseous phase in the pore, and υ is the molar volume of the condensed phase. For the sake of simplicity, let us consider only the excess grand potential, ΔΩ(x, ζ) = Ω(x, ζ) − Ω_liq(ζ), over the one found if the pore was completely filled with the capillary-condensed phase at this particular pressure. ΔΩ(x, ζ) per unit pore volume in a cylindrical pore can be represented as

where t is the thickness of the adsorbed film on the pore walls. In Eq. (3), the contribution of the two menisci has been neglected. The latter is justified if the channel length l exceeds significantly the channel diameter x. Equation (3) can further be rewritten as

where the new variables y = x − 2t, γ_a = γ_gl, and E_a = R_gT/υ are used.

Figure 2 exemplifies the transition behavior and ΔΩ(x, ζ) as resulting from Eq. (4) in long cylindrical pores with x = 12 nm. It is instructive to outline the computation procedure for ΔΩ(x, ζ) in Fig. 2(b) in order to understand how these computations were performed in what follows for the disordered chains. The computations of ΔΩ(x, ζ) were based on the phase states known for each point along the isotherms in Fig. 2(a). Indeed, as it has been mentioned in the Sec. II, for each point along the isotherms the two functions ψ(x) and ψ′(x) are exactly known. Notably, because ΔΩ(x, ζ) in Eq. (4) depends only on the surface area and volume of the gaseous domains, only ψ′(x) is of relevance. Thus, with ψ′(x) known, the entire surface area and volume of the gaseous domains can be computed and ΔΩ(x, ζ) evaluated. Because the behavior of ψ′(x) is trivial for one ideal pore, anything can easily be traced to get a flavor of the computation procedure.

The free energies shown in Fig. 2(b) reveal that ΔΩ(x, ζ) first remains equal to zero with decreasing pressure. This corresponds with the state where the capillary-condensed phase fills the entire pore. Obviously, ψ′(x) = 0 and there is no change in the grand free energy due to the contribution of the gaseous phase. As a certain pressure is reached, evaporation from the entire channel occurs. Because we are considering one single pore, the surface area and volume of the gaseous phase in the core part of the pore are easily found by properly taking account of the adsorbed film thickness obtained from the adsorption kernel. For more complex phase configurations, the surface areas and volumes need to be obtained by properly averaging over ψ′(x). According to the data in Fig. 2, evaporation from the channel is accompanied by no change in the grand free energy. This signals about phase transition occurring at equilibrium between two phases and is in agreement with the well-known behavior for this configuration for an open-ended channel.⁴² The respective relative transition pressure will be referred to as ζ_eq.

Further pressure decreases decrease lowers ΔΩ(x, ζ). In this pressure range, the second term on RHS of Eq. (3) dominates over the first one. If now pressure is reversed, ΔΩ(x, ζ) changes reversibly with increasing pressure until ζ_eq is reached. With further increase of pressure, ΔΩ(x, ζ) continues to grow. Notably, the computations of ΔΩ(x, ζ) here are based on the behavior of the adsorption isotherm in Fig. 2(a) revealing the occurrence of the gaseous phase. Altogether, this indicates that the system is found in a metastable state. The condensation transition occurs at ζ_n by nucleating a liquid bridge triggering the filling of the entire channel with the capillary-condensed phase. Once again, this renders ψ′(x) = 0, hence ΔΩ(x, ζ) = 0.

If there is a distribution of the channel sizes, but the channels remain disconnected (the IPM case), the qualitative behavior remains unchanged as demonstrated in Fig. 3. In addition to the boundary transitions, the scanning desorption behavior is also shown. The most important observation is that the free energies along the desorption boundary curve represent the lowest free energy states in the entire hysteresis region.

In order to find the free energies in connected pores, the results obtained earlier for SDCM will be used.¹⁴ In particular, the excess energies ΔΩ(ζ) can easily be found with the help of the distribution functions ψ′(x). The latter includes only the part of the entire PSD in which the adsorbed film on the pore walls coexists with the gaseous phase in the pore core. Indeed, only the pore section containing the gaseous phase contributes to ΔΩ(ζ); hence,

Let us analyze the results obtained with Eq. (5) separately for two cases of weak and strong disorder, in which the phase transition behaviors differ significantly. In particular, for the weak disorder case, the phase transition in the pore chains is dominated by avalanches spanning over almost the entire chains. In contrast, in the strong disorder case, the alternating domains with different phases form along the chains and exist over the entire pressure range in the hysteresis region. Classification of different disorder from the perspective of phase transitions in random pore systems may be based on the disorder parameter ξ as introduced in Ref. 4,

It takes account of both thermodynamic [the first term in the brackets on the RHS of Eq. (6), where x_n and x_pg are the critical pore diameters for nucleation and phase growth, respectively] and structural [the second terms in the brackets on the RHS of Eq. (6), where $x̄$ is the average pore diameter and σ is the distribution width assuming normal PSD] properties of the fluid-porous solid systems to predict the qualitative phase transition behavior. Thus, ξ > 1 corresponds with the weak disorder case, while ξ < 1 signals about the strong disorder. As an example, for ξ = 1 nucleation of the capillary-condensed phase in only 2% of the pore sections leads to the filling of 98% of all pore sections due to phase growth.

The excess free energies obtained for the weak disorder case are illustrated in Fig. 4 representing a system with ξ = 3.8. They reveal that, in contrast to the IPM, the system encounters the metastable states now along both boundary isotherms. If metastability on adsorption is caused by the delayed nucleation of the capillary-condensed phase in line with the findings made for IPM, metastability on desorption is caused by the pore blocking effect. ΔΩ(ζ) calculated along two boundary isotherms coincide at a certain relative pressure ζ_x. Remarkably, this crossing point is common also for all desorption scanning isotherms. The fact that the density distributions and phase compositions resulting at ζ_x along different isotherms are different suggests that, at this pressure, two phases are found in thermodynamic equilibrium with each other. The origin of the latter observation will be explained in Sec. IV.

In the strong disorder case, as exemplified in Fig. 5 for a system with ξ = 0.7, there are several points distinct from that observed for the weak disorder case. First of all, metastability along the adsorption branch is now less pronounced than on desorption. Second, the crossing point, as characterized by the respective pressure ζ_x is shifted notably toward the pressures associated with the boundary adsorption transition. Moreover, this crossing point no longer coincides with the crossing point for the ΔΩ(ζ) curves associated with the desorption scanning isotherms.

In order to show that the phenomena reported are not simply attributes of the capillary model, we present the computations of the grand free energies by using the mean field theory (MFT) for the lattice gas.⁴³ The details of the respective theory and its implementation for disordered channels can be found elsewhere.²⁸ The calculations for an ideal channel and a collection of independent pores are in qualitative agreement with the results obtained for the capillary model and are not shown. Here, we present only the results obtained for the disordered chains composed of slit-like pore sections. They are demonstrated in Figs. 6 and 7 for the weak and strong disorder case, respectively. In full agreement with the results of the capillary model, in the former case, one finds a common crossing point for the grand potential, while in the latter case it is absent. On considering all other aspects, the qualitative behavior is found to be identical. This is evidence that the capillary model captures accurately the essential physics of the phase coexistence and, due to its simplicity and analytical character, can be used to explore the system behavior further without loss of generality and without need for demanding computations.

The phenomena discussed will be general for any transitions. To illustrate this, we also consider solid–liquid coexistence in confined spaces. The grand potential Ω of the substances confined in a single cylindrical channel of size x at temperature ζ = T/T₀ can be expressed as (see for the derivation  Appendix B),

where Ω_liq(x, ζ) is the grand potential of the supercooled liquid completely filling the pore space, A and V are the surface area and volume of the ice core in the pore, respectively, γ_sl is the surface energy per unit area of the ice–liquid interface, T₀ is the bulk transition temperature, and ΔH_f is the heat of fusion per unit volume. Because Ω_liq(x, ζ) is generally unknown, in line with the gas–liquid equilibria considered earlier, we evaluate ΔΩ(x, ζ) = Ω(x, T) − Ω_liq(x, ζ), i.e., the excess grand potential over the supercooled liquid phase per unit volume. In a spirit of Eq. (4), ΔΩ(x, ζ) can be expressed as

where y = x − 2τ, τ is the thickness of non-frozen layers formed between the ice core and pore walls,^41,44,45γ_a = γ_sl, and E_a = ΔH_f. It turns out that Eqs. (4) and (8) have similar functional forms [note that, for ζ not much less than 1, ln(ζ) ≈ ζ − 1]. Thus, two phenomena, gas sorption, and melting/freezing will differ only quantitatively and can be discussed mechanistically from the same perspective.

For the collections of pores with different pore sizes, the excess grand potential per unit volume can be obtained as

ΔΩ(T) obtained for single channels and collections of individual channels with different channel diameters behave similar to the behaviors shown in Figs. 2 and 3. Here, Figs. 8 and 9 demonstrate the grand potentials obtained for the SDCM in the weak and strong disorder cases, respectively. In full agreement with the scenarios observed for the gas–liquid coexistences, for the weak disorder case, there is one temperature T_x at which the free energies are irrespective of the history of the system preparation. In the strong disorder case, however, this common temperature is absent.

Let us discuss the finding presented in Sec. III from the perspective of the statistical theory for SDCM. As revealed by the respective equations and their results shown in Figs. Figs. 4 and 7–9, the solid–liquid and gas–liquid equilibria differ only qualitatively. Hence, in what follows, the melting/adsorption and freezing/desorption phenomena may be discussed interchangeably.

We start the discussion with the puzzling observation made for the weak disorder case of the coinciding excess grand potentials at ζ_x irrespective of how this temperature or pressure was approached. As it has been discussed earlier, in the weak disorder case, all transitions occur in the form of avalanches spanning over entire chains.^4,8 For example, on heating the entire pore chain melts in almost one step as soon as the liquid phase is nucleated in the smallest pore section in the chain. Thus, due to the statistical distribution of the sections with the smallest pore sizes from one chain to another, the phase coexistence established is given by the collections of the chains containing either liquid or ice phases in the entire chains. In this case, ψ′(x, ζ), i.e., the distribution function describing PSD of the pore sections containing the ice or gas phases, obtained along the boundary melting or adsorption curves to a very good approximation differs from the PSD φ(x) only by a numerical factor. The latter factor simply reflects the relative fraction of the chains in the entire ensemble containing the capillary-condensed, liquid phase. This is evident in Fig. 10(a), exemplifying these distribution functions resulting along the boundary melting curve as predicted by the statistical theory. Note that the freezing and desorption processes occur similarly, where the freezing and desorption processes occur in the entire chains in one step as soon as ice or gas can invade the smallest pores responsible for pore blocking. Thus, for the weak disorder, also along any freezing and desorption transitions, including the scanning freezing and scanning desorption transitions, the shape of ψ′(x, ζ) will remain almost identical to φ(x). This can also be judged from Fig. 10(b) showing ψ′(x, ζ) obtained at the relative temperature ζ_x where all free energies obtained via the boundary freezing and scanning freezing transitions coincide.

By noting that, for this specific case of weak disorder, the common relative temperature ζ_x is very close to that given by ΔΩ(ζ_x) = 0, Eq. (9) yields the following equation for ζ_x,

While in the weak disorder case ψ′(x, ζ) = const × φ(x), ζ_x is irrespective of where the freezing scan has been started. Finally, Eq. (10) can be solved for ζ_x to yield

where ⟨y⟩ and ⟨y²⟩ are the mean and mean square crystal sizes obtained upon averaging over the number PSD φ(x)/x², respectively. Interestingly, Eq. (11) coincides formally with the Gibbs–Thomson equation predicting transition at equilibrium, i.e., without free energy change, for a cylindrical ice crystal formed in a cylindrical uniform pore and having the crystal diameter,

ζ_x = 0.95 obtained using Eq. (11) is shown in the inset of Fig. 9(b) and is found to perfectly coincide with the crossing point for ΔΩ(ζ).

Similar conclusions can be drawn for gas sorption yielding

which is the Kelvin equation applied for gas invasion into uniform cylindrical pores of size x_c = y_c + 2t with y_c given by Eq. (12).

With increasing structural disorder, the phase coexistences between the ice and liquid phases or the gas and liquid phases become possible also along the individual chains. Under these conditions, the shape similarity between φ(x) and ψ′(x, ζ) becomes lost. This is proven by, e.g., the data of Fig. 11(a) obtained for the strong disorder case. It is noted that ψ′(x, ζ) obtained along the boundary melting (or adsorption) transition are found to differ significantly from φ(x). Not surprisingly, these differences remain or even amplify along the subsequent scanning freezing (or desorption) transitions as demonstrated in Fig. 11(b). In this case, clearly, no common ζ exists where all free energy coincides due to different distribution functions in Eq. (10).

We have seen earlier that, for the collections of independent pores (see, e.g., Fig. 3), all states within the hysteresis region approached along the melting or adsorption branches are less energetically favorable than those obtained along the freezing or desorption branch. This behavior changes cardinally when disordered chains are considered. Here, the metastable states emerge along both boundary transitions. In the case of weak disorder, see Figs. 4 and 8, ζ_x divides all states into two regions in which either boundary melting/adsorption or boundary freezing/desorption yield the lowest energy states. Consequently, any other experimental realizations of different thermodynamic states at these temperatures or pressures will result as metastable, higher grand free energy states.

The phenomena become more complicated with increasing disorder strength, as shown in Figs. 5 and 9. The grand free energy profiles obtained along the freezing and desorption scanning transitions reveal that both phase configurations obtained along the boundary freezing/desorption and boundary melting/adsorption transitions can exhibit metastabilities. Indeed, it is seen that in the vicinity of the crossing point between ΔΩ(ζ) obtained along the boundary melting/adsorption and freezing/desorption transitions, the states resulting along the scanning transitions are found to yield lower grand free energies. Thus, just slightly more complex experiments, such as scanning experiments, allow for the creation of the thermodynamic states with grand free energies lower than those obtained along the boundary transitions. This observation poses an important question: what are the states with the lowest grand free energies that one may create by performing the freezing and desorption scanning experiments?

To answer this question, we have compared the grand free energies in the hysteresis region as obtained along the boundary transitions and a large family of the scanning experiments for each ζ. For a given ζ, we have compared the grand free energies obtained along different scanning transitions and selected the one yielding the lowest energy. The respective θ yielding the lowest energy attainable via the scanning experiments are shown in Figs. 8(a) and 9(a) by the dashed-dotted lines. It turns out that these lowest grand free energy states are always those obtained along the freezing and desorption scans. The melting and adsorption scanning experiments yield higher grand free energies.

As revealed by the dashed-dotted line shown in the weak disorder case, Fig. 8, ζ_x separates two different thermodynamic states in the pores yielding the lowest energies for ζ above and below ζ_x. Indeed, in a close vicinity of ζ_x there is a sharp, but continuous transition separating the states created along the freezing/desorption scanning experiments and yielding the lowest energy states. The situation is different for the case of strong disorder. Here, there is no sharp separation between the states obtained along the boundary transitions. Rather, the lowest energy states result from the scanning experiments and are associated with the coexisting ice/gas and liquid phases over a broad range of ζ.

The next important question concerns the states, i.e., geometric phase configurations, yielding the absolute minima in the grand free energy at various pressures and temperatures.⁶ It turns out that these states can be predicted analytically based on the structure of Eqs. (2) and (7). They reveals that, for the given ζ and average phase composition, among all possible geometric configurations the lowest grand free energy state will be associated with the one yielding the lowest surface-to-volume ratio A/V for the ice/gas core. The latter is achieved with the phase configuration ψ′(x, ζ) resulting in the IPM limit L = 1 along the equilibrium freezing or desorption transitions. In this case, the volume V of the ice/gas core is contributed by the largest pores having simultaneously the lowest specific surface area. In this way, the lowest surface-to-volume ratio A/V is attained, which, consequently, yields the lowest free energy state. Note, that this conclusion is quite general and is applied to networks with arbitrary interconnectivities rather than only pore chains considered earlier.

The respective IPM-like boundary transitions, corresponding with the global grand free energy states, are as well shown in Figs. 8 and 9 for two disorder cases. They show the average phase composition, as expressed here by the relative fraction of the liquid phase, needed at given ζ to realize the state with the lowest energy. The same phase compositions obtained, e.g., via scanning experiments will yield, thus, states with the higher grand free energies. Showing the three-dimensional plots by complementing the phase composition—temperature plots with the third dimension associated with the free energies are outside of the scope of this paper and will be presented in a subsequent publication.

The global grand free energy line divides the hysteresis region into two regions. All states prepared along any arbitrary variation of the externally controlled ζ and having the final average phase composition in the region above this line may minimize their energies by a process in which the fraction of the liquid phase is lowered. For this, either freezing or desorption is needed. Vice versa, the states with the resulting liquid fractions below the global free energy line will necessitate melting or adsorption as the process to lower the free energy. This finding allows us to make a global prediction for the intrinsic relaxation processes occurring at constant chemical potential, either temperature or pressure.^9,11 Their direction will be such that the average phase composition will change toward the one corresponding with the global free energy minimum. As just argued, the latter can easily be established by computing the boundary freezing or desorption transitions as soon as the real PSD of the material under study is known. Luckily, the most accurate PSD can be obtained by an analysis of the melting or adsorption boundary transitions by applying the statistical theory for SDCM. As it has been discussed elsewhere, these transitions are only marginally subject to percolation constraints and predict accurately the phase configuration and composition irrespective of the pore space interconnectivity.^4,14 Hence, irrespective of its one-dimensional character, the statistical theory for SCDM can accurately predict the melting and adsorption transitions in three-dimensional pore networks and, thus, allows for obtaining of the most accurate PSD.²⁸ 

The main objective of the present work was to study the grand free energy landscapes for fluids and solids confined in mesoporous materials with the structural disorder. The latter gives rise to very complex density distributions along the pore spaces depending on how the system chemical potential was varied. If these distributions or geometries of the phase configurations are known, the free energies can be evaluated based on some modeling approaches. Because these distributions are unknown, the problem was approached previously in a few studies with the help of computer simulations of selected models mimicking structural disorder. Even though these studies shed light on non-equilibrium phenomena occurring in these systems, in-depth understanding has not previously been possible.

Here, we used the advantage of the statistically disordered chain model (SDCM), allowing for analytical solutions in terms of the density distributions, to explore the system grand free energy the along boundary and scanning loci. Because SDCM captures all essential physics of pore networks, the results obtained are qualitatively valid for porous solids with arbitrary connectivities of the pore networks. With the analytically available phase configurations, the free energies could be efficiently computed for a broad spectrum of the system states along different pathways of the variation of thermodynamic parameters. In particular, we have performed the respective calculations for gas sorption in the frame of the capillary approximation. By using the mean-field-theory of lattice gas in similar pore systems, we have proven that the results obtained are irrespective of the model used. Furthermore, by considering freezing and melting in disordered pore chains, we have shown their complete similarity to capillary condensation phenomena.

In the first part, we demonstrated how the pore connectivity results in the emergence of the metastable thermodynamic states along both boundary transitions beyond that observed with the collections of independent pores. It turned out that the system behavior for the pore chains depends strongly on the disorder strength. In particular, we consider two situations of (i) weak disorder, where the phase transformations in the pores occur in the form of avalanches spanning the entire pore networks, and (ii) strong disorder, where the two phases coexists in the pore networks.

We find that, in the former case, there exists a specific pressure or temperature where the system grand free energy is unchanged irrespective of how the states were approached. To demonstrate this, we have explored the boundary and a family of scanning transitions showing a common crossing point for the free energies. This common point is found to separate all experimentally realizable thermodynamic states into two regions where either boundary adsorption or boundary desorption yields the lowest grand free energy state. By means of theoretical analysis and by direct evaluation of the phase configurations resulting in the weak disorder case, we interrelate this finding with the specific behavior of the density distribution functions, which remain almost intact when a large collection of the pore chains is considered. From the experimental point of view, this corresponds with non-monolithic porous solids made out of many grains with different disorder realization in each grain. Indeed, in this case of weak disorder considered, the phase transitions in each grain occur in a form of avalanches spanning over the entire pore network in the grain. For statistical reasons, there will be from grain to grain distributions of the smallest pores triggering these avalanches at slightly different temperatures or pressures. However, because PSDs in the different grains differ marginally, the PSD in the grains, in which the phase transition has not yet occurred, will remain identical with the original PSD for the entire sample.

In the second case of strong disorder, the free energy landscape is determined by particular realizations of the phase configurations. The most instructive example concerns monolithic porous solids, where the variation of chemical potential results in the variation of the density distributions in the material. In this case, a common crossing point is no longer found. As an interesting observation, however, it is found that the desorption or freezing scanning experiments allow for experimental realization of the thermodynamic states with lower free energies than those obtained along the boundary transitions. In this way, we find the lowest grand free energy states attainable with the scanning experiments.

One of the main findings of this work is identifying the global grand free energy minimum and providing with a simple method for computing it. By a rigorous analysis, we show that this coincides with the limit of independent pores prepared in the thermodynamic states approached along the desorption or freezing boundary transition. Even though the global free energy minimum state is not attainable experimentally in the interconnected pore systems, finding its location has a paramount importance. For example, it becomes now possible to predict the direction of the global relaxation process for any state realized in the experiment. Indeed, at a given chemical potential, the average phase composition of the experimentally prepared state will relax toward the one corresponding with the global free energy minimum state. As an experimental means to find these global free states, the monolithic material containing the pore network can be ground to very fine powder, so that the single pore limit is approached. According to the analysis presented in this work, the desorption or freezing boundary transitions measured with this powdered sample will correspond with the states of global free energy.

In future studies, it will be interesting to compute the entire spectrum of the free energies that can be realized in the system and to identify among them the family of states accessible in the experiments.

The German Science Foundation (DFG) is gratefully acknowledged for the financial support (Project No. 411771259).

The authors have no conflicts to disclose.

A. Alzaidi: Data curation (equal); Formal analysis (supporting); Investigation (equal); Software (equal); Validation (equal); Visualization (lead). E. S. Kikkinides: Conceptualization (supporting); Data curation (lead); Formal analysis (supporting); Methodology (supporting); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (supporting). D. Schneider: Data curation (supporting); Formal analysis (supporting); Methodology (equal); Software (supporting). P. A. Monson: Conceptualization (supporting); Formal analysis (supporting); Methodology (supporting); Writing – review & editing (supporting). R. Valiullin: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Supervision (lead); Validation (equal); Writing – original draft (lead); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

To evaluate the grand free energy for the gas–liquid coexistence in a pore, an illustration of the experimental system in Fig. 12 can be used. The chemical potential and pressure in the system are fixed by the gas reservoir. Volumetric changes of the porous solids during adsorption/desorption are minor, hence the pore volume V_pore can be considered to be constant. Temperature is fixed along all experiments and only pressure is varied. Because the particle number in the pore may vary, the grand canonical ensemble is used to evaluate the system behavior.

The grand potential for the liquid-filled pore [Fig. 12(a)] is

where p_l is pressure in the capillary-condensed liquid, γ_lw is the surface energy of the liquid–pore wall interface, and A_w is the surface area of the pore walls. If the gaseous phase is formed in the core part of the pore [Fig. 12(b)], the grand potential Ω_e is

where p_core is the pressure in the pore core filled with the gaseous phase, γ_gl is the gas–liquid interfacial energy, A_core is the surface area of the gas phase in the pore. Pressure p_lf in the adsorbed film is reasonable to approximate by that in the capillary-condensed liquid p_l. With V_pore = V_lf + V_core, Ω_e results as

In what follows, we evaluate the difference ΔΩ = Ω_e − Ω of the grand free energies,

Quite generally the differential of the grand potential is

Let us consider the equilibrium condition for evaluating p_l − p_core in Eq. (A4), i.e., Ω_core + Ω_lf = Ω_l. This means that in the vicinity of this state

For keeping the equations compact, the following subscripts will be used in what follows: g for the gaseous phase in the pore core, l for liquid, lf for the adsorbed film. With Eq. (A5), Eq. (A6) for isothermal conditions expands to

With dV_l = dV_g + dV_lf and by assuming reasonably that the pressures in the capillary-condensed liquid and in the adsorbed liquid film are identical, this equation simplifies to

By noting that N_i = n_iV_i, where n is concentration, Eq. (A8) modifies to

The pressure difference is now found as

Using

and by assuming the ideal gas behavior, i.e., $(∂μ/∂pg)T=Vg/Ng$⁠, Eq. (A10) results in

By substituting this into Eq. (A4), we finally obtain

Let us now fix the difference in the grand potential ΔΩ and find the pressure, which results in this ΔΩ. This can be done by integrating with the reference point taken at p₀, i.e., the pressure of saturated vapor,

The integrals are readily evaluated to yield

At equilibrium, i.e., for ΔΩ = 0, we obtain the Kelvin equation describing the decrease of the equilibrium transition pressure by confinement in the pore,

Typical experimental conditions of thermoporometry experiments, in which freezing and melting transitions are measured in porous solids, are illustrated in Fig. 13. The porous solid under study is imbibed with liquid (for example, water) and excess liquid is provided. With an initial freezing/melting run, the system is prepared in a state where the bulk phase is frozen (ice), while the intrapore phase still remains in the liquid state. In this way, the boundary conditions at the pore openings are fixed, and the subsequent freezing or melting experiments are performed with the well-defined initial state. It is important to note that, upon freezing of the intra-pore liquid, a liquid-like layer (NFL) is formed between the ice core and pore walls. Because ice has a lower density than water and the liquid is only slightly compressible, upon freezing some part of liquid is pushed out of the pores and crystallizes somewhere as bulk ice. It is further assumed that, upon melting, the missing water density in the pores is compensated by melting bulk ice.

Because the porous solid is always surrounded by bulk ice, thermodynamic conditions correspond with the fixed chemical potential of bulk ice μ_ice. We assume that the porous solid does not expand upon freezing of the confined liquid, hence the pore volume is fixed. The number of water molecules in the pore is not conserved. Hence, the appropriate ensemble to consider is the μVT-ensemble and the respective thermodynamic potential to evaluate is the grand canonical potential Ω.⁴⁶ 

The grand potential for the liquid-filled pore [Fig. 13(a)] is

where p_l is the pressure in a liquid, V_pore is the pore volume, γ_lw is the surface energy of the liquid–pore wall interface, and A_w is the surface area of the pore walls. If the ice crystal is formed in the core part of the pore [Fig. 13(b)], the grand potential Ω_f is

where p_core is the pressure in the ice crystal, γ_ls is the liquid–ice interfacial energy, and A_core is the surface area of the ice crystal. Pressure in the NFLs is reasonable to approximate by that in liquid p_l. With V_pore = V_NFL + V_core, Ω_f results as

In Eq. (B3), the long-range interaction between the ice core and pore walls was neglected.

In what follows, we evaluate the difference ΔΩ = Ω_f − Ω of the grand free energies,

Quite generally the differential of the grand potential is

Let us consider the equilibrium condition for evaluating p_l − p_core in Eq. (B4), i.e., Ω_core + Ω_NFL = Ω_l. This means that in the vicinity of this state

For keeping the equations compact, the following subscripts will be used in what follows: s for ice core, l for liquid, n for NFLs. With Eq. (B6), Eq. (B6) expands to

With dV_l = dV_s + dV_n and by assuming reasonably that the pressures in the liquid and in NFLs are identical, this equation simplifies to

Let us denote with the lowercase s entropy per particle, so that S_i = N_is_i, and use that N_i = ρ_iV_i/m₀, where ρ is the mass density and m₀ is the mass of one particle. Then, equation Eq. (B8) modifies to

The pressure difference is now found as

dμ/dV_s can be found as

As the next step, we need to find how the chemical potential of the bulk ice varies with temperature. Here, we need to refer to the experimental realization of common thermoporometry experiments in which bulk frozen ice (surrounding the pore system) is found in contact with its vapor phase in a sample tube or chamber (the formation of a non-frozen layer between the gas and ice phases will not change the result). Hence, Ω_s = Ω_g and p_s = p_g along the coexistence line, where the sub-script g refers to the gas phase. An analog of Eq. (B7) in this case is

which can further be evaluated to

The last term on RHS introduces a small correction and can reasonably be neglected. Finally,

Taking account of

where ΔH_f is the latent heat of fusion per unit volume of ice and T₀ is the equilibrium transition temperature, Eq. (B14) becomes

Substituting it into Eq. (B4), we finally obtain

Let us now fix the difference in the grand potential ΔΩ and find the temperature, which results in this ΔΩ. This can be done by integrating with the reference point taken at T₀,

which leads to

At equilibrium, i.e., for ΔΩ = 0, the Gibbs–Thomson equation describing lowering of the equilibrium transition temperature under confinement,

is obtained.