E Kierlik

J Puibasset

G Tarjus

Effect of the reservoir size on gas adsorption in inhomogeneous porous media

Keywords:

We study the influence of the relative size of the reservoir on the adsorption isotherms of a fluid in disordered or inhomogeneous mesoporous solids. We consider both an atomistic model of a fluid in a simple, yet structured pore, whose adsorption isotherms are computed by molecular simulation, and a coarse-grained model for adsorption in a disordered mesoporous material, studied by a density functional approach in a local mean-field approximation. In both cases, the fluid inside the porous solid exchanges matter with a reservoir of gas that is at the same temperature and chemical potential and whose relative size can be varied, and the control parameter is the total number of molecules present in the porous sample and in the reservoir. Varying the relative sizes of the reservoir and the sample within experimental range may change the shape of the hysteretic isotherms, leading to a "reentrant" behavior compared to the grand-canonical isotherm when the latter displays a jump in density. We relate these phenomena to the organization of the metastable states that are accessible for the adsorbed fluid at a given chemical potential or density.

I. INTRODUCTION

Under confinement in disordered mesoporous materials, the characteristic time scale for relaxation of a fluid can become extremely long. As a result, and although not always appreciated, equilibrium is often not attained and, accordingly, bona fide thermodynamic transitions such as the liquid-gas transition in one-component fluids or macroscopic phase separation in mixtures, are unobservable. "Capillary condensation" in disordered solids is an out-of-equilibrium phenomenon, as illustrated by the irreversibility and hysteresis effects found in experiments. One typically observes a hysteresis loop that describes the isothermal evolution of the amount of fluid adsorbed in the porous solid as a function of the applied pressure, with branches that differ on adsorption (filling) and on desorption (draining). This hysteresis loop appears rate-independent, and its size and shape vary with temperature as well as with the characteristics of the solid (e.g. its porosity) or those of the solid-fluid interaction potential. Such a phenomenon is related to the existence of a large number of "metastable states" in which the system can be trapped on the experimental time scale; evolution from one metastable state to another then only occurs as a result of the action of the applied pressure (or equivalently, chemical potential) and it proceeds through a sequence of irreversible cooperative condensation (or evaporation) events, generically denoted as "avalanches". [1] The fact that the location and the shape of the hysteresis loop are reproducible in experiments indicates that the observation time is smaller than the time to reach the global equilibrium state, but is larger than local equilibration processes by which the system settles in one metastable state. As a result, the behavior of a fluid during filling or draining in disordered mesoporous materials can be rationalized by envisaging the evolution of the system in a free-energy landscape characterized by many local minima, i.e. metastable states. [2][3][4] The above picture of gas adsorption in disordered porous media brings in a strong analogy with the out-ofequilibrium response of systems driven by an external force in the presence of impurities or other types of quenched disorder. This is for example the case of magnetization cycles in ferromagnetic materials when a magnetic field is ramped up and down and of hysteretic martensitic transformations in alloys; [START_REF] Bertotti | The Science of Hysteresis[END_REF] in both examples, avalanches can be detected through some "crackling noise", [START_REF] Sethna | The Science of Hysteresis[END_REF] magnetic Barkhausen noise in the former, acoustic emission in the latter. In such driven disordered systems, one expects the occurrence of out-of-equilibrium phase transitions as one changes, on top of the driving force, some external parameters such as the temperature or the characteristics of the intrinsic disorder (e.g. the porosity in a porous solid). [START_REF] Sethna | [END_REF] The branches of the hysteresis loop, in particular the adsorption and the desorption isotherms, may then display jumps (discontinuities): indications for such behavior are for instance seen in the adsorption of helium in very light aerogels. [8][9][10] Such discontinuities and out-of-equilibrium phase transitions, however, are theoretically predicted on the basis of a grand-canonical set-up (for gas adsorption) in which the gas reservoir is infinite. In real experiments, the reservoir has a finite size which may not always be large enough for considering that the fluid inside the porous material is in a grand-canonical ensemble with fixed chemical potential. What should be expected in such situations ? The evolution of the system among metastable states may depend on the specific experimental set-up. As far as we know, there has been no systematic experimental study of the effect of changing the size of the gas reservoir (relative to that of the porous medium) and no estimate of the condition under which a grand-canonical situation is approximately reached. Answering these questions is the primary goal of the present paper.

The question of the dependence of the hysteretic response on the chosen control parameter has been a little more studied in other driven disordered systems. Two extreme cases have been considered: the response of an extensive quantity to a change in the (conjugate) "force" and the response of an intensive quantity (a force) to a change in the conjugate extensive quantity; the former is akin to the grand-canonical situation for a fluid in which the chemical potential (or actually, the pressure in the gas reservoir) is controlled and the latter to a canonical situation in which the number of adsorbed fluid molecules is controlled. Examples of such studies are found in the context of martensitic transformations in which either stress or strain is controlled [11] and in that of Barkhausen noise in which either magnetic field or, via some feed-back mechanism, magnetization is controlled. [START_REF] Bertotti | The Science of Hysteresis[END_REF] In all cases, the loop obtained with the extensive variable as control parameter appears as reentrant when compared to that obtained with the force as control parameter. (For a theoretical study, see Ref. [12])

One anticipates that the influence of the relative size of the reservoir[13] on the adsorption isotherms of a fluid in a disordered porous material is intimately connected to the organization of the metastable states in the adsorbeddensity/chemical-potential plane. Indeed, when the isotherms are smooth both on adsorption and on desorption, one may experimentally probe metastable states located inside the main hysteresis loop by studying ascending and descending "scanning curves", which involve partial filling or draining. [START_REF] Bertotti | The Science of Hysteresis[END_REF] These curves lead to a variety of hysteretic subloops that provide direct evidence for the presence of metastable states inside the main loop. [14][15][16] Actually, metastable states are expected everywhere inside the latter. The situation is quite different, however, when the adsorption or the desorption branch displays a jump (in a grand-canonical setting with a very large reservoir). A discontinuity prevents the realization of scanning curves in some portion of the isotherm. As a result, part of the interior of the main hysteresis loop is now inaccessible.

In the present work we address the above questions, namely the effect of the relative size of the gas reservoir on the adsorption isotherms of a fluid in a disordered or inhomogeous porous solid and the connection to the distribution of metastable states inside the hysteresis loop. We show that even when the grand-canonical isotherms display discontinuous jumps, there are branches of metastable states inside the main hysteresis loop and that many of these states can be reached by varying the relative size of the reservoir, the smaller the reservoir the larger the extent to which the branches of metastable states are probed. In particular, the extent is maximal when the reservoir is so small that the fluid inside the porous solid behaves as in a canonical situation of fixed number of adsorbed molecules. When compared to the loop obtained for an infinite reservoir (grand-canonical situation) and plotted as the amount adsorbed versus the chemical potential, the hysteresis loop for a finite reservoir size then appears as reentrant. On the contrary, when the isotherms are smooth (continuous), there is no influence of the size of the reservoir.

The rest of the paper is organized as follows:

In section II, we introduce the two models that have been investigated and give some details on the methods used to compute the adsorption isotherms. We have considered: (i) an atomistic model of a fluid in simple, yet structured pore, whose adsorption isotherms are computed by molecular simulation, and (ii) a coarse-grained model for adsorption in a disordered mesoporous material, studied by a density functional approach in a local mean-field approximation. In both cases, the fluid inside the porous solid exchanges matter with a reservoir of gas that is at the same temperature and chemical potential and whose relative size can be varied. The overall system composed of the sample plus the reservoir is taken in the canonical ensemble, with the total number of molecules as control parameter.

In section III, we present the results for the atomistic model. The simplicity of the system allows one to get a clear interpretation of the physical nature of the branches of metastable states, whose number is limited. The existence of metastable states is a direct consequence of the intrinsic inhomogeneity of the pore space. As the relative size of the gas reservoir is varied, we find that the exploration of the branches changes, the isotherm displaying larger jumps in the adsorbed density as the size increases.

Section IV is devoted to the coarse-grained model of a disordered porous material. The number of metastable states is now enormous, since it is likely to increase exponentially with the size of the sample. We find a drastic change of behavior between the "strong-disorder" regime for which the isotherm is continuous and the "weak-disorder" one for which the isotherm, here the adsorption isotherm on which we focus our study, displays a jump corresponding to a macroscopic avalanche. In the former regime, the isotherm does not display any dependence on the relative size of the reservoir whereas in the latter, it is reentrant for finite reservoir sizes, the more so as one decreases the reservoir size, and it shows increasing jumps in adsorbed density with increasing reservoir size.

Finally, we summarize our main results and give our conclusions in section V. In particular, we discuss the relevance of our study to experimental situations and we stress the important role of the intrinsic inhomogeneity induced by the solid matrix.

