Numerical characterization of the density of metastable states within the hysteresis loop in disordered systems

Joël Puibasset

1

Numerical characterization of the density of metastable states within the hysteresis loop in disordered systems

Joël Puibasset

Centre de Recherche sur la Matière Divisée, CNRS-Université d'Orléans, 1b rue de la Férollerie, 45071 Orléans cedex 02, France E-mail : puibasset@cnrs-orleans.fr Abstract: An improved approach is proposed to analyze the density of metastable states within any hysteresis loop, as those observed in magnetic materials or for adsorption in porous materials. Except for few analytically tractable models, most calculations have to be performed numerically on finite systems. The main points to be addressed thus concern the average over various material samples (the socalled realizations of the disorder), and the finite size analysis to estimate the thermodynamic limit. As an improvement of previously existing methods, it is proposed to introduce the Fourier transform of the density of metastable states (characteristic function). Its logarithm is shown to be additive and can straightforwardly be averaged over disorder. This procedure leads to a new definition of the complexity in finite size, giving the usual quenched complexity in the thermodynamic limit, while being better suited to perform finite size analysis.

Introduction

Hysteresis appears in many situations, the most frequent being magnetism. It is also found for fluids confined in nanopores: for some experimental conditions, the fluid adsorption and desorption isotherms may differ, and exhibit a hysteresis. Puzzling experimental observations

have not yet been accounted for by a molecular approach. For instance, one observes large hysteresis in regular mesopores even when they are closed at one end. [1][2][3] Wallacher and collaborators 4 suggested that the hysteresis emerges as a consequence of intrinsic disorder within the pores. The implication of such a disorder to explain hysteresis was shown in structures like aerogels. [5][6][7] Taking into account the effect of physico-chemical heterogeneities on fluid adsorption, 8 and the interconnections between the elementary domains delimited by these heterogeneities, we arrived at the conclusion that these heterogeneities within pores give rise to a large hysteresis loop with a triangular shape. 9,10 A simplified lattice gas approach gives similar conclusions. 11 The boundary of the hysteresis corresponds to the main adsorption and desorption branches.

As for magnetism, the inner region of the loop can be explored by performing incomplete adsorption/desorption cycles (scanning curves). This shows the existence of a plethora of metastable states and gives some information regarding their organization. The same procedure could be used in our numerical experiment. However, it misses many metastable states, as can be shown by changing the size of the reservoir. 12 We thus perform a systematic exploration of all metastable states within the hysteresis. These metastable states arise as a consequence of the complex free energy landscape which is ubiquitously observed in systems exhibiting disorder. In our case, this disorder is a consequence of the strong heterogeneity in the physico-chemical properties of the inner surface of the porous material.

An interesting issue is to characterize the density of metastable states within the hysteresis loop. Such studies have a long history in magnetic systems, but recent works have shown their relevance to adsorption hysteresis. 13,14 Previous works in both fields have shown that this density grows exponentially with system size inside the loop. One thus introduces the complexity, i.e. the logarithm of the density normalized to the system size. This quantity depends on the spatial distribution of the heterogeneities in the pore, and one has to somehow perform an average over all possible pore realizations. This is however not immediate since the average may be performed either on the density of metastable states, leading to the socalled annealed complexity, or on the complexity itself, leading to the so-called quenched complexity. The annealed complexity, which characterizes the average density of metastable states, is unsuited to describe the hysteresis loop in the thermodynamic limit. 13 On the other hand, the quenched complexity, which characterizes the typical density of metastable states, is expected to lead to the correct thermodynamic limit, and has been extensively used to characterize the adsorption/desorption hysteresis loop. [13][14][15] Another approach may however be proposed, based on the characteristic function associated to the density of metastable states. We show that its logarithm is straightforwardly averaged over disorder, providing a natural characterization of the distribution of metastable states in terms of average over disorder. This approach is shown to be in close connection with the previous one, leading to the usual quenched complexity in the thermodynamic limit. The calculations are more involved and more sensitive to numerical inaccuracies, which may explain why it has never been introduced in this context. This method is however shown to exhibit a remarkably faster convergence in terms of finite size analysis, because it relies on an exact description of the growth of the density of metastable states with system size when subsystems become essentially independent. From this point of view, it proves to be particularly suited if one is interested in numerical extrapolations to the thermodynamic limit.

These results have been obtained in the context of fluid adsorption in inhomogeneous linear pores to be described hereafter, but are expected to be quite general.

The paper is organized as follows: After a presentation of the molecular model and the standard methods, we introduce the approach exploiting the properties of the characteristic function associated to the distribution of metastable states. We finally show that this method is able to provide quantitative information regarding the quenched complexity in the thermodynamic limit. The last section summarizes and gives some conclusions.

Model and methods

Our model is designed to catch the main features of simple fluid adsorption in heterogeneous siliceous linear nanopores (3 nm in diameter). It can be seen as a simplified model for MCM-41 or oxidized porous silicon, if one admits that these materials exhibit corrugated pores. This is however a highly debated issue and alternative models have been proposed. [16][17][18] The results are however expected to be transferable to any other system exhibiting heterogeneous nanopores, like for instance the nanoporous alumina, 2,3 or pores modulated at mesoscopic scale. 19 All interactions are supposed to follow the Lennard-Jones (12, 6) potential. The parameters are those for argon-argon and argon-oxygen species to describe fluid-fluid and fluid-silica interactions. 9 Heterogeneities arise as a modulation of the fluid/wall interaction along the pore axis: this results in a string of more or less attractive domains (see Fig. 1). In our model, we chose domains to be 7 nm long. Due to short range interactions (less than the domain length) and rapid decrease of correlations, the fluid properties in a given domain are essentially determined by the imposed chemical potential  and the fluid state (density) in the two neighboring domains (the first and last domains are in contact with the bulk gas state since pores are supposed to be open at both ends). For a given , the number of such fluid states (densities) is actually quite small, which allows considering all of them to calculate all "conditional" isotherms of the domain under consideration. The fluid confined in a domain may either be in the gaslike state, denoted G, corresponding to the fluid adsorbed only at the walls, or in the liquidlike state, denoted L, corresponding to the fluid filling the whole domain. This gives at most four conditional isotherms for each domain. After these calculations are performed for all domains, the fluid properties can be determined in the whole pore for any chemical potential path. This approach thus allows a decoupling between the atomistic level calculations performed by molecular simulation in each domain (Grand Canonical Monte Carlo), and the large-scale calculation of the fluid state in the heterogeneous pore of micrometric length. It is emphasized that this is not a simple independent domain approach since the strong interaction between adjacent domains are taken into account. More details on the model can be found in refs [ 9 ] and [ 10 ]. For instance the main adsorption isotherm, from vacuum to saturation, is calculated as follows: at low , all domains are in the gaslike state and follow the gaslike branch of its G-G conditional isotherm until one domain reaches its limit of stability and jumps to the liquidlike state. This obviously impacts the fluid state in the neighboring domains with a possible avalanche, until a new self-consistent state is reached. This avalanche corresponds to the propagation of a liquid meniscus. The chemical potential is further increased, and the same procedure is applied until saturation is reached.

Desorption and scanning isotherms are obtained in a similar way.

As explained in the Introduction, we need to determine all the possible metastable states for the confined fluid. Let us consider a pore with N domains. At chemical potential , the fluid may be in one of two states (G or L) for each domain. One then expects at most 2 N possible fluid states, which allows a systematic analysis for small enough N. Most of these configurations are actually impossible due to the interdependence between neighboring domains. Our algorithm selects the few ones that are consistent with the conditional isotherms. These configurations are stable against thermal fluctuations (taken into account in the pre-calculated isotherms) and correspond to the metastable states (see Fig. 2). It was shown in previous work that the number of metastable states grows exponentially with system size but remains much smaller than 2 N . The density

  N ; ;   N of metastable states in the -  plane is defined so that      d N ; ;
N is, for a given chemical potential , the number of metastable states in the interval d around , where  denotes the amount adsorbed normalized to its maximum. The number of metastable states at a given  is

          d N N ; ; ; N N .
These quantities actually strongly depend on the particular distribution of heterogeneities in the porous material. In our case, it corresponds to a given string of domains characterizing the linear pore. Any such "realization of the disorder" will be denoted by . As usual in disordered systems, it is interesting to consider both the average and typical number of

