G Rickayzen
email: gerald.rickayzen@physics.org

D M Heyes
email: dheyes@bournemouth.ac.uk

Bounded Soft-sphere fluids by Molecular Dynamics and Integral Equation Techniques

Keywords: bounded soft spheres, integral equation, molecular dynamics, gaussian core model

The structural properties are determined for fluids composed of particles interacting with the Bounded Soft-Sphere (BSS) potential, V(r)=εσ 2n /(a 2 +r 2 ) n , where a is a variable parameter, r is the separation between the particle centres, n is an integer exponent, and ε and σ, respectively, set the energy and length scales of the potential. The density and temperature dependence of the radial distribution function, g(r), of the BSS fluid have been calculated using Molecular Dynamics (MD) simulation and by numerical solution of the Ornstein-Zernike, OZ, integral equation with various closures. Comparisons with the Gaussian Core Model (GCM) potential fluid are also made. As previously found for the GCM, the Hypernetted Chain (HNC) closure of the OZ equation reproduces the MD g(r) very well over essentially the complete fluid range. Another integral equation approach, which we call the Mean Field Approximation (MFA), is applied and gives comparable accuracy to the HNC OZ closure for not too high values of the potential at the origin.

Part of the equation of state of these fluids is computed and the exact low and high density limits determined. It is shown that the BSS and GCM fluids in the high density limit are far less compressible than the ideal gas at the same density and temperature.

I. Introduction

The physical properties and phase behaviour of systems whose particles interact via pair potentials that are positive and finite at zero separation is a subject which has attracted much interest in recent years. Such potentials can be used as coarse-grained representations of certain soft matter systems, notably polymeric molecules in solution [1,2], and are useful as interaction spreading kernels in various mesoscale dynamical (fluid mechanics-related) modelling techniques such as Dissipative Particle Dynamics [3] and Smooth Particle Hydrodynamics [4]. The 'bounded' nature of the potential introduces some new features into the phase diagram which are not seen in potential systems whose interaction potentials diverge at the origin. For example, the 'bounded' system can pass through a fluid-to-solid phase transition with increasing density and sufficiently low temperature, and then melt again at even higher density (so-called 're-entrant' melting) [2].

Alternatively, depending on the analytic form of the potential, multiple occupancy of the particles on the lattice sites can occur at sufficiently high density (thereby avoiding re-entrant melting) [5].

Recent developments include new duality relations, which link the energy of states of soft pair interactions to the corresponding energy of the Fourier-transformed ('dual') potential [6,7].

Perhaps the most widely studied example of a fluid generated from particles interacting via a purely repulsive bounded potential is that formed from a Gaussian pair interaction. The static properties and phase behaviour of this fluid were first considered by Stillinger and co-workers [8][9][10] and in recent years by many researchers (see for example Refs. [11][12][13][14][15]), in the context of modelling soft matter. Archer and Evans applied a mean-field equation of state to calculate the density-concentration phase diagrams of repulsive Gaussian particles over a range of size ratios [16]. The Gaussian interaction or 'Gaussian Core Model' (GCM) can be written in the following analytic form,

V r V - = (1)
where V(0)=Mε for M>0, ε=k B T , k B is Boltzmann's constant and T is the temperature. The reduced temperature T * ≡ k B T/ε is conveniently set to unity, as for this study. Energy is expressed in terms of ε and distance in σ. A convenient (non-dimensional) measure of the density of the system is the nominal volume or packing fraction, ζ, defined as

3 6 ρσ π ς = (2)
where ρ is the number of particles per unit volume. Both ρ and ζ can be used as measures of the particle number density. The Gaussian potential tends to the hard potential system in the limit of zero temperature but with a diverging effective hard sphere diameter tending to infinity in that limit, so there is no simple path from this potential to those typically used to represent the class of 'simple liquids' dominated by excluded volume effects (such as the soft sphere or inverse power potentials [17]). A recent study of the Gaussian potential has considered its phase behaviour in high Euclidean dimensions [18].

An alternative simple analytic form for a bounded repulsive potential is,

