Andrew J Archer
email: a.j.archer@lboro.ac.uk

Alexandr Malijevský

On the interplay between sedimentation and phase separation phenomena in two-dimensional colloidal fluids

Keywords: colloids, sedimentation, phase transitions, dynamical density functional theory, Brownian dynamics colloids, sedimentation, phase transitions, dynamical density functional theory, Brownian dynamics

Colloidal particles that are confined to an interface effectively form a two-dimensional fluid. We examine the dynamics of such colloids when they are subject to a constant external force, which drives them in a particular direction over the surface. Such a situation occurs, for example, for colloidal particles that have settled to the bottom of their container, when the container is tilted at an angle, so that they 'sediment' to the lower edge of the surface. We focus in particular on the case when there are attractive forces between the colloids which causes them to phase separate into regions of high density and low density and we study the influence of this phase separation on the sedimentation process. We model the colloids as Brownian particles and use both Brownian dynamics computer simulations and dynamical density functional theory (DDFT) to obtain the time evolution of the ensemble average onebody density profiles of the colloids. We consider situations where the external potential varies only in one direction so that the ensemble average density profiles vary only in this direction. We solve the DDFT in one-dimension, by assuming that the density profile only varies in one direction. However, we also solve the DDFT in two-dimensions, allowing the fluid density profile to vary in both the x-and y-directions. We find that in certain situations the twodimensional DDFT is clearly superior to its one-dimensional counterpart when compared with the simulations and we discuss this issue.

Introduction

Whilst colloidal fluids are, of course, of interest because of their importance to industry and every-day-life -for example the blood in our arteries and veins is a complex colloidal suspension -well controlled model colloidal suspension are also studied in order to further our fundamental understanding of fluid behaviour in general. In contrast to atoms, colloids are often large enough that they can be observed under the microscope and so much important knowledge has been gleaned concerning the microscopic properties of liquids in general from studying colloidal suspensions. In fact, it has been found that much from the theory of simple liquids [START_REF] Hansen | Theory of Simple Liquids[END_REF][START_REF] Barrat | Basic Concepts for Simple and Complex Liquids[END_REF] can be utilised for their description. On the left we display a side view of the set-up and on the right we display the view from above. When the surface is tilted at an angle θ to the horizontal, the colloids experience a force m B g sin θ parallel to the surface, which causes them to 'sediment' down to the bottom edge of the surface.

In recent years, the non-equilibrium properties of (colloidal) fluids on microscopic length and time scales has attracted considerable interest [START_REF] Hansen | Theory of Simple Liquids[END_REF][START_REF] Dhont | An Introduction to Dynamics of Colloids[END_REF]. Dynamical density functional theory (DDFT) [START_REF] Marconi | [END_REF][5][6][7] presents a particularly promising approach for such study. DDFT provides a theory for the time evolution of Brownian particles with Langevin stochastic equations of motion and therefore colloidal models serve as the natural subject for the theory. Recent work in this direction includes studies of colloidal sedimentation such as that described in Ref. [8], in which DDFT was found to accurately describe the time evolution of the density profile of sedimenting colloids measured using confocal microscopy, down to the scale of the individual colloids. In this system, the effective interactions between the colloids are purely repulsive and well-modelled by the hard-sphere pair-potential and therefore the sole mechanism for the aggregation of the colloids at the bottom of the sample is the gravity driven sedimentation of the colloids. However, in suspensions where the effective interactions between the colloids includes attraction between the particles, the aggregation can also be driven by these attractions. Such attraction can be generated by the addition of non-adsorbing polymer to the solution, which results in attraction due to the so-called depletion mechanism [9,10], and can lead to the suspension phase separating into a colloid rich 'liquid' phase and a colloid poor 'gas' phase. Such a system has, for example, been used to study how droplets of the colloidal liquid phase coalesce with the bulk of this phase that has already sedimented to the bottom of the container [11]. The benefit of studying this phenomenon in a colloidal fluid is that it occurs on length and time scales accessible to optical microscopy. The authors of Ref. [11] were able to study in detail the various stages of the drop coalescing with the bulk fluid.

There is also a great deal of interest in two-dimensional (2D) colloidal fluids. Such systems occur when the colloids become confined to an interface. This can be at a fluid interface, such as the air-water interface, where the colloids often become pinned to the interface by capillary forces [12,13]. Alternatively, a 2D system may be obtained when the colloids settle under gravity to the bottom of their container. Once on the bottom, the colloids are still mobile and able to diffuse over the surface, but because the height above the interface h th that they are able to reach due to the thermal fluctuations is much less than the colloid diameter σ at common temperatures, they can effectively be considered to be a 2D fluid of particles constrained to move on the bottom surface of the container. The height h th ∼ k B T /m B g, where k B is Boltzmann's constant, T is the temperature, m B is the buoyant mass of the colloids and g is the acceleration due to gravity. In Refs. [14,15] the authors studied the flow behaviour of such a 2D fluid down a narrow channel with a width of only a few times the diameter of the particles. In this set up, the colloids are ordered into layers within the channel and a reduction in the number of layers due to a density gradient along the channel was observed. In this paper, we present results obtained for a model 2D system that is related to the system studied in Refs. [14,15]. In Fig. 1 we display a sketch of our system. We consider cases where the colloids are initially uniformly dispersed over the surface, with a fairly low average surface coverage. However, due to the fact that there is a uniform external force on the colloids towards one side of the surface, they begin to 'sediment'. In addition to this, there are attractive forces between the colloids in our system. As the colloids move over the surface, both due to the driving and also due to their Brownian motion, they come into contact with one another and due to the attractive interactions between the particles, they aggregate together. We study the interplay between these two effects on the time evolution of the density profile of the colloids using both Brownian Dynamics (BD) computer simulations and DDFT. We solve the DDFT in both one-dimension (i.e. assuming that the density profiles only varies in the y-direction, which is the only direction in which the external potential varies) and we also solve the DDFT in 2D, allowing the fluid density profiles to vary in both the x-and y-directions. Perhaps the most striking result in this paper is that we find that in certain situations the results from solving the DDFT in 2D and then subsequently averaging over the x-direction agree much better than the 1D DDFT results with what we find from the BD simulations.

This paper is structured out as follows: In Sec. 2 we describe our model system and give the governing stochastic equations of motion for the particles. In Sec. 3 we present BD simulation results showing the typical behaviour of the system. Then, in Sec. 4 we present our theory for the system and give a brief discussion of the equilibrium properties predicted by the theory. In Sec. 5 we present our DDFT results and make comparison with our BD simulation results. Finally, in Sec. 6 we make a number of concluding remarks.

Model System

We model the colloids by assuming that they interact via the following pair potential:

v(r) = v hd (r) + v at (r), ( 1
)
where r is the distance between the centres of the particles. The hard-disk pair potential is given by

v hd (r) = ∞ r ≤ σ 0 r > σ, (2)
and the attractive contribution to the pair-potential is modelled by the following simple form:

v at (r) = -ǫ exp(-r/l), ( 3
)
where ǫ is a parameter that determines the strength of the attraction and l sets the range of the potential. For simplicity we set l = σ, the diameter of the hard-core of the particles. Note that at contact v(r = σ+) = -ǫ/e. We model the external potential acting on the colloids as follows:

u ext (x, y) =      ∞ y ≤ 0 ay 0 < y < L, ∞ y ≥ L (4)
where y is the coordinate direction perpendicular to the hard parallel walls that are located at y = 0 and y = L along the edge of the surface on which the particles are located. The y-axis is also parallel to the uniform external potential ay, where a = m B g sin θ, for the set-up pictured in Fig. 1. We consider a system of N particles and we denote the set of position coordinates of the particles r N = {r i ; i = 1, . . . , N }. The potential energy of the system is given by

U N (r N , t) = N i=1 u ext (r i , t) + 1 2 N i=1 j =i v(|r i -r j |). (5)
We assume that the dynamics of the colloids, with positions r i (t), is governed by the following set of over-damped stochastic equations of motion [START_REF] Dhont | An Introduction to Dynamics of Colloids[END_REF]:

Γ -1 dr i (t) dt = -∇ i U N (r N , t) + ξ i (t), (6)
where Γ -1 is a friction coefficient characterising the one-body drag of the solvent on the particles and ξ i (t) is a stochastic white noise term [START_REF] Dhont | An Introduction to Dynamics of Colloids[END_REF][START_REF] Marconi | [END_REF][5][6][7]. Note that by modelling the equations of motion via Eq. ( 6) the hydrodynamic interactions between the colloids are neglected [START_REF] Dhont | An Introduction to Dynamics of Colloids[END_REF].

Brownian Dynamics Simulation Results

We have performed standard Brownian Dynamics computer simulations [START_REF] Allen | Computer simulation of liquids[END_REF] of the time evolution of the particles in our system. This corresponds to numerically integrating Eqs. (6). To make this integration more straight-forward, we replace the discontinuous hard-disk potential (2) with the following continuous potential: Choosing the amplitude βα = 20, where β = 1/k B T is the inverse temperature, together with a large value for the exponent n = 50 in v * hd (r), ensures that using this potential instead of the hard-disk-potential (2) leads to no change in the structure of the fluid. We initiate our system with the non-overlapping particles randomly distributed over the surface, which corresponds to an equilibrium configuration for a system where there is no attraction between the particles (ǫ = 0) and the drive force amplitude a = 0, corresponding to the particles being on a flat horizontal surface.

v * hd (r) = α exp(-(r/σ) n ). (7)
In Fig. 2 we display snapshots from a simulation for the case when L = 43.5σ, βa = 0.1, βǫ = 4.8 and there are N = 200 particles in the system, corresponding to an average density ρσ 2 = 0.072. Note that when βǫ = 4.8, this corresponds to a potential well depth βv(r = σ+) = -βǫ/e = -1.77. This amounts to a relatively strong attraction between the particles and over time it leads to the particles aggregating into clusters. We also see in Fig. 2 the influence of the external potential which causes the particles to sediment to the bottom. Due to the fact that there is no attraction between the particles and the wall, the 'drops' of the colloidal liquid phase do not spread over the lower wall, so that it remains 'dry'. From the snapshots in Fig. 2 we infer that in this situation the contact angle is close to 180 o . Note however, that the contact angle decreases as the strength of the driving force a due to the external potential is increased. When a is increased to the value βa = 1, then as can be observed from the BD simulation results displayed in Fig. 3, we find that the force on the particles towards the lower edge is strong enough to cause them to spread to form a uniform layer there. Thus, increasing the parameter a, which corresponds to increasing the angle θ of the surface to the horizontal (c.f. Fig. 1), enhances the coalescence effect.

In Fig. 4 we display the fluid one-body density profiles ρ(y), which are calculated by averaging the full 2D density profile ρ(x, y) over the x-direction. These are obtained for the same parameters as for the results displayed in Fig. 2 by averaging over 40 different simulation runs, each with a different initial configuration and realisation of the random noise terms in Eq. ( 6). We see that the contact density ρ(y = 0) at the lower wall, is much less than the density just above the wall, at around y ∼ 3σ. This low contact density comes as a consequence of the fact that the amplitude a of the driving force is small and also that there are no attractive forces between the particles and the wall and so the particles seek to gather to themselves away from the wall. The surface tension (excess free energy) between the wall and the liquid is therefore very similar in value to the liquid-gas surface tension between the colloidal liquid and colloidal gas phases. As a consequence of this and the fact that the surface tension between the vapour and the wall is low, the drops on the surface displayed in Fig. 2 have a high contact angle, as one would expect from Young's equation. DFT shows that for a given external potential u ext (r), there is a unique equilibrium one-body density profile ρ(r) which is obtained as the minimum of the following grand potential functional [START_REF] Hansen | Theory of Simple Liquids[END_REF][START_REF] Evans | [END_REF][START_REF] Evans | Fundamentals of Inhomogeneous Fluids[END_REF]:

Ω[ρ(r)] = F [ρ(r)] -µ drρ(r), ( 8
)
where µ is the chemical potential and F [ρ(r)] is the Helmholtz free energy functional. Minimization leads to the following Euler-Lagrange equation that may be solved for the equilibrium fluid density profile:

δF [ρ(r)] δρ(r) = µ. (9)
The Helmholtz free energy functional for a 2D fluid is:

F [ρ(r)] = k B T drρ(r)[ln(Λ 2 ρ(r)) -1] +F ex [ρ(r)] + dr u ext (r)ρ(r). ( 10
)
The first term on the right hand side is the ideal-gas Helmholtz free energy; Λ is the thermal de Broglie wavelength. F ex [ρ(r)] is the excess Helmholtz free energy, which is the contribution to the free energy from the particle-particle interactions. This quantity is not known exactly, however, there are various well-trodden paths available in the literature for obtaining a good approximation for this quantity [START_REF] Hansen | Theory of Simple Liquids[END_REF][START_REF] Evans | [END_REF][START_REF] Evans | Fundamentals of Inhomogeneous Fluids[END_REF]. The final term in Eq. ( 10) is the contribution to the free energy due to the external potential.

For the work described in this paper, we approximate the Helmholtz free energy functional as follows:

F [ρ(r)] = drf hd (ρ(r)) + 1 2 dr dr ′ ρ(r)ρ(r ′ )v at (|r -r ′ |) + drρ(r)u ext (r), (11)
where f hd (ρ) is the Helmholtz free energy per unit area of a uniform fluid of harddisks with bulk density ρ. We use the scaled particle approximation [START_REF] Rosenfeld | [END_REF],

βf hd (ρ) = ln(Λ 2 ρ) -2 -ln(1 -η) + (1 -η) -1 , (12)
which includes the exact ideal-gas contribution to the free energy and where η = πρσ 2 /4 is the packing fraction. The second term in Eq. ( 11) is a simple meanfield approximation for the contribution to the free energy from the attractive interactions between the particles. Note that the approximation for F [ρ] in Eq. ( 11) does not give reliable results for the fluid density profile when ǫ = 0, or when the external potential due to the wall includes an attractive tail. In this case there are oscillations in the density profile near the wall, due to the effects from packing of the hard cores of the particles against the wall. This is an effect which . Equilibrium (t → ∞) density profiles obtained using DFT for the case when the parameters in the external potential are L = 43.5σ, βa = 0.1 and when the average density in the system is ρσ 2 = 0.2. We display the (exact) ideal-gas density profile, together with the density profiles obtained using Eq. ( 11) when there is no attraction between the particles βǫ = 0 (pure hard-disks) and when the attraction amplitude is βǫ = 2.5.

Eq. ( 11) is completely unable to describe. However, for the present model fluid, the approximation in Eq. ( 11) is good, as long as a is small (βa ∼ 0.1), and ǫ is sufficiently large (βǫ 1), owing to the fact that as ǫ is increased, the density of the particles in contact with the wall is low, as discussed above in the context of the results displayed in Figs. 2 and4.

To further illustrate this point, in Fig. 5 we display the particle density profiles obtained using Eq. ( 11) for the case when the parameters in the external potential are L = 43.5σ and βa = 0.11 . In our calculations we have set the value of the chemical potential µ so that in all cases the average density in the system is ρσ 2 = 1 L L 0 dyρ(y)σ 2 = 0.2. We display in Fig. 5 the density profile for the case when βǫ = 2.5, which corresponds to a strong enough attraction between the particles for the system to exhibit liquid-gas phase separation and due to the external driving force a towards the bottom, we find the high density liquid phase is pushed to the bottom of the system. However, the density near to contact with the wall (at small values of y) is low, which justifies our use of the local density approximation (LDA) for the hard-disk contribution to the free energy in Eq. (11). We also display in Fig. 5 the (exact) ideal-gas density profile and the density profile obtained using Eq. ( 11) when there is no attraction between the particles βǫ = 0 (i.e. pure hard-disks). We see that in this case the density profile near the wall is monotonic, whereas in reality we should expect to see oscillations due to packing near the wall [START_REF] Evans | Fundamentals of Inhomogeneous Fluids[END_REF].