metastable states:   N ; ;    N and  ] ; ; ln exp[ N    N
where the bar denotes the average over the disorder . Since these quantities are both expected to scale exponentially with N within the hysteresis, one introduces the corresponding so-called annealed and quenched complexities:

    N N N ; ; ln 1 lim ; A      N      (1)     N N N ; ; ln 1 lim ; Q      N      (2)
The calculation of the finite-size annealed complexity

    N N N ; ; ln 1 ; ; A      N    (3)
is straightforward. However, the finite-size complexity

  N N ; ; ln 1    N  diverges each
time the density of metastable states reaches zero, which may always occur for some configurations of the heterogeneities in the pore, leading to an undefined quenched complexity. To circumvent the problem, Pérez-Reche 14 and collaborators proposed a method inspired by an analytical calculation, 13 which consists in introducing the Legendre-Fenchel transform

  N g ; ;    of the finite size disorder dependent complexity   N N ; ; ln 1    N  :                     g N N N g ; ; ln 1 max ; ; N (4)
This quantity is always defined for any real value g, even if

  N ; ;    N reaches zero. Let us denote   N g ; ;   
the -value which maximizes the quantity between parentheses in Eq.

4. One has:

            g N N N g    ; ; ln 1 ; ; N (5)
Averaging over disorder, one gets:

      N g g N N N g ; ; ; ; ln 1 ; ;           N (6) where     N g N g ; ; ; ;      
. The quenched complexity can thus be defined for a finite system as:

    ) ; ; ( ; ); ; ; ( ; ; F - L Q N g N N g N             (7) where ) ; ; ( N g   is implicitly defined by   N g ; ;   .   N g ; ;    is expected
        ; ; ; lim Q F - L Q       N N (8)

The characteristic function

As any other quantity, the annealed and quenched complexities are expected to depend on the system size N. A general form was proposed:

1 1 ln       bN N aN
involving three free parameters which depend on  and . It has proven to be a fairly good estimate of the system size dependence, which can be used in a standard finite-size analysis for instance. 14   is the expected thermodynamic limit. Figure 3 shows the annealed complexity as a function of  for  = -8.55, where the hysteresis is the largest, and for various system sizes N = 15, 20, 25, and 30. Noticeable finite size effects are visible. The evolution of the maximum can be described using the above mentioned ansatz with the parameters:   = 0.25, a = 0.65, and b = -0.64.

Similar results are obtained for the quenched complexity. Can this system size dependence be interpreted in terms of interdependence between domains? The answer is no, since a straightforward calculation shows that the annealed and quenched complexities are expected to exhibit finite size effects even if the domains are independent. A question thus arises: what contribution to system size dependence is specifically due to the interconnections between domains? To illustrate this point, we have considered pores of size N = 30 made of the union of two independent smaller pores of size N = 15. The calculated annealed complexity is given in the inset of Fig. 3. As can be seen the results are close to the fully interconnected pores with N = 30. This suggests that the main contribution to the observed system size dependence is not due to interactions between subsystems. The aim of this section is to propose an analysis which exhibits solely the finite size effects due to domain interdependence.

The main idea underlying the average procedure over disorder is that a large system can be seen as the ensemble of many subsystems which tend to be essentially independent in the thermodynamic limit. As a consequence any additive quantity is straightforwardly averaged. This is for instance the case for the envelope of the hysteresis (main adsorption/desorption loop). The case of the distribution of metastable states within the hysteresis loop is more

involved. If we denote       N ; ;
N the distribution of metastable states for the subsystem  made of N domains, then the distribution for the union of two independent subsystems  and  is given by:

                                                    d d N N N N N N N N ; ; ; ; ; ; , N N N (9)
where the integral is performed over all possible values of adsorbed densities in  and  subsystems, and  is the Dirac measure. This formula can easily be generalized to more that two subsystems. Introducing the Fourier transform (characteristic function) of the distribution:

                 d iN N N e ; ; ; ; ~N N (10) one has:                       N N N N ; ; ; ; ; ; , N N N   (11)
As previously, this formula can be straightforwardly generalized for more than two subsystems. This shows that the logarithm of the characteristic function is additive. We thus introduce

             N N N ; ; ln ; ; 1 N   
. This quantity is expected to be selfaveraging, with fluctuations decreasing as N / 1

. The rule for averaging over disorder realizations follows immediately from the additivity relation for

      N ; ;  :           N N ; ; ; ;    ( 12
)
where the bar denotes the average over disorder realizations  with identical number of domains N = N.

Some subtleties need to be fixed regarding the numerical implementation of this average. The ;    appears as the best quantity to perform a quantitative finite size analysis. We used the same ansatz

1 1 ln       bN N aN , applied to the maximum   N ; ; 0    
. A best fit is obtained for   = 0.2266, a = -0.0120, and b = 0.0735. The two parameters characterizing the N-dependence are two orders of magnitude smaller than for the maximum of the annealed or the quenched complexity (see previous section), leading to a much better estimate of the thermodynamic limit. The fact that the curves for N = 15 and 30 are quite similar means that pores of length 30 can be seen as the union of independent sub-pores of size 15. This is a quantitative confirmation of the results anticipated in Fig. 3. More precisely, one notes that the curve for N = 15 slightly differs from the other ones, but all curves for N = 20 to 30 essentially coincide. This shows that the intrinsic system size dependence due to subsystem interconnections is quite small above N = 20.

The metastable state distribution of the average system is recovered by inverse Fourier transform (the normalization is N/2). The associated complexity

                          
                d z N N N z e ; ; ; ; N LN (14)
with z any complex number (the integral is always defined). Following a route similar to the previously introduced Fourier transform, it is straightforward to show that the natural quantity to be averaged over disorder is  can be easily calculated using the Laplace method (the integral of the exponentially large term is dominated by its largest value):

           N z N
      1 o ; ; ln 1 max ; ;                     g N N N g N (15)
One recognizes in the right hand side the Legendre transform

     N g ; ;  of        N N ; ; ln 1 N  (Eq.
        ; ; ; lim Q F Q       N N 
The two approaches however differ for finite system size. The Fourier transform (or equivalently the Laplace transform) is exact, whereas the Legendre transform appears as an approximation for large N. The second method is easier to handle, and closely follows the definition of the quenched complexity in finite systems. It however introduces important contributions of order N -1 lnN and N -1 making the finite size analysis less informative. Note that in principle these contributions could be calculated. On the other hand, the Fourier transform is more involved and sensitive to numerical inaccuracies (especially when inversion back to the density of metastable states is concerned), but is perfectly suited to study finite size dependence.

Extrapolation to the thermodynamic limit

Since the intrinsic system size dependence due to subsystem interconnections has proven to be quite small, it is possible to propose a good approximation for the thermodynamic limit.

Let us consider a large system made of n independent subsystems of size N (n is for instance the number of disorder realizations used in the previous averages, but we actually can take it as large as we want to reach the thermodynamic limit). The total number of domains is nN.

Following the additivity rule for

  N ;

;     , the complexity associated to the distribution of metastable states in the thermodynamic limit is given by:

                              d i N nN nN nN N e n ; ; 2 ln 1 lim ; ; F, Q (16)
Note that this definition differs from Eq. 13:

  N ; ; F, Q    
is the thermodynamic limit complexity of a system made of the union of independent subsystems of size N. Since An analytical approach is expected to be much more efficient and precise. This can be performed for instance for the simple binomial curves obtained for the model with independent domains (ID). These curves can then be used as test functions to fit the data, and deduce the thermodynamic limit. The binomial curves are defined as follows: For N independent domains at a given chemical potential , if we suppose that each domain can be in one of two states (empty or filled), and that each filled domain contributes 1/N of the total adsorbed density (equal to 1 at saturation), then the corresponding distribution

  N ; ;    is almost independent of N,   N ; ; F, Q     is expected
  N ; ; ID  

N

is simply given by the binomial coefficients as:

       N k N k k N k N N 0 ID ) ( )! ( ! ! ; ;     N ( 17
)
where  is Dirac measure. Note that it is possible to smear out the delta functions over -bins, of width 1/N for instance, which is simply equivalent to take into account a dispersion in the amount adsorbed per domain at saturation. In this ID model, the total number of metastable states at the given  is

  N N 2 ; ID   N
. A more realistic model should however take into account the fact that some of the domains may actually be in a single state (due to the particular value of the fluid/wall interactions in these domains), and that for the domains which can be in two states (bi-stable domains), the gaslike state is generally not completely empty, and the liquidlike state has a density generally less than 1. We are thus lead to consider the phenomenological distribution

                     N N           ; ; 1 ; ; ID phen. N N (18)
where  is the proportion of bi-stable domains and  

    ,
is the interval where the density of metastable states is not zero. The normalization follows from the Dirac measure. This phenomenological distribution is thus defined by 3 independent parameters:   ,   and  .

After calculating the Fourier transform, one gets:

  2 2 cos 2 ln ; ; phen.                                      i N (19)
The three free parameters are determined by least squares method and given in Table 1. The corresponding curves for  = 30 are given in Fig. 5 and6 (dotted lines). Note that in Fig. 6, for clarity, the discrete peaks associated to the delta functions in Eq. 17-18 have been smeared

out over the bins in density. As can be seen in Fig. 5, a fairly good agreement is obtained around the maximum of   N ;

;    , while the large  behavior is not correctly reproduced.

These large  contributions are actually associated to sharp variations of the density of metastable states and neglecting them is equivalent to smoothing the data. It can be seen on Fig. 6 that these contributions are not relevant since the agreement for the density of metastable states is very good. The advantage of this treatment is that one can now give the analytic expression for the thermodynamic limit (using Eqs 16, 18 and Stirling formula):

                                                        ln ln ; ; F, Q N ( 20
)
where the N dependence originates from the 3 parameters. The results for the various system sizes are shown in Fig. 8. The remarkable feature is that the residual finite size effect is now quite small compared to the one observed for the Legendre-Fenchel quenched complexity (Eq. 7, shown as dotted lines for comparison) and properly measures the effect of interconnections between domains. In particular for  < 0.6 the thermodynamic limit behavior seems to be reached even for N = 15, in close agreement with previous conclusion. The maximum is also quite stable, very close to its theoretical thermodynamic limit 0.22557 (the difference is 1.5 percent). 20 Beside these successes, one however notes failures around

  0 ; ; F, Q    N   . Both complexities   N ; ; F - L Q    and   N ; ; F, Q    
do not seem to reach zero for the same limiting -values. This is clearly visible for the low  bounds since they are essentially independent of system size, but remains true for the high  bound. In these regions the complexities vary rapidly, and one does not expect

  N ; ; F, Q    
to be very accurate because the quality of the fit is questionable for large  values (see Fig. 5).

The same analysis and fitting procedure has been performed for all chemical potential values.

The result is given in Fig. 9. The main averaged adsorption/desorption isotherm is also given.

As can be seen,

  N ; ; F, Q    
is positive within the hysteresis and reaches its maximum at the center of the adsorption/desorption loop, where the density of metastable states is maximal. As previously shown for  = -8.55, the complexity, as given in Fig. 9, is quantitative within few percent around the line of maxima. The uncertainty grows larger when approaching the region where

  0 ; ; F, Q    N  
. This can be seen on Fig. 9 since

  N ; ; F, Q    
does not reach zero exactly on the main adsorption/desorption isotherm. 13 It is also possible to impose that   and   be explicitly the expected main average adsorption and desorption densities. In this case, one gets a quantitative determination of the quenched complexity in the vicinity of the hysteresis loop but also around its maximum. We however expect to improve the numerical treatment accuracy in order to be able to show that

  N ; ; F, Q    
vanishes exactly on the main adsorption loop. This could be achieved using better trial functions able to fit the characteristic function data on a larger range.

Conclusion

The disorder-induced adsorption/desorption hysteresis has been studied in the particular case of a simple Lennard-Jones fluid confined in a linear inhomogeneous pore. The randomly distributed physico-chemical heterogeneities, associated to domains, induce a complex free energy landscape, and consequently metastable states. Our approach proceeds in two steps: we first calculate "conditional isotherms" by molecular simulation, which describe the fluid state in a given domain as a function of the chemical potential and for the various possible fluid states in the neighboring domains. Once these data are known, algorithms are used to calculate the adsorption/desorption isotherms for the entire pore, or enumerate all possible metastable states. Their number is shown to grow exponentially with system size (number of heterogeneities), and one thus introduces the complexities (logarithm normalized to system size). The difficult issue regarding the average over disorder is discussed. One generally introduces the annealed and quenched complexities, but we show that the characteristic function, associated to the distribution of metastable states, can straightforwardly be averaged over disorder. The associated complexity yields the quenched complexity in the thermodynamic limit.

As usual in numerical approaches, a finite size analysis is required to address the thermodynamic limit. The annealed and quenched complexities exhibit substantial finite size effects, most of which being irrelevant regarding the interconnections between domains: most of these effects would arise even if the domains were independent. To disentangle these contributions, the characteristic function associated to the distribution of metastable states is shown to be the best approach. Its main property is that, by construction, it is expected to be independent of system size for N > N0, if subsystems of size N0 are independent. Its dependence with N thus solely measures the effect of the interconnections between domains of size N. The calculations for our linear pore model clearly show that the annealed and quenched complexities exhibit finite size effect, while the characteristic function is essentially independent of system size above 20 domains. This shows that the interconnections can be neglected between pore sub-regions containing more than 20 domains. This important information can be exploited to determine the thermodynamic limit, since it is possible to calculate exactly the distribution of metastable states for a system made of an arbitrarily large number of independent sub-regions of size 20, 25 or 30 domains (for which exact enumeration was previously performed). The maximum of the distribution, which can straightforwardly be obtained, is very close to the expected theoretical limit (within 1.5%).

The other cases however require to perform an inverse Fourier transform with an arbitrarily large exponent, which introduces numerical noise. We thus introduced a collection of Ndependent, analytically tractable test-functions used to perform best fits with the simulation 20 data. The obtained quenched complexity (in the thermodynamic limit) takes maximal values in the center of the hysteresis loop, in quantitative agreement (within few percent) with previous theoretical calculations performed for some particular points. The uncertainties increase in the boarder region where the complexity reaches zero. The results are however consistent with the usual picture that the quenched complexity reaches zero on the main adsorption/desorption hysteresis loop. Noticeable improvement is expected with better trial functions able to fit the characteristic function data on a larger range. This approach offers a helpful framework to perform quantitative finite size study and extrapolation to the thermodynamic limit, relevant for many hysteresis phenomena, including adsorption and magnetic materials. Table 1: Values for ,  -, and  + for various N used to fit   N ;

Tables

;    (see Fig. 5) with Eq. complexity derived from Eq. 18 with the best-fit parameters given in Table 1 for N = 30. The delta peaks are smeared over the density bins.



  This has obviously no consequences for recovering the initial distribution for a given disorder realization (inverse Fourier transform), but might introduce spurious effects when performing the average if no care is taken to get the smooth (continuous) numerical solution for the logarithm of function of . Note that this is almost always possible except when origin, which is quite rare and occurs only for large values of . Figure4gives the Fourier transforms for two disorder realizations. The low  part contains the main information on the shape of the distribution, while the high frequency part, which generally approaches the origin while turning around, is due to fast variations in the distributions, and is actually irrelevant (its average over disorder vanishes). The results for the average   Fig.5for various N and for  = -8.55. The maximum of    = 0 and corresponds to the finite-size quenched complexity associated to the total number of metastable states at  = -8.55: The imaginary part is essentially linear, the slope giving the density  which maximizes the number of metastable states (maximum of complexity).Let us now focus on the dependence of   on N if the domains are independent. More generally, if the system can be seen as the union of independent subsystems of typical size N0, then   N above N0. In that sense,   N ;

  Fig.6for N = 15 to 30. The subscript Q is justified later. As can be seen, for low values of N the numerical fluctuations are quite important (due to fluctuating values of decrease for larger N (as these contributions become negligible in Eq.13).What is the connection, in the thermodynamic limit, between introduced complexities? Figure7shows a superimposition of the various complexities for N = 30. As can be seen, quenched complexity in the thermodynamic limit.More formally, the link between quenched complexity can be rationalized as follows. We introduce the generalized Laplace Transform  LN of  N as:

.

  4), which allowed to define The logarithm of the Laplace transform being self-averaging, one expects the Legendre transform and the quenched complexity to be self-averaging too. In the thermodynamic limit, and real axes respectively. The associated complexities, defined as inverse transformations,

..

  to be very close to the true thermodynamic limit Let us examine Eq. 16 for the particular value of  corresponding to the slope of the argument of   This choice is expected to give the maximum of the complexity. For this choice obtained for  = 0 and equals the quenched complexity the expected limit for the maximum of the complexity. For other values of  the calculations are more involved and cannot be performed numerically straightforwardly.

Figure 1 :

 1 Figure 1: Left: schematic representation of a pore exhibiting disorder. The pore is made of N

Figure 2 :

 2 Figure 2: Typical example of adsorption/desorption hysteresis for a simple fluid adsorbed in

Figure 3 :

 3 Figure 3: Symbols: annealed complexity (Eq. 1) versus amount adsorbed for various system

Figure 4 :

 4 Figure 4: Upper panel: Fourier transform

Figure 5 :

 5 Figure 5: Lines: real (upper panel) and imaginary (lower panel) parts of   N ; ;   

Figure 6

 6 Figure 6:

Figure 7 :

 7 Figure 7: Superimposition of the various complexities (Eqs. 3, 7 and 13) versus  for N = 30

Figure 8 :

 8 Figure 8: Solid lines:

Figure 9

 9 Figure 9:

  to deviate from its average value  by terms of order

												1	/	N	. As a consequence,	ln	N		  	;		;	N	
	approaches	ln	N				;		;	N		in the thermodynamic limit, and one has:

Acknowledgement: Dr. E. Kierlik and Dr. G. Tarjus are gratefully acknowledged for helpful discussions. The simulations were performed thanks to grant from the Institut de Développement des Ressources en Informatique Scientifique IDRIS-CNRS, Orsay, France.


Joël Puibasset
email: puibasset@cnrs-orleans.fr

Counting metastable states within the adsorption / desorption hysteresis loop: a molecular simulation study of confinement in heterogeneous pores

Counting metastable states within the adsorption/desorption hysteresis loop: A molecular simulation study of confinement in heterogeneous pores

Joël Puibasset

I Introduction

Despite its apparent simplicity, adsorption of fluids in regular mesoporous solids is still challenging. Many experimental efforts have been made to produce such solids presenting almost perfectly cylindrical and non-interconnected mesopores, like MCM-41, 1,2 and porous silicon. 3 Due to their importance for application and fundamental studies, they have been widely characterized, and many experimental efforts have produced accurate measurements of adsorption/desorption hysteresis in such materials for a large variety of fluids. [4][5][6][7][8][9] Extensive theoretical work has provided deep understanding of many issues, but some points remain unclear even for the above mentioned regular mesopores. 10 One of the puzzling behaviors is the appearance of large hysteresis in regular mesopores even when they are closed at one end. 7,11 Wallacher and collaborators 8 have performed extensive studies showing that the pores of porous silicon most probably present an intrinsic disorder which is more efficient to produce hysteresis than the enforced modulation of the pore size during the material production. Other authors arrived at the same conclusion that the inner surfaces of MCM-41 or porous silicon most probably exhibit some roughness and irregularity. 8,[12][13][START_REF] Berenguer-Murcia | Studies in Surface Science and Catalysis[END_REF][15][16][17][18] The quenched or frozen disorder, originating in the irregular structure of the solid or the solidfluid interaction, has been shown to explain the hysteresis in the case of gas adsorption in 3dimensional structures like aerogels. [19][20][21] It was however not yet clear if such disorder could explain hysteresis in non-interconnected mesopores like MCM-41 or porous silicon (1dimensional pore structure). To answer the question, we focused on a molecular approach to catch the main features of fluid adsorption. We introduced a simplified model of surface heterogeneities based on preliminary studies showing that the modulation of the fluid/wall interaction has strong consequences on fluid adsorption in simple pore geometry, while purely geometric deformation of the pore without changing the intensity of the fluid/wall potential has minor effects. [22][23][24][25] Note that a constriction in a pore is generally associated to both a geometric deformation and a modulation of the fluid/wall interaction (the smaller the diameter of the constriction, the larger the local fluid/wall interaction), and thus plays an important role in gas adsorption or meniscus displacement. 26,27 Simple models of inhomogeneous linear pores could then be built by modulating the fluid/wall potential along the cylindrical pore axis. Molecular simulations in the grand canonical ensemble allowed to study the cases of few heterogeneities, or domains as introduced by Everett. 28 New metastable states where shown to be stabilized by the most attractive portions of the pore, and complex coexistence diagrams of the various metastable states were proposed. 29 The complexity is growing faster than the number of domains, and the time-scale of the simulation did not allow extensive studies on larger systems. A multiscale approach could however be introduced because it is possible to study each domain of a tubular pore if one knows the fluid state in the two neighboring domains (the pores are non-interconnected, and thus present a 1D chain structure). This approach allowed to study pores of micrometric length, presenting few hundred domains, much beyond computer capabilities for a brute force molecular simulation. It is emphasized that this approach completely differs from pore network approaches where simple rules are postulated for filling and emptying (pore blocking etc.). 30 Within this approach it has been shown that the presence of disorder in tubular pores induces the appearance of a large hysteresis in adsorption/desorption cycles, very similar to what is experimentally observed in disordered materials like aerogel or Vycor. 31 This behavior arises in each pore, as a consequence of disorder, and does not require any connection between pores. Naumov and co-workers arrived at similar conclusions in a simplified lattice gas approach. 32 The effect of an additional disorder between pores (like a pore size distribution) has also been considered and shown to be negligible when the pores already present an intrinsic disorder. 33 The main scanning curves obtained in the grand canonical ensemble, i.e. desorption or adsorption branches started from any given point of the main hysteresis loop, are qualitatively comparable to experimental data: they reveal the presence of a large number of metastable states which densely fill up the whole hysteresis since all scanning curves are continuous and can be close parallel for slightly different initial filling of the porous material. Studying these curves constitutes a preliminary exploration of the underlying organization of the metastable states which might contain some information on the mesoporous solid structure. However, only few metastable states can be attained by first scanning curves, 34 and a good understanding of hysteresis requires a more systematic search for all metastable states. Such studies are generally extremely difficult, especially if one wants to infer conclusions in the thermodynamic limit, i.e. the experimental situation. Such questions have been slightly more studied in other fields like hysteresis in magnetic materials, but recent works have shown that their relevance to adsorption hysteresis cannot be underestimated. 35,36 As shown in previous papers, the fluid may be in essentially two states in each domain. This allows a systematic study of all possible states of the confined fluid, at least for small enough systems. An algorithm (given below) determines which states are metastable and their location in the - plane (chemical potential versus amount adsorbed). However, due to the exponential increase of the number of possible states with the number of domains, the study is limited to quite small systems with typically less than 35 domains. Two issues need then to be addressed: 1) the results actually strongly depend on the particular distribution of heterogeneities in the porous material. One thus has to somehow perform an average over disorder, i.e. over many possible pore samples or "realizations". 2) the thermodynamic limit of infinitely long pore has to be discussed. These issues cannot be addressed rigorously in our numerical approach, except when exact calculations are possible (used as benchmark).

