Andrew O Parry

Nelson Rei Bernardino

Carlos Rascón

Scaling of correlation functions near capillary condensation

Keywords: Capillary Condensation, Wetting, Correlation Function, Effective Hamiltonian, Scaling

published or not. The documents may come

Scaling of correlation functions near capillary condensation

I. INTRODUCTION

Capillary condensation [1][START_REF] Derjaguin | Fluid Interfacial Phenomena[END_REF][START_REF] Schoen | [END_REF][4][5][6][7] refers to the shift in the value of the chemical potential µ cap (or pressure), relative to that of bulk saturation µ sat , at which a vapour confined in a slit of width L condenses to a liquid-like phase, see Fig. 1. As is well known, for macroscopically wide slits, the shift in the coexistence curve, δµ cap ≡ µ sat -µ cap , is well described by the classical Kelvin equation

δµ cap = 2σ cos θ ∆ρL + • • • (1)
where θ is the contact angle, σ is the surface tension of the liquid/vapour interface and ∆ρ = ρ l -ρ v is the difference in bulk number density. For thin-film uni-axial (Ising-like) magnets there is an analogous shift in the value of the external field h, at which up-spin and down-spin-like phases coexist [6,8]. In this magnetic notation, "condensation" between phases with net positive and negative magnetisation occurs when h cap ≈ -σ cos θ/m 0 L where m 0 is the spontaneous magnetisation and σ is the surface tension of the free up-spin/down-spin interface. Modern density functional theory (DFT) studies and computer simulations have shown that, at low temperatures, far from the capillary critical point, the Kelvin equation remains accurate for slit widths as small as 10 molecular diametersprovided the walls are only partially wet by the liquid, corresponding to θ > 0. For complete wetting, θ = 0, (or indeed complete drying θ = π) on the other hand corrections to the Kelvin equation arise due to the presence of wetting layers of liquid (vapour for the case of drying) at each wall and are particularly important for systems with long-ranged (dispersion-like) interactions in the wetting layers [9]. Similar (fluctuation-induced) modifications to the Kelvin equation arise in lower dimensions due to the long-ranged entropic repulsion of the interface from the wall(s) [10].

In this paper, we point out that, for systems with shortranged forces, complete wetting layers also strongly influence the decay of (density-density) correlations near capillary condensation. While there has been some discussion of the moments of the correlation function and, in particular, their relation to thermodynamic observables such as the force of solvation [1,11], we are not aware of any detailed studies of how wetting influences the asymptotic decay of the correlation functions in the capillary slit geometry. We ask the following question: what is the correlation (function) between particles situated on opposite side of the capillary? On the "condensed" side of the capillary condensation transition, µ = µ + cap , it is reasonable to assume the correlation function across the slit is similar to that of a bulk liquid. On the "vapour" side of the transition, µ = µ - cap , however it is possible that the correlation is enhanced due to the interaction between the fluctuating liquid-vapour interfaces associated with the wetting films. We shall show that the full quantitative form of this enhancement is particularly revealing and reflects the non-local nature of interfacial interactions. Indeed, this leads to universal scaling behaviour and a novel reinterpretation of the Kelvin equation.

Our article is arranged as follows: we begin by recalling some simple background theory about complete wetting transitions and the Kelvin equation. We then turn attention to the nature of correlation functions, and use effective interfacial Hamiltonian descriptions to predict how the spatial decay across the capillary depends on the slit width. We focus primarily on three-dimensions (3D) in systems with shortranged forces, and stress some subtle quantitative differences between predictions based on a local and non-local description of interfacial phenomena. Only the latter is fully consistent with calculations based on a more microscopic squaregradient density-functional theory. We finish with a summary of our main predictions and also make some remarks about the decay of correlations in 3D with long-ranged (dispersion) forces, and also in 2D concentrating on Ising-like systems.

II. BACKGROUND THEORY

Consider the interface between a planar substrate (wall) of infinite area and a bulk fluid at chemical potential µ and temperature T , below that of bulk criticality. If µ > µ sat , the bulk fluid is a high density liquid while for µ < µ sat it is a low density gas. Exactly at bulk two-phase co-existence the surface tensions of the respective wall-vapour (σ wv ) and wall-liquid (σ wl ) interfaces are related through Young's equation for the contact angle of a macroscopic drop:

σ wv = σ wl + σ cos θ. (2)
Complete wetting corresponds to θ = 0, in which case a macroscopically thick layer of liquid intervenes between the wall and bulk vapour when µ = µ - sat [START_REF] Dietrich | Phase Transitions and Critical Phenomena[END_REF][START_REF] Schick | Liquids at interfaces[END_REF][START_REF] Forgacs | Phase Transitions and Critical Phenomena[END_REF][START_REF] Bonn | [END_REF]. The complete wetting transition refers to the continuous divergence of the thickness π ∼ δµ -β co s of this layer as saturation is approached (δµ = µ sat -µ → 0). Here, the subscript on the thickness refers to the semi-infinite nature of the planar wall-fluid geometry. Associated with this transition are the vanishing of a singular contribution to the (wall-vapour) surface tension, σ sing ∼ δµ 2-α co s , and the divergence of a correlation length parallel to the unbinding liquid-vapour interface, ξ ∼ δµ -ν co (and an associated interfacial roughness ξ ⊥ ) arising from its fluctuations. Only one of the exponents is free since 1 -α co s = -β co s and α co s = 2ν co which follow from the Gibbs adsorption equation and a compressibility sum-rule respectively. The phenomenon of complete drying (θ = π) is completely analogous to this but instead refers to the divergence of a macroscopic layer of vapour at a wall-liquid interface as µ → µ + sat .

Complete wetting transitions can also be seen in semiinfinite Ising-like magnets [16][17][18]. Suppose for example that the spins are subject to a bulk field h ≤ 0, (which ensures that the bulk phase has negative magnetisation) but also that a strong positive surface field h 1 orders the spins at the surface (wall). As h → 0 -a wetting layer of up-spins grows and intrudes between the surface and the bulk down-spin state. Provided one is above the roughening temperature (for models based on lattices) the unbinding (and associated fluctuations) of the up-spin/down-spin interface are the same as those occurring for the liquid-vapour interface in fluid systems.

The modern theory of complete wetting, which is in excellent agreement with experiment [START_REF] Dietrich | Phase Transitions and Critical Phenomena[END_REF][START_REF] Schick | Liquids at interfaces[END_REF][START_REF] Forgacs | Phase Transitions and Critical Phenomena[END_REF][START_REF] Bonn | [END_REF], builds on the pioneering work of Frumkin and Derjaguin in the 1940s [9,19]. Theorists have two main approaches; Density functional models and effective interfacial Hamiltonians which each have complementary strengths (and weaknesses). Within DFT the equilibrium one-body density profile is found from minimisation of a grand-potential functional. The advantage of density functional theory is that it is microscopic, taking as its starting point a molecular description of short-ranged repulsive and long-ranged attractive intermolecular forces. Thus modern DFTs based on fundamental measure theory can describe very well packing effects which give rise to density oscillations when a high density liquid is close to a wall [20].

The disadvantage of the mean-field DFT approach is that it does not include fluctuation effects arising from the thermally induced wandering of the unbinding interface. To model these one must usually rely on interfacial Hamiltonian models which describe long wavelength fluctuations in a collectivecoordinate (x), representing the local height of the liquidvapour interface (or up-spin/down-spin interface for Ising-like magnets) at position x along the wall. The standard local interfacial model, which is partly phenomenological in character, has the form [START_REF] Forgacs | Phase Transitions and Critical Phenomena[END_REF] H

[ ] ≡ dx σ 2 (∇ ) 2 + W ( ) (3)
where W ( ) is the binding potential function describing the direct interaction of the interface with the wall (arising from intermolecular forces) which can be determined from some more microscopic theory -most commonly a DFT for fluids or a Landau-Ginzburg-Wilson (LGW) model for Ising-like magnets [21]. By definition, the binding potential is identified as the excess-free energy per unit area of a wetting film which is constrained to be at a uniform height . There is therefore an implicit assumption of locality in the standard interfacial model in that the interface interacts with the wall via its local height only. This is something we will return to later. Note that in writing (3) we have specialised to isotropic fluid interfaces in which the coefficient of the gradient term (the stiffness) is identified with the surface tension. This is also appropriate for Ising-like systems close to the bulk-critical point.

Thus the first term in (3) arises simply due to the increase in surface area induced by a fluctuation. Implicit in the Hamiltonian is a high-momentum cut-off Λ of the order of an inverse bulk correlation length. A mean-field analysis of fluctuation effects at complete wetting in three-dimensional systems suffices for almost all purposes, implying that the film thickness π and excess freeenergy σ sing can be determined from the minimisation conditions

W ( π ) = 0 σ sing = W ( π ) (4)
respectively. The length scales ξ and ξ ⊥ can be determined from the structure factor (the parallel Fourier transform of the height-height correlation)

S(q) ≡ dx 12 δ (x 1 )δ (x 2 ) e iq.x12 (5)
where δ (x) ≡ (x)π . At mean-field, S(q) = 1/ C(q) where C(q) is the Fourier transform of the "direct" correlation function C(x 12 ) ≡ β

δ 2 H δ (x1)δ (x2) where β = 1/k B T . Clearly C(q) = β (W ( π ) + σq 2 ) (6)
so that

