Dynamic mean field theory for lattice gas models of fluids confined in porous materials: higher order theory based on the Bethe-Peierls and path probability method approximations.

Dynamic mean field theory for lattice gas models of fluids confined in porous materials: Higher order theory based on the Bethe-Peierls and path probability method approximations John R. Edison and Peter A. Monson Citation: The Journal of Chemical Physics 141, 024706 (2014); doi: 10.1063/1.4884456 View online: http://dx.doi.org/10.1063/1.4884456 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/2?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Dynamics of capillary condensation in lattice gas models of confined fluids: A comparison of dynamic mean field theory with dynamic Monte Carlo simulations J. Chem. Phys. 138, 234709 (2013); 10.1063/1.4811111 Fluids in porous media. III. Scaled particle theory J. Chem. Phys. 134, 074503 (2011); 10.1063/1.3532546 Publisher’s Note: “Mean field kinetic theory for a lattice gas model of fluids confined in porous materials” [J. Chem. Phys.128, 084701 (2008)] J. Chem. Phys. 128, 149902 (2008); 10.1063/1.2904477 Mean field kinetic theory for a lattice gas model of fluids confined in porous materials J. Chem. Phys. 128, 084701 (2008); 10.1063/1.2837287 Role of elasticity forces in thermodynamics of intercalation compounds: Self-consistent mean-field theory and Monte Carlo simulations J. Chem. Phys. 116, 3083 (2002); 10.1063/1.1436472

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 129.22.67.107 On: Mon, 24 Nov 2014 05:44:41

THE JOURNAL OF CHEMICAL PHYSICS 141, 024706 (2014)

Dynamic mean field theory for lattice gas models of fluids confined in porous materials: Higher order theory based on the Bethe-Peierls and path probability method approximations John R. Edison and Peter A. Monson Department of Chemical Engineering, University of Massachusetts, Amherst, Massachusetts 01003-9303, USA

(Received 28 April 2014; accepted 9 June 2014; published online 10 July 2014) Recently we have developed a dynamic mean field theory (DMFT) for lattice gas models of fluids in porous materials [P. A. Monson, J. Chem. Phys. 128(8), 084701 (2008)]. The theory can be used to describe the relaxation processes in the approach to equilibrium or metastable states for fluids in pores and is especially useful for studying system exhibiting adsorption/desorption hysteresis. In this paper we discuss the extension of the theory to higher order by means of the path probability method (PPM) of Kikuchi and co-workers. We show that this leads to a treatment of the dynamics that is consistent with thermodynamics coming from the Bethe-Peierls or Quasi-Chemical approximation for the equilibrium or metastable equilibrium states of the lattice model. We compare the results from the PPM with those from DMFT and from dynamic Monte Carlo simulations. We find that the predictions from PPM are qualitatively similar to those from DMFT but give somewhat improved quantitative accuracy, in part due to the superior treatment of the underlying thermodynamics. This comes at the cost of greater computational expense associated with the larger number of equations that must be solved. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4884456] I. INTRODUCTION

The lattice gas model of fluids has been extensively used to investigate the behavior of fluids under confinement.1, 6–8, 19, 20, 35, 40, 42, 43 The chief advantage of this model lies in the ability to perform an extensive investigations with simulations and theoretical approaches in a computationally efficient manner. One theoretical approach in particular the lattice gas Mean Field Theory (MFT) has been used by us and others to study fluids confined in ordered and disordered porous materials. The advantage of the MFT over off lattice classical density functional theories is the ability to study three-dimensional systems of large sizes with less computational effort. A recent review discusses the use of lattice models and MFT to describe the equilibrium behavior of ordered and disordered porous materials.35 Also interesting and important is understanding the relaxation dynamics of a system approaching equilibrium. The motivation behind this is to develop a better understanding of the transport resistances to equilibration in adsorption isotherms and how these depend on the structure and surface properties of the porous material, as well as the dynamic behavior in the hysteresis region of adsorption/desorption isotherms. Although molecular dynamics is the obvious candidate to investigate this problem, the system sizes and the time scales required make it infeasible. In the past Kawasaki dynamics or vacancy hopping dynamics has been used to describe the relaxation dynamics of lattice gas fluids, including the nucleation and dynamics associated with capillary condensation, wetting and dying phenomena.23, 24, 38, 41, 42 These simulations can capture the essential physics of the system but they can also be computationally expensive.

0021-9606/2014/141(2)/024706/10/$30.00

With this in mind we recently developed a theory of relaxation dynamics for confined lattice fluids which we refer to as dynamic mean field theory (DMFT) or mean field kinetic theory.33 The evolution of the molecular density distribution of the system is obtained by applying a mean field approximation to the master equation that describes Kawasaki dynamics. Solution of this equation yields a time dependent density distribution approximating an average over an ensemble of Kawasaki dynamics Monte Carlo trajectories. We have also used the DMFT used to study the relaxation dynamics of fluids in pore networks, partially wetting systems, and partially/completely drying situations9–12, 33, 34 and applications of the approach have been explored by others as well.18, 22, 32 DMFT predicts the evolution of a system between two stable/metastable states described by equilibrium MFT after a change in the chemical potential. The ability to investigate static and dynamic behavior consistently in a single framework is a major main advantage of this approach. The DMFT is similar in spirit to the Cahn-Hilliard approach, where a square gradient free energy functional2–4 was incorporated into a diffusion equation and used the resulting equation was used describe the early stages of spinodal decomposition. The DMFT is based on the dynamic mean field theories of the Ising37 and binary alloy28 models first developed almost 20 years ago and reviewed recently by Gouyet et al.15 Prior to our work Aranovich et al.29–31 used it extensively to model diffusion in porous membranes and bulk fluids. In a recent paper we showed through comparisons with equilibrium and dynamic Monte Carlo simulations that the quantitative accuracy of DMFT as a theory of lattice model dynamics is to a significant extent limited by the underlying