Reasonable approximations allow however to determine appropriate methods to perform average over disorder, and a systematic finite size analysis allows to infer trends for the thermodynamic limit. Note that in a typical MCM-41, pores are several micrometers long, corresponding to approximately 10 2 -10 3 domains. The thermodynamic limit might thus not be the only relevant case. In this first study of the metastable state structure we focus on their density distribution inside the hysteresis loop, and the associated complexity, i.e. the logarithm of the density normalized to the system size. Two examples are studied, one with usual continuous adsorption/desorption isotherms, and one exhibiting a sharp jump during adsorption. They are shown to present important differences in their distribution of metastable states, suggesting the existence of an underlying transition.

The paper is organized as follows: In section II, we introduce the pore model, the algorithms used to determine the isotherms and metastable states, and briefly introduce general results regarding disordered systems, in particular the annealed and quenched averages, where the variables associated to the source of disorder, originating in the pore irregularities in our case, are either free to fluctuate or frozen. In section III we present the isotherms and all metastable states for various pore samples. The distributions of the number of metastable sates over disorder realizations and for a given chemical potential are shown and discussed. The quenched complexity is shown to appear as self-averaging. In section IV and V we properly define and calculate, in finite size, the annealed and quenched complexities as a function of the chemical potential , and as a function of  and the amount adsorbed . Finally, in section VI, the case of a porous material exhibiting a sharp jump during adsorption is studied. We finally conclude by comparing the two systems and the qualitative differences observed for the quenched complexity.

