A Numerical Test of a High-Penetrability Approximation for the One-Dimensional Penetrable-Square-Well Model

A Numerical Test of a High-Penetrability Approximation for the One-Dimensional Penetrable-Square-Well Model

Riccardo Fantoni rfantoni@ts.infn.it National Institute for Theoretical Physics (NITheP), Stellenbosch 7600, South Africa    Achille Giacometti achille@unive.it Dipartimento di Chimica Fisica, Università Ca’ Foscari Venezia, Calle Larga S. Marta DD2137, I-30123 Venezia, Italy    Alexandr Malijevský a.malijevsky@imperial.ac.uk E. Hála Laboratory of Thermodynamics, Institute of Chemical Process Fundamentals
of the ASCR, v. v. i., 165 02 Prague 6, Czech Republic
Department of Physical Chemistry, Institute of Chemical Technology, Prague, 166 28 Praha 6, Czech Republic
   Andrés Santos andres@unex.es; URL:http://www.unex.es/fisteor/andres/ Departamento de Física, Universidad de Extremadura, E-06071 Badajoz, Spain
July 16, 2019

The one-dimensional penetrable-square-well fluid is studied using both analytical tools and specialized Monte Carlo simulations. The model consists of a penetrable core characterized by a finite repulsive energy combined with a short-range attractive well. This is a many-body one-dimensional problem, lacking an exact analytical solution, for which the usual van Hove theorem on the absence of phase transition does not apply. We determine a high-penetrability approximation complementing a similar low-penetrability approximation presented in previous work. This is shown to be equivalent to the usual Debye–Hückel theory for simple charged fluids for which the virial and energy routes are identical. The internal thermodynamic consistency with the compressibility route and the validity of the approximation in describing the radial distribution function is assessed by a comparison against numerical simulations. The Fisher–Widom line separating the oscillatory and monotonic large-distance behavior of the radial distribution function is computed within the high-penetrability approximation and compared with the opposite regime, thus providing a strong indication of the location of the line in all possible regimes. The high-penetrability approximation predicts the existence of a critical point and a spinodal line, but this occurs outside the applicability domain of the theory. We investigate the possibility of a fluid-fluid transition by Gibbs ensemble Monte Carlo techniques, not finding any evidence of such a transition. Additional analytical arguments are given to support this claim. Finally, we find a clustering transition when Ruelle’s stability criterion is not fulfilled. The consequences of these findings on the three-dimensional phase diagrams are also discussed.

I Introduction

Recent advances in chemical synthesis have unveiled more and more the importance of soft-matter systems, such as dispersions of colloidal particles, polymers, and their combinations. Besides their practical interest, these new developments have opened up new theoretical avenues in (at least) two instances. Firstly, it is possible to experimentally fine-tune the details of interactions (range, strength, …), making these systems a unique laboratory for testing highly simplified models within an effective-interaction approach where the microscopic degrees of freedom are integrated out in favor of renormalized macroparticle interactions. Secondly, they offer the possibility of exploring new types of equilibrium phase behaviors not present in the simple-fluid paradigm.

As early as in 1989, Marquest and WittenMarquest89 () suggested that the experimentally observed crystallization in some copolymer micellar systems could be rationalized on the basis of a bounded interaction, that is, an interaction that does not diverge at the origin. Successive theoretical work showed that this class of bounded or ultra-soft potentials naturally arises as effective interactions between the centers of mass of many soft and flexible macromolecules, such as polymer chains, dendrimers, star polymers, etc. (see, e.g., Ref. Likos01, for a reference on the subject). Two well-studied cases belonging to the above class are the Gaussian core model (GCM) introduced by StillingerStillinger76 () and the penetrable-sphere (PS) model introduced in Refs. Marquest89, ; Likos98, , whose freezing transition turns out to display rather exotic features with no analogue in the atomistic fluid realm.

In the present paper, we shall consider a close relative of the PS model, first introduced in Ref. Santos08, , denoted as the penetrable-square-well (PSW) model, where a short-range attractive tail is added to the PS model just outside the core region. In the limit of infinite repulsive energy, the PS and PSW models reduce to the usual hard-sphere (HS) and square-well (SW) models, respectively.

An additional interesting feature common to both PS and PSW, as well as to all bounded potentials, is the fact that, even when confined to one-dimensional systems, they may exhibit a non-trivial phase diagram due to the penetrability which prevents an exact analytical solution.

This is because particles cannot be lined up on a line with a well-defined and fixed ordering in view of the possibility of reciprocal interpenetration (with some positive energy cost), thus lacking an essential ingredient allowing for an exact solution in the respective hard-core counterparts (HS and SW). It is then particularly useful to discuss some motivated approximations in the one-dimensional model which can then be benchmarked against numerical simulations and subsequently exploited in the much more complex three-dimensional case.