141, 024706-1

© 2014 AIP Publishing LLC


024706-2

J. R. Edison and P. A. Monson

J. Chem. Phys. 141, 024706 (2014)

mean field approximation for the thermodynamics.13 This makes it worthwhile to investigate the accuracy of extensions of DMFT involving higher order approximations. In this paper we explore the use of the path probability method (PPM) that was introduced by Kikuchi and co-workers17, 21 to study relaxation dynamics in a binary alloy model, which is isomorphic with the lattice model considered here. The pair approximation in the path probability method provides a theory of the dynamics that in the equilibrium limit yields a description of the thermodynamics that is identical to the Bethe-Peierls approximation (BPA) or quasichemical approximation. The BPA was recently applied to lattice gas models of confined fluids by Salazar and Gelb.39 In this paper we apply the PPM to lattice fluids confined in a slit pore that is in contact with the bulk fluid and we study responses of the confined fluid to changes in the bulk fluid chemical potential. We find that the qualitative predictions from the PPM are similar to those from DMFT but the quantitative accuracy is improved. This improved accuracy comes at a cost due to the larger number of equations that must be solved in the PPM. The rest of the paper is outlined as follows. In Sec. II we describe the model, review the thermodynamic formulation in MFT and the BPA, and describe the theories of the dynamics represented by DMFT and the PPM, including a derivation of the steady state limit of the PPM. In Sec. III we describe some applications to slit pores, comparing DMFT and PPM with dynamics Monte Carlo simulations, and in Sec. IV we summarize our results and conclusions.

where ρ i is the average occupancy (density) of site i. The equilibrium density distribution of the system is obtained by minimizing the free energy with respect to ρ i , and this yields a set of nonlinear equations that can be solved for the ρ i ’s, ρi − kT ln ρi+a + φi − μ = 0 ∀i. (3) 1 − ρi a One can also study the system in the canonical N V T ensemble, by fixing the total density of the system. The implementation of the MFT in the canonical ensemble is discussed in our earlier papers.33 The BPA treats first nearest neighbor interactions via absolute occupancies, and the remainder via a mean field. A detailed derivation can be found in the work of Salazar and Gelb.39 Every site in the lattice is associated with both an occupancy probability or density, ρ i , and z pair occupancy probabilities pi, i + a , where z is the coordination number of the lattice. The pair occupancy probability pi, i + a , is the probability that a site i and its first neighbor i + a are both occupied. The grand free energy in the BPA is given by39 (1 − z)ρi ln ρi + (1 − ρi ) ln(1 − ρi ) = kT i

1 ln(1 − ρi − ρi+a + pi,i+a ) 2 a 1 − ρi − ρi+a + pi,i+a − ρi ln ρi − pi,i+a 1 − ρi − ρi+a + pi,i+a − ρi+a ln ρi+a − pi,i+a

pi,i+a (1 − ρi − ρi+a + pi,i+a ) + pi,i+a ln (ρi − pi,i+a )(ρi+a − pi,i+a ) 1 − pi,i+a + ρi φi − μ ρi . (4) 2 a i i +

II. MODEL AND METHODS A. Thermodynamics in MFT and the BPA

The Hamiltonian of a single occupancy lattice gas, with nearest neighbor interactions, in the presence of an external field is given as ni ni+a + ni φi , (1) H =− 2 i a i where is the nearest neighbor interaction strength, ni is the occupancy (0 or 1) of site with i denoting a set of lattice vectors, and a is the vector that denotes the set of nearest neighbors. The field imposed on site i by the confining solid is given by φ i . We use a simple cubic lattice in our calculations and the fluid molecules interact with the walls via a nearest neighbor interaction of strength −α. The simplest mean field approximation is the point approximation. In this approximation a particle present in any site i interacts with an effective field created by the neighboring sites rather than interacting with them directly.5 The ensemble most relevant to study the equilibrium behavior of confined fluids is the grand ensemble (μ,V ,T). The grand free energy of the system in MFT is given by ρi ln ρi + (1 − ρi ) ln(1 − ρi ) = kT i

−

ρi ρi+a + ρi φi − μ ρi , 2 a i

i

i

(2)

The density distribution at equilibrium is determined by minimizing the grand free energy ∂ =0 ∀ ∂ρi ∂ =0 ∂pij

i, (5)

∀ i, j.

This yields the following set of coupled nonlinear equations: (1 − ρi )(ρi − pi,i+a ) ρi kT ln + kT ln (1 − ρi ) ρi (1 − ρi − ρi+a + pi,i+a ) a +φi − μi = 0 ∀

i

(6)

and pi,i+a (1 − ρi − ρi+a + pi,i+a ) − e/kT = 0 ∀ i. (ρi − pi,i+a )(ρi+a − pi,i+a )

(7)

The latter equation can be solved explicitly for the pair occupancy probability in terms of the average site occupancies so


024706-3


J. Chem. Phys. 141, 024706 (2014)

that pi,i+a =

ρi + ρi+a + 2

1−

1 + 2(e/kT − 1)(ρi + ρi+a − 2ρi ρi+a ) + (e/kT − 1)2 (ρi − ρi+a )2