II Model and methods

Our model is designed to catch the main features of simple fluid adsorption in heterogeneous tubular pores like MCM-41 or oxidized porous silicon. The pores, of diameter around 3 nm, are supposed to be made of pure silica. The interactions are supposed to follow the Lennard-Jones potential. The parameters are those for argon-argon and argon-oxygen species to describe fluid-fluid and fluid-silica interactions. 31 All thermodynamic quantities are reduced using fluid-fluid interaction parameters. Heterogeneities are introduced by modulation of the fluid/wall interaction along the pore axis, within +/-30%. The modulation factor, denoted h, is the ratio between the local fluid/wall interaction intensity to its reference value. For instance, h = 0.85 means that the interaction is 15 percent lower. It is chosen to be a step function along the pore axis. Each step corresponds to a domain as defined by Everett 28 possibilities need to be considered. Figure 1 shows a typical set of 4 such "conditional" isotherms obtained in a cylindrical domain. After these calculations are performed for all possible domains, the fluid properties can be easily determined in the whole pore, according to the algorithms described below. This approach allows a decoupling between the atomistic level calculations performed by Grand Canonical Monte Carlo molecular simulation in each domain, and the large-scale calculation of the fluid state in pores of micrometric length. More details on the model can be found in refs [ 31 ] and [ 33 ]. With this approach we are able to study fluid properties in very long pores at the atomistic level, or perform extensive studies like counting metastable states in systems of reasonable sizes, which is impossible otherwise. In the rest of the paper, we consider only pores open at both ends: as a consequence, the first and last domains are always in contact with the bulk gas state.

