The van Hove distribution function for Brownian hard spheres: dynamical test particle theory and computer simulations for bulk dynamics
Paul Hopkins, Andrea Fortini, Andrew Archer, Matthias Schmidt
aa r X i v : . [ c ond - m a t . s o f t ] O c t The van Hove distribution function for Brownian hard spheres:dynamical test particle theory and computer simulations for bulk dynamics
Paul Hopkins, ∗ Andrea Fortini, Andrew J. Archer, † and Matthias Schmidt
1, 2 H.H. Wills Physics Laboratory, University of Bristol, Tyndall Avenue, Bristol BS8 1TL, UK Theoretische Physik II, Physikalisches Institut Universit¨at Bayreuth, D-95440 Bayreuth, Germany Department of Mathematical Sciences, Loughborough University, Loughborough LE11 3TU, UK (Dated: October 12, 2010)We describe a test particle approach based on dynamical density functional theory (DDFT) forstudying the correlated time evolution of the particles that constitute a fluid. Our theory providesa means of calculating the van Hove distribution function by treating its self and distinct parts asthe two components of a binary fluid mixture, with the ‘self’ component having only one particle,the ‘distinct’ component consisting of all the other particles, and using DDFT to calculate the timeevolution of the density profiles for the two components. We apply this approach to a bulk fluidof Brownian hard spheres and compare to results for the van Hove function and the intermediatescattering function from Brownian dynamics computer simulations. We find good agreement atlow and intermediate densities using the very simple Ramakrishnan-Yussouff [Phys. Rev. B ,2775 (1979)] approximation for the excess free energy functional. Since the DDFT is based on theequilibrium Helmholtz free energy functional, we can probe a free energy landscape that underliesthe dynamics. Within the mean-field approximation we find that as the particle density increases,this landscape develops a minimum, while an exact treatment of a model confined situation showsthat for an ergodic fluid this landscape should be monotonic. We discuss possible implications forslow, glassy and arrested dynamics at high densities. PACS numbers: 05.20.Jj, 61.20.-p, 64.70.Q-
I. INTRODUCTION
The structure of condensed matter is commonlyprobed by X-ray or neutron scattering techniques, thatyield quantities such as S ( k ), the static structure factor,and its dynamical counterpart F ( k, t ), the intermediatescattering function . However, in recent years, with theadvent of modern confocal microscopes, which are ableto characterise the structure of colloidal suspensions inreal-space, there is an equally great emphasis on deter-mining the radial distribution function g ( r ) and its dy-namical counterpart G ( r, t ), the van Hove distributionfunction . The van Hove distribution function G ( r, t )is a real-space dynamical correlation function for char-acterising the spatial and temporal distributions of pairsof particles in a fluid. It gives the probability of find-ing a particle at position r at time t , where | r | = r ,given that one of the particles was located at the ori-gin at time t = 0. The intermediate scattering function F ( k, t ) is simply obtained from G ( r, t ) via spatial (three-dimensional) Fourier transform. Pair correlation func-tions are important because of the significant amountof information that they contain: Transport coefficientscan be calculated via Kubo formulae – for example thediffusion coefficient D can be obtained from G ( r, t ) –and thermodynamic quantities such as the internal en-ergy and the pressure can be related to spatial integrals involving g ( r ). Whether or not a liquid is near to freez-ing can often be discerned from inspecting the heightof the principal peak in S ( k ): it was first noticed byHansen and Verlet that many simple liquids freeze whenthe principal peak in S ( k ) at k = k m , takes the value S ( k m ) ≃ .
85. Whether a system is a glass (i.e. an amor-phous solid) rather than a fluid may be determined fromthe long time limit value of F ( k, t ), because in a fluid the t → ∞ limit of F ( k, t ) is zero, whereas for a glass thislimit takes non-zero values. This brief (and incomplete)survey is intended to demonstrate that both dynamicaland static pair correlation functions are fundamental forcharacterising and understanding liquids.In the history of liquid state physics, fluids of hardspheres have proved to be an important model systemfor developing new techniques and theories. The hardsphere model is composed of particles interacting via thepair potential v hs ( r ) = (cid:26) ∞ r < σ , (1)where r is the distance between the centers of the par-ticles and σ is the hard sphere diameter. Hard spheresplay an important role in describing real systems, be-cause attractive interactions such as those present in theLennard-Jones potential, can often be treated as a per-turbation to the hard sphere system . Hence a theorythat can successfully describe the properties of the hardsphere fluid forms a good candidate to work for morerealistic systems. The hard sphere model has furthergrown in importance in recent decades due to the factthat Eq. (1) provides a good model for the effective in-teraction potential between colloidal particles in suspen-sion, in the case when the charges on the colloids aresmall or well screened – see e.g. Ref. 5 for an example ofsuch a system. As the density of a hard sphere fluid isincreased, the system freezes to a crystalline state, andthe hard sphere model has played an important role indeveloping our understanding of this phase transition.In contrast, although the glass transition has attractedmuch interest in recent years, it is still not completely un-derstood. An introduction to the vast literature on thissubject can be found in Refs. 1,6 and references therein.A number of universal processes have been discovered, in-cluding dynamical heterogeneity , stretched exponentialdecay of correlation functions , and two-stage relaxationtimes . In order to understand the processes involvedin structural arrest, a number of different theoretical ap-proaches have been used. In particular, mode-couplingtheory (MCT) has been successful in describing the bulkglass transition for hard sphere colloids , and has beenapplied e.g. to suspension rheology . Nevertheless,alternatives to MCT have been developed . What isclear from the many studies of arrested systems, is thatkey signatures of the slow dynamics are manifest in dy-namical pair correlation functions.In our previous Rapid Communication , a theory tocalculate the van Hove function was proposed. For abulk fluid of particles interacting via Gaussian pair po-tentials, comparison with Brownian dynamics (BD) com-puter simulation results showed that the theory is veryreliable for determining G ( r, t ) for this particular modelsystem. The theory is formulated for inhomogeneoussystems and hence was also applied to investigate thedynamics of hard spheres confined between two parallelhard walls . This approach has since been applied toinvestigate dynamics in liquid crystalline systems . Inthe current paper we explore the theory further, and ap-ply it to study a bulk fluid of Brownian hard spheres.We present results for the self and distinct parts of thevan Hove distribution function, G s ( r, t ) and G d ( r, t ) re-spectively, and by Fourier transforming, for the inter-mediate scattering function F ( k, t ). We also displayresults for the scaled intermediate scattering function φ ( k, t ) ≡ F s ( k, t ) /F s ( k, t = 0) evaluated at the wavenumber kσ = 2 π . This function is often the centralobject of focus of MCT . The two stage relaxation of φ ( k, t ), that MCT predicts close to the glass transition,is also present in our theory. We also determine G ( r, t )and F ( k, t ) using BD computer simulations and comparethese results with those from the theory.Our starting point is a dynamical generalisation of Per-cus’ test-particle approach for determining the radialdistribution function g ( r ). Percus showed that for a fluidof classical particles interacting via the pair interactionpotential v ( r ), that if one sets the external potential act-ing on the fluid u ( r ) = v ( r ), then the one body den-sity distribution ρ ( r ) of the fluid around the fixed ‘test’particle is equal to the radial distribution function, mul-tiplied by the bulk fluid density ρ ; i.e. Percus showedthat ρ ( r ) = ρg ( r ). When using equilibrium density func-tional theory (DFT) to study a fluid, the test-particlemethod provides a useful route to obtain g ( r ) because u ( r ) (and hence v ( r )) enters the framework explicitly.We also describe an alternative ‘zero-dimensionality’ ap- proach for calculating g ( r ). This forms a stepping stonein the development of the dynamical theory.We apply a dynamical extension of Percus’ idea,together with dynamical density functional theory(DDFT) in order to calculate the van Hove function G ( r, t ) in general inhomogeneous situations. We imple-ment the method using the very simple Ramakrishnan-Youssouff (RY) approximation for the Helmholtz free en-ergy functional . We find that the results from the the-ory agree well with those from BD computer simulationswhen the fluid density ρσ . .
6. At higher densitiesthe free energy underlying the dynamics develops a min-imum, corresponding to the appearance of a free energybarrier that must be traversed for a particle to escapefrom the cage formed by the neighbouring particles.In addition, we compare our results for G ( r, t ) to thoseobtained by assuming that G s ( r, t ) takes a simple Gaus-sian form for all times t , together with the Vineyardapproximation , which sets G d ( r, t ) to be a simple con-volution of G s ( r, t ) and the radial distribution function g ( r ), as described in detail below. We find that in con-trast to the received wisdom , the simple Vineyard ap-proximation is actually a fairly good approximation forthe van Hove function for Brownian hard sphere fluids atlow and intermediate densities.We also compare to an equilibrium DFT based ap-proach with which we are able to calculate a series of den-sity profiles that agree well with those from the DDFT.This is done by performing a constrained minimisationof the free energy through the judicious choice of a suit-able external potential to confine the test particle. Thisapproach is easier to implement than the full DDFT andallows for the free energy landscape underlying the dy-namics to be mapped out and examined in detail. How-ever, this approach does not give any of the time infor-mation that the full DDFT gives – i.e. it yields the vanHove function with the time labels ‘removed’. One ofthe advantages of this approach is that for a particular(parabolic) choice of external potentials, we are able tocalculate exactly the fluid density profiles, which are pre-cisely those predicted by Vineyard’s theory . We dis-cuss the significance of this result below, after we havelaid out the general structure of the theory and shownthe results.This paper is structured as follows: In Sec. II we out-line the necessary theoretical background, including thedefinition of the van Hove function, the Vineyard approx-imation, DDFT, and the static test particle limit. Mostof this section may be safely skipped by expert readers.In Sec. III the dynamical test particle limit is introduced.Sec. IV summarises the model used and describes thesimulation details. In Sec. V we describe results from thedifferent dynamical approaches, the corresponding equi-librium approaches, and the free energy landscape. InSec. VI we make some concluding remarks. Appendix Apresents an exact solution of a corresponding equilibriumsituation. II. BACKGROUNDA. The van Hove function
We first recall the definition of the van Hove functionand some of its properties; for a more detailed accountsee Refs. 1,2. Consider a set of N particles with timedependent position coordinates r i ( t ), where i = 1 , .., N is the particle index, and t is time. The van Hove cor-relation function is defined as the probability of findinga particle at position r at time t , given that there was aparticle at the origin at time t = 0: G ( r, t ) = 1 N * N X i =1 N X j =1 δ ( r + r j (0) − r i ( t )) + , (2)where h·i represents an ensemble average and δ ( · ) is thethree-dimensional Dirac delta function. G ( r, t ) can benaturally separated into two terms, conventionally re-ferred to as its “self” and “distinct” part, by discrimi-nating between the cases i = j and i = j , respectively.So G ( r, t ) = 1 N * N X i =1 δ ( r + r i (0) − r i ( t )) + + 1 N * N X i = j δ ( r + r j (0) − r i ( t )) + ≡ G s ( r, t ) + G d ( r, t ) , (3)where the self part, G s ( r, t ), describes the average motionof the particle that was initially at the origin, whereasthe distinct part, G d ( r, t ), describes the behavior of theremaining N − t = 0, Eq. (3) reducesto the static particle-particle auto-correlation function,which is defined as G ( r ,
0) = δ ( r ) + 1 N * N X i = j δ ( r + r j (0) − r i (0)) + = δ ( r ) + ρg ( r ) , (4)where g ( r ) is the (static) pair distribution function. Forthe homogeneous bulk fluid ρ ( r ) = ρ ; isotropy impliesthat the dependence is only on r = | r | . Thus, at t = 0: G s ( r,
0) = δ ( r ) (5) G d ( r,
0) = ρg ( r ) . (6)From the definitions of G s ( r, t ) and G d ( r, t ), Eq. (3), itis clear that the volume integral of these functions mustbe a conserved quantity for all times t : Z d r G s ( r, t ) = 1 , (7) Z d r G d ( r, t ) = N − . (8) The asymptotic behaviour of G ( r, t ) in bulk in the ther-modynamic limit is obtained by considering N → ∞ andvolume V → ∞ such that N/V = ρ is finite:lim r →∞ G s ( r, t ) = lim t →∞ G s ( r, t ) = 0 , (9)lim r →∞ G d ( r, t ) = lim t →∞ G d ( r, t ) = ρ. (10)A key quantity that we use below to characterise G s ( r, t )is its width w ( t ) defined via( w ( t )) = 4 π Z ∞ d r r G s ( r, t ) , (11)the second moment of G s ( r, t ). It is often convenientto consider the intermediate scattering function which isrelated to the van Hove function via a spatial Fouriertransform, F ( k, t ) = Z d r G ( r, t ) exp( − i k · r ) . (12)This quantity is directly accessible in light and neutronscattering experiments . B. Approximating G s ( r, t ) A commonly used approximation for the self part ofthe van Hove function is to assume a Gaussian shape : G s ( r, t ) = 1 π / W ( t ) exp (cid:18) − r W ( t ) (cid:19) (13)where the width W ( t ) = q w ( t ), when w ( t ) is calculatedvia Eq. (11). The form (13) is exact in the limits t → t → ∞ for all densities when the system is fluid .It is also exact for all times t in the low density limit ρ →
0, where interactions between the particles can beneglected. There are a number of approximations for W ( t ). In molecular dynamics, at very short times t ≪ τ c ,where τ c is the mean collision time, particles in a fluid donot experience collisions and therefore move freely at aconstant velocity. This is akin to an ideal gas where theparticle velocities follow a simple Maxwellian (Gaussian)distribution, giving W ( t ) = t p /βm, (14)where m is the particle mass. Over longer times t ≫ τ c the particles in the fluid undergo many collisions withneighbouring particles, so that the trajectory of a givenparticle is a random walk and thus its probability distri-bution G s ( r, t ) is the solution of the diffusion equation ∂G s ( r , t ) ∂t = D l ∇ G s ( r , t ) , (15)where D l is the long time self diffusion coefficient. For theDirac delta initial condition (5), the solution of Eq. (15)corresponds to the Gaussian form (13), with W ( t ) = 2 p D l t. (16)For colloidal particles, the collisions with the solventatoms happen so frequently that the time scale τ c is muchsmaller than all other time scales relevant for the dynam-ics, such as the Brownian time scale τ B which is roughlythe time for a particle to travel a distance equal to itsown diameter, and also the typical time scale betweencollisions of pairs of colloids, τ col . This means that wemay set τ c → t . Thus, we maycombine Eqs. (13) and (16) to obtain G s ( r, t ) = (4 πD l t ) − / exp (cid:18) − r D l t (cid:19) . (17)We find below that for Brownian hard spheres this ap-proximation is not only reliable in the low density limit,but is also fairly good up to intermediate densities ρσ . . C. Vineyard approximation for G d ( r, t ) Vineyard suggested that one may rewrite the distinctpart of the van Hove function as G d ( r, t ) = Z d r ′ g ( r ′ ) H ( r , r ′ , t ) , (18)which is merely a redefinition of G d ( r, t ) in terms of theunknown function H ( r , r ′ , t ), which is the probabilitythat if there was a particle at the origin at time t = 0and a second particle located at r ′ , this second particleis later located at r at time t . Vineyard’s approximationis to replace H ( r , r ′ , t ) by G s ( r − r ′ , t ), giving G d ( r, t ) = Z d r ′ g ( r ′ ) G s ( r − r ′ , t ) . (19)Some comments in the literature state that the Vine-yard approximation ignores important correlations thatinhibit the rate at which the structure of the liquidbreaks up and it therefore predicts too rapid decay ofthis structure . This may indeed be the case for fluidsundergoing molecular dynamics, but for the system withBrownian dynamics (over-damped stochastic equationsof motion) that we consider here, we find that taking Eq.(17) together with Eq. (19), is actually fairly reliable –in particular when the fluid is at low and intermediatedensities ρσ . .
6. We will henceforth refer to Eqs.(17) and (19) as the ‘Vineyard approximation’ for thevan Hove function.
D. DDFT and equilibrium DFT
The dynamics of a system of N Brownian (colloidal)particles with positions r i ( t ) can be modelled with thefollowing set of (over-damped) stochastic equations ofmotion : Γ − d r i ( t )d t = −∇ i U N ( r N , t ) + ζ i ( t ) , (20) where r N = { r i ; i = 1 , . . . , N } is the set of particle co-ordinates, Γ − is a friction constant characterizing theone-body drag of the solvent on the particles, ζ i ( t ) is astochastic white noise term and the total inter-particlepotential energy is U N ( r N , t ) = N X i =1 u ( r i , t ) + 12 N X i =1 X j = i v ( | r i − r j | ) , (21)which is composed of a one-body contribution due tothe external potential u (which may or may not betime dependent), and a sum of contributions from thepair interactions between the particles. The time evo-lution of the probability density for the particle coor-dinates P ( N ) ( r N , t ) is described by the Smoluchowskiequation : ∂P ( N ) ∂t = Γ N X i =1 ∇ i · [ k B T ∇ i P ( N ) + ∇ i U N P ( N ) ] . (22)The one-body density is obtained by integrating over theposition coordinates of all but one particle: ρ ( r , t ) = N Z d r ... Z d r N P ( r N , t ) . (23)Integrating the Smoluchowski equation (22) we obtainthe key equation of DDFT : ∂ρ ( r , t ) ∂t = Γ ∇ · (cid:20) ρ ( r , t ) ∇ δF [ ρ ( r , t )] δρ ( r , t ) (cid:21) , (24)where F [ ρ ] is taken to be the equilibrium total Helmholtzfree energy functional: F [ ρ ( r )] = k B T Z d r ρ ( r )[ln(Λ ρ ( r )) − F ex [ ρ ( r )] + Z d r u ( r ) ρ ( r ) , (25)where the first term on the right hand side is the ideal-gas contribution to the free energy, Λ is the thermal deBroglie wavelength, F ex [ ρ ( r )] is the excess (over ideal gas)part of the free energy, which is in general unknown ex-actly, and we have suppressed the dependence on temper-ature T and volume V in the notation. In obtaining (24)we have made the approximation that equal-time two-body correlations at each time t in the non-equilibriumsituation are the same as those of an equilibrium fluidwith the same one-body density profile ρ ( r , t ), generatedby an appropriate external potential . It has beenshown in a variety of cases that the DDFT (24) is reliablein predicting the time-evolution of ρ ( r , t ), when solvedin conjunction with a sufficiently accurate approxima-tion for the equilibrium Helmholtz free energy functional F [ ρ ( r )] – see for example the results presented in Refs.21,28–36.Although in the following we will not go beyond dy-namics that are local in time, it is worth mentioning thatmore generally, going beyond the case of particles withstochastic over-damped equations of motion (20), Chanand Finken established a rigorous DDFT for classicalfluids, showing that the time evolution of the one bodydensity ρ ( r , t ) is obtained from the solution of ∂ρ ( r , t ) ∂t = −∇ · j ( r , t ) , (26) ∂ j ( r , t ) ∂t = P [ ρ ( r , t )] , (27)where j ( r, t ) is the particle current density, and Eq. (26)represents a continuity equation for the one-body density ρ ( r , t ). One should, of course, expect on general groundsfor the dynamical equations to be of this form –recall that Eq. (27) is the continuity equation. However,the functional P [ ρ ( r , t )] that governs the time evolutionof ρ ( r , t ), takes a form that depends on the equationsof motion of the particles – i.e. it depends on whetherthe particles evolve under Newtonian dynamics or havestochastic equations of motion such as (20).Due to the fact that in general the functional P [ ρ ( r , t )]in Eq. (27) is an unknown quantity, one is preventedfrom applying our DDFT approach for calculating dy-namic correlation functions, and we are restricted to thethe Brownian case (20) outlined above. The particularapproximation used in Eqs. (26) and (27) to obtain Eq.(24), is to assume the one particle current density to beof the form j ( r , t ) = − Γ ρ ( r , t ) ∇ δF [ ρ ( r , t )] δρ ( r , t ) . (28)Nevertheless, there is much active research aimed at go-ing beyond the simple overdamped case .In what follows we will relate the van Hove function tothe time evolution of the one-body density profiles of abinary mixture; we therefore require the multicomponentgeneralization of Eq. (24): ∂ρ i ( r , t ) ∂t = Γ ∇ · (cid:20) ρ i ( r , t ) ∇ δF [ { ρ i } ] δρ i ( r , t ) (cid:21) , (29)where F [ { ρ i } ] has the following form (cf. Eq. (25)) : F [ { ρ i } ] = k B T X i Z d r ρ i ( r )[ln Λ ρ i ( r ) − F e x [ { ρ i } ] + X i Z d r u i ( r ) ρ i ( r ) . (30)where the summations run over all species i . Given aninitial set of density profiles, { ρ i ( r , t = 0) } , we mayemploy the DDFT equations (29) and (30) to calculatethe full time evolution of the one-body density profiles ρ i ( r , t ).For completeness, we recall some of the key resultsfrom equilibrium . For a given set of (one-body) ex-ternal potentials { u i ( r ) } , the unique set of one-bodydensity profiles { ρ i ( r ) } are those which minimize the Helmholtz free energy of the system F [ { ρ i } ], subject tothe constraint that the average number of particles ofeach species R d r ρ i ( r ) = N i , is fixed. This is equivalentto an unconstrained minimization of the grand potentialfunctionalΩ[ { ρ i } ] = F [ { ρ i } ] − X i µ i Z d r ρ i ( r ) , (31)where the Lagrange multiplier µ i is the chemical poten-tial of species i . Minimization with respect to variationsin the density profiles yields the following set of Euler-Lagrange equations : δF [ { ρ i } ] δρ i ( r ) = µ i . (32)The Euler-Lagrange equations can be rewritten as ρ i ( r ) = Λ − exp (cid:20) βµ i − βu i ( r ) + c (1) i ( r [ { ρ j } ]) (cid:21) , (33)where c (1) i ( r ; [ { ρ j } ]) = − β δF ex [ { ρ j } ] δρ i ( r ) , (34)is the one-body direct correlation functional. The set ofdensity profiles that satisfy (33) minimize the free energyand are the equilibrium density profiles. When the equi-librium set of profiles { ρ i ( r ) } are substituted into (31),the grand potential Ω of the system is obtained. E. Percus’ test particle limit
Here we give a derivation of Percus’ (static) test par-ticle limit closely following Ref. 48. Consider a one com-ponent system such that the Helmholtz free energy, F ,is given by (25). We are interested in the change in ρ ( r )when the external potential is changed from the poten-tial u ′ ( r ) to the potential u ( r ). To this end we per-form a functional Taylor expansion of F ex [ ρ ] in powersof ∆ ρ ( r ) = ρ ( r ) − ρ ′ ( r ). For the sake of simplicity weconsider the change in going from u ′ ( r ) = 0 to a spheri-cally symmetric external potential u ( r ). The variable inthe Taylor expansion is then ∆ ρ ( r ) = ρ ( r ) − ρ b , where ρ b is the bulk density and the expansion of F ex [ ρ ] to secondorder in ∆ ρ ( r ) is F ex [ ρ ] = F ex [ ρ b ] + Z d r δF ex [ ρ ] δρ ( r ) (cid:12)(cid:12)(cid:12)(cid:12) ρ b ∆ ρ ( r )+ 12 Z d r Z d r ′ δ F ex [ ρ ] δρ ( r ) δρ ( r ′ ) (cid:12)(cid:12)(cid:12)(cid:12) ρ b ∆ ρ ( r )∆ ρ ( r ′ )+ O (∆ ρ ) . (35)Although the form of F ex is not specified, the functionalderivatives are related to identifiable properties of thesystem, so that evaluating them at ρ ( r ) = ρ b gives δF ex [ ρ ] δρ ( r ) (cid:12)(cid:12)(cid:12)(cid:12) ρ b = − k B T µ ex ,δ F ex [ ρ ] δρ ( r ) δρ ( r ′ ) (cid:12)(cid:12)(cid:12)(cid:12) ρ b = − k B T c (2) ( | r − r ′ | ) , (36) δO (∆ ρ ) δρ ( r ) = B ( r ) , where µ ex is the excess chemical potential, c (2) ( r ) is the(pair) direct correlation function, and B ( r ) is an un-known function that contains the higher order terms ofthe Taylor expansion. Substituting Eqs. (35) and (36)into (25), and then minimizing the functional with re-spect to variations in ρ ( r ), we obtain the following Euler-Lagrange equation [c.f. Eq. (32)]: δFδρ ( r ) = µ = k B T ln(Λ ρ b ) + µ ex , (37)where we have separated the chemical potential µ into anideal-gas and an excess (over ideal-gas) contribution µ ex .In the case of a spherically symmetric external potential,Eq. (37) may be rewritten as: ρ ( r ) ρ b = exp h − βu ( r ) + Z d r ′ c ( | r − r ′ | )∆ ρ ( r ′ )+ B ( r ) i . (38)For the same one-component system the bulk Ornstein-Zernike (OZ) equation for the total correlation function, h ( r ) = g ( r ) −
1, reads as follows: h ( r ) = c ( r ) + ρ b Z d r ′ h ( r ′ ) c ( | r − r ′ | ) . (39)It can be shown through diagrammatic-methods thatthe OZ equation has the general solution h ( r ) = c ( r ) + ln( g ( r )) + βv ( r ) + b ( r ) (40)where v ( r ) is the inter-particle pair potential, and b ( r )is the bridge function composed of the sum of all the so-called ‘bridge’ diagrams . Substituting (40) into (39) weobtain: g ( r ) = exp h − βv ( r )+ Z d r ′ c ( | r − r ′ | ) ρ b h ( r ′ )+ b ( r ) i . (41)If we compare Eqs. (38) and (41) we find they have thesame structure, and that they may be formally identified.If we set u ( r ) = v ( r ) in Eq. (38), it can be shown throughdiagrammatic methods that b ( r ) = B ( r ) and that g ( r ) = ρ ( r ) /ρ b , (42)or alternatively, ρh ( r ) = ∆ ρ ( r ). Thus when u ( r ) = v ( r ),Eqs. (38) and (41) become identical. Therefore we notethat not only can the OZ relationship be derived from the free energy functional , but that the equilibriumone-body density profile in the presence of an externalpotential u ( r ) = v ( r ) is related to the (two-body) ra-dial distribution function via Eq. (42). We should recallat this point that although many formal statements canbe made about the bridge function b ( r ), in practice it isan unknown function, and all theories for g ( r ) constitutesome form of approximation for b ( r ) . For example, ifwe set b ( r ) = B ( r ) = 0 then (38) is equivalent to usingthe hyper-netted chain (HNC) approximation to the OZequation. Furthermore, Percus showed that by Taylorexpanding with different functions of ρ ( r ) one may re-trieve the Percus-Yevick and other closures to the OZequations. This result may also be generalized to inho-mogeneous systems . F. Zero-dimensionality route to g ( r ) We present an alternative method for calculating g ( r ),although its basis is the same key idea that underpinsPercus’ test particle limit described above: that g ( r ) canbe obtained from the density profile of a fluid around afixed test particle. The key difference is that instead oftreating the test particle as a fixed external potential,in the zero-dimensionality route we treat the test parti-cle via its density distribution. The density profile of aparticle fixed at a point (i.e. in zero-dimensional space –hence our choice of name for this limit) takes the form ofa Dirac delta function. Having fixed this contribution tothe density distribution [c.f. Eq. (4)], one can then cal-culate the density distribution of the remaining particlesin the presence of the test particle. Specifically, we canwrite the grand potential functional as:Ω ∗ [ ρg ( r )] = F id [ ρg ( r )] + F ex [ δ ( r ) + ρg ( r ))] − µ Z d r ρg ( r ) , (43)where ρg ( r ) is the density distribution of the remainingparticles – the quantity we wish to calculate. Note thathere, and in what follows, ρ is the bulk density. The idealgas term F id does not contain the Dirac delta contribu-tion – we have crossed over to a system with N − ρ is necessarily fixed,we must simply minimise Ω ∗ with respect to variationsin g ( r ), giving the following Euler-Lagrange equation tobe solved for g ( r ): δ Ω ∗ δg ( r ) = 0 . (44)An alternative means of calculating g ( r ) is to treat thesystem as a binary mixture. The test particle (which welabel ‘ s ’), with density distribution ρ s ( r ) = δ ( r ), is onespecies and then we use the DFT for a binary mixtureto calculate the density profile of the remain particles ρ d ( r ) in the presence of the density profile ρ s ( r ) for thefixed particle, treating the remaining particles as a secondspecies ‘ d ’ in the mixture. ρ d ( r ) is the solution of δ Ω † δρ d ( r ) = 0 , (45)where Ω † is a modified version of Eq. (31) where ρ s ( r ) = δ ( r ) is fixed and therefore the ideal Helmholtz free energy, F id [ ρ d ], does not depend on ρ s ( r ), c.f. Eq. (43).When using an approximate free energy functional,there is a difference between the the zero-dimensionallimit and Percus’ limit for calculating g ( r ). This is be-cause in the zero-dimensionality limit, in contrast to Per-cus’ method, the test particle at the origin does not inter-act with the fluid via an external potential u ( r ) = v ( r ),that is identical to the pair potential, but rather via anapproximation u ∗ ( r ) to the pair potential, generated bythe approximate density functional. We calculate be-low in Sec. IV C an explicit expression for u ∗ ( r ) in thecase of hard-spheres treated using the RY approxima-tion for the free energy. We will also discuss further therelation between Percus’ test particle limit and the zero-dimensionality limit. III. DYNAMIC TEST PARTICLE LIMITA. Definition
We next extend the static test particle limit and con-sider the dynamical situation which allows us to calculatethe van Hove function G ( r, t ). The key is the followingobservation: Consider a fixed test particle of species ‘ s ’(self) located at the origin; due to Percus we know thatin this situation the density distribution of the remainingparticles ρ ( r ) = ρg ( r ). Now consider releasing the testparticle at time t = 0 and allowing it to move throughthe fluid. When this happens its probability (density)distribution ρ s ( r, t ) changes from a Dirac delta function(at t = 0) to a distribution with a non-zero value for somepoints away from the origin. If we now recall the defini-tion of the function G s ( r, t ) in Eq. (3), we see that G s ( r, t )gives the probability that a particle initially located atthe origin has moved a distance r away from the origin af-ter time t . Therefore, in this situation, ρ s ( r, t ) ≡ G s ( r, t )for all times t ≥
0. Similarly, if we consider how theremaining particles redistribute themselves as the testparticle moves away from the origin, we see from Eq. (3)that the probability of finding any one of these parti-cles a distance r from the origin at time t is given by G d ( r, t ). We label the remaining particles as being parti-cles of species ‘ d ’ (distinct) and having the density profile ρ d ( r, t ). Thus, as in the static test particle case, we mayconnect the two parts of the van Hove function with thedensity profiles of a two-component system of species s and d : G s ( r, t ) ≡ ρ s ( r, t ) ,G d ( r, t ) ≡ ρ d ( r, t ) , (46) where species s is composed of only one particle, the testparticle, and R d r ρ s ( r, t ) = 1, so that Eq. (7) is satisfied.We may therefore formally set the pair potential for inter-actions between species s particles v ss ( r ) = 0. The den-sity profile for species d must satisfy the normalizationconstraint (8), and the self-distinct and distinct-distinctpair potentials must be identical, v sd ( r ) = v dd ( r ) = v ( r ).This is equivalent to modelling a one-component system,but treating one particle separately from the rest. Recallthat at time t = 0 we know the test particle’s positionexactly from Eq. (5) and combining this with Eqs. (6),(42) and (46) we obtain G s ( r, t = 0) ≡ ρ s ( r, t = 0) = δ ( r ) ,G d ( r, t = 0) ≡ ρ d ( r, t = 0) = ρg ( r ) . (47)The connections made in Eq. (46) between the self anddistinct parts of the van Hove function and the densityprofiles ρ s ( r, t ) and ρ d ( r, t ) in the dynamical test particlelimit described above are conceptually important. How-ever, we have merely shifted the problem of how to de-termine G ( r, t ) onto the problem of how to determine thetime evolution of the two coupled density profiles ρ s ( r, t )and ρ d ( r, t ). The solution that we use in this paper is touse DDFT, i.e. equations (29) are integrated forward intime with Eqs. (47) providing the initial time, t = 0, den-sity profiles. The resulting time series of density profilesgives the self and distinct parts of the van Hove functionthrough Eq. (46). Henceforth, we refer to this as the‘dynamical test particle’ theory. B. Approximate solution
Before proceeding to the results that we obtain fromfollowing the calculation scheme described above, it isworth examining an approximate analytical solution thatmay be obtained as follows: From Eqs. (29), (30) and (34)we may write the DDFT equations for the two densityprofiles ρ s ( r, t ) and ρ d ( r, t ) as: ∂ρ i ( r , t ) ∂t = D ∇ ρ i ( r , t ) + Γ ∇ · h ρ i ( r , t ) ∇ c (1) i ( r , t ) i , (48)where i = s, d and the diffusion coefficient D = k B T Γ.If we set the second term on the right hand side to zeroand we set D = D l then we obtain Eq. (15) for ρ s ( r, t ) = G s ( r, t ) and thus the solution to the DDFT for the selfpart of the van Hove function in this limit is the Gaussianform in Eq. (17). Similarly, for species d , when we assumethat the second term on the right hand side of Eq. (48)can be neglected, then we obtain: ∂ρ d ( r , t ) ∂t = D ∇ ρ d ( r , t ) . (49)If we now assume the Vineyard approximation ρ d ( r, t ) = Z d r ′ g ( r ′ ) ρ s ( | r − r ′ | , t ) (50)[c.f. Eq. (19)], then after Fourier transforming Eq. (49),together with Eq. (50) we obtainˆ g ( k ) ∂ ˆ ρ d ( k, t ) ∂t = − k D ˆ g ( k )ˆ ρ d ( k, t ) , (51)where ˆ g ( k ) is the Fourier transform of the radial distri-bution function g ( r ) and ˆ ρ d ( k, t ) is the Fourier transformof ρ d ( r, t ). Dividing both sides of Eq. (51) by ˆ g ( k ) andthen taking the inverse Fourier transform we obtain Eq.(15) for ρ s ( r, t ) = G s ( r, t ). Thus the Gaussian form inEq. (17) together with the Vineyard approximation (50)for ρ d ( r, t ), together form a self consistent solution to theDDFT equations in the dynamical test particle limit, inthe case where we can neglect the contribution from thesecond term on the right hand side of Eq. (48). Thisterm is zero in the ideal-gas limit when the excess contri-bution to the free energy F ex = 0, in Eq. (30), or when ρ →
0. However, we find below for hard spheres thatthis approximation is reliable well beyond the ideal gasregime, which suggests that in the test particle limit c (1) s and c (1) d in Eq. (48) must both be slowly varying (almostconstant) functions, so that their gradients are small. IV. MODEL FOR HARD SPHERESA. Simulation method
In order to provide benchmark results, we calculate thevan Hove function by integrating the equations of motion(20) using standard BD computer simulations . In orderto apply the algorithm we model the hard spheres witha steep continuous pair potential: βv ( r ) = (cid:26) ( σ/r ) − r < σ, δt = 1 × − τ B ; recall that theBrownian time τ B = σ /D , where D = Γ k B T is theStokes-Einstein diffusion coefficient. The random forces ζ i in Eq. (20) mimic the interaction between particlesand solvent, and are sampled from a Gaussian distribu-tion with zero mean and variance 2 Dδt . The simula-tions are carried out using N =1728 particles at densities ρσ = N ( σ/L ) = 0 . , . , . , .
8, and 1 in a cubic boxof volume L .After an equilibration time of 50 τ B , we sampled thedistribution functions G s ( r, t ) and G d ( r, t ) at the times t/τ B =0.01, 0.1, and 1. The distribution functions areaveraged over all possible time intervals t/τ B of a singlesimulation run. The total simulated times are 2 τ B , 20 τ B ,and 200 τ B for the short, medium and long time intervals,respectively. The scaled intermediate scattering functionis calculated from the density autocorrelation functionin Fourier space, φ ( k, t ) = h n k ( t ) n − k (0) i / h n k (0) n − k (0) i ,where n k ( t ) = P Ni =1 exp( − i k · r i ( t )) are the Fourier com-ponents of the local number density. B. The excess free energy functional
In order to implement the dynamical test particle limitwe must (as is almost always the case in density func-tional theory calculations) select an approximation forthe excess part of the Helmholtz free energy functional,Eq. (30). We use the RY functional , which is obtainedfrom the two-component generalisation of Eqs. (35) and(36) by neglecting all terms of order O (∆ ρ ) and higher.We obtain: F ex [ ρ s , ρ d ] = V f ex ( ρ )+ f ′ ex ( ρ ) (cid:26)Z d r ( ρ d ( r ) − ρ ) + Z d r ρ s ( r ) (cid:27) − Z d r Z d r ′ c ( | r − r ′ | ) (cid:26) ( ρ d ( r ) − ρ )( ρ d ( r ′ ) − ρ )+ ( ρ d ( r ) − ρ ) ρ s ( r ′ ) + ρ s ( r )( ρ d ( r ′ ) − ρ ) (cid:27) (53)where f ex ( ρ ) is the bulk excess free energy per unit vol-ume, V is the volume of the system, f ′ ex = ∂f ex /∂ρ and c ( r ) is the bulk pair direct correlation function of thehard sphere fluid with bulk density ρ . We use f ex and c ( r ) as given by the Percus-Yevick approximation .Eq. (53) is perhaps the simplest functional that onecould use to calculate hard sphere fluid density profiles.Our reasons for using this functional are: (i) The struc-ture of the functional is relatively simple (and as a con-sequence is widely used within liquid state approaches).(ii) Within the RY functional it is straightforward to ne-glect the species s intra-species interactions which is nec-essary to ensure that ρ s ( r, t ) represents a single particle.(iii) Finally, the RY functional was the first functional tocorrectly reproduce freezing phenomena in hard spheres. C. Static structure of the fluid
In Fig. 1 we display the radial distribution function andstatic structure factor for a bulk fluid of hard spheres forthe densities ρσ = 0 .
2, 0.4, 0.6, 0.8 and 1. We showresults obtained from BD simulations, together with theresults from the static test particle limit using the RYapproximation for the excess free energy (53). One canobserve that for low and intermediate densities ρσ . . g ( r = σ + ). This in turn leads to the discrepan-cies in the static structure factor at high densities in Fig.1(d); S ( k ) is obtained by Fourier transforming the datain (b). The overall conclusion to be drawn from Fig. 1 isthat the test particle method combined with the RY func-tional provides a reliable description of the fluid structurefor low and intermediate densities, but at higher densities ρσ > .
6, the theory is only qualitatively correct. g ( r ) r / s (a) (b)(c) (d)BD rs =0.2 rs =0.4 rs =0.6 rs =0.8 rs =1 0 1 2 3 4 5 6 7 0 1 2 3 4 5 r / s DFT rs =0.2 rs =0.4 rs =0.6 rs =0.8 rs =1 0 0.5 1 1.5 2 2.5 3 3.5 0 5 10 15 20 25 30 S ( k ) k s k s FIG. 1: (a) and (b) display radial distribution functions, and(c) and (d) static structure factors for a bulk fluid of hardspheres at densities ρσ = 0 .
2, 0.4, 0.6, 0.8 and 1. Parts (a)and (c) are obtained from BD simulations and (b) and (d) viaPercus’ static test particle limit using the RY functional.
We return now to the discussion of the somewhat sub-tle issues concerning the relation between Percus’ testparticle limit and the zero-dimensionality limit. The re-sults from these two calculations are not the same whenone uses an approximate expression for the free energy,such as the RY functional (53). Combining Eqs. (53) and(33), we obtain the following: ρ d ( r ) = Λ − exp (cid:20) βµ d − βf ′ ex ( ρ )+ Z d r ′ c ( | r − r ′ | ) ρ s ( r ′ )+ Z d r ′ c ( | r − r ′ | )( ρ d ( r ′ ) − ρ ) (cid:21) . (54)Now, recall that a test-particle calculation involves fixingone of the particles at the origin, treating it as an exter-nal potential, and then determining the density profile ofthe fluid (species d ) under the influence of this externalpotential. Doing this, using the RY approximation for F ex [ ρ s , ρ d ], Eq. (53), we obtain: ρ d ( r ) = Λ − exp (cid:20) βµ d − βf ′ ex ( ρ ) − βu ( r )+ Z d r ′ c ( | r − r ′ | )( ρ d ( r ′ ) − ρ ) (cid:21) . (55)Comparing equations (54) and (55), we see that Eq. (54)is merely Eq. (55) with the external potential βu ( r ) = βv ( r ) replaced by the effective potential: βu ∗ ( r ) = − Z d r ′ c ( | r − r ′ | ) ρ s ( r ′ ) . (56)In the w = 0 limit, when ρ s ( r ) = δ ( r ), we then obtain: βu ∗ ( r ) = − c ( r ) . (57)One consequence of this random phase-like approxima-tion is that the core condition is violated. The degree towhich the core condition is violated could be used as anindicator towards the reliability of any approximate freeenergy functional.In the remainder of this paper we will display resultsand distribution functions that are derived from Percus’test particle results used as initial condition, though wewill draw attention to results from the zero dimensional-ity route where appropriate. V. RESULTSA. Dynamic approaches
In Fig. 2 we display the two parts of the van Hovefunctions, G s ( r, t ) and G d ( r, t ), measured in the BD sim-ulations for fluid densities ρσ = 0 .
2, 0.4, 0.6, 0.8 and 1.The different curves correspond to the times t/τ B = 0 . G s ( r, t ) appears to have a near-Gaussian formfor all times and densities. For short times G d ( r, t ) ex-hibits a correlation hole for r < σ and it is highly struc-tured for (larger) r > σ . At later times the structurein G d ( r, t ) diminishes and the correlation hole becomes‘filled in’. Recall that Eq. (10) defines the long time limit.Increasing the density beyond ρσ = 1, we find that thesimulated system crystallises onto a regular lattice andthere is no evidence of glass forming behaviour.In Fig. 3 we display the one-body density profiles, ρ s ( r, t ) and ρ d ( r, t ), from the DDFT dynamical test par-ticle method for bulk fluid densities ρσ = 0 .
2, 0.4, 0.6and 0.8. As initial condition, ρ d ( r, t = 0) = g ( r ), wehave used Percus’ test particle method for calculating g ( r ), as shown in Fig. 1. The results in Fig. 3 corre-spond to the same times as the BD curves displayed inFig. 2, namely t/τ B = 0 .
01, 0.1, and 1. Comparing theBD results in Fig. 2 with the DDFT results in Fig. 3 weobserve that for densities ρσ = 0 .
2, 0.4, and 0.6, there isgood qualitative agreement between the simulation andDDFT results. The ρ d ( r, t ) results show a similar amount0 -4 -2 (a)(b)(c)(d)(e) G s ( r , t ) s G d ( r , t ) s t / t B =0.01 t / t B =0.1 t / t B =1 0.00.10.20.310 -4 -2 -4 -2 -4 -2 -4 -2 r / s r / s FIG. 2: The self and distinct parts of the van Hove distribu-tion function, G s ( r, t ) and G d ( r, t ), for a hard sphere fluid,measured in BD simulations at densities: (a) ρσ = 0 . ρσ = 0 . ρσ = 0 . ρσ = 0 . ρσ = 1. The resultsare plotted for times t/τ B =0.01 (solid line), 0.1 (dashed line)and 1 (dotted line). In the semi-logarithmic scale of G s ( r, t )versus r a Gaussian appears as a parabola. The G d ( r, t ) re-sults are shown on a linear scale. of structure as the G d ( r, t ) results, and ρ s ( r, t ) has a verysimilar magnitude and range as G s ( r, t ), although for ρσ = 0 . ρ s ( r, t ) shows some departures from the al-most Gaussian shape observed in the simulation results,particularly at t/τ B = 1. For ρσ = 0 . t/τ B ∼ . ρ s ( r, t ) and ρ d ( r, t )cease to change with time and that the system becomes‘arrested’. One could interpret this state as the taggedparticle remaining localised within the cage formed bythe neighbouring fluid particles. We discuss the signifi- -4 -2 (a)(b)(c)(d) r s ( r , t ) s r d ( r , t ) s t / t B =0.01 t / t B =0.1 t / t B =1 0.00.10.20.310 -4 -2 -4 -2 -4 -2 r / s r / s FIG. 3: The ‘ s ’ and ‘ d ’ density profiles, ρ s ( r, t ) and ρ d ( r, t ),obtained from the dynamical test particle theory, for densi-ties: (a) ρσ = 0 . ρσ = 0 . ρσ = 0 . ρσ = 0 . t/τ B =0.01 (solid line), 0.1(dashed line) and 1 (dotted line). In (d), after a short time,the system reaches an ‘arrested state’, where the density pro-files no longer evolve in time and the width of ρ s ( r, t → ∞ )remains finite. cance of this phenomenon in Sec. VI.In Fig. 4 we show the van Hove functions calculatedusing the Vineyard approximation, Eqs. (17) and (19)with D l = D = k B T Γ, for fluid densities ρσ = 0 . t/τ B = 0 .
01, 0.1, and 1.As in the DDFT we have used g ( r ) calculated using theRY functional and Percus’ test particle method, althoughone could use g ( r ) obtained from any reasonable method,including g ( r ) measured in the BD simulations. Compar-ing the Vineyard results to the BD simulation results inFig. 2, we find that there is reasonably good agreementbetween the two. The form of G s ( r, t ) is fixed to beGaussian, so there is good agreement with G s ( r, t ) fromthe BD simulations, though it is clear that the width of G s ( r, t ) increases more rapidly in the Vineyard approxi-mation. For densities ρσ ≤ . G d ( r, t ) for r > σ in the Vineyardapproximation as in the simulation results. However, for ρσ = 1 the Vineyard approximation does not exhibit the1 -4 -2 (a)(b)(c)(d)(e) G s ( r , t ) s G d ( r , t ) s t / t B =0.01 t / t B =0.1 t / t B =1 0.00.10.20.310 -4 -2 -4 -2 -4 -2 -4 -2 r / s r / s FIG. 4: The self and distinct parts of the van Hove distribu-tion function, G s ( r, t ) and G d ( r, t ), calculated using the Vine-yard approximation for densities: (a) ρσ = 0 . ρσ = 0 . ρσ = 0 . ρσ = 0 . ρσ = 1. The results areplotted for times t/τ B =0.01 (solid line), 0.1 (dashed line) and1 (dotted line). same degree of structure that is present in the simulationresults.In Fig. 5 we compare the width, w ( t ), of the self partof the van Hove function, G s ( r, t ), obtained from (i) BDsimulation results, (ii) dynamical test particle limit, and(iii) the Vineyard approximation. In the Vineyard ap-proximation the time dependence of w ( t ) is defined byEq. (16) and does not depend on density, so there is onlyone master curve. This is because in the same way as inthe dynamical test particle limit, we set D l = D , where D = k B T Γ is the short time diffusion coefficient, whichis strictly only equal to the long time self diffusion coef-ficient D l in the limit ρ →
0. Since w ( t ) ∝ √ t , on the -1 -3 -2 -1 w / s t / t B Increasing r -1 -3 -2 -1 w / s t / t B Increasing r -1 -3 -2 -1 w / s t / t B Increasing r FIG. 5: The width w of the self part of the van Hove function,defined in Eq. (11), measured in the BD simulations for ρσ =0 . × ), 0.6 ( ∗ ), 0.8 ( (cid:3) ), and 1 ( ◦ ), and from thedynamical test particle theory, for ρσ = 0 . double logarithmic scale in Fig. 5 this is represented bya straight line with gradient 1 /
2. We find that the sim-ulation results are also approximately linear in this rep-resentation for all densities considered, but that there isa slowing down effect as density is increased, due to thefact that D l decreases as the fluid density is increasedand is no longer equal to D . For ρσ = 0 .
2, the simula-tion results are close to the Vineyard result. As the bulkdensity is increased, the BD results move away from thisline.The dynamical test particle results for w ( t ) in Fig. 5exhibit a much stronger dependence on density. At lowdensities the curves are similar to the Vineyard and theBD results, but as the density is increased the w ( t ) curvesshow a slowing down, and then (unphysical) speedingup of the dynamics, unlike that seen in the BD results.This slowing down is greatly exaggerated so that theDDFT curve for ρσ = 0 . ρσ = 0 .
8. Furthermore, the DDFT curves for ρσ = 0 . ρσ = 0 . w ( t ) curve shows extremely exagger-ated slowing down and speeding up. We believe that theunphysical speeding up for t/τ B & is due to the factthat the DDFT incorrectly sets the long-time diffusioncoefficient D l equal to the short time diffusion coefficient D , so that as the particle escapes the cage of neighboringparticles, it is forced to “catch-up” to give the incorrectlong time behavior. Note that from the Smoluchowskiequation (22) it can be shown that w ( t ) must be sub-diffusive for intermediate times, and that the long timediffusion coefficient must be smaller than the short timeone, a feature which is well established in Brownian dy-2 F s ( k , t ) F d ( k , t ) t / t B =0.01 t / t B =0.1 t / t B =1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 k s k s FIG. 6: Intermediate scattering function F ( k, t ) as a functionof the scaled wave vector kσ , obtained by a spatial Fouriertransform, Eq. (12), of the BD simulation results for the vanHove functions displayed in Fig. 2. namics simulations and experiments . The curve for ρσ = 0 . w ( t → ∞ ) is fi-nite, as one would infer from the density profiles shown inFig. 3. We postpone a discussion of the possible physicalimplications to Sec. VI.For completeness, we also plot the intermediate scat-tering function F ( k, t ). In Fig. 6 we display results fromBD computer simulations, and in Fig. 7 the results fromthe DDFT. We find that for ρσ = 0 .
2, 0.4 and 0.6, theresults from both approaches exhibit very similar struc-ture. However, since in the DDFT the dynamics becomearrested at ρσ = 0 .
8, so F ( k, t ) becomes arrested aftera very short time, unlike the BD simulations result. In F s ( k , t ) F d ( k , t ) t / t B =0.01 t / t B =0.1 t / t B =1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 k s k s FIG. 7: Intermediate scattering functions F ( k, t ), obtained bya spatial Fourier transform, Eq. (12), of the dynamic test par-ticle density profiles, ρ s ( r, t ) and ρ d ( r, t ), displayed in Fig. 3. Fig. 8 we plot the scaled intermediate scattering function, φ ( k, t ) = F s ( k, t ) F s ( k, t = 0) , (58)for fixed kσ = 2 π , obtained from the Vineyard approxi-mation and compare to the BD simulation results (Fig.8(a)) and the DDFT results (Fig. 8(b)). At the lower den-sities the BD simulation results and the DDFT results areclose to the Vineyard approximation and both show someslowing down with density. At the higher densities theBD results continue to show a steady decay. However, inthe DDFT results for ρσ = 0 . φ ( kσ = 2 π, t ) de-cays in two stages over a much longer time. For ρσ = 0 . φ ( kσ = 2 π, t ) to remain fi-nite in the limit t → ∞ . B. Relating dynamic to static density profiles
One may connect DDFT and equilibrium DFT by find-ing the unique set of effective external potentials { u i ( r ) } that in equilibrium generate the same set of density pro-files as obtained in the dynamic approach at a particular3 -4 -3 -2 -1 f ( k s = p , t ) t / t B (a) Vineyard approx.0.00.20.40.60.81.010 -4 -3 -2 -1 f ( k s = p , t ) t / t B (a) BD rs =0.2 rs =0.4 rs =0.6 rs =0.8 rs =1 0.00.20.40.60.81.010 -4 -3 -2 -1 f ( k s = p , t ) t / t B (b) Vineyard approx.0.00.20.40.60.81.010 -4 -3 -2 -1 f ( k s = p , t ) t / t B (b) DDFT rs =0.2 rs =0.4 rs =0.6 rs =0.7 rs =0.8 FIG. 8: Scaled intermediate scattering function φ ( kσ = 2 π, t )as a function of time t/τ B , calculated from the Vineyard ap-proximation (thick gray line), compared to (a) BD simulationresults, and (b) the dynamical test particle method. Thedynamical test particle results exhibit slowing and arresteddynamics for ρσ = 0 . time t . These potentials represent the net effect of fi-nite time and limited diffusion preventing the fluid fromfinding the structure that minimises the system free en-ergy. For the two-component system considered here, thetwo external potentials βu i ( r, t ) may be recovered, up toan overall time-dependent additive constant βµ i ( r, t ), byrearranging Eq. (33): βu i ( r, t ) − βµ i = c (1) i ( r ; [ ρ s , ρ d ]) − ln[Λ ρ i ( r, t )] , (59)where ρ s ( r, t ) and ρ d ( r, t ) are the solution of the DDFTat time t . In Fig. 9 we plot these external potentials cor-responding to the density profiles from DDFT displayedin Fig. 3. We find that at the lowest densities ρσ = 0 . .
4, the shape of u s ( r, t ) is approximately parabolicfor all times and distances r . As the fluid density is in-creased, u s ( r, t ) departs from the parabolic form. For ρσ = 0 . r , butat small r they become distorted. We find that u d ( r, t )does not vary significantly with density. At short timesit is dominated by strong repulsion within the hard-core -4 -2 (a)(b)(c)(d) b u s ( r , t ) b u d ( r , t ) t / t B =0.01 t / t B =0.1 t / t B =1 r -4 -2 -4 -2 -4 -2 r / s -1 -1.00.01.02.03.04.0-1.00.01.02.03.04.0-1.00.01.02.03.04.0-1.00.01.02.03.04.0 r / s FIG. 9: External potentials, u s ( r, t ) and u d ( r, t ), requiredto generate equilibrium density profiles, ρ s ( r, t ) and ρ d ( r, t ),identical to those obtained from the dynamical test particleapproach – see Fig. 3. The external potentials also includean overall time-dependent additive constant which is not in-dicated. Note the logarithmic x -axis for u s ( r, t ). At low den-sities u s ( r, t ) is approximately parabolic, as indicated by thestraight line (dashed-dotted). At high densities u s ( r, t ) is dis-torted at small r , but still parabolic at large r . diameter, r < σ . Recall that in order to calculate g ( r ),which corresponds to t = 0, we chose to use Percus’ testparticle method, and hence have introduced an externalpotential equal to the hard sphere potential. The strongrepulsion found for short times is a remnant of this exter-nal potential. As t increases, the strength of this repul-sion decreases and becomes almost zero for t/τ B = 1 . C. Corresponding equilibrium approach
Having established the form of the external poten-tials necessary to create equilibrium fluid density profilesequal to the profiles calculated using the dynamical testparticle method, we now seek a simple approximationfor this set of external potentials, to allow us to easilycalculate equilibrium density profiles that mimic the dy-4namic profiles. In other words, we seek to determine thefull van Hove function when G s ( r, t ) has a given width w , without calculating the entire preceding time series ofprofiles. In doing this we lose time t as a function ar-gument and instead we must ‘label’ the density profileswith w . In what follows, we will disregard the associatedproblem of relating w to time t .We parametrise the external potentials using a sim-ple functional form. Firstly, we assume that u s ( r, w ) isparabolic for all widths: βu s ( r, λ ) = λr , (60)where λ is the strength of the confining potential, and w is now an unknown function of λ . For the externalpotential that acts on species d we consider two options.The first is to assume that u d ( r, w ) = 0. In this case, it ispossible to solve exactly for the equilibrium distributionfunctions and the free energy, as outlined in Appendix A.We find that the species s density profile is a Gaussian, ρ s ( r, λ ) = exp( − βλr )( π/βλ ) / , (61)where the dependence of the width w on λ is, w ( λ ) = (2 λ/ − / , (62)and the d profile is given by a convolution of the radialdistribution function g ( r ), together with the Gaussianprofile ρ s ( r ), and multiplied by ρ : ρ d ( r ) = ρ Z d r ′ ρ s ( r ′ ) g ( | r − r ′ | ) . (63)These distribution functions are identical to those fromthe (dynamic) Vineyard approach.The second approach that we consider is to calculatethe density profiles without defining the external poten-tial u d ( r, w ) at the outset of the calculation. Instead, wedetermine this potential self consistently ‘on-the-fly’ aspart of our iterative numerical solution routine, based onthe following considerations: Firstly, recall the normali-sation constraints on the van Hove function in Eqs. (7)and (8). In order to satisfy the normalisation constraint(7) on the density profile for the single tagged s particle,we introduce a Lagrange multiplier µ s . One may alsoconsider λ to be a Lagrange multiplier that enforces thewidth constraint (11) on the profile ρ s ( r ). In our calcula-tions the value of µ s is determined on-the-fly by enforcingEq. (7) at each step of our iterative routine. However,we are not able to do the same for the density profile ofthe remaining d particles, because we also must have ρ d ( r, w ) → ρ, as r → ∞ . (64)Multiplying ρ d ( r ) by a single factor breaks this condition,so we are not able to simply enforce (8) at each stepof our iterative routine in the same way as we do for ρ s ( r ). The condition in Eq. (64) implies that we require -4 -2 (a)(b)(c)(d) r s ( r , w ) s r d ( r , w ) s -4 -2 -4 -2 -4 -2 r / s r / s FIG. 10: The density profiles ρ s ( r, w ) and ρ d ( r, w ), calculatedusing the equilibrium DFT. The curves are calculated at thedensities: (a) ρσ = 0 .
2, (b) ρσ = 0 .
4, (c) ρσ = 0 . ρσ = 0 .
8. The curves are chosen so that the widths w of ρ s ( r, w ) correspond to the same widths of G s ( r, t ) at times t/τ B =0.01 (solid line), 0.1 (dashed line) and 1 (dotted line)displayed in Fig. 2. an a priori unknown inhomogeneous external potential, u d ( r, w ), with the property that u d ( r, w ) → r → ∞ .This may be achieved by scaling the quantity ρ d ( r ) − ρ (instead of scaling ρ d ( r ) itself) at each step, so that ρ d ( r ) satisfies both (8) and (64). Once convergence of thenumerical procedure is achieved, one may then inspectthe effective external potential u d ( r, w ) by substitutingthe resulting density profiles into Eq. (59).The density profiles calculated using this equilibriummethod are shown in Fig. 10 where we plot ρ d ( r, w ) and ρ s ( r, w ) having widths identical to those of the van Hovefunctions from simulation, displayed in Fig. 2. Note thatwe consider only the densities, ρσ = 0 .
2, 0.4, 0.6 and 0.8.These equilibrium profiles have been calculated using anormalisation constant taken from the approximation for g ( r ) calculated using Percus’ test particle method. Wefind that there is reasonable qualitative agreement be-tween the equilibrium DFT density profiles displayed inFig. 10 and the BD simulation results displayed in Fig. 2.5 -10.0-5.00.05.010.00.0 0.5 1.0 1.5 2.0 b F ( w ) w / s Exact result-10.0-5.00.05.010.00.0 0.5 1.0 1.5 2.0 b F ( w ) w / s DFT rs =0.2 rs =0.4 rs =0.6 rs =0.7 rs =0.8-10.0-5.00.05.010.00.0 0.5 1.0 1.5 2.0 b F ( w ) w / s DDFT rs =0.2 rs =0.4 rs =0.6 rs =0.7 rs =0.8 FIG. 11: The free energy landscape F ( w ) plotted as a func-tion of w , the width of ρ s ( r, w ), calculated using both theDDFT and equilibrium DFT approaches, compared to the ex-act result, Eq. (65). For ρσ = 0 . F ( w ). However, the profiles predict too much infilling in the re-gion close to the origin r < σ , particularly at higherdensities, which in turn results in an underestimate inthe structure of the profiles at larger r > σ . This erroroccurs for the reasons discussed in the previous subsec-tion IV C; i.e. that the RY functional does not exert astrong enough interaction from the test particle onto therest of the fluid.For ρσ = 0 .
8, shown in Fig. 10(d), we are able tocalculate density profiles for all values of w , even thoughin the DDFT the profiles became ‘trapped’ at small w .The most striking aspect of these density profiles is thatfor intermediate values of w we see that ρ s ( r ) exhibitsa plateau and a long tail. These features were not ob-served in the BD simulation results. However, similarfeatures are present in G s ( r, t ) at intermediate times forcolloidal spheres at densities close to the glass transition,where they are interpreted as the signature of dynamicalheterogeneity in the system . D. Free energy landscape
Since we are able to convert the dynamic density pro-files into their equilibrium equivalents via a set of externalpotentials, c.f. Eq. (59), we are also able to calculate theequilibrium Helmholtz free energy for this correspond-ing equilibrium situation. Although this free energy isstrictly an equilibrium quantity, since it underlies thetime evolution of our dynamic approach we believe thatit plays a relevant role. Therefore, by substituting thedensity profiles calculated using the DDFT into the freeenergy functional, Eqs. (30) and (53), we are able to mapout a ‘free energy landscape’ as a function of t or w . Fig.11 plots this free energy landscape, F ( w ), for densities ρσ = 0 .
2, 0.4, 0.6, 0.7 and 0.8. For ρσ = 0 .
2, 0.4, -10.0-5.00.05.010.00.0 0.5 1.0 1.5 2.0 b F ( w ) w / s Exact resultVineyard rs =0.2 rs =0.4 rs =0.6 rs =0.7 rs =0.8 FIG. 12: The free energy landscape F ( w ) plotted as a functionof the width w/σ , calculated by substituting the density pro-files from the exact equilibrium solution [Eqs. (61) and (63)]into the RY functional. Here the deviation from the exactfree energy and the emergence of the minimum are entirelydue to the RY functional. -2 -1 bl w / s Exact result10 -2 -1 bl w / s DFT rs =0.2 rs =0.4 rs =0.6 rs =0.7 rs =0.8 10 t* / t B =2 t* / t B =200 FIG. 13: The (scaled) strength of the parabolic potential, βλ ,versus the resulting width, w/σ , of ρ s ( r ). Results are calcu-lated using the equilibrium DFT and compared to the exactresult. For the bulk fluid densities ρσ = 0 .
2, 0.4, 0.6 and 0.7, λ decreases monotonically with w . However, for ρσ = 0 . t ∗ is the simulation run time. and 0 . F ( w ) decreases monotonically with w . This decrease is initially steep and then the gradientbegins to reduce as w increases. For ρσ = 0 .
7, we findthat after the initial steep descent the landscape devel-ops an almost constant plateau, but there is still a verysmall negative gradient. For ρσ = 0 . u d ( r, w ) = 0, one hasan exact solution (see Appendix A) for the free energy6landscape as a function of the width, F ( w ) = F id − ln( Z ′ N ) −
32 ln (cid:18) πw β (cid:19) , (65)where F id is the ideal Helmholtz free energy, and Z ′ N isan irrelevant constant representing the partition functionof the fluid when the test particle is located at the ori-gin. We plot Eq. (65) alongside the landscapes from theDDFT approach in Fig. 11. We find that this curve islocated close to the DDFT landscape for ρσ = 0 . ρσ = 0 .
2, 0.4 and 0.6,we find that there is good agreement between the DDFTand equilibrium DFT approaches. For ρσ = 0 . w , but around w/σ = 0 . F ( w ) in the equi-librium DFT results that is not present in the DDFT re-sults. For ρσ = 0 . w .However, since we can calculate the density profiles usingour equilibrium approach for all widths we can thereforecalculate F ( w ) for all w . We find a free energy landscapewith two disconnected branches. Therefore, for a rangeof values of λ we find two solutions to the Euler-Lagrangeequations with different widths. Whether this is an in-dication of dynamic heterogeneity or an artefact of thefunctional is an interesting question.If we take the density profiles calculated using the ex-act route [Eqs. (61) and (63)] and substitute these intothe RY functional, we find that the free energy curves,shown in Fig. 12, do not follow the exact result (65), butare very similar to those from the DDFT and DFT ap-proaches and even exhibit minima for ρσ = 0 . λ , againstthe width, w , and compare it to the results from theequilibrium DFT approach. At the lowest densities theequilibrium DFT approach closely follows the exact re-sult, but as the density is increased the width decreasesfor a given λ . For ρσ = 0 . ρσ = 0 . t ∗ is the simulation run time. The results for thelonger simulation run time agree well with the exact re-sult, but results over the shorter time underestimate thewidth, particularly at the smaller values of λ . A similareffect may exist in the DFT approach, where the non-ergodicity arises from not including the states where theparticle is far from the origin. Recall that formally thedensity profile from DFT is an average over all possiblestates of the system. Using an approximate functionalsome of these states may be neglected . VI. CONCLUDING REMARKS
On the basis of the theoretical and simulation resultsthat we presented in this paper for the dynamics of thebulk hard sphere fluid, we conclude that the dynamicaltest particle limit, combined with DDFT, provides a re-liable method for calculating the van Hove (and otherrelated) dynamical pair correlation functions at low andat intermediate densities ρσ . .
6. In the previouspublication we have shown that the theory may beapplied in a fairly straightforward manner also to inho-mogeneous situations, hence we conclude that the dy-namic test particle theory may indeed be used to calcu-late the van Hove function for fluids at interfaces and un-der confinement. Furthermore we have shown that, quitesurprisingly, at low and intermediate densities the verysimple Vineyard approximation (reviewed in Sec. II C) isactually quantitatively fairly good. This approximationonly requires as input the radial distribution function g ( r ) and the diffusivity D l and therefore provides a veryuseful quick approach for obtaining an approximation forthe van Hove function for colloidal fluids.The possible conclusions about the performance andeven about the qualitative status of the predictions ofthe theory at higher densities are far more intricatethough. In this regime, the theory in its current formis clearly quantitatively unreliable – compare for exam-ple, the results from BD simulations in Fig. 2(d) tothose from the dynamic test particle theory shown inFig. 3(d) for ρσ = 0 .
8, where the theory predicts thatthe system jams, whereas in simulations the system isan ergodic fluid. Moreover, at even higher densities, themonodisperse hard sphere system crystallises in simula-tions rather than undergoing a glass transition; polydis-persity would be required to suppress freezing . We be-lieve that the behaviour of the theory at higher densitiescan primarily be ascribed to our choice of approxima-tion for the free energy functional – see Sec. IV B for thereasons for using the RY functional (53) in the presentstudy. We believe that if we had used a more accu-rate functional, such as Rosenfeld’s fundamental measuretheory , quantitatively more accurate results couldbe obtained at these higher densities. Nevertheless, thetheory yields a clear prediction that there is a dynamic(glass) transition, whereby the tagged (self) particle be-7comes trapped in the cage formed by the surroundingparticles. We can offer three possible interpretations ofthis result.Firstly, one may conclude that this glass transitionstems simply from the use of the approximate RY func-tional and since the theory predicts the glass transitionto be at a density value that is well below where theglass transition is believed to occur, our results at higherdensities should be disregarded.An alternative conclusion that one may draw is thatthe theory is correctly describing (some of) the physicsof the glass transition, but that the predictions are onlyqualitative in nature and occur in reality at higher densi-ties and in polydisperse systems. Support for this pointof view comes in particular from results such as thosein Fig. 8, where we display results for φ ( kσ = 2 π, t ).For ρσ = 0 .
7, which is near to where the theory pre-dicts the glass transition to be, we find a plateau in φ ( k, t ) and the clear presence of two-stage relaxation,which indeed has been observed in hard sphere colloidalsuspensions . Hence our results are (qualitatively) simi-lar to those from MCT . Further support for the aboveinterpretation comes from the results in Figs. 11 and 12for the behaviour of the free energy landscape that un-derlies the DDFT. The appearance of a minimum in thefree energy as a function of displacement, correspond-ing to a particle becoming trapped in the cage formedby its neighbours, is one central prediction of the the-ory by Schweizer and Saltzmann , who combined ele-ments of MCT, DFT, and activated rate theory, in orderto describe localization and transport in glassy fluids.Furthermore, given some of the work in the literaturebased on DFT to study the glass transition, our predic-tion that particles become localized should not come asa surprise: Wolynes and coworkers developed a suc-cessful model of hard sphere vitrification, which is similarto the DFT treatment of crystallization . Using a ran-dom close-packed, non-periodic lattice they found a fluid-glass transition where the fluid “crystallizes” onto thislattice. The success of this method, along with its abilityto model the freezing transition (onto a regular lattice),has provoked a number of further developments .Other approaches have investigated dense Browniansystems through modelling via certain stochastic differen-tial equations and found that the system exhibits glassybehaviour. Thus, overall our results seem to be qualita-tively consistent with other DFT based theories and withMCT for the glass transition. Nevertheless, the densitywhere the glass transition occurs, as predicted by the the-ory in its present form, is far too low. Furthermore, itcould be the case that the similarity between our resultsand those from the MCT are somehow a mathematical(rather than physical) coincidence, since an essential fea-ture of MCT is the presence of memory in the dynamicalequations. This important feature is absent from thepresent DDFT.The third possible conclusion that one may draw con-cerns the question whether the minimum in the free en- ergy and the localization of the tagged self particle aremerely a signature of freezing in the theory. The RY func-tional is well-known to predict the freezing transition tooccur at a density below that where it occurs in reality.It could simply be the case that this functional overlyfavours freezing, so that when it is applied in the waywe use it here, where we constrain all density profiles, ρ s ( r ) and ρ d ( r ), to be spherically symmetric, a signatureof freezing shows up as the tagged self particle becominglocalized.Some merit can be found in all of the arguments out-lined above and we find ourselves unable to judge whichone(s) are correct. Indeed further work is required toprovide a clear assessment of these issues. In particu-lar, the dynamical test particle theory should be imple-mented with a more sophisticated approximation for thefree energy functional than we have used here.As we have shown, our approach is based on integrat-ing the Smoluchowski equation (22) over all except oneof the position coordinates, in order to derive an equa-tion for the one-body density distribution. An alterna-tive approach is to integrate over all but two of the po-sition coordinates, in order to obtain (23) an equationfor the two-body distribution function ρ (2) ( r , r , t ) = N ( N − R R d r ... R d r N P ( r N , t ). The resulting dy-namical equation depends on the three body distributionfunction ρ (3) ( r , r , r , t ). On making a suitable closureapproximation, this provides a different starting point forstudying the pair correlations in a colloidal fluid – see e.g.Refs. 52,53 and references therein, which also considerthe effect of the hydrodynamic interactions between thecolloids. Developing the theory for the dynamical paircorrelation functions in this way is very natural. How-ever, we believe that the strength of our method, wherewe use the dynamical test particle approach allowing usto work at the one-body level, is that we are able to useDFT to close our equations and therefore we are able todescribe the fluid spatial correlations very accurately.Finally, we mention other possible directions for de-veloping the theory in the future. One important aspectin the dynamics of colloidal dispersions, that we haveentirely neglected here, are the hydrodynamic interac-tions between the particles. Rex and L¨owen haveshown how to include the hydrodynamic interactions ina DDFT treatment and so it would be worthwhile to usetheir DDFT formulation together with the present dy-namical test particle limit, in order to calculate the vanHove function under the influence of hydrodynamic in-teractions.A further aspect of our work, that offers possible exten-sions of the theory, concerns the question how to modelthe diffusivity of the tagged particle in a better way: Inthe dynamical test particle calculation one could replacethe (constant) diffusion coefficient in Eq. (48) with a dif-fusion coefficient that depends on time; i.e. to replace D → D ( t ). In doing this one could ensure that D ( t )takes the correct values at both short and long times.However, doing this still does not treat memory effects8in the dynamics. As MCT demonstrates, memory effectsare key for a system to exhibit the ideal glass transitionscenario . Thus, we believe that including memoryinto our theory would be a crucial step in future work.This could possibly be done along the lines of the in-teresting work of Medina-Noyola and coworkers . Toinclude memory in our theory one could replace Eq. (29)with : ∂ρ i ( r , t ) ∂t = ∇ · Z t dt ′ Z d r ′ Γ( r − r ′ , t − t ′ ) × (cid:20) ρ i ( r ′ , t ′ ) ∇ δF [ { ρ i } ] δρ i ( r ′ , t ′ ) (cid:21) , (66)where the mobility coefficient Γ has been replaced byone that is non-local in time and space. However, thiswould result in a considerable increase in computationalcomplexity as within DFT the correlations in space arealready treated in a complex manner and these wouldneed to be coupled to the correlations in time. Whethersuch non-locality helps to cure some of the deficiencies ofour approach is an open question. Appendix A: Exact Results
We consider a fluid of N particles with positions r i ,momenta p i and mass m in the presence of an arbitraryexternal field that acts only on particle i = 1, u ( r ) = λ r . Assuming that we are in the classical limit, theHamiltonian is given by H N = K + V + U where thecontributions are due to the (classical) kinetic energy,the total inter-particle potential (not necessarily pairwiseadditive), and the external potential, respectively; K = N X i =1 p i m (A1) V = v ( r , .., r N ) (A2) U = λ r . (A3)The canonical partition function , Q N ( V, T ), is given by Q N ( V, T ) = h − N ( N − Z Z d r N d p N exp[ − βH N ( r N , p N )](A4)where h is Planck’s constant, and the ( N − i = 1, the re-maining particles are indistinguishable. The integrationsover momenta in Eq. (A4) can be carried out explicitly,leaving a configuration integral over positional degrees offreedom: Z N = Z d r ... d r N exp( − β ( V + U )) . (A5)Note that for Brownian particles Z N is also the quantitythat characterises the structure of the fluid. Substituting our external potential (A3) into (A5) we obtain Z N = Z d r exp( − βλ r ) Z d r ..N exp( − βV ( r N ))= Z d r exp( − βλ r ) Z d r ′ ..N exp( − βV ( r ′ ..N )) , where in the second step we have made the substitution r ′ i = r i − r , for i = 2 ..N , so that we can do the integra-tions over the positions r ′ ... r ′ N . This gives, Z N = R d r exp( − βλ r ) Z ′ N = ( π/βλ ) / Z ′ N , where Z ′ N is the configuration integral for N particleswhere one particle is located at the origin. The Helmholtzfree energy is then given by, F = − β − ln( Q N ( V, T )) = − β − ln( Q id N Z N ( V,T ) V N ) where V N is the volume occupiedby the particles, which yields βF = F id − ln( Z ′ N /V N ) −
32 ln (cid:18) πβλ (cid:19) . (A6)Therefore, the Helmholtz free energy only depends on theconfining potential in a simple way.One can also obtain the one body density profiles. Ingeneral, for a system of N particles, the one body densityprofile, ρ (1) N ( r ) can be obtained from , ρ (1) N ( r ) = N ! Z N ( N − Z d r ( N − exp (cid:2) − β ( V ( r N ) + Φ( r N )) (cid:3) , (A7)where the N ! / ( N − ρ (1)1 ( r ) = exp( − βλ r ) Z N Z d r ..N exp( − βV ( r N )) , = exp( − βλ r )( π/βλ ) / Z ′ N Z d r ′ ..N exp( − βV ( r ′ ..N )) , = exp( − βλ r )( π/βλ ) / , which is a normalised Gaussian. It can be shown thatsince ρ s ( r ) is a Gaussian, then w and λ are simply relatedby w = (2 λ/ − / , (A8)and we can rewrite Eq. (A6) as F = F id − ln( Z ′ N ) −
32 ln (cid:18) πw β (cid:19) . (A9)We now seek the density profile of the remaining parti-9cles: ρ (1)2 ( r ) = ( N − Z N ( N − Z d r exp( − βλ r ) × Z d r ..N exp( − βV ( r N )) , = R d r exp( − βλ r )( π/βλ ) / (A10) × ( N − Z ′ N Z d r ..N exp( − βV ( r N )) . To progress we make use of the formal relationship be-tween g ( r ) and the two body density profile, ρ (2) N ( r , r ),which for a homogeneous fluid can be shown to be , g (2) N ( r − r ) = 1 ρ ρ (2) N ( r − r ) , = 1 ρ N ! Z N ( N − × Z d r ..N exp( − βV ( r N )) , = VρN N ( N − Z ′ N V × Z d r ..N exp( − βV ( r N )) , where we have made the substitutions, ρ = N/V and Z N = Z ′ N V . Cancelling terms and rearranging we get, ρg (2) N ( r − r ) = ( N − Z ′ N Z d r ..N exp( − βV ( r N )) . (A11)Substituting (A11) into (A10) gives, ρ (1)2 ( r ) = R d r exp( − βλ r )( π/βλ ) / ρg (2) N ( r − r )= ρ ( π/βλ ) / Z d r exp( − βλ r ) g (2) N ( r − r ) , which is the normalised Gaussian convolved with ρg ( r ). Acknowledgements
PH thanks the EPSRC for funding under grantEP/E065619/1 and AJA gratefully acknowledges finan-cial support from RCUK. MS and AF thank DFG forsupport via SFB840/A3. ∗ [email protected] † [email protected] J.-P. Hansen and I. R. McDonald,
Theory of Simple Liquids (Academic Press, London, 2006), 3rd ed. L. van Hove, Phys. Rev. , 249 (1954). W. K. Kegel and A. van Blaaderen, Science , 290(2000). J.-P. Hansen and L. Verlet, Phys. Rev. , 151 (1969). P. N. Pusey and W. van Megen, Nature , 340 (1986). W. G¨otze,
Liquid, Freezing and Glass Transition, Proceed-ings of the Les Houches Summer School, ed. J.-P. Hansen,D. Levesque, and J. Zinn-Justin (North-Holland, Amster-dam, 1989). E. R. Weeks and D. A. Weitz, Phys. Rev. Lett. , 095704(2002). L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, D. E.Masri, D. L’Hote, F. Ladieu, and M. Pierno, Science ,1797 (2005). W. van Megen and S. M. Underwood, Phys. Rev. Lett. ,2766 (1993). J. M. Brader, T. Voigtmann, M. E. Cates, and M. Fuchs,Phys. Rev. Lett. , 058301 (2007). J. M. Brader, M. E. Cates, and M. Fuchs, Phys. Rev. Lett. , 138301 (2008). L. Yeomans-Reyna and M. Medina-Noyola, Phys. Rev. E , 3382 (2000). L. Yeomans-Reyna and M. Medina-Noyola, Phys. Rev. E , 066114 (2001). L. Yeomans-Reyna, H. Acu˜na Campa, F. Guevara-Rodr´ıguez, and M. Medina-Noyola, Phys. Rev. E , 021108 (2003). M. A. Ch´avez-Rojo and M. Medina-Noyola, Physica A , 55 (2006). M. Medina-Noyola and P. Ramirez-Gonzalez, J. Phys.:Condens. Matter , 504103 (2009). A. J. Archer, P. Hopkins, and M. Schmidt, Phys. Rev. E , 40501 (2007). M. Bier, R. van Roij, M. Dijkstra, and P. van der Schoot,Phys. Rev. Lett. , 215901 (2008). J. K. Percus, Phys. Rev. Lett. , 462 (1962). R. Evans,
Fundamentals of Inhomogeneous Fluids (NewYork: Dekker, 1992). U. M. B. Marconi and P. Tarazona, J. Chem. Phys. ,8032 (1999). U. M. B. Marconi and P. Tarazona, J. Phys.: Condens.Matter , A413 (2000). A. J. Archer and R. Evans, J. Chem. Phys. , 4246(2004). T. V. Ramakrishnan and M. Yussouff, Phys. Rev. B ,2775 (1979). G. H. Vineyard, Phys. Rev. , 999 (1958). J.-P. Hansen and I. R. McDonald,
Theory of Simple Liquids (Academic Press, 1986), 2nd ed., although, interestingly,the third edition has this comment removed. J. K. G. Dhont,
An Introduction to Dynamics of Colloids (Elsevier, New York, 1996). J. Dzubiella and C. N. Likos, Phys. Rev. , L147 (2003). F. Penna, J. Dzubiella, and P. Tarazona, Phys. Rev. E ,61407 (2003). A. J. Archer, J. Phys.: Condens. Matter , 1405 (2005). A. J. Archer, J. Phys.: Condens. Matter , S3253 (2005). M. Rex, H. L¨owen, and C. N. Likos, Phys. Rev. E ,21404 (2005). M. Rex, C. N. Likos, H. L¨owen, and J. Dzubiella, Mol.Phys. , 527 (2006). C. P. Royall, J. Dzubiella, M. Schmidt, and A. vanBlaaderen, Phys. Rev. Lett. , 188304 (2007). M. Rex, H. H. Wensink, and H. L¨owen, Phys. Rev. E ,21403 (2007). M. Rauscher, A. Dominguez, M. Kr¨uger, and F. Penna, J.Chem. Phys. , 244906 (2007). G. K. L. Chan and R. Finken, Phys. Rev. Lett. , 183001(2005). H. J. Kreuzer,
Nonequilibrium Thermodynamics and itsStatistical Foundations (Oxford University Press, Oxford,1981). P. M. Chaikin and T. C. Lubensky,
Principles of Con-densed Matter Physics (Cambridge University Press, Cam-bridge, 2000). W. B. Russel, D. A. Saville, and W. R. Schowalter,
Colloidal Dispersions (Cambridge University Press, Cam-bridge, 1992). U. M. B. Marconi and P. Tarazona, J. Chem. Phys. ,164901 (2006). A. J. Archer, J. Phys.: Condens. Matter , 5617 (2006). U. M. B. Marconi and S. Melchionna, J. Chem. Phys. ,184109 (2007). U. M. B. Marconi, P. Tarazona, F. Cecconi, and S. Mel-chionna, J. Phys.-Condes. Matter (2008). U. M. B. Marconi and S. Melchionna, J. Chem. Phys. (2009). A. J. Archer, J. Chem. Phys. , 014509 (2009). Note that the free energy in Eq. (30) is, strictly speak-ing, a grand canonical quantity. This raises the issue ofthe validity of Eq. (30) for describing the free energy ofspecies s , which has only one particle and should strictlybe treated in the microcanonical ensemble. However, if onedoes consider the statistical mechanics of a single particlein a trap (external potential) then one finds that one canwrite the ideal-gas contribution to the free energy of thisparticle as F id [ ρ s ] = k B T R d r ρ s ( r ) ln Λ ρ s ( r ), i.e. differingfrom the expression in Eq. (30) by the term k B T R d r ρ s ( r ).However, due to the fact that the density profile ρ s ( r ) isnormalized (only one s particle), then this term only con-tributes an additional (irrelevant) constant to the free en-ergy (30), and this difference can be ignored. See Ref. 78for further discussion of this issue. M. Oettel, J. Phys.: Condens. Matter , 429 (2005). P. Attard,
Thermodynamics and Statistical Mechanics:Equilibrium by Entropy Maximisation (Academic Press, London, 2002). M. P. Allen and D. J. Tildesley,
Computer Simulation ofLiquids (Oxford University Press, Oxford, 1987). G. N¨agele, J. Bergenholtz, and J. K. G. Dhont, J. Chem.Phys. , 7037 (1999). J. F. Brady, J. Fluid Mech. , 109 (1994). J. F. Brady, J. Chem. Phys. , 567 (1993). D. R. Foss and J. F. Brady, J. Rheol. , 629 (2000). E. R. Weeks, J. C. Crocker, A. C. Levitt, A. Schofield, andD. A. Weitz, Science , 627 (2000). We have also calculated the external potentials correspond-ing to using the zero-dimensionality route for calculating g ( r ) and found that the external potentials are largely un-changed. However, in the zero dimensionality route themagnitude of the repulsion in u d ( r, t ) is much smaller whichresults in more infilling in the core region. D. Reguera and H. Reiss, J. Chem. Phys. , 2558 (2004). S. Auer and D. Frenkel, Nature , 711 (2001). Y. Rosenfeld, Phys. Rev. Lett. , 980 (1989). P. Tarazona, J. A. Cuesta, and Y. Martinez-Raton, Lect.Notes Phys. , 247 (2008). R. Roth, J. Phys.: Condens. Matter , 063102 (2010). E. J. Saltzman and K. S. Schweizer, J. Chem. Phys. ,1197 (2003). K. S. Schweizer, J. Chem. Phys. , 244501 (2005). J. P. Stoessel and P. G. Wolynes, J. Chem. Phys. , 4502(1984). Y. Singh, J. P. Stoessel, and P. G. Wolynes, Phys. Rev.Lett. , 1059 (1985). X. Xia and P. G. Wolynes, Phys. Rev. Lett. , 5526(2001). K. Kim and T. Munakata, Phys. Rev. E , 21502 (2003). C. Kaur and S. P. Das, Phys. Rev. Lett. , 2062 (2001). C. Kaur and S. P. Das, Phys. Rev. E , 26123 (2002). M. Baus and J. L. Colot, J. Phys. C: Solid State Phys. ,L135 (1986). H. L¨owen, J. Phys.: Condens. Matter , 8477 (1990). L. M. Lust, O. T. Valls, and C. Dasgupta, Phys. Rev. E , 1787 (1993). M. Rex and H. L¨owen, Phys. Rev. Lett. , 148302(2008). M. Rex and H. L¨owen, Euro. Phys. J. E , 139 (2009). B. G¨otzelmann and S. Dietrich, Phys. Rev. E , 2993(1997). P. Ram´ırez-Gonz´alez and M. Medina-Noyola (2010), pri-vate communication. T. Koide, G. Krein, and R. O. Ramos, Phys. Lett. B ,96 (2006).78