From Fig. 5 we can also observe the influence of the nature of the particle interactions on the fluid density profiles: on comparing the ideal gas density profile with the ǫ = 0 pure hard-disk case, we see that the effect of the excluded area (hardcore) interactions between the particles is to push some of them further up into the system. However, the effect of adding attraction between the particles is to reverse this effect, so that when βǫ = 2.5 the particles are on average much closer to the bottom of the system. This effect, of course, is also due to the external potential (a > 0) driving the particles to the bottom of the system. (11). Within the binodal curve the uniform fluid is predicted to be unstable. Note that our simple theory is not able to predict the existence of the hexatic or solid phases.

On the nature of the density profiles

There is an important issue that we must raise at this point in our discussion: When a fluid is confined in an external potential such as that given in Eq. ( 4), that only varies in one direction (in this case the y-direction), then we know that the fluid density profile ρ(r) will only vary in this direction. This is because ρ(r) is an ensemble average quantity, which means that even if particular configurations of the system are inhomogeneous in the x-direction, such as in the configuration displayed in the final panel of Fig. 2 (which is an equilibrium configuration), the average over the full ensemble of possible configurations will not vary in the xdirection. Thus, in calculating the density profiles displayed in Fig. 5, we have used this fact to reduce solving Eq. ( 9), to that of solving for a one-dimensional (1D) function ρ(y). However, although formally ρ(r) does not vary in the x-direction, as we show below, when one uses an approximate free energy functional, such as that in Eq. ( 11), one can in fact obtain density profiles that vary in both the x-and in the y-directions. We will further discuss this issue below.

Bulk phase diagram

From our DFT (11) we are able to calculate the bulk (uniform) fluid phase diagram for our system, by setting the density in Eq. ( 11) to be a constant ρ(r) = ρ, which gives the bulk fluid Helmholtz free energy. From this the phase diagram may be calculated in the usual way [START_REF] Chandler | Introduction to Modern Statistical Mechanics[END_REF]. The resulting phase diagram is displayed in Fig. 6. We see that as the temperature T is decreased or, equivalently, as ǫ is increased, the theory predicts that the system exhibits phase coexistence between a low density suspension and a high density suspension of particles (i.e. gas-liquid phase coexistence). The pairs of coexisting state points define the locus of the binodal. The line in the phase diagram along which the compressibility is predicted to be zero defines the spinodal. Inside this curve the DDFT (see Eq. ( 15) below) predicts that the uniform fluid is linearly unstable. Note that our simple theory is not able to predict the transition to either the hexatic phase nor the solid phase, that we expect to find at higher densities. Note also that the phase diagram in Fig. 6 is only qualitatively correct in describing the gas-liquid phase separation. For example, when βǫ = 4.8, corresponding to the results displayed in Figs. 2 and4, we find that the density of the liquid at coexistence is ρσ 2 ≃ 0.85. However, from Fig. 6 we see that we must select βǫ = 3.5 (i.e. k B T /ǫ = 0.29) in order for the density of the liquid at coexistence to equal this value. As a consequence of this, we find below in Sec. 5 better agreement between the DDFT results for βǫ = 3.5 with the BD simulation results for βǫ = 4.8, than we do for the simulation results for βǫ = 3.5.

DDFT

DDFT [START_REF] Marconi | [END_REF][5][6][7] was developed to describe the dynamics of the one-body density of a fluid of Brownian particles with equations of motion given by Eq. ( 6). The Smoluchowski equation [START_REF] Dhont | An Introduction to Dynamics of Colloids[END_REF]6]:

∂P (N ) ∂t = Γ N i=1 ∇ i • [k B T ∇ i P (N ) + ∇ i U N P (N ) ] (13)
governs the time evolution of the N -particle probability density P (N ) (r N , t) for systems with dynamics governed by Eq. ( 6). The one-body density is obtained by integrating over this function, as follows:

ρ(r 1 , t) = N dr 2 ... dr N P (r N , t). (14)
This non-equilibrium density profile represents an ensemble average over all realisations of the stochastic noise [7]. On integrating the Smoluchowski equation ( 13) we obtain the central equation of DDFT [START_REF] Marconi | [END_REF][5][6][7]:

∂ρ(r, t) ∂t = Γ∇ • ρ(r, t)∇ δF [ρ(r, t)] δρ(r, t) . (15)
Note that this depends on the equilibrium Helmholtz free energy functional F [ρ], given in Eq. (10). In order to derive Eq. ( 15), we have made the approximation that equal-time two-body correlations at time t in the non-equilibrium fluid are the same as those of an equilibrium fluid with the same one-body density profile ρ(r, t) [START_REF] Marconi | [END_REF][5][6]. As we see below, this closure approximation is good for the situation of interest here. Other situations where DDFT (15) has proved to be reliable are discussed in Refs. [START_REF] Marconi | [END_REF]8,[START_REF] Dzubiella | [END_REF][22][23][24][25][26][27][28]. We should emphasise that in addition to the closure approximation used in deriving the DDFT (15), there is the additional approximation that one must make in selecting a suitable approximation for the Helmholtz free energy functional F [ρ(r)]. The approximation that we use here is that in Eq. (11). We solve the DDFT by discretizing the equations on a cartesian grid, with grid spacing ∆x = ∆y = 0.5σ, and we set the initial time t = 0 density profile to be that of the ideal-gas in the equilibrium situation when the drive in the external potential a = 0 (i.e. ρ(x, y) = ρ, for 0 < y < L), which has added to it a small random number at each grid point. This random number is drawn from a uniform distribution between ±ρ/20, where ρ is the average fluid density in the system. The results do not depend on the amplitude of this noise. The random noise is added to the initial time t = 0 density profile in order to break the symmetry of the profile which subsequently in some situations allows the system to develop density profiles that vary in both the x-and in the y-directions, when the DDFT is solved in 2D. Note that the density profiles from the DDFT do not vary in the x-direction in all situations; they only vary in this direction in certain specific situations -see below.

One can also assume that the density profiles only vary in the y-direction (which is strictly true, since the external potential only varies in the y-direction) and use this fact to reduce solving Eq. ( 15), to that of solving for a one-dimensional (1D) time series of density profiles ρ(y, t). Below, we compare results from doing this, to the results obtained from the full 2D solution to the DDFT.

Results

All the results displayed in this section are for the case when the parameters in the external potential (4) take the values L = 43.5σ and βa = 0.1. In Fig. 7 we display results for the case when there are 200 particles in the system, which corresponds to an average density ρσ 2 = 0.072, and when βǫ = 1. This level of attraction between the particles is not strong enough for the system to exhibit liquid-gas phase separation as can be inferred from Fig. 6. The system takes a time t * ∼ 400 to reach equilibrium, where t * ≡ t/τ B , and where τ B = βσ 2 /Γ is the Brownian time scale, which is the time it takes for a particle to diffuse a distance of order its own diameter. Since the amplitude of the driving force a is rather small, we see that the density near the top of the system ρσ 2 ∼ 0.01 is not insignificant. This also comes as a result of the attraction between the particles being quite weak. We see that the agreement between the BD simulation results and the DDFT results is fairly good. Note that in this situation, the DDFT produces exactly the same result whether it is solved in 1D or in 2D. This is in contrast to the results which we display below for cases where βǫ is sufficiently large for phase separation to be observed.

In Fig. 8 we display density profiles obtained from the 2D DDFT for the case when βǫ = 3.5 and ρσ 2 = 0.072 at various times during the sedimentation process. These results should be compared with the BD simulation results (which are for βǫ = 4.8) displayed in Figs. 2 and4. We see that the theory predicts that the system prefers to form a number of droplets, rather than form a layer of uniform thickness along the bottom. Strictly speaking, since ρ(r, t) is an ensemble average quantity and since the external potential only varies in the y-direction, these profiles should also only vary in the y-direction. It is rather striking, however, that the resulting density profiles in-fact resemble more closely what we observe in the simulation snapshots in Fig. 2.

We also calculate the density profiles in this situation using the 1D DDFT. We display these results in Fig. 9, together with the results from the 2D DDFT that have subsequently been averaged over the x-direction to produce a 1D density profile. Initially, for t * 50 there is no difference between the results from these two approaches. However, for latter times, when the 2D DDFT predicts that the system breaks up into a number of separate drops, the results from the two approaches are, of course, different. We also display in Fig. 9 the BD simulation results (c.f. Figs. 2 and4), and remarkably the results from the 2D DDFT that are subsequently averaged over the x-direction are in much better agreement with the BD simulation results then the 1D DDFT results. This suggests that the additional degree of freedom in the 2D DDFT helps to improve the theory, by 'allowing' it to better describe the density inhomogeneities that are present in the x-direction.