The aim of the present paper is to complete a study on the one-dimensional PSW model started in Refs. Santos08, ; Fantoni09, , as well as the general results presented in Ref. Santos08, which are particularly relevant in the present context. In the first paper of the series,Santos08 () we introduced the model and discussed the range of stability in terms of the attractive versus repulsive energy scale. We provided, in addition, exact analytical results in the low-density limit (second order in the radial distribution function and fourth order in the virial expansion) and a detailed study of the Percus–Yevick and hypernetted-chain integral equations. These were used in the following paperFantoni09 () to propose a low-penetrability approximation at finite density which was then tested against numerical simulation. This low-penetrability approximation is expected to break down in the opposite regime, namely when temperatures and densities are such that particles easily interpenetrate each other. In the present paper, we address this latter regime by proposing a complementary approximation (the high-penetrability approximation) and discussing its range of validity and the relationship with the low-penetrability regime. Note that a similar matching of the low- and high-penetrability approximation has already been carried out by two of the present authors in the framework of the PS model.Malijevsky06 (); MYS07 () It turns out that the high-penetrability approximation in the context of bounded potentials coincides with the linearized Debye–Hückel classical approximation originally introduced in the framework of the Coulomb potential.Hansen86 () It has been recently shownSantos09 () that two of the three standard routes to thermodynamics (the energy and the virial routes) are automatically consistent within the linearized Debye–Hückel approximation, for any potential and dimensionality. This means that a deviation from the third standard route to thermodynamics (the compressibility route) can be exploited to assess the degree of reliability of the high-penetrability (or linearized Debye–Hückel) approximation. This is indeed discussed in the present paper, where we also discuss the full hierarchy of approximations ranging from the full Debye–Hückel approximation to the simplest mean-spherical approximation.

In view of the boundness of the potential, the usual van Hove no-go theoremVanHove50 (); Cuesta04 () on the absence of phase transitions in certain one-dimensional fluids does not hold. It is then natural to ask whether a phase transition occurs in the one-dimensional PSW fluid by noting that the addition of an attractive tail to the pair potential of the PS model extends the question to the fluid-fluid transition, in addition to the fluid-solid transition possible even within the PS model. In the present paper we confine our attention to the fluid-fluid case only and discuss this possibility using both analytical arguments and state-of-the-art numerical simulations.Panagiotopoulos87 (); Panagiotopoulos88 (); Smit89a (); Smit89b () Our results are compatible with the absence of such a transition, as we shall discuss. This is also supported by recent analytical resultsFantoni10 () using a methodology devised for one-dimensional models with long-range interactions.Kastner08 () We discuss possible reasons for this and a plausible scenario for the three-dimensional case.

Finally, we note that the approach to a critical point is frequently anticipated by the so-called Fisher–Widom lineFisher69 () marking the borderline between a region with oscillatory behavior in the long-range domain of the correlation function (above the Fisher–Widom line) and a region of exponential decay. We discuss the location of this line within the high-penetrability approximation and again the matching of this result with that stemming from the low-penetrability approximation.

The structure of the paper is as follows. We define the PSW model in Sec. II. We then construct the high-penetrability approximation in Sec. III and in Sec. IV we discuss some approximations related to it. Section V contains a discussion on the routes to thermodynamics, as predicted by the high-penetrability approximation. The structure predicted by the approximation is compared with Monte Carlo data in Sec. VI. The Fisher–Widom line and the possibility of a fluid-fluid transition are discussed in Sec. VII. The paper ends with some concluding remarks in Sec. VIII.

Ii The penetrable-square-well (PSW) model

The penetrable-square-well (PSW) model is defined by the following pair potentialSantos08 (); Fantoni09 ()


where is the well width and and are two positive constants accounting for the repulsive and attractive parts of the potential, respectively. In the following, we shall restrict our analysis to the case and , where we know the one-dimensional model to be stable with a well defined thermodynamic limit.Santos08 () It is shown in Appendix A that, more generally, the one-dimensional PSW model is guaranteed to be stable if , where is the integer part of . For lower values of the model may or may not be stable and we will come back to this point in Section VII.

An important role in the following is played by the corresponding Mayer function


where is the parameter measuring the degree of penetrability varying between 0 (free penetrability) and 1 (impenetrability), while measures the strength of the well depth. Here with the temperature and the Boltzmann constant, is the Mayer function for hard spheres, and is the Heaviside step function.

A detailed discussion of the limiting cases of the PSW model can be found in Ref. Santos08, . Here we merely note that the PSW Mayer function is immediately related to the usual SW Mayer function by


where we have introduced the ratio . At a given value of , increases quasi-linearly with , its minimum value corresponding to .

Iii A High-Penetrability Approximation (HPA)

In Ref. Fantoni09, we discussed a low-penetrability approximation (LPA) to the PSW model. Within this approximation, one assumes so that the repulsive barrier is sufficiently higher than the thermal energy , penetrability is small, and the system is almost a hard-core one. The advantage of this theoretical scheme is that one can use the general recipe leading to the exact solution for the one-dimensional SW problem —in fact, valid for any potential with a hard core and short-range attractions— and perform some ad hoc adjustments to ensure that some basic physical conditions on the radial distribution function are satisfied. Comparison with Monte Carlo (MC) simulations showed a good behavior of the LPA even for (), provided the density was moderate ().