The equations from MFT and the BPA are solved numerically by iteration of Eqs. (3) and (6), respectively, starting from a suitable initial guess. Due to the analytical relationship between the single and pair occupancy probabilities in Eq. (8) the computation effort to solve the BPA and MFT is comparable,39 although this does not carry over into the dynamics. Typically we start at low chemical potential and assume a uniform low density in the system. We then obtain a series of solutions on an isotherm by starting each new state from the solution at the previous chemical potential in the sequence. We study sequences of states with both increasing chemical potential (adsorption) and decreasing chemical potential (desorption) to obtain the branches of the hysteresis loops in the isotherms. The system we consider is an ideal slit pore in contact with the bulk fluid and with α = 3.0. The geometry is shown in Figure 1. Since the density is uniform in the y-direction we can incorporate that into the theoretical calculations so that the density distribution is two-dimensional.

2(e/kT − 1)

where EijMF T =

⎧ ⎨0

EjMF T < EiMF T

⎩ E MF T − E MF T j i

EjMF T > EiMF T

and EiMF T = −

1. Dynamic mean field theory

∂ρi = − wi,i+a ({n})ni (1−ni+a )−wi+a,i ({n})ni+a (1−ni )t , ∂t a (9) where wi,i+a ({n}) is the transition probability for transitions from site i to site i + a for a configuration {n}. The occupancy

z

y

(8)

factors ni (1 − ni+a ) and ni+a (1 − ni ) impose the requirement that in a hopping move from site i to site j, site i must be occupied and site j unoccupied, and vice versa. In the mean field approximation Eq. (9) can be written as ∂ρi = − [wi,i+a ({ρ})ρi (1−ρi+a )−wi+a,i ({ρ})ρi+a (1−ρi )]. ∂t a (10) Given expressions for the transition probabilities, Eq. (10) can be solved to obtain {ρi } as a function of t, and is equivalent to that given by Matuszak et al.29 The mean field approximation for the Metropolis transition probabilities in Kawasaki dynamics yields29

wij ({ρ}) = w0 exp − EijMF T /kT , (11)

B. Dynamic behavior

The DMFT33 is based on earlier work on the Ising model by Penrose37 and the binary alloy model by Martin28 as well as later work by Aranovich et al.29–31 Gouyet et al.15 have presented a detailed review of using this approach to study relaxation dynamics in binary alloys. The DMFT gives an approximation to the time evolution of the density distribution averaged over an ensemble of Kawasaki dynamics Monte Carlo trajectories. So we begin with the master equation that describes Kawasaki dynamics and apply mean field approximation to it. The evolution of the ensemble average density at site i can then be expressed exactly in terms of the net fluxes from site i to its nearest neighbor sites i + a via

.

ρi+a + φi

(12)

(13)

a

is the energy associated with site i in mean field. w0 is the jump or transition rate in the absence of interactions and can be viewed as defining the timescale. There are two ways in which the mean field approximation appears in this dynamic theory. The first is the expression for the energies in the transition probabilities and is similar to the mean field approximation made in the free energy. The second approximation is in the occupancy factors ρi (1 − ρi+a ) and ρi+a (1 − ρi ), which imply that there are no occupancy correlations in the system dynamics. Using Eqs. (3) and (10) we obtain ∂ρi =− wi,i+a ({ρ})ρi (1−ρi+a )[1−exp{−(μi −μi+a )/kT }]. ∂t a (14) This expression shows that a limit of the DMFT equations where flux approaches zero is associated with uniform chemical potential throughout the system. For the applications studied here and in our other works DMFT is used to study the evolution in the density distribution.

H x

L FIG. 1. Schematic representation of a slit pore of length L and height H, in contact with the bulk phase. Periodic boundaries are applied in the y-direction.


024706-4


J. Chem. Phys. 141, 024706 (2014)

equation as

2. Dynamics via the path probability method

Again our approach closely follows the work of Gouyet et al.15 who presented the PPM to describe the relaxation dynamics for a model of binary alloys with transition probabilities given by an egg box potential. We extend their work to a lattice gas model with transition probabilities given by the Metropolis criterion. In this case we have to consider the evolution of both the local density, ρ i , and the pair occupancy probability, pi, j . We begin with the evolution of the local density from Eq. (9). The flux directed from site i to l is given by Jildir = wil ({n})ni (1 − nl )t ,

(15)

where wij ({n}) is the transition probability for transitions from site i to site l for a configuration {n}. The angled brackets in the above equation imply an ensemble average over independent trajectories at an instant of time t. In DMFT we replace the average in the right-hand side, by the product of three averages, i.e., Jildir = wil ({ρ})ρi (1 − ρl ). Moreover, the ansatz employed to compute the average transition probability in DMFT is to simply replace the actual occupancy ni , with the average occupancy ρ i in the expression for the Metropolis criterion. In the pair approximation the above factorization cannot be done, since it implies that sites which are nearest neighbors are uncorrelated. Following Gouyet we employ the following ansatz to factorize the product in Eq. (15): Jildir = wil ({n})ni (1 − nl ).

(16)

In the above equation the term within the second set of angled brackets on the right-hand side is the probability of observing a pair of sites il, with site i being occupied and site l being vacant. The first term within the angled brackets on the right-hand side is the Metropolis transition probability given by wil ({n}) = wo exp(−Eil /kT ).

(17)