S(q) = k B T W ( π ) + σq 2 (7)
and identifies the parallel correlation length as ξ = σ/W ( π ). Integrating S(q) over all wave-vectors determines the (squared) interfacial roughness ξ 2 ⊥ ≡ δ 2 in three dimensions as

ξ 2 ⊥ = k B T 2π Λ 0 dq q W ( π ) + σq 2 (8)
The presence of the high momentum, or short-distance, cutoff is essential to regularise this integral and leads to

ξ 2 ⊥ ≈ k B T 2πσ ln Λξ (9)
which diverges (weakly) as one approaches complete wetting [22,23]. It is this length scale that is missing in current meanfield DFT descriptions of wetting transitions.

For systems with short-ranged forces, the appropriate choice of binding potential, is

W ( ) ≡ h + ae -κ + • • • (10)
where the coefficient of the exponential repulsion a > 0 and κ ≡ 1/ξ b is the inverse bulk correlation length ξ b of the phase that wets the wall. Here h is proportional to the bulk-ordering field; specifically h = ∆ρ δµ for fluids and h = 2m 0 |h| for Ising-like magnets, where m 0 is the absolute value of the magnetization of the bulk phase. The physical origin of this binding potential will be discussed later. Using this binding potential it follows that

κ π ∼ -ln h (11)
and the associated free-energy singularity σ sing ∼ -ξ b h ln h.

Similarly for the parallel correlation length we obtain

ξ = σ κ h ( 12
)
The values of the critical exponents are therefore β co s = 0(ln), α co s = 1 and ν co = 1/2. Even though the upper critical dimension for complete wetting with short-ranged forces is d * = 3, explicit renormalisation group analysis of the interfacial Hamiltonians shows that these mean-field predictions are largely unaltered by fluctuation effects. Only the coefficient of the logarithmic divergence is (slightly) altered while the expression for ξ is unchanged [24].

Let us now consider a fluid confined between two identical planar walls, of infinite area A, which are located in the planes z = 0 and z = L respectively, and which is also in contact with a bulk reservoir at chemical potential µ < µ sat and temperature T . The Kelvin equation ( 1) follows from consideration of the total grand potential Ω for vapour-like (v) and liquid-like (l) phases that may coexist in the slit at temperatures far below the capillary critical point [1]. For the gas-like phase we write

Ω v A = -pL + σ(L; µ) ( 13
)
where p is the pressure of the bulk vapour and σ(L; µ) is the total surface tension of the finite-width capillary. To first approximation we write this as σ(L; µ) ≈ 2σ wv corresponding to the tensions of two semi-infinite wall-gas interfaces. For the liquid-like phase on the other hand

Ω l A ≈ -p * L + 2σ wl ( 14
)
where p * is the pressure of a metastable bulk liquid and σ wl is the surface tension of a semi-infinite wall-liquid interface.

Phase co-existence in the slit occurs when Ω l = Ω v implying that p -p * = 2(σ wv -σ wl )/L. This reduces to eq. ( 1) on using Young's equation σ wv = σ wl + σ cos θ and the bulk relation p -p * = (ρ l -ρ v ) δµ. The Derjaguin correction to the Kelvin equation for complete wetting arises from the singular contribution to the surface tension σ(L; µ) due to the wetting films. Neglecting the interaction between the wetting layers on either side of the strip we now approximate σ(L; µ) ≈ 2 (σ wl + σ + σ sing ) where, as above, σ sing ∼ h2-αs . On making use of the Gibbs adsorption equation this leads to an effective reduction in the slit width L eff appearing in the Kelvin equation. As first shown by Derjaguin [9] for dispersion forces L eff = L -3 π where recall π ≈ δµ -1/3 . For short-ranged forces on the other hand the correction has a more obvious geometrical interpretation, L eff = L -2 π , and is not quite so important owing to the slow logarithmic divergence of π .

III. CORRELATION FUNCTIONS AT CAPILLARY CONDENSATION IN SYSTEMS WITH SHORT-RANGED FORCES

A. The Zeroth Moment in systems with short-ranged forces. To do this, we will need to account for the interaction between the wetting films on either side of the slit. To this end, we write σ(L; µ) = 2σ wl + 2σ + W sing [START_REF] Bonn | [END_REF] and seek to determine the width dependence of the singular contribution to the excess-free energy. At mean-field level, which is accurate for almost all purposes in 3D, we can determine W sing as the minimum value of a binding potential function which allows for the interaction of the interfaces with the walls and also each other. We write

W L ( 1 , 2 ) = W ( 1 ) + W (L -2 ) + ∆W ( 1 , 2 ) (16)
where 1 and 2 denote the positions of two planar constrained liquid-vapour interfaces representing the thickness of the wetting layers (see Fig. 1). The interaction of each interface with the (closest) wall is modelled using the binding potential W ( ) defined for a semi-infinite geometry (10) while for the interface-interface interaction we use

∆W ( 1 , 2 ) = -ce -κ( 2-1) + • • • (17)
where c ∝ σ is the coefficient of an exponential attraction between the interfaces. Strictly speaking the inverse length scale κ now corresponds to the inverse bulk correlation length of the gas phase (or down-spin phase in a magnetic context), but for simplicity we shall assume an Ising-like symmetry for the bulk phases. The form of ∆W ( 1 , 2 ) may be easily derived from an underlying square-gradient DFT (or LGW Hamiltonian) by considering the excess free-energy of a constrained interfacial configuration. We will comment more on this later when we consider a non-local treatment of interfacial interactions. Since, at equilibrium, the wetting films have the same thickness at each wall we can define W ( ; L) ≡ W L ( , L -) and determine W sing as the minimum of

W ( ; L) = 2( h + ae -κ ) -ce -κ(L-2 ) + • • • (18)
Provided the bulk ordering field is not too small, the interaction between the wetting films has a negligible influence on the equilibrium thickness which is extremely well approximated by the semi-infinite expression (11). This is certainly the case close to capillary condensation where, according to the Kelvin equation the film thickness ≈ ξ b ln L L. We note in passing that the interaction term does become important for much smaller values of δµ and, indeed, eventually induces a spinodal when h ≈ e -2κL/3 corresponding to thick wetting films of order ≈ L/3. However, this metastable region is of no relevance to our present discussion of correlation functions near capillary condensation. It follows that near capillary condensation W sing ≈ 2 h π (δµ co ), so that equating [START_REF] Schick | Liquids at interfaces[END_REF] and [START_REF] Forgacs | Phase Transitions and Critical Phenomena[END_REF]

leads to hcap = 2σ (L -2 π ) + • • • (19)
which is the modified Kelvin equation for complete wetting with short-ranged forces.

We now return our attention to our main concern and introduce the connected density-density correlation function G(r 1 , r 2 ) = δρ(r 1 )δρ(r 2 ) where δρ(r) = ρ(r) -ρ(r) is the fluctuation (operator). It is convenient to exploit the translational invariance in the directions parallel to the wall and define a transverse or parallel Fourier transform

G(z 1 , z 2 ; q) ≡ dx 12 G(r 1 , r 2 ) e iq.x12 (20)
where x 12 is the parallel displacement vector between the two points at heights z 1 and z 2 from one of the walls. As a starting point for interpreting the correlation function we first consider the zeroth moment G0 (z 1 , z 2 ) ≡ G(z 1 , z 2 ; 0) corresponding to the zero wave-number limit. Indeed, exact sum-rules and also density-functional calculations shows that, for particle positions on either side of the slit, this function is related formally to the derivative of the total surface tension via [START_REF] Henderson | Fundamentals of Inhomogeneous Fluids[END_REF] G0 (0, L) ∝ -

∂ 2 σ(L) ∂L 2 . ( 21
)
Therefore, for the vapour-like phase, provided δµ ≥ δµ cap , we can expect that

Gv 0 (0, L) ∝ e -κ(L-2 π ) . (22)
As with the modified Kelvin equation, this expression depends on the reduced effective slit width L -2 π . Substituting for π , we have

Gv 0 (0, L) ∝ e -κL h2 (23)
so that, exactly at capillary condensation,

Gv 0 (0, L) ∝ L 2 e -κL . ( 24
)
This result is notable because of the polynomial prefactor ∝ L 2 (arising from the reduction in the effective slit width) and contrasts with the expected result for the liquid-like phase

Gl 0 (0, L) ∝ e -κL . (25)
Therefore correlations are enhanced due to the presence of wetting films. Note that the last expression ( 25) has two consistent interpretations. Firstly, model calculations for systems with short-ranged forces show that for the liquid-like phase the grand-potential ( 14) contains higher-order exponential corrections ∝ e -κL which lead directly to (25) using the sum-rule (21). More physically we know that away from the immediate vicinity of the walls where there are packing effects, the density profile is almost constant, in the liquid-like phase. Hence, we can anticipate that the pair-correlation function decays as it does in the bulk, provided of course that we are far below the capillary critical temperature. In other words, above the Fisher-Widom line, where the correlation function and associated inhomogeneous density profiles decay monotonically [START_REF] Fisher | [END_REF]27], we anticipate a simple Ornstein-Zernike-like decay where r ≡ |r 1 -r 2 | is the particle separation. The zeroth parallel moment of this, G0 ∝ e -κz , decays as a pure exponential, so that, for particles on either side of the slit, we can expect [START_REF] Henderson | Fundamentals of Inhomogeneous Fluids[END_REF]. We emphasize again that the predictions ( 24) and ( 25) are only valid far below the capillary critical point. For wide capillary slits, this is equivalent to the condition that L ξ b .