The opposite limit is also inherently interesting for several reasons. From a physical viewpoint this amounts to starting from the ideal gas limit (one of the common reference systems for simple fluids) and progressively building up interactions by increasing . An additional mathematical advantage stems from the simple observationMalijevsky06 (); MYS07 (); AS04 () that in the (exact) cluster expansion of only the dominant chain diagrams need to be retained at all orders, thus leading to the possibility of an exact summation of those leading contributions. As we shall see shortly, this is in fact a procedure known as the Debye–Hückel approximation in the context of charged fluids.Santos09 (); Hansen86 ()

Our main goal is the computation of the cavity function , from which one can immediately compute the radial distribution function (RDF) . In the PSW case one then has from Eq. (3)


where and, in conformity with previous work,Santos08 (); Fantoni09 () we have redefined all lengths in units of so we set in most of the following equations.

As shown in Ref. AS04, for the PS case, the exact form of the PSW cavity function in the limit at finite and is


where the function is defined through its Fourier transform


being the Fourier transform of . Note that in the limit one has and .

Generalizing an analogous approximation in the context of the PS model,Malijevsky06 (); MYS07 () our high-penetrability approximation (HPA) consists of assuming Eqs. (6) and (7) for finite, but small, values of . It is worth noting that the combination of the two expressions (6) and (7) defines what is usually referred to, in a different context, as the linearized Debye–Hückel (LDH) approximation,Santos09 (); Hansen86 () and this will be further elaborated below.

Equations (6) and (7) hold for any dimensionality. In the specific one-dimensional case, and taking into account Eq. (4), we have


The function can be numerically evaluated in real space by Fourier inversion as


An explicit expression for the density expansion of within the HPA is reported in Eqs. (43) and (44) of Appendix B, where the radius of convergence of the expansion is also analyzed.

From Eqs. (3), (5), and (6) the total correlation function, , within the HPA is easily obtained as


or in Fourier space,


From this equation it is straightforward to get the structure factor


and the Fourier transform of the direct correlation function . The zero wavenumber value of the structure factor is


where we have taken into account that . This completes the calculation of the correlation functions within the HPA.

Iv Approximations related to the HPA

As anticipated in Section III, the HPA is the exact equivalent to the well-known LDH approximation, which is widely used in the context of charged simple fluids.Hansen86 () The latter is actually an intermediate step of a hierarchy of successive approximations ranging from the simplest mean-spherical approximation (MSA) to the full non-linear version of the Debye–Hückel approximation (see Ref. [11] in Ref. Santos09, for a discussion on this point). For the PSW model —and more generally for any bounded potential— the LDH approximation (the HPA in the present language) is particularly relevant in view of the fact that one can make arbitrarily small by letting , thus justifying the approximation of neglecting non-chain diagrams. It is then interesting to check the performance of the other approximations included in the aforementioned class, which will be translated in the present context for simplicity.

On top of the hierarchy of approximations there is the non-linear HPA (nlHPA)


which is equivalent to the non-linear Debye–Hückel approximation, as remarked. The HPA, Eq. (6), is obtained upon linearizing the exponential, an approximation valid again in the limit . An additional approximation —denoted here as the modified HPA (mHPA)— can be considered with the help of Eq. (10) by neglecting the quadratic term in . This yields


This is equivalent to keeping only the first term on the right-hand side of Eq. (11), which implies or, in real space,


The lowest rank in the hierarchy is occupied by the MSA, which is obtained from Eq. (16) upon linearization of the Mayer function ,


Since is a convolution, it must be continuous at and . It follows that the approximations with a continuous cavity function at and are nlHPA and HPA. For instance, in the mHPA (15) the cavity function is , so that one has and .

It has been shown in Ref. Santos09, that the virial and energy routes to thermodynamics (to be discussed below) are consistent one another within the HPA, for any potential and any dimensionality. A similar statement holds true for soft potentials within the MSA.Santos07 () This clearly includes the PSW potential in both cases.

It is interesting to make contact with previous work carried out by Likos et al. Likos01b (); Likos07 () on a general class of unbounded potentials which are free of attractive parts, thus resulting particularly useful in the context of the fluid-solid transition. In Ref. Likos01b, , the MSA given in Eq. (17) along with the spinodal instability to be discussed in detail in Section VII, have been introduced for a general class of unbounded potentials including the PS as a particular case. This has been further elaborated and extended to include the GCM in Ref. Likos07, . In both cases, the authors discuss directly the three-dimensional case, so that a direct comparison with the present work cannot be drawn at the present stage, but they also provide a detailed discussion of various approximations, within the general framework of density functional theory, that provide a unified framework where even the present model could be included.

V Equation of state

Given an approximate solution of a fluid model there are several routes to the equation of state which, in general, give different results. The most common are three:Hansen86 () the virial, the compressibility, and the energy route. The consistency of the outcome of these different routes can be regarded as an assessment on the soundness of the approximation. For some particular approximations it may also happen that the consistency of two of the three routes are automatically enforced (see Ref. Santos09, and references therein for a detailed discussion on this point). This is the case of the HPA, where the virial and energy routes coincide, as anticipated. Hence, the consistency with the compressibility route will provide a rough estimate of the regime of validity of the HPA within a phase diagram for the PSW potential.

Let us briefly recallHansen86 () the methods to compute the compressibility factor, , associated with the three different routes. The virial route is defined by