Here wo is the jump rate in the absence of any interactions (infinite temperature) and sets the time scale. The energy term Eil which is the change in energy associated with a translation of a particle to one of its neighboring sites is given by 0 El < Ei Eil = . (18) El − Ei El > Ei Here Ei is the interaction energy experienced by a particle at site i, with site l being empty, and El is the interaction energy experienced by the particle at site l with site i being empty. We only need to consider the case El > Ei since if El < Ei then wil ({n}) = wo . In the PPM we compute the average of the transition probability as follows: Ei − El exp(Ei /kT )P P M exp = . (19) kT exp(El /kT )P P M PPM Next we express the above equation in terms of point and pair occupancy probabilities. We can, for example, rewrite the term in the numerator of the right-hand side of the above

exp(Ei /kT )P P M =

ni exp − j=l nj + φi /kT ni

. (20)

The form of this equation guarantees that only contributions where ni = 1 are considered. The index j represents the nearest neighbors of site i excluding site l, which is empty. The exponential then is expanded as a weighted average to each of the two occupancy possibilities of site j, 0 (vacant) and 1 (occupied) and this gives exp(Ei /kT )P P M exp(−/kT )ni nj + ni (1 − nj ) φi = exp , (21) kT ρi j=l

which simplifies to

φi exp(−/kT )pij +ρi −pij exp(Ei /kT )P P M = exp . kT ρi

j=l

(22) In summary the Metropolis transition probability within the PPM can be written as ⎧ e−/kT p +ρ −p φ /kT ⎪ ij i ij i ⎪ e ρ ⎪ j=l i ⎨ w0 φ /kT e−/kT p +ρ −p wil < w0 . (23) wil = lk l lk l e ρ k=i ⎪ l ⎪ ⎪ ⎩ w0 wil ≥ w0 Note that the way we factorize Eq. (19) is not unique. One other possibility would be to write Ei − El exp(−El /kT )P P M exp = . (24) kT exp(−Ei /kT )P P M PPM Though Eqs. (19) and (24) appear equivalent, only Eq. (19) guarantees that the long time limit of the evolution equations to approach the equilibrium solutions obtained from the BPA. Now we turn to the pair occupancy probability. The evolution of the pair occupancy probabilities is described by the following exact equation (detailed derivation is given in Appendix A.2.2 of Ref. 15): ∂p ∂nl nk = lk ∂t ∂t = wil nk ni (1 − nl ) − wli nk nl (1 − ni ) i=k

+

i=k

wjk nl nj (1 − nk ) −

j=l

wkj nk nl (1 − nj ). j=l

(25) Following Gouyet the above equation which describes the evolution of pair occupancy probability of a pair of sites lk can be represented graphically as ∂ •l •i → ◦l ◦i ← •l − = •k •k ∂t •k i=k •l •l − . (26) + •j → ◦k ◦j ← •k j=l



J. Chem. Phys. 141, 024706 (2014)

wil nk ni (1 − nl ) =

nk (1 − nl ) wil ni (1 − nl ). (1 − nl )

(27)

The first term on the right-hand side is the conditional probability of observing site l in the pair lk vacant. The average wil ni (1 − nl ) is simply the flux from site i to site l, with site l described by (16), wil nk ni (1 − nl ) =

ρk − plk dir Jil ({n}) t . (1 − ρl )

(28)

The presence of a particle in site k is accounted for by carefully evaluating the Metropolis transition probability. The term e−/kT plk + ρ l − plk which is the contribution to the average Metropolis transition probability (see Eq. (23)) from the pair lk is replaced with e−/kT plk since site k is occupied. The third term in Eq. (26) can be computed in the same spirit and the summation runs over sites j which represents the nearest neighbors of site k excluding site l. The second (fourth) term in Eq. (26) is merely the flux from site l to site i (site k to site j). The presence of a particle in site k (site l) is accounted by evaluating the contribution to the average Metropolis criterion in the manner described above. In the long time limit the net flux between any two pair of sites vanishes, i.e., Jildir = Jlidir . Using Eqs. (16) and (19) and upon further simplification it can be shown that −/kT e pij +ρi −pij eφi /kT j=l ρi (ρ − pli ) wil . (29) = l = e−/kT plk +ρl −plk wli (ρi − pil ) eφl /kT k=i ρ l

Using Eqs. (6) and (7) it can be shown that the above equation reduces to exp[(μi − μl )/kT ] = 1,

1.5

1

1.25

0.75

1

MFT BPA Series approx.

0.75 0.5

Tr

A filled dot represents an occupied site while an open dot represents a vacant site. Terms 1 and 3 in the above Eq. (26) are pair creation terms, and terms 2 and 4 are pair annihilation terms. Consider term 1 on the right-hand side of Eq. (26). This term describes the creation of a pair lk, with the empty site l filled via a hop from one of its neighboring sites i. One way to factorize the above equation, in terms of pair and point occupation probabilities, is as follows:

T

024706-5

MFT BPA Series approx.

0.5 0.25

0.25 0

0 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

ρ

ρ

(a)

(b)

0.8

1

FIG. 2. Phase coexistence diagram of the bulk lattice gas (a) T vs ρ and (b) Tr = T/Tc vs ρ, computed the MFT, BPA methods and compared with the series approximation results of Essam and Fisher.14

probabilities at the interface between the system and the control volume evolve with time. For the DMFT we have found Euler’s method to be of acceptable accuracy for time steps less than about w0 t = 0.1. The PPM equations are stiffer and require smaller time steps, so we have used w0 t = 0.02 for the cases presented here.

C. Monte Carlo simulations