G l (r 1 , r 2 ) ∝ e -κr r (26)

B. Capillary-wave theory

Let us now extend the analysis to calculate the paircorrelation function across the slit in the presence of wetting films assuming that δµ ≥ δµ cap . The simplest approach to this is based on a local interfacial Hamiltonian

H[ 1 , 2 ] = dx σ 2 (∇ 1 ) 2 + σ 2 (∇ 2 ) 2 + W L ( 1 , 2 ) (27)
in which fluctuations are resisted by the surface tension and the interface interaction is modelled using the binding potential (16) with the expressions (10) and (17). To determine G for particle positions close to the walls on opposite sides of the slit we need to determine the height-height correlations between the two different interfaces. To this end we define a (2 × 2) matrix S(q) with elements

Smn (q) ≡ dx 12 δ m (x 1 )δ n (x 2 ) e iq.x12 (28)
where δ m (x) = m (x)m . Recall that only two elements are independent since by symmetry S11 = S22 and S12 = S21 . At mean-field level, the matrix S(q) is obtained as the inverse of the matrix

C(q) = β      ∂ 2 W L ∂ 2 1 + σq 2 ∂ 2 W L ∂ 1 ∂ 2 ∂ 2 W L ∂ 2 ∂ 1 ∂ 2 W L ∂ 2 2 + σq 2      (29)
the elements of which are simply the Fourier transforms of the second functional derivatives δ 2 H δ mδ n evaluated at equilibrium. Close to capillary condensation, C11 (q) is essentially unchanged from the expression (6) for complete wetting in the semi-infinite geometry, while for the off-diagonal element is determined as the second derivative of ∆W ( 1 , 2 ). Thus,

C(q) ≈ β   W ( π ) + σq 2 -κ 2 ce -κ(L-2 π ) -κ 2 ce -κ(L-2 π ) W ( π ) + σq 2   (30)
This leads to the two structure factors

β S11 (q) ≈ 1 κ h + σq 2 (31)
and

β S12 (q) ≈ cκ 2 e -κ(L-2 π ) (κ h + σq 2 ) 2 . ( 32
)
The first result is the same expression (or rather very nearly the same) as the structure factor for complete wetting at a single wall, eq. ( 7). The weak correlation between the two interfaces is quantified by the off-diagonal term S12 which is exponentially small in the slit width. Note also that the denominator is the square of that appearing in S11 (q). To determine the density-density correlation function we note, that in local effective Hamiltonian theories, a local fluctuation in the interfacial height 1 (x) produces an analogous local density fluctuation at the wall δρ ≈ e -κ 1 δ 1 (x). This can be thought of as arising from the tail of the exponentially decaying density profile from the liquid-vapour interface. It follows that the Fourier transform G(0, L; q) of the densitydensity correlation function for particle positions on either side of the slit is given by

Gv (0, L; q) ∝ e -κL (κ h + σq 2 ) 2 . ( 33
)
Exactly at capillary condensation, this reads

Gv (0, L; q) ∝ L 2 e -κL (1 + L 2κ q 2 ) 2 . ( 34
)
Setting q = 0, this recovers our earlier sum-rule result for the zeroth moment (24) demonstrating the consistency of the calculation. Note that the large length scale present in the wavevector dependence is simply the parallel correlation length ξ ∝ h-1/2 describing the interfacial fluctuations in the wetting film evaluated at the chemical potential corresponding to capillary condensation. Substituting the Kelvin equation expression δµ ≈ 2σ/(∆ρL) into (12) yields

ξ = L 2κ (35)
which is the same as that appearing in (34).

Finally we need to perform the inverse Fourier transform to determine the desired decay in real space:

G v (0, L) = C L 2 e -κL 2π Λ 0 dq q (1 + L 2κ q 2 ) 2 (36)
where C is an unimportant constant. It is interesting that this result is essentially independent of the cut-off Λ because the high momentum behaviour is integrable. One can see this explicitly if one introduces a simple change of variable to rewrite the expression as

G v (0, L) = C Le -κL 2πξ b I CW (37)
where

I CW = Λ 2 ξ b L/2 0 dz 1 (1 + z) 2 (38)
is a dimensionless integral. Obviously this is non-singular and approaches unity in the limit L → ∞. Thus within the local effective Hamiltonian theory we are led to the prediction 6 which shows explicitly the enhancement of correlations across the slit compared to correlations in the condensed, liquid phase. The non-universal constant C is the same as that appearing in the zeroth moment, and can be eliminated by defining a dimensionless amplitude ratio Lξ b G v (0, L)/ Gv 0 (0, L) or, more accurately, the limit of this as L → ∞. Thus, according to the present local Hamiltonian, this dimensionless amplitude ratio is given by

G v (0, L) = C Le -κL 2πξ b (39)
ξ b L G v (0, L) Gv 0 (0, L) ≈ 1 2π (40)
the value of which we shall return to shortly.

C. Non-local theory

The analysis described above, based on a local interfacial Hamiltonian, is reliable regarding most of its predictions. However, when we compare it to the results of a more microscopic DFT we find that it is not quite correct (see Appendix). In this section, we show that this can be traced directly to non-local effects associated with the damping of interfacial interactions and has a similar interpretation to recent studies of non-locality at short-ranged wetting transitions [28][29][30][31].

One may improve on the basic phenomenological interfacial Hamiltonian of wetting by starting from a more microscopic LGW theory and integrating out degrees of freedom systematically. This has been done recently for (short-ranged) wetting in a semi-infinite geometry and has provided new insights into a number of long-standing problems in the theory of critical wetting transitions [32]. We recall some details of this approach, concentrating only on what we need for complete wetting transitions, before describing how the analysis generalises to the parallel-plate geometry. The systematic derivation of an interfacial model was first considered by Fisher and Jin [33] who used a perturbative approach (as an expansion in the gradient ∇ ). However as described in [29], it is possible to do this non-perturbatively using a Green's function method closely related to multiple reflection expansion techniques [34] . This has a number of advantages and allows one to derive an interfacial model for wetting at nonplanar walls, the position of which we denote by a function ψ(x).

The non-local (NL) Hamiltonian can be written

H NL [ , ψ] = σA lv + hV + W NL [ , ψ] (41)
and is valid for length scales larger than the bulk correlation length. Curvature corrections related to rigidity-like terms can be calculated but are not considered here. The interfacial area is denoted A lv while V = dx( (x)-ψ(x)) is the volume of the wetting film. The binding potential functional W NL [ , ψ] describes the interaction of the interface and wall and has an expansion

W N L [ , ψ] = aΩ 1 1 [ , ψ] + • • • (42)
with a geometry independent coefficient a which is the same as the binding potential (10). Higher order contributions to W NL exist but are not important for complete wetting transitions [35]. Each of these contributions has a diagrammatic representation in which a zigzagging straight line (each zig and zag representing a kernel ) connects points on the wall (lower wavy-line) to points on the interface (upper wavy-line).

The dominant diagram has the algebraic expression

Ω 1 1 [ , ψ] = = ds ψ ds K(r ψ , r ). ( 43
)
Here r ψ and r are points on the wall and interface respectively, which are integrated over using the appropriate infinitesimal areas ds ψ = dx 1 + (∇ψ) 2 and ds = dx 1 + (∇ ) 2 . These surface integrations are implied by the black circles in the diagrams. The kernel, denoted by each straight line, is

K(r 1 , r 2 ) = κ 2π|r 1 -r 2 | e -κ |r1-r2| (44)
and is simply a normalised version of the bulk correlation function.

Specialising to the case of wetting at a planar wall, the NL Hamiltonian reduces to

H NL [ ] = dx σ 2 (∇ ) 2 + h + W NL [ ] (45)
where we have taken the usual gradient expansion for the bending term. If the wall is flat the integral over it can be done exactly in which case the wetting diagram reduces to

= ds e -κ (46)
where ds = dx 1 + (∇ ) 2 . One can see immediately that this contribution is essentially identical to that appearing in the standard interfacial Hamiltonian arising from the (leading) exponential decay of the binding potential (10). That is, for a planar wall, the leading-order contribution to the binding potential functional is local in character. This is not the case for higher order terms, which are important at critical wetting transitions, but for complete wetting the form of the local Hamiltonian ( 3) is supported by the NL theory. Consequently the predictions for the parallel correlation ξ and structure factor S(q) are the same as described in section III A. Despite this the description of density-density correlation functions is better within the NL theory because the relation between fluctuations in the interfacial height and density fluctuations is controlled explicitly by the kernel K. Indeed this is an additional source of non-local effects with the result that the NL description is in precise agreement with an exact sumrule [31]. The singular (interfacial) contribution to the pair correlation function is given by where r 1 and r 2 are the correlated points, while r 3 = (x 3 , π ) and r 4 = (x 4 , π ) lie on the interface. This has the short-hand diagrammatic representation

G(r 1 , r 2 ) ∝ ∂ 2 ∂z 1 ∂z 2 dx 3 dx 4 K(r 1 , r 3 )S(x 34 )K(r 4 , r 2 ) (47)
G(r 1 , r 2 ) ∝ ∂ 2 ∂z 1 ∂z 2 (48)
where the wiggly line represents the height-height correlation function S(x 12 ). In Fourier space, this reads G(z 1 , z 2 ; q) ∝ K (z 1π ; q) K (z 2π ; q) S(q) (49) where