which, using standard manipulations,Hansen86 () yields


Thus the problem is reduced to the computation of the cavity function which, in the present context, follows from Eqs. (6) and (9).

Next, we consider the compressibility route,


where the integral can be readily evaluated with the help of Eq. (13).

Regarding the energy route, we start from the internal energy per particle


In the above equation, the expressions given by Eqs. (1) and (5) have been used. In order to obtain from we exploit the following standard thermodynamic identitySantos09 ()


thus leading to


We have used the exact consistency between the virial and energy routes within the HPA as a test of the numerical calculations.

Figure 1: Plot of the compressibility factor as computed from the virial(-energy) route to the pressure for and . Results for three different reduced temperatures are displayed.

Figure 1 depicts the results of the virial —and hence energy, as remarked— route under the condition , which constitutes the borderline range of stability of the one-dimensional PSW model with .Santos08 () Under this demanding condition, we have considered three reduced temperatures from to , whereas the width of the well has been fixed to the value . We remark that is not limited in values from above due to the boundness of the potential. The clear downturn of all three curves for sufficiently large reduced density is a consequence of the existence of a maximum density [see Eq. (46)], beyond which the HPA breaks down, as described at the end of Appendix B. In particular, the values of the maximum density for and are , , and at , , and , respectively, in agreement with Fig. 1.

Figure 2: Comparison between the compressibility factor from the virial(-energy) and the compressibility route to the pressure. The circles represent MC simulation results. Chosen parameters are , , and .

We compare in Fig. 2 the results from the virial(-energy) and the compressibility routes with MC simulationsnote1 () for an intermediate value of the reduced temperature () and other parameters as before. Rather interestingly, the virial(-energy) route appears to reproduce rather well the numerical simulation results up to the region where the artificial downward behavior shows up, whereas the compressibility route begins to deviate for densities .

Figure 3: Rough estimate of the region of reliability for the HPA, on the basis of the consistency between the virial(-energy) and compressibility routes, in the reduced density versus reduced temperature plane. Here the curves (a)–(d) correspond to , , , and , respectively. The region below each curve represents states where the relative deviation between the virial(-energy) route and the compressibility one is smaller than and hence regarded as reliable. The inset shows the curves in the versus plane for (dashed line) and (solid line).

As the temperature increases the HPA theory clearly remains a good approximation for a larger range of densities. We can naturally measure this by the requirement that virial(-energy) and compressibility routes are consistent within a few percent. This is indeed shown in Fig. 3, where we depict a transition line separating a “reliable” from an “unreliable” regime, as measured by the relative deviation of the two routes (here taken to be ), for four choices of the model parameters: , , , and . The value is frequently used in the SW counterpart.Vega95 () We observe that the region where the HPA is reliable is hardly dependent on [compare curves (a)&(b) and (c)&(d) in Fig. 3]. On the other hand, at given values of and , the range decreases with increasing [compare curves (a)&(d) and (b)&(c) in Fig. 3], as expected. However, this effect is much less important if the increase of takes place at fixed (see inset of Fig. 3). It is interesting to note that, as illustrated by Fig. 2, the HPA virial route keeps being reliable up to a certain density higher than .

Figure 4: Regions of reliability for the LPA and the HPA in the reduced density versus reduced temperature plane. Here =5 and . The labels LPA, HPA, and LPA+HPA indicate the regions where only the LPA, only the HPA, or both approximations are reliable, respectively.

As said above, in Ref. Fantoni09, we introduced a LPA that was accurate for states where the penetrability effects were low or moderate. The LPA is complemented by the HPA presented in this paper. It is then interesting to compare the regions where each approximation can be considered reliable according to the same criterion as in Fig. 3. This is shown in Fig.  4 for the case and . The two transition lines split the plane into four regions: a region where only the LPA is reliable, a region where only the HPA is reliable, a region where both approximations are reliable (and provide equivalent results), and a region where none of them is sufficiently good. The latter region shrinks as decreases (thanks to the HPA) or decreases (thanks to the LPA).

Vi Structure

As an additional test of the soundness of the HPA, we also study the RDF , which can easily be obtained from Eqs. (5) and (6), or equivalently from Eq. (10), once the auxiliary function has been determined. For a sufficiently high temperature (and hence high penetrability), the HPA is clearly well performing, as can be inferred from Fig. 5, when compared with standard NVT MC results. Here we have considered the same parameters as in the preceding section ( and ) at a corresponding high-temperature value and a density where overlappings are unavoidable. Under these conditions, there is no visible difference among the various approximations considered in Section IV. The excellent performance of the HPA observed in Fig. 5 agrees with the reliability criterion of Fig. 3 since the state and is well below the curve (c) corresponding to and .

Figure 5: Results for the RDF as a function of with , , , and . Predictions from the HPA (solid line) are compared with MC results (circles).
Figure 6: Comparison of different approximations in the results for the RDF vs . Chosen parameters are , , and . MC results (circles) are compared with the HPA (solid line), the nlHPA (dashed line), the mHPA (short dashed line), and the MSA (long dashed line). The two different panels refer to different reduced temperatures: (top panel) and (bottom panel). Note that both panels are drawn within the same scale.

