Anna Oleksy

Jean-Pierre Hansen

O n l y Wetting and drying scenarios of ionic solutions

Keywords: wetting, drying, ionic solutes, density functional theory

Wetting and drying scenarios of ionic solutions.

Solutes can profoundly modify the interfacial properties, in particular the surface tension γ, of solvents. While surfactants radically reduce the surface tension of water/oil interfaces, ions, on the contrary, lead to a significant increase of the surface tension of the water/air interface [1], which may be traced back to strong Coulombic interactions [2][3][4]. Since γ, which is the surface excess contribution to the free energy of an inhomogeneous fluid, controls the wetting of a solid (or liquid) substrate by a saturated vapour, or the conjugate drying behaviour of the coexisting liquid phase, it is to be expected that any solute will affect the wetting and drying transitions of the pure solvent. However, while the experimental investigation [5], and theoretical understanding [6] of wetting transitions are well advanced for pure fluids and mixtures, much less attention has been devoted to the study of wetting and drying phenomena of ionic solutions. The reason for this may lie in the complication linked to the additional length scale, the Debye screening length λ D , which competes with the usual length scales, namely the film thickness ζ, the width δ of the liquid/vapour interface, and the range of the substrate/fluid attraction. λ D controls the strength of the electrostatic repulsion between the electric double-layers which form at the substrate/fluid and liquid/vapour interfaces, and hence affects the thickness of the liquid wetting film [START_REF] Langmuir | [END_REF][8][9][10]. However, the influence of such effects on the wetting and drying scenarios of ionic solutions has only recently been investigated theoretically, within a mean-field framework, using the tools of density functional theory of inhomogeneous fluids [START_REF] Evans | Fundamentals of Inhomogeneous Fluids[END_REF][START_REF] Hansen | Theory of Simple Liquids[END_REF].

In the present paper we review this recent work [START_REF] Denesyuk | [END_REF][14][15][16][17], and present some new results on the variation of the substrate/fluid and liquid/vapour surface tensions with ion concentration. After a brief reminder of mean field theories of wetting and drying transitions in simple fluids (Section 2), we introduce the models and approximations used to describe inhomogeneous ionic solutions (Section 3). In Sections 4-6 we present the DFT results obtained with increasingly "realistic" representations of ionic solutions, namely the "primitive" model (Section 4), the "semi-primitive" model (Section 5) and the "civilised" model (Section 6). Our new data concerning the ion concentration dependence of the surface tension are presented in Section 7, while conclusions are drawn in Section 8.

Mean field theories of surface phase transitions.

The present discussion will be restricted to three-phase equilibrium between a "spectator phase" a, taken here to be a planar solid substrate, and two fluid phases (liquid b and vapour c). The competition, as a function of temperature and proximity of liquid/vapour coexistence, of the interfacial free energies, or surface tensions γ ab , γ ac and γ bc gives rise to surface phase transitions, namely wetting and drying. If a saturated vapour is put into contact with a solid substrate, which exerts an attractive force on the vapour molecules, a microscopically thin "liquid" film forms at contact. This partial wetting situation is characterised by the thermodynamic inequality:

γ ac < γ ab + γ bc (1)
and is observed well below the liquid/vapour critical temperature T c . However, as the temperature is increased along liquid/vapour coexistence, a simple argument due to Cahn [18] shows that at a wetting temperature T w , the inequality ( 1) becomes an equality, and a macroscopically thick liquid film intrudes between the substrate and the vapour. This wetting transition is usually first order, i. e. the liquid film thickness ζ jumps discontinuously from a finite (microscopic) value below T w (partial wetting) to infinity (complete wetting). In some rare cases a second order scenario has been observed [5], where the film thickness diverges continuously as T w is approached from below (second order or continuous wetting transition). In the case of a first order wetting transition at T w , there is an additional twist, predicted theoretically by Cahn [18], and by Ebner and Saam [19] in 1977, namely a prewetting transition as vapour/liquid coexistence is approached along an isotherm T > T w from the undersaturated vapour side. The film thickness undergoes a finite discontinuity (jump) as a prewetting line is crossed; after that ζ diverges continuously as the coexistence line is approached. The prewetting line starts on the coexistence curve at T = T w , where the discontinuity is infinite, and moves into the undersaturated vapour region where the discontinuity remains finite. Its amplitude decreases with increasing T and vanishes at a prewetting critical point T pwc , above which the film thickness varies continuously.

On the liquid side of the phase diagram, conjugate drying may be observed, where a gas layer intrudes between the bulk liquid and a solid substrate. The thickness of this "vapour film" increases continuously as coexistence is approached along an isotherm, and diverges at coexistence [20,21].

The film thickness, or equivalently the adsorption Γ, is the obvious order parameter characterising the surface phase transitions. If the structureless, impenetrable substrate is taken to be the xy plane (z=0), the local fluid density, or density profile ρ(z) is the fundamental quantity within DFT of classical inhomogeneous fluids [START_REF] Evans | Fundamentals of Inhomogeneous Fluids[END_REF][START_REF] Hansen | Theory of Simple Liquids[END_REF]. The adsorption is defined as:

Γ = ∞ 0 ρ(z) -ρ (0) dz (2)
where ρ (0) is the bulk density of the fluid, far from the substrate (i. e. as z → ∞). The film thickness ζ may be defined as ζ = Γ = Γ/ρ (0) .

Within DFT the equilibrium density profile is obtained by minimizing an approximate expression for the grand potential functional (per unit area of the substrate) Ω [ρ(z)] with respect to ρ(z), for a fixed chemical potential µ, namely:

Ω [ρ(z)] = F [ρ(z)] + ∞ 0 [V(z) -µ] ρ(z)dz (3)
where F [ρ(z)] is an approximate free energy functional (per unit area) and V(z) is the potential of the force exerted by the substrate on the fluid molecules. The Euler-Lagrange equation leads to the following generic expression for the density profile [START_REF] Evans | Fundamentals of Inhomogeneous Fluids[END_REF][START_REF] Hansen | Theory of Simple Liquids[END_REF]:

ρ(z) = ρ (0) exp - 1 k B T µ[ρ(z)] -µ(ρ (0) ) + V(z) (4)
where is the "local" chemical potential, and µ(ρ (0) ) = µ is the overall (bulk) chemical potential. For a given free energy functional F, eq. ( 4) must be solved iteratively. Substitution of the solution into the surface excess grand potential yields the surface tension γ:

ρ (0) = ρ(z → ∞) is the bulk density, µ [ρ(z)] = δF [ρ(z)] δρ(z) (5)
γ = Ω ex [ρ(z)] = Ω [ρ(z)] -Ω(ρ (0) ) = F [ρ(z)] -F(ρ (0) ) + ∞ 0 V(z)ρ(z)dz -µΓ (6)
Note that the adsorption ( 2) and the surface tension are related by the Gibbs adsorption equation [START_REF] Evans | Fundamentals of Inhomogeneous Fluids[END_REF][START_REF] Hansen | Theory of Simple Liquids[END_REF]:

Γ = - ∂γ ∂µ T (7)
which provides a useful test of the thermodynamic consistency of results based on approximate functionals [22].

At this stage, one may schematically distinguish between two classes of DFTs of wetting and drying. Following the pioneering work of Cahn [18], phenomenological theories are based on generalisations of van der Waals classic "square gradient" theory of the liquid/vapour interface [START_REF] Hansen | Theory of Simple Liquids[END_REF]23,[START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF]. Microscopic theories, starting from a molecular description of the fluid, were pioneered by Ebner and Saam [19], and by Sullivan [START_REF] Sullivan | [END_REF]26], and culminated in the decisive work of Evans and collaborators [21,[27][28][29]. Following van der Waals, Cahn assumes the density profile ρ(z) to vary slowly near contact. The substrate-fluid interaction is assumed to act only at contact:

V(z) = V 0 (ρ(z)) δ(z) (8a)
with:

V 0 (ρ) k B T = γ 0 -γ 1 ρ + 1 2 γ 2 ρ 2 (8b)
The square gradient expression for the surface excess grand potential hence reads [6,18]:

Ω ex [ρ(z)] k B T = ∞ 0 W (ρ(z)) + B 2 dρ(z) dz 2 dz + γ 0 -γ 1 ρ s + 1 2 γ 2 ρ 2 s ( 9
)
where ρ s = ρ(z = 0) is the contact value of the local density. W(ρ) is the reduced bulk grand potential density (local density approximation), while the coefficient B/2 of the non-local square gradient correction is assumed to be independent of density. In the vicinity of the critical point W(ρ) may be approximated by the double well form:

W(ρ) = c(ρ -ρ l ) 2 (ρ -ρ v ) 2 (10)
where ρ l and ρ v are the temperature-dependent densities of the coexisting liquid and vapour phases, and c is another phenomenological coefficient. The equilibrium profile is easily calculated to be: where δ = B/2c(ρ l -ρ v ) -1 is the van der Waals width of the liquid/vapour interface [START_REF] Hansen | Theory of Simple Liquids[END_REF]23,[START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF], and ζ = -δ ln C = -δ ln [(ρ l -ρ s )/(ρ s -ρ v )] is the film thickness. Substituion of ( 11) into ( 9) yields the reduced surface excess grand potential in the form [START_REF] Denesyuk | [END_REF]14]:

