A finite-size scaling study of wedge filling transitions in the 3D Ising model

Depósito de investigación de la Universidad de Sevilla https://idus.us.es/ “This is an Accepted Manuscript of an article published in Soft Matter on May 2013 available at https://doi.org/10.1039/C3SM50207D.” A finite-size scaling study of wedge filling transitions in the 3D Ising model Jos´ e Manuel Romero-Enrique,∗a,´ Alvaro Rodr´ ıguez-Rivas a, Luis F. Rull aand Andrew O. Parry b Received Xth XXXXXXXXXX 20XX, Accepted Xth XXXXXXXXX 20XX First published on the web Xth XXXXXXXXXX 200X DOI: 10.1039/b000000x We propose a fully quantitative theory for the finite-size scaling of the filling transition in a three dimensional double wedge geometry, based on the exact transfer-matrix solution of a phenomenological interfacial model. Antisymmetric fields act at the top and bottom wedges, so each one favours a different bulk phase under coexistence conditions, i.e. gas and liquid phases in fluid models such as the lattice gas, or equivalently ferromagnetic domains of opposed magnetization in the Ising model. From this formalism we obtain an analytical form for the magnetization probability distribution function at critical filling which is valid for any aspect ratio. To test our predictions we revisit and perform new simulation studies of filling in the Ising model double-wedge geometry and use our finite-size scaling theory to locate accurately the critical filling transition. 1 Introduction Fluid adsorbed on micropatterned and sculpted solid substrates are known to exhibit phase transitions which differ from those observed at planar, homogeneous walls1–3. The simple 3D wedge geometry has been extensively studied in the past decade theoretically4–17, experimentally3,18–20 and by computer simulation21–23. It is well established by thermodynamic arguments24–26 that the wedge is completely filled with liquid provided that the contact angle θ is less than the tilt angle α . This transition may be first-order or continuous (critical filling). For the latter, the characteristic lengthscales are the averaged interfacial height ℓW, the correlation length ξ yalong the wedge axis and the interfacial roughness ξ ⊥. In an infinite wedge these all diverge continuously as θ → α according to the power-laws ℓw≡∼( θ − α )− β W, ξ y∼( θ − α )− ν y, ξ ⊥∼( θ − α )− ζ W ν y(1) where, for systems with short-ranged forces, the critical exponents take the values β W=1/4, ν y=3/4 and ζ W=1/3. The last exponent identifies the interfacial wandering exponent for 3D critical filling which is much larger than for 3D wetting transitions at planar substrates. This arises because close to critical filling the dominant fluctuations, which correspond to local translations in the position of the interface (breather mode excitations), are strongly anistropic owing to the substrate geometry. This effective dimensional reduction aDepartamento de F´ ısica At´ omica, Molecular y Nuclear, ´ Area de F´ ısica Te´ orica, Universidad de Sevilla, Apartado de Correos 1065, 41080 Sevilla (Spain) bDepartment of Mathematics, Imperial College London,180 Queen’s Gate, London SW7 2BZ, United Kingdom strongly enhances interfacial fluctuation effects5,6. In addition the conditions for critical or first-order filling can be quite different to critical or first-order wetting. In a shallow wedge, corresponding to small α , the conditions for observing continuous wedge filling transition are less restrictive for small α than for critical wetting at planar walls, and critical filling may occur even if the walls of the wedge show first-order wetting5,6. In more acute wedges the situation is more complicated and the order of the filling transition may be influenced by a number of factors at mean-field level and beyond. For example recently it has been proposed that if the wedge is acute enough, the filling transition may become first-order even if the wetting at the planar walls is continuous; an effect due to the self-interaction of the interface16,17. Furthermore, a phenomenological theory which takes into account the breather mode interfacial fluctuations shows that the filling transition may be driven first-order by modifying the fluid-solid interactions close to the wedge bottom13–15. Critical filling transitions have been observed unambiguously in computer simulation studies of the 3D Ising model21–23. These authors considered a special geometry (the double wedge) in which they were able to make a systematic finite-size characterization of the critical behaviour of the filling transition, confirming the theoretical predictions for the critical exponents (1) obtained from interfacial Hamiltonian models5,6. In this paper we generalize the phenomenological interfacial model of wedge filling, previously applied to the infinite wedge geometry, to the finite-size double-wedge geometry. The importance of this is threefold. Firstly we want to test the phenomenological theory with existing computer simulations. Secondly, if the theory is accurate, it will provide us with a means of accurately locating the critical filling transition. Fi- 1–8 |1 nally, the theory may guide us to determing with greater accuracy when the filling transition becomes first-order by, for example, changing the tilt angle α or modifying the fluid-solid interactions close to the wedge bottom. Our paper is organized as follows. In Section 2 we revisit the breather-mode interfacial model for wedge filling, and generalize this to the finitesize double-wedge geometry. Section 3 is devoted to the Ising model computer simulations, and the comparison between the simulation results and the phenomenological theory. Finally, in Section 4 we will present our conclusions. 2 The phenomenological model In order to propose a phenomenological theory for the fluid adsorption in the double wedge, we first review the breathermode model of filling in an infinite wedge. A mean-field analysis of the standard interfacial Hamiltonian model applied to the wedge geometry shows that near bulk coexistence the interface is locally flat across the wedge and the fluctuation effects are dominated by pseudo-one-dimensional local translations in the height of the filled region along it5,6. This suggests that filling phenomenona can be understood using an effective pseudo-one-dimensional wedge Hamiltonian which accounts only for such breather-mode excitations5,6: HW[ℓ] = Zdy(Λ(ℓ) 2dℓ dy2 +VW(ℓ))(2) where ℓ(y)>0 is the local height of the interface above the wedge bottom. The effective bending term Λ(ℓ)resisting fluctuations along the wedge can be expressed as5,6 Λ(ℓ)≈2Σℓ tan α (3) where Σis the interfacial tension (more properly the stiffness for lattice models). Similarly for large ℓthe effective wedge binding potential VW(ℓ)takes the form5,6: VW(ℓ)≈hℓ2 tan α +2Σ(cos α −cos θ )ℓ sin α +AF (p−1)ℓ1−p+... (4) where his the bulk ordering field (h=0 at vapour-liquid coexistence) and θ is the contact angle of a liquid drop at the planar wall-vapour interface. The first two terms reflect the bulk and surface thermodynamic costs of forming a filled region7 and are written here as appropriate for general tilt angle generalizing earlier expressions reported for a shallow wedge5,6. The last contribution arises from the fluid-solid interactions so that AFis the effective Hamaker constant associated with the attractive part of the wall-fluid interface binding potential. The exponent pspecifies the range of the attractive part of the intermolecular forces. For three-dimensional wedges and for (non-retarded) van der Waals forces the appropriate value is p=2. For systems with short-ranged forces, such as the Ising model, this term is exponentially small, of order exp(− κ ℓ), where κ −1≡ ξ bis the bulk correlation length. Unlike 3D critical wetting where the form of this exponential interaction is crucial owing to the marginal nature of the phase transition, the strong nature of fluctuation efects at 3D wedge filling mean that we can further simply and model this as a contact potential. This can be considered to be an effective line tension associated with the wedge bottom13–15. The quasi-one-dimensional character of the Hamiltonian (2) means it is ameniable to an analytic solution. The partition function corresponds to path integral27: Z(ℓb,ℓa,Ly) = ZDℓexp(−HW[ℓ]) (5) where we set kBT=1 for convenience and ℓaand ℓbare the (fixed) values of the interfacial height at the end points aand b, separated a distance Lyalong the wedge axis. The presence of a position-dependent stiffness coefficient makes that some care is required in the definition of the partition function and its measure8. This issue is similar to the ”factor-ordering” problem in solid-state quantum mechanics when there is an effective position-dependent mass28,29. We follow the prescription proposed by us in earlier work which provides a solution which is mathematical consistent and agrees with necssary thermodynamic requirements13–15. There it is shown that the partition function Z(ℓb,ℓa,Ly)verifies the differential equation: HWZ(ℓb,ℓa,Ly) = − ∂ Z(ℓb,ℓa,Ly) ∂ Ly(6) with the initial condition Z(ℓb,ℓa,Ly)→ δ (ℓb−ℓa)as Ly→0. The operator HWis defined as29: HW≡−1 2 ∂ ∂ ℓb1 Λ(ℓb) ∂ ∂ ℓb+VW(ℓb)+ ˜ VW(ℓb)(7) where ˜ VW(ℓ)is given by ˜ VW(ℓ) = −1 2Λ(ℓ)"3 4Λ′(ℓ) Λ(ℓ)2 −Λ′′(ℓ) 2Λ(ℓ)#(8) and the prime denotes the derivative with respect to its argument. Substituting Eq. (3) into Eq. (8), we obtain ˜ VW(ℓ) = −3tan α 16Σℓ3(9) The solution of the equation (6) can be expressed as27: Z(ℓb,ℓa,Ly) = ∑ αψα (ℓb) ψ ∗ α (ℓa)e−E α Ly(10) 2|1–8 α l(y) x y z L Ly Fig. 1 Schematic illustration of a typical interfacial configuration in the 3D double-wedge geometry. where { ψα (ℓ)}is a complete orthonormal set of eigenfunctions of the Hamiltonian operator HW, with associated eigenvalues E α . In the case of a continuous spectrum, the sum is replaced by an integral. From this result we can obtain all relevant structural properties. For example in the thermodynamic limit of an infinitely long wedge, the probability distribution function (PDF), PW(ℓ), for the interface height at position Y, follows as PW(ℓ) = lim Ly→∞ Z(ℓ,ℓ−Ly/2,Y+Ly/2)Z(ℓLy/2,ℓ,Ly/2−Y) Z(ℓLy/2,ℓ−Ly/2,Ly)(11) and reduces simply to PW(ℓ) = | ψ 0(ℓ)|2, and is independent of the value of Yconsidered. For the case of short-ranged forces, the contact potential leads to a boundary condition for Zat the wedge bottom, ℓ=0, which is controlled by the potential (line tension) strength −U13–15. Along the filling transition boundary h=0, θ = α , two different situations may arise depending on the value of U13–15. Above some value Utc, the operator HWhas an isolated bound ground state together with a continuous spectrum with E>0. This means that the interface remains bound to the wedge bottom as θ → α +implying that the filling transition is first-order. Below Utc, there is only a continuous spectrum with E>0 and no bound states, implying the interface is delocalized and that the filling transition, occuring as θ → α +, is critical. Under these conditions, the PDF shows a short-distance expansion PDF ∝ℓ3, consistent with scaling arguments7. This is the necessary dependence that ensures thermodynamic consistency with exact sum-rules for filling fluctuation regime, which dictate that the local density at the wedge bottom ρ W(0)is non-singular, i.e. ρ W(0)− ρ l∼T−Tf, where ρ lis the bulk liquid density, and Tfis the filling transition temperature30. The borderline situation U=Utc corresponds to a tricritical point13–15. Now we turn to the finite-size geometry. As we want to compare with our phenomenological theory with Ising model simulations we will consider the same geometry as in the Monte Carlo simulations used for the latter. Thus, the fluid will be confined in a double-wedge geometry pore, as shown in Fig. 1. The fluid is confined in a one-dimensional pore of lateral length Ly, with a rhombic cross-section formed by a double wedge characterized by a tilt angle α and a side length L. Periodic boundary conditions are imposed along the wedge axis, and antisymmetric (opposing) surface fields are applied at the upper and lower confining surfaces, favouring gas and liquid adsorption respectively. Assuming that, as for the infinite wedge, the relevant interfacial fluctuations are breather modes, and that the interface is locally flat across the wedge, the interfacial height above the wedge bottom can be regarded as a fuction of just one co-ordinate: ℓ(y)with 0< ℓ(y)<2Lsin α . Consequently the effective interfacial Hamiltonian is still given by Eq. 2, but the effective bending coefficient is now given by: Λ(ℓ) =      2Σℓ tan α 0< ℓ < Lsin α 2Σ(2Lsin α −ℓ) tan α Lsin α < ℓ < 2Lsin α (12) Similarly for systems with short-ranged forces the only longranged part-of the effective wedge binding potential VW(ℓ) arises from the surface energy cost of forming a interfacial configuration, and at bulk coexistence, h=0, we find VW(ℓ) =      2Σ(cos α −cos θ )ℓ sin α ;0 < ℓ < Lsin α 2Σ(cos α −cos θ )(2Lsin α −ℓ) sin α ;Lsin α < ℓ < 2Lsin α (13) Analogous to the case of the infinite wedge, there are contact potentials with equal strengths at the wedge top ℓ=2Lsin α and bottom ℓ=0. For the case of periodic boundary conditions the interfacial partition function in the double wedge can be obtained as Zperiodic(L,Ly) = Z2Lsin α 0Z(ℓ,ℓ,Ly)dℓ(14) where Z(ℓa,ℓb,Ly)satisfies the differential equation Eq. (6) with Λ(ℓ)and VW(ℓ)given by Eqs. (12) and (13). Since the spectrum of HWis discrete, we can expand Zperiodic(L,Ly)in terms of the set of eigenvalues {E α }with associated eigenfunctions ψα , satisfying HW ψα =E αψα . Using the spectral expansion Eq. (10) and the orthonormal character of the set of eigenfunctions of HW, we arrive at: Zperiodic =∑ α e−E α Ly(15) From this we can define the PDF for the interfacial height as PW(ℓ,L,Ly) = Z(ℓ,ℓ,Ly) Zperiodic(L,Ly)(16) 1–8 |3 which reduces to PW(ℓ,L,Ly) = ∑ α | ψα (ℓ)|2e−E α Ly ∑ α e−E α Ly(17) In order to proceed, we now need to solve the eigenvalue problem HW ψ (ℓ) = E ψ (ℓ)and determine the spectrum. As the operator HWis symmetric with respect to ℓ=Lsin α , the corresponding eigenfunctions will be either symmetric or antisymmetric. Consequently we only need to solve the eigenvalue problem in the interval 0 < ℓ < Lsin α and impose the appropriate symmetry at the mid-point. At ℓ=0, we force the eigenfunctions to have the same short-distance expansion as in the infinite wedge situation15 lim ℓ→0+ √ℓ ψ (ℓ) d( ψ (ℓ)/√ℓ) dℓ=±Γ[−1/3]3−2/3 Γ[1/3] ξ −1 u(18) where ξ u∼ |U−Utc|−1is a characteristic length associated with the contact potential at the wedge bottom and the positive and negative signs corresponds to U>Utc (first-order filling) and U<Utc (critical filling), respectively15. Thus, the eigenfunctions have a short-distance expansion: ψ (ℓ)∼√ℓ±Γ[−1/3]3−2/3 Γ[1/3] ℓ3/2 ξ u(19) For critical filling, ξ u→0 , so ψ (ℓ)∼ℓ3/2. At the mid-point height ℓ=Lsin α , we have to take into account that Λ(ℓ)given by Eq. (12) has a cusp singularity so that Λ′(ℓ) = ∓2Σ/sin α as ℓ→Lsin α ±. Equivalently there is Dirac delta contribution to Λ′′(ℓ). Consequently, we can integrate the eigenvalue equation between ℓ=Lsin α −and ℓ=Lsin α +, and obtain ∂ψ ∂ ℓℓ=Lsin α +− ∂ψ ∂ ℓℓ=Lsin α − + ψ (Lsin α ) Lsin α =0 (20) We are now in a position to use the symmetry properties of the eigenfunctions. For odd eigenfunctions, ( ∂ψ / ∂ ℓ)at ℓ=Lsin α +is equal to ( ∂ψ / ∂ ℓ)at ℓ=Lsin α −. Consequently Eq. (20) reduces to ψ (Lsin α ) = 0. On the other hand, for even eigenfunctions, ( ∂ψ / ∂ ℓ)at ℓ=Lsin α +is equal to −( ∂ψ / ∂ ℓ)at ℓ=Lsin α −, in which case the boundary condition at ℓ=Lsin α simplifies to: ∂ψ ∂ ℓℓ=Lsin α − = ψ (Lsin α ) 2Lsin α (21) To proceed, the general solution of the eigenvalue problem HW ψ =E ψ for θ > α in the range 0 < ℓ < Lsin α is ψ (ℓ)∝√ℓeℓ ξθ ˜ ε −ℓ ξθ "H˜ ε 2 4−1 2√2ℓ ξθ −˜ ε √2 +C1F11 4−˜ ε 2 8,1 2,(˜ ε −2ℓ/ ξθ )2 2#(22) where Hs(x)is the s-order Hermite function and 1F1(a,b,x)is the Kummer confluent hypergeometric function31. Here the characteristic lengthscale ξθ and the reduced eigenvalue ˜ ε are defined as ξθ =√sin α [2Σ2cos α (cos α −cos θ )]1/4˜ ε =Σ ξ 3 θ tan α E(23) The constant Cand the reduced eigenvalues ˜ ε are obtained after imposing the boundary conditions Eqs. (18) and (20). Solutions (22) may be analytically continued to θ < α . Exactly at the filling transition, θ = α , the eigenfunctions take the form ψ (ℓ)∝√ℓhAi(−4 ε )1/3ℓ+CBi(−4 ε )1/3ℓi (24) where Ai(x)and Bi(x)are Airy functions, and the reduced eigenvalue ε =ΣE/tan α . As for the case θ > α , the constant Cand the reduced eigenvalues are obtained from the boundary conditions Eqs. (18) and (20). In particular, the boundary condition Eq. (18) reduces to: ± ξ −1 u= (−4 ε )1/31 √3−C 1 √3+C(25) Now consider the finite-size behavior of the fluid adsorption occurring exactly at a filling critical transition, i.e. θ = α and ξ u=0. In this case the boundary condition Eq. (25) reduces simply to the condition C=−1/√3. Application of the boundary condition Eq. (20) then leads to the following relationships for the values of x≡(−4 ε )1/3Lsin α : 1 √3=(Ai′(x) Bi′(x) ψ (ℓ)even Ai(x) Bi(x) ψ (ℓ)odd (26) These equations can be solved numerically or graphically (see Fig. 2), leading to an infinite series of solutions at xi, with i=0,1,2,.... Table 1 shows the values corresponding to the first few solutions. The corresponding eigenvalues Eican be expressed as Ei=−tan α x3 i/4Σ(Lsin α )3. We see that the ground state corresponds to an even eigenfunction, and that the parity of the eigenfunctions, ordered by the eigenvalue, is alternating even and odd. It follows that the PDF for the interfacial height in the region 0 < ℓ < Lsin α is given by PW(ℓ,L,Ly) = ∑∞ i=0NiℓhAi−xiℓ Lsin α −1 √3Bi−xiℓ Lsin α i2e κ x3 i ∑∞ i=0e κ x3 i (27) where Niare the normalization factors of the eigenfunctions and κ = (tan α /4Σsin3 α )(Ly/L3). Note that the PDF exactly at critical filling does not depend on Land Lyindependently, 4|1–8 -10 -5 05 x -3 -2 -1 0 1 2 3 Fig. 2 Graphical solution of the Eq. (26). The continuous line curve corresponds to the function Ai′(x)/Bi′(x)−1/√3, and the dashed line one to Ai(x)/Bi(x)−1/√3. The values of xiare the zeroes of these functions. Table 1 First solutions of Eq. (26). i xiEigenfunction parity 0 -1.5149 Even 1 -2.6664 Odd 2 -3.5341 Even 3 -4.3425 Odd 4 -5.0562 Even 5 -5.7410 Odd 6 -6.3726 Even 7 -6.9861 Odd 8 -7.5639 Even but via the scaling combination Ly/L3, i.e. PW(ℓ,L,Ly) = PW(ℓ,Ly/L3), in agreement with the assumptions of previous scaling arguments21–23. In order to compare these predictions with the existing Ising model computer simulation results, we must covert the dependence on the interfacial height into an appropriate microscopic observable. For Ising model this is the local magnetization in the vertical plane at position Yalong the wedge. In the breather-mode picture this change of variables is straightforward32. In the case of the double wedge we can write m mb=(ℓ Lsin α 2−1 ;0 < ℓ < Lsin α 1−2Lsin α −ℓ Lsin α 2;Lsin α < ℓ < 2Lsin α (28) where mb>0 is the bulk magnetization. That is, in standard lattice-gas language, the bulk “liquid” phase has an associated magnetization density +mb, while for the bulk “gas” the magnetization density is −mb. The magnetization PDF is then related to the interfacial height PDF via PW(m) = -0.5 00.5 m 0 0.2 0.4 0.6 0.8 1 PW(m) Ly/L3=0.0046 Ly/L3 infinity Fig. 3 Magnetization PDFs at the critical filling transition obtained from the phenomenological theory. The parameters considered are α = π /4, Σ=0.0981 and mb=0.75 to mimic the existing Ising model studies. The continuous curve corresponds to Ly/L3=0.0046, and the dashed line to Ly/L3→∞. PW(ℓ(m))|dℓ/dm|. Thus the PDF for the magnetization exactly at the critical filling transition is given by Pw(m,Ly/L3) = ∑∞ i=0˜ ψ 2(m) R+mb −mbdm ˜ ψ 2(m)e κ x3 i ∑∞ i=0e κ x3 i(29) where κ and xiare defined as above, and the functions ˜ ψ i(m) are given by: ˜ ψ i(m) = Ai −xis1−|m| mb!−1 √3Bi −xis1−|m| mb! (30) Fig.3 shows the magnetization PDFs predicted by this phenomenological theory. For Ly→∞, the PDF is unimodal, corresponding to the situation PW(m) = ˜ ψ 2 0(m). As Lydecreases, the magnetization PDF becomes bimodal, where the mostprobable magnetization density in each section depend on the value of the scaling variable Ly/L3, but in absolute value is always smaller than mb. 3 Computer simulation results In order to check the prediction of the phenomenological theory, we will resort to computer simulations of the Ising model. In particular, we revisit the Ising model simulations in the double-wedge geometry considered by Milchev et al.21–23. They considered Monte Carlo simulations of the nearest neighbour Ising model on a simple cubic lattice with linear dimensions L×L×Ly. Periodic boundary conditions are applied along the ydirection, and in the remaining bound- 1–8 |5 aries they considered free boundary conditions, but with applied surface fields +Hsin two of the surfaces which form a wedge (W1), and −Hsin the other ones (W2). The Hamiltonian of the Ising model in this geometry is β H=−1 4∑ hi,jibulk SiSj−1 8∑ hi,ji∈W1∪W2 SiSj − β Hs∑ i∈W1 Si+ β Hs∑ i∈W2 Si(31) where β =1/kBT, β J=1/4 and the surface exchange constant Js=J/221–23. We will take Jas the energy unit, so J=1. In this way, a double wedge is created, characterized by a tilt angle α = π /4. Under these conditions, the value of mb≈0.75, and Σ= βσ ≈0.0981, with σ being the interfacial tension of the Ising model22. The filling transition is found to be critical under these conditions. By using finite-size techniques, Milchev et al. estimated the value of Hs(which determines the contact angle θ ) at the filling transition to be Hs≈0.7221–23. These studies focus on properties related to the volume averaged magnetization density, while our theory predicts the behaviour of the magnetization density at every section perpendicular to the wedge axis. Thus we performed Ising model simulations under the same conditions, in order to test our phenomenological theory. These Monte Carlo simulations were performed by using the standard Metropolis algorithm33. As the evaluation of the quantities we are interested in are quite slow to equilibrate, we considered runs of order of 108sweeps, where a sweep is L2×Lyattempted updates of a spin chosen at random. In addition to the usual quantities such as the energy or total magnetization, during the simulation we also evaluated the magnetization associated with the L2spins in each slice perpendicular to the ydirection, which we will denote as ˜ Rj, where j=1,...,Ly. In order to minimize the effect of the enhanced order of the magnetization close to the boundaries, we excluded the spins on sites which are nearest or next-to- nearest neighbours to the walls. The set of (L−4)2sites which are in ˜ Rjbut are not either nearest or next-to-nearest neighbours to the walls will be denoted as Rj. The magnetization in each slice Rjis defined as M(j) = ∑i∈RjSi, with j=1,...,Ly, and the corresponding magnetization density m(j) = M(j)/(L−4)2. We also evaluated the surface magnetization at each slice, i.e. Ms(j) = ∑i∈Mj∩(W1∪W2)Si. The magnetization PDF can be obtained from the computer simulation as PW(m) = (L−4)2 2Ly*Ly ∑ j=1 δ M(j),(L−4)2m+(32) Here δ i,jis the Kronecker symbol and we have taken advantage of the translational symmetry along the yaxis due to the periodic boundary conditions, and h...icorresponds to the average. This quantity can be compared with the theoretical prediction from the phenomenological theory. We also considered the joint PDF for mand Ms, which can be obtained as: PW(m,¯ Ms) = (L−4)2 4Ly*Ly ∑ j=1 δ M(j),(L−4)2m δ Ms(j),¯ Ms+(33) Note that PW(m) = 2∑MsPW(m,Ms)≈RdMsPW(m,Ms). This quantity is particularly useful as it is possible to predict the magnetization PDF at a value of Hswhich is different from the value used in the simulation by using reweighting techniques34,35. If ¯ Hsis the simulation value of the surface magnetic field (the conjugate field to Ms), and PW(m,Ms;¯ Hs)is the joint PDF obtained from the simulation, we can estimate the joint PDF at a different surface field Hsas: PW(m,Ms;Hs)≈PW(m,Ms;¯ Hs)e− β (Hs−¯ Hs)Ms RdmdMsPW(m,Ms;¯ Hs)e− β (Hs−¯ Hs)Ms(34) This is practical only if Hsdoes not differ too much from the simulation value. We performed computer simulations for different box sizes: 19×19×16, 24×24×37 and 34×34×122. The simulation box aspect ratio Ly/(L−4)3is approximately the same for the three cases, and equal to 0.0046. The simulations are done for Hs=0.723, and by using the reweighting technique we tune the value of Hsin order to match with the simulations for each size. Ideally, at exactly the critical filling transition, the simulated magnetization PDF should be independent of the size, and the matching with the theoretical value would lead to the same set of magnetization moments (in particular, the Binder cumulant). The results are shown in Fig. 4. We see that the simulated magnetization PDFs show tails for |m|>mb. So, in order to make the aforementioned, we have to multiply the simulated PDFs by a constant factor (always around unity) in order to get good agreement. This factor will not affect the evaluation of quantities such as the Binder cumulant. In any case, the matching is reasonable for the most-probable magnetizations in the range of sizes considered. The apparent values of Hsat the transition are not exactly the same (Hs= 0.7355 for L=19, Hs=0.7286 for L=24 and Hs=0.7238 for L=34), indicating that corrections to scaling still play a role in the location of the critical filling value of Hs. 4 Discussion and conclusions In this paper we have generalized the breather-mode effective Hamiltonian theory of filling transitions to finite-size systems13–15 and have in particular, applied it to fluid adsorption in a double wedge geometry used in previous Ising 6|1–8 -1 -0.5 00.5 1 m 0 0.2 0.4 0.6 0.8 PW(m) 19x19x16 (x1.05) , H s= 0.7355 24x24x37 (x1.02) , H s =0.7286 34x34x122 (x1.01) , H s =0.7238 Phenomenological theory Fig. 4 Plot of the magnetization probability distribution functions for: Hs=0.7355, 19×19×16 (squares), Hs=0.7286, 24×24×37 (circles) and Hs=0.7238, 34×34×122 (triangles). Under these conditions an optimal match with the theoretical prediction (wide line) is obtained. The magnetization PDF was scaled by a factor 1.05, 1.02 and 1.01, respectively. model simulations21–23. We found that, at the critical filling phase boundary, the probability distribution function for the local magnetization density, in square section perpendicular to the wedge axis, has the scaling form PW(m,L,Ly) = PW(m,Ly/L3), in agreement with previous assumptions and provided an analytical expression for the scaling function. We have tested our theory by performing further Monte Carlo simulations of the Ising model under the same conditions as that reported earlier in the literature21–23, and have found reasonable quantitative agreement between theoretical predictions and simulation results for the precise form of the magnetization PDF. This study shows that the local magnetization density in each section is a better order parameter for the 3D wedge filling than the overall instantaneous magnetization, and it clearly confirms that the breather modes are the most relevant interfacial fluctuations in this phenomenon. Our work opens a number of questions which will be addressed in future studies. Although it has been established that the phenomenological model is an accurate description of the universal scaling properties for the adsorption in an Ising double wedge, it would be desirable to identify what are the most important corrections to scaling. A possible candidate would be the bulk magnetization fluctuations (responsible for the the magnetization PDF tails when for m>|mb|). Their characterization would provide a means of extrapolating the apparent finite-size critical filling parameters to the infinite wedge thermodynamic limit. In addition we hope that the present observed agreement between theory and simulation regarding critical filling will helps us explore recent theoretical predictions that the filling transition may be driven first-order. This can be studied by for example introducing suitable modifications of the spin-wall interactions close to the wedges. Another interesting problem would be to change the tilt angle to see if the filling transition becomes first-order in more acute wedges16,17. For these goals it is important to know as accurately as possible the phase boundary of the wedge filling transition. The procedure outlined in this paper is a step in that direction. We hope to address these problems in the future. Acknowledgments We aknowledge financial support from the Spanish Ministerio de Econom´ ıa y Competitividad through grants no. FIS2009- 09326 and FIS2012-32455, and Junta de Andaluc´ ıa through grant no. P09-FQM-4938, all co-funded by the EU FEDER. References 1 H. Gau, S. Herminghaus, P. Lenz and R. Lipowsky, Science, 1999, 283, 46. 2 C. Rasc´ on and A. O. Parry, Nature, 2000, 407, 986. 3 L. Bruschi, A. Carlin and G. Mistura, Phys. Rev. Lett., 2002, 89, 166101. 4 K. Rejmer, S. Dietrich and M. Napi´ orkowski, Phys. Rev. E, 1999, 60, 4027. 5 A. O. Parry, C. Rasc´ on and A. J. Wood, Phys. Rev. Lett., 2000, 85, 345. 6 A. O. Parry, A. J. Wood and C. Rasc´ on, J. Phys.: Condens. Matter, 2001, 13, 4591. 7 M. J. Greenall, A. O. Parry and J. M. Romero-Enrique, J. Phys.: Condens. Matter, 2004, 16, 2515. 8 A. Bednorz and M. Napi´ orkowski, J. Phys. A: Math. Gen., 2000, 33, L353. 9 J. R. Henderson, J. Chem. Phys., 2004, 120, 1535. 10 J. R. Henderson, Phys. Rev. E, 2004, 69, 061613. 11 J. R. Henderson, Mol. Simul., 2005, 31, 435. 12 C. Rasc´ on and A. O. Parry, Phys. Rev. Lett., 2005, 94, 096103. 13 J. M. Romero-Enrique and A. O. Parry, J. Phys.: Condens. Matter, 2005, 17, S3487. 14 J. M. Romero-Enrique and A. O. Parry, Europhys. Lett., 2005, 72, 1004. 15 J. M. Romero-Enrique and A. O. Parry, New J. Phys., 2007, 9, 167. 16 A. O. Parry and C. Rasc´ on, J. Phys.: Condens. Matter, 2011, 23, 015004. 17 N. R. Bernardino, A. O. Parry and J. M. Romero-Enrique, J. Phys.: Condens. Matter, 2012, 24, 182202. 18 L. Bruschi, A. Carlin and G. Mistura, J. Chem. Phys., 2001, 115, 6200. 19 L. Bruschi, E. Carlin and G. Mistura, J. Phys.: Condens. Matter, 2003, 15, S315. 20 L. Bruschi, A. Carlin, A. O. Parry and G. Mistura, Phys. Rev. E, 2003, 68, 021606. 21 A. Milchev, M. Muller, K. Binder and D. P. Landau, Phys. Rev. Lett., 2003, 90, 136101. 22 A. Milchev, M. Muller, K. Binder and D. P. Landau, Phys. Rev., 2003, 68, 031601. 23 K. Binder, M. Muller, A. Milchev and D. P. Landau, Comput. Phys. Commun., 2005, 169, 226. 24 P. Concus and R. Finn, Proc. Natl Acad. Sci. USA, 1969, 63, 292. 25 E. H. Hauge, J. Colloid Interface Sci., 1986, 113, 5. 26 Y. Pomeau, Phys. Rev. A, 1992, 46, 4994. 27 T. W. Burkhardt, Phys. Rev. B, 1989, 40, 6987. 1–8 |7 28 J. Thomsen, G. T. Einevoll and P. C. Hemmer, Phys. Rev. B, 1989, 39, 12783. 29 L. Chetouani, L. Dekar and T. F. Hammann, Phys. Rev. A, 1995, 52, 82. 30 A. O. Parry, M. J. Greenall and A. J. Wood, J. Phys.: Condens. Matter, 2002, 14, 1169. 31 N. N. Lebedev, Special Functions and their Applications, Dover, New York, 1972. 32 E. V. Albano, A. D. Virgiliis, M. M¨ uller and K. Binder, J. Phys.: Condens. Matter, 2003, 15, 333. 33 N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. N. Teller and E. Teller, J. Chem. Phys., 1953, 21, 1087. 34 A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett., 1988, 61, 2635. 35 A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett., 1989, 63, 1195. 8|1–8