As we cool down, significant differences among various approximations (HPA, nlHPA, mHPA, and MSA) begin to appear, as depicted in Fig. 6, where results corresponding to temperatures (top panel) and (bottom panel) are reported within the same scale. The state and is still lying in the reliable region of Fig. 3, but close to the boundary line (c), while the state and is clearly outside that region. In the case the HPA and its three variants are practically indistinguishable, except in the region , which is very important to describe the correct thermodynamic behavior, where the best agreement with MC data corresponds to the nlHPA, followed by the HPA. The two approximations that do not preserve the continuity of the cavity functions (mHPA and MSA) overestimate the jump at . In the lower temperature case all the approximations overestimate the oscillations of the RDF. Interestingly, the HPA captures quite well the values of near the origin. The worst overall behavior corresponds again to the MSA, which even predicts negative values of for .

Figure 7: An additional comparison of different approximations in the results for the RDF vs . Here fixed parameters are and . MC results (circles) are compared with the HPA (solid line), the nlHPA (dashed line), the mHPA (short dashed line), and the MSA (long dashed line). The top panel refers to the state and , whereas the bottom panel refers to and . In the bottom panel the MSA is not depicted. Again both panels are on the same scale.
Figure 8: Comparison of different approximations in the results for the RDF vs . Chosen parameters are , , , and . MC results (circles) are compared with the HPA (solid line) and the LPA (dashed line).

Additional insights can be obtained by decreasing the range of interactions, in close analogy with what we considered in previous work for the complementary LPA.Fantoni09 () This is reported in Fig. 7 for the case and with (top panel) and (bottom panel). Again we stress that the same scale is used for both panels in order to emphasize the effect of lowering the temperature. Clearly this is a more demanding situation. In fact, both states are above curve (a) of Fig. 3 and thus outside the corresponding reliability region. Therefore, clear deviations from MC results appear in all considered approximations, especially in the lower temperature case (bottom panel). Yet, the HPA is still a reasonably good approximation that follows the main qualitative features of the correct . In the higher temperature case the only noticeable limitations of the HPA practically take place near the origin, this deficiency being largely corrected by the mHPA.

To conclude this section, it is worthwhile comparing the two complementary approaches HPA and LPA at a case where both are expected to be reliable, according to the diagram of Fig. 4. This is done in Fig. 8 for , , , and . We observe that both approximations agree well each other and with MC data, except in the region , where the HPA RDF presents an artificial curvature.

Vii Fisher–Widom line and fluid-fluid transition

We now turn to an interesting point raised in previous work,Fantoni09 () namely the question of whether the model can display a phase transition in spite of its one-dimensional character.

The existence of general theorems —all essentially based on the original van Hove’s resultVanHove50 ()— on the absence of phase transitions for a large class of one-dimensional models with short-range interactions is well established.Cuesta04 () PSW and PS models, however, do not belong to the class for which these general theorems hold. This is because boundness allows multiple partial (or even total) overlapping at some energy cost, thus rendering invalid the arguments used in the aforementioned theorems.

On the other hand, none of these theorems provide a general guideline to understand whether a one-dimensional model may or may not display a non-trivial phase transition, and one has then to rely on the specificity of each model.

As discussed in our previous work,Fantoni09 () it is instructive to first address the simpler question of the location of the Fisher–Widom (FW) line. This is a line separating two different regimes for the large-distance behavior of the RDF in the presence of competing repulsive/attractive interactions.Fisher69 () The rationale behind the FW line is that on approaching the critical points where attractions become more and more effective, the behavior of correlation functions must switch from oscillatory (characteristic of repulsive interactions) to exponential with a well defined correlation length . In the previous work, we analyzed the location of this line for PSW within the LPA. Here we extend this analysis to the HPA regime and discuss the compatibility of the two results.

Let us first briefly recall the main points of the analysis, referring to Ref. Fantoni09, for details. From Eq. (5) we note that the asymptotic behavior of is the same as that of . In view of Eq. (6), this is hence related to , whose asymptotic behavior is governed by the pair of conjugate poles of with an imaginary part closest to the origin. If the real part of the pair is zero, the decay is monotonic and oscillatory otherwise.

According to Eq. (7), the poles of are given by


Let be the imaginary pole and be the pole with the imaginary part closest to the origin. The FW line is determined by the condition . This gives, at a given temperature, three equations in the three unknowns , , and .Fantoni09 () More specifically, after some algebra, one gets


The inverse of the parameter represents the correlation length . From a practical point of view it is more convenient to use rather than as a free parameter to construct the FW line. In that way, Eq. (25) becomes a transcendental equation that gives as a function of ; once is known, the solution to Eq. (26) gives as a function of ; finally, insertion of and into Eq. (27) provides . The corresponding values of the pressure are obtained from either Eq. (19) (virial-energy route) or Eq. (20) (compressibility route).

We observe that decreases as decreases, until a critical value is found in the limit . In that limit, Eqs. (25)–(27) simplify to