In Fig. 10 we display density profiles obtained from the 2D DDFT for the case when βǫ = 3.5 and ρσ 2 = 0.2, at various times during the sedimentation process. The average density in this case is higher than in the case considered above in Fig. 8. We see that the initially uniform fluid breaks up into droplets of the liquid phase, due to the attraction between the particles. Due to the external potential, these droplets sediment to the bottom of the system and coalesce to form a single uniform layer along the bottom of the system. This is in contrast to the case in Fig. 8, and is due to the fact that in the case in Fig. 10 there are more particles in the system. This difference is best understood by considering the fact that as the density of the fluid is increased, the area of the surface that is covered by the two-dimensional liquid, which we denote A, also increases. As discussed above in Sec. 3, the interfacial tension between the liquid and the surface is roughly equal to the interfacial tension γ between the liquid and the gas. Thus, as the number of particles (average density) in the system is increased, the ratio l/A decreases, where l is the length of the interface between the liquid and the gas plus the length of the interface between the liquid and the wall (i.e. it is the total length of the interface(s) around the 2D liquid). As the system seeks to minimise the free energy, it therefore seeks to minimise the contribution due to these interfaces ∼ γl. Thus, the final equilibrium configuration is essentially determined by minimising l subject to the constraint that A is fixed. When the ratio l/A is large (low ρ), the system must form drops to minimise l; this is the case in Fig. 8. However, when l/A is small (i.e. larger ρ), l is minimised by the system forming a layer of uniform thickness along the bottom; this is the case in Fig. 10. Note, however, that this argument only applies in the limit when a is small. As the amplitude of the external potential a is increased, then the external potential energy contribution to the free energy overcomes the interfacial tension contribution lγ and the particles are squeezed to the bottom of the system, as can be seen in Fig. 3. This argument should also apply to three-dimensional (3D) systems in a similar situation, as long as one instead considers the ratio A/V , where A is the surface area of the liquid and V is its volume. In this case, the system will seek to minimise the free energy by balancing the interfacial contribution ∼ γA, against the external potential energy contribution. We have shown that the behaviour during sedimentation of a 2D model colloidal suspension with attraction between the particles is rather rich due to the interplay between the aggregation stemming from the particle attractions and the sedimentation due to the external potential. We have used DDFT to describe the dynamics of the system, and as the results in Figs. 7 and9 show, the theory is in quite good agreement with the simulation results. The agreement is particularly good when the colloids are undergoing pure sedimentation (Fig. 7). However, even for the case displayed in Fig. 9, where in addition to the sedimentation, the particles are also aggregating due to the strong attractions between them, we still observe good qualitative agreement between the DDFT and the simulations.

The results from our work raise a key question: why does the DDFT when it is solved in 2D, and then afterwards averaged over the x-direction, give better agreement with the simulations than when the DDFT is solved in 1D? As discussed above in Sec. 4.2, solving equilibrium DFT is in principle equivalent to evaluating the partition function for the system, and therefore since our external potential varies in only the y-direction, the density profile obtained by averaging over the ensemble of all configurations of the system, must also only vary in the y-direction. This should at least be the case for the final equilibrium density profile obtained from the DDFT, since recall that the equilibrium t → ∞ solution of the DDFT is the solution of Eq. ( 9) -i.e. it is a minimum of the free energy [6] and is equivalent to solving equilibrium DFT. However, as can be seen in the final frame for t * = 300 in Fig. 9, this is not the case. The final frame of Fig. 9 shows that the results from the 2D calculation are in very good agreement with the simulation results, whereas the results from the 1D calculation do not agree with the simulation results.

When our DDFT (DFT) is solved in 2D, we believe the reason that it does not always give a density profile that only varies in the y-direction, is due to the fact that we are not using the (unknown) exact free energy functional, but instead, we are using an approximation for this quantity -specifically that given in Eq. ( 11). This approximate free energy functional neglects some of the fluctuation contributions to the free energy. This manifestation has been known for many years now, in particular, in the context of using DFT to calculate the density profile for the free interface between a liquid and its vapour. It is known that DFT fails to include some of the interfacial fluctuations [START_REF] Evans | [END_REF]. See also the more recent work in Refs. [29][30][31][32][33] and references therein. This issue has also been discussed in a very illuminating manner by Reguera and Reiss [34]. In Ref. [34] the behaviour of drops of a 3D liquid, with a fixed number of particles and confined within a spherical cavity were studied using both DFT and simulation. In terms of the physics involved, the situation studied in Ref. [34] is somewhat akin to the situation we study here, but with the external potential parameter a in our model set to zero. Using an approximation for the Helmholtz free energy functional that is the 3D analogue of the functional in Eq. ( 11), the authors of Ref. [34] showed that the DFT is unable to describe the fluctuations corresponding to drop translations within the cavity and also fluctuations corresponding to variations in the number of particles in a liquid drop. By comparing their DFT results with Monte-Carlo simulation results for a system where the centre of mass is constrained to coincide with the centre of the cavity (i.e. so that the translational fluctuations are suppressed) they found better agreement between the DFT and the simulation results. In other words, the DFT describes well the density distribution within the drop(s) but it is unable to describe the translational fluctuations of the drop(s). This is also the conclusion that we draw on the basis of the results in Fig. 9 and leads us to surmise that when DFT is implemented together with an approximate free energy functional, such as that in Eq. (11), it leads to one obtaining an average over a subsection of phase space, rather than averaging over all of phase space.

In Ref. [34] the authors also discuss the relationship between DFT and field theory (FT). In FT, the central quantity is a Hamiltonian E[φ(r)], which is a functional of the order-parameter φ(r), which may be interpreted as being a coarsegrained density profile (See Ref. [35] for a possible definition of this quantity). When solving the FT to obtain the order parameter profile φ(r) which corresponds to the most probable configuration of the system, one is required to solve the following Euler-Lagrange equation [34]:

δE[φ(r)] δφ(r) = µ, (16)
[c.f. Eq. ( 9)] where µ is the Lagrange-multiplier that originates from the constraint from fixing the total number of particles in the system, which mathematically plays the same role as the chemical potential. Now, of course, the functional E[φ(r)] is not known exactly and so one is required to make an approximation for this quantity based on physical insight and any other knowledge we have about the system of interest. In DFT, there is a similar process used to construct an approximation for the free energy functional F [ρ] [START_REF] Hansen | Theory of Simple Liquids[END_REF][START_REF] Evans | [END_REF][START_REF] Evans | Fundamentals of Inhomogeneous Fluids[END_REF]. The end result of these two processes can be very similar or even the same; i.e. one considers ρ(r) and φ(r) to be the same quantity (strictly speaking, they are not the same, and ρ(r) = φ(r) , i.e. ρ(r) is a statistical average over all possible configurations of the coarse-grained density φ(r) [7,35]) and therefore one is tempted to make the same approximations in constructing E[φ] as are made in constructing F [ρ]. See Ref. [35] for an example of a FT constructed along these lines. With these arguments in mind, it is now clear why our 2D DDFT results agree better with the simulation results than the 1D DDFT results: Our approximation for the free energy ( 11) is unable to describe the translational fluctuations of drops in the system. However, by solving in 2D we obtain density profiles which correspond to particular likely configurations of the system. By subsequently averaging this density profile over the x-direction, we are averaging over several droplet profiles, and so we get better agreement with the simulations than the 1D DDFT because we are averaging over more 'typical' configurations. In this respect, we are treating the DFT more like a FT. If we fully take the FT point of view, then the natural extension of our approach is to perform several different DDFT calculations for each situation, with each calculation having a different realisation of the random noise that we add to the initial density profile (see Sec. 4.4). One would then average over all of these different density profiles and presumably therefore average over more of phase space (i.e. averaging over more possible fluid configurations). We believe that doing this may give even better agreement with the simulation results than our present 2D DDFT results. We should add, however, that following this approach would be rather computationally expensive, particularly if one were investigating a 3D system -it may be quicker to just simulate the system.