The adsorption/desorption isotherms in a complete pore are calculated as follows: For low , the adsorption in each domain is single valued, and does not depend on the fluid state in the rest of the pore. The amount adsorbed is simply the sum over all domains. Note that in this case, all domains are in the gaslike state and follow the conditional isotherm with G-G boundaries. The chemical potential is then increased by small amounts until one domain reaches its limit of stability. This domain then jumps to the liquidlike branch. This has consequences for the two neighboring domains, which boundaries are now G-L and L-G.

Their state has then to be reconsidered according to their G-L isotherms. They may stay in the gaslike state, meaning that the fluid condenses only in the initial domain, or switch to the liquidlike state, meaning that the meniscus created in the initial domain propagates to neighboring domains. In this case, the next domains need now to be reconsidered too, and so on until the meniscus stops, corresponding to the end of the avalanche. The amount adsorbed is given by the sum over all domains. The chemical potential is then further increased, and the same procedure is applied until saturation is reached. Desorption and scanning isotherms are obtained in a similar way.

We now present the method used to determine all metastable states for any given realization of the disorder in the pore and for a given chemical potential . As previously explained, the fluid properties are entirely determined by the chemical potential  and the fluid state (G or L) in each domain. One then expects at the most 2 N possible fluid states, where N is the number of domains. However, due to the interdependence between neighboring domains, many configurations are actually impossible. For instance, a domain in between a gaslike and a liquidlike domain, actually has only one (instead of two) possible state entirely determined by the chemical potential since G-L isotherms are single valued. But many other situations may appear, depending on the exact forms of the pre-calculated isotherms. From a physical viewpoint, these configurations correspond to situations where the fluid state is unstable and has to evolve towards another state, by filling or emptying some domains. On the other hand, the configurations which are self-consistent are stable against thermal fluctuations, which are taken into account in the pre-calculated isotherms, and correspond to the metastable states.

For any given , the algorithm explores all possible states, and extracts the few ones that are self-consistent, i.e. metastable. It is observed that any given metastable state can exist on a range of , which can easily be determined by slowly decreasing (increasing)  until the lower (upper) stability limits min (max) are reached. The corresponding () curves can be drawn in the - plane to represent each metastable state.

The number of possible fluid states to be examined increases exponentially with system size.

As a consequence, our algorithm cannot consider pores with more than 35 domains for few disorder realizations, 30 domains for 5×10 3 disorder realizations, 25 domains for 5×10 4 disorder realizations, and 20 domains or less for 10 5 disorder realizations. It will however be shown below that the results indicate that the thermodynamic limit can be reasonably extrapolated from these system sizes. The last difficulty is how to perform averages over disorder realizations? Each realization produces its own hysteresis. How could we manage to get an "average" hysteresis? For a given chemical potential far from the hysteresis region, the fluid state is single-valued and given by the amount of fluid adsorbed. This quantity being additive, one expects that the "average" fluid state over disorder be precisely the averaged amount adsorbed. Note that the additivity hypothesis relies on the fact that for large enough systems, the amount adsorbed in any macroscopic region of the porous solid is not very different from what would be obtained if this region was isolated from the rest of the material.

However, in the hysteresis region, the amount adsorbed is not uniquely defined, and can take any value between the so-called main adsorption and desorption branches. One can easily average these main branches over disorder to define the "average main hysteresis loop". However, characterizing the distribution of metastable states within the hysteresis and performing the average is a much more difficult task, to be explained now.

Let us define the density of metastable states in the - plane ;   N as a density is necessary in our case (continuous values for ) and differs slightly from the usual definition used in Ising models where  is discrete. This has no consequences for the definition of the corresponding complexities (see Eqs. 1 and 2 below) since the differences vanish in the thermodynamic limit. We also define the number of metastable states at a given  as

          d N N ; ; ; N N . The total number of metastable states is denoted   N N . Note that despite similar notation "
N ", these quantities are not identical and the parentheses will help to differentiate them.

Furthermore, one has to keep in mind that

  N N is not the integral of   N ;  N over 
because all metastable states exist over a finite range of , and thus each metastable state would be counted several times.

The density of metastable states

  N ;

;   N not only depends on the statistical properties of the disorder in the pore, but is also strongly sensitive to the particular disorder realization in the pore sample. It is then natural to introduce an average over disorder. As usual in disordered systems, it is interesting to consider both the average and typical number of

metastable states:   N ; ;   N and  ] ; ; ln exp[ N  

N

where the bar denotes the average over disorder. 35,36 Since these quantities are both expected to scale exponentially with N, one introduces the corresponding annealed and quenched complexities:

    N N N ; ; ln 1 lim ; A     N      (1)     N N N ; ; ln 1 lim ; Q     N      (2)
As will be shown in next section, the quenched complexity is of particular interest because it is expected to be self-averaging (its fluctuations scale as

N / 1
), which is not the case for the annealed complexity.

III Isotherms and metastable states

The main adsorption/desorption isotherms and all metastable states are calculated for various pore length and disorder realizations. Some examples are shown in Fig. 2 for N = 5, 10 and 15. All isotherms present two reversible portions at low and high , and a hysteresis in between. The saturation branch is essentially identical in all cases, whereas the initial adsorption branch exhibits some differences which reflect the variations in the average fluid/wall interaction in the whole pore. These variations are expected to decrease with the system size. As can be seen, the hysteresis is also strongly dependent on the particular realization of the disorder. One can notice differences in width, height, and slope of the adsorption and desorption branches.

The metastable states appear as portions of lines corresponding to their metastability range.

These lines can be seen as the projection in the - plane of the local minima in the free energy landscape. As expected, the slopes are always positive due to increasing film thickness with . These lines are almost straight and essentially parallel, meaning that despite possibly large differences in fluid configuration between two adjacent metastable states, the response of the system to a chemical potential variation is essentially the same in both cases. The small difference appears as a balance between the number of domains in the gaslike and the liquidlike state, expected to be of order 1/N. However, large variations can be observed in their distribution in the - plane depending on the disorder realization. Nevertheless, in all cases, all metastable states fall within the main hysteresis loop. It is actually observed that the main hysteresis constitutes the external envelope of the metastable states. This feature is observed in other systems and seems to be quite general. 35,36 .

Another interesting quantity is the total number of metastable states   N N , given in each panel of Fig. 2. The first point to be noticed is that this number grows faster than N, but is significantly smaller than the total number of fluid states 2 N , especially for large N. It also takes a wide range of values depending on the pore disorder realization. Its distribution has been numerically investigated for various N. As can be seen in Fig. 3 (first and second panels)

the distributions are quite large, with long tails. As previously noticed in undulated pores, 29 the typical number of metastable states grows quite fast with the number of domains. This is why the distributions cannot be drawn on the same panel. The average and standard deviation were also calculated for 30  N

, and shown to be compatible with an exponential behavior.

Furthermore the width of the distribution remains of the same order as its average value. On the other hand, the distributions for

N

. As a consequence, the quenched complexity appears as the natural self-averaging quantity to characterize the number of metastable states.

IV Calculation of

 

N ; ln  N and   N ; ln  N
We now focus on the number of metastable states for a given chemical potential . Since numerical simulations are performed in finite systems, we define finite size annealed and quenched complexities for a given  as:

    N N N ; ln 1 , A   N   (3)     N N N ; ln 1 , Q   N   (4)
where N is the pore length. Their corresponding thermodynamic limit  

  N are denoted as      A and      Q .
Let us first focus on the particular choice  = -8.55 corresponding to the central region of the hysteresis (see Fig. 2). The finite size complexities have been calculated for N = 15, 20, 25 and 30 with respectively 10 5 , 10 5 , 5×10 4 , and 5×10 3 disorder realizations (circles in Fig. 4). The data exhibit a small dependence with the system size, which is well described by the ansatz

        1 A/Q 1 A/Q A/Q A/Q ln ,         N b N N a N    
involving three free parameters which depend on . 36 A least-squares fit to the data for  = -

yields

  A = 0.2555, A a = 0.0313, A b = -7.538×10 -3 and   Q = 0.2266, Q a = -0.0120, Q b = 0
.0735 (the best fit are shown as dashed lines in Fig. 4).

For the particular choice  = -8.55 an analytical calculation is possible and can be used to check the accuracy of the previous system-size dependence analysis (see Appendix A). The results are shown as red solid lines in Fig. 4. As can be seen, the data obtained by exact enumeration of all metastable states are in close agreement with the expected values, showing the good convergence in terms of average over disorder. We are now in a position to address the thermodynamic limit. One can first perform the numerical calculations for the largest possible value of N, which gives good results due to the weak dependence with N, but the safest is to perform the analytical derivation (see Appendix B). The results are

  25715 . 0 55 . 8 A       and   22557 . 0 55 . 8 Q      
. They are given as horizontal straight lines in Fig. 4. These exact calculations show that the numerical finite size analysis was correct within 1%, and is likely to be reliable even when analytical calculations are not possible.

This systematic enumeration of metastable states has been extended to other values of the chemical potential. The results are displayed in Fig. 5 (symbols). The previous theoretical derivation of the annealed and quenched complexities can be extended to any < -8.53 (see Appendices). The corresponding theoretical thermodynamic limits are given as lines in Fig. 5.

As can be seen, the simulation results exhibit a weak dependence with system size and are quite reliable, even in regions where the complexities vary rapidly. The data show that both complexities are not zero for  between approximately -8.67 and -8.25, and exhibit a maximum around  = -8.58. As expected,

           A Q
. A remarkable feature is that both complexities reach zero for the same chemical potential. It corresponds to the situation where all domains are of the same type. The theoretical derivation shows that both complexities reach zero when the probability pA reaches zero at  = -8.6757: straightforward calculations show that at first order

     1 1 A A        p and     1 A Q ln   p    , for   0 A   p . A
similar behavior is expected around  = -8.25. Note that these points correspond to the lower and upper closure points of the main hysteresis loop. The same is observed in Ising chain models. 35 V Calculation of

  N ; ; ln   N and   N ; ; ln  

N

As previously, we introduce the annealed and quenched complexities for finite systems. The annealed complexity associated to the distribution of metastable states versus  and  is defined as follows:

    N N N ; ; ln 1 , ; A     N   (5)
where the bar denotes average over disorder. This quantity has been calculated on a discrete grid in the - plane, for the same -values as previously, and for 24 -values between 0 and 1. Our definition of the complexities based on the density of metastable states thus differs from the one used in Ising models by

 ln 1   N
where the bin width  = 1/24 in our case.

Figure 6 shows the results for N = 25. This system size allows a good compromise between large N and a large number of disorder realizations (5×10 4 ). For comparison, the average over disorder of the main hysteresis loop for N = 25 has been superimposed in the top view. As can be seen, the annealed complexity presents a typical shape with maximal values within the region of hysteresis. It can also be noted that the region where the annealed complexity is positive extends well beyond the main adsorption/desorption loop. It also exhibits a spurious finite region where a single fluid state is expected whatever the disorder (for instance for  < -8.7 and  > -8.2). In between these values, visual inspection shows that the hysteresis defined by the main loop qualitatively corresponds to

     ln , ; 1 A     N N
(green to red region in Fig. 6), i.e. to an average number of metastable states larger than 1 per bin. How does this change in the thermodynamic limit? The system size dependence is illustrated in Fig. 7 for  = -8.55. The first point to be noticed is that the system size dependence is significantly more

pronounced for   N , ; A    than for   N ;
A    symbols are superimposed in Fig. 5, not in Fig. 7. However, as previously, it is possible to extrapolate the thermodynamic limit using the

N-dependent ansatz 1 1 ln      bN N aN
, where the three parameters depend on  and .

When ;   N is exponentially large in the thermodynamic limit and thus its integral over , giving      A , is dominated by its largest term. A thermodynamic limit curve can thus be proposed (Fig. 7, dashed line) exploiting the system size dependence. In particular, one observes that the densities for which

  0 , 55 . 8 ; A     N  
are essentially stable with N. Note that, in the thermodynamic limit, the regions defined by

     ln , ; 1 A     N N or   0 , ; A   N  
are identical and clearly extend well beyond the hysteresis defined by the main adsorption/desorption loop. This conclusion is supported by theoretical works. 35,36 The finite-size quenched complexity associated to the distribution of metastable states versus  and  may be introduced as follows:

    N N N ; ; ln 1 ; ; Ch Q     N   (6)
where the bar denotes average over disorder, and the logarithm is modified as proposed by Cherrier and collaborators: 37 since for finite systems, the number of metastable states can locally be equal to zero, the logarithm diverges. It is thus arbitrarily set to zero in such situation. This introduces non-controllable errors in the evaluation of the quenched complexity. However, the results are expected to be reliable in the central region of the hysteresis where the density of metastable states almost never vanishes. The results are shown in Fig. 8 (for N = 25) as a function of  and , for

     ln ; ; 1 Ch Q     N N
so that the typical number of metastable states is always larger than 1 per bin. The average main hysteresis loop for N = 25 has been superimposed in the top view. As can be seen, the largest values for the quenched complexity are reached within the hysteresis loop, but, as previously, the region determined by

     ln ; ; 1 Ch Q     N N
extends beyond the main adsorption/desorption loop. Note however a significant qualitative improvement for  < -8.68

and  > -8.28. The system size dependence is illustrated in Fig. 9 for  = -8.55. As previously, a best fit can be obtained with the ansatz

1 1 ln      bN N aN
, and the limiting value of the maximum for  = -8.55 is compatible with

  55 . 8 Q     
. A putative thermodynamic limit can be drawn (Fig. 9, dashed line). The figures show that the quenched complexity reaches its largest values inside the hysteresis in the thermodynamic limit. However, since

  N ; ; Ch Q   
is not reliable when it approaches zero, it is impossible to discuss whether the quenched complexity exactly vanishes on the main adsorption/desorption loop in the thermodynamic limit. Another approach is required.

To circumvent the difficulty raised previously, it has been proposed 36

  N N ; ; ln 1    N  :                     g N N N g ; ; ln 1 max ; ; N (7)
where  denotes a given disorder realization. Note that this quantity is always defined for any real value g, even if ;    the -value which maximizes the quantity between parentheses in Eq. 7. One has:

            g N N N g    ; ; ln 1 ; ; N (8)
Averaging over disorder, one gets:

      N g g N N N g ; ; ; ; ln 1 ; ;           N (9) where     N g N g ; ; ; ;      
. A quenched complexity can thus be defined for a finite system as: ) selects the smallest (respectively largest)  for which the density of metastable states is not zero, corresponding to the main adsorption (respectively desorption) branches for each realization of disorder. As a consequence, the envelope of the projection in Fig. 10 also corresponds to the average over disorder of the main hysteresis loop. The agreement is quite good within the uncertainties of our numerical approach: for each of the 24 -values used to perform the calculations in the range [-8.8;-8.1],

    ) ; ; ( ; ); ; ; ( ; ; F - L Q N g N N g N             (10)
  0 25 ; ; F - L Q    N  
superimpose with the main hysteresis within few percent; the visible discrepancies or roughness in Fig. 10 arise as a consequence of interpolation between the data. In order to discuss the thermodynamic limit, the system size dependence is shown in Fig. 11 for  = -8.55. It is slightly more pronounced than in previous cases, but still well fitted with the ansatz

VI Systems with macroscopic transitions

When the disorder is moderate, the hysteresis loop associated to the capillary condensation exhibits an unusual behavior. One observes vertical branches for the main adsorption or desorption curve, generally associated to a true out-of-equilibrium transition characterized by a macroscopic jump in the hysteresis. 19 Metastable states are expected to be found in the hysteresis loop, but distributed in such a way that some regions are totally depleted. 34,36,38 20 Such macroscopic transitions are observed in our system, for instance, when the Gaussian distribution for the heterogeneities is replaced by a flat and truncated distribution. 33 We suspect that the underlying mechanism is quite different in both situations, for the simple reason that a true out-of-equilibrium transition is not expected in materials exhibiting a onedimensional pore network. The important point is most probably the truncation of the distribution. Gaussian distribution being strictly positive, all heterogeneity values between 0.7 and 1.3 may be found in a given pore. This is not the case if one considers a distribution which is strictly zero for heterogeneity out of a given interval. For instance, we now consider a flat distribution over [0.9,1.1]. In this case, the extreme values of the heterogeneity are not observed at all in the pore, while they could be found in our previous example with a Gaussian distribution. From this point of view, the level of disorder is lower.

For any realization of the disorder, a macroscopic transition is observed during adsorption (3 examples are given in Fig. 12 for N = 25 domains), while desorption occurs in several steps. In all cases, the sharp adsorption is around  = -8.45, which corresponds to the gas to liquid transition in elementary pores with heterogeneity around its median value. 31 The small variations around this value yield a rounding of the jump in adsorption after averaging over disorder. However, the larger the system size, the steeper the jump (see Fig. 13). A macroscopic transition is thus suspected in the thermodynamic limit.

All metastable states have been determined and shown in Fig. 12. All of them fall within the main hysteresis loop, but does not distribute uniformly: no metastable state could be found in the right hand part of the hysteresis, except for the gaslike and liquidlike branches. This is explained by noticing that when transition occurs somewhere in the pore for  around -8.5, the meniscus propagates freely in the whole pore (absence of strongly "repulsive" domain to stop the meniscus). The system is thus either in the gaslike or the liquidlike state. On the other hand, an exponentially large number of metastable states is found in the left half-region of the hysteresis. Their density is largest around  = 0.6, and decreases monotonously while approaching the gaslike and liquidlike branches.

A quantitative analysis including the average over disorder has been performed by calculating the quenched complexity

  N ; ; F - L Q   
. The result is shown in Fig. 14 for the largest system size N = 30 and 5×10 3 disorder realizations. A smoothly varying complexity is obtained, with a single maximum for any given chemical potential. As previously noticed by visual inspection on various disorder realizations, the quenched complexity reaches monotonously zero while approaching the main hysteresis loop for  < -8.52 (Fig. 14, second panel). The interesting behavior is for  > -8.52: the quenched complexity is not defined since no metastable states were observed in the 5000 disorder realizations used to perform the average, except for the two gaslike and liquidlike branches. In other words, in this system, the quenched complexity is shown to vanish strictly inside the main hysteresis loop. Before going farther, let us discuss the thermodynamic limit. We have performed the same finite size analysis as in the previous system. The maximum values reached by the quenched complexity slightly decrease while approaching the thermodynamic limit, but, as previously, the region where

  0 ; ; F - L Q   N   is essentially independent of N. In particular, one has   0 ; ; F - L Q   N  
for  > -8.52 even in the thermodynamic limit, and thus the quenched complexity vanishes strictly inside the main hysteresis loop.

Such a behavior has been reported in the random field Ising model (RFIM), 35,36 with however noticeable differences. For instance, in our case, the depleted region extends to all density values between the gaslike and the liquidlike states, while, in the RFIM, metastable states could be found in regions close to the main gaslike and/or liquidlike branches (materialized by an ellipse in Fig. 14). The origin for this difference may arise due to the treatment of the interdependence between adjacent domains, or to the particular choice for the fluid/wall interactions and the associated distribution of heterogeneities. It is however suspected that the difference is more fundamental. Remember that the RFIM with the same 1-dimensional topology as our pores predicts a continuous curve whatever the disorder strength. The appearance of a jump in adsorption in our system is closely connected to the truncation of the distribution of heterogeneities. As previously mentioned, the region of chemical potential  > -8.52 corresponds to a situation where a liquid/gas meniscus can propagate freely through any heterogeneity in the pore, the truncation of the distribution resulting in the absence of strong enough heterogeneities to pin the meniscus. Furthermore, the gaslike state is metastable while the liquidlike state is stable. As a consequence, any state of intermediate density is unstable and immediately turns into liquidlike state. However, as for the RFIM, the main adsorption loop corresponds to the grand-canonical envelope of the metastable states, where

  0 ; ; F - L Q   N  
; it is reminded that along the gaslike and liquidlike branches, the quenched complexity equals zero in the thermodynamic limit.

If the pores are closed at one end, the sharp transition to liquidlike state shifts towards lower value 33 ( = -8.52), which corresponds precisely to the lower limit of the depleted region. This is in agreement with our previous interpretation: the pore end helps to nucleate the liquid/gas meniscus and suppresses the metastable gaslike branch above  = -8.52. In this situation, the main adsorption/desorption loop is still the envelope of the metastable states, except that one observes a sharp transition during adsorption without depletion of the hysteresis loop (i.e. many intermediate states can be found closely along the sharp transition). This is still different from the RFIM results. It is however emphasized that our numerical study cannot exclude a very small depleted region around the transition. Using a finite reservoir may help to investigate these questions. 34,39

Conclusion

We have studied numerically the density of metastable states embedded in the main adsorption/desorption hysteresis loop for a simple fluid confined in a single disordered tubular pore. Very few works can be found on such investigations, essentially in lattice gas models.

The originality of this study is that the model is based on a molecular simulation approach of the adsorption of a Lennard-Jones fluid in a heterogeneous linear pore. The substrate is divided into domains, and the fluid adsorption properties are pre-determined for each domain by molecular simulation as a function of the fluid state in the neighboring domains, thus explicitly taking into account the effect of interconnections between domains. This approach completely differs from pore network models where ad-hoc rules for filling or emptying are postulated (pore blocking, etc.). The adsorption properties of the entire pore can be then calculated straightforwardly in our approach, which allows to characterize all the numerous metastable states. The annealed and quenched averages over the disorder originating in the physico-chemical heterogeneities in the pore are introduced to calculate the corresponding complexities. Finite size effects have been studied extensively. Interestingly, they are shown to be quite small concerning the total number of metastable state at a given chemical potential, and the data can be safely extrapolated to the thermodynamic limit. The annealed and quenched complexities, characterizing the density of metastable states, present more pronounced finite size effects; the thermodynamic limit is inferred from the expected 1/N and lnN/N behavior. When possible, the numerical results have been compared to exact analytical calculations, validating the average over disorder and the finite size analysis.

One observes that the density of metastable states grows exponentially with system size and presents a pronounced maximum in the central region of the hysteresis. The annealed complexity remains positive in a region which extends beyond the main average hysteresis loop, confirming that the average density of metastable states is unsuited to characterize the hysteresis in the thermodynamic limit. On the other hand, the quenched complexity, characterizing the typical number of metastable states, is shown to be strictly positive only inside the hysteresis and to vanish on the main grand canonical adsorption and desorption branches. Equivalently, the main adsorption/desorption hysteresis averaged over disorder corresponds to the envelope of the typical distribution of metastable states. An interesting point concerns the comparison of two systems with different distribution of heterogeneities, one Gaussian, the other one flat and truncated. The main difference is that in the second one the high values of disorder are strictly suppressed. In the first case, a smooth hysteresis is observed, and an exponentially large density of metastable states is found everywhere in the loop. In the second case, a macroscopic jump is observed during adsorption, and a depleted region with no metastable states appears in the vicinity of the transition within the hysteresis loop. The interpretation is as follows: The distribution of heterogeneities being truncated, the less attractive domains are not strong enough to impede the propagation of a meniscus. As a consequence, the possible states are just liquidlike and gaslike. The first one is actually the stable state. The gaslike state is metastable due to the truncation of the distribution: the most attractive regions in the pore are not strong enough to nucleate a liquid meniscus. When the pores are closed at one end, the depleted region and the hysteresis shrink, as a consequence of the facile nucleation of a meniscus in the bottom, while the transition remains sharp and macroscopic. Further investigation is needed to assess whether the depleted region completely disappears or remains as a thin area close to the transition. Significantly larger systems are required to minimize finite size effects and reach the required resolution for the quenched complexity.

Appendices:

These appendices present the main steps in deriving the analytic formulae for the annealed and quenched complexities when the domains can be separated in two types (the so-called Atype and B-type, see below); i.e. for  < -8.53.

A: Annealed and quenched complexities for finite size system

Examination of the various possible conditional isotherms induced by disorder shows that for  < -8.53, all domains necessarily fall into one of two types (see Ref [ 31 ] and Fig. 1). Type A corresponds to domains for which

2 coex     
occurring when the disorder parameter h > 0.945 (h is the modulation factor of the fluid/wall interaction). 31

      R A L R A G R A N N N N N N   (A1)       R A L R A G R A G 1 N N N N N N    (A2)       R A L R A G R A L 2 1 N N N N N N    (A3)
and For N domains, the probability to have R A N A-regions is given by:

    1 1 R A L R A G    N N N N = 1.

Straightforward calculations show that

  R A R A 2 2 1 1 R A N N N       N (A4)
                        1 2 0 1 B A 1 1 2 1 1 ; R A e B R A e B R A A A A e B e B A R A A R A R A N N N N N N N N N N N N N N N N N N N N p p P (A5)
where  

1 1 A R A   N N
is the number of ways to split the A N domains into R A N A-regions, and

  1 2 e B A R A     N N N N
is the number of ways to split the

5 . 0 0 R A   n , N n n / 1 2 1 0 R A e B    
, and

N n n n n / 1 1 R A e B A R A      . ( B1
)
Since now on, all quantities negligible in the thermodynamic limit are dropped. Using Stirling formula:

        ) ( o 1 ln 1 ln M x x x x M M xM e M        (B2)
and replacing the discrete sums by integrals, the probability   N N ; R A P can be written as:

    A e B A e B R A e B R A R A ) ( o , , , , ; dn dn N n n n p p N e n N N n N B A N       D P  (B3)
where the integral is performed over the domain D defined by Eq. B1 (the terms 1/N are neglected) and

        e B A e B A A A A A e B R A 1 ln 1 ln ln 1 ln , , , , n n n n n n p n p n n n n p p A B A B A                   R A R A R A e B A R A e B A R A A R A A ln 2 1 ln 1 ln n n n n n n n n n n n n            (B4)
The dominant term of the integral is given by the maximum of  over D .  is strictly decreasing with e B n and thus the maximum is for e B n = 0, i.e. when the pore has A-domains at

     B A R A * A * A * A R A * A 1 1 p p n n n n n n      (B5)
For a given R A n , the l.h.s is an increasing function from 0 to infinity for

    ) ( o , , ; R A max R A R A N n p p N e N N n N B A N       P (B7)   R A max , , n p p B A  , seen as a function of R A n exhibits a maximum for * R A n (the corresponding value of * n A is denoted ** n A ) such that     0 , 0 , , , , , * * A A * R A R A * R A R A ; R A A e B R A R A R A max         n n n n B A n n B A n n n n p p dn n p p d   (B8)
The first equality holds since by definition

  0 , 0 , , , * * A A A A e B R A     n n B A n n n n p p  . The result is      2 * R A * R A A * A A 1 n n n n n ** **    
, which gives, using Eq. B5:

B A * R A p p n  and A ** p n  A .
Inserting these values into Eq. B6 shows that this maximum is actually zero. Hence,

  0 , , R A max  n p p B A 
and the leading term of the probability given by Eq. B7 is correctly normalized. Note that the probability

  N N ; R A P is thus essentially Gaussian around B A * R A p p n  :     ) ( o 2 ; 2 B A R A R A R A N p p n N e N N n N N         P (B9) with     2 B A 2 R A R A max 2 * R A R A / , ,        p p n n p p n n B A   > 0.

Annealed complexity:

The annealed complexity is given by

                          ;max A A 0 R A R A A ; ln 1 lim ; ; ln 1 lim ; R N q N N N N N N N N N P N     (B10)
Turning the discrete sum into integral gives:

        R A 1 R A R A max 1 0 R A R A ) ( o ln , , ; ;max A A n d N n n p p N e N N N B A N q R              N P (B11)
The integral of the exponential is again dominated by its largest term, obtained here for

Q R; A R A n n  (the corresponding value of * n A is denoted Q ; A * n ) such that:          0 1 ln ln , , , , 2 Q R; A Q R; A Q *; A Q R; A Q *; A 1 1 R A R A max R A R A max Q R; A R A Q R; A R A                     n n n n n n n p p dn n p p d n n B A n n B A     (
                       B A Q R; A 5 1 2 1 1 1 4 5 1 p p n (B13)                 B A B A B A Q *; A 5 1 2 1 4

Quenched complexity:

The quenched complexity is given by:

                  ;max A A 0 R A R A Q ln ; 1 lim ; ; ln 1 lim ; R N q N N N N N N N N N P N     (B15)
Using Eqs. A4 and B7 and the fact that the integral of the exponential is dominated by its largest term in the thermodynamic limit (obtained for

* R A A R; A R A n n n  
), one gets.

    1 B A 0 R A R A 1 Q ln ; 1 lim ln ; ;max A A     p p N N N N R N q N         P (B16)

of approximately 7

 7 nm. The values of the steps are uncorrelated and follow a Gaussian distribution centered on 1.0 and half-width 0.2 in the interval [0.7; 1.3]. Interactions being short range, typically less than a domain length, the fluid properties in a given domain are entirely determined by the fluid configurations in the two neighboring domains. Since fluidfluid correlations decrease very rapidly, the relevant quantity in the neighboring domains is, on average, the fluid density. This is entirely determined by the chemical potential , except if the adsorbed fluid may exist in several states. This is precisely the case for a fluid confined in a tubular domain, which, for some values of , may be in the gaslike (G) or liquidlike (L) state, corresponding respectively to fluid adsorbed only at the wall or filling the whole domain. As a consequence, the fluid properties in a given domain and for a given  are entirely determined by the fluid state (G or L) in the two neighboring domains. At the most 4

  a given chemical potential , the number of metastable states with adsorbed density in the interval d around .Note that this definition of   N ;

  to Gaussians (see Fig. 3, third panel). The average value over disorder   N N ln is linear with N, and the standard deviation of

  performed for the maxima, occurring around  = 0.53 for  = -8.55, one gets a limiting value   very close to the previously calculated  

  size, disorder dependent complexity

N

  is zero for some values of . To see the connection with the quenched complexity, let us denote   N g ;

  Fig.10as a function of  and , for N = 25. As previously, the largest values are reached inside the hysteresis loop. But, contrary to previous cases, the quenched complexity exactly vanishes on the main adsorption/desorption hysteresis loop, shown as a solid line in the second panel of Fig.10. This is expected: it can be seen by inspection of Eq. 7 that the limit   g (respectively   g

.

  Fig.11, dashed line). An interesting feature is that the locations where

= 2 .

 2 Such domains can be in two states (G or L) only if both neighbors are gaslike. Otherwise, these domains are necessarily in the liquidlike state. Type B corresponds to domains for which 945. They can be in two states (G or L) only if both neighbors are liquidlike. Otherwise, these domains are necessarily in the gaslike state. For our particular choice of , any pore can thus be seen as a chain of A's and B's, and the total number of domains is N = NA + NB where NA and NB denote the number of domains A and B. Knowing the relation between the domain type and h as well as the distribution of h, one can deduce that the probability of having a domain A is pA = 0.625 (and pB = 0.375). Considering the properties of domains A and B, if two A's or two B's are in contact they are necessarily in the same fluid state. This can be generalized to any group of A's or B's. The pore can thus be seen as a chain of Aand Bregions of various lengths. R A N and R B N denote the number of A-regions and B-regions, and of regions. Any metastable fluid state in the pore can be described by a chain of R N fluid states corresponding to the R N regions of the pore. It is also noted that any B-region at one pore end is necessarily in the gaslike state. e B N denotes the number of domains B belonging to B-regions which are at the pore ends. If the pore has an Astates for a chain of R A N A-regions (the number of B-regions is necessarily R A N -1, R A N , or R A N +1 depending of the nature of pore ends. Simple inspectionshows that Let us consider the last A-region of the pore fluid state in this region is either G or L. We denote liquidlike in this region. One has the following relations:

  disorder configuration. We now calculate their probability.

  the number of ways to append the e B N B-domains at both pore ends. The summations are performed over all possible choices for e B obtained numerically, which allows to calculate the quenched and annealed complexities as a function of N for  = -8.55. The results are shown as red solid lines in Fig. 4

  Figure captions:

Figure 1 :

 1 Figure 1: Schematic representation of the typical set of 4 "conditional" isotherms

Figure 2 :

 2 Figure 2: Black curves: main adsorption/desorption hysteresis. Red lines: all metastable states obtained by systematic exploration. The total number of metastable states N is given in each figure. From left to right: increasing number of domains (N = 5, 10, and 15). For each N,

Figure 3 :

 3 Figure 3: Upper panels: Distribution of the number of states N found in pores made of N =

Figure 4 :

 4 Figure 4: System size dependence of the annealed and quenched complexities for  = -8.55.

Figure 5 :

 5 Figure 5: Quenched and annealed complexities (Eq. 3 and 4) versus chemical potential. Symbols: exact enumeration. Circles: N = 15 (10 5 disorder realizations); squares: N = 20 (10 5 disorder realizations); up triangles: N = 25 (5×10 4 disorder realizations); down triangles: N = 30 (5×10 3 disorder realizations). Solid lines: theoretical thermodynamic limit (see text and

Figure 6

 6 Figure 6: Annealed complexity

Figure 7 :

 7 Figure 7: Symbols: annealed complexity versus amount adsorbed  for N = 15, 20, 25, 30 and

Figure 8

 8 Figure 8: Quenched complexity

Figure 9

 9 Figure 9: Quenched complexity

Figure 10

 10 Figure 10: Quenched complexity

Figure 11

 11 Figure 11: Quenched complexity

Figure 12 :

 12 Figure 12: Main adsorption/desorption hysteresis (black curves) and all metastable states (red lines) obtained by systematic exploration. The total number of metastable states N is given in each panel, corresponding to three realizations of the disorder. The number of domains is N

Figure 13 :

 13 Figure 13: Enlargement of the average adsorption isotherm around the macroscopic jump, for various system sizes: N = 10, 20 and 25 (respectively 10 5 , 10 5 , and 5×10 4 disorder

Figure 14

 14 Figure 14: Quenched complexity

3 

 3 Figure 4

Acknowledgement: Dr. E. Kierlik and Dr. G. Tarjus are gratefully acknowledged for helpful discussions. The simulations were performed thanks to grant from the Institut de Développement des Ressources en Informatique Scientifique IDRIS-CNRS, Orsay, France.

B: Annealed and quenched complexities in the thermodynamic limit

The thermodynamic limit of the quenched and annealed complexities (Eqs. 3