n n a r r V ) ( ) ( 2 2 2 + = εσ (3)
which for |a| > 0 is bounded (i.e., 0 < V(0) < ∞). Assuming a>0, equation (3) has the useful feature that it goes over smoothly to the soft-sphere system as a→0+; it is also longer ranged than the Gaussian potential. We refer to this potential as the Bounded Soft-Sphere or BSS potential. Unlike the GCM, the BSS therefore can be 'transformed' to a standard model simple fluid potential (i.e., soft-sphere potential) at arbitrary temperature through a continuous and differentiable transformation process for a, and as a→0. The static and dynamical properties of the soft-sphere fluid itself (i.e., a=0) have been mapped out in some detail in recent years (e.g. ref. [17]). In fact, the soft-sphere fluid could be used as a reference or limiting case of the BSS fluid for a given n and arbitrary a. The static and dynamical properties of the BSS fluid as a function of density are still largely unknown. To ensure thermodynamic convergence in the a→0 limit, for integer n values, we require 2 n ≥ .

The GCM is a special case of a more general class of potentials with the form, V(r)=Mε exp(-r m ), which we refer to as the Generalised Exponential Potential (GEP). In figure 1 the V(r) for m=1/2, 2 and 4 are compared with the BSS potential, taking n=2 and n=6. The GCM potential is the GEP with m=2. For notational convenience we set ε≡1 and σ≡1 in the equations below. The curvature of the GCM potential at the origin is -2V(0), and for m>2 it is zero (the function is flat at the origin), and for the BSS potential it is -2nV(0)/a 2 . The GCM is therefore 'flatter' close to the origin than the n > 2 examples of the BSS potential with the same V(0) and a 2 < 1 . The point of inflection, r i , of a bounded potential with these analytic forms is also the value of r where the force is a maximum. For the GEP, r i,= [(m-1)/m] 1/m and for the BSS potential, r i = a/(2n+1) 1/2 which indicates inter alia that the soft-sphere potential (a=0) does not have a point of inflection. The value of a is set by the desired value of the potential at r=0, i.e., a=V(0) -1/2n . The range of ζ accessible to the bounded potentials can be essentially infinite in contrast to those fluids formed if the interaction potential has an excluded volume dominance at short range (i.e., r<σ), the usual 'simple liquid' case. The bounded potential fluids can exist in states at high density which are structurally fluid-like, almost ideal gas-like, but whose thermodynamic behaviour can be quite different from an ideal gas or fluids formed at more usual densities (i.e., ζ < 1), see below. The soft-sphere potential is widely used in many branches of liquid state physics, with applications to, for example, liquid metals [19], and hence its extension to include a bounded form could be of potential interest in various practical areas. For example, some classes of polymeric system might be better represented by the BSS potential than the GCM, as the 'steepness' of the potential can be 'tuned' with the parameter 'n' in the BSS analytic form. The BSS potential has a longer tail than the GCM, and not all polymeric system effective potentials are necessarily Gaussian in form. Some of these complex particles will have longer tails than the Gaussian which dies very sharply and it is therefore worth exploring whether this results in somewhat different properties. In fact, the soft-sphere potential has been used to represent microgels [20], and the BSS form might find application to represent microgels with low cross-linked densities. The thermodynamic stability of bounded Lennard-Jones fluids has been investigated by the authors, [21][22][23] and the behaviour of a purely repulsive inverse power system is a natural development of these studies.

The GCM system is fluid above a certain temperature, with a solid region at intermediate densities and lower temperatures. The present definition of reduced temperature is different to that found in previous works, where the reduced value of the potential at the origin was set to 1, and various temperatures, t*<<1 were typically studied in order to investigate the fluid-solid coexistence and surrounding regions. In the present work the reduced temperature was set to T*=1 (the usual softsphere choice) for all cases studied, and the value of the potential at the origin, V(0), was allowed to vary. The largest t* where there is solid-fluid coexistence for the GCM in 3D is t*=0.01, [15] which is equivalent to V(0) =100 in the present study. Most of the reported GCM and BSS simulations here were for V(0)=2 and a few for V(0)=20, a regime most appropriate to polymer solutions. The structure factor, S(k), for large V(0) (small t*) values has been calculated using the HNC closure of the OZ equation. It was found that at ρ=0.225 value for the Gaussian potential (where the highest melting temperature is found for this potential) the Hansen-Verlet freezing criterion (the height of the first maximum in S(k)=2.85) occurs at t*=0.007 (or equivalently V(0)=143) which is about 70% of the literature value. The corresponding calculations for the BSS potential indicate that the maximum t* is even lower. While this clearly is a highly approximate analysis of potential phase boundaries, it does indicate that the BSS states considered in the paper are even further into the fluid state than even those of the GCM. Also there is no evidence from the radial distribution function from the Molecular Dynamics simulations that we have entered a solid or pre-solid part of the phase diagram. For most of the figures the Gaussian and BSS system properties were plotted.

The trends in the behaviour revealed in the figures are qualitatively very similar for the BSS and Gaussian potentials, and these Gaussian state points are known to be in the fluid phase, at an effective a t* well above those where there is any solid in the phase diagram.

Theory

In this Section, the theoretical and computational aspects of this study are described.

2a Ornstein-Zernike

The radial distribution function, g(r), is a fundamental function of any fluid, whose density and temperature dependence gives the equation of state and other static properties. The total correlation function, h(r) = g(r)-1, can be obtained from a given potential form from the direct correlation function, c(r), using the Ornstein-Zernike (OZ) equation,

∫ ′ ′ - ′ + = ) ( |) (| ) ( ) ( r h r r c r d r c r h ρ , (4)
which in reciprocal space is

) ( 1 ) ( ) ( ~k c k c k h ρ - = , ( 5
)
where k is the reciprocal space inverse length. For a generic bounded function, f(r)

∫ ∞ = 0 ) ( ) sin( 4 ) ( ~dr r f kr r k k f π . (6)
h(r) and hence g(r) can be obtained either analytically or numerically from equation (5) using,

∫ ∞ = 0 2 ) ( ) sin( 2 1 ) ( dk k f kr k r r f π . ( 7
)
Iterative solutions for h(r) from equation ( 5) using equations ( 6) and ( 7) have been performed here closures for c(r) (for the numerical procedure, see ref. [24]).

The exact asymptotic form of c(r) forms the basis of the Mean Spherical Approximation (MSA) closure, first introduced in liquid state theory by Lebowitz and Percus in 1966 [25]. Although for simple liquid non-bounded potentials, the MSA is not accurate under normal liquid state conditions, for the GCM bounded potential it has found more success. For the GCM it can be solved analytically (in reciprocal space) and over a reasonably large density and temperature range proves to be quite accurate, when compared with molecular simulation results [2]. In the MSA the closure of the OZ equation is

) ( ) ( r V r c β - = , (8)
where, β=1/ k B T. The diminished importance of the excluded volume term in the bounded potentials is presumably the reason the MSA approximation is more successful for the bounded potential fluids considered here.

Substituting equation (8) in equation ( 6) yields V ~(k), and hence from equation (5),

) ( 1 ) ( ) ( ~k V k V k h ρβ β + - = (9)
where

) ( ) ( ~k V k c β - =
. For the GCM of equation (1),

) 4 / exp( ) ( ~2 2 3 2 / 3 σ εβσ π k k c - - = (10)
Equation ( 7) is needed to obtain h(r) and hence the radial distribution function, g(r) from equation (9). We also note that,

) ( 1 1 ) ( 1 1 ) ( k V k c k S βρ ρ + = - = (11)
the last expression applying only to the MSA approximation. Within the MSA framework, if ) ( ~k V is positive for all k then at high density re-entrant melting rather than clustering on a crystalline lattice is predicted [5]. For the BSS potential, ( ) ( ) ( )

0 2 2 sin 4 n r kr V k dr k r a π ∞ = + ∫ % (12)
This can be evaluated analytically when n is a positive integer. One way to do this is to use the generating function (

)

2 2 1 2 2 2 2 0 2 2 1 , 1 1 n n r a t F a t t r a r a r a t ∞ +   + -   = =     + + +     = + + ∑ (13)
Then the integral in equation ( 12) is the coefficient of ( )

1 n t - -
in the Taylor expansion of ( ) ( )

2 2 2 2 0 sin 4 2 exp r kr dr k a t k r a t k π π ∞   = - +     + + ∫ (14)
The expressions for

) ( ~k V
with the values n=2 to 6 are (from the symbolic mathematics software package, MAPLE [26]) are, ( ) ( ) As these are positive for all k we expect the BSS fluids to exhibit a re-entrant melting transition [5],

at least for n up to 6, with any given closure.

An alternative to the OZ integral equation approach is proposed in the next Section.

2b Mean Field Approximation (MFA)

Consider the self-consistent potential, V S (r), in a fluid at a position r when there is an additional similar molecule centred at the origin. This is made up of the potential V(r) directly due to the molecule at the origin together with the change in the potential due to the change in the local density of all the other molecules as a result of the additional molecule at the origin. Hence,

( ) ( ) ( ) ( ) 0 S V V d V ρ ρ ′ ′ ′ = + - -     ∫ r r r r r r (16)
where ( ) ρ r is the density at position r when the additional molecule is at the origin and 0 ρ is the equilibrium density. The quantity in the square bracket is the local change in the fluid density due to the additional molecule.

If we assume that the distribution of molecules is determined by this potential, V s (r), then If this expression for ( ) ρ r is substituted in equation ( 16) we obtain a self-consistent equation for the potential V S (r), namely

( ) ( ) ( ) ( ) ( ) ( ) 0 0 0 exp exp S B S g V k T V ρ ρ ρ ρ β = = - = - r r r r (17) .
( ) ( ) ( ) ( ) { } 0 exp 1 S S V V d V V ρ β ′ ′ ′ = + - - -     ∫ r r r r r r (18)
.

From this and equation ( 17), one obtains the self-consistent equation for g(r),

( ) ( ) ( ) ( ) 0 ( ) exp 1 g r V r d V g β βρ   ′ ′ ′ = - - - -   ∫ r r r r (19) Note that 1 1 ) ( << - ′ r g
for r′ > 1 and therefore the integral is dominated by these small values of r′ (i.e., r′ <1) which means that one only needs to integrate in effect over small values of r′ , helping convergence. At high temperatures one can expand the exponential in equation ( 18) in powers of β . Retaining the first term only leads to

( ) ( ) ( ) ( ) 0 S S V V d V V βρ ′ ′ ′ = - - ∫ r r r r r r (20)
or its Fourier transform

) ( ) ( ) ( ) ( ~0 k V k V k V k V s s βρ - = (21)
with the solution

) ( 1 ) ( ) ( ~0 k V k V k V s βρ + = (22)
If we use the same high temperature linear approximation in the pair distribution function in equation ( 17), we find

) ( 1 )] ( exp[ 1 ) ( ) ( r V r V r g r h s s β β - ≈ - - = - = (23) Hence, as ) ( ) ( ~k V k h s β - =
, to first order in β from equation (22),

F o r P e e r R e v i e w O n l y 13 ) ( 1 ) ( ) ( ~0 k V k V k h βρ β + - = (24)
. This is just the MSA solution of the OZ equation. From this point of view, the MSA can be seen to be the high temperature approximation to the MFA. Is the MFA a better approximation than the MSA? In order to decide this, we need to solve equation ( 18) and determine g(r) from equation ( 17) (setting V s ≡V to initiate the iteration sequence). One advantage of the MFA over the MSA is that g(r), because it is the exponential of a real function, will be positive definite, as it should be.

The MFA has already been used to study the behaviour of Gaussian-type soft molecule fluids by Likos and co-workers [2] although they derive it from a different point of view. They begin from the HNC closure of the OZ equation, in the form,

( ) ( ) ( ) ( ) 0 ( ) exp 1 g r V r d c g β ρ   ′ ′ ′ = - + - -   ∫ r r r r (25)
but, instead of combining it with the OZ equation to obtain both c(r) and g(r), they replace c(r) by the mean spherical approximation, equation (23). The result is the MFA equation ( 19) for g(r).

Actually one should also include the possibility that if there is a molecule at r there will also be another one at ′ r . This requires a further factor g(r-′ r ) in the integrand of equation ( 16). This equation then becomes

( ) ( ) ( ) ( ) ( ) 0 S V V d V g ρ ρ ′ ′ ′ ′ = + - - -     ∫ r r r r r r r r (26)
. This version of the MFA we call MFA2 but it proved to show convergence problems probably owing to there being three functions in the integral rather than the usual two (as for MFA and the OZ equation). MFA2 will not be considered further here. Convergence of the integral equation was achieved by the standard Picard method of mixing in a fraction of the new solution with the old (see ref. [24]), using V s (r) as the function to be converged. Even for MFA, convergence for typically ζ>10 could be problematic by this method. It could be that for MFA and indeed MFA2, more sophisticated convergence procedures would resolve this problem [27].

2c Molecular Dynamics (MD) simulations

The radial distribution function was computed using the histogram method, by partitioning r into intervals, r k ± ∆r/2, where r k =k∆r. The average value of g(r) in this interval is

k k ij k V r N r g ∆ > < = ρ ) ( ) ( (27)
where <N ij (r k )> is the average number of particles j around a given particle, i, in the interval r k ± ∆r/2. ρ is the number density and the volume element,

∆V k is [ ] 3 3 ) 2 / ( ) 2 / ( 3 4 r r r r V k k k ∆ - - ∆ + = ∆ π (28)
Because the volume element decreases rapidly in magnitude as r tends to zero, the statistical uncertainty in the MD-generated g(r) increases markedly in this limit, and quite long MD simulations are required to obtain reasonable precision for r→ 0. Molecular Dynamics (MD) simulations were carried out for typically 1 million time steps of magnitude 0.002 reduced units for N=1372 particles in the simulation cell [START_REF] Heyes | The Liquid State[END_REF]. Packing fractions in the range 0.1 < ζ < 10 were simulated.

Results and Discussion

In this Section the structural and static properties of the BSS fluid are investigated as a function of density and temperature using several analytic and computational techniques. These results are compared with those obtained for the GCM potential, which has been studied extensively in previous publications. The following figures explore and compare the static structural and thermodynamic properties of GCM and BSS fluids as a function of density and potential height at the origin, V(0), for T * =1 . The available parameter space is rather large for the BSS systems, and only selected but representative state points can be covered in such an initial study. In figures 2-8, the performance of the various integral equation solutions at selected densities are compared with 'exact' simulation results. The MSA undershoots the MD data to a significant extent. At the much higher density of ζ=10, for the Gaussian potential (V(0)=2), the HNC and MFA solutions are equally very good at reproducing the MD data, and now the MSA is also very good also (see figure 4). The PY and RY OZ closures give very poor agreement with the MD data as r tends to zero. The effective pair potential function, V s (r), in the MFA solution is given in figure 6 for the Gaussian potential with V(0)=2. V s (r) follows V(r) in a qualitative way, in having its maximum value at r=0 and decaying monotonically to zero with increasing distance. In the low density region (see the ζ=0.01 curve) V s (r) tends to be very close to V(r). The magnitude and range of V s (r) decreases with increasing fluid density, which is not surprising as the differences represent the effects of the other molecules (i.e., the fluid as a whole). Qualitatively the same trends are evident for the BSS potential (n=3), as revealed in figure 7. In the low density limit V s (r) tends to V(r) more rapidly (see the ζ=0.01 curve, which is hardly distinguishable from V(r)). Also on the figure is shown the HNC c(r) for ζ=0.01 , which can be seen to converge to the MSA form, -βV(r) for large r, as might be expected from this form for the closure, but not for all r as is assumed in the MSA approximation. is not as good with the MD data at this high value of V(0). MFA performs much better at V(0)=2. A wider range of packing fractions is explored in figure 9, which shows this time, g(0) against ζ for the Gaussian and BSS potentials with various V(0) values. For both potential forms, the HNC closure predictions agree well with the corresponding MD-derived values, even for V(0)=20. As density increases to very high values (>>1), the radial distribution function tends to unity for all but the smallest of distance ranges (compared to σ). Therefore the total correlation function is zero for most of the structurally important distance range, r. In this limit, the MSA, MFA and HNC approximations g(r) converge (for not too high V(0) in the MFA case). For the HNC and MSA closures of the OZ equation, the direct correlation function, c(r), is related to h(r) through, (see ref. [START_REF] Hansen | Theory of Simple Liquids[END_REF] p.89). The HNC closure is,

] 1 ) ( ln[ ) ( ) ( ) ( + - + - = r h r h r V r c β (29)
When h(r) tends to zero, the HNC approximation converges to the MSA solution. 18 very high density, so the system cannot be described as being completely 'ideal' structurally (i.e., g(r)=1 for all r) even at these extreme densities. The short range structure in g(r) is somewhat more accentuated and oscillatory in the BSS case, possibly because this potential is less flat near the origin than the GCM example.

Figure 12 shows the density dependence of ) 0 ( h ρ , which is seen to tend to -1 for volume fractions in excess of about 10 for the GCM and 100 for the BSS systems. This limit follows automatically from equation ( 9), taking the high density limit (ρ>>1), and hence is a characteristic feature of bounded potential fluids which is not shared by the 'excluded volume' class of 'simple' liquids (e.g, Lennard-Jones or soft sphere with reasonably high exponent) which have a narrower accessible ρrange for the equilibrium fluid.

The ζ-dependence of the compressibility factor, Z, is shown in figure 13, where (Z-1)/ζ is plotted as function of ζ for the GCM and BSS fluids with V(0)=2. Z was computed from the radial distribution function numerically evaluated and included in the virial expression for the pressure (see ref. [START_REF] Hansen | Theory of Simple Liquids[END_REF], p. 32). Values obtained directly from the MD computations using the virial expression directly are shown also on the figure to give good agreement with the HNC curves. The low density limit is known exactly from the second virial coefficient, and the high density limit can be derived from the virial expression assuming g(r)=1 for all r. These limits are referred to as b 2 and c 2 , respectively, and are defined through the equations, These exact lower and upper bounds for the equation of state were evaluated by numerical integration independently and are given on the figure as horizontal lines (annotated as B2 and CB2, respectively). The good agreement found between the HNC integral equation and MD simulation routes to the equation of state at intermediate densities suggests a relatively simple accurate closedform equation of state for the fluid phase of these bounded potential systems could be obtained.

∞ → + + = → + + = ζ ζ ζ ζ ζ ζ ) ( 1 , 0 ) ( 1 2 2 2 2 O c Z O b Z (30) where, ( ) ∫ ∞
In the large ζ limit, figure 13 reveals that for both classes of bounded potential, Z increases linearly with ζ to a good approximation. As the isothermal compressibility, χ T , is given by,

T B T T Z Z T k P V V             + ∂ ∂ =       ∂ ∂ - = ρ ρ ρ χ 1 1 1 (33)
Then for the bounded potential fluids, χ T decreases as ∼ ρ -2 at high density while the ideal gas (i.e., Z=1) decreases more slowly with density (∼ρ -1 ). The behaviour of χ T at high density has consequences for the zero wave vector component of the structure factor, S(0), which is related by S(0)= ρk B T χ T . At high densities, therefore, it is expected that S(0)∼ ρ -1 . But as, ). 0 potential fluids become more rapidly incompressible in the high density limit than an ideal gas, even though the radial distribution functions of the two classes of system become increasingly similar, apart from a small region close to r=0 which progressively diminishes in size, as revealed in figure 11.

A possible reason why the particles are uncorrelated is that, at high density, the soft particles can overlap. At very high density the density outweighs the interaction; the repulsive interaction is sufficient to ensure that the particles do not simply pile up, the density becomes uniform, and the particles are essentially uncorrelated. This contrasts with the case of hard particles where the hard centres ensure that centres cannot get very close and the particles are correlated to high order distribution function level. This might explain why the variance in the number for large spherical volumes would scale to a power of the sphere radius lower than r 3 , leading to S(k)→0 for k→0. This argument would hold for the Gaussian and BSS potentials. So-called, Hyperuniform systems consist of point patterns that do not possess infinite-wavelength fluctuations, also satisfy the limit of S(0) vanishing in the limit of zero wave vector [START_REF] Torquato | [END_REF]31] . The GCM and BSS clearly share the same feature in the high density limit previously also found to be the case for the GCM [2]. A Mean Field Approximation (MFA) integral equation proposed here is as good at reproducing the radial distribution function as the HNC for not too high values of the bounded potential at zero particle separation. Despite this limitation, the MFA formulation has the physically appealing feature that the local structure in the fluid can be expressed in terms of a 'self-consistent' potential, V s (r), which reduces to the pair potential in the low density limit. The equation of state of the GCM and BSS fluids at high very density reflects a decrease in the isothermal compressibility which appears to vary approximately as the inverse second power of the density. This is a more rapid decrease with increasing density than that of the ideal (gas) fluid, where the compressibility decreases only as the inverse power of density. On the figure this can be seen to be a good approximation for ca. r ≥ 0.9 (taking the HNC to be nearly exact here). We thank the referee for reviewing the revised manuscript and adding further clarification of the original points, which we try our best to address below and a further revision.