At the critical point and , one has or, equivalently, . Therefore, at this point does not decay for long distances and in Fourier space one has and for short wave numbers. The condition is also satisfied for if , in agreement with Eq. (13). This defines in the - plane a spinodal line or locus of points of infinite isothermal compressibility (within the compressibility route). The spinodal line cannot be extended to temperatures larger than the critical value because if , where is the maximum density beyond which the HPA is unphysical at a given temperature (see Appendix B). We further note that the spinodal line only has a lower density (or vapor-like) branch, thus hampering the interpretation of as a conventional critical point.

The above features are already suggestive of considering the HPA spinodal line as an artifact of the theory when used in a region of parameter space where the approximation is invalid. Additional support to this view stems from the fact that the HPA keeps predicting a spinodal line and a critical point even in the SW case (, ), a clearly incorrect feature. As we shall discuss further below, specialized numerical simulations coupled with a recent analytical study,Fantoni10 () strongly support the absence of any phase transition in the present one-dimensional PSW model.

Figure 9: Plot of the FW transition line in the vs plane with and (top panel) and (bottom panel). The long dashed curves are the spinodal lines predicted by the HPA. The FW and spinodal liness meet at the critical point (denoted by a circle).

In Fig. 9 we report the comparison between the FW lines, as predicted by the LPA and the HPA, in the vs plane for and two different energy ratios and . The spinodal line predicted by the HPA is also included. As said above, the FW and spinodal lines meet at the critical point. While at high temperatures (above ) there is a remarkable agreement between the two approximations, deviations occur at lower temperatures.

Figure 10: Same as in Fig. 9, but in the vs plane. Note that while in the LPA the three routes to the pressure are not distinguishable one from the other on the graph scale, for the HPA the difference between the virial and the compressibility route is noticeable at low temperatures.

As in the original work by Fisher and Widom,Fisher69 () we also report the FW line in the vs plane (see Fig. 10), where again we compare the lines as derived from the HPA and LPA schemes. For the energy ratio (top panel) we see that the HPA and the LPA give qualitative similar forms of the FW line with a significant deviation at low temperatures, where again the HPA FW line is interrupted at . Note that, while all three standard routes give practically identical results within the LPA, the virial(-energy) route in the HPA differs from the compressibility one at low temperatures (more than for ). For consistency, the HPA spinodal line is obtained via the compressibility route only. Similar features occur for the second lower value of the energy ratio, namely (bottom panel). Again, the LPA and HPA lines are qualitatively similar, the three routes in the LPA provide indistinguishable results, and the virial and energy routes in the HPA deviate more than for . The main distinctive feature in this case is a marked upswing of the tail of the FW line, absent the previous case . This means that, on increasing penetrability —that is, on decreasing — the transition from oscillatory (above the line) to monotonic (below the line) behaviors occurs at a higher pressure and a higher density for a fixed temperature .

Despite the important differences in the steps followed to derive the LPA and the HPA, it is noteworthy that they agree in the qualitative shape of the FW lines (even though the HPA predicts a spurious spinodal line). It is then reasonable to expect that the true FW line should interpolate the LPA line at low temperatures with the HPA line at high temperatures.

Finally, we now tackle the issue of the existence of a phase transition for the PSW one-dimensional model. In view of the HPA results on the seemingly existence of a spinodal line (and hence of a critical point), we here consider the fluid-fluid transition. As we shall see, our numerics is compatible with the absence of such a transition, thus supporting the view that the above findings of a spinodal line is indeed a consequence of the application of the HPA to a regime where the theory is not valid.

As the FW line always anticipates the critical point, as remarked, we can then look for the existence of a fluid-fluid coexistence line in the region predicted by the interpolation of the LPA and HPA Fisher–Widom lines. We have carried out extensive simulations of the PSW fluid using Gibbs ensemble Monte Carlo techniques and employing all standard improvements suggested in the literature.Panagiotopoulos87 (); Panagiotopoulos88 (); Smit89a (); Smit89b () In order to validate our code, we tested it against the case of the one-dimensional SW potential, where exact analytical predictions for all thermodynamic quantities are available.

We have used up to particles and carefully scanned the temperature range and the density range , as suggested by the FW line (see Fig. 9). We have also considered different values of and for cases giving a significant overlapping probability. In all the cases we have not found any signature of a fluid-fluid phase separation.

Although the absence of a critical transition is always much more difficult to assess as opposed to its presence, the first scenario is consistent with more than one indication. The first indication stems from a lattice model counterpart of the one-dimensional PSW model. This is discussed in Appendix C, where the lattice version of the PSW model is constructed following standard manipulations with the result that no phase transitions are present for finite occupancy. An additional evidence supporting the absence of any fluid-fluid or freezing transition stems from the very recent exact analytical work alluded to earlier,Fantoni10 () that, using the methodology presented in Ref. Kastner08, , concludes that no phase transitions are present for the PSW and PS models in one dimension.

Figure 11: Plot of the RDF obtained from MC simulations for , , , and (top panel) and (bottom panel).