We have made Monte Carlo simulations of the systems using the grand ensemble (GCMC) method for static properties and dynamic Monte Carlo (DMC) simulations for the dynamics. These calculations follow the methods presented in our recent paper13 and we just highlight some of the essential features. The dimensions of the system are as follows. The wall spacing in the slit is H = 6 (in units of the lattice constant) and the length of the slit was either L = 20 or L = 40. The slit pore has periodic boundary conditions in the y-direction and the pore width in the y-direction is 12 lattice constants. The dimensions of the bulk region on each side of the slit are 20 (length) × 8 × 12 lattice constants. To study the dynamics of the model we use a dual control volume simulation technique16, 25, 36 with the bulk regions at each end of the pore divided into two regions of length 10 lattice units. The region furthest from the pore is the control volume. The simulation protocols for GCMC and DMC are identical to those we used earlier.13 The DMC simulation results presented in

(30)

which states that the long time limit of these equations corresponds to uniform chemical potential throughout the system. Also at the steady state the first two terms in Eq. (26) can be equated to zero. Using Eq. (29) we can show that the first two terms at the steady state reduce to Eq. (7). Thus the states predicted by the PPM in the long time limit are consistent with the states predicted by the BPA. We implement the DMFT and the PPM as follows. For a given pore geometry we add a layer of sites at the periphery of the system (control volume), where the density is fixed at the value associated with the activity of the state to which we want the system to evolve. The fixed density layer or the control volume acts as a source/sink of fluid during the dynamics. The system is then evolved by numerically integrating the evolution equations using an explicit integrator. In case of the PPM after every single integration step the BPA equations are resolved in the control volume, since the pair occupancy

2 MFT BPA GCMC 1.5 ρ

1

0.5

0 0

0.25

0.5 λ/λ0

0.75

1

FIG. 3. Isotherms of density vs. relative activity for a slit pore of width H = 6 and length L = 40 computed with MFT (black), BPA (blue), and GCMC (red). For clarity the curves from the BPA and GCMC have been shifted in density by 0.5 and 1.0, respectively.



J. Chem. Phys. 141, 024706 (2014)

1

1

0.75

0.75 ρ

ρ

024706-6

0.5

0.25

DMFT PPM DMC

0.25

10000

20000

PPM DMC

0

0 0

0.5

30000

0

20000

30000

-t

-t FIG. 4. Density evolution curves for a step change in relative activity from λ/λ0 = 0.001 to λ/λ0 = 0.95 in a L = 20, H = 6 slit pore at the same reduced temperature T/Tc = 2/3 predicted by the DMFT (black lines), PPM (blue lines), and DMC (red line). Density of the entire pore: solid line. Density of a layer of sites in the middle of the pore: dashed line.

10000

FIG. 6. Density evolution curves for a step change in relative activity from λ/λ0 = 0.001 to λ/λ0 = 0.95 at the same temperature T∗ = 0.822 for a L = 20, H = 6 slit pore predicted by the PPM (blue lines) and DMC (black line). We do not show the DMFT results as MFT does not predict a capillary condensation transition below saturation at this temperature.

B. Static behavior

this work were obtained by averaging data from 1536 independent DMC trajectories.

III. RESULTS AND DISCUSSION A. Bulk phase behavior

A study of relaxation dynamics of confined fluids requires knowledge of the static behavior of the bulk and the confined fluid and so first we show the bulk phase diagram of the lattice gas model computed with the MFT and BPA, compared with the highly accurate results from series expansions by Essam and Fisher.14 As shown in Figure 2 inclusion of correlations at just the pair level, decreases the critical temperature of the system from Tc∗ = 1.5 (MFT) to Tc∗ = 1.233 (BPA), much closer to the critical temperature of the lattice gas which is Tc∗ = 1.128. Figure 2(b) shows the three binodals scaled with their respective critical temperatures. With these differences in mind, in most cases we will compare the results of the DMFT and PPM methods with DMC at the same reduced temperature Tr = T/Tc , in order to reduce the impact of differences in the thermodynamic properties upon the results.

The comparison is done in the first instance at a fixed reduced temperature T/Tc = 0.66. In Figure 3 we plot isotherms of a slit pore of width H = 6 computed with the MFT and the BPA methods together with results from GCMC simulations. The fluid in the pore evaporates at a higher activity than would happen for an infinite pore without contact with the bulk.33 Contact with the bulk reservoir means there is no nucleation barrier to evaporation at the vapor-liquid interface and the desorption occurs at the equilibrium transition between the vapor and liquid phases in the pore, an observation first made by Marconi and van Swol.26, 27 The isotherms predicted by all three methods agree qualitatively, with the BPA isotherm showing the monolayer and the capillary condensation transition at a lower relative activity than the MFT isotherm. The hysteresis loops are wider for the two theories than for the GCMC results due to the neglect of fluctuations in the theories. C. Dynamic behavior

1. Pore filling in an L = 20 pore

We first show the mechanism of pore filling in an L = 20 slit pore. We initialize the system at a very low relative

FIG. 5. Visualizations of states emergent in the dynamics during a step in change in relative activity from λ/λ0 = 0.001 to λ/λ0 = 0.95 at the same reduced temperature T/Tc = 2/3 for a L = 20 slit pore predicted by the DMFT (left panel), the PPM (center panel), and DMC (right panel).


024706-7


J. Chem. Phys. 141, 024706 (2014) 1

ρ

0.75 0.5 DMFT PPM DMC

0.25 0 0

20000

40000

60000

80000