Referee(s)' Comments to Author:

Referee: 1 Comments to the Author I appreciate the authors' attempt to incorporate additional references from the literature. However, the discussion of these references is a bit cursory. In particular, the mention of hyperuniformity in the manuscript is supported by no physical arguments to suggest why these systems are becoming hyperuniform in the high-density regime, particularly since the pair correlation functions are essentially unity in this regime (as pointed out by the authors on pages [18][19]. Why does the pair correlation function become decorrelated at high densities? Why is this being accompanied by hyperuniformity? I can find no satisfactory answers to these questions. The authors need to do a better job here.

As we mentioned in the revised manuscript, our route to the conclusion that S(k)→0 for k→0 relied on the fact that the density of these systems is effectively unbounded so that from equation (9), in the very high density limit ) 0 ( h ρ tend to -1. Hence S(0) tends to 0. Our understanding of the hyperuniformity treatment is that density is not considered as a parameter. The success of the HNC closure may also be another reflection of the lack of high order structural correlations. We have added the text below on page 20 of the new manuscript:

A possible reason why the particles are uncorrelated is that, at high density, the soft particles can overlap. At very high density the density outweighs the interaction; the repulsive interaction is sufficient to ensure that the particles do not simply pile up, the density becomes uniform, and the particles are essentially uncorrelated. This contrasts with the case of hard particles where the hard centres ensure that centres cannot get very close and the particles are correlated to high order distribution function level. This might explain why the variance in the number for large spherical volumes would scale to a power of the sphere radius lower than r 3 , leading to S(k)→0 for k→0. This argument would hold for the Gaussian and BSS potentials.

The authors also mention the duality papers but do not, based on their comments, find the results relevant to their study. I disagree with their statements. Although it is true that they are concerned in particular with fluid-phase behavior, it is essential to have some idea of where the fluid-solid boundary is on the phase diagram to provide the appropriate context to understand to the structural properties of the fluid. Are any of their results preempting a fluid-solid phase transition (particularly with the possibility of re-entrant melting for this system)? One cannot answer this question without knowing at least qualitatively the shape of the phase boundary, including the potential ground states. The duality relations are invaluable toward this end. Also, it is not clear to me that the authors realize that it is not essential for the pair potential to be self-similar to apply the duality relations; in particular, in their response to comment 1, they claim that the only useful application of the duality relations has been for the Gaussian core model, a statement which is clearly incorrect.

We did add some extra text in the revised manuscript to argue that our state points are at an effective temperature well above where one would expect the solid part of the phase diagram for these systems. We repeat this text below: "The GCM system is fluid above a certain temperature, with a solid region at intermediate densities and lower temperatures. The present definition of reduced temperature is different to that found in previous works, where the reduced value of the potential at the origin was set to 1, and various temperatures, t*<<1 were typically studied in order to investigate the fluid-solid coexistence and surrounding regions. In the present work the reduced temperature was set to T*=1 for all cases studied and the value of the potential at the origin, V(0), was allowed to vary. The largest t* where there is solid-fluid coexistence for the GCM in 3D is t*=0.01, [15] which is equivalent to V(0) =100 in the present study. Most of the reported GCM and BSS simulations here were for V(0)=2 and a few for V(0)=20. Bearing in mind the similarity between the behaviour of the GCM and BSS systems which is revealed below, the systems reported here should lie well within the fluid part of the phase diagram, which is incidentally the regime most appropriate to polymer solutions."

In response to the referee's second reply we calculated the structure factor, S(k), for large V(0) (small t*) values. Using the HNC closure of the OZ equation we found that at ρ=0.225 value for the Gaussian potential (where the highest melting temperature is found for this potential) the Hansen-Verlet freezing criterion (The height of the first maximum in S(k)=2.85) occurs at t*=0.007 (or equivalently V(0)=143) which is about 70% of the literature value. The corresponding calculations for the BSS potential indicate that the maximum t* is even lower. While this clearly is an approximate method, it does give support for the view expressed above that the BSS states considered in the paper are even further into the fluid state than even those of the GCM. Also there is no evidence from the radial distribution function from the Molecular Dynamics simulations that we have entered a solid or pre-solid part of the phase diagram. For most of the figures the Gaussian and BSS system properties were plotted. The trends in the behaviour revealed in the figures are qualitatively very similar for the BSS and Gaussian potentials. These Gaussian state points are known to be in the fluid phase, at an effective a t* well above those where there is any solid in the phase diagram. So we see no reasonable evidence in these results to expect any interference from an underlying solid phase or pretransitional region. We agree that the phase diagram is definitely worthy of calculating but we feel this needs to be the subject of a separate systematic study. The solid part of the phase is not relevant to the isotherms or iso-V(0) lines reported in the manuscript. We have modified the above additional piece of text a little to strengthen this point: "The GCM system is fluid above a certain temperature, with a solid region at intermediate densities and lower temperatures. The present definition of reduced temperature is different to that found in previous works, where the reduced value of the potential at the origin was set to 1, and various temperatures, t*<<1 were typically studied in order to investigate the fluid-solid coexistence and surrounding regions. In the present work the reduced temperature was set to T*=1 (the usual soft-sphere choice) for all cases studied, and the value of the potential at the origin, V(0), was allowed to vary. The largest t* where there is solid-fluid coexistence for the GCM in 3D is t*=0.01, [15] which is equivalent to V(0) =100 in the present study. Most of the reported GCM and BSS simulations here were for V(0)=2 and a few for V(0)=20, a regime most appropriate to polymer solutions. The structure factor, S(k), for large V(0) (small t*) values has been calculated using the HNC closure of the OZ equation. It was found that at ρ=0.225 value for the Gaussian potential (where the highest melting temperature is found for this potential) the Hansen-Verlet freezing criterion (the height of the first maximum in S(k)=2.85) occurs at t*=0.007 (or equivalently V(0)=143) which is about 70% of the literature value. The corresponding calculations for the BSS potential indicate that the maximum t* is even lower. While this clearly is a highly approximate analysis of potential phase boundaries, it does indicate that the BSS states considered in the paper are even further into the fluid state than even those of the GCM. Also there is no evidence from the radial distribution function from the Molecular Dynamics simulations that we have entered a solid or pre-solid part of the phase diagram. For most of the figures the Gaussian and BSS system properties were plotted. The trends in the behaviour revealed in the figures are qualitatively very similar for the BSS and Gaussian potentials, and these Gaussian state points are known to be in the fluid phase, at an effective a t* well above those where there is any solid in the phase diagram.

We do not see any evidence that the "results preempting a fluid-solid phase transition (particularly with the possibility of re-entrant melting for this system)".

The point about 'duality' only applies at t*=0, and is therefore only relevant if one is interested in that part of the phase diagram encompassing a solid phase, which we argue above, we are not. As far as we understand the papers on this topic and those mentioned by the referee, the only area where duality has been used practically for simple real space potentials is for the Gaussian system where the real space and Fourier transform (FT) potentials are essentially the same. This is not the case for the BSS potential.

I also disagree with the authors' statement that identifying the class of orthogonal polynomials that appear in the dual potential is not relevant to the study. On page 10, the authors use their results to justify re-entrant melting as in the GCM, but only up to n = 6 in the BSS family of potentials. A broader analysis could determine if this behavior holds generally for the class of potentials and could provide additional insight into the phase behavior of these systems (e.g., via duality relations). I agree that it's a minor point and perhaps not necessary for publication, but given that the authors have emphasized the We thank the referee for this suggestion and further clarification of this point. The dual potential was emphasized in the potential only because it is required to solve the OZ equation. The fact that all the coefficients in the polynomial coefficient of the exponential are positive (at least for n up to 6) was used as an aside to comment that it satisfies the accepted literature criterion for reentrant melting. We were not claiming that for all higher n this is the case. We have not been able to prove that this is in general the case for all higher n.

Given that the authors are concerned with mapping the equations of state for the BSS fluids, it is perhaps not essential for them to address all of these issues before publication, and I do agree that the literature on the GCM has a bit of a "head start" compared to this new class of BSS potentials. There are some interesting results in this paper, and future work on these potentials might be warranted. However, I still can find no mention of how physical these potentials are, other than the trivial case where the parameter 'a' vanishes (i.e., the "simple" soft-sphere fluid). Are these results all mathematical curiosities, or do they have applications to actual physical systems? I'm not convinced that it's enough to claim that the BSS potentials are prototypical bounded soft potentials since any number of similar functions could be constructed, most of which would have little or no physical relevance. Can the authors identify an analog of the Flory-Krigbaum effective potential from the GCM literature that would justify this study of BSS potentials? I would expect a complete answer to this question before publication.

Actually, we did address this point in the revised manuscript, which had the following text inserted.

The soft-sphere potential is widely used in many branches of liquid state physics, with applications to, for example, liquid metals [19], and hence its extension to include a bounded form could be of potential interest in various practical areas. For example, some classes of polymeric system might be better represented by the BSS potential than the GCM, as the 'steepness' of the potential can be 'tuned' with the parameter 'n' in the BSS analytic form. In fact, the soft-sphere potential has been used to represent microgels [20], and the BSS form might find application to represent microgels with low cross-linked densities. The thermodynamic stability of bounded Lennard-Jones fluids has been investigated by the authors, [21][22][23] and the behaviour of a purely repulsive inverse power system is a natural development of these studies.

We think the BSS potential is highly relevant to real systems. All linear polymer systems are not necessarily Gaussian in effective form. The BSS potential has a longer, tunable, 'tail' which might apply to some real systems, which could in principle be synthesized. The square mound potential is widely studied in the literature and that is even further from representing a known real system. One does not know in advance where a particular potential form could be useful in the future. In fact, the bounded Lennard-Jones potential We have slightly revised the above text in the new version of the manuscript:

The soft-sphere potential is widely used in many branches of liquid state physics, with applications to, for example, liquid metals [19], and hence its extension to include a bounded form could be of potential interest in various practical areas. For example, some classes of polymeric system might be better represented by the BSS potential than the GCM, as the 'steepness' of the potential can be 'tuned' with the parameter 'n' in the BSS analytic form. The BSS potential has a longer tail than the GCM, and not all polymeric system effective potentials are necessarily Gaussian in form. Some of these complex particles will have longer tails than the Gaussian which dies very sharply and it is therefore worth exploring whether this results in somewhat different properties. In fact, the soft-sphere potential has been used to represent microgels [20], and the BSS form might find application to represent microgels with low cross-linked densities. The thermodynamic stability of bounded Lennard-Jones fluids has been investigated by the authors, [21][22][23] and the behaviour of a purely repulsive inverse power system is a natural development of these studies.

  -known Percus-Yevick (PY), Hypernetted Chain (HNC) and Rogers-Young (RY)

  figures 12-14 and the associated discussion.

Figure 2

 2 figure reveals that the MFA and HNC curves are in very good agreement with the MD data at all r.

Figure 5

 5 Figure5presents the radial distribution function for the BSS potential with V(0)=2 at ζ=1 for three

Figure 8

 8 Figure 8 presents the radial distribution function for the Gaussian and BSS (n=3) fluids over a

Figure 10

 10 Figure 10 presents 1-g(0) for MD, and the MSA and HNC solutions of OZ for volume fractions up to, ζ >>1 using the Gaussian and BSS potentials with V(0)=2 in both cases. This figure

Figure 11

 11 Figure 11 presents g(r) by MSA and HNC for ζ=1000 and V(0)=2. Note the high resolution and

.

  This is consistent with the results displayed in figure 12. The function S(k) for the HNC closure at ζ=1000 and for V(0)=2 is shown for the GCM and BSS fluids in figure 14. This figure shows that S(k) decays monotonically to zero as k→ 0.Therefore, to summarise, as ζ→∞, the limiting density dependence of the isothermal compressibility of these bounded potential fluids departs from that of the ideal gas. The bounded

  It has been shown that for model fluids formed from repulsive particles interacting via the soft Bounded Soft Sphere (BSS) potential, that despite some quantitative differences, they display very similar qualitative trends in the radial distribution function as a function of density to those of the well-studied Gaussian Core Model (GCM) potential fluid. One may therefore expect similarities in static properties expressible in terms of the radial distribution function. When compared with essentially exact Molecular Dynamics simulation data, the radial distribution functions are well represented by the Hypernetted Chain (HNC) closure of the Ornstein-Zernike equation, which was

Figure 1 :

 1 Figure 1: The generalised exponential potential, V(r)=exp(-r m ), where m=1/2, 2 and 4 (top four

Figure 2 :

 2 Figure 2: Radial distribution function for the Gaussian potential with V(0)=2 and VF denotes the volume fraction or ζ. The g(r) from Molecular Dynamics (MD) data (symbols) and various OZ

  The volume fraction is equal to 0.5. To help clarify the trends, the MFA data are shifted upwards by 0.2 and the MSA, PY and RY are shifted upwards by 0.5. The MD data are shifted upwards appropriately to compare with the theoretical curves.

Figure 3 :

 3 Figure 3: As for figure 2, except data for the BSS (n=3) potential are shown. The MFA data are

Figure 4 :

 4 Figure 4: As for figure 2 except that the VF (or ζ) is 10. The MFA data are shifted upwards by 0.2

Figure 5 :

 5 Figure 5: Comparison of the radial distribution function generated by the HNC OZ closure for

Figure 6 :

 6 Figure 6: A comparison between the pair potential, V(r), and MFA self-consistent potential, V s (r) for four values of ζ ('VF' in the figure annotation): 0.01, 0.5, 1.0, and 10 for the Gaussian potential

Figure 7 :

 7 Figure 7: As for figure 6 except that the BSS potential fluid is considered, with n=3. Also shown in this figure is c(r) for the HNC solution. In the MSA case, c(r) = -V(r), as β=1 in the present units.

Figure 8 :

 8 Figure 8: Comparison between the MD and HNC OZ closure g(r) for the Gaussian potential ('G')

Figure 9 :

 9 Figure 9: Comparison between g(0) values as a function of volume fraction for the Gaussian and

Figure 10 :Figure 11 :

 1011 Figure 10: The quantity 1-g(0) against ζ (volume fraction) for the Gaussian and BSS (n=3)

Figure 12 :

 12 Figure 12: This figure shows the volume fraction dependence of ) 0 ( h ρ

Figure 13

 13 Figure 13 This compares (Z-1)/ζ as a function of volume fraction, ζ, for the Gaussian and BSS

Figure 14 S

 14 Figure 14 S(k) from the HNC equations as a function of wavevector, k, for the Gaussian and BSS V(0)=2 fluids at ζ=1000.

  dual potential in the paper, one might have expected a more thorough analysis.

  use in the literature by the thermodynamic integration (free energy) community.

Acknowledgements

It is a great pleasure to contribute to this Special Issue of Molecular Physics in honour of Prof. R. Evans. DMH would like to thank S. Bendall of the School of Design, Engineering and Computing at Bournemouth University for providing invaluable technical assistance in support of these Molecular Dynamics simulations carried out on Bournemouth University computers. The authors would like to thank a referee for bringing their attention to a number of relevant publications for