We believe that if we used a more accurate DFT such as Rosenfeld's fundamental measure theory [START_REF] Rosenfeld | [END_REF] in our 2D DDFT calculations, we would still see the same 'symmetry breaking' (i.e. the theory would predict that the density profiles have variations in the x-direction). Using this alternative approximation for the free energy functional would undoubtably improve the description of the fluid structure by doing a better job of describing the short range correlations between the particles. However, what is really required for the problem here is a DFT that includes and averages over the large scale density fluctuations; in particular those corresponding to modes such as a droplet of particles translating along the bottom of the system. Constructing such a DFT requires developing a theory which is capable of describing fluctuations on several different length scales, which is a notoriously difficult problem in physics.

A final point that we should emphasise is that in both our simulations and in our theory we have not included the hydrodynamic interactions that play an important role in the dynamics of real colloids [START_REF] Dhont | An Introduction to Dynamics of Colloids[END_REF]. It is possible to include these interactions between the particles within the DDFT approach [36,37]. However, we leave this as future work.

Figure 1 .

 1 Figure1. The model system that we consider, consisting of colloidal particles distributed over a surface. On the left we display a side view of the set-up and on the right we display the view from above. When the surface is tilted at an angle θ to the horizontal, the colloids experience a force m B g sin θ parallel to the surface, which causes them to 'sediment' down to the bottom edge of the surface.

Figure 2 .

 2 Figure 2. Snapshots from our BD simulations for the case when L = 43.5σ, βa = 0.1, βǫ = 4.8 and there are N = 200 particles in the system, corresponding to an average density ρσ 2 = 0.072. Results are displayed at the times t/τ B ≡ t * = 0, 50, 150, and 400, where τ B = βσ 2 /Γ is the Brownian time-scale.

Figure 3 .Figure 4 . 2 .

 342 Figure 3. Snapshots from our BD simulations for the case when all the parameters (including the initial time t = 0 starting positions of the particles) are the same as in Fig.2. the only difference is that the driving force is one order of magnitude larger -i.e. βa = 1. Results are displayed at the times t/τ B ≡ t * = 10, 20, 30, and 70, where τ B = βσ 2 /Γ is the Brownian time-scale.

Figure 5

 5 Figure 5. Equilibrium (t → ∞) density profiles obtained using DFT for the case when the parameters in the external potential are L = 43.5σ, βa = 0.1 and when the average density in the system is ρσ 2 = 0.2. We display the (exact) ideal-gas density profile, together with the density profiles obtained using Eq. (11) when there is no attraction between the particles βǫ = 0 (pure hard-disks) and when the attraction amplitude is βǫ = 2.5.

Figure 6 .

 6 Figure 6. Phase diagram in the reduced temperature k B T /ǫ versus density ρ plane, obtained from our simple approximation for the free energy in Eq.(11). Within the binodal curve the uniform fluid is predicted to be unstable. Note that our simple theory is not able to predict the existence of the hexatic or solid phases.

Figure 7 .Figure 8 .Figure 9 .

 789 Figure 7. Time evolution of the fluid density profiles ρ(y, t). The solid lines are obtained from BD simulations, for L = 43.5σ, βa = 0.1, βǫ = 1 and there are N = 200 particles in the system, corresponding to an average density ρσ 2 = 0.072. We also display results from DDFT (dashed line). Results are displayed at the times t/τ B ≡ t * = 0, 100 and 400, where τ B = βσ 2 /Γ is the Brownian time-scale.

Page 11 of 16 URLFigure 10 .

 1610 Figure 10. Fluid density profiles obtained from DDFT for the case when βa = 0.1, βǫ = 3.5 and ρσ 2 = 0.2 at the times t * = 1, 20, 40, 60, 100, 160, 200, 260 and ∞ (starting top left and finishing bottom right).

Note that in our DFT and DDFT calculations we replace the upper hard wall at y = L in Eq. (4) by a softened purely repulsive form. This is done to improve the stability of our numerical algorithm for solving the DDFT. This, of course, has no effect on the behaviour of the system at the lower wall, which is the primary focus of our work described here.

Theory for the system

We use dynamical density functional theory [4][5][6][7], in order to develop a theory to further understand the behaviour of our system. Since DDFT is based on equilibrium density functional theory (DFT) [1,17,18], which is a theory for the equilibrium one-body density profile ρ(r) of an inhomogeneous fluid, we start this section by introducing DFT and the approximation we use for the free energy functional. Once set, DFT also allows one to calculate thermodynamic quantities and phase behaviour for the system.

Acknowledgements

It is a pleasure to dedicate this paper to Professor Bob Evans, whose contributions to liquid state theory are immense. He has inspired several generations of scientists in the field and it has been a privilege and a joy to have worked with him. We wish him all the best in the years to come.

AJA gratefully acknowledges financial support from RCUK and AM acknowledges financial support of the MSMT of the Czech Republic under Project No. LC512 and the GAAS of the Czech Republic (Grant No. IAA400720710).


Roland Roth

Andrew O Parry

Drying in a Capped Capillary

Keywords: wetting and drying, capillary condensation, density functional theory

We use a model Density Functional, based on the White-Bear version of Fundamental Measure theory, to test recent predictions, due to Evans and co-workers, that capillary condensation in a capped capillary-slit is a continuous interfacial critical phenomenal related to the complete wetting transition. Using a model with a square-well intermolecular fluid-fluid attraction we first determine accurately the location of the first-order capillary evaporation transition in an infinite (open) hard-wall capillary slit. Extending the density functional model to allow for a two dimensional order-parameter profile, we then study the adsorption as the chemical potential is reduced to capillary evaporation but now in a capillary-slit that is capped at one end. The equilibrium density profiles obtained show that, sufficiently close to the phase boundary, a meniscus separating liquid-like and vapour-like phases forms near the capped end, and that as capillary evaporation is approached, continuously unbinds from the capped end. Our numerical results indicate that the divergence of the adsorption due to the unbinding of the meniscus is logarithmic and is the same as for the complete wetting transition in systems with short-ranged forces.

Introduction