K(z; q) = κ κ q e -κq|z| (50)
is the parallel Fourier transform of the Green's function (44), the prime denotes the derivative w.r.t. z and we have abbreviated κ q = κ 2 + q 2 . Note that the length scale controlling the exponential decay of correlations as we move away from the interface is κ q and not κ (as is often assumed in local theories). Thus for particle positions very close to the wall, the NL theory predicts that G contains a singular term

Gsing (0, 0; q) ∝ h 1 + ξ 2 q 2 e -π q 2 /κ . (51)
Note the presence of a damping factor e -π q 2 /κ which shows that, in addition to ξ , there is a diverging length scale ξ NL ∝ √ π which strongly damps correlations. It is the allowance for this length in the NL theory which gives consistency with the exact sum-rule requirement that the second moment G2 (0, 0) sing ∝ σ + σ sing [30,31].

The generalisation of the NL theory to wetting films in a parallel-slit geometry is straightforward. The interaction of each interface with the (closest) wall is the same as in the semi-infinite geometry and is essentially local in character. The difference with the local Hamiltonian approach however is that the interaction between the interfaces is now strongly non-local. Thus, we have to write

H NL [ 1 , 2 ] = dx σ 2 (∇ 1 ) 2 + σ 2 (∇ 2 ) 2 + W L [ 1 , 2 ]
(52) where

W L [ 1 , 2 ] ≡ dx {W ( 1 ) + W (L -2 )} + ∆W [ 1 , 2 ]
(53) analogous to (16). The interaction between the interfaces is now modelled by the non-local expression

∆W [ 1 , 2 ] ≡ -cΩ 1 1 [ 1 , 2 ] (54)
where, in an obvious notation,

Ω 1 1 [ 1 , 2 ] ≡ ds 1 ds 2 K(r 1 , r 2 ) (55)
Again, we can visualise this interaction diagrammatically as

Ω 1 1 [ 1 , 2 ] ≡ (56)
where the two horizontal lines represent the capillary walls. If the two interfaces are planar then this integral reduces to

= A e -κ( 2-1) (57)
where A is the area. This exactly reproduces the interaction term (17) in the binding potential function (16) implying that, for planar interfaces, the NL Hamiltonian is (per unit area) identical to the binding potential function W L ( 1 , 2 ) appearing in the local theory. Hence, at mean-field level, the excess free-energy and therefore the zeroth moment of the pair correlation function G0 (0, L) are the same as calculated earlier.

The determination of the structure factor S12 (q) follows along the same lines described in III A. The difference between the local and NL theories arises because now the offdiagonal terms for the direct correlation function matrix is wave-number dependent

C12 (q) ≈ -κ 2 ce -κq(L-2 π ) (58)
and recall κ q = κ 2 + q 2 . This feeds directly into the structure factor

β S12 (q) ≈ cκ 2 e -κq(L-2 π ) (κ h + σq 2 ) 2 (59)
which only matches with the corresponding expression of the local theory when q = 0. To calculate the pair-correlation function, we use (49) which is equivalent to

Gv (0, L; q) ∝ e -2κq π S12 (q) (60)
where we have assumed that the wetting layer thickness is the same as at a single wall. Hence, for h ≥ hcap , we have

Gv (0, L; q) ∝ e - √ κ 2 +q 2 L (κ h + σq 2 ) 2 . ( 61
)
This is similar, but not quite identical, to the result (34) we arrived at using a local capillary-wave theory because of the wave-vector dependence in the numerator. Since the NL Hamiltonian is only valid for length scales longer than the bulk correlation length we can expand κ q ≈ κ + q 2 /2κ + • • • and re-write this as

Gv (0, L; q) ∝ e -κL e -q 2 Lξ b /2 (κ h + σq 2 ) 2 (62)
which identifies an additional damping term in the numerator ∝ e -Lq 2 ξ b /2 compared with the corresponding result of the local theory. Thus both approaches are in agreement with regard to the value of the zeroth moment but differ for larger wave-vectors. It is apparent that the structure factor across the slit is characterised by two large length scales: the usual parallel correlation length ξ which, for h ≥ hcap , is well approximated by the semi-infinite result [START_REF] Dietrich | Phase Transitions and Critical Phenomena[END_REF], and a second length

ξ damp = Lξ b 2 (63)
which emerges from the damping term e -Lq 2 ξ b /2 ≡ e -q 2 ξ 2 damp . It is clear that this damping effect is very similar to that which occurs for the correlation function G(0, 0; q) at complete wetting discussed earlier. The extra length scale controlling the damping in this case was ξ NL = √ π ξ b . The damping across the capillary is similar but now occurs on the scale of the slit width rather than the wetting film thickness.

A remarkable feature of the length scales ξ damp and ξ is that at capillary condensation, hcap ≈ 2σ/L, they are equal. In other words the Kelvin equation ( 1) is equivalent to the statement that

ξ damp = ξ ( 64
)
where ξ is given by [START_REF] Dietrich | Phase Transitions and Critical Phenomena[END_REF]. Expressed in terms of the slit width, we find for the Fourier transform Gv (0, L; q) ≈ CL 2 e -κL e -Lq 2 /2κ (1

+ L 2κ q 2 ) 2 (65)
which is valid for wave-numbers q Λ ≈ κ, and L ξ b . The scaling of the correlation lengths (64) therefore implies that

G v (0, L) = C Le -κL 2πκ I NL ( 66
)
where

I NL = Λ 2 ξ b L/2 0 dz e -z (1 + z) 2 (67)
is different to the corresponding expression I CW of the local theory. Again this integral is very weakly dependent on L (which can be set to ∞) leading to

I NL = 1 -e E 1 (1) ≈ 0.404 (68)
where E 1 (x) is the exponential integral

E 1 (x) ≡ ∞ x dt e -t t (69)
Hence the dimensionless amplitude ratio

ξ b L G v (0, L) Gv 0 (0, L) ≈ 1 -e E 1 (1) 2π (70)
and is reduced compared to the prediction of the local interfacial Hamiltonian. In the Appendix, we show that (65) is the result obtained from directly solving the Ornstein-Zernike equation for a square gradient DFT. That is, only the nonlocal description of interfacial interactions reproduces the full quantitative features of an underlying microscopic theory.

Finally, we mention that within the non-local theory there is a simple diagrammatic expression for the correlation function across the slit, analogous to that for a single wall, eq. (48). For particle positions at distances z 1 ≈ 0 and z 2 ≈ L (and arbitrary parallel displacement), the expression is given by

G(r 1 , r 2 ) ∝ ∂ 2 ∂z 1 ∂z 2 (71)
where the lowest and highest horizontal lines represent the capillary walls, and the two middle horizontal lines represent the two liquid-vapour interfaces. The diagram has the same interpretation as earlier, in which the wavy lines denote height-height correlations S 11 (= S 22 ), and the thick line denotes the kernel K. In this way, we can see rather clearly that the total correlation across the slit arises from both bulk-like decays, represented by the thick lines, and capillary-wave excitations at both interfaces (wavy lines).

IV. DISCUSSION

In this paper, we have discussed how the presence of wetting films influences the density (or magnetisation) correlation function G(0, L) between particle positions on either side of a capillary slit (or Ising thin-film). Concentrating on systems in three dimensions with short-ranged forces, our central predictions are as follows:

1) The amplitude of the exponential decay of correlations is enhanced due to the presence of the wetting films. Near capillary condensation, the correlation function in the vapour phase decays as G v (0, L) ≈ Le -κL compared to G l (0, L) ≈ e -κL /L in the liquid phase.

2) The enhancement of correlations arises from the interaction between the fluctuating interfaces on either side of the slit and, unlike the interfacial roughness, does not depend on the high-momentum cut-off in the effective Hamiltonian.

3) The attractive interaction between the wetting films must be modelled using a non-local interaction potential which we can represent diagrammatically as

∆W [ 1 , 2 ] = -c . ( 72
)
For planar interfacial configurations this recovers the usual binding potential function present in local descriptions. This non-local description correctly accounts for the damping of interfacial interactions controlled by a length scale

ξ damp = Lξ b 2 (73)
4) At capillary condensation hcap ≈ 2σ/L the damping length scale and the parallel correlation length are equal: This scaling of the characteristic lengths leads to the following amplitude ratio for the correlation function to its zeroth moment

ξ damp = ξ . (74)
ξ b L G v (0, L) Gv 0 (0, L) ≈ 1 -e E 1 (1) 2π . (75)
This may be regarded as a signature of non-locality. Our analysis has only considered Gaussian fluctuations about the mean-field approximation. However, we believe that our quantitative predictions remain valid even if one did a more thorough renormalisation group analysis based on the interfacial models. As remarked earlier, even though the upper critical dimension for complete wetting with short-ranged forces is d * = 3, it is known that fluctuations do not alter the mean-field expression for the parallel correlation length [START_REF] Dietrich | Phase Transitions and Critical Phenomena[END_REF]. Since both the form of the exponential decay of correlations and the expression for ξ are not altered by interfacial fluctuation effects in three dimensions, we believe all our predictions are exact (for sufficiently large slit widths L).

We finish our article by making some brief remarks about systems with long-ranged forces and also correlation functions near (pseudo) capillary condensation in two dimensions:

Long-Ranged Forces -Let us imagine that instead of being short-ranged, the attractive component of the fluid-fluid (and wall-fluid) intermolecular potential decays as φ(r) ≈ -r -γ where, of course, γ = 6 for dispersion forces. In the bulk, this long-ranged decay, also describes the large-distance behaviour of the correlation function. It is possible to generalise the above interfacial Hamiltonian analysis and determine how the correlation function behaves in the slit geometry in the presence of wetting films. Consider a local treatment first using the interfacial Hamiltonian (27) but now with a binding potential appropriate for long-ranged forces. For the wallinterface interaction we use the standard binding potential [START_REF] Dietrich | Phase Transitions and Critical Phenomena[END_REF][START_REF] Schick | Liquids at interfaces[END_REF][START_REF] Forgacs | Phase Transitions and Critical Phenomena[END_REF][START_REF] Bonn | [END_REF] W

( ) = h + a -p + • • • (76)
where p = γ -4 (in three dimensions). To describe the interinterfacial interaction we now use

∆W ( 2 , 1 ) ∝ ( 2 -1 ) -p + • • • (77)
which is simply the excess free-energy per unit area of two planar parallel interfaces separated by a distance 2 -1 .

Repeating our analysis, we find that this local treatment predicts that, at capillary condensation, the correlation decays as G(0, L) ∼ L -q where the exponent q = (γ-2) 2 γ-3 , i.e., an algebraic decay which is different to the intermolecular potential! For example, q = 16/3 < 6 for dispersion forces. This result is, in fact, wrong owing to the incorrect local description of the interface-interface interaction. Recall that, for simple liquids, it is often assumed that the contribution to the Helmholtz functional F [ρ], arising from attractive intermolecular forces is

F att = 1 2 dr 1 dr 2 ρ(r 1 )ρ(r 2 )φ(r 12 ). (78)
Applying a simple sharp-kink approximation to the density profile in this expression leads directly to an interaction term

∆W [ 1 , 2 ] ∝ (79)
where, following the diagrammatic notation of [35], the upper and lower black dots (which now lie off a surface) indicate an integral over the volume of the upper and lower wetting layer. The dashed line now represents the intermolecular potential φ(r) rather than the Yukawa function K(r). If both 1 and 2 are planar ∆W [ 1 , 2 ] = A∆W ( 1 , 2 ) (where A is the area), so that the local theory is recovered in the long-wavelength limit. It is straightforward to show that the non-local interaction (79) corrects the erroneous prediction of the local theory and leads to

G(0, L) ≈ L -γ , (80)
i.e., the pair correlation decays algebraically with the same exponent as the intermolecular potential. Wetting films therefore do not strongly influence the decay of correlations as they do for systems with short-ranged forces.

Two-Dimensional Systems -Finally, we mention that correlation functions at two-dimensional capillary condensation are also of interest. Consider, for example, a two dimensional Ising strip of infinite length, finite width L and surface spins fixed to +1 (infinite surface field). Since the system is only infinite in one direction, it is pseudo-one dimensional and the capillary condensation "transition" is rounded. However, for wide slits, the cross-over from up-spin to down-spin phases is sharp and occurs over an exponentially narrow region centred around hcap [36]. It is still possible, therefore, to compare how the correlation function between spins on either side of the slit decays compared to that in the bulk.

Recall that exactly in zero bulk field, and at sub-critical temperatures, the spin-spin correlation function decays as G(r) ≈ e -κr /r 2 where the inverse length scale is related to the surface tension κ = 2βσ [37]. These features of the bulk correlation function are very beautifully explained by Abraham's bubble model of correlations [38]. In this picture, correlations in a phase of net down-spins (say) arise from bubbles of up-spins which enclose the two points being correlated. The sum over all Boltzmann weighted bubble configurations determines the correlation function G(r) ∝ Σe -βE . In Fig. 2 we schematically illustrate the lowest energy bubble excitation in the down-spin and up-spin-like phases close to condensation. For the down-spin phase, the bubble that connects spins on either side of the slit has the form of a capillary bridge made from two semi-circular menisci. For the upspin phase, on the other hand, a circular bubble of down-spins stretches across the slit connecting the two correlated spins. The energy cost of both these configurations is determined by simple macroscopic considerations involving surface tensions (conjugate to the length of the domain wall) and the coupling to the bulk magnetic field (conjugate to the area of the excitation). Precisely at capillary condensation the energy cost of In 2D systems, spins on either side of the strip are correlated when excitations corresponding to bubbles of the appropriate metastable phase extend between the two spins. Left) The elementary excitation on the down (-) spin phase is a capillary bridge of up (+) spins. Right) The elementary excitation on the up (+) spin phase is a circular bubble of down spins.

both the bridge and circular bubble excitations are the same and a simple calculation leads to the prediction G(0, L) ≈ e -πκL/4 (81)

for both the up-spin and down-spin phase. Here, we have ignored a possible multiplicative algebraic prefactor which depends on L and is unimportant compared to the exponential term. It is remarkable therefore that near to capillary condensation the correlation function decays with an effective (inverse) length scale κ eff = π 4 κ. Since κ eff < κ, it follows that correlations are therefore enhanced in the slit geometry. We believe this prediction can be checked numerically in the Ising model using density-matrix renormalisation group methods [START_REF]Density-Matrix Renormalization[END_REF]40] or perhaps simulations.

  FIG.2. In 2D systems, spins on either side of the strip are correlated when excitations corresponding to bubbles of the appropriate metastable phase extend between the two spins. Left) The elementary excitation on the down (-) spin phase is a capillary bridge of up (+) spins. Right) The elementary excitation on the up (+) spin phase is a circular bubble of down spins.

  FIG.1. Schematic illustration of coexisting vapour-like and liquidlike phases in a capillary slit geometry of width L. On the vapour side of the condensation transition, thick wetting films of liquid form at the substrates if the contact angle is zero.

					2
				L
	substrate	vapour	substrate	Condensation substrate _	substrate
	liquid	liquid	liquid
		z	o F	z
				r
				P
				e e r
					R e v i e w
					O n
					l y

URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics

ACKNOWLEDGMENTS

It is our pleasure to dedicate this paper to Bob Evans on the occasion of his (unexpected) 65th birthday. Throughout his academic life, he has been nothing other than an inspiration. A.O.P. is particularly grateful for the wonderful memories of his time in Bristol.

NRB acknowledges support from the Portuguese Foundation for Science and Technology (SFRH/BPD/63183/2009). CR acknowledges support from grants MOSAICO (Ministerio de Educación y Ciencia) and MODELICO (Comunidad de Madrid).

Appendix A: Correlation Functions within a Square-Gradient DFT

In this appendix, we determine the pair-correlation function across a capillary slit using a simple square gradient DFT. We adopt a magnetic notation and consider the LGW free-energy functional

with two bounding walls in the z = 0 and z = L planes. The local order parameter (magnetization) is denoted m(r), and the potential φ(m(r)) has an appropriate double-well structure modelling bulk two-phase coexistence at subcritical temperatures. We use the simplest possible choice, corresponding to a double-parabola approximation

where the bulk "liquid" and "gas" phases have order parameters m 0 and -m 0 at zero field, respectively. We fix the magnetization at the walls to be m 1 , and note that the condition m 1 > m 0 corresponds to complete wetting, θ = 0. The equilibrium magnetization profile is determined by the Euler-Lagrange equation

whose solution, for the double-parabola approximation, is readily expressed in terms of elementary exponential functions. The correlation function

is determined from solution of the Ornstein-Zernike equation [41]. In Fourier space, the parallel transform G satisfies

where we have set k B T = 1. Again within double-parabola approximation, this equation is readily solved since φ (m) = κ 2 , provided m = 0 [31]. For the "vapour-like" phase, the full solution for G is

) where

Here, denotes the equilibrium position of the interface (where m = 0) from the wall at z = 0. We now consider the case of a very wide capillary slit, L l 1/κ, and small magnetic fields, |h| κ 2 m 0 . We can then approximate X 1 ≈ X 2 ≈ X 3 ≈ 1 and assume that the interface position is located in the same position as in the semi-infinite geometry. Also, setting z 1 ≈ 1/κ and z 2 ≈ L -1/κ, the correlation function across the slit reduces to

Expanding κ q to order q 4 in the denominator and q 2 in the numerator, this reduces to Gv (0, L; q) ≈ 2 κ 3 σ 2 e -κL e -q 2 Lξb/2 (2|h|m 0 κ + σ q 2 ) 2 . (A7) where σ ≡ κm 2 0 is the surface tension.


Jose ; Tavares

J M Tavares

P I C Teixeira

Patching up dipoles: Can dipolar particles be viewed as patchy colloids?

Keywords: Liquid-vapour transition, Self-assembly, Magnetic fluids PACS numbers: 64.70.F-, 64.75.Yz, 47.65.Cb