-t FIG. 8. Density evolution curves for a step change in relative activity from λ/λ0 = 0.001 to λ/λ0 = 0.95 at the same reduced temperature T/Tc = 2/3 for a L = 40, H = 6 slit pore predicted by the DMFT (black lines), PPM (blue lines), and DMC (red line). Density of the entire pore: solid line. Density of a layer of sites in the middle of the pore: dashed line. FIG. 7. Visualizations of states emergent in the dynamics during a step in change in relative activity from λ/λ0 = 0.001 to λ/λ0 = 0.95 at the same temperature T∗ = 0.822 for a L = 20, H = 6 slit pore predicted by the PPM (left panel) and DMC (right panel).

activity of λ/λ0 = 0.001 and make a step change in the relative activity of the control volume to λ/λ0 = 0.95. The relative activity of λ/λ0 = 0.95 corresponds to a condensed state in the pore, and hence the dynamics involves a condensation process. The density evolution curves from the DMFT, the PPM, and DMC are shown in Figure 4 and show the nucleation of the liquid phase via the formation of a liquid bridge as was discussed in detail in our first DMFT paper.33 The full line is the density averaged over the entire pore and the dashed line is the density averaged over two layers of sites in the x-z plane equidistant from the ends of the pore. Both theories show a cusp in the density evolution of the pore which coincides with a sharp rise in density in the middle of pore. Note that the DMFT predicts an enhancement of density in the middle of the pore prior to the formation of a liquid bridge, corresponding to formation of an undulate and this is less clear in the results of the PPM. The density uptake curve computed from DMC simulations is qualitatively similar to those from DMFT and the PPM but the timescale is longer than both. This is in part due to the fact that while the values of T/Tc are the same in the three cases the values of T∗ = kT/ are in the sequence ∗ ∗ ∗ < TP∗P M < TDMF TDMC T . For lower values of T the dynam13 ics will be slower. Note that pore averaged density curve in

Figure 4 does not show a cusp like that seen for DMFT and the PPM. This is due to the distribution of nucleation times in the DMC rather than the single nucleation time seen in the theory. This issue was discussed in detail in our earlier paper comparing DMFT with DMC simulations.13 The visualizations of the states encountered in the dynamics are shown in Figure 5. Overall it is evident that there is good qualitative agreement between the three sets of results. The DMC results show a more diffuse distribution reflecting the broader temporal and spatial distributions of the bridge nucleation events.13 We have been able to compare results from the PPM and DMC at the same value of T∗ = 0.822. (At this temperature the hysteresis is much wider in DMFT and we have to take the bulk fluid to supersaturated states in order to see condensation.) The results for the density vs. time are shown in Figure 6 and visualizations shown in Figure 7. We see that the agreement is now somewhat closer although the DMC results are again on a somewhat longer timescale. 2. Pore filling in an L = 40 pore

Next we compare the filling of a slit pore of length L = 40 from a dilute state to a completely filled state. The density evolution curves predicted by the DMFT and PPM are shown in Figure 8 together with that from DMC. The visualizations of the states encountered in the dynamics are shown in Figure 9. Dynamics at short times is associated with the

FIG. 9. Visualizations of states emergent in the dynamics during a step in change in relative activity from λ/λ0 = 0.001 to λ/λ0 = 0.95 at the same reduced temperature T/Tc = 2/3 for a L = 40, H = 6 slit pore predicted by the DMFT (left panel), the PPM (center panel), and DMC (right panel).


024706-8


J. Chem. Phys. 141, 024706 (2014)

1

1 DMFT PPM DMC

0.5 0.25

DMFT PPM DMC

0.75 ρ

ρ

0.75

0.5 0.25

0

0 0

15000

30000

-t

0

4000

8000

-t

FIG. 10. Density evolution curves for a step change in relative activity from λ/λ0 = 0.95 to λ/λ0 = 0.001 at the same reduced temperature T/Tc = 2/3 for a L = 20, H = 6 slit pore predicted by the DMFT (black lines), PPM (blue lines), and DMC (red line).

FIG. 12. Density evolution curves for a step change in relative activity from λ/λ0 = 0.95 to λ/λ0 = 0.001 at the same temperature T∗ = 1.0 for a L = 20, H = 6 slit pore predicted by the DMFT (black lines), PPM (blue lines), and DMC (red line).

formation of a liquid-like mono-layer at the pore walls. It is then followed by the appearance of liquid bridges at each end of the pore which then proceed to fill the pore. The appearance of two bridges rather than one is related to the choice of the final state chosen on the liquid branch of the isotherm and also the difference in chemical potential between the initial and final states.13, 33 The visualizations clearly indicate that the predictions of the DMFT qualitatively agree with the PPM technique. The density evolution curve also agrees qualitatively with each other. Again there are subtle differences when we compare the density evolution of the middle of the pore (dashed curves) shown in Figure 8. In case of the DMFT this curve separates from the overall density curve before the cusp corresponding to the formation of the liquid bridge. This indicates that the undulates that precede the liquid bridges are formed by temporary redistribution of density from the middle of the pore to the ends of the pore. The results of the PPM, however, show that this separation happens at the cusp or when the liquid bridges are formed in the pore. The loss of density in the middle of the pore during the instant the bridge is formed (cusp of the overall pore density curve) is predicted by both the methods. This is due to the bridge attracting fluid from elsewhere in the pore as well as from the bulk. In the middle of the pore there is a period of time when the flux is reversed relative to the overall flux of fluid into the