In 1979 Bob Evans wrote a pioneering review article [1] on the application of classical density functional theory (DFT) to the statistical mechanics of inhomogeneous fluids and, with exceptional clarity, explained how free-energies, density profiles and correlation functions can all be calculated from a model functional. It is no exaggeration to say that this brought "calculation to the masses" paving the way for the development of ever more sophisticated model microscopic functionals [2][3][4][5][6][7] together with their application to phase transitions such as wetting, layering and capillary condensation [START_REF] Dietrich | Phase Transitions and Critical Phenomena[END_REF][START_REF] Evans | [END_REF][10][11][12]. Together with effective interfacial Hamiltonian models and computer simulations, DFT is part of the triumvirate of powerful tools theoreticians have used to revolutionise our understanding of interfacial phenomena, which has grown enormously since the 1980s. Interest in the last decade has been further renewed by technological developments which allow one to pattern and sculpt solid surfaces on the nanometre and micrometre scale. Thus as well as considering adsorption at flat substrates and in capillary slits and pores one can consider corrugated surfaces and wedges/grooves of different cross-section [13][14][15][16][17], and also heterogeneous surfaces patterned into domains with different wetting properties. As well as being of huge practical importance this is of fundamental interest to statistical mechanics since different substrate geometry can induce novel types of interfacial phase transition which lie "between" wetting and capillarycondensation. For example there has been a great deal of interest in the physics of filling transitions for fluids adsorbed in linear wedges which can exhibit very large interfacial fluctuation effects compared to wetting at flat substrates [18,19]. It has also been noticed, from studies based primarily on effective interfacial Hamilto- nian models, that there appear to be precise connections between absorptions in different geometries and between apparently different phase transitions [20,21]. The purpose of the present article is to test some of these predictions using more microscopic density functional models. By doing so we follow in spirit an approach pioneered by Bob Evans. Capillary-condensation refers to the shift of the liquid-vapour coexistence curve when a fluid is confined in a capillary slit or pore. Thus vapour confined in a slit of width L, as shown in Fig. 1(a), will condense to a liquid-like phase, at a pressure p C (L, T ), which is different to the value p sat (T ) at bulk coexistence, given by the macroscopic Kelvin equation

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¢ ¡¢
¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ ¦ ¡¦ § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § L L (a) (b)
p sat (T ) -p C (T, L) = 2σ cos θ L + • • • (1)
where σ is the liquid-vapour surface tension and θ is the contact angle. Thus p C (L, T ) < p sat (T ) for contacts angles 0 ≤ θ < π/2 corresponding to capillary condensation, while p C (L, T ) > p sat (T ) for contact angles π/2 < θ ≤ π corresponding to capillary evaporation. The vast majority of studies of capillary condensation based on Landau-like square gradient models and density-functional theories, assume that the walls are of infinite area and the density is translationally invariant parallel to them. Thus the equilibrium one-body density profile ρ(r) = ρ(x) where x is the co-ordinate perpendicular to the wall(s). In this case capillary condensation is certainly a first-order phase transitions and mean-field treatments of it display adsorption isotherms with a van der Waals loop. However capping the capillary at one end [see Fig. 1(b)] may have a profound effect on this transition. Building on earlier numerical studies, Parry, Rascon, Wilding and Evans [22], pointed out that the presence of the capped end necessarily induces the formation of a meniscus, separating liquid-like and vapour-like phases. For walls that are completely wet by liquid interfacial Hamiltonian models indicate that the meniscus is bound near to the capped end for pressures p < p C (T, L). As the pressure is increased to p C (T, L) the meniscus continuously unbinds from the capped end thus continuously filling the capillary with liquid and is directly analogous to the complete wetting transition occurring at flat walls. However instead of the unbinding of a liquid-vapour interface as bulk coexistence is approached, the unbinding is that of a meniscus as the shifted phase boundary (corresponding to capillary condensation) is approached. Thus condensation in a capped capillary is actually a second-order interfacial phase transition. For systems with short-ranged forces, and in the ab- We wish to test this prediction with a microscopic DFT model. This is a nontrivial numerical task of course because the equilibrium density profile ρ(r) = ρ(x, z) is no longer translationally invariant (parallel to the walls) and depends explicitly on z the distance from the capped end. We choose a model fluid with a square-well intermolecular potential and a capillary-slit with purely hard walls. These are completely wet by the vapour phase so our attention is on the capillary evaporation which occurs as the pressure is reduced to p C (T, L). In this case the meniscus separates vapour (towards the capped end) from liquid but otherwise the phenomenology is identical. We begin with a brief recap of the density functional method before accurately determining the location of the capillary-evaporation pressure in an infinite slit with open ends (no cap). This precision is necessary in order to test the logarithmic growth predicted for the capped system. Our numerical results for the two dimensional DFT corresponding to the capped capillary are then presented and discussed.

Density Functional Theory

Within the framework of density functional theory for classical systems [1,[START_REF] Evans | Fundamentals of In homogeneous Fluids[END_REF] there exists the grand canonical potential functional of the form

Ω[ρ] = F[ρ] + d 3 r ρ(r)(V ext (r) -µ), (2)
where F[ρ] is the functional of the intrinsic Helmholtz free energy, V ext (r) is the external and µ the chemical potential. It can be proven that the functional Ω[ρ] possesses the property that it is minimal for the equilibrium density distribution ρ o (r), for which it reduces to the grand potential Ω = Ω[ρ 0 ] of the system. It follows that the equilibrium density distribution can be obtained from the variational principle

δΩ[ρ] δρ(r)) ρ(r)=ρ0(r) = 0. ( 3
)
Once ρ 0 (r), i.e. the inhomogeneous structure, is known the thermodynamics is accessible through the grand potential Ω. The functional of the intrinsic Helmholtz free energy

F[ρ] = F id [ρ] + F ex [ρ] (4)
can be split into an exactly known ideal gas part F id [ρ] and an excess (over the ideal gas) free energy F ex [ρ], which for most systems of interest is known only approximately. The excess free energy contains all the information about the inter-particle interactions. Here we wish to study a simple fluid with short-ranged interactions, that can phase separate into a liquid and a vapour phase. A well studied case is the square-well fluid [START_REF] Evans | [END_REF]25] with the inter-particle interaction given by As usual β = 1/(k B T ), with the Boltzmann constant k B and the absolute temperature T . Within DFT we treat the square-well interaction by a hard-sphere reference system, described by the White-Bear version [5][6][7] of fundamental-measure theory [4,26], and a square-well attraction as a perturbation:

βV sw (r) =    ∞ r < 2R -βε 2R < r < 2R sw 0 otherwise. (5)
F ex [ρ] = F HS ex [ρ] + 1 2 d 3 r ρ(r) d 3 r ρ(r ) φ sw (|r -r |). (6)
The perturbation term of this form underestimates the correlation in the system. To compensate this effect, one usually introduces an empirical modified square-well potential

βφ sw (r) = -βε r < 2R sw 0 otherwise, (7)
where the attraction is extended into the core, i.e. to r → 0.

In the absence of any external field, V ext (r) = 0, the density profile of the fluid reduces to the constant bulk density ρ(r) = ρ and the density functional provides simple expressions for the chemical potential µ(ρ, T ) and the pressure p(ρ, T ) of the fluid. It is straightforward to construct the fluid phase diagram by demanding mechanical-and chemical-equilibrium between two fluid phases with densities ρ I and ρ II , respectively, i.e.

µ(ρ I , T ) = µ(ρ II , T ) and p(ρ I , T ) = p(ρ II , T ). ( 8
)
In the following we chose the range of the square-well attraction to be R sw = 3σ, where σ is the hard-sphere diameter. The resulting phase diagram of a square-well bulk fluid is shown in Fig. 2 in the (dimensionless) density-temperature representation, where η = ρ(4π/3)R 3 is the packing fraction of the bulk fluid. Note that for a confined fluid, as studied in the following the packing fraction is a measure for the density in the particle reservoir. The DFT calculations presented here are performed at a fixed temperature, corresponding to βε = 1.1, as indicated by the dotted line in Fig. 2. The dashed line in Fig. 2 denotes the Fisher-Widom line [27], which separates region in the phase diagram in which the bulk correlation function decays monotonic from those in which it decays oscillatory. Within the framework of DFT Bob Evans and co-workers could show [28][29][30][31][32] that the character of the bulk correlations, i.e. monotonic or oscillatory decay, is also reflected by inhomogeneous density distributions. It is this insight that made us calculate and plot the Fisher-Widom line for our system in Fig. 2. The fluid state points we consider are clearly on the oscillatory side of the Fisher-Widom line.

Before studying the square-well fluid confined by a capped capillary, for which the density profile ρ(r) = ρ(x, z) depends on two coordinates, it is of great importance to establish the accuracy and internal consistency of the implementation of the DFT in order to ensure that the numerical results are reliable. For our study we require consistency between the structure of the fluid, i.e. the density profile ρ(r), and thermodynamic quantities that can be derived from the grand potential Ω = Ω[ρ(r)]. To this end we verify that our numerical results at a single planar hard wall satisfy the contact theorem and Gibbs' adsorption theorem. Note that we have employed the two-dimensional DFT implementation even at a single planar wall, although the problem could be reduced to an effectively one-dimensional problem. Furthermore, we have verified that our implementation satisfies Gibbs' adsorption theorem in a slit geometry and in a capped capillary. Using the two dimensional implementation with 10 grid points per hard-sphere radius we have established at the fixed temperature shown in Fig. 2 for various densities that the relative error of the sum rules is less than 10 -5 .

Infinite Slit

It is well known and studied that for all values of T between the critical temperature T c and the triple point temperature T t at a single planar hard wall complete drying occurs for a square-well liquid as bulk coexistence with the vapour phase is approached, i.e. as the chemical potential µ → µ sat (T ) approaches from above its value at coexistence, µ sat (T ), a macroscopically thick vapour film wets the interface between the bulk liquid and the planar hard wall. It is also well established that a square-well liquid off coexistence, µ > µ sat (T ), confined by two parallel planar hard wall (slit) at separation L undergoes the phase transition of capillary evaporation (CE) as the chemical potential µ → µ CE (T, L) approaches the value µ CE (T, L) > µ sat (T ) from above [12]. At µ = µ CE (T, L) a high density liquid and a low density vapour, that is meta-stable in the bulk, can coexist in a slit of width L.