We investigate whether the liquid-vapour phase transition of strongly dipolar fluids can be understood using a model of patchy colloids. These consist of hard spherical particles with three short-ranged attractive sites (patches) on their surfaces. Two of the patches are of type A and one is of type B. Patches A on a particle may bond either to a patch A or to a patch B on another particle. Formation of an AA (AB) bond lowers the energy by ǫ AA (ǫ AB ). In the limit ǫ AB /ǫ AA ≪ 1, this patchy model exhibits condensation driven by AB-bonds (Y-junctions). Y-junctions are also present in low-density, strongly dipolar fluids, and have been conjectured to play a key role in determining their critical behaviour. We map the dipolar Yukawa hard-sphere (DYHS) fluid onto this 2A + 1B patchy model by requiring that the latter reproduce the correct DYHS critical point as a function of the isotropic interaction strength ǫ Y . This is achieved for sensible values of ǫ AB and the bond volumes. Results for the internal energy and the particle coordination number are in qualitative agreement with simulations of DYHSs. Finally, by taking the limit ǫ Y → 0, we arrive at a new estimate for the critical point of the dipolar hard-sphere fluid, which agrees with extrapolations from simulation.

Even after some forty years of efforts we still do not know whether dipolar forces alone suffice to condense a vapour into a liquid. We shall not, of course, attempt to solve this difficult problem here. Rather, we will show how some of our results for a different model system, where liquid-vapour coexistence is also elusive, may be relevant to the criticality of dipolar fluids. But first we give a brief review of the current state of our knowledge.

In the classic Van der Waals picture of phase transitions, simple fluids interact via isotropic intermolecular potentials that comprise a short-ranged repulsion and a longerranged attraction. Condensation is then driven by the free energy balance between the high-entropy vapour phase (where attraction dominates) and the low-energy liquid phase.

(where repulsion prevails).

Most molecular species, however, exhibit permanent dipole moments. The simplest model of a dipolar fluid is the dipolar hard sphere (DHS) fluid, which consists of hard spheres of diamater σ, having at their centre a dipole momnent of strength µ. The dipole-dipole interaction is strongly anisotropic and long-ranged, hence notoriously difficult to treat theoretically. For this reaeon, a number of authors [1][2][3][4][5] have sought to derive isotropic effective interaction potentials that nevertheless capture the essentials of the true dipolar fluid, for the purpose of easier computation of its phase behaviour. In all cases, the leading term of this effective potential is attractive and has the same distance dependence as, e.g., the long-range part of the Lennard-Jones (LJ) potential. de Gennes and Pincus [6] conjectured that the phase diagram of DHSs should be similar to that of a Van der Waals fluid, with vapour, liquid and solid phases. Surprisingly, however, numerical simulations of DHS [7] and also of Stockmayer fluids [8] have revealed that the anisotropy of the dipolar potential promotes the formation of selfassembled aggregates (chains, rings and more complex clusters) if the dipolar interaction strength is of the order of the thermal energy. This is in sharp contrast with the isotropic compact clusters observed in simple fluids. Moreover, unlike in simple fluids the pair correlation function of DHS is strongly peaked at contact and the internal energy is nearly independent of the density [9]. It is at present unclear whether strong association pre-empts condensation, or any other kind of fluid-fluid phase separation, as failure to observe it may be an artifact of the simulation techniques [10][11][12][13]. New methods are needed to allow a re- assessment of the dipolar and related condensation problems in the light of recent theoretical results [14].

Thus dipolar particles may or may not condense, but they certainly associate. The simplest theories of association as applied to this system [15][16][17][18][19] assume that the only effect of the strongly anisotropic interparticle interaction is to drive cluster formation. Consequently, the fluid is described as an ideal mixture of self-assembling clusters, failing to exhibit any phase transitions unless direct or indirect interactions between the clusters are added. Nevertheless, these theories reproduce well the slow variation of the internal energy with the density and cluster size (or mass) distribution. Alternative treatments of the competition between phase separation and association in the Stockmayer and DHS fluids include variants of the Flory-Huggins (FH) model of equilibrium polymerisation [20,21] and a thermodynamic perturbation theory for associating fluids [22,23]. The former rely on casting the free energy of the Stockmayer fluid in FH form; the latter includes a contribution from association to the free energy, but does not distinguish between chaining and branching, the two types of self-assembly that appear to influence the phase behaviour of dipolar fluids.

A different approach to the phase behaviour of the dipolar fluid was proposed in [24].

Following the results for the structure suggested by simulations (long chains that may form branches, see, e.g. [18]), it was assumed that the dipolar fluid at low densities and temperatures can be described as a perturbation of a ground state that consists of infinitely long chains. This perturbation consists of the appearance of two types of thermally excited defects: end defects, which correspond to the formation of shorter chains, and Y-junction defects, which correspond to the merging of three chains. Chain ends and junctions carry energy penalties ǫ e and ǫ j , respectively. The entropic gain of defect forming is related to the lattice model from which the free energy was derived [25]. It was shown that a critical point would exist only if the energy cost of forming a junction is small enough, ǫ j < ǫe 3 . When such a critical point exists, phase coexistence is obtained between a low-density gas of short chains (i.e., rich in end defects) and a high-density networked liquid (i.e., rich in Y-junctions). More strikingly, it was found that the phase diagram displays a pinched, or reentrant, behaviour: as temperature is decreased, the densities of the coexisting liquid and vapour phases become more similar, rather than more different, as is the norm in most fluids. This model thus introduces a new perspective on association in dipolar fluids: the anisotropy of the dipolar potential promotes competition between chaining and branching, which determines the peculiar phase behaviour of DHSs. Safran's model can be quantitatively related to the dipolar fluid by expressing the two energy scales ǫ j , ǫ e , in terms of the dipole moment. The shortcomings of this description of the DHS fluid are: (i) the difficulty of relating ǫ j and the dipole moment; (ii) its origin in a lattice model, which hinders quantitative comparion with continuous models like DHSs; and (iii) the asymptotic character of its thermodynamics, which is only valid at extremely low densities and temperatures.

In recent work [26], we have established an analogy between Tlusty and Safran's (TS) model [24] and patchy particles with three bonding sites (patches), two of type A and one of type B. In this model, henceforth denoted 2A + 1B model, particles can bind to each other through their patches, thereby forming AA bonds that lower the energy by ǫ AA , and AB bonds that lower the energy by ǫ AB . If ǫ AA > 2ǫ AB , the ground state consists of long chains of particles connected by AA bonds. Two end defects are created when a long chain breaks into two, with an energy cost of ǫ AA . A Y-junction is formed when, after breaking a long chain into two, one of the ends (a non-bonded A patch) bonds to an interior particle of another chain (to a non-bonded B patch), creating a AB bond. The energy cost of this Y-junction is thus -ǫ AB + ǫ AA /2. Moreover, it was also shown that Wertheim's first order perturbation theory [27][28][29][30], as applied to the patchy model predicts the same phase behaviour as TS: a critical point can only exist if ǫ AB > ǫ AA /3 (which is equivalent to ǫ j > ǫ e /3 in the TS model); coexistence, when present, is between a low-density phase with few AA and AB bonds, and a high-density phase rich in AB and AA bonds; and pinching is also observed [31].

Despite the similiraties between the TS model and the patchy model, there are advantages in using the latter. These come from the fact that the patchy model, combined with Wertheim's perturbation theory, yields a thermodynamics that is able to describe a fluid.

In fact, TS relies on a lattice model and, in its thermodynamic description, the entropy of the ends and junctions is solely related to the properties of the underlying lattice [25]. By contrast, the entropy of the bonded patches is related to their size and to the properties of the reference system [26,32]. Therefore, when comparing both models with simulation results, better agreement is expected for the patchy model.

Our purpose is to provide a simple analogy between the 2A + 1B model and DHSs. This Yukawa hard-sphere (DYHS) fluid [13,23], for decreasing values of the isotropic Yukawa interaction strength. This paper is organized as follows: in section II we detail the analogy between the patchy model and the DYHS (II A), and use Wertheim's perturbation theory to derive the thermodynamics (II B). Results for the critical parameters, phase diagrams, internal energy and particle coordination numbers are presented and discussed in section III. Finally in section IV we conclude and highlight further directions for research.

II. THEORY A. Analogy between DYHSs and patchy particles

The DYHS fluid studied in [13,23] consists of hard spheres (HSs) of diameter σ with a central point dipole of strength µ and interacting through a potential that comprises two contributions: the anisotropic dipole-dipole interaction and a Yukawa isotropic attraction of energy scale ǫ Y and decay parameter z (we have set z = 1.8 for consistency with [13,23]).

The 2A + 1B model of patchy particles will then consist of HSs of diameter σ, which attract isotropically through the same Yukawa potential as the DYHS fluid, but where the central dipole has been replaced by two patches of type A and one patch of type B, placed on the HS surface [26]. These patches may form AA and AB bonds, which lower the energy by ǫ AA and ǫ AB , respectively. To proceed with this analogy, ǫ AA and ǫ AB must be related to the dipolar energy scale µ 2 /σ 3 . A sketch of this analogy is given in figure 1.

When ǫ AA > ǫ AB /2 (and ǫ Y = 0) the ground state of the 2A+1B patchy particles consists of long chains connected by AA bonds. Likewise, the ground state of DHSs consists of long chains (or rings) of particles aligned head to tail. Because the absolute minimum of the DHS pair potential (the energy of two dipolar spheres at contact and aligned head to tail) is 2µ 2 /σ 3 , we shall take