In our simulations we have also investigated values of and which violate Ruelle’s stability criterion (see Appendix A) and thus the one-dimensional PSW model is not necessarily stable in the thermodynamic limit. Here the phenomenology turns out to be much more interesting. For sufficiently low temperatures and sufficiently high densities we observe the formation of a “blob” of many-particle clusters (each made of a large number of overlapping particles) having a well defined and regular distribution on the axis, and occupying only a portion of the system length. In the “blob” phases the energy per particle grows with the number of particles, thus revealing the absence of a thermodynamic limit. The transition from a “normal” phase to a “blob” phase if , where is the integer part of , is illustrated in Fig. 11. This figure shows the RDF at (top panel) and (bottom panel) for , , and . At the lower density the structure of the PSW fluid is qualitatively not much different from that expected if (compare, for instance, with Fig. 6). However, the structure at the higher density is reminiscent of that of a solid, except that the distribution of particles does not span the whole length (here ). Instead, the particles distribute into a few clusters regularly spaced with distance , so that the particles of a given cluster interact attractively with all the particles of the nearest- and next-nearest-neighbor clusters. Note that the number of particles within any given cluster is not necessarily identical. It is also worth noticing that, in spite of the huge difference in the vertical scales of both panels in Fig. 11, they are consistent with the condition .

The decay of the peaks of is mainly due to the lack of translational invariance, i.e., the first and last clusters have only one nearest-neighbor cluster, the first, second, next-to-last, and last clusters have only one next-nearest-neighbor cluster, and so on.

We stress that the above phenomenon is specific of bounded potentials, such as PSW, and has no counterpart in the hard-core domain. It is then plausible to expect their appearance even in the corresponding three-dimensional versions of these models where freezing transition (and phase separation for SW) are present but could both be hampered by the presence of this clustering phenomenon in the region of parameter space where Ruelle’s stability criterion is violated. This would extend the interesting phenomenology already established for the PS case Likos01b (). Work along these lines on the three-dimensional case is underway and will be reported elsewhere.

Viii Conclusions

In this paper we have completed the study initiated in previous workSantos08 (); Fantoni09 () on the penetrable-square-well model in one-dimension. This is a model combining the three main ingredients present in many colloidal-polymer solutions, namely repulsions, attractions, and penetrability. While the first two are ubiquitous even in simple fluids, the latter is a peculiarity of complex fluids where there exist many examples of colloid-polymer systems which are penetrable (with some energy cost) to some extent, and they involve both steric repulsions and short-range attractions. This model then captures all these crucial features at the simplest level of description within an implicit solvent description.

The main new point of this paper was to present an additional and complementary approximation, denoted as the HPA, valid in regimes complementary to those valid for the low-penetrability scheme LPA discussed in Ref. Fantoni09, . While the idea behind the LPA was to modify the exact relations valid for the one-dimensional SW fluid —and in fact for any fluid with a hard core and short-range attraction— to allow penetrability within some reasonable approximation, the driving force behind the HPA is the fact that for bounded potential the Mayer function can be made arbitrarily small by considering sufficiently high temperatures. As a consequence, only the linear chain diagrams need to be retained at each order in the cluster expansion. As it turns out,Santos09 () this is tantamount to considering the celebrated Debye–Hückel theory for charged fluids, and we have here considered the soundness of this approximations at various regimes as compared to specialized MC simulations. The latter were also compared with other approximations which parallel the entire hierarchy of approximations in the framework of charged fluids, ranging from the most sophisticated non-linear Debye–Hückel to the simplest mean-spherical approximations. We have assessed the regime of reliability of these approximations both for thermodynamic and correlation functions by comparison with MC simulations and by internal consistency between different routes to thermodynamics.

Next we have also discussed the location of the FW line, separating oscillatory from monotonic behavior in the correlation function, within the HPA and compared with that obtained from the LPA introduced in previous work.Fantoni09 () In agreement with previous findings, we find that penetrability enhances the region where correlation functions have a monotonic regime. The FW derived from the HPA and LPA schemes are found to be in qualitative agreement thus making it possible the drawing of a line interpolating the high- and low-temperature regimes.

As a final point, we investigated the possibility of a fluid-fluid transition. This possibility arises because the boundness of the potential renders the van Hove theorem on the absence of phase transition for one-dimensional model with short-range interactions non-applicable. In fact, the HPA is seen to predict a critical point where the FW line meets a spinodal line. However, this prediction takes place in a region of densities and temperatures where the HPA is not reliable. A careful investigation using both NVT and Gibbs ensemble Monte Carlo techniques akin to those exploited in the investigations of the analog problem for the three-dimensional SW model yields negative results. These findings are also supported by analytical arguments based on the lattice gas counterpart where the absence of transition can be motivated by the absence of an infinite occupancy of each site, as well as by an exact analytical proofFantoni10 () of a no-go theorem proving the absence of any phase transition in this model, which is in agreement with, and beautifully complements, our work.

In our quest for a possible thermodynamic transition in the one-dimensional PSW model we have explored values of the energy ratio and well width for which the stability of the system in the thermodynamic limit is not guaranteed by Ruelle’s criterion.Fisher66 (); Ruelle69 () We have found that, as the temperature decreases and/or the density increases, a transition from a normal fluid phase to a peculiar solid-like phase takes place. The latter phase is characterized by the formation of clusters of overlapping particles occupying a small fraction of the available space and with non-extensive properties. This clustering transition preempts both the fluid-fluid and fluid-solid transitions.