ρ(z) = ρ v + (ρ l -ρ v ) 1 + exp [(z -ζ)/δ] (11)
Ω ex k B T ≡ ω ex = ω ex 0 + 6γ 1 2(1 + C) 2 - 1 3(1 + C) 3 - 1 6 + (p 1 -p 2 ) C 1 + C + 1 2 p 2 C 1 + C 2 (12)
where

ω ex 0 = γ + γ 0 -γ 1 ρ l + γ 2 2 ρ 2 l , γ = √ 2Bc(ρ l -ρ v ) 3 /6
is the van der Waals surface tension of the liquid/vapour interface, and p 1 and p 2 are the two control parameters

p 1 = γ 1 √ 2Bc (ρ l -ρ v ) 2 ; p 2 = γ 2 √ 2Bc (ρ l -ρ v ) -1 (13)
ω ex is finally minimised with respect to C (or equivalently to the contact density ρ s ), which determines the film thickness ζ corresponding to the lowest minimum of ω ex for any given temperature (i. e. ρ l and ρ v ). Depending on the values of the two control parameters p 1 and p 2 , a first (discontinuous) or second order wetting transition is predicted along the liquid/vapour coexistence curve. The Cahn theory also predicts the prewetting scenario as the coexistence curve is approached from the undersaturated vapour side, provided the grand potential density ( 10) is generalised to incorporate an additional contribution linear in (ρ -ρ v ) proportional to the degree of undersaturation [14,18]. Cahn's theory is very attractive because of its fully analytic nature, but it requires a large number of thermodynamic or phenomenological parameters as input (ρ l (T), ρ v (T), B, c, γ 1 , γ 2 ) which must be taken from a separate microscopic theoretical description, or from experiment.

Microscopic density functional theories are based on a molecular model of the fluid, which specifies the Hamiltonian of the system, i. e. the interactions between the molecules [19,21,[START_REF] Sullivan | [END_REF][26][27][28][29]. Following van der Waals, the pair potential is split into a short-range repulsion and long-range (LR) contributions (dispersion and Coulombic interactions). The former is conveniently taken to be of the simple hard sphere (HS) form. Accordingly, the free energy functional F[ρ] may be split into ideal, HS and LR contributions:

F [ρ(z)] = F id [ρ(z)] + F HS [ρ(z)] + F LR [ρ(z)] (14)
While F id [ρ(z)] is known exactly [START_REF] Evans | Fundamentals of Inhomogeneous Fluids[END_REF][START_REF] Hansen | Theory of Simple Liquids[END_REF]

F id [ρ(z)] = k B T ∞ 0 ρ(z) log Λ 3 ρ(z) -1 dz (15)
approximate expressions must be used for the contributions F HS and F LR due to intermolecular forces. Weighted density approximations have been used for F HS by the Bristol group [21,[27][28][29], but the work on ionic solutions reported in the following Sections is based on the more accurate "fundamental measure" functional introduced by Rosenfeld [30] in its most refined form [31,32]. The long-range contribution to the free energy is almost invariably taken to be of the non-local mean-field form: where w is the long-range part of the pair potential, suitably integrated over the transverse (x,y) coordinates.

F LR [ρ(z)] = 1 2 ∞ 0 dz ∞ 0 dz ρ(z)w(|z -z |)ρ(z ) (16)
Minimisation of the corresponding grand potential ( 3) with respect to ρ(z) can no longer be carried out analytically but requires careful iterative numerical solutions [21,[27][28][29]. The various wetting scenarios predicted by Cahn's phenomenological theory are confirmed by the "first principles" DFT approach. A key finding of the Bristol group is that the order of the wetting transition is very sensitive to the range of the substrate/fluid attraction, relative to that of the attractive potential w(z) between fluid molecules. Short-range substrate/fluid attractions favour a continuous wetting transition, while first order wetting and prewetting transitions are found to occur when the range of the attraction exerted by the substrate on the fluid molecules exceeds a threshold [27].

In the following Sections the phenomenological and microscopic DFT approaches are generalised to the case of simple models of ionic solutions [START_REF] Denesyuk | [END_REF][14][15][16][17]. It is important to stress at the outset that for any theory of wetting or drying to be consistent, the interfacial and bulk thermodynamic properties must be calculated throughout from the SAME free energy functional, which reduces to a function of the constant densities in the bulk.

3 Simple models of ionic solutions.

Ionic solutions are (at least) three-component systems involving the solvent, cations and anions. To investigate their interfacial properties next to a neutral or charged planar substrate, one must introduce three density profiles ρ 0 (z) (solvent), ρ + (z) and ρ -(z). A fundamental combination of the latter two profiles is the local charge density:

ρ c (z) = α=+,- q α ρ α (z) (17)
where q α is the charge carried by ionic species α. The profiles determine the partial adsorptions of the three species:

Γ α = ∞ 0 dz ρ α (z) -ρ (0) α ( 18
)
which are in turn related to the interfacial tension γ by the generalisation of the Gibbs adsorption equation [START_REF] Langmuir | [END_REF]:

Γ α = - ∂γ ∂µ α T,{µ β } (19)
Through their long-range Coulombic interaction and the requirement of charge neutrality, both globally, and locally through the screening mechanism, the ions introduce qualitatively different interfacial behaviour (compared to the pure solvent), essentially through the formation of electric double-layers at the substrate/fluid and liquid/vapour interfaces. We have studied how the ions modify the wetting and drying behaviour of the pure solvent, using three different models of the ionic solution:

a) The "primitive" model, where ions and solvent molecules are completely decoupled, the solvent playing only the role of a dielectric continuum of fixed permittivity ε relative to the ions. The phenomenological theory combines Cahn's theory for the solvent as sketched in Section 2, with Poisson-Boltzmann theory for the inhomogeneous "fluid" of ions. This Cahn-Poisson-Boltzmann approach [START_REF] Denesyuk | [END_REF]14] will be presented in Section 4.

b) The "semi-primitive" model (SPM) takes the molecular granularity of the solvent explicitly into account, by representing the solvent molecules by neutral hard spheres with attractive interactions between all solvent-solvent, solvent-ion and ion-ion pairs; the ions are charged hard spheres. Since the model ignores the polar nature of the solvent molecules, it must be supplemented by the introduction of a local dielectric permittivity ε(z) depending on the local solvent density to allow for solute dissociation [15,16]. This model will be considered in Section 5.

c) The "civilized" model [33,34] introduces an explicitly polar solvent, by representing its molecules by dipolar hard spheres. Apart from this crucial refinement, the model is similar to the SPM, but does not require the phenomenological introduction of a local permittivity ε(z), since the dielectric properties of the solvent are now controlled by its molecular dipolar moment. The wetting and drying behaviour of the "civilized" model [17] will be the subject of Section 6.

The primitive model: Cahn-Poisson-Boltzmann theory.

Within the primitive model, the surface excess grand potential is assumed to be the sum of two independent contributions, both functions of the film thickness ζ; the solvent contribution Ω s , derived from Cahn's square gradient theory (cf. eq. ( 9)), and an electrostatic contribution Ω el due to the electric double-layers formed by the ions at the substrate/liquid and liquid/vapour interfaces; the latter is calculated within the non-linear Poisson-Boltzmann (PB) theory of ions in a continuous solvent of dielectric permittivity ε.

Ω ex (ζ) = Ω s (ζ) + Ω el (ζ) (20)
The substrate carries a uniform surface charge chosen to be negative for clarity, σ = -q/a, where a is the area per unit charge. The counterions carry a charge +q, while the coions have a charge -q (symmetric electrolyte). The generalisation to ions of different valences is straightforward. The cation and anion density profiles within the wetting film satisfy Poisson's equation:

d 2 Φ(z) dz 2 = -4πl B [ρ + (z) -ρ -(z)] (21)
where Φ(z) is the dimensionless electrostatic potential qΨ(z)/k B T and l B = q 2 /(εk B T) is the Bjerrum length.