ǫ AA = 2µ 2 /σ 3 .
(

) 1
This first relation between the DHS and 2A + 1B models is therefore obtained by matching their ground sates. together with the breakup of short chains, is expected to appear as a low-energy, thermallyexcited perturbation of the ground state of infinite chains [24]. A Y-shaped arrangement of N DHSs that results from bringing together the ends of three chains of aligned dipoles (see figure 1b) has energy -2(N -3)µ 2 /σ 3 -u 3 µ 2 /σ 3 , where -u 3 µ 2 /σ 3 is the dipolar interaction energy between the three dipoles that make up the junction. On the other hand, the Yshaped arrangement of N patchy particles that results from bonding an A patch at a chain end with a B patch on an interior particle of another chain, has energy -(N -2)ǫ AA -ǫ AB .

Matching these two results and using equation ( 1), we obtain

ǫ AB = (u 3 -2)µ 2 /σ 3 . (2)
u 3 can be evaluated by considering three DHSs in contact (see figure 2). The configurations depicted show that a simple approximation for u 3 should be avoided. In fact, both configurations seem reasonable for Y-arrangements of dipoles: (a) results from slightly bending one chain and attaching another chain at the bend; (b) is a highly symmetric (ring-like) configuration that corresponds to the lowest possible energy. Choosing (a) or (b) has completely different consequences for the phase behaviour [26,33]: (a) corresponds to ǫ AB /ǫ AA ≈ 0.105 (i.e., ǫ AB /ǫ AA < 1/3) and thus to the absence of a critical point, whereas (b) corresponds to ǫ AB /ǫ AA = 0.875 (i.e., ǫ AB /ǫ AA > 1/2) and thus to the existence of a critical point that is not of the TS type, since the ground state is not a gas of long chains [31]. In order to avoid this somewhat arbitrary (and possibly biased) choice, we will let u 3 (or, equivalently, ǫ AB ) be a free parameter of the patchy model, and determine it by fitting the DYHS thermodynamics as determined by simulations [13] to that of our 2A + 1B model. Notice that this apparent shortcoming has one obvious advantage: we cannot know in advance whether the DHS fluid will exhibit a critical point. This conclusion will depend on the value obtained for ǫ AB , which is determined not only by the analogy we propose, but also by the simulation results.

The bonding of A and B patches is characterised not only by two energy scales, but also by the bond volumes, v AA and v AB , [26,32], which are related to the patch sizes. We will let these two volumes be free parameters of the patchy model, and they will also be determined by fitting the thermodynamics of the DYHS fluid to that of 2A + 1B patchy particles. The free energy per particle f of the 2A + 1B patchy model with an additional Yukawa attracton will be approximated by:

βf = βf HS + βf b + βf Y , (3)
where β ≡ 1/(k B T ), T is the temperature, k B is the Boltzmann constant, f HS is the HS contribution (which we approximate by the the Carnahan-Starling CS expression [34]), f b is the bonding contribution due to the patches, and f Y is the Yukawa contribution.

For HSs with two A and one B patches, Wertheim's theory [27,28,36] gives the following expression for the bonding contribution to the free energy [26,35]:

βf b = 2 ln X A + ln X B -X A - X B 2 + 3 2 , ( 4
)
where X α is the probability of having a patch of type α not bonded. The variables X α are related to the density and temperature through the laws of mass action that are derived by treating bond formation as a chemical reaction [27,28,36]. For the case where only AA and AB bonds can form, these are [26,35],

X A + 2η∆ AA X 2 A + η∆ AB X A X B = 1, (5)
X B + 2η∆ AB X A X B = 1, (6)
where η ≡ (N/V )v s is the packing fraction, and

∆ αβ = 1 v s v αβ g ref (r) [exp(βǫ αβ ) -1] dr. ( 7
)
This integral is calculated over v αβ , the volume of bond αβ, and g ref is the pair correlation function (PCF) of the HS reference system. These equations yield the probability p α = 1 -X α that a site of type α is bonded, as a function of density, temperature, and the interaction strengths. As we shall see, this enables us to calculate structural properties, in addition to thermodynamic properties.

The integral in equation ( 7) is evaluated by replacing the PCF by its value at contact in the CS approximation [34]: This is expected to be reasonable provided that the typical size of the patches (and, consequently, that of the bonds) is ≪ σ. Within this approximation, equation (7) becomes

g ref (r) = g ref (σ + ) = 1 -η 2 (1 -η) 3 . (8
∆ αβ = v αβ v s [exp(βǫ αβ ) -1] 1 -η 2 (1 -η) 3 . (9)
The contribution to the free energy coming from the Yukawa attraction, f Y , will be given by a simple mean-field (MF) approximation, as used, e.g., in [36]:

βf Y = - 1 2 βη φ Y (r 12 )dr 12 = -2πβ z -1 + z -2 ǫ Y η, (10)
with, as in [12], z = 1.8. We expect equation (10) to become increasingly inaccurate as attractive Yukawa forces become dominant over dipolar forces and the density of the coexisting phases (as well as that of the critical point) increases. A slightly more sophisticated alternative consists in setting

βf Y = -βB Y 2 η ≈ 2π 3 × 3.097 (βǫ Y ) 1.0724 η ( 11
)
where B Y 2 is the second-virial coefficient of the HS Yukawa fluid, interpolated from the numerical results in [37] after subtracting the HS contribution (which is already contained in βf hs ). Below we shall assess the relative merits of equations ( 10) and (11); it turns out that there is very little difference between them in the ranges of parameters (ρ, T, ǫ Y ) of interest.

We are now in possession of a free energy, equation (3), which is a function of, besides (η, T ), ǫ AB , v AA , v AB (through equation ( 4)) and ǫ Y (through equation (10) or equation (11)). We shall first derive the formal expressions for the critical density and temperature, and then search for the values of ǫ AB , v AA and v AB that best reproduce the critical densities and temperatures found in the DYHS simulations of [12] for several values of ǫ Y . Once these parameters are known, all thermodynamic quantities of interest can be calculated, as well as some structural properties.

III. RESULTS

We introduce reduced units defined as usual: reduced temperature 1 of [13], for a range of Yukawa interaction strengths. This was achieved using the least-squares routine LMDIF from MINPACK library, coupled with NETLIB routine HYBRD for solving the critical point equations, which yielded the parameters collected in table I. The quality of fits, as measured by their root mean square deviation, does not depend strongly on how many points we fit, or how we treat the attractive Yukawa term. Note that, for all fits, one has ǫ AB /ǫ AA < 1/2, i.e., results are consistent with the hypothesis that the ground state of the system is long chains. These results thererfore indicate that competition between chaining and branching plays a prominent role in the phase behaviour. Note also that, for all fits but one, ǫ AB /ǫ AA > 1/3, which strongly suggests that there exists a critical point in the DHS limit. Consequently, we predict that the critical point of the DHS fluid (corresponding to ǫ * Y = 0) will occur at T * c ≈ 0.14 -0.16, ρ * c ≈ 0.03 -0.05, which is in line with earlier estimates reviewed in [13]. In what follows we have used the set of parameters shown in bold, chosen on the basis of (i) simplicity of the MF treatment of attractions; (ii) sensible values for v AA and v AB . The other sets of parameters (except the very first one, for which ǫ AA /ǫ AB < 1/3 ) nevertheless yield very similar results. Finally, one apparent shortcoming of these fits is that the values of v αβ are larger than those usually adopted in patchy particle models (see, e.g., [33]), which must not violate Wertheim's theory assumption that no more than two particles connect at any given bond. Further work is needed to clarify if the bond volumes in table I satisfy this requirement.

T * = k B T σ 3 /µ 2 ; reduced Yukawa energy parameter ǫ * Y = ǫ Y σ 3 /µ 2 , reduced dipolar interaction parameter (µ * ) 2 = µ/(ǫ Y σ 3 ) = 1/ǫ
In figure 3 we plot the reciprocal of the critical temperature in Yukawa units (figure 3a), the critical temperature (figure 3b) and the critical density (figure 3c). Agreement, which is in general quite good, deteriorates slightly as we go to larger ǫ * Y , i.e., as we approach the non-polar Yukawa limit. This is a consequence of the fact that we fitted only to the six data points with smallest ǫ * Y (see above), coupled with our rather simplistic treatment of Yukawa attractions at the MF level.

We found the phase diagram by equating the pressures and chemical potentials of the coexisting liquid and vapour phases, together with the mass-action law equations ( 5) and (6) for the two phases. The resulting set of six coupled non-linear equations was solved numerically using NETLIB routine HYBRD. Results are shown in figure 4, for different values of the Yukawa interaction parameter ǫ * Y . As ǫ * Y → 0 and the DHS limit is approached, the phase diagram develops a 'pinch', or re-entrance, at low temperatures: the densities of the coexisting liquid and vapour phases become more similar, rather than more different, as is the norm in most fluids. This is a consequence of the fact that 1/3 ≤ ǫ AB ≤ 1/2, and was previously observed in a Monte Carlo simulation of a variant of our patchy fluid with two A and nine B patches [31]. In the present model this pinching only appears at the smallerst ǫ Y , for which the chaining/branching competition dominates over the effect of isotropic interactions.

The change in internal energy as a function of density at T * = 0.15 and ǫ * Y = 0.025, can bee seen in figure 5a (to be compared with figure 4 of [13]). Whereas the total energy decreases with density, in agreement with the classical picture of the liquid (higher-density) phase as low-energy and low-entropy, the patchy contribution saturates, thus reproducing a caractheristic footprint of DHSs. The latter contribution actually equals the internal energy of the patchy fluid at the same density and temperature, since our theory treats bonding and Yukawa attractive interactions separately.

In figure 5b we plot several structural quantities at the critical point, in order to describe the evolution of the structure of the fluid when ǫ * Y is varied. For extremely small ǫ * Y (i.e., close to the DHS limit), ≈ 90% of A patches and ≈ 15% of B patches are bonded. The number of AA bonds per particle is very large, hence there must be many long chains, forming Y-junctions whenever a B patch is bonded. This scenario seems to hold up to ǫ * Y ≈ 0.1, with slightly shorter chains (a decreasing fraction of AA bonds) and a larger number of Y-junctions (an increasing fraction of AB bonds). However, for ǫ * Y > 0.1 there is a dramatic change in the fluid structure at the critical point: the number n of chains becomes negligible (there is a sharp decrease in the fraction of AA bonds) and the number of AB bonds, although decreasing, is larger than the number of AA bonds. Since the critical temperature and the critical density are both increasing, AB bonds become favoured relative to AA the Yukawa attraction and not by the competition between chaining and branching. This result shows that, as already discussed in [13], the DHS limit obtained from any model that approaches it, must be taken with care, since dramatic changes in the strucure and thermodynamics are expected.

Finally, in figure 6 we show the variation with density at T * = 0.15 of the fractions of particles bonded to zero (n 0 ), one (n 1 ), two (n 2 ) or three (n 3 ) other particles, corresponding to, respectively, isolated particles, chain ends, particles inside chains, or Y-junctions. These quantities are calculated using the probabilities p A = 1-X A and p B = 1-X B , and assuming that bonds are independent. As the density increases n 0 +n 1 decreases rather sharply whereas n 3 gently increases, implying that the lower-density (vapour) phase is richer in chain ends and unattached particles, whereas the higher-density (liquid) phase is richer in junctions.

n 2 , on the other hand, first rises sharply, indicating that chains are increasing in length, then reaches a maximum at ρ * ∼ 0.04, and finally drops very slowly, with increasing ρ * .

This slow decrease of n 2 (together with the increase in n 3 ) is of entropic origin: AA bonds are turning into AB bonds because v AB ≫ v AA . Except for the initial increase in n 2 , the qualitative trends are the same as in the nearest-neighbour histograms of figure 6 of [13];

note, however, that the results in the latter figure do not extend below ρ * = 0.025, and that n 2 appears to be peaking around that value (see middle panel of figure 6 in [13]).

IV. CONCLUSIONS

We have been able to map the DYHS fluid onto a model of patchy colloids with three interaction sites, two of type A and one of type B, where B sites do not interact. This was achieved by requiring that this model reproduce the critical densities and temperatures found from simulation, for a range of values of the Yukawa interaction parameter ǫ * Y . Wertheim's association theory as applied to this model is able to reproduce the correct trends of the particle coordination number: on increasing the density at constant temperature, the fraction of particles with zero (isolated particles), one (chain ends) or two (inner particles in chains) nearest neighbours decreases, whereas that of particles with three or more nearest neighbours (junctions or particles inside compact clusters) increases. Moreover, we predict that the dipolar contribution to the internal DYHS energy should saturate at high densities, again in agreement with simulations. A more interesting consequence of using the 2A+1B patchy model is that, if our mapping is valid, then the phase diagram of DYHSs should be rentrant, or 'pinched' in the ǫ * Y → 0 limit: for T * < ∼ (2/3)T * c , the density difference between the coexisting liquid and vapour phases decreases as the temperature is further lowered. This is in stark contrast to the behaviour of most fluids, but has actually been confirmned by simulation of 2A + 9B patchy colloids [31].

Observation of the same effect for actual DYHS would be a stringent test of our approach.

In future work we plan to apply our mapping to other models that, while remaining simple to simulate and study theoretically, converge to the DHS on varying some parameter, e.g., mixtures of polar and non-polar HSs in the limit of infinite dilution of the latter. In the DHS fluid this has energy -2(N -3)µ 2 /σ 3 -u 3 µ 2 /σ 3 (-u 3 µ 2 /σ 3 is the dipolar interaction energy of the three particles that have come together). In the 2A + 1B model, this forms when a free A patch A at a chain end bonds to a B patch on an interior particle of another chain; its As in figure 4 of [13], the dipolar energy saturates at higher densities. Notice that in our theory, the dipolar (i.e., patchy) contribution to the internal energy of the DYHS fluid is identical to the internal energy of the DHS fluid, since the Yukawa interaction does not affect bonding. (b) Degrees of non-association at the cricial point, X Ac and X Bc .: both approach 1 (the fully non-associated limit) as ǫ * Y grows and the fluid becomes less DHS-like and more Yukawa-like. We also plot the number of AA and AB bonds per particle at the critical point: n AA /N ≡ 1+X Bc

F

energy is -(N -2)ǫ AA -ǫ AB . 15

3

 3

4 Page

 4 analogy will be quantitatively tested by comparing our theoretical results for the thermodynamics and the structure of the patchy model, with those of simulations of the dipolar

5 Page

 5 The relation between ǫ AB and the dipolar energy is established by considering branched Y-shaped arrangements (see figure1b) containing N particles. It is this structure that,

) 7 Page

 7