pore. The failure to observe pronounced undulate states in the pore in the PPM might be a consequence of the difference in temperatures used in the two methods. The average density throughout the pore from DMC is also shown in Figure 8 and again we see qualitatively similar behavior to the results from theory but the timescale is again longer. The visualizations from DMC are also shown in Figure 9 and we see that multiple bridging is also seen in this case. Due to the presence of fluctuations in the DMC we also have the possibility that additional bridges could appear in the system and we have explored that possibility in other recent work.13 What we found in that work is that for even very long pores such additional bridges are quickly absorbed into the two bridges near the pore ends. 3. Pore emptying in an L = 20 pore

We have also studied the dynamics of pore emptying in an L = 20 slit pore. We initialize the system at λ/λ0 = 0.95, and make a step change in relative activity of the bulk to λ/λ0 = 0.001. The pore emptying proceeds via the meniscus receding to the middle of the pore and then the snapping of the liquid bridge followed by desorption of the fluid layers near the pore walls. The density evolution curves from the DMFT, PPM, and DMC are shown in Figure 10. Interestingly the

FIG. 11. Visualizations of states emergent in the dynamics during a step in change in relative activity from λ/λ0 = 0.95 to λ/λ0 = 0.001 at the same reduced temperature T/Tc = 2/3 for a L = 20 slit pore predicted by the DMFT (left panel), the PPM (center panel), and DMC (right panel).


024706-9


J. Chem. Phys. 141, 024706 (2014)

FIG. 13. Visualizations of states emergent in the dynamics during a step in change in relative activity from λ/λ0 = 0.95 to λ/λ0 = 0.001 at the same temperature T∗ = 1.0 for a L = 20, H = 6 slit pore predicted by the DMFT (left panel), the PPM (center panel), and DMC (right panel).

differences in timescales are apparently much greater in this case even though the mechanisms are similar as shown in the visualizations in Figure 11. For desorption we have been able to compare all three methods at the same value of T∗ and this is shown in Figures 12 and 13. The agreement is now very close, suggesting that the earlier shown differences are simply due to the effect of using a different T∗ but the same T/Tc value upon the dynamics.

IV. SUMMARY AND CONCLUSIONS

We have presented an extension of the DMFT for lattice gas models of confined fluids to higher order based on the PPM approach. The PPM theory generates a theory of the dynamics that is consistent with the Bethe-Peierls approximation for the thermodynamics, a more accurate thermodynamics theory than the mean field approximation underlying DMFT. We have applied the PPM to a slit pore in contact with a bulk gas as studied in earlier work.13, 33 We have studied the dynamics of uptake from an empty to a filled pore for two pore lengths as well as the dynamics of desorption. Overall the results from the PPM theory are qualitatively similar to those from DMFT and the description of the condensation and evaporation phenomena in the pore is also consistent with DMFT. The agreement with DMC is somewhat better for the PPM theory in part due to there more accurate treatment of the thermodynamics (BPA vs. MFT). In assessing respective utilities of the PPM and DMFT theories of the lattice gas model dynamics we can say that they both provide a good qualitative description of the dynamics. The PPM provides somewhat improved quantitative predictions. However, unlike the case of the MFT and BPA where the computational costs are similar, for the dynamics the PPM is a much more expensive calculations. This is due to two features. The first is that we must explicitly track the dynamics of the pair occupancy probabilities and there are z such functions for each location in the system, where z is the coordination number of the lattice. This leads to (z + 1) times as many equations to solve. In addition we find that the evolution equations for the PPM are somewhat stiffer than

those for the DMFT and stability dictates the use of a smaller time step. As a consequence we find that PPM calculations typically require in the range of 10–50 times the computation time required for the DMFT. It is gratifying that the two theories give similar qualitative predictions and the much lower computational costs of the DMFT suggest it as the method of choice in these studies. On the other hand there may well be situations where observed phenomena depend on a close balance of different interaction effects and the PPM may prove superior in such cases. Tables of data from Figures 2–4, 6, 8, 10, and 12 of this paper and gnuplot scripts for plotting them are available as the supplementary material.44 ACKNOWLEDGMENTS

This work was supported by the National Science Foundation (Grant Nos. CBET-0853068 and CBET-1158790). 1 G.

Aranovich and M. Donohue, J. Colloid Interface Sci. 200(2), 273 (1998). 2 J. W. Cahn, J. Chem. Phys. 30(5), 1121 (1959). 3 J. W. Cahn and J. E. Hilliard, J. Chem. Phys. 28(2), 258 (1958). 4 J. W. Cahn and J. E. Hilliard, J. Chem. Phys. 31(3), 688 (1959). 5 D. Chandler, Introduction to Modern Statistical Mechanics (Oxford University Press, 1987). 6 M. J. De Oliveira and R. B. Griffiths, Surf. Sci. 71(3), 687 (1978). 7 M. D. Donohue and G. L. Aranovich, J. Colloid Interface Sci. 205(1), 121 (1998). 8 C. Ebner, Phys. Rev. A 23(4), 1925 (1981). 9 J. Edison, M. Ganz, B. Novello, and P. Monson, Adsorption 17, 769 (2011). 10 J. Edison and P. Monson, Microporous Mesoporous Mater. 154(0), 7 (2012). 11 J. R. Edison and P. A. Monson, J. Low Temp. Phys. 157(3-4), 395 (2009). 12 J. R. Edison and P. A. Monson, Faraday Discuss. 146, 167 (2010). 13 J. R. Edison and P. A. Monson, J. Chem. Phys. 138(23), 234709 (2013). 14 J. W. Essam and M. E. Fisher, J. Chem. Phys. 38(4), 802 (1963). 15 J. F. Gouyet, M. Plapp, W. Dieterich, and P. Maass, Adv. Phys. 52(6), 523 (2003). 16 G. Heffelfinger and F. van Swol, J. Chem. Phys. 100(10), 7548 (1994). 17 T. Ishii, H. Sato, and R. Kikuchi, Phys. Rev. B 34, 8335 (1986). 18 E. Kierlik, F. Leoni, M. L. Rosinberg, and G. Tarjus, Mol. Phys. 109, 1143 (2011). 19 E. Kierlik, P. A. Monson, M. L. Rosinberg, L. Sarkisov, and G. Tarjus, Phys. Rev. Lett. 87(5), 055701 (2001).