In view of the results presented here, it would be very interesting to discuss the phase diagram of the corresponding three-dimensional PSW model. The phase diagram of the SW model () is indeed well established and includes both a fluid-solid transition —present even in the HS counterpart— and a fluid-fluid transition line. The latter is present for any value of and but is stable against freezing only for , the depth of the well being irrelevant. Kranendonk88 () The results presented here strongly suggest the importance of the additional parameter . A first interesting issue would be the Ruelle instability in three-dimensions. A straightforward extension of the arguments presented in Appendix A predicts a guaranteed stability for (if ), but the actual onset of the instability cannot be assessed through these arguments. One could expect that for sufficiently high penetrability (i.e., and high density) a phenomenon akin to the “clustering” transition found here could be present. In the case , where the clustering transition is not expected, the interesting point is to assess the influence of the ratio on the location of the fluid-fluid critical point and coexistence line. All of this opens the possibility of a rich and interesting phase diagram which would complement that already present in a general class of bounded potentials with no attractive tails. Likos01b (); Likos07 () We note that the high-penetrability regime is indeed the realm of the HPA presented here, which can be obviously extended to three-dimensions.

Work on the three-dimensional PSW model including the above points and other aspects is underway and will be reported elsewhere.

One of us (RF) wishes to thank M. Kastner for useful discussions and ongoing collaboration. He also gratefully acknowledges the support of the NITheP of South Africa. The support of a PRIN-COFIN 2007B58EAB grant is acknowledged. The research of A.S. was supported by the Spanish government through grant No. FIS2007-60977, partially financed by FEDER funds. AM would like to acknowledge the financial support of the MSMT of the Czech Republic under Project No. LC512 and the GAAS of the Czech Republic (Grant No. IAA400720710).

Appendix A Ruelle’s stability criterion

Let us consider the 1D PSW model characterized by and . Let us call the integer part of , i.e., . According to Ruelle’s criterion, a sufficient condition of thermodynamic stability is Fisher66 (); Ruelle69 ()


for all configurations , where is a fixed bound.

Given the number of particles , we want to obtain the configuration with the minimum potential energy . Without loss of generality we can see any given configuration as a set of clusters , each cluster being made of overlapping particles (i.e., any pair of particles of a given cluster are separated a distance smaller than ). In Ref. Santos08, we proved that for a fixed value of the minimum energy corresponds to , all the particles of each cluster being located at the same point and the centers of two adjacent clusters being separated a distance . Therefore, we can restrict ourselves to this class of ordered configurations and use as the variational variable.

The repulsive contribution to the potential energy is


To compute the attractive contribution we need to take into account that all the particles of a given cluster interact attractively with the particles of the nearest clusters. The total number of pairs of interacting clusters are . Therefore,


The total potential energy is


We then see that the value that minimizes is


This value is only meaningful if . Otherwise, . In summary, the absolute minimum value of is


Therefore, if the potential energy is bounded from below by with and thus the system is stable in the thermodynamic limit. On the other hand, if there exist configurations that violates Ruelle’s criterion and so the thermodynamic stability of the system is not guaranteed.

Appendix B Density expansion of within the HPA

Starting from Eq. (7) and for , the Fourier transform can be expanded in power series as


Upon inverse Fourier transform one then has




Equation (8) can be rewritten as




The origin () is a regular point of and hence of [but not of each separate term in Eq. (41)], so we can choose to save the point in Eq. (39) either from above or from below. Here we do it from above with the result


where the path in the complex plane goes from to and closes itself on the upper plane if and in the lower one if . In Eq. (39) we then find


It is interesting to note that if . Thus Eq. (38) can be rewritten as


where is the integer part of .

Figure 12: Plot of the radius of convergence and of the maximum density versus for and .

The radius of convergence of the series (44) depends on temperature and can be obtained by the same arguments as in the PS case.Malijevsky06 () From the denominator of Eq. (7), it follows that the series converges provided that , where


Here denotes the absolute maximum value of . From Eq. (8) if that maximum corresponds to , i.e., , and so . On the other hand, if the maximum value takes place at and so must be obtained numerically. For sufficiently large values of one has , so that and this coincides with the maximum physical density (see below). Figure 12 shows as a function of for two values of . In the PS limit () one has . As the strength of the attractive part of the potential (measured by the product ) increases, the radius of convergence first grows, reaches a maximum, and then decays.

Except when is so large that , the maximum corresponds to a negative value of and so the singularity responsible for the radius of convergence is located on the negative real axis. Therefore, is still well defined beyond the radius of convergence, i.e., for . On the other hand, analogously to the PS case,AS04 (); GK77 () the HPA for the PSW fluid becomes unphysical, at a given temperature, for densities larger than a certain value given by the condition


where is the absolute maximum value of . Since , it is obvious that . For sufficiently large values of (actually, for temperatures below the critical value defined in Section VII) one has , and so . In that case, the line of maximum density becomes a spinodal line, as discussed in Section VII. Figure 12 also includes a plot of as a function of for the same two values of . Note the kink of the curve for at the critical point (