The slit enters our DFT calculation via the external field

βV slit ext (x) = ∞ |x| > L 2 0 otherwise. (9)
In order to determine the location of the capillary evaporation we perform a series of DFT calculations for a fixed temperature, which corresponds to a dimensionless square-well depth of βε = 1.1, as indicated by the dotted line in Fig. µ > µ CE (T, L) we observe a high density liquid in the slit, while for µ < µ CE (T, L) the density profile jumps from a high density liquid profile to a low density vapour one [12]. However, the jump in the density profile is not sufficient to locate the capillary evaporation phase transition, because of hysteresis effects [12]. Within DFT it is possible to calculate the grand potential of the system in addition to the density profile. In Fig. 3 we show the grand potential Ω = Ω[ρ(x)] as a function of the reservoir packing fraction obtained from DFT for a square-well fluid in a slit with fixed width. As can be seen, there are two distinct branches, one corresponding to a liquid (squares, full line) and one corresponding to a vapour (circles, dashed line). The grand potential of the liquid depends much stronger on the reservoir packing fraction than that of the vapour. The two branches of the grand potential intersect at the capillary evaporation phase transition, which can be described by η CE (T, L) or equivalently µ CE (T, L). We find for our choice of parameters the capillary evaporation to correspond to a reservoir packing fraction of η CE (T, L) = 0.360641, which is marked by the circle in Fig. 2. The temperature is chosen so that βε = 1.1 (see dotted line in Fig. 2). The grand potential of the liquid branch (squares) and that of the vapour branch (circles) intersect at the point of capillary evaporation, corresponding to a reservoir packing fraction of η CE (T, L) = 0.360641 (circle in Fig. 2).

In the inset of Fig. 3 the grand potential of both the liquid and the vapour branch is shown in the vicinity of capillary evaporation. It is clear that there are hysteresis effects and meta-stable liquids or vapours can be observed in the slit, as is typical for a first order phase transition. Note that the accurate location of capillary evaporation can only be determined, if the theory describes the structure (density profiles) and thermodynamics (grand potential) of the system consistently. It is for this reason that we have verified that the contact theorem at a single wall and the Gibbs adsorption theorem at a single wall and in the slit are satisfied.

For a reservoir packing fraction of η = η CE (T, L) = 0.360641 a liquid and a vapour can coexist in the slit for the given parameters, L = 19.6 and βε = 1.1. We plot he corresponding density profiles of the coexisting phases in Fig. 4. The density profile of the liquid possesses an oscillatory high density part in the middle of the slit, reflecting the fact that its state point is located in the oscillatory region of the phase diagram (see Fig. 2), and has a contact value at the walls, |x| = L/2, which is significantly lower than the reservoir density [START_REF] Evans | [END_REF]. Note that the state point of the high density liquid lies on the oscillatory side of the Fisher-Widom line, while that of the low density vapour is located on the monotonic side. In agreement with this observation we find an oscillatory structure only in the liquid profile. As the capillary evaporation is crossed the stable state of the fluid in the slit jumps between these two profiles.

Since we have chosen a temperature well below the critical temperature the difference between the high density liquid and the low density gas at the same chemical potential is large. As a consequence, the change in the density profile as the point of capillary evaporation is crossed is big.

Capped Capillary

In order to study the behaviour of the square-well fluid in a capped capillary we change the external potential from the slit potential, Eq. ( 9), to that of a capped capillary

βV cap ext (x, z) = ∞ |x| > L 2 or z < 0 0 otherwise. ( 10
)
The geometry of the capped capillary is shown in Fig. 1(b). The length of the capillary, D is assumed to be infinite. Clearly, the resulting density profiles are now functions of x and z. The reservoir packing fractions considered here correspond to chemical potentials that are larger than its value at capillary evaporation, i.e. µ = µ CE (T, L) + δµ with δµ > 0. This choice is denoted by the full line at fixed temperature in Fig. 2, which ends at the location of capillary evaporation, marked by the circle.

A typical example of the two dimensional density profile of a square well fluid with a reservoir packing fraction of η = 0.38 in the capped capillary is shown in Fig. 5. The chemical potential for the liquid at η = 0.38 is clearly away from capillary evaporation, as is measured by βδµ = 0.3719.

For z → ∞, the profile recovers that of an infinite slit with a low contact value at the walls, |x| = L/2, and an oscillatory high liquid density part in the middle of the capillary. The oscillations represent packing effects in the liquid.

For z → 0 the influence of the third wall becomes apparent. wedges at z = 0 and |x| = L/2 is the smallest, as expected, but for the density used for Fig. 5 the wedges are still separated by a region along the third wall at z = 0 with higher density.

x/R z/R ρ(x, z)R 3 ρ(x, z)R 3
As we decrease the reservoir packing fraction η or the corresponding distance in chemical potential δµ from capillary evaporation, the density profile ρ(x, z) changes considerably. The first important change takes place at the third wall at z = 0. The low density regions at the wedges at |x| = L/2 connect to each other, as can be seen in Fig. 6(a)-(b). At the same time the oscillatory structure in the middle of the slit for z → ∞ becomes weaker. Note that also the contact value of the liquid at the walls at |x| = L/2 for z → ∞ decreases, but not nearly as fast as the contact value at z = 0 and x = 0. In Fig. 6(a) and (b) we show density profiles for reservoir packing fractions of η = 0.375 and 0.370, which corresponds to βδµ = 0.2659 and 0.1670, respectively.

As we further approach the location of capillary evaporation, one can observe that a vapour film wets the third wall at z = 0 [22], as is shown in Fig. 6(c) and (d), for reservoir packing fractions of η = 0.365 and 0.36067. The corresponding values of βδµ are 0.0749 and 0.00049, respectively.

While it is very interesting to get an impression of the behaviour of the density profile as capillary evaporation is approached, it is important to quantify this phenomena by thermodynamic quantities. Here we consider

Ω ex = Ω[ρ(x, z)] + pV, (11)
the excess (over a bulk liquid) grand potential. The grand potential of the corresponding bulk system is simply Ω bulk = -pV , which assumes as a reference system a liquid with constant density that occupies a volume V . Note that while the total grand potential of the system is well defined and unique, the excess grand potential is not because the volume V is not uniquely defined. Here we use the volume that is accessible to the centrers of fluid spheres. Strictly speaking, the thermodynamic limit of Eq. ( 11), in which the system size becomes infinite, is problematic because both Ω and Ω ex are infinite in this limit. Since we are using a fixed, and finite system size in our DFT calculations, we do not face this problem in practice. The finite system considered here is extended by the behaviour of the liquid in an infinite slit, i.e. we make use of the fact that ρ(x, z → ∞) → ρ slit (x), where ρ slit (x) is the density profile of the liquid at the same state point in a slit geometry. We have chosen the system size sufficiently large that this asymptotic behaviour is realised.

In order to isolate the effects of the third wall at z = 0 it would be cleaner to subtract from Ω the grand potential of a liquid confined by a slit of volume V , rather than that of a bulk liquid. However, the influence of the parallel walls for z → ∞ is rather small and for the following discussion of no importance so that we can stick to the definition of the excess grand potential given in Eq. (11). The line is a guide to the eye. Γ diverges for δµ → 0 like ln δµ.

Related to the excess grand potential is the excess adsorption defined by

Γ = V d 3 r (ρ(r) -ρ) = - ∂Ω ex ∂µ T,V (12)
which quantifies the excess (over a bulk system) amount of liquid found in the capped capillary. As indicated in Eq. ( 12) the excess adsorption can be either calculated from the inhomogeneous density profile ρ(x, z), or using Gibbs' adsorption theorem as a derivative of the excess grand potential, Eq. ( 11), w.r.t. the chemical potential µ. We have verified that for our results both routes are equivalent. Since the excess grand potential and the excess adsorption are related via Gibbs' adsorption theorem, the same remarks about the thermodynamic limit and the reference system apply.