024706-10 20 E.


Kierlik, P. A. Monson, M. L. Rosinberg, and G. Tarjus, J. Phys.: Condens. Matter 14(40), 9295 (2002). 21 R. Kikuchi, Prog. Theor. Phys. Suppl. 35, 1 (1966). 22 F. Leoni, E. Kierlik, M. L. Rosinberg, and G. Tarjus, Langmuir 27(13), 8160 (2011). 23 K. Leung and A. Luzar, J. Chem. Phys. 113(14), 5845 (2000). 24 A. Luzar and K. Leung, J. Chem. Phys. 113(14), 5836 (2000). 25 J. MacElroy, J. Chem. Phys. 101(6), 5274 (1994). 26 U. M. B. Marconi and F. van Swol, Phys. Rev. A 39(8), 4109 (1989). 27 U. M. B. Marconi and F. van Swol, Europhys. Lett. 8(6), 531 (1989). 28 G. Martin, Phys. Rev. B 41(4), 2279 (1990). 29 D. Matuszak, G. L. Aranovich, and M. D. Donohue, J. Chem. Phys. 121(1), 426 (2004). 30 D. Matuszak, G. L. Aranovich, and M. D. Donohue, Phys. Chem. Chem. Phys. 8(14), 1663 (2006). 31 D. Matuszak, G. L. Aranovich, and M. D. Donohue, Ind. Eng. Chem. Res. 45(16), 5501 (2006). 32 Y. Men, X. Zhang, and W. Wang, J. Chem. Phys. 134(12), 124704 (2011). 33 P. A. Monson, J. Chem. Phys. 128(8), 084701 (2008).

J. Chem. Phys. 141, 024706 (2014) 34 P.

A. Monson, Proceedings of the 8th International Symposium on the Characterization of Porous Solids, Special Publications No. 318, edited by S. Kaskel, P. Llewellyn, F. Rodriguez-Reinoso, and N. A. Seaton (Royal Society of Chemistry, 2009), pp. 103–110. 35 P. A. Monson, Microporous Mesoporous Mater. 160, 47 (2012). 36 A. Papadopoulou, E. D. Becker, M. Lupkowski, and F. van Swol, J. Chem. Phys. 98(6), 4897 (1993). 37 O. Penrose, J. Stat. Phys. 63(5–6), 975 (1991). 38 F. Porcheron and P. A. Monson, Langmuir 21(7), 3179 (2005). 39 R. Salazar and L. Gelb, Phys. Rev. E 71(4), 041502 (2005). 40 L. Sarkisov and P. A. Monson, Langmuir 17(24), 7600 (2001). 41 R. Valiullin, S. Naumov, P. Galvosas, J. Karger, H. J. Woo, F. Porcheron, and P. A. Monson, Nature (London) 443, 965 (2006). 42 H. J. Woo and P. A. Monson, Phys. Rev. E 67(4), 041207 (2003). 43 H. J. Woo, L. Sarkisov, and P. A. Monson, Langmuir 17(24), 7472 (2001). 44 See supplementary material at http://dx.doi.org/10.1063/1.4884456 for tables of data from Figures 2, 3, 4, 6, 8, 10, and 12 and Gnuplot scripts for plotting them. For information on the supplementary material, see http://www.aip.org/pubservs/epaps.html.


Beyond mean field theory: statistical field theory for neural networks.

Embedded mean-field theory.

Effective field theory for lattice nuclei.

Cluster mean-field theory study of J1-J2 Heisenberg model on a square lattice.

Field-based dynamic light scattering microscopy: theory and numerical analysis.

Real space mean-field theory of a spin-1 Bose gas in synthetic dimensions.

Empowerment theory: clarifying the nature of higher-order multidimensional constructs.

Symmetry-preserving mean field theory for electrostatics at interfaces.

Comments on quantum probability theory.

Electron inelastic mean free path theory and density functional theory resolving discrepancies for low-energy electrons in copper.

A transition state theory for calculating hopping times and diffusion in highly confined fluids.

Constructor theory of probability.

Calorimetric glass transition in a mean-field theory approach.

The mobile order theory versus UNIFAC and regular solution theory-derived models for predicting the solubility of solid substances.

An analysis of item response theory and Rasch models based on the most probable distribution method.

Zero-Temperature, Mean-Field Theory of Atomic Bose-Einstein Condensates.

Probability state modeling theory.

Transition path theory.

A bond-order theory on the phonon scattering by vacancies in two-dimensional materials.

A mean-field theory on the differential capacitance of asymmetric ionic liquid electrolytes.

Adapting SAFT-γ perturbation theory to site-based molecular dynamics simulation. II. Confined fluids and vapor-liquid interfaces.

Intrinsic flexibility of porous materials; theory, modelling and the flexibility window of the EMT zeolite framework.

Classifying linearly shielded modified gravity models in effective field theory.

Landau-like theory for universality of critical exponents in quasistationary states of isolated mean-field systems.