The principles of the method are the following. One considers a starting situation where N particles are placed inside two cells which are in thermal equilibrium with an infinite heat bath at temperature T . One of the cells represents a porous solid sample of volume V P and the other a reservoir of volume V R . The total volume is then equal to V = V R + V P . We allow mass exchange between the cells whose volumes are kept fixed so that at (exchange) equilibrium the chemical potential µ R of particles inside the reservoir is equal to the chemical potential µ P of particles inside the porous solid: µ = µ R = µ P . We measure the average amount of particles N P present in the porous solid. Varying the total number of particles N in small steps, one changes both the chemical potential of the cells, µ(N ), and the number of adsorbed particles, N P (N ). We then monitor the adsorption isotherm ρ P = N P /V P (for simplicity, we omit the subscript P in the fluid density adsorbed in the porous sample) as a function of µ. Introducing the ratio α = V P /V and the average reservoir density ρ R = N R /V R , where N R = N -N P , one has to satisfy for each µ the constraint (N/V ) = αρ P + (1 -α)ρ R with the condition µ = µ R (ρ R , T ) = µ P (ρ P , T ). Taking the limit α → 0 corresponds to the grand-canonical ensemble for the adsorbed fluid (with a controlled chemical potential) whereas taking the limit to 1 corresponds to the canonical ensemble for the adsorbed fluid (with a controlled number of adsorbed particles). Between these two extremes, one has a mixed ensemble, with a reservoir of variable size compared to the sample, somewhat reminiscent of the mesoscopic canonical ensemble considered by Neimark and coworkers. [17] Low values of α correspond to a large reservoir with respect to the porous material size : experimental values are typically lower than 10 -2 which corresponds to a volume of the reservoir which is few hundreds the volume of the sample. More details about the algorithms used in our two models are given below.

B. Atomistic model of a fluid in a single structured pore An atomistic model of confined fluid was considered because it allows a realistic description of the various (metastable) states the system can adopt and, possibly, of the transition mechanisms at a molecular scale. The chosen system is a simple atomic fluid (Lennard-Jones like) confined in a mesoporous substrate. (The model more specifically corresponds to Argon adsorbed in nanoporous solid carbon dioxide.) [18,19] A fully realistic model of the substrate should take into account the surface roughness, pore morphology (pore size distribution), and interconnections between pores. The number of metastable states would however be too large for a systematic study. The model was therefore designed to exhibit only a few metastable states.

The main potential sources for generating metastable states are nanometer-scale heterogeneities, due to pore-size distribution, and variations in surface chemistry. In both cases, the prominent effect on the adsorbed fluid comes from the modulation of the effective fluid/wall interaction. [20,21] We then chose to investigate a cylindrical pore containing few domains of variable fluid/wall interaction. [22] The diameter of the nanopore is 8 σ f f (σ f f is the fluidfluid Lennard-Jones diameter). In this way, the smooth-wall approximation can be applied: the external potential seen by a fluid particle in the pore is calculated by integrating the fluid-wall (6 -12) Lennard-Jones potential (ǫ sf = 1.277 ǫ f f , σ sf = 1.094 σ f f ) over a uniform distribution of substrate interacting sites of density 0.8265 σ 3 f f . The calculated reduced external potential Ψ * cyl (r) = Ψ cyl (r)/ǫ f f in this perfectly cylindrical pore is given in Fig. 1. The heterogeneity is introduced by modulating this external potential along the axial direction z:

Ψ * pore (r, z) = 1 + a(z) sin 4πz L Ψ * cyl (r) (1)
where a = 0.3 for z < 0 and a = 0.2 for z > 0, and L = 24σ f f is the simulation box length. This external potential exhibits four domains of spatial extension of a few molecular diameters along the axial direction (see Fig. 2). Periodic boundary conditions are applied along the axial direction z. The interactions are truncated at half the simulation box size (minimal image convention). The fluid adsorption properties are calculated by Monte-Carlo simulations. Thermalization of the adsorbed fluid is performed by particle displacement trials. Chemical equilibration between the adsorbed fluid and the reservoir is done by particle exchange trials. The acceptance probabilities are given by the Metropolis algorithm. The fluid in the reservoir is assumed to be ideal, which is a good approximation since the simulations are performed well below the critical point. As a consequence, the gas in the reservoir does not need to be treated explicitly, which speeds up the calculations. In the limit of infinite reservoir size, one recovers the usual Grand Canonical Monte Carlo algorithm (GCMC). The adsorption/desorption curves are determined as follows. The initial configuration is the empty pore. Few molecules are then added in one step to the system (pore+reservoir). In order to mimic the experimental situation,

r / σ σ σ σ f-f 0 1 2 3 4 Ψ Ψ Ψ Ψ cyl / / / / ε ε ε ε f-f -4 -2 0 2
FIG. 1: External potential Ψ * cyl (r) for the perfectly cylindrical pore as a function of radial distance r. the molecules are actually added to the reservoir cell, which results in an increase in µ R , and, accordingly, to a lack of equilibrium between the pore and its reservoir. The whole system is then relaxed for a while. Molecules exchange between the cells, which produces a small flow from the reservoir to the pore in the first steps. This flow comes to zero when stationarity is reached (µ R = µ P ). We emphasize that this matter flow does not describe correctly the mass transport, since the Monte Carlo algorithm allows transfer of matter everywhere in the porous substrate. 10 6 Monte Carlo trials per particle are then performed to acquire statistics for computing averaged quantities. This gives the first point of the isotherm. The subsequent points are obtained according to the same procedure, i.e. a small increase in the number of molecules in the reservoir, followed by a relaxation run, and, finally, by a long run for acquisition.

Ψ Ψ Ψ Ψ pore / / / / ε ε ε ε f-f r / σ f-f z / σ f-f 1
After complete adsorption has been achieved, the desorption isotherm can be calculated by decreasing the amount of particles in the reservoir step by step, as in a real experiment. This algorithm allows us to obtain the scanning curves embedded within the main hysteresis loop. Note that for the case of an infinite reservoir, it is the chemical potential which is increased stepwise instead of the amount of particles in the reservoir, which is then infinite.

C. Coarse-grained model for adsorption in a disordered porous material

As discussed in previous papers, [2,3,23] our approach to fluid adsorption in disordered mesoporous materials is based on a coarse-grained lattice-gas description which incorporates the essential physical ingredients of the solid-fluid system. We consider a three-dimensional BCC-lattice of linear size L in which each of the N SP = 2L 3 sites may be occupied by a fluid or by a solid particle (L is measured in units of the lattice spacing a and we set a ≡ 1; for a lattice, the volume V P is then simply equal to the number of sites N SP times the volume of the elementary cell a 3 ≡ 1, so that we shall speak of volume and number of sites indistinctly). Multiple occupancy of a site is forbidden and only nearest-neighbor (nn) interactions are taken into account. The fluid particles can equilibrate, as explained below, with a finite-size reservoir at fixed temperature whereas the solid particles are "quenched" and distributed according to a specific choice of the porous material structure. For the sake of simplicity, we study a random matrix with a porosity φ = 75%. The relevant correlation length of the solid is around one lattice spacing (it is a purely random matrix) so that even the smallest systems studied (with a linear size L = 25) can correctly describe the collective effects occurring inside the matrix on long length scale (such as a sharp condensation event in the whole pore space).

The starting point of our theoretical analysis is the following expression of the free-energy functional in the local mean-field approximation (hereafter this analysis will be called LMFT) :

F P [{ρ i }] = k B T i [ρ i ln ρ i + (η i -ρ i ) ln(η i -ρ i )] -w f f <ij> ρ i ρ j -w sf <ij> [ρ i (1 -η j ) + ρ j (1 -η i )] (2)
where ρ i is the thermally averaged fluid density at site i and η i = 0, 1 is a quenched variable describing the occupation of the lattice by solid particles (η i = 0 if the site i is occupied by the solid and η i = 1 otherwise); w f f and w sf denote the fluid-fluid and fluid-solid attractive interactions, respectively, and the double summations run over all distinct pairs of nn sites. We first start with the grand canonical situation (α = 0) where fluid particles can equilibrate with an infinite reservoir that fixes the chemical potential µ. Minimizing the grand-potential functional Ω P [{ρ i }] = F P [{ρ i }] -µ i ρ i with respect to {ρ i } at fixed T and µ for a given realization of the solid yields a set of coupled equations,

ρ i = η i 1 + exp[-β(µ + w f f j/i ρ j + w sf j/i (1 -η j ))] , (3)
where the sums run over the c = 8 nearest neighbors of site i. By using a simple iterative method to solve these equations, one finds solutions that are only minima of the grand potential surface, i.e. metastable states. For a given realization of the solid, the adsorption isotherm is obtained by increasing the chemical potential in small steps δµ. At each subsequent µ, the converged solution at µ -δµ is used to start the iterations. What happens in the mixed situation (α = 0)? The fluid particles can equilibrate between the solid sample and a finite-reservoir so that their total number N is fixed. The isotherm is then obtained by increasing N in small steps δN . In this case, the system tries to minimize its total Helmholtz free-energy