10 Page 10 of 20 URL

 1020 bonds, because the associated entropy is larger on account of v AB ≫ v AA . Nevertheless, the total number of bonds decreases strongly and the fluid becomes less and less associated as ǫ * Y increases, as expected. Therefore, above ǫ * Y ≈ 0.1 criticaliy is already determined by

11

 11

14 Page 14 of 20 URLFIG. 1 :

 14201 FIG. 1: Sketch of the mapping of the DHS fluid onto the 2A + 1B patchy model. The small black (white) circles represent patches of type A (B). (a) A linear chain of N DHSs, in contact and aligned head-to-tail, has energy -2(N -1)µ 2 /σ 3 , whereas a row of N patchy particles has N -1 AA bonds and thus energy -(N -1)ǫ AA . (b) A Y -shaped arrangement, made up of N particles.

Page 15 of 20 URLFIG. 2 : 15 4 µ 2 /σ 3 . 16 Page 16 of 20 URLFIG. 3 :

 20215216203 FIG. 2: Two examples of orientational configurations of three DHSs at contact. The dipolar energy of configuration (a) is (-2 √ 3 + 5 4 )µ 2 /σ 3 ≈ -2.21µ 2 /σ 3 , and that of the ring configuration (b) is -15 4 µ 2 /σ 3 .

17

 17

Page 17 of 20 URLFIG. 4 : 18 Page 18 of 20 URLFIG. 5 :

 20418205 FIG. 4: Temperature-density phase diagram for different Yukawa interaction parameters ǫ * Y as shown, and ǫ AB , v vAA and v AB as in figure 3. As ǫ * Y → 0 the diagrams 'pinch' (become re-entrant)at low temperatures, as seen in[31].

2 - 19 Page 19 of 20 URL 3 FIG. 6 : 20 Page 20 of 20 URL

 21920362020 FIG.6: Fractions of particles either not bonded or bonded to one other particle (n 0 + n 1 ), bonded to two other particles (n 2 ), or to three other particles (n 3 ), vs density at T * = 0.15, for ǫ AB , v AA and v AB as in figure3. In our theory n i do not depend on ǫ * Y . Compare with figure6of[13].

  HS diameter, as our unit of length, and σ 3 /µ 2 as our unit of energy. Note that this choice of energy scale is neutral from the point of view of criticality, as AA bonds can only lead to chaining, and not phase coexistence, at any finite temperature. The DHS limit corresponds to ǫ * Y = 0. We started by estimating the model parameters ǫ AB , v AA and v AB , by requiring that our 2A + B patchy model reproduce the critical temperatures and densities listed in Table

	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y
	8

* Y ; reduced number density ρ * = ρσ 3 . This amounts to taking σ, the

TABLE I :

 I Fit parameters for the DYHS. n is the number of points from Table1in[13] included, starting at the lowest ǫ * Y . For each set we also give the estimates for the critical temperature and density of the DHS fluid (corresponding to ǫ * Y = 0).

	n Yukawa ǫ AB /ǫ AA	v AA /σ 3	v AB /σ 3 RMS deviation T * c (DHS) ρ * c (DHS)
	8 eq. (10) 0.284805 0.00222266 0.459357	0.00691166	-	-
	7 eq. (10) 0.347745 0.00234394 0.231612	0.00564450	0.142	0.033
	6 eq. (10) 0.37328 0.00220629 0.173476 0.00565144	0.150	0.040
	8. eq. (11) 0.390595 0.00259764 0.145093	0.00697690	0.151	0.040
	7. eq. (11) 0.46689 0.00282981 0.0665826 0.00695441	0.162	0.049
	6. eq. (11) 0.46331 0.00247637 0.0668395 0.00657444	0.161	0.049
	o F o
	r r
	P P
	e e r e e r
			R e v i e w R e v i e w
				O n O n
					l y l y

Table 1

 1

in

[13]

, the lines are fits from our theory, using the values boldfaced in table I for ǫ AB , v AA and v AB . See the text for details.

URL: http://mc.manuscriptcentral.com/tandf/tmph

Acknowledgement

Financial support from the Foundation of the University of Lisbon and the Portuguese Foundation for Science and Technology (FCT) under Contracts nos. POCI/FIS/55592/2004, POCTI/ISFL/2/618 and PTDC/FIS/098254/2008, as well as through Pluriannual contracts with CFTC, is gratefully acknowledged.