Assuming the ion concentration in the vapour phase above the liquid film to be negligible (i. e. ions do not leak outside the film), the potential must satisfy the boundary conditions:

dΦ(z) dz z=0 = - 4πl B a ; dΦ(z) dz z=ζ = 0 (22)
which are equivalent to the overall charge neutrality condition: The PB free energy functional is the sum of ideal and mean field electrostatic contributions:

ζ 0 [ρ + (z) -ρ -(z)] dz = 1 a (23)
f e = F e Ak B T = α=± ζ 0 ρ α (z) log Λ 3 α ρ α (z) -1 dz + 1 2 ζ 0 ρ c (z)Φ(z)dz (24)
where ρ c (z) = ρ + (z) -ρ -(z) -δ(z)/a is the total charge density (in units of q). Imposing chemical equilibrium throughout the wetting film:

µ α (z) = δf e [{ρ β (z)}] δρ α (z) = µ (0) (25)
where µ (0) is the bulk chemical potential of anions and cations, the Boltzmann form of the density profiles is recovered:

ρ ± (z) = ρ (0) exp {∓Φ(z)} (26)
Substitution of ( 26) back into eq. ( 21) results in the PB equation:

d 2 Φ(z) dz 2 = 8πl B ρ (0) sinh (Φ(z)) (27)
which must be solved, for a given film thickness ζ, subject to the boundary conditions ( 22). This can be achieved analytically, in parametric form [START_REF] Denesyuk | [END_REF]14]. The solution may be substituted in the free energy functional [START_REF] Rowlinson | Molecular Theory of Capillarity[END_REF] to yield an analytic expression for the electrolyte contribution to the surface excess grand potential:

ω e = Ω e Ak B T = f e -µ (0) ζ 0 [ρ + (z) + ρ -(z)] dz (28)
as a function of the film thickness ζ. Gathering the results for the pure solvent (eq. ( 12)) and for the electrolyte, one obtains the total surface excess grand potential as a funcction of ζ. The following scenarios emerge: a) In the simple case where only counterions (released by the substrate into the film upon ionisation) are present (i. e. ρ -(z) ≡ 0), the electrostatic contribution ω e is always positive, and decays like 1/ζ due to the lack of screening, while the solvent contribution ω s decays exponentially with ζ.

Consequently the variation of ω ex = ω s + ω e with film thickness is dominated by ω e as z → ∞, and the wetting transition will always be first order, irrespective of the discontinuous (first order) or continuous (second order) nature of the wetting transition of the pure solvent [START_REF] Denesyuk | [END_REF]14].

b) In the more general case of added salt (i. e. the presence of counterions and coions), the wetting behaviour is controlled by three dimensionless parameters, in addition to the parameters p 1 and p 2 , defined in eq. ( 13), which control the wetting scenarios of the pure solvent; these three parameters are determined by the two additional physical variables, namely the reservior chemical potential µ (0) , or equivalently the reservoir ion concentrations ρ

(0) + = ρ (0)
-= ρ (0) (symmetric electrolyte) and the surface charge density σ = q/a (the three parameters are hence not fully independent). They are: where δ and γ are the liquid/vapour interfacial width and reduced surface tension of the pure solvent defined after eqs. ( 11) and ( 12), while κ D = 8πρ (0) l B is the inverse Debye screening length.

δ = κ D δ; η = ρ (0) 3κ D γ ; K = 4πl b aκ D ( 29
)
When salt is present (ρ (0) = 0), the electrostatic interactions are screened, and the wetting transition of the ionic solution may be either continuous or discontinuous, depending on the specific values of ρ (0) and σ. When ρ (0) = 0, ω e decays exponentially with ζ = κ D ζ. The wetting transition is always discontinuous (as for the salt-free case) as long as ρ (0) < ρ(0) 0.05M [14]. When ρ (0) > ρ(0) , three different wetting scenarios may occur [14]:

(i) Upon increasing the solvent control parameter p 1 , a discontinuous wetting transition occurs at some p tr 1 > p 2 . (ii) As p 1 is increased, a continuous wetting transition occurs at p tr 1 = p 2 . (iii) A succession of two wetting transitions is observed: first a discontinuous jump between a smaller value ζ tr 1 and a larger value ζ tr 2 of the film thickness occurs at some p tr 1 < p 2 ; this discontinuous transition is followed by a continuous divergence of ζ as p 1 → p 2 . Unlike in the case of prewetting, this sequence of wetting transitions is observed along the vapour/liquid coexistence curve, with the continuous transition occuring close to the critical point; this scenario is a direct consequence of the appearance of two minima in the total surface excess grand potential ω as a function of ζ.

For high salt concentrations ρ (0) the electrostatic forces are strongly screened, and do not change the order of the wetting transition of the pure solvent.

The various scenarios may be understood in terms of the competition of two different length scales in the ionic solution -the bulk correlation length δ of the solvent and the Debye screening length κ -1 D . Maps of the occurence of the wetting scenarios in parameter space are given in ref. [14]. In the case of aqueous ionic solutions the wetting behaviour is expected to be very similar to that of pure water due to its very high surface tension γ. Interesting effects may be observable for weakly polar solvents (with γ << γ H 2 O ) or in the presence of organic surfactants.

The semi-primitive model

The phenomenological Cahn-PB theory of wetting by ionic solutions has some obvious limitations: the underlying "primitive" model ignores the coupling between ions and solvent; the square-gradient and PB approximations cannot account for the strong layering of solvent molecules and ions near the substrate, due to excluded volume effects; the free energy functional is not directly related to an underlying microscopic model of the solution, and requires the input of a number of phenomenological parameters, which must be taken from experiment or independent theoretical calculations, e. g. of the liquid/vapour phase diagram. A more microscopic approach, based on a well defined Hamiltonian, is clearly desirable. This would be a formidable task for realistic models of water, but simplified models of the solvent may be expected to provide at least a qualitative picture of the influence of dissolved ions on the wetting and drying behaviour of the solvent, which may be confronted with the trends predicted in Section 4 on the basis of Cahn-PB theory. The "semi-primitive" model (SPM) briefly introduced in Section 3 is a step in that direction. Solvent molecules are taken to be neutral hard spheres of diameter d 0 , while cations and anions are charged hard spheres of diameters d + and d -; the investigations in reference [16] are restricted to monovalent ions (q + = +e, q -= -e), and the diameters taken equal to their Pauling values. This is the basic model examined in ref. [15]. The Coulombic interactions between ions are reduced by a constant permittivity ε associated with the solvent, hence the name of "semi-primitive" model; it takes into account the granularity of the solvent (unlike the "primitive" model), but its polar nature is accounted for only on a macroscopic level, as is the case with the dielectric continuum representation of the "primitive" model.

Ion and solvent density profiles and related interfacial properties of the SPM near a neutral or charged hard wall were calculated in ref. [15] on the basis of a free energy functional detailed below including an approximate Coulomb correlation term [35], for NaCl, KCl and CaCl 2 . For the monovalent ion cases and moderate surface charges σ the correlation term does not significantly modify the predictions based on the mean-field density functional, which yields density profiles in reasonable agreement with Monte Carlo data for a symmetric SPM [15,36].