We show the DFT results for Ω ex as a function of δµ in Fig. 7. The excess grand potential of the liquid in a capped capillary decreases as the limit δµ → 0 is approached. In the inset of Fig. 7 we show Ω ex for very small values of δµ. In the whole range we considered the excess grand potential is well behaved.

The situation is clearly different for the excess adsorption, shown in Fig. 8. Γ is negative and decreases in the limit δµ → 0, however, in contrast to the excess grand potential, the decrease in Γ becomes stronger with smaller values of δµ. In fact for δµ → 0 one can observe a divergence of the excess adsorption. This is highlighted in the inset of Fig. 8. In order to better understand the divergence we fit the DFT data. Clearly, a polynomial in δµ is not sufficient to capture the behaviour of Γ in the range of δµ shown in Fig. 8. However, we find, in accordance with the prediction by Rascon and Parry [13,22], that the excess adsorption can very nicely be fitted by a polynomial in δµ plus a logarithmic term, i.e. in the limit δµ → 0 we find that Γ ∝ -ln δµ. The corresponding term in the excess grand potential in the same limit is given by Ω ex ∝ δµ ln δµ, which is non-analytic but remains finite.

It is also interesting to consider cuts ρ(x = 0, z) of the two-dimensional density profiles at the middle of the capped capillary. The predictions of Rascon and Parry Molecular Physics cc [13,22] suggest that these cuts should essentially behave like wetting (or drying) profiles of the fluid at a single planar wall, when bulk coexistence is approach as δµ → 0. Some cuts are shown in Fig. 9. . The density ρ(x = 0, z) in the middle of the capped capillary, as capillary evaporation is approached with δµ = µµ CE (T, L) → 0. These profiles behave essentially like those of a fluid at a single planar wall, as bulk coexistence is approach with δµ = µµsat(T ) → 0.

Due to the close connection between the excess grand potential and the excess adsorption, we can test the consistency of the fits. The fit parameters for the excess adsorption can be used, up to a constant, to describe the data for the excess grand potential, including the non-analytic contribution.

It is important to realise that the fact that the divergence in the excess adsorption is (numerically) observed at δµ = 0 is a consequence of the high level of consistency between the structure and thermodynamic quantities in our implementation of DFT. Any slight discrepancy could result in an apparent shift of the divergence.

Discussion

In this paper we have used the White-Bear version of fundamental measure DFT to investigate the nature of capillary evaporation in a capped hard-wall capillary slit. Our numerically results support the prediction made using effective interfacial Hamiltonian models predict that in this geometry capillary evaporation is a secondorder transition related to complete wetting. This is perhaps most clearly seen in our final diagram (Fig. 9) which shows various density profiles (measured at the mid-point between the walls) as a function of the distance from the capped end on approaching capillary evaporation. These profiles resemble, qualitatively, those for drying films at a planar wall and show the growth of a thick "drying" layer as "coexistence" is approached. However in this case the "wall" is the capped end, coexistence refers to capillary evaporation, the unbinding "interface" is actually a cross-section of the meniscus, and the coexisting phases are the "liquid-like" and "vapour-like" phases of the shifted bulk coexistence curve. Using an accurate and independently determined value for µ CE , we have shown that the growth and divergence of the thickness of the vapour-like phase (the distance of the meniscus to the cap) as δµ = µ -µ CE is -ln δµ -the same as for complete wetting with short-ranged forces in accordance with the predictions of Evans and co-workers. Our present model DFT calculation is of course only mean-field like and does not include interfacial fluctuation effects. For the present case of systems with shortranged forces these could be large and eventually lead to cross-over behaviour in which the adsorption changes from its mean-field expression Γ ∝ -ln δµ to Γ ∝ δµ -1/3 . The latter is the universal result for complete wetting (in bulk two dimensions). Such cross-over behaviour would be very interesting to study as a function of the slit width. Our intuition is that the wider the slit, the narrower the true asymptotic regime and the broader the regime when the mean-field calculation is valid. Unfortunately none of this is as yet possible using DFT and one must still resort to coarse-grained interfacial models and simulation studies.

However there are still many interesting questions one can ask using the present DFT and substrate geometry. For example it would be interesting, and relatively simple, to add an attractive potential to the wall-fluid interaction in order to alter the contact angle. Indeed one may use a different material to cap the capillary so that the interaction potential with the side walls and cap is different. By doing this we might be able to alter the manner in which the meniscus unbinds from the capped end similar to studies of tri-criticality filling transitions in linear wedges. Similarly it would be interesting to see what happens as one moves along the capillary-evaporation line towards the capillary-critical point. In this case we can expect cross-over to critical-adsorption-like behaviour in which the decay of the profile from the capped end follows a slow algebraic power law. Finally we hope that our present work serves to encourage further density functional model studies of adsorption on structured surfaces including more complex fluids and their mixtures.

DFigure 1 .

 1 Figure 1. The geometry of an infinite slit (a) of two parallel hard walls which are separated by the width L and of a capped capillary (b) with width L. In case (b) the length D of the capillary is assumed to be infinite.

  sence of fluctuation effects, this implies that the adsorption in the capped capillary should diverge as Γ ≈ -ln |p -p C (T, L)|.

Figure 2 .

 2 Figure 2. Bulk phase diagram of a square-well fluid with well depth ε and square-well range Rsw = 3σ.The dashed line is the Fisher-Widom line that separates the phase diagram into a region in which the bulk correlation function decays monotonic and one in which it decays oscillatory. We perform our DFT calculations at a fixed temperature, corresponding to βε = 1.1, as indicated by the dotted line, i.e. in the oscillatory region of the phase diagram. For this temperature and the employed slit width, L = 19.6R, capillary evaporation occurs at η CE (T, L) = 0.360641 (marked by the circle), well away from bulk phase coexistence. The full line indicates the path in the phase diagram taken for our calculations in the capped capillary.

Figure 3 .

 3 Figure 3. Grand potential per unit area Ω of a square-well fluid in an infinite slit. The temperature is chosen so that βε = 1.1 (see dotted line in Fig.2). The grand potential of the liquid branch (squares) and that of the vapour branch (circles) intersect at the point of capillary evaporation, corresponding to a reservoir packing fraction of η CE (T, L) = 0.360641 (circle in Fig.2).

Figure 4 .

 4 Figure 4. Density profiles of a square-well fluid in an infinite slit at the coexisting reservoir packing fraction of η = η CE (T, L) = 0.360641. The difference between the liquid and the vapour profile is well pronounced. Note that the state point of the high density liquid lies on the oscillatory side of the Fisher-Widom line, while that of the low density vapour is located on the monotonic side. In agreement with this observation we find an oscillatory structure only in the liquid profile. As the capillary evaporation is crossed the stable state of the fluid in the slit jumps between these two profiles.

Figure 5 .

 5 Figure 5. Density profile of a square-well fluid in a capped capillary. The packing fraction is η = 0.380. Corresponding to this packing fraction βδµ = 0.3719.

Figure 6 .Figure 7 .

 67 Figure 6. Density profiles of a square-well fluid in a capped capillary. The packing fractions are (a) η = 0.375, (b) 0.370, (c) 0.365, and (d) 0.36067, respectively. These values of η correspond to (a) βδµ = 0.2659, (b) 0.1670, (c) 0.0749, and (d) 0.00049, respectively.

16 Figure 8 .

 168 Figure 8. The excess adsorption Γ of a square-well fluid in a capped capillary as a function of δµ = µ-µ CE .The line is a guide to the eye. Γ diverges for δµ → 0 like ln δµ.

Figure 9

 9 Figure 9. The density ρ(x = 0, z) in the middle of the capped capillary, as capillary evaporation is approached with δµ = µµ CE (T, L) → 0. These profiles behave essentially like those of a fluid at a single planar wall, as bulk coexistence is approach with δµ = µµsat(T ) → 0.

  ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¡£ ¤£ ¡£ ¡£ ¡£ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥ ¡¥ ¤¥ ¡¥ ¡¥ ¡¥

	¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢	£

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

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

Acknowledgments

We are just two of the masses who were taught how to calculate by Bob (although any mistakes in this article are entirely our own). It is our heartfelt pleasure to thank him for his deep physical insight and the inexhaustible encouragement he has always given us.