Table Of Content

Capillary condensation in a square geometry with surface fields M. Zubaszewska,1 A. Gendiar,2 and A. Drzewin´ski1 1Institute of Physics, University of Zielona Góra, ul. Prof. Z. Szafrana 4a, 65-516 Zielona Góra, Poland and 2Institute of Physics, Slovak Academy of Sciences, SK-845 11 Bratislava, Slovakia (Dated: January 3, 2013) We study the influence of wetting on the capillary condensation for a simple fluid in a square geometry with surface fields, where the reference system is an infinitely long slit. The Corner 3 Transfer Matrix Renormalization Group method has been extended to study a two-dimensional 1 Ising model confined in an L×L geometry with equal surface fields. Our results have confirmed 0 that in both geometries thecoexistence lineshift is governed bythesame scaling powers, buttheir 2 prefactors are different. n a PACSnumbers: 64.60.an,64.60.De,68.35.Rh,05.10.-a J 1 Porousmaterialsaresolidsconsistingofinterconnected out spins is identified with the wall. The above equation ] network of pores. In recent years, microporous and is known in the literature as the Kelvin equation7. h mesoporous materials have been a focus of nanoscience c andnanotechnologyastheirpropertiesdiffersignificantly e m from the same bulk materials1,2. Both the pore size and 2.5 its shape, as well as the chemical nature of its surface, - t i.e., whether hydrophilic or hydrophobic, determine the a 2 t properties of porous materials3. In a mesopore, the cu- s mulativeeffectofthewallsbecomesimportant. Afterthe . e t formation of adsorbate layers with 2-3 molecular thick- ur1.5 a at h (T) m ness on walls, further adsorption induces attractions be- per hcw (T) d- tswateieonnaodfslioqrubiadt-elikmeoaledcsuolrebsa,tleeamdoinlegcutolesaisnusdiddeenthceopnodreens-. Tem 1 hhw1 == 00..120705 1 n This effect is analogous to the capillary condensation h1 = 0.371 o phenomena4. 0.5 h1 = 0.500 c The Ising square of a finite size L with the field h [ 1 acting on all four surfaces (boundaries) can be used as 0 0 0.2 0.4 0.6 0.8 1 3 an idealized representation of a simple fluid in a pore Surface field v or between finely divided colloidal particles. The choice 0 of the model entails the assumption that the intermolec- FIG. 1. The critical wetting lines for semi-infinite systems 3 ular forces are short-ranged in character, there are no at T = 2.1: the corner wetting (hcw(T)) on the two planar .28 dleinspteorfsitohne-lbikuelkfomracegsn.etIicnfifleuldidHexdpeesrcirmibeenststhtehdeeevqiautiivoan- spularnfaacressuforframcein.gWaistthrariegshptecatngtoletahnedpwlaentatirngwe(thtwin(gT)h)1o=nt0h.1e 0 and 0.275 values correspond to the dry regime, whereas the of the critical chemical potential H ∼µ−µ that is de- 1 0 h1 =0.371 and 0.5 valueto thewet one. terminedbythe densityofthe fluidinthereservoir. The 2 1 phenomenonequivalenttothecapillarycondensationcan Wetting occurs in systems close to the line of phase : be studied in magnetic systems5. Generally, in the bulk, v coexistencewhenone phasemayadsorbpreferentiallyat the phase coexistence occurs for temperatures T < T i c asolidsubstrate. Typically,wemodelthesubstrateusing X and for vanishing bulk magnetic field H. In a slit L×∞ a planar surface where the critical line h (T), presented ar with identical surface fields at the boundaries, the com- in Fig. (1), is known exactly12 w bined effect of surface fields and confinement shifts the phase coexistence to a non-zero value of the bulk magnetic field H =H (L), which for large L scales as coe exp(2J/k T)(cosh(2J/k T)−cosh(2h J/k T)) B B w B =sinh(2J/k T). (2) B σ cosθ 1 0 H (L)= , (1) coe For this semi-infinite system if we approach the co- m L b existence line from the gas side along a given isotherm where σ , m , and θ are the surface tension of the free (T > T > T ) the amount of liquid adsorbed on the 0 b c w up-spin/down-spin interface, bulk spontaneous magne- surface l diverges l ∼ H−βsco. Moreover, the phase tran- tization, and contact angle given by Young’s equation, sition, called the complete wetting, is characterized by respectively6. In our case the liquid phase is represented the presence of the singular part of the excess surface byup-spins,gasphasebydown-spins,whereasareawith- free-energy fsing = H2−αcso. For the two-dimensional 2 Ising model, the values of the critical exponents are boundary fields with the following Hamiltonian: 2−αco =2/3 and βco =1/3. s s Consequently, when we consider the (pseudo)two-   dimensional Ising system in a slit geometry and the in- H=−J S S −h S −H S , X i,j k,ℓ 1 X i,j X i,j teraction between the liquid and walls is strong (above   ijkℓ surface all  the h (T) line), the Kelvin equation fails and some cor- spins spins w (4) rections are necessary. The reason is that, although the withJ >0andS =±1. Thefirstsumistakenoverthe i,j macroscopicallythick layerforms only for a semi-infinite nearest neighbours, while the second sum is performed system,anoticeableliquidlayerintervenesbetweenagas on spins at the surface only. The surface field h corre- 1 and the wall for a finite-size system as well. Therefore, sponding to direct short-range interactions between the Foster pointed out8 that if adsorbed layers were formed walls and spins is related to the preferential adsorption prior to condensation, the slit width L in Eq. (1) should on the surface for one of the two phases. The uniform be correctedby the layerthickness l. Derjaguinshowed9 bulk magnetic field H acts over all spins. that if solid-fluid forces decayed exponentially or had a finiterange,theporewidthLcouldbereplacedbyL−2l. Next, Evans et al.4 showed that the effects of wetting layerswere of quantitative rather than ofqualitative im- portance for capillary condensation. Albano et al.10 and Parry and Evans11 analyzed the next-order correction 2.2 termto the Kelvinequationfor the semi-infinite system. e ur Both studies, using scaling and thermodynamics argu- at er ments, concluded that for temperatures below the wet- mp 2 tingtemperatureTw (the dryregime)theleadingcorrec- Te tion to scalingterm was of type L−2. Above the wetting CTMRG DMRG (the wet regime) the correction is expected to be non- 1.8 Tc analytic due to a singularity of the surface free energy. For the two-dimensional Ising model10,11 the predicted -0.004 -0.003 -0.002 -0.001 0 correction term is proportional to L−5/3. Bulk field Inthesubsequentnumericalinvestigation,thedensity- FIG. 2. Pseudo-coexistence lines for the square (triangles) matrix renormalization techniques were employed13,14. and slit (circles) systems for L = 500 and h1 = 0.8. L is Foralargerangeofsurfacefieldsandtemperaturehigher- measuredinunitsoflatticeconstant,whereash1isinJ units. order corrections were not compatible with L−5/3, but they were of type L−4/3. It has been shown that this The origin of applied numerical method, called apparent disagreement was due to the fact that even for the Corner Transfer Matrix Renormalization Group thelargesizesconsidered(L∼150)thewettinglayerhas (CTMRG), came from Baxter20. Next Nishino and a limited thickness, so that the singular part of the sur- Okunishi21 combined his corner transfer matrix method facefreeenergythatdeterminesthecorrection-to-scaling with the ideas from the density matrix renormalization behavior is dominated by the contacts with the walls. group method (DMRG) approach. The last technique In real systems, where properties of both pure fluids wasdevelopedbyWhite22,23forthestudyofgroundstate and fluid mixtures confined to nanoporous and microp- properties of quantum spin chains and next extended orousmaterialsareunderconsideration15,thesquare-like by Nishino to two-dimensional classical systems in a slit geometryismorecommonthantheslitone. Inthiscase, geometry24. thegeometryofthesystemsignificantlyaffectsthecourse The generalideais to find arepresentationofthe con- of wetting phenomena because close to a corner the im- figurational space in a restricted space which is much pactoftheindividualwallsisstrong,whichshouldleadto smaller than the original one m≪2L2. This truncation moreintensiveformationofthewettinglayer. Therefore, is done through the construction of a reduced density the corner wetting transition should also be taken into matrix whose eigenstates provide the optimal basis set account16,17. Astwosidesofthesquareformthestraight m. Of course, the larger m, the better accuracy, so in angle, the corresponding (L → ∞) corner wetting line the present case we keep this parameter up to m=400. h (T), presented in Fig. (1), is known exactly18,19 cw Although in the original CTMRG algorithm the full transfermatrixis never constructed,wehavemodified it to determine the two eigenvectors related to the largest cosh(2J/kBT)−exp(−2J/kBT)sinh2(2J/kBT) eigenvalues. Because each of these vectors dominates on =cosh(2Jh /k T). (3) the opposite side of the coexistence line, using both vec- cw B tors for the construction of the density matrix guaran- In order to model the influence of wetting phenomena tees that the Hamiltonian is properly projected on the on the capillary condensation in a square geometry we subspace of most probable states. Our results have not consider a square Ising ferromagnet subject to identical shown any ambiguities of the calculated free energy and 3 by increasing the number of states kept m, our results Furthermore,because for a weak surfacefield (below the (the free energy) converged. To have a point of refer- wetting line) there remains a thin liquid layer between ence, we compared the results for the square geometry a gas and the wall, the gas-wall interface formally never with the results for the slit geometry, where the DMRG occurs. For the same reason, below the wetting line the technique was applied13,14. isolateddropletsofliquidcannotbeobservedonthewall, Due to the finiteness of L and to the nonvanishing so the contact angle cannot be drawn. surface field h the (pseudo)coexistence lines (for both As a useful tool for the analysis of the higher terms of 1 geometries) are shifted with respect to the bulk coexis- theKelvinequation,weintroducethelogarithmicderiva- tenceline(H =0). AsonecanseeinFig.(2)thiseffectis tive which acts as an effective dominant exponent muchstrongerforthe square-likesystem. This is easyto ln[H (L+∆L)]−ln[H (L)] understand, if one remembers that the influence of walls α(L)=− coe coe . (5) on the individual spins is here reinforcedwith respect to ln(L+∆L)−ln(L) the slit geometry. When the following expansion of the Kelvin equation is A schematic draw of the fluid shape adopted in a assumed square geometry below (h = 0.1) and above (h = 0.6) 1 1 wetting is demonstrated in Fig. (3). Both upper magnetization profiles correspond to the phase when a liquid A B H = + , (6) fills the square-like pore. For a weak field the magneti- coe Lα Lγ zation is lower nearby walls, whereas for a strong field combining both formulae, we obtain the first-order ex- the magnetization is essentially uniformly high through- pansion for the effective exponent out the square. The bottom magnetizations illustrate the situation on the other side of the coexistence line, when the gas phase occupies the middle of the square. B 1 α(L)=α+(γ−α) . (7) While the weak-field magnetization profile only slightly ALγ−α increases at the square edges, the strong-field magnetization profile increases considerably, which corresponds TheKelvinequationisexpectedtobevalidtothefirst to creation of the liquid layer. Note that the surface im- order for all T < Tc, but since only a limited size L is pact is always enhanced in the square corners, although available for the numerical computation, it is preferable limited to a relatively small area. to consider only temperatures not too close to the bulk critical temperature, where the scaling of the capillary critical point H (L) ∼ L15/8 is present6. Therefore, crit our results refer to temperature T = 2.1, but we have checked that they do not qualitatively change for other temperatures. 1.25 1.2 1.15 L) α( 1.1 hh11 == 00..530701 ((ssqquuaarree)) h1 = 0.275 (square) h1 = 0.100 (square) 1.05 h1 = 0.500 (slit) h1 = 0.371 (slit) h1 = 0.275 (slit) 1 h1 = 0.100 (slit) FIG.3. Two-dimensional magnetization profiles hsi,ji on the 0 0.02 0.04 0.06 0.08 0.1 0.12 square geometry 31 ×31 for the two values of the surface 1/L field h1 = 0.1 (left graphs) and h1 = 0.6 (right graphs) calculated at fixed temperature T = 2.0. The profiles change FIG. 4. Plots of the local exponent for both geometries at theirshapesandspinpolarizationswhilevaryingthebulkfield T = 2.1. Lines correspond to the square geometry, whereas value H below and above the wetting temperature Tw(h1 = thesymbols are related tothe slit one. 0.1)=2.257 and Tw(h1 =0.6)=1.814, respectively. Figure (4) shows that for a wide range of the surface The real phase transition (the complete wetting) oc- fields, when the system grows, the value α(L) goes to curs only in a semi-infinite system. Since we deal with α = 1. It confirms that to the accuracy of the first- a finite-size system, the thickness of the wetting layer is order expansion, the shift of the phase coexistence and limited, and a sharp liquid-gas interface is not observed. the system size are reciprocal in both geometries. 4 For the wet phase (θ = 0) in the slit geometry, the regime with respect to both wetting curves (see Fig. 1). exact formula for the coefficient is known (see Eq.(1)) Althoughthepresentedcurvescorrespondtoinfinitesys- giving A(T = 2.1) = 0.335 which very well agrees tems, we are dealing with systems large enough to guar- with our estimated value 0.345. Generally, our numer- antee that both points are localised on the appropriate ical fittings for curves H (L) according to Eq.(6) show side of the wetting curves. coe that the ratio of the main prefactors for the slit and To sum up, the Corner Transfer Matrix Renormal- square geometries is A /A ∼ 1/2 for each value ization Group method has been extended to study the slit square of the surface field. For example, below the wetting equilibrium statistical mechanics of simple fluids con- A (h = 0.1) = 0.121 and A (h = 0.1) = 0.242, fined in the square geometry providing some input for slit 1 square 1 whereas above the wetting A (h = 0.5) = 0.345 and understanding experiments in porous materials. As the slit 1 A (h =0.5)=0.669. This fact canbe explainedin method is not perturbative, it can be applied to arbi- square 1 the following way: when the size of the square increases, traryvaluesofthemodelparametersandyieldsaccurate the impact ofthe cornersremains moreor less the same, results, provided that a convergence of the free energy but the influence of the sides increases. Thus, when L is reached. Accuracy of the CTMRG method is com- goes to infinity, the system begins to resemble a set of pletely controlled by varying the number of the states two slits perpendicular to each other. kept. Moreover,the method does notsuffer fromgetting stuck in a local minimum of the free energy instead of the true global free energy minimum what is a typical problem in Monte Carlo simulations. There are neither 1.2 α = 1 + 5.2 L-2/3 - 8.6 L-1 metastabilities nor hysteresis effects provided that m is settobesufficientlylargewhichisourcaseconcludedby 1.15 α = 1 + 7 L-2/3 - 18.7 L-1 the convergedfree energy. h = 0.5 (slit) Our results confirm that in a more realistic (square) 1 L) h = 0.5 (square) geometry the coexistence line shift is inversely propor- α( 1.1 h11 = 0.1 (slit) tional to the system size L and the main pre-factor is h = 0.1 (square) 1 two times larger than for the case of an infinitely long 1.05 α = 1 + 2.7 L-1 - 15.5 L-2 α = 1 - 2.4 L-1 + 228 L-2 slit of width L. Moreover, for the square geometry the leading correction term to the Kelvin equation becomes negative in the dry regime. 1 0 0.005 0.01 0.015 0.02 We have also found that similarly to the slit geometry 1/L the leadingcorrectionsto scaling inthe squaregeometry areof type 1/L2 in the dry regimeand oftype 1/L5/3 in FIG.5. Theeffectiveexponentsα(L) andtheirfittingcurves thewetregime. Thelastcorrectionisagainnon-analytic denotedbythethicklinesfortwosamplevaluesofthesurface field: h1 = 0.1 for the dry system and h1 = 0.5 for the wet due to a singularity of the surface free energy. system (see Fig. (1)). Curve fitting was carriedout for the values of the sur- facefieldsthatcorrespondtothedry/wetregimewithre- spect to both wetting lines defined by Eqs.2 and 3. The AsonecanseeinFigs.4and5,theshapeofthecurves area between the wetting lines was not analysed in de- indicatesthatwearedealingwithvariousexpressionsofa tail, althoughthe strongestcompetitioncanbe expected higherorder. Firstofall,the B coefficientthatis always between the two types of wetting. This requires an ex- positiveintheslitgeometry,isnegativeinthedryregime amination of larger systems, which implies the need for for the square geometry. This is manifested by the α(L) a more precise calculation. In our next studies, we aim function minima for small 1/L, where the leading correction B/Lγ dominates. Obviously, the next correction to increase the precision of the numerical results by ad- ditional improvement of the CTMRG algorithm. has to be positive. Moreover, in the wet regime, where the B coefficients are positive, as one cansee in Fig.(4), thecurvescorrespondingtothesamevalueofthesurface field start to overlap when L becomes enough large. I. ACKNOWLEDGEMENTS Inthesquaregeometry,asfarastheleadingcorrection termsareconcerned,thepreciseanalysisofournumerical results shows, in both the dry and wet regimes, they are This workwasdone under projects VEGA-2/0074/12, the same as in the slit geometry (see Fig. 5). APVV-0646010 (COQI) and POKL.04.01.01-00-041/09- It is worth adding that the value of the surface field 00. Numerical calculations were performed in WCSS h = 0.1 (h = 0.5) was chosen to be in the dry (wet) Wrocl aw (Poland, grant 82). 1 1 1 J.N.Israelachvili,Intermolecular and surface forces, (Aca- 2 Handbuch of Nanophysics. Functional Nanomaterials, demic Press, New York,2nd ed., 1992). ed.KlausD.Sattler,(CRCPress,Taylor&FrancisGroup 5 2011). 16 A. O. Parry, A. J. Wood, and C. Rascón, J. Phys.: Con- 3 Lev D. Gelb, K. E. Gubbins, R. Radhakrishnan, and dens.Matter12,7671(2000);M.J.Greenall,A.O.Parry, M. Sliwinska-Bartkowiak, Rep. Prog. Phys. 62, 1573 andJ.M.Romero-Enrique,J.Phys.: Condens.Matter16, (1999). 2515 (2004). 4 R. Evans, U. Marini Bettolo Marconi, and P. Tarazona, 17 A. Milchev, M. Mueller, K. Binder, and D. P. Landau, J. Chem. Phys.84, 2376 (1986). Phys. Rev. E 68, 031601 (2003); A. Milchev, A. De Vir- 5 Kurt Binder and David Landau, J. Chem. Phys. 96, 1444 giliis, and K Binder, J. Phys.: Condens. Matter 17, 6783 (1992); Kurt Binder, David Landau and Marcus Mu¨ller, (2005). J. Stat.Phys. 110, 1411 (2003). 18 A. O. Parry, A. J. Wood, E. Carlon, and A. Drzewinski, 6 M. E. Fisher and H. Nakanishi, J. Chem. Phys. 75, 5857 Phys. Rev.Lett. 87, 196103 (2001). (1981);H.NakanishiandM.E.Fisher,J.Chem.Phys.78, 19 D. B. Abraham, and A. Maciolek, Phys. Rev. Lett. 89, 3279 (1983). 286101 (2002). 7 W. Thomson, Philos. Mag. 42, 448 (1871). 20 R.J. Baxter, J. Stat.Phys. 19, 461 (1978). 8 A. G. Foster, Trans. Faraday Soc. 28, 645 (1932). 21 T. Nishino and K. Okunishi, J. Phys. Soc. Jpn. 65, 891 9 B. V.Derjaguin, Zh.Fiz. Khim.14, 137 (1940). (1996); K. Ueda, R.Krˇcmár, A. Gendiar, and T. Nishino, 10 E. V.Albano, K. Binder, D. W. Heermann, and W. Paul, J. Phys. Soc. Jpn. 76, 084004 (2007). J. Chem. Phys.91, 3700 (1989). 22 S.R.White,Phys.Rev.Lett.69,2863(1992);S.R.White, 11 A. O.Parry and R. Evans, J. Phys. A 25, 275 (1992). Phys. Rev.B 48, 10345 (1993). 12 D. B. Abraham, Phys.Rev. Lett.44, 1165 (1980). 23 U. Schollwoeck, Rev. Mod. Phys. 77, 259 (2005); K. 13 E. Carlon, A. Drzewin´ski, and J. Rogiers, Phys. Rev. B Hallberg, Adv. Phys. 55, 477 (2006); U. Schollwoeck, 58, 5070 (1998). Ann.Phys.(NY) 326, 96 (2011). 14 A. Drzewin´ski, and K. Szota, Phys. Rev. E 71, 056110 24 T. Nishino, J. Phys. Soc. Jpn. 64, 3598 (1995). (2005). 15 I.BrovchenkoandA.Oleinikova,Interfacial and Confined Water, (Elsevier, Amsterdam, 2008).