Wetting and drying phenomena being controlled by cohesive substrate/fluid and fluid/fluid forces, the above SPM must be augmented by the addition of attractive interactions between the solvent molecules and the ions, and between the substrate and the solution. Following the work of Sullivan [START_REF] Sullivan | [END_REF]26] and of Evans and coworkers [21,[27][28][29], Yukawa-like attractions between all species are assumed. The pair potentials are hence chosen of the generic form (α, β = 0, +, -):

v αβ (r) =    ∞; r < d αβ = (d α + d β )/2 qαq β εr -u r exp[-k f (r/d αβ -1)]; r > d αβ ( 30
)
where ε is the (local) dielectric permittivity of the solvent, u is the energy scale (multiplied by unit length) of the attractive interaction, and k f controls its range. For simplicity u and k f are assumed to be the same for all αβ pairs. The external potential exerted by the planar substrate with surface charge σ on particles of species α is of the form:

V α (z) =    ∞; z < 1 2 d α -u w exp {-k w (z/d α -1)} -2πσqα ε(z) z; z > 1 2 d α ( 31
)
where z is the distance from the substrate (xy plane); u w is the energy scale of the non-Coulombic substrate/solution attraction, while k w controls its range. Convenient reduced variables are T * = k B Td 0 /u, ρ * α = ρ α d 3 0 and σ * = σd 2 0 /e. We know from the work of Tarazona nad Evans [27] that the ratio k f /k w controls the order of the wetting transition of the pure solvent.

The total free energy functional per unit area, generalising the pure solvent form ( 14), may be conveniently split into the following contributions:

F[{ρ α (z)}] = F id [{ρ α (z)}] + F HS [{ρ α (z)}] + F Y [{ρ α (z)}] + F el [{ρ α (z)}] (32)
The ideal term F id is the sum of three contributions of the form ( 15). For the hard sphere term we have used the most accurate fundamental measure form adapted to a three-component HS mixture [15,32]. For the Yukawa attraction term F Y and the electrostatic term F el the mean-field form is adopted:

F Y [{ρ α (z)}] = 1 2 α,β ∞ 0 dz ∞ 0 dz ρ α (z)v Y αβ (|z -z |)ρ β (z ) (33)
where, using cylindrical coordinates:

v Y αβ (|z -z |) = 2π 0 dφ ∞ d αβ rdrv Y αβ (| r -r |) (34)
and

F el [{ρ α (z)}] = 1 2 ∞ 0 dzρ c (z)Ψ(z) (35)
where ρ c (z) is the charge density ( 17), while Ψ(z) is the local electrostatic potential, which satisfies Poisson's equation:

d dz ε(z) d dz Ψ(z) = -4πρ c (z) (36)
Short-range correlations due to excluded volume effects are properly included in the functional ( 32) through the fundamental measure HS term. Correlations due to the attractive Yukawa interactions and to the long-range Coulombic forces are neglected in the mean-field approximations ( 33) and ( 35), but it was shown in ref. [15] that this is a reasonable approximation for the Coulombic term, at least for monovalent ions and moderate surface charges σ. Three different choices for the local dielectric permittivity ε(z) have been made in ref. [16]:

• A constant value, equal to the bulk permittivity of the solvent.

• A local Clausius-Mossotti (CM) expression

ε(z) = 1 + 8π 9k B T m 2 ρ0 (z) 1 -4π 9k B T m 2 ρ0 (z) ( 37
)
where m is the assumed dipole moment of the solvent molecules, while ρ0 (z) is the local density of the solvent, averaged over a sphere of diameter d 0 (weighted density).

• A phenomenological, sigmoidal function of the weighted local density, mimicking the permittivity of water [16].

The CM approximation ( 37) results from a mean-field description of the genuine polar HS solvent used in the "civilized" model considered in Section 6. It is afflicted by the familiar "dielectric catastrophe" at high solvent density, and must hence be restricted to sufficiently small dipole moment m. In the calculations of reference [16] the reduced dipole moment m * = m/ (ud 2 0 ) 1/2 was set at m * = 0.75, thus avoiding the "dielectric catastrophe" throughout the fluid range, and yielding dielectric permittivities ε between 5 and 10 in the dense liquid phase (weakly polar solvent). For a given ε(z) Poisson's equation ( 36) was solved numerically by a Runge-Kutta method, subject to the boundary conditions: The equilibrium density profiles are determined by minimizing the grand potential per unit area:

dΨ(z) dz z=0 = - 4πσ ε(0) , dΨ(z) dz z→∞ ---→ 0 (38)
Ω[{ρ α (z)}] = F[{ρ α (z)}] + α ∞ 0 dz [V α (z) -µ α ] ρ α (z) (39)
with respect to the ρ α (z) for fixed chemical potentials µ α . The surface tension γ is finally given by the three-component generalisation of eq. ( 6). The liquid/vapour binodal line of the pure solvent and at finite salt concentrations is calculated from the bulk limit ρ α (z) → ρ ((0) α of the free energy ( 32). Note that in the mean-field approximation F el vanishes in that limit due to global charge neutrality. The effect of the ions on the phase diagram is hence only due to their contributions to the other terms in the free energy. Up to ion concentrations of 2M this effect is small, and only leads to a modest decrease (1%) of the critical temperature T c [16]. Inclusion of correlation terms in the Coulombic contributions to the free energy, e.g. within the mean spherical approximation [37], might lead to more substantial changes in the phase diagram, but would spoil thermodynamic self-consistency at the level of the Gibbs adsorption equation ( 19).

The equilibrium density profiles of the liquid film which forms upon approaching the binodal from the vapour side along an isotherm exhibits strong layering, as illustrated in Fig. 1. As expected the film thickness grows as coexistence is approached, and diverges above a wetting temperature T w . A second order scenario is observed when the ranges of the wall/fluid and fluid/fluid attractions are comparable. The case k w = k f = 1.8 is illustrated in Fig. 2. The distance from coexistence is measured as a difference in the Gibbs function g = α x α µ α (where x α is the molar fraction of species α) between undersaturated and saturated vapour. When the range of the wall/fluid attraction increases relative to the fluid/fluid attraction, the wetting transition goes over to a discontinuous scenario: prewetting is observed off coexistence, followed by a continuous divergence of the film thickness as coexistence is approached; the case with k w = 1.2, k f = 1.8 is illustrated in Fig. 3.

Concerning the key question of the present investigation, namely the change in the wetting and drying behaviour of the pure solvent induced by the ionic solute, the main conclusions of ref. [16] may be summarised as follows, in terms of two additional physical variables, the salt concentration c and the surface charge density σ. Higher c results in stronger screening of the electrostatic forces, rendering wetting a priori less favourable. The addition of salt is indeed found to generally increase the wetting temperature T w compared to that of the pure solvent; the effect increases with c and can be as large as 4% for c=1M, depending on the assumed permittivity ε. On the other hand, for a given c, T w is found to decrease linearly with |σ| if the wetting transition of the pure solvent is first order. Conversely, the variation of T w is non-monotonic, and depends moreover on the sign of σ when the pure solvent undergoes a second order wetting transition.

The order of the wetting transition in most ionic solutions is predicted to be the same as in the corresponding solute-free fluids. However, in concentrated solutions (short Debye length), where the excess of counterions extends only over roughly one layer, a transition which was continuous in the ion-free solvent may become first order in contact with a highly charged substrate. By contrast, at low salt concentrations (large Debye length), where charge separation near the wall can span the entire width of the liquid film, discontinuous wetting is favoured even for moderate values of |σ|. These predictions are in overall qualitative agreement with those of the Cahn-PB theory of Section 4. Quantitative comparison is difficult because of the phenomenological nature of the latter, and because the primitive model cannot account for layering effects, which play a subtle role in determining the interfacial behaviour predicted by the explicit solvent model. The work based on the SPM found no evidence of the sequence of two wetting transitions along the coexistence curve predicted in ref. [14] (Section 4) for some ionic solutions in the region of parameter space adjacent to that where the crossover from a discontinuous to a continuous transition takes place. The scenario cannot be completely ruled out, since the fraction of parameter space explored in ref. [16] is much smaller than in ref. [14] due to the much larger computational effort required with the explicit solvent model.

When the liquid/vapour coexistence curve is approached from the liquid side along an isotherm, drying (or dewetting) is commonly observed near a solid substrate if u w = 0. The solute-free HCY solvent was shown to dry completely in contact with a hard wall, i. e. an infinitely thick vapour layer intrudes between the liquid and the substrate for any temperature between the triple and critical points [16]. When ions and surface charge on the substrate are added, the drying behaviour of the system changes radically due to the Coulombic interaction between the electric double layers which form at the fluid/substrate and liquid/vapour interfaces. The monovalent counterions are strongly attracted to the oppositely charged substrate, where they form an adsorbed layer, which only partially neutralizes the surface charge, while the Cl -ions outnumber the Na + cations at the liquid/vapour interface, which is hence negatively charged, in agreement with experimental observations [38] and molecular dynamics simulations [39,40] of aqueous solutions. The situation is illustrated by the density profiles calculated at liquid/vapour coexistence for σ > 0 (Fig. 4) and for σ < 0 (Fig. 5). When σ > 0 the two interfaces attract each other because the electric double-layers carry opposite charges; this attraction hinders complete drying which is only observed when σ → 0. When σ < 0, cations are preferentially adsorbed on the substrate, but the corresponding double-layer charge remains globally negative, thus repelling the negatively charged liquid/vapour interface, and favouring drying. Complete drying is indeed predicted for |σ| less than a threshold value (of the order of 1e/nm 2 , the precise value depends on the concentration and temperature of the solution). But drying can be hindered if |σ| is sufficiently large for the Yukawa attraction between the adsorbed cations and the solvent to outweight the solvent-solvent cohesion forces and the electrostatic repulsion between the two interfaces. For a fixed moderate value of |σ| the calculations of ref. [16] predict a first order drying transition between a partially wet wall (finite film thickness) and a completely dry wall above a drying temperature T d along the liquid/vapour coexistence curve. Contrary to the conjugate wetting temperature T w , T d is found to increase with |σ|. For very low values of |σ|, T d lies below the triple point, and the wall is completely dry for any thermodynamic state on the binodal. Conversely, if the substrate carries a very large negative charge density T d > T c and the wall is always partially wet. Analogously to wetting, pre-drying is observed upon approaching the liquid/vapour coexistence curve along an isotherm above T d . The first order drying transition is predicted to occur in the absence of any wall/liquid dispersion attraction. It would be interesting to investigate the competition between substrate/fluid Coulombic and dispersion forces.

6 The "civilized" model.

In a further step towards physical reality the "civilized" model, featuring a dipolar solvent was used in ref. [17] to investigate the wetting and drying behaviour of ionic solutions. The solvent molecules are HCY spheres carrying a dipole moment m 0 . Charge-dipole and dipole-dipole contributions must be added to the pair potentials ( 30), namely:

v d αβ ( r, m α , m β ) = 1 r 3 q α m β • r - m α m β r 3 [3( mα • r)( mβ • r) -mα • mβ ] (40)
where m α = 0 for α = 0, while q α = 0 for α = 0. The external potential acting on the polar molecules is now:

V 0 (z, θ) =    ∞, z < 1 2 d 0 -u w exp {-k w (z/d 0 -1)} -2πσm 0 cos(θ) (41)
where θ is the polar angle of m 0 relative to the z axis (normal to the substrate).

The density profile of the solvent depends on the local orientation of the dipole:

ρ(z, θ) = ρ(z)α(z, θ) (42)
where α(z, θ) may be expanded in Legendre polynomials:

α(z, θ) = ∞ l=0 α l (z)P l (cos(θ)) (43)
with α 0 (z) = 1/4π for proper normalization. Consider first the solute-free pure solvent, placed inside a condenser of oppositely charged plates (±σ, to ensure overall charge neutrality in the absence of ions). The free energy functional F [ρ 0 (z, θ)] may be split into ideal, hard sphere, Yukawa and electrostatic contributions, as in eq. ( 32), with (temporarily dropping the subscript 0)

F id [ρ(z, θ)] k B T = ∞ 0 dzρ(z) log Λ 3 ρ(z) -1 + 2π ∞ 0 dzρ(z) π 0 dθ sin(θ)α(z, θ) log [4πα(z, θ)] (44)
The hard sphere contribution F HS is given by the one-component version of the Rosenfeld functional [30,31], the Yukawa contribution is given by the one-component version of ( 33), while the electrostatic contribution is, within the mean-field approximation:

F el [ρ(z, θ)] = - 1 2 ∞ 0 P(z)E dd (z) dz (45)
where P(z) is the local polarization density:

P(z) = 2π π 0 m cos(θ)ρ(z, θ) sin(θ)dθ = 4πm 3 ρ(z)α 1 (z) (46)
while E dd (z) is the local electric field generated by the dipolar particles: The orientation-independent density ρ(z) and the orientational order parameter α 1 (z) (and hence the polarization ( 46)) are determined self-consistently from the Euler-Lagrange equation derived from the minimization of the grand potential functional ( 3) of ρ(z, θ) [17]. Due to the mean-field approximation for F el the orientational order parameter takes the classic Langevin form:

E dd (z) = 2π ∞ 0 dz |z -z | 2 d 3 - 1 d P(z )Θ H (d -|z -z |) (47)
α 1 (z) = 3 4π coth(βmE(z)) - 1 βmE(z) (48)
where E(z) = E 0 + E dd (z) = 2πσ + E dd (z) is the total electric field acting on the dipoles. For the same reason the dielectric permittivity of the bulk dipolar solvent reduces to the classic Clausius-Mossotti form, the local version ( 37) of which was used in the SPM description of Section 5. To avoid the "dielectric catastrophe", all DFT calculations are restricted to a reduced dipole moment m * = m/ud 3 ≤ 0.9.

For σ = 0 the model reduces to the HCY solvent model of Section 5; when σ = 0 the electric field partly orients the dipole moments, and this leads to a drop of the critical temperature T c proportional to m * [17]. Note that if dipole-dipole correlations are taken into account [41,42], T c of the nonelectrified (σ = 0) polar solvent increases with m * . Detailed DFT calculations [17] show that the wetting scenarios of the dipolar solvent are very similar to those of the HCY solvent (m * = 0). The wetting temperature T w drops as m * increases from 0 (HCY limit) for fixed values of σ, but the ratio T w /T c is practically independent of σ and m * , both for first and second order transitions.

When ions are added the DFT calculations very much procede as for the SPM in Section 5, except that the local electric potential Ψ(z) now satisfies Poisson's equation:

d 2 Ψ(z) dz 2 = -4πρ c (z) + 4π dP(z) dz (49)
rather than eq. ( 36). As illustrated in Fig. 6, the density profiles calculated for the "civilized" model are practically indistinguishable from those obtained for the SPM under identical physical conditions (T, c, m * ), provided the local CM permittivity ( 37) is used in the SPM calculations. This remarkable agreement seems to be the first justification of the use of a local dielectric permittivity ε(z) in the description of inhomogeneous ionic solutions (in particular in explicit solvent descriptions of electric double layers), at least within a mean-field context. A direct consequence of this close agreement between the density profiles, and hence the adsorptions Γ α predicted by the two models, is that the wetting and drying scenarios of the "civilized" model are qualitatively and quantitatively identical to those described in Section 5 for the SPM. In particular the novel drying transitions predicted in Section 5 are confirmed by the DFT calculations based on the "civilized" model [17].

Surface tension.

In this section we focus on the surface tension of liquid/vapour and solid/fluid interfaces of the threecomponent explicit solvent models of ionic solutions discussed above. In particular, we investigate the influence of ion concentration and distribution (electric double-layers) on interfacial tension. Within DFT the surface tension is a functional of the equilibrium density profiles of the species present in the system. A generalization of eq.( 6) appriopriate to describe the free liquid/vapour interface is:

γ lv = F [{ρ α (z)}] -F ρ SK α (z) - α ∞ -∞ µ α ρ α (z) -ρ SK α (z) dz (50)
where ρ SK α (z) are sharp-kink density profiles:

ρ SK α (z) = ρ SK α (z) = ρ v α , z < 0 ρ SK α (z) = ρ l α , z > 0 (51)
while ρ v α and ρ l α are the densities of species α in the vapour and in the liquid, respectively, which coexist at a given temperature. The three sharp-kink profiles form a reference set against which the excess grand potential per unit area is measured. The same method was used by Groh and coworkers [43] to calculate the liquid/vapour surface tension of the restricted primitive model of an ionic solution. Note that due to the discontinuity of the auxiliary sharp kink profiles in eq. ( 51), the associated free energy F ρ SK α (z) is discontinuous at z=0, and takes on bulk values for z<0 (vapour) and for z>0 (liquid). We have calculated the variation of the surface tension of the free liquid/vapour interface with ion concentration for several temperatures. A representative sample of the results we have obtained is shown in Fig. 7. γ lv of the model solutions is significantly lower than that of real aqueous solutions at room temperature (about 70 mN/m). Such a result is not surprising, as our model does not account for the complex network of hydrogen bonds which form between water molecules, and which introduce a substantial contribution to its liquid/vapour surface tension. Contrary to the behaviour observed for alkali halides in water, our calculations show a slow but systematic decrease of γ lv with ion concentration (as seen in the inset of Fig. 7). This trend was found to be true for both "semi-primitive" and "civilized" models. The main reason behind this disagreement is the poor decription of electrostatic interactions in the bulk fluid within the mean-field approximation. The electrostatic contribution to the uniform fluid's free energy vanishes identically, resulting in a phase diagram where the coexisting phases' total densities are essentially independent of ion concentration (see Fig. 1 of ref. [16]), unlike in real ionic solutions, where the total density of the liquid phase increases with c, leading to an increase in γ lv . The actual decrease of γ lv (c) for our model solution is most likely related to the accumulation of chlorine anions (the largest species in the system) in the interfacial region. Such an entropic effect has been observed in dilute electrolytes, where the change in the total density of the solution is minor. We have also investigated the solid/liquid surface tension γ sl . This quantity is related to the inhomogeneous density profiles of an ionic solution under coexistence conditions in contact with a substrate by the following equation:

γ sl = F [{ρ α (z)}] -F ρ l α + α ∞ 0 V α (z)ρ α (z) -µ α ρ α (z) -ρ l α dz (52)
As for the free liquid/vapour interface, we analysed the influence of ion concentration on surface tension, and compared γ sl (c) curves obtained for different dielectric permittivity models and substratefluid interaction parameters. Unlike γ lv (c), for the vast majority of systems investigated γ sl (c) is not a monotonic function and has a minimum. The location of this minimum depends on both the dielectric properties of the solvent and the substrate-fluid potential. The typical variation of γ sl with c is illustrated in Fig. 8. However, note that even if no minimum was found within the concentration range investigated, the system in question most likely has one at a lower c, since for very dilute solutions the amount of counterions needed to neutralise the substrate's charge is very high in comparison with this species' bulk density, which introduces a large contribution to the surface tension. Complete wetting occurs if the condition that γ sv = γ lv + γ sl is fulfilled. Therefore, our latest results for γ lv and γ sl indicate that the wetting transition temperature of an ionic solution near a charged substrate is generally a non-monotonic function of solute concentration.

Conclusions.

We have reviewed and supplemented recent theoretical work investigating the influence of dissolved monovalent ions on the wetting and drying scenarios of model solvents near a charged solid substrate. The earlier approach is based on the primitive model representation of ionic solutions, and combines Cahn's square gradient DFT for the solvent with non-linear Poisson-Boltzmann theory of the ionic double-layers. Subsequent approaches use explicit solvent representations, and build on the classic DFT theory of Evans and collaborators for the wetting and drying behaviour of simple fluids. Ions have a small but significant effect on the wetting behaviour of ionic solutions, leading in particular to a shift of the wetting temperature. Close to the physical conditions where a cross-over from a first to a second order wetting transition is observed for the pure solvent, the addition of ions can change the order of the transition.

The predictions based on the "semi-primitive" model (with a local Clausius-Mossotti dielectric permittivity) are very close to those obtained with the "civilized" model using a mean-field description of the electrostatic interactions between ions and dipoles. Both models predict a novel drying transition which is completely absent in the pure solvent.

The new results for the surface tension γ lv of the liquid-vapour interfaces predict a slow decrease of γ lv with increasing ion concentration c, and a sharp drop with increasing temperature. The surface tension of the substrate-liquid interface, on the other hand, is found to have a minimum at a certain c min . The variation pattern of the surface tension associated with these two types of interfaces indicates that wetting temperatures of ionic solutions are typically non-monotonic functions of the solute concentration c. The decrease of γ lv as c increases is opposite to the trend observed experimentally for aqueous solutions, pointing to the expected failure of the "semi-primitive" and "civilized" solvent models to provide a realistic description of water. A DFT of wetting by aqueous ionic solutions, based on a realistic representation of water, e.g. the SPC/E [44] or the TIP4P [45] models, is not yet available, and should be a goal for future investigations. Another important objective for future work is an extension of the present DFT and of time-dependent density-functional theory (TDFT) [46] to explore the technologically important electro-wetting phenomena which are well documented experimentally [47].

Page 6 of 27 URL

 27

Θ

  H being the Heaviside step function.

Figure 1 :Figure 2 : 20 PageFigure 3 :

 12203 Figure 1: Density profiles ρ α (z) formed by NaCl vapour in contact with a charged substrate, illustrating the liquid film growth as coexistence conditions are approached along an isotherm. The saturated vapour coexists with a 0.1M solution. The dashed lines correspond to ρ (0) α = 0.99ρ co α , dotted lines to ρ (0) α = 0.999ρ co α and solid lines to ρ (0) α = 0.9999ρ co α .Na + (counterion) profiles are marked with triangles, Cl -(coion) profiles with circles, and H 2 O profiles with squares. ε(z) is approximated using a sigmoidal function mimicking the permittivity of water (eq. (34) from ref.[16]), k w = k f = 1.8, σ * = -0.022, T * = 0.736.

Figure 4 :

 4 Figure 4: Density profiles of 1M NaCl at liquid-vapour coexistence near a hard wall carrying a surface charge σ * = +0.0087. T * = 0.75, ε(z)is the CM function defined by( 37). The inset shows the total amount of charge per unit area I c (z) contained between 0 and z (including the substrate charge).

Figure 5 : 23 Page 24 of 27 URLFigure 6 :

 523276 Figure 5: As in Fig. 4, but for σ * = -0.072.

Figure 7 :

 7 Figure 7: The liquid-vapour surface tension of SPM NaCl solutions as a function of ion concentration. The inset shows a single curve corresponding to T * = 0.75 to illustrate the steady decrease of γ lv with c. ε(z) is approximated using a sigmoidal function mimicking the permittivity of water.

Figure 8 :

 8 Figure 8: The liquid-vapour surface tension of SPM NaCl solutions as a function of ion concentration. For all systems illustrated in this figure k w = 1.2, and T * = 0.77. The line marked with circles corresponds to σ * = -0.022 and sigmoidal ε(z) as in Fig.1(values of ε in the liquid are close to 80, as in real water), the line marked with diamonds to identical σ * , but ε(z) given by the CM formula (with m * = 0.75, which yields liquid phase permittivity about 10 times smaller than the sigmoidal formula), and the line marked with squares (shown in the inset for clarity) to σ * = -0.043 and sigmoidal ε(z).

Page 26 of 27 URL

				Molecular Physics
		28.3
		28.25
		28.2
	γ sl [mN/m]	28.05 28.1 28.15				24.68 24.7 24.72 24.74
		28				24.66
		27.95				24.62 24.64
		27.9					0.1 0.2 0.3 0.4 0.5 0.6
	27.85 r F o	0.1	0.2	0.3	c [M]	0.4	0.5	0.6
		P
		e
			e r
				R
				e
					v i e
							w
								O n l
									y

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

Page

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

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

Acknowledgements

We are indebted to Bob Evans and Roland Roth for their interest in our work. JPH would like to thank Bob Evans for constant inspiration over more than 30 years.


Richard James Matthews
email: r.matthews1@physics.ox.ac.uk

Ard A Louis

Julia M Yeomans

Confinement of knotted polymers in a slit

Keywords: knotted polymers, confinement, Langevin dynamics knotted polymers, confinement, Langevin dynamics

come

Introduction

Whilst knots are generic, and any polymer of sufficient length in solution is very likely to be knotted [1], they have proven to have particular relevance in biological systems. The DNA contained in the cells of living organisms is typically both long and confined: the genome of E. coli, a well studied bacterium [2], is about 3 × 10 4 persistence lengths [2] (∼ 1.5mm) and is contained in a cell with a typical volume ∼ µm 3 [3]. Although the DNA is organized using supercoiling and proteins (particularly in eukaryotes) [4], confinement makes entanglements likely [5,6], particularly during processes such as replication [7]. Knots in DNA can prevent transcription [8], a process necessary for the production of proteins. Indeed, there is a family of enzymes, the topoisomerases [9], one of whose functions is to control knotting. Evidence of highly confined DNA knotting has also been found in bacteriophages [START_REF] Arsuaga | Proc. Nat. Acad. Sci[END_REF], viruses that infect bacteria. These knots may affect the ejection speed of a bacteriophage's DNA [START_REF] Matthews | [END_REF][START_REF] Marenduzzo | Proc. Nat. Acad. Sci[END_REF]. Given its relevance in biology, it is important to understand the interplay of knotting and confinement.

The confinement of linear polymers in a slit, between two infinite, parallel walls, has been well studied. The simplest model is to consider a single chain with walls that provide only geometric constraints. For D > R, where D and R are slit and polymer size respectively, the chain will maintain a shape similar to that seen without confinement. For D < R, however, significant deformation occurs. A scaling form for the free energy of such a chain for D < R was given by Daoud and de Gennes [START_REF] Daoud | [END_REF] using a blob picture, in which the polymer is taken to form a series of independents blobs of size determined by D, within each of which the polymer behaves as if it is free. They predicted

F conf k B T ≅ N a D 1/ν , ( 1
)
where N is the number of monomers in the chain, a is monomer size and ν is a three dimensional exponent ≈ 3/5. Predictions made by Daoud and de Gennes for the polymer size were verified by Webman, Lebowitz and Kalos using simulations [14]. Later, it was predicted that the ratio between the monomer density by the wall and the force exerted on the wall is universal, independent of microscopic details [15]. Again, attempts have been made to verify this prediction by simulation, see ref. [16] and references therein. Whilst the majority of simulation work has used Monte Carlo methods, a recent study has also compared the results of Molecular Dynamics simulations to scaling predictions [17]. Other work has considered issues such as attractive walls [18]. By comparison to the case of linear polymers, there has been little work on knotted polymers in slits and basic issues remain open. Tesi et al. [6] looked at the probability that a knot will be found in a slit-confined polymer. They found that for relatively wide slits the probability increased as the slit width decreased. However, for narrower slits, the probability peaked and eventually became less than the value for wide slits. Janse van Rensburg [19] applied both analytical calculations and a Monte Carlo approach that samples configurations with the same topology but varying N . The change of the average value of N with slit width was investigated. Knotted polymers were found to expand as the width increased, in contrast to unknotted polymers whose size showed a plateau after a certain width was reached. In ref. [20], the ends of knotted and unknotted polymers were attached to parallel plates, separated so as to stretch the chains. It was found that the forces exerted on knotted polymers deviated from scaling predictions. The results were used to deduce an estimate of how the knot size scales with polymer length.

Motivated by the common occurrence of confined knotted polymers in biology, we use a coarse-grained model to gain insight into the behaviour of knotted ring polymers in a slit, concentrating on the force exerted on the walls. We then go on to look at other quantities, such as the radius of gyration and the monomer density distribution, that help us to interpret the trends seen.

We find that the effect of knotting differs depending on the degree of confinement of the polymer. We find that for narrow slits, more complex knots exert higher forces on the walls, whereas for relatively wide slits, the opposite is true. We relate the forces to the monomer densities near the walls and interpret the results in terms of the effect of knots on the ability of the polymer to spread out in the slit.

Model

We use a standard, bead-spring, polymer model [21]. Excluded volume is included through a truncated Lennard-Jones potential: where σ is the size of a bead, ǫ is an energy scale and x is the bead separation.

V EV (x) = 4ǫ σ x 12 -σ x 6 + 1 4 , x ≤ 2 1/6 σ, V EV (x) = 0, x > 2 1/6 σ, (2)
Neighbouring beads are connected with FENE springs [22] V

F EN E = - kX 2 0 2 ln 1 - x X 0 2 , (3)
where X 0 , the maximum extension, is set to 1.5σ and k is chosen as 30ǫ/σ 2 . The polymer motion is simulated using the Langevin equation [23]:

mẍ i = f i -ζ ẋi + 2ζk B T r i (t), (4)
where x i is the position of the i th bead, f i is the total force acting on it, m is the bead mass, ζ is a friction constant and k B T is the temperature, which we set to ǫ. r i (t) are random vectors that satisfy

r i (t) = 0, r i (t)r j (t ′ ) = δ ij δ(t -t ′ )I, ( 5
)
where I is the identity matrix. The second and third terms in Eq. ( 4) represent the drag due to the solvent and Brownian noise respectively. The parameters may be combined to give a convenient unit of simulation time, t 0 = (σ 2 m/ǫ). We integrate the motion using a velocity Verlet algorithm [START_REF] Frenkel | Understanding Molecular Simulation: from Algorithms to Applications[END_REF] with a timestep of ∆t = 0.01. The polymer is placed in a slit between two parallel plane walls, modelled by the same excluded volume potential as for the bead-bead interactions, Eq. ( 2). The surfaces of the walls are taken as lying at the point at which the potential is truncated. We consider slit widths from D = 4σ to D = 18σ. This range was chosen, through initial exploratory simulations, to show the crossover between the regimes where more complex knots are harder or easier to confine. The simulations resemble those of Dimitrov et al. [17], who considered confined, linear polymers.

We choose a polymer length, N = 300, which is sufficiently short that the knots will significantly affect the polymer properties. It should be borne in mind that, since knots in polymers in a good solvent are weakly localized [START_REF] Orlandini | [END_REF], different results would be obtained for other values of N . The results do, however, give a qualitative picture of how a knot, which is of signif icant length compared to N , will affect a polymer in a slit.

We consider linear and unknotted ring polymers (denoted 0 1 ) and also polymers with knot types 3 1 , 6 1 and 9 1 . C p is a standard notation [START_REF] Orlandini | [END_REF], where C gives the minimum number of crossings in any projection onto a plane, the essential crossings. p is used to distinguish knots with equal C. 3 1 and 9 1 are torus knots; they both belong to a group of knots with similar topology. 6 1 is an even-twist knot. In addition, we consider two other knots. Firstly, we include a twelve-crossing even-twist knot. We label it 12 1 , although the standard notation only extends up to C = 10. Secondly, we also look at 8 18 , which is chosen because in initial simulations it was found to have an unusually spherical average instantaneous shape [28]. In Fig. 1. we show a knot diagram [START_REF] Livingston | Knot Theory[END_REF] for 8 18 . The polymer model used should prevent chain crossings [21] and we verified that the knot type remained constant by calculating the Alexander polynomial [START_REF] Orlandini | [END_REF], A k (t), at t = -2.

The radius of gyration of the unconfined linear polymer, R G (linear) = 14.2±0.1σ. For the ring polymers with different knots, R G (0 1 ) = 10.65 ± 0.03σ, R G (3 1 ) = 9.01±0.02σ, R G (6 1 ) = 7.78±0.01σ, R G (8 18 ) = 6.86±0.01σ, R G (9 1 ) = 7.28±0.01σ and R G (12 1 ) = 6.90 ± 0.01σ. However, it should be kept in mind that the instantaneous configurations of linear polymers are not spherical, but rather one can define 3 axes that characterize a prolate shape [28]. We investigated the anisotropic instantaneous distribution of knotted polymer shapes, which tend to be more spherical than linear polymers [29], although here there is no universal behaviour because it depends on the size of the knot compared to the length of the polymer.

We note that for some of the more complex knots, an increase in the internal energy of the chain larger than the standard error was observed: for narrow slits and complex knots the probability of collision increases and so the amount of overlap between beads is higher, leading to an internal energy increase due to the excluded volume interaction. The largest such increase, ∆U = 0.56 ± 0.06k B T , was seen for 12 1 between D = 8σ and D = 4σ. In the middle of this range, the force on the confining walls was 7.734 ± 0.003k B T /σ. This suggests that the internal energy change accounts for 1% of the free energy increase of the chains, with the rest being due to entropy. We are thus confident that this effect, which depends on the specific potential chosen, does not affect our conclusions, which are qualitative.

The starting configurations for the polymers were created by performing runs in which D was reduced at a rate of 10 -4 σ/t 0 from a wide slit, D = 25σ, to the narrowest slit we considered, D = 4σ. The final configuration was then taken as the starting configuration for all the simulations at different D for that particular knot type. By performing extra compression runs, we checked that the configurations used were not unusual: all those seen at end of the additional compression runs were qualitatively similar. When viewed from the direction across the slit, they showed a number of chain crossings that was typically equal to, or only slightly larger than, the number of essential crossings of the knot.

Simulations for each D were performed separately. Starting from a polymer in the configuration obtained from the compression run, an initialisation period of 10 5 t 0 was allowed before results were recorded for 10 7 t 0 . Relaxation was slower in the two free directions than the confined direction and also depended on D and knot type [30]: depending on the particular system, the timescale may vary with knot complexity in different ways [31]. The relaxation time in the free directions, taken as the time for the auto-correlation function of the relevant component of the radius of gyration to decay to 1/e, was typically on the order of 10 2 t 0 for the knotted polymers and even for the slowest case (D = 4σ, linear polymer) was clearly less than 10 4 t 0 .

Errors were estimated from the variance of the relevant quantities, initially assuming that all data points were independent. During the simulations, data was recorded every 100t 0 . If the estimate of the relevant correlation time, t C , was seen to be significantly greater than 100t 0 then the estimated error was increased by a factor of 2t C /100t 0 [START_REF] Frenkel | Understanding Molecular Simulation: from Algorithms to Applications[END_REF].

Results

We study the average force, f , exerted on the confining walls by the polymer as function of D for different knot types. The average force for knotted and unknotted rings is plotted as a ratio to the value for a linear chain in Fig. 2.

We first consider the forces exerted by linear and unknotted ring polymers, which are very similar over the range of slit widths investigated. The scaling picture of a confined polymer is based on it forming a series of blobs of dimension ∼ D [START_REF] Daoud | [END_REF]. For small widths, both linear and ring polymers should form the same number of blobs. However, since linear polymers have a larger radius of gyration, it would be expected that, for wide enough slits (D ≫ R G ), the forces exerted by the unknotted ring would decrease below those for linear polymers. At D = 30σ, the component of the radius of gyration in the direction of confinement was 5.24 ± 0.01σ for a linear polymer, compared to 5.14 ± 0.01σ for an unknotted ring, and the force exerted by the unknotted polymer was 0.886 ± 0.006 of the linear polymer value. The difference reflects the fact that R G (0 1 ) < R G (linear).

The topological constraints of a knot that is relatively large compared to the polymer prevent the formation of independent blobs so it is expected that results for knotted rings will deviate from those for linear and unknotted polymers. Fig. 2 shows that, for narrow slits, polymers with more complex knots exert a higher force, whilst for wide slits the opposite is true. The crossover occurs at smaller slit widths for more complex knots, likely because they are more compact [30]. 8 18 , picked because it is unusually spherical, shows, uniquely amongst the knots studied, non-monotonic behaviour for the force ratio as a function of D. The force itself always showed monotonic behaviour for each knot type. We conjecture that the non-monotonicity in the ratio for 8 18 reflects a qualitatively different behaviour for this knot compared to the others due to its more spherical shape. A polymer which typically has an extended shape may respond to increasing confinement by losing rotational degrees of freedom for its overall configuration. A highly spherical polymer, such as one with an 8 18 knot, on the other hand, cannot do this and must deform. This leads to a more rapid increase in the force at relatively large D.

Although the force exerted by 8 18 is the highest at intermediate slit widths, for narrower D it is very close to 6 1 , a knot with less essential crossings. Once the polymer is significantly deformed, its unconfined shape would not necessarily be expected to further influence its behaviour. Polymers in narrower slits should show behaviour like two-dimensional, absorbed polymers and, as shown in ref. [START_REF] Orlandini | [END_REF] for example, knots in such polymers show quite different behaviour to those in three-dimensional polymers.

We next compare the radius of gyration of the polymers, focussing on the com- , where z is the distance from the centre of the slit and δ is a correction necessary to achieve agreement with scaling formulae when D is comparable to a [16,17]. (a) Monomer density across the slit for D = 4σ ( ), D = 12σ ( ) and D = 18σ (×) for linear polymers. The solid line is a scaling formula [16,17]. We use a scaling exponent ν = 0.588 [17] ponent in the confined direction, R perp . In Fig. 3, results are plotted as ratios to the values for a linear polymer. The results for R perp show very similar trends to those seen for the forces: as may be appreciated by comparing Figs. 2 and3, the ordering of the different knot types is almost identical for all D.

To help interpret the results for f and R perp , we examine the monomer density as a function of position across the slit for various D. The results are shown in Fig. 4. Firstly, in Fig. 4(a), we compare the results for linear polymers to the following scaling formula for monomer density [16,17]:

ρ(ζ) = Γ(2 + 2/ν) (Γ(1 + 1/ν)) 2 (1/4 -ζ 2 ) 1/ν , (6)
where ζ is the position across the slit divided by D + δa. δ is a correction, necessary when D is comparable to a [16,17]. Assuming a scaling exponent, ν = 0.588 [17], we find the best fit to the data using δ = 1.38. The results show reasonable agreement with the scaling formula, although clear differences are seen, particularly for D = 4σ. Given that, for this width, D is only a few times bigger than the bead size, it is unsurprising that such deviations are seen. We also checked our results by comparing the forces for linear polymers to a scaling formula derived from Eq. ( 1).

As in the work of Dimitrov et al. [17], we found good agreement using the same value of δ as for the monomer density. In Figs. 4(b)-(d), we compare monomer densities across the slit for different knot types for D = 4σ, 12σ and 18σ, by plotting them as ratios to the results for linear polymers. These simulations show that knots change how the polymer is distributed across the slit in a different way for different slit widths. For wide slits, D 14σ, the constraints of more complex knots make the polymer more compact [30] leading to a lower monomer density at the walls. For narrower slits, D 8σ, however, the more complex the knot the higher the wall density. The value of the monomer density is expected to control the force exerted by the polymer and comparison of Figs. 2 and4 show a close correlation between the two quantities.

Similarly, the trends for the radius of gyration across the slit shown in Fig. 3 can be explained in terms of the results in Fig. 4: for narrower slits the monomers tend to lie closer to the walls for more complex knots leading to an increase in R perp , whereas for wider slits the reverse is true.

Although it is intuitively reasonable that for wider slits more complex knots lead to a more compact polymer and hence a lower monomer density at the walls, it is less obvious why the monomer density at the walls should increase with knot complexity in narrower slits. A suggestion as to why this may be the case is found by looking at the distribution of (instantaneous) R para , the radius of gyration in the two unconfined directions, for a slit with D = 18σ, the widest we consider. As shown in Fig. 5, for more complex knots the distribution is narrower and has a maximum at lower R para indicating that the more complex knots prevent the polymer from exploring configurations where they are spread out along the slit. Hence this suggests that when the slit becomes narrow and significantly distorts the polymer, the polymer cannot so easily distribute its monomers further out in the free directions and instead their density near the wall increases. As shown in Fig. 6 the increase in the average R para as D is changed from 18σ to 4σ is less for more complex knots.

Summary and Discussion

To summarize, we have investigated the effect of knot type on the properties of a ring polymer confined in a slit of width D. We found that, for wide slits, more complex knots exert lower forces than a linear polymer on the walls. For more narrow slits the opposite occurs: knotting leads to larger forces. The forces are seen to correlate with the monomer densities near the walls. The crossover occurs at different slit widths for different knots. For larger slit width, when the polymers are not significantly deformed, those with more complex knots have a smaller radius of gyration [30], reducing the monomer density by the walls and correspondingly the force exerted on them. For smaller slits, the knot prevents the polymer from deforming as much in the free directions as an unknotted polymer can, and so the monomer densities, and consequently forces, are increased. The example of 8 18 , which has an unusual spherical shape and shows nonmonotonic behaviour of the force ratio (to that for a linear polymer) as a function of D, demonstrates that the particular topology of the knot plays a role. We conjecture that the maximum for 8 18 occurs because, due to their shape, polymers with 8 18 cannot respond to decreasing D by losing rotational degrees of freedom for their overall configuration, as more elongated polymers can. Instead, they begin to deform at a point, relative to their size, where other polymers still maintain their preferred shape, if with reduced orientational freedom.

Although we exclusively consider one polymer length, the trends we have uncovered should be generic: except in the large N limit, knots make polymers more compact and will decrease their ability to spread out as they are confined. Therefore we expect that there is a wide range of polymer lengths for which similar qualitative trends will be observed. An interesting issue to be explored in the future is the crossover from weak to strong knot localization as the polymer is confined: polymers in a good solvent in three-dimensions are weakly localized [START_REF] Orlandini | [END_REF] whilst those in two-dimensions are strongly localized. It would also be useful to look at the crossover from the type of behaviour we observe to the large N limit, where, due to localization, knot type will not affect polymer properties.

  , A.A. Louis and J.M. Yeomans

3 Figure 1 .

 31 Figure 1. Knot topology for 8 18 , based on a diagram from ref. [26].

Figure 2 .

 2 Figure 2. The average force, f , exerted on the walls of a confining slit by polymers with different knot types -0 1 (×), 3 1 ( ), 6 1 (△), 8 18 ( ), 9 1 (+) and 12 1 (♦) -divided by the values for a linear chain, as a function of slit width, D. Errorbars are plotted but for most points they are much smaller than the symbols.

Figure 3 .Figure 4 .

 34 Figure 3. The radius of gyration, in the confined direction, Rperp, of polymers in a slit with different knot types -0 1 (×), 3 1 ( ), 6 1 (△), 8 18 ( ), 9 1 (+) and 12 1 (♦) -divided by the values for a linear chain, as a function of slit width, D.

  Figure 4. Monomer densities as a function of ζ = z/(D + δa), where z is the distance from the centre of the slit and δ is a correction necessary to achieve agreement with scaling formulae when D is comparable to a[16,17]. (a) Monomer density across the slit for D = 4σ ( ), D = 12σ ( ) and D = 18σ (×) for linear polymers. The solid line is a scaling formula[16,17]. We use a scaling exponent ν = 0.588[17] and a correction, δ = 1.38, found to give the best fit to the data. (b)-(d) Monomer density as function of the position across the slit for various knot types: 0 1 (×), 3 1 ( ), 6 1 (△), 8 18 ( ), 9 1 (+) and 12 1 (♦). Results, for (b) D = 4σ, (c) D = 12σ and (d) D = 18σ, are plotted as ratios to the values for linear polymers.

Figure 5 .

 5 Figure 5. Distribution of the (instantaneous) radius of gyration in the two free directions, Rpara, for ring polymers in a slit of D = 18σ with different knots: 0 1 (×), 3 1 ( ), 6 1 (△), 8 18 ( ) and 12 1 (♦). Results for 9 1 are omitted to improve the clarity of the plot but are seen to fit with the general trend.

Figure 6 .

 6 Figure 6. Change in the radius of gyration in the two free directions, Rpara, between slits of width D = 18σ and D = 4σ for different knot types. Knot type is indicated on the horizontal axis by C, the essential crossings.

Acknowledgements

We are very pleased to contribute to this Special Issue of Molecular Physics dedicated to Professor Bob Evans in recognition of his important contributions to the field. We wish him much success in the coming years. We would also like to thank Dirk Aarts and Friederike Schmid for helpful discussions.