F T [{ρ i }, ρ R , T ] = F P [{ρ i }, T ] + F R [ρ R , T ],
where F R is the Helmholtz free-energy of the gas reservoir

(F R [ρ R , T ] = V R {k B T [ρ R ln ρ R + (1 -ρ R ) ln(1 -ρ R )] - w f f cρ 2
R /2} with V R the volume of the reservoir, i.e., again up to a 3 ≡ 1, the number of sites of the reservoir, while satisfying the global constraint

N = N P + N R = i ρ i + V R ρ R .
This can be solved in a natural way by the method of Lagrange multipliers. We consider the function

ΩT [{ρ i }, ρ R , λ, T ] = F T [{ρ i }, ρ R , T ] + λ{N - i ρ i -V R ρ R }, (4)
where λ is a Lagrange multiplier that has the meaning of the chemical potential coupled to the densities. Minimizing F T with the constraint on densities amounts to solve simultaneously the coupled equations ∂ΩT ∂ρi = 0, ∂ΩT ∂ρR = 0 and

∂ΩT ∂λ = 0 or equivalently k B T ln[ ρi ηi-ρi ] -λ -w f f j/i ρ j -w sf j/i (1 -η j ) = 0, 1 ≤ i ≤ N SP (5) k B T ln[ ρR 1-ρR ] -λ -w f f ρ R = 0, (6)
N -i ρ i -V R ρ R = 0. ( 7
)
One has then to define an iterative scheme that specifies how the system goes from one converged solution to another as the total number of particles is slowly changed. The details are given in the appendix. We are searching for configurations for the local densities ρ i that are local extrema of the grand-potential functional Ω P for the special value of the chemical potential µ = λ. However, it is not fully guaranteed that the above algorithm necessarily converges to a local minimum, nor even to an extremum, since the constraint could in principle stabilize an unstable state. Therefore, we regularly ascertained that configurations obtained in the mixed ensemble are stable in the grand-canonical ensemble, i.e. that they indeed are metastable states, by starting grand-canonical calculations with the converged solution at chemical potential λ. In addition, it must be emphasized that the global constraint is not satisfied until convergence is reached, in the spirit of the Lagrange method. It is therefore doubtful that one can attribute any physical meaning to the intermediate stages of the iterative process. Recall also that one does not take into account mass transport: local densities can change everywhere and instantaneously inside the porous sample.

Moreover, we focus on the adsorption (filling) process. As shown in previous papers, [3,24] fluid desorption may crucially depend on the presence of an external surface for the porous solid: the system then includes a real interface between the adsorbed fluid and the external vapor, and, during desorption, the vapor domain may penetrate and drain the solid from the outer surface (the so-called percolation and depinning transitions discussed in Ref. [24]. These mechanisms are no longer bulk phenomenon and their, a priori more subtle, analysis will be carried out elsewhere. We have therefore used periodic boundary conditions for the sample and the reservoir (separately). The ratio y = w sf /w f f has been fixed to 0.8 so that the adsorption isotherms exhibit the interesting range of phenomena as the temperature changes (see Ref. [3]).

III. RESULTS ON RESERVOIR-SIZE DEPENDENCE FOR THE ATOMISTIC MODEL A. Grand-canonical isotherms

We first focus on the adsorption/desorption isotherms obtained in the limit of infinite reservoir size (α = 0), i.e. the GCMC data. The results are given in Fig. 3 for various temperatures. As can be seen, the highest temperature isotherm (reduced temperature T * = k B T /ǫ f f = 1.0) is perfectly reversible, i.e. the adsorption and desorption curves are superimposed. On the other hand, for lower temperatures, the adsorption and desorption curves differ and exhibit hysteresis with vertical steps. We emphasize that the vertical lines are guides to the eye and do not correspond to GCMC data. The isotherm jumps in one step, and cannot be stabilized in between.

Quite noticeably, the simulation points may be grouped into branches. For instance, for a reduced temperature of 0.9, three branches are present: the gas-like branch (with a fluid layer adsorbed at the wall, the rest of the pore being filled with gas), which exists down to very low µ, the liquid-like branch (pore filled by dense fluid) at high µ, and a branch of intermediate density which exists only for a limited range of µ. We have found that adsorption and desorption are reversible along any given branch. These branches are associated to deep local minima in the free energy landscape describing the system (the grand potential). The energetic barriers separating these minima are generally large compared to the thermal fluctuations sampled by the Monte-Carlo algorithm: the system then remains trapped in these local minima, which explains the reversibility of adsorption/desorption along the branches. However, for some particular values of µ the barriers become sufficiently small for allowing thermal fluctuations to make the system jump into an adjacent local minimum. These values define the limit of stability of the branches (in the grand canonical ensemble). The vertical lines in Fig. 3 indicate the new local minimum (branch) reached by the system. Note that the µ-range of the various branches overlap, which means that for some µ, the system may be stabilized in GCMC at various degrees of pore filling ρ P . This is the main origin of hysteresis. As can be seen, the lower the temperature, the larger the µ-range of existence of the branches, and the wider the hysteresis. The number of branches increases for decreasing temperature, at least down to T * = 0.8 (Fig. 3). We shall show later that this is actually still true down to T * = 0.6 if one takes into account metastable states that are hidden in the GCMC. For T * = 0.9, the three branches are visited during both adsorption and desorption. They belong to the main adsorption and desorption curves. For T * = 0.8, five metastable states belong to the main adsorption curve, and three to the main desorption curve. The branch of intermediate density reached during desorption (second branch, ρ * P = ρ P σ 3 f f ∼ 0.4) looks like a continuation of the third branch visited during adsorption (ρ * P ∼ 0.4). This point was checked by showing the reversibility of the complete branch by GCMC (by increasing and decreasing µ). The total number of metastable states is then five, and their range of existence was determined in a systematic way by µ-variations. The results are shown in Fig. 4. Here again, lines are guides to the eye, and the vertical lines show the new metastable state reached by the system beyond the stability limit. The case T * = 0.6 deserves special attention. The procedure that consists in following simple ascending and descending µ-paths does not allow to reach the two branches with reduced densities around 0.3 and 0.45. These states were obtained from the corresponding states at higher temperature T * = 0.8 by following a particular (µ,T ) path along which their stability is preserved. This path essentially consists of slowly decreasing the temperature and increasing the chemical potential in order to keep the average number of particles roughly constant. When T * = 0.6 has been reached, the µ-range stability can be determined for both states (see Fig. 4). It is possible that more metastable states actually exist, but the procedures previously explained did not allow finding more than five branches in our simple system.

B. Metastable states

Figure 5 displays an enlargement of the lowest temperature (T * = 0.6) GCMC results, without the vertical lines. As previously mentioned, the simulation points group into reversible branches (solid lines are guides to the eye) which correspond to local minima in the free-energy landscape. We now give a molecular-level description of these local minima. The lowest density branch corresponds to gas-like fluid filling the pore with adsorbed molecules at the wall. The highest density branch corresponds to liquidlike fluid saturating the pore. The three other states are stabilized by the chemical corrugation. Visual inspection of molecular configurations (see lower panel of Fig. 5) shows that the three branches correspond respectively to: one liquid-like bridge in the most attractive domain of the pore; two liquid-like bridges in the two attractive regions of the pore; one single gas-like "bubble" in the least attractive domain of the pore. Lines are guides to the eye indicating the path followed by the system beyond the stability limits of each metastable state. The hysteresis loops found in Fig. 3 are the envelopes of all metastable branches.

C. Mixed-ensemble results

We now focus on the adsorption/desorption results obtained for a finite-size reservoir. Two ratios V R /V P = 500 and 2000 have been considered (α ≃ 2.×10 -3 and 5.×10 -4 ). The procedure previously used to obtain all metastable states will be applied here in its simplest form, i.e. isothermal µ-paths, because it corresponds to what is actually feasible for realistic mesoporous substrates and in real experiments (one essentially measures the main adsorption/desorption hysteresis and scanning curves). The results for the main hysteresis loops are given in Fig. 6, where the grand canonical case (α = 0) has also been plotted for comparison. In all cases the symbols are the simulation data, and the color (thick) lines are guides to the eye that connect the simulation points in the same order as they were obtained by slowly increasing (adsorption) or decreasing (desorption) the total number of particles in the system. The thin (black) lines are the five branches previously obtained by extensive GCMC study (Fig. 5). In the limit of infinite reservoir size (α = 0), starting from the empty system and gradually increasing µ results in the continuous filling of the system up to µ * = µ/ǫ f f = -9.43 where the system jumps into the second local minimum. Increasing further the chemical potential makes the system jump directly to the fifth local minimum. Upon desorption, the system remains in the saturated state for a while (hysteresis) and finally jumps (for µ * = -9.94) directly into the first local minimum (gas-like branch) without visiting any intermediate state. The vertical lines represent the constraint imposed by the reservoir (µ * = cte) but not necessarily the path actually followed by the system during the transition (transient phenomenon are not properly described by the Monte-Carlo algorithm).

In the cases of finite reservoir, it is the total number of particles in the system plus its reservoir which is increased by small steps as in a real experiment. As can be seen from Fig. 6, the amount of adsorbed fluid initially increases along the gas-like branch in a continous manner until it reaches its stability limit. (Note that this limit occurs at slightly larger µ for smaller reservoir: µ * = -9.43 for α = 0, µ * = -9.41 for α ≃ 5. × 10 -4 and µ * = -9.36 for α ≃ 2. × 10 -3 .) At this stability limit, a small addition of extra molecules in the reservoir destabilizes the adsorbed fluid which then jumps into the second state. During the transition, particles leave the reservoir and adsorb into the pore in the most attractive region to form a liquid-like bridge (second metastable state), which results in a decrease of the chemical potential. Here again, our algorithm is not meant to describe this transition. The path actually followed by the system is not known. The (curved) line shown in the figure corresponds to the constraint of conservation of the total number of particles in the system plus its reservoir. After equilibration, the system reaches the second branch. Further increase in the total number of particles makes the system follow this branch until it reaches a new stability limit. (As previously, this limit slightly increases for a finite reservoir: µ * = -9.27 for α = 0, µ * = -9.23 for α ≃ 5. × 10 -4 and µ * = -9.24 for α ≃ 2. × 10 -3 .) The system then jumps onto the fifth branch for α = 0, onto the fourth branch for α ≃ 5. × 10 -4 , and onto the third branch for α ≃ 2. × 10 -3 .

The above results illustrate one of the main effects of the reservoir: the system does not necessarily visit the same metastable states as one changes the size of its reservoir. From Fig. 6, this effect may be tentatively interpreted as follows: the reservoir imposes a constraint on the relation between the chemical potential and the amount adsorbed, materialized by the vertical or inclined curves. After a branch stability limit has been reached, the system approximately follows this line of constraint during the transition, and finally meets a new branch at a location that depends on the ratio between the reservoir and the system size. To illustrate this point, starting from the highest µ GCMC point of the second branch (limit of stability), we have determined the new metastable state reached by the system after a small increase in the total number of molecules for various reservoir sizes. The results are given in Fig. 7. This picture is however somewhat oversimplified and does not apply in all situations since it does not take into account the complexity of the underlying free-energy (grand-potential) landscape that determines which path is actually followed by the system (being in configurational space, the free energy landscape is an object of very high dimension). At a molecular level, this means that complex nucleation and fluid adsorption processes determine the path actually followed by the system. One of the consequences is that the system does not necessarily stop on the first branch crossed by the constraint line. This is illustrated in our simple system for α ≃ 2. × 10 -3 during desorption from the fourth branch: as can be seen, the constraint line crosses the third metastable branch, but the system avoids this state and actually reaches the second branch. In the grand-potential landscape picture, this means that, starting from the fourth local minimum, and following the steepest slope, the system has not crossed the basin of attraction associated to the third minimum and has followed its way until being finally trapped by the basin associated to the second local minimum. In the molecular-level description, the fourth branch corresponds to a large liquid domain with a single gas-like "bubble" in the least attractive region. Upon desorption, the system would have to nucleate a bubble in between the attractive regions in order to reach the third branch corresponding to two liquid-like bridges in the attractive regions (see Fig. 8, upper path). This nucleation barrier is probably too high, and it is more favorable for the system to desorb by recession of the liquid bridge in the less attractive region until it is left with one single liquid bridge in the most attractive region corresponding to the second state (see Fig. 8, lower path).

As for the grand-canonical case, it is possible to perform a systematic search for all possible metastable states in the system in contact with a finite reservoir. The five branches already described are found in all cases. Figure 9 shows the superimposition of all simulation points (stable over long runs) obtained for the various ratios between the reservoir and system sizes already presented. As can be seen, the points group onto branches corresponding to the five metastable states found in the system. It is however reminded that the limits of stability of the branches slightly vary with the reservoir size. This may be related to the fact that transitions proceed via energetic barriers and that the amplitude of the fluctuations allowed by the reservoir depends on its size. To summarize the results, reducing the size of the reservoir allows one to explore metastable states upon adsorption/desorption that are not accessible in the grand-canonical situation. It is also found that the number of metastable states visited by the fluid increases with decreasing reservoir size.

IV. RESULTS FOR THE COARSE-GRAINED MODEL A. Grand-canonical isotherms

Typical results for adsorption isotherms obtained at various temperatures in the grand-canonical ensemble are shown in Fig. 10. When the temperature decreases, the isotherm shape changes from smooth to steep. As discussed in previous papers, [24] this corresponds to a true out-of-equilibrium phase transition, the so-called "avalanche transition", [START_REF] Sethna | [END_REF] with the sudden appearance of a macroscopic, connected liquid domain in the whole porous sample. A previous extensive scaling study [24] showed that it exists a critical value of the temperature for which the isotherm changes from continuous to discontinuous in the thermodynamic limit. At higher temperature, isotherms look gradual but at reduced temperatures T * = 1.40, they consist of little steps of varying sizes (see insert in Fig. 11).

The above results are related to the characteristics of the grand-potential landscape. At low temperatures, this multidimensional landscape (recall that it is a function of the local fluid densities, with here N SP ∼ 10 4 -10 6 ) is characterized by a large number of local minima, the metastable states in which the system may be trapped. Since thermally activated processes are neglected in the LMFT (an approximation that, as already discussed, finds its justification in the experimental reproducibility of the adsorption isotherms on the time scale of most experiments), the evolution of the system is only due to a variation of the external driving (here the chemical potential). As µ varies, the system either follows the minimum in which it is trapped as this minimum deforms gradually (the flat portions in the insert of Fig. 11), or it falls instantaneously into another minimum when the former reaches its stability limit. This later move is a discontinuous and irreversible process, an avalanche, which is at the origin of the history-dependent behavior of the system, e.g. the hysteresis. The avalanche corresponds to a collective condensation event inside the porous sample that manifests itself by a jump in the adsorption isotherm [1]. The size of the avalanche may be macroscopic as for T * = 0.8, and the adsorption isotherm is discontinuous in the thermodynamic limit; or it may be microscopic as for T * = 1.40, and the adsorption isotherm remains smooth in the thermodynamic limit. (More calculations should be necessary to conclude on what happens at T * = 1.10.) The hysteresis loop encloses all the metastable states of the system: this is illustrated in Fig. 11 by the scanning curves obtained by performing incomplete filling of the matrix and then decreasing the chemical potential to drain the adsorbed fluid. Of course, the number of metastable states is very large, expected to be exponential in the system size as in disordered magnetic systems (see e.g. [25]), and only few of them are revealed with this simple procedure.

B. Canonical and mixed-ensemble isotherms

What happens when the porous sample is coupled to a finite reservoir? First, we begin with a vanishingly small reservoir, which is the canonical situation where one controls the number of adsorbed particles. The main results are summarized in Fig. 12 where we compare the adsorption isotherms obtained with the "µ-driven" (grand-canonical) and "ρ P -driven" (canonical) procedures. The behavior is quite different in the low and the high temperature regimes, respectively characterized by the absence and the presence of a macroscopic µ-driven avalanche. At the highest temperature studied (T * = 1.80, not shown in the figure), there exists only one stable state, the equilibrium one, and controlling either the adsorbed density or the chemical potential yields the same result. This is also true for all temperatures at very low adsorbed densities (not shown in the figure). At lower temperatures, avalanches appear in the grand-canonical isotherm. Collective localized events, which we define as "avalanches" irrespective of the control variable, also show up in the canonical isotherm in the form of small jumps in the chemical potential toward a lower value: with the disappearance of the initial minimum, the system has to find another metastable state with the required adsorbed density; as this state cannot be found at a higher value of the chemical potential since the grandcanonical adsorption/desorption isotherms have been shown to represent the extremal curves that encompass all the metastable states of the system, [3] the chemical potential must decrease. (Note that in a grand-canonical setting "avalanches" appear as jumps in the adsorbed density ρ P at constant chemical potential µ whereas in a canonical one, they appear as jumps in µ at constant ρ P .) When a metastable state is found and fluid is further added, the system smoothly follows this state (see the quasi-linear portions in Fig. 13) until it reaches the corresponding stability limit. Then, there is a new jump in the chemical potential to a smaller value, and the evolution proceeds in this way until the porous sample is completely filled with liquid. The contrast between µ-driven and ρ P -driven isotherms that is illustrated in Fig. 12 is very reminiscent of what was found in ferromagnetic systems when comparing magnetizationdriven and magnetic-field-driven protocols. [12] As in the latter case, the canonical (ρ P -driven) isotherms are closely related to the distribution of metastable states inside the (grand-canonical) hysteresis loops.

We now discuss in more detail the canonical (ρ P -driven) isotherms. At T * = 1.40 with a sample of linear size L = 100, the jumps in µ(ρ P ) are very small and one needs to zoom in on the isotherm to see small differences between the canonical and grand-canonical protocols, as illustrated in Fig. 13. This smallness is related to the presence of many metastable states in the close vicinity of the grand-canonical isotherm. Since both the grand-canonical and the canonical avalanches remain of microscopic size, we expect that the small jumps in either ρ P or µ become infinitesimally small in the thermodynamic limit. This is indeed what is observed in Fig. 13 where isotherms obtained for different system sizes are compared. Therefore, we predict that the canonical and grand-canonical curves should become identical in the thermodynamic limit.

A similar behavior is observed at T * = 1.10, as illustrated in Fig. 14. However, the jumps in ρ P are larger than at T * = 1.40 and are clearly associated with a much more fluctuating chemical potential µ(ρ P ): each jump in ρ P generates a re-entrant behavior in µ(ρ P ). It is difficult to conclude on how the isotherm will evolve in the thermodynamic limit: as the size of the system increases (see Fig. 14), there are more and more intermediate points in the steepest part of the grand-canonical isotherm and, at the same time, the canonical isotherm is less and less fluctuating. The temperature T * = 1.10 is probably barely above the critical temperature of the (grand-canonical) avalanche transition and we expect that both the canonical and the grand-canonical isotherms would become continuous and would coincide in the thermodynamic limit.

On the other hand, T * = 0.80 is undoubtedly below the critical temperature of the avalanche transition: all samples studied, whatever their size, display a large jump in their grand-canonical isotherms (see Fig. 15). (In a previous work on a random matrix, [24] we concluded through an extensive finite-size study that a macroscopic jump exists FIG. 13: A close-up of the canonical (symbols, α = 1) and grand-canonical (lines, α = 0) adsorption isotherms obtained by LMFT at T * = 1.40 for different linear sizes of the system (colors on line).

in the thermodynamic limit at T * = 0.8 for y = 0.9; since the value y = 0.8 used in this paper corresponds to a weaker disorder than y = 0.9, this is necessarily true here as well.) As shown in Fig. 15, the canonical curves show a pronounced reentrant behavior, with a large difference with the grand-canonical ones: a whole region void of metastable states then appears. There are always finite-size effects, which affect both the canonical and the grandcanonical isotherms. In the latter case, the rare events that trigger the macroscopic avalanches are very sensitive on details about the structure of the matrix, and the corresponding chemical potentials vary strongly with the system size at small L. We have nonetheless checked by performing grand-canonical isotherms on very large systems of linear size L = 200 that the position of the jump stays very close to that found in the isotherm for L = 100 (displayed in Fig. 15), which indicates that one is close to the thermodynamic limit. For the canonical isotherms, the chemical potentials fluctuate less and less as the size of the system grows (especially at the end of the adsorption process, when the porous sample is nearly filled with liquid) and the overall location of the isotherm in the ρ P -µ plane does not shift significantly with L.

To discuss the fluctuations of the chemical potential with the density in the canonical protocol, we have performed finite-size studies of the first moments of their distribution, after averaging over disorder (i.e. porous sample) realizations. More precisely, we compute the mean µ(ρ P ) (hereafter an overbar denotes the average over disorder) and the variance ∆ µ (ρ P ) = µ(ρ P ) 2 -µ(ρ P )

2 for three sizes of the system. As explained in Ref. [12], the standard argument concerning self-averaging quantities, whose value in a macroscopic sample is equal to the average over all disorder realizations, cannot be applied: as a consequence, one cannot be sure that µ(ρ P ) is a self-averaging quantity. To delve more into its behavior, it is interesting to investigate the dependence on system size of the standard deviation σ µ (ρ P ) = ∆ µ (ρ P ). The result is shown in Fig. 16, where σ µ (ρ P ) is plotted as a function of the adsorbed density (computed in bins of width 0.01). Except at low and high density, the standard deviation remains nearly constant with ρ P (typically, between ρ P = 0.2 and ρ P = 0.6). It decreases with L but quite slowly: we find that in the region between ρ P = 0.2 and ρ P = 0.6, ∆ µ (L) = σ µ (L) 2 ∼ L -γ with a finite-size scaling exponent γ around 1.4, i.e. significantly smaller than 3. This indicates that µ(ρ P ) is self-averaging, but only weakly so. As in Ref. [12], we have checked that the corresponding histograms are roughly Gaussian (see Fig. 16). Although the situation with respect to self-averaging remains somewhat unclear as far as the whole isotherm is concerned, [12] it is nonetheless instructive to study the mean µ(ρ P ) as a function of system size. This is displayed in Fig. 17, where one can see that the µ(ρ P )-isotherms vary significantly with the size of the system, becoming steeper as L increases. We cannot conclude whether or not the isotherm becomes vertical in the thermodynamic limit, but the data suggests that it will be distinct in this limit from the grand-canonical isotherm: the µ(ρ P )-isotherms for all studied system sizes intersect near the point µ * = -4.47 and ρ P = 0.3, at a significantly lower chemical potential than those found for the grand-canonical jumps (µ * ∼ -4.43 -4.44) for L = 200 (see Fig. 17). We now comment on the two types of behavior found along the ρ P -driven adsorption isotherm. For increasing adsorbed density ρ P , the chemical potential may either increase continuously or decrease by discontinuous jumps. In the former case, there are only slight swellings of the liquid domain, with the liquid-gas interfaces retaining the same shapes locally and all the ρ i 's increasing. In the latter case, we find that each jump corresponds to a single condensation event, a ρ P -driven "avalanche", in which all sites of a compact region become liquid. This is illustrated in Fig. 18. Considering the two consecutive configurations with average densities ρ P and ρ P + ∆ρ P (before and after the jump), sites are considered as turning liquid (respectively, gas) when the variation of the local fluid density is larger than 0.3 (respectively, lower than -0.3) and are shown in red (respectively, grey) in the figure. The size of the condensed region varies for each jump along the isotherm, but the number of particles involved in the local condensation can go much beyond the controlled increment V P ∆ρ P : in Fig. 18, the condensation event (in red) corresponds to a local increase in the number of particles that is 110 times larger than the initial increment. Therefore, others regions of the porous sample must be drained: as shown in the figure, the most important evaporation events (in grey) remain smaller and less compact than the condensation event and appear at random throughout the material. (Whereas the definition of "turning liquid" is well characterized since the histogram of positive variations of the local density is always peaked at a value close to 1, this is much less so for its counterpart "turning gas": the histogram of negative variations is continuously decreasing between 0 and -1, so that the decrease of the average density that balances the condensation event mostly comes from small but widespread local variations.) This suggests the following interpretation: adding a small amount of fluid in the sample may trigger a condensation "avalanche" whose size is larger than the added fluid amount; a decrease of the chemical potential allows a slight draining of the rest of the sample and provides the required amount of fluid to "feed" the avalanche.

Finally, Fig. 19 presents our results in the mixed-ensemble for different sizes of the reservoir. As the size increases, the µ(ρ P )-isotherms exhibit flat segments with decreasing length and jumps in the µ-ρ P plane that are more and more pronounced. This behavior is similar, albeit in a more complex systems and in the presence of a much larger number of metastable states, to that found in the atomistic model and illustrated in Fig. 6. When the volume of the reservoir becomes 2000 times greater than the volume of the porous material, the isotherm displays a unique jump between the grand-canonical low-density branch and the grand-canonical high-density branch, but it still comes with a decrease of the chemical potential (lower right panel of Fig. 19). For larger reservoir sizes, the jump gets closer and closer to the grand-canonical isotherm. Convergence however is rather slow. On the other hand, for the smallest reservoir sizes (a reservoir volume 500 times the volume of the porous material or less: see the two upper panels of Fig. 19), the isotherms tend to superimpose on the canonical isotherm. It can be however noticed that the explored metastable states are not necessarily identical. As explained in section III, one indeed has to be very careful with the representation of states in the µ-ρ P plane which is only a projection from the highly multi-dimensional configurational space: it is possible to cross branches as the mixed-ensemble constraint allows the system to explore other areas of the underlying free-energy landscape.

FIG. 19: Adsorption isotherms at T * = 0.80 predicted by LMFT for different relative sizes of the reservoir for a sample of linear size L = 50. In the four panels, the canonical isotherm is shown in light grey and the grand-canonical isotherm as a dashed curve. Top left: α = 10 -2 , top right: α = 2. × 10 -3 , bottom left: α = 10 -3 , bottom right: α = 5. × 10 -4 , where α = VP /(VP + VR). Note that both the slope and the length of the inclined segments that represent jumps in the µ-ρP plane (the lines are therefore a guide to the eye) increase as α decreases: the slope is directly related to α (compare with Fig. 7).

V. CONCLUSION

In this paper, we have presented a study of the out-of-equilibrium, hysteretic response of a fluid adsorbed in inhomogeneous porous materials when one couples the porous sample with a finite-size reservoir and controls the total number of particles. Varying the relative size of the sample and the reservoir allows one to interpolate between a canonical situation with a controled adsorbed density and a grand-canonical situation with a controled chemical potential. We have considered both an atomistic model of a fluid in simple, yet structured pore and a coarse-grained model for adsorption in a disordered mesoporous material. The adsorption isotherms have been computed by molecular Monte Carlo simulation in the former case and by a density functional approach in a local mean-field approximation in the latter. The Monte Carlo simulations give us a clear picture of what occurs at the molecular scale in a small and weakly disordered system whereas the density functional approach provides insights at the mesoscopic scale for a fully disordered system. In both cases, we have found metastable states that appear as branches of finite extent in the µ -ρ P plane. The metastable states correspond to inhomogeneous configurations of the fluid and the number of branches increases rapidly with the complexity of the material. It is worth stressing two points: first, these states are not unstable, i.e. not purely stabilized by the constraint on the total number of particles [26]; second, metastability is induced by the intrinsic inhomogeneity of the solid. the picture of metastability in such systems is quite different than that encountered in mean-field-like descriptions of homogeneous (bulk) systems undergoing first-order transitions.

We have shown that the way the system evolves between these metastable states may depend on the protocol, controlled here by the relative size of the reservoir. In particular, our results suggest that a discontinuity in the grand-canonical adsorption isotherm (an out-of-equilibrium "avalanche transition") is associated with the absence of metastable states in a whole region of the µ -ρ P plane and that the corresponding canonical adsorption isotherm (and, more generally, isotherms performed with a small enough reservoir size) differs from the grand-canonical one and displays a reentrance, even in the thermodynamic limit.

What is the relevance of our results for experimental set-ups ? For illustration, we consider two examples. In order to determine the adsorption isotherm with the volumetric method using for instance nitrogen, known amounts of fluid are admitted stepwise in the sample cell. The amount of adsorbed fluid is the difference between the admitted gas and the amount of gas that fills the "dead volume", i.e. the free space in the sample cell including connections: this is equivalent to the procedure discussed in this paper, with the dead volume playing the role of a finite-size reservoir for the sample. The amount of gas in this reservoir is calculated from the fluid equation of state and from measurements of the pressure, the temperature and the dead volume. In an experiment with a sample size of 20 mm 3 , a cell of a few cm 3 could be used. In such a situation, the ratio of volumes (porous medium versus total) is around 10 -2 and the ratio of adsorbed amount on the total amount around 1. On the other hand, to measure Helium adsorption in aerogels, B. Cross et al [27] use a different experimental set-up in which the adsorbed amount is controlled through the temperature of a helium gas reservoir connected to the experimental cell. The total amount of Helium is fixed, but varying the temperature of the reservoir transfers atoms from the reservoir to the cell, or conversely. In the recent experiments, [10] with temperatures ranging from 4K to 5K, the adsorbed amount at saturation corresponds to around 20% of the total amount.

In both above examples, the ratio between the adsorbed amount and the total amount of fluid is fully in the range of the parameters of our study. However, in most experiments in disordered porous materials the adsorption isotherms display hysteresis but are smooth and continuous: in consequence, as predicted by the present study, no effects of the reservoir size are expected in this case. We have shown here that for such effects to be observable, the temperature should be low enough for the grand-canonical adsorption isotherm to exhibit a true discontinuity and not only a very steep variation. It is unlikely that this true discontinuity could be observed in disordered porous materials such as Vycor or xerogels: the disorder and the confinement are too strong and prevent the appearance of the avalanche transition at temperatures higher than the triple point. The case of helium adsorption in aerogels is however more promising since at low temperature (below 4K), jumps have been predicted [28,29]and may have already been observed in some experiment [8]. It would therefore be interesting in this case to perform controlled experiments for different relative sizes of the reservoir to see if the predicted reentrance of the isotherm as one moves away from the grand-canonical situation is indeed encountered, as observed in hysteretic martensitic transformations. [11]

2 3 4 0FIG. 2 :

 42 FIG.2: External potential Ψ * pore (r, z) for the heterogenous pore as a function of coordinates r and z.

FIG. 3 :FIG. 4 :

 34 FIG.3: Adsorption/desorption isotherms obtained by GCMC (α = 0) simulation for various temperatures (symbols). Lines are guides to the eye. Note the hysteresis present at low enough temperature.

6 FIG. 5 :

 65 FIG. 5: Enlargement of the low-temperature metastable states found by GCMC (α = 0). Solid lines indicate reversible branches. The lower panel shows at µ * = -9.52 (dotted line), for each of the five metastable states sketched in the upper panel, one typical molecular configuration.

3 FIG. 6 :

 36 FIG.6: Adsorption/desorption isotherms for various values of the relative size of the gas reservoir: α = 0 (GCMC), α ∼ 5.×10 -4 and 2. × 10 -3 , where α = VP /(VP + VR).

FIG. 7 :FIG. 8 :

 78 FIG. 7: Circles: GCMC results for T * = 0.60. Lines are guides to the eye. The arrows indicate the jump made by the representative point of the adsorbed fluid for various reservoir sizes, after destabilization of the initial point A by a small increase in the total number of particles in the system. The curved arrow indicates increasing values of the relative reservoir size: α = 0 (GCMC), 2.2 × 10 -4 , 2.5 × 10 -4 , 5.0 × 10 -4 , 1.0 × 10 -3 , 2.0 × 10 -3 , 2.4 × 10 -2 .

8 FIG. 9 :

 89 FIG. 9: Metastable states of the adsorbed fluid obtained for all the different sizes of the reservoir studied here (shown with the different colors used in Fig. 6; colors on line).

FIG. 11 :

 11 FIG. 11: Grand canonical (α = 0) descending scanning curves calculated by LMFT in a sample of linear size L = 50 at T * = 1.40. A close-up of the adsorption isotherm is shown in the inset.

FIG. 14 :

 14 FIG. 14: Canonical (symbols, α = 1) and grand-canonical (lines, α = 0) adsorption isotherms obtained by LMFT at T * = 1.10 for samples of different linear sizes: L = 25 (left), L = 50 (middle), L = 100 (right).

FIG. 17 :

 17 FIG. 17: Average canonical (α = 1) adsorption isotherms µ(ρP ) calculated by LMFT at T * = 0.80 for different system sizes; from left to right in the upper part of the plot (colors on line): L = 25 (black), L = 50 (red), L = 100 (blue). The corresponding numbers of realizations are 256, 64 and 8, respectively. Typical grand-canonical (α = 0) adsorption isotherms for L = 200 samples are also shown (dashed curves).

Acknowledgments

We thank M.-L. Rosinberg for many fruitful discussions.

Appendix: DFT Algorithms

As a starting point, it is convenient to rewrite the condition of minimization of the grand potential, ∂ΩT ∂ρi = 0, as exp(-βλ)

and sum over i so as to express λ as a function of the densities exp(-βλ

Our algorithm is then the following: changing the total density ρ by a small step ∆ρ, i.e. ρ new = ρ old + ∆ρ, one first supposes that the supplementary amount is adsorbed in the porous material: ρ new P = ρ old P + ∆ρ P . One then computes the new Lagrange parameter from exp(-βλ new ) = i {(η i -ρ old i ) exp[β(w f f j/i ρ old j + w sf j/i (1 -η j ))]} V P ρ new P (A.3) (with the old local densities) and the new local densities from

The new reservoir density ρ new R is obtained by inverting

within a given precision and, if not,

until convergence is reached. In practice, to improve the convergence of the iteration procedure, we use a mixing scheme, retaining a part of the previous iteration for the subsequent iteration. It appears that the configurations visited do not depend on the mixing parameter when the convergence criterions are strong enough.

Note that other equivalent algorithms can be devised as well and, interestingly, the output of the calculation appears to be quite robust to the choice of the algorithm. For instance, changing Eq. A.2 by the equivalent equation,

yields the same trajectory for the converged states even if the intermediate stages and the speed of convergence greatly differ. One can also add the increments ∆ρ in the reservoir as in the molecular simulation: ρ new R = ρ old R + ∆ρ; then, one computes the new Lagrange multiplier with Eq. A.5 and the new local densities with Eq. A.4, checks the constraint, and iterates using ∆ρ = 1 1-α [ρ -α

. This does not change the isotherm when the evolution is adiabatic. This algorithm seems simpler and possibly closer to the experimental protocol. However, it becomes problematic in the limit of the canonical ensemble (α=1). Therefore, we have preferred the first algorithm described above.

Our calculations were performed with samples of linear sizes varying from L = 25 to L = 100 and convergence was assumed when for the n th iteration, max

< ǫ and ρ (n-1) -

< ǫ with ǫ = 10 -6 for L = 25 and ǫ = 10 -8 for L = 50 and L = 100. In addition, the steps ∆ρ in ρ were taken as small as 10 -5 so that most of the avalanches could be resolved (see Ref. [1] for more details about identification of avalanches). Enforcing further the convergence criterions (e.g. with ǫ = 10 -9 for L = 50) does not change the path followed by the system: adiabatic regime has been reached. However, for the smallest system studied, with linear size L = 25, diminishing too much ∆ρ could prevent convergence as the system is not able to find a metastable state with the required density.


S Badiou

A Boudet

H Leray-Moragues

A Rodriguez

A S Bargnouxa

A M Dupuy

J P Cristol
email: jp-cristol@chu-montpellier.fr

Monthly reference change value of cardiac troponin in hemodialysis patients as a useful tool for long-term cardiovascular management☆

Keywords:

Purpose of research: Circulating cardiac troponin (cTn) has been identified as a risk factor for cardiovascular and overall mortality in patients undergoing hemodialysis. However, its interpretation remains difficult due to the high prevalence of patients with cTn level beyond the 99th percentile. Determining the cTn reference change value (RCV) may help in assessing a clinically significant change of cTn during regular follow-up of patients. We aimed to determine the long-term RCV of cTn in such patients and to calculate the perdialytic reduction rate of cTn. Design and methods: To calculate RCV, high-sensitivity (hs)-cTnT (Roche), hs-cTnI (Abbott), and cTnI-ultra (Siemens)were determined every month before the midweek dialysis session over a 3-month period in 36 stable hemodialysis patients. cTn was also measured after the midweek dialysis session to calculate the cTn removal rate. Results: The mean RCV (95% confidence interval) was 22% (18-26) for hs-cTnT versus 53% (34-73) for hs-cTnI versus 65% (45-84) for cTnI-ultra. Log-normal RCV (%) was -19/+25 for hs-cTnT, -33/+96 for hs-cTnI, and-39/+115 for cTnI-ultra. The index of individuality was b0.6 regardless of the cTn assay used. A significantly greater reduction rate was observed for hs-cTnT (48%) than for hs-cTnI (30%, p b 0.001) and cTnI-ultra (29%, p b 0.05). Conclusions: These results underline the need to use the RCV approach rather than cutoff points to identify the critical change in long-term serial cTn levels. In addition, RCV should be determined for each available assay due to significant differences between assays. Removal of cTn during hemodialysis sessions should also be considered if acute coronary syndrome is suspected during a session.

Introduction

Circulating cardiac troponin (cTn) level is commonly used as a biomarker for diagnosing myocardial infarction (MI) using the 99th percentile as a cutoff and considering kinetic measures [START_REF] Roffi | ESC guidelines for themanagement of acute coronary syndromes in patients presenting without persistent ST-segment elevation: task force for the Management of Acute Coronary Syndromes in patients presenting without persistent ST-segment elevation of the European Society of Cardiology (ESC)[END_REF]. With the emergence of high-sensitivity (hs)-cTn, circulating levels below the 99th percentile are detected in N50% of presumably healthy volunteers, thus extending its field of interpretation. Indeed, circulating cTn level has been suggested as a prognostic factor in different populations including the elderly [START_REF] Kuster | Estimation of age-and comorbidities-adjusted percentiles of high-sensitivity cardiac troponin T levels in the elderly[END_REF], those with chronic heart failure [START_REF] Jungbauer | Highsensitive troponin T in chronic heart failure correlates with severity of symptoms, left ventricular dysfunction and prognosis independently from N-terminal pro-btype natriuretic peptide[END_REF], and those undergoing hemodialysis (HD) [START_REF] Bargnoux | Cardiovascular risk stratification in hemodialysis patients in the era of highly sensitive troponins: should we choose between hs-troponin I and hs-troponin T?[END_REF]. It has been previously reported that HD patients evidenced chronic elevation of cTn levels [START_REF] Wayand | Cardiac troponin T and I in end-stage renal failure[END_REF], leading to inconclusive interpretation of a single cTn value above the 99th percentile.

Therefore, serial measures of circulating cTn are needed to distinguish between acute coronary syndrome (ACS) and chronic elevation of cTn in HD. However, the magnitude of cTn change in HD that is clinically relevant both for MI diagnosis and for the worsening of chronic cardiac disease remains to be defined. Indeed, parameters including renal removal impairment, adsorption onto dialyzer membranes [START_REF] Gaze | Cardiac troponin I but not cardiac troponin T adheres to polysulfone dialyser membranes in an in vitro haemodialysis model: explanation for lower serum cTnI concentrations following dialysis[END_REF], or change in troponin molecular structure [START_REF] Ziebig | Renal elimination of troponin T and troponin I[END_REF] could influence cTn levels, particularly cTnI. In the general population, algorithms based on the 99th percentile and/or relative and absolute changes of cTn during the kinetics have been proposed for MI diagnosis [START_REF] Roffi | ESC guidelines for themanagement of acute coronary syndromes in patients presenting without persistent ST-segment elevation: task force for the Management of Acute Coronary Syndromes in patients presenting without persistent ST-segment elevation of the European Society of Cardiology (ESC)[END_REF]. However, there are no recommendations to define the clinical change of long-term cTn in patients with high cardiovascular risk such as HD patients.

In this context, knowledge of cTn reference change value (RCV) is necessary to better interpret cTn levels in the regular follow-up of HD patients. The concept of RCV considers both the analytical and biological variations of a biomarker in order to determine which percentage change between two consecutive measurements may be considered as a critical difference [START_REF] Fraser | Generation and application of data on biological variation in clinical chemistry[END_REF]. Biological variations depend on intra-and inter individual variations assessed through coefficients of variation (CV) within (CVi) and between (CVg) subjects, respectively. The determination of RCV for a biomarker in the healthy population is not always transposable to specific clinical conditions, as recently demonstrated in HD patients [START_REF] Corte | Biological variation of cardiac troponin T in patientswith end-stage renal disease and in healthy individuals[END_REF]. To date, only limited data are available for longterm RCV of cTn in this population.

The aim of the present study was therefore to calculate long-term RCV of different cTn assays including cTnT and cTnI in patients with end-stage renal disease on stable HD maintenance, as well as to assess the reduction rate of cTn during an HD session.

Patients and methods

This longitudinal study included stable patients with end-stage chronic kidney disease undergoing hemodialysis thrice a week at the GCS HELP Dialysis Center at Montpellier (France). All patients enrolled were monitoring for their cTn every month for a 3-month period after receiving oral information and providing written consent. The study was approved by the local ethics committee and was declared under the number DC-2008-417.

Patients were included if they were clinically stable, without cardiovascular events or change in dialysis conditions during the study period. Peripheral blood samples were collected before dialysis during a midweek session into evacuated tubes containing heparin as an anticoagulant. In addition, samples were drawn before and after dialysis in the last month of the study, to assess the removal rate of cTn. The reduction rate was corrected for hemo concentration using Bergson's formula: [c]before -([c]after × (1 + (weight before -Weight after) / 0.2 × weight after)) / [c]before). Here, [c] corresponds to the concentration of the molecule [START_REF] Bergström | No change in corrected B2-microglobulin concentration after cuprophane haemodialysis[END_REF].

Hs-cTnT was measured by the chemiluminescent method using a Cobas8000/e602 analyzer and an hs-cTnT STAT reagent (Roche®) corresponding to our routine method. The samples were measured over the 3-month study period using multiple calibration curves. The manufacturer's claim was a CV10% at 13 ng/L, the 99th percentile at 14 ng/L, and a CVa at 2.7% at 17.7 ng/L. A previous study in our laboratory revealed a 2.4% CVa for hs-cTnT using a pool of plasma at 11 ng/L on a Cobas8000/e602 analyzer [START_REF] Dupuy | Biological variability of hs-cardiac troponin T on the Roche Cobas 8000/e602(R) immunoanalyzer[END_REF].

STAT hs-cTnI was measured on frozen aliquots with the chemiluminescent method on an Architect i1000 (Abbott®), which is used in our laboratory for other routine methods. All samples were measured using a single calibration curve. The manufacturer's claim was a CV 10% at 4.7 ng/L, the 99th percentile at 26.2 ng/L, and CVa of 3.4% for a level of 43 ng/L. cTnI-ultra was determined on frozen aliquots using an acridinium ester-based assay on an ADVIA Centaur (Siemens®) device used in our laboratory for other routine methods. All samples were measured using a single calibration curve. The manufacturer's claim was a CV10% at 30 ng/L, the 99th percentile at 40 ng/L, and CVa at 5.3% for a level of 80 ng/L.

RCV was calculated as RCV=21/2 * Z * CVt, with a Z-score of 1.96 for p b 0.05 and a CVt 2 = Cvi 2 + CVa 2. A log-normal RCV was calculated according to the formula proposed by Fokkema et al. [START_REF] Fokkema | Reference change values for brain natriuretic peptides revisited[END_REF], as log transformation is recommended for biomarkers with skewed distribution. The total CV (CVt) was individually calculated using the results of serial measures. CVi was obtained after removing the analytical variation to CVt according to the formula CVi=(CVt 2-CVa 2)1/2. The index of individuality (II) was assessed through the ratio CVi/CVg, with CVg calculated as the between-subject variation.

Results

Thirty-six HD patients fulfilled the inclusion criteria. The median age was 75.6 ± 17.2 years, and 27.8% (n = 10) were women. All patients were treated thrice a week by hemodiafiltration with a Fresenius 5008/5008Mix dialysis system (Fresenius Medical Care, Germany), using polysulfone (91.5%) (Fresenius FX 800 and 100 (86%), Asahi APS (5.5%)) and AN69 (Evodial, Hospal, 8.5%) membranes. The duration of the dialysis session was 220 ± 25 min with an intradialytic weight loss of 1.82 ± 1.02 kg and a Kt/Veq of 1.63 ± 0.25 according to the Daugiras equation.

At the first measurement, the median and interquartile (Q1;Q3) values for cTn (ng/L) were 76 (56;100), 23 (13;34), and 23 (15; 38) for the Roche, Abbott, and Siemens assays, respectively. The mean of cTn was higher than the median values for each cTn assay, revealing a rightskewed distribution. The cTn values were below the 99th percentile value provided by manufacturers (i.e., b14 ng/L for the Roche hs-cTnT, b26.2 ng/L for the Abbott hs-cTnI, and b40 ng/L for the Siemens cTnI-ultra) in, respectively, 0%, 62%, and 77% of patients. The level of cTn did not increase over the period study, regardless of the assay used. The individual kinetics of cTn for each assay are presented in Fig. 1. The mean reduction rate of cTn during a dialysis session calculated using a correction for hemo-concentration was significantly higher for cTnT (48% (43-51)) than for hs-cTnI (30% (22-37), p b 0.001) and for cTnI-ultra (29% (25-44), p b 0.05). We noted a return to value below the 99th percentile for four patients with hs-cTnT, five patients with hs-cTnI, and none for cTnI-ultra. The median cTn and interquartile of cTnI before and after the last dialysis session are presented in Fig. 1.

The results for the calculation of RCV, CVi, CVg, and II are presented in Table 1. The CVi values were significantly higher for both cTnI compared with cTnT, whereas II was b0.6 for all assays. The RCV for both cTnI (upper 50%) was significantly higher than the RCV for cTnT (22%). Asymmetric RCV using the log-normal approach confirmed the significantly higher RCV for cTnI compared to cTnT for both decrease and increase.

Discussion

This study highlights the difficulties of interpreting cTn levels in the HD population. Indeed, as previously reported, the cTn levels are elevated and frequently exceed the 99th percentile, which is generally recognized as an alert threshold value. In addition, RCV values and low II strongly suggest that serial measurements are needed to obtain clinically relevant information.

Recently, the European Society of Cardiology (ESC) proposed an alternative 0/1-h rapid ruleout and rule-in algorithm for non-ST elevation ACS based on basal levels or very short (1-h) kinetic studies. A rule-in cutoff value of 52 ng/L has been validated for the Roche hscTnT and been proposed for the Abbott hs-cTnI assays, based on a single hs-cTn determination [START_REF] Roffi | ESC guidelines for themanagement of acute coronary syndromes in patients presenting without persistent ST-segment elevation: task force for the Management of Acute Coronary Syndromes in patients presenting without persistent ST-segment elevation of the European Society of Cardiology (ESC)[END_REF]. For hs-cTnT, 81% of our HD cohort exceeded the threshold established in the general population, whereas 15% of our patients exceeded the threshold for hs-cTnI. These repartitions highlight the complexity of cTn interpretation in HD. Nevertheless, considering the complexity of their condition, HD patients should be considered for further workup. Furthermore, clinicians must particularly focus on the laboratory method used for cTn measurement. This is because its interpretation depends on the characteristics of the test used, including the cutoff point, RCV, and additional clearance rates for HD patients.

Beyond the diagnosis of MI, cTn was evidenced as a predictive factor of mortality in asymptomatic HD patients, even with classical cTnT [START_REF] Apple | Predictive value of cardiac troponin I and T for subsequent death in end-stage renal disease[END_REF] and hs-cTnT or cTnI [START_REF] Bargnoux | Cardiovascular risk stratification in hemodialysis patients in the era of highly sensitive troponins: should we choose between hs-troponin I and hs-troponin T?[END_REF][START_REF] Mcgill | Over time, highsensitivity TnT replacesNT-proBNP as themost powerful predictor of death in patients with dialysis-dependent chronic renal failure[END_REF]. Cardiovascular disease is known to be the first cause of mortality in HD, although often asymptomatic. Thus, regular follow-up of cTn may improve the clinical management of these high-risk patients. In this approach, determining long-term RCV of cTn is useful in assessing the significant change indicating worsening or occurrence of chronic cardiac disease. In the present study, RCV was significantly higher for both cTnI (Roche and Abott assays) compared with cTnT (Siemens assay), highlighting the need for determining RCV for each available cTn assay. Indeed, our results support that a change around 20% for hs-cTnT measured every month can be considered as a non-critical variation, whereas this change could be around 50% for hs-cTnI and 60% for cTnI-ultra. A more fine-tuned approach can be used with the log-transformed RCV, which allows distinguishing between variation of cTn in rise or fall.

The CVi and RCV observed are close to the weekly (hs-cTnT and hs-cTnI) and monthly (hs cTnT and sensitive cTnI) data previously reported in HD patients [START_REF] Aakre | Weekly and 90-minute biological variations in cardiac troponin T and cardiac troponin I in hemodialysis patients and healthy controls[END_REF][START_REF] Mbagaya | Biological variation of cardiac troponin in stable haemodialysis patients[END_REF]. They support the need to define different magnitude changes of cTn between HD patients and the general population, especially in absolute values. The use of cTn during HD sessions will be even more difficult if ACS is suspected. Indeed, the variation of cTn during the HD session may interfere in kinetic interpretation [START_REF] Wayand | Cardiac troponin T and I in end-stage renal failure[END_REF]. Although previous studies measured cTn before and after HD [START_REF] Wayand | Cardiac troponin T and I in end-stage renal failure[END_REF], only a few used hs-cTn [START_REF] Wolley | Variation in and prognostic importance of troponin T measured using a high-sensitivity assay in clinically stable haemodialysis patients[END_REF]. Contrary to previous data reporting no change in cTnT, the present study showed a mean removal of cTnT at 48% using hemo-diafiltration. Of note, this reduction rate, corrected for hemo-concentration, appeared higher than the first reported rate using the same hs-cTnT assay [START_REF] Wolley | Variation in and prognostic importance of troponin T measured using a high-sensitivity assay in clinically stable haemodialysis patients[END_REF]. Further, it is worth noting that lower removal of cTnI was observed compared to hs-cTnT. This finding was in discordance with previous data reporting higher adsorption onto polysulfone membrane dialysis for cTnI compared to cTnT [START_REF] Gaze | Cardiac troponin I but not cardiac troponin T adheres to polysulfone dialyser membranes in an in vitro haemodialysis model: explanation for lower serum cTnI concentrations following dialysis[END_REF]. ___________________________________________________________________________ Fig. 1. Individual evolution of cTn over the 3-month period (M1, M2, and M3) and effect of a dialysis session on cTn level for each assay. The box plot represents median, first and third quartile, as well as minimum and maximum values.

___________________________________________________________________________

The limitations of this study were the number of patients included and the relatively short follow-up period. However, this study was designed to provide information on three cTn assays in parallel. It was also difficult to ensure the cardiac stability of our population outside acute events. This was established through an overall clinical and biological stability. Table 1 Intra individual coefficient of variation (CVi, %), inter individual coefficient of variation (CVg, %)), index of individuality (II), and reference change value (RCV, %) for cardiac troponin (cTn) in hemodialysis patients. CVa was claimed by the manufacturer.

In conclusion, this study determined the long-term RCV of several cTn assays in patients on regular maintenance HD for improved use of cTn in a context of chronic elevation of cTn. Our data evidenced higher biological variations with the cTnI assay compared to cTnT, highlighting the need to determine RCV for all available assays for identifying the magnitude of a critical change. These results strongly suggest that a serial monitoring of cTn, for instance, on a monthly or quarterly basis, could be included in the follow-up of HD patients. The long-term evolution of cTn in this population could be interpreted according to the RCV determined in a stable HD population.