Liquid-Liquid Phase Transitions for Soft-Core Attractive Potentials

Using event driven molecular dynamics simulations, we study a three dimensional one-component system of spherical particles interacting via a discontinuous potential combining a repulsive square soft core and an attractive square well. In the case of a narrow attractive well, it has been shown that this potential has two metastable gas-liquid critical points. Here we systematically investigate how the changes of the parameters of this potential affect the phase diagram of the system. We find a broad range of potential parameters for which the system has both a gas-liquid critical point and a liquid-liquid critical point. For the liquid-gas critical point we find that the derivatives of the critical temperature and pressure, with respect to the parameters of the potential, have the same signs: they are positive for increasing width of the attractive well and negative for increasing width and repulsive energy of the soft core. This result resembles the behavior of the liquid-gas critical point for standard liquids. In contrast, for the liquid-liquid critical point the critical pressure decreases as the critical temperature increases. As a consequence, the liquid-liquid critical point exists at positive pressures only in a finite range of parameters. We present a modified van der Waals equation which qualitatively reproduces the behavior of both critical points within some range of parameters, and give us insight on the mechanisms ruling the dependence of the two critical points on the potential's parameters. The soft core potential studied here resembles model potentials used for colloids, proteins, and potentials that have been related to liquid metals, raising an interesting possibility that a liquid-liquid phase transition may be present in some systems where it has not yet been observed.


I. INTRODUCTION
The discovery and investigation of liquid-liquid phase transitions in a one-component system is of current interest, since recent experiments for phosphorus [1,2] show a firstorder phase transition between two stable liquids in the experimentally accessible region of the phase diagram. A liquid-liquid phase transition, ending in a critical point was initially proposed to explain the anomalous behavior of network-forming liquids such as H 2 O [3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18]. In particular, the density anomaly, consisting in the expansion under isobaric cooling of these systems, has been related to the possible existence of a phase transition between low density liquid (LDL) and high density liquid (HDL). Simulation results and experimental studies of water predict a LDL-HDL phase transition in an experimentally inaccessible region of the phase diagram [9,12,15,19,20]. Computer simulations of realistic models of carbon [21], phosphorus [22], SiO 2 [23], and Si [24,25] strongly suggest the existence of first-order LDL-HDL phase transitions in these substances. Recently the step changes of the viscosity of liquid metal, such as Co, have been theoretically interpreted as evidence of liquid-liquid phase transitions [26].
The presence of the first-order phase transitions in solids and solid-solid critical points, determined experimentally [27] and with simulations [28,29,30,31], have suggested the possibility of the existence of liquid-liquid critical points and polymorphism in the amorphous state [32,33,34]. It has been proposed that systems with solid polymorphism may exhibit several liquid phases with local structures similar to the local structures of various crystals. Experimental evidence of sharp structural transitions between liquid polymorphs of Se, S, Bi, P, I 2 , Sn, Sb, As 2 Se 3 , As 2 S 3 , and Mg 3 Bi 2 are consistent with phase diagrams with firstorder liquid-liquid phase transitions [33,35], analogous to the liquid-liquid phase transition seen in rare earth aluminate liquids [36,37].
These results call for a general interpretation of the basic mechanisms underlying the liquid-liquid phase transition. Here we aim to delineate the conditions ruling the accessibility of the two liquid phases. A first step in this direction was taken in Refs. [38,39], where we have shown that a specific isotropic soft-core attractive potential, for a one-component system, has a phase diagram with LDL-HDL phase transition, with two fluid-fluid critical points, and with no density anomaly.
Here we extend this analysis by varying the parameters of this potential (Fig. 1). We find that, for a wide range of parameters, this potential has a phase diagram with a liquid-liquid critical point, and we show how the phase diagram depends on the parameters. We develop a modified van der Waals equation (MVDWE) able to describe the behavior of the two critical points as a function of the potential parameters, elucidating a mechanism for the liquid-liquid phase transition and the conditions under which the liquid-liquid critical point occurs at positive pressure.
In Sec. II we introduce the isotropic soft core potential; in Sec. III we describe the two different molecular dynamics techniques we use; in Sec. IV we present our results for different combinations of parameters that give rise to a liquid-liquid phase transition ending in a liquid-liquid critical point; in Sec. V we construct a modified van der Waals equation which can qualitatively reproduce the behavior of the two critical points; in Sec. VI we discuss the role of potential parameters in changing the position of the critical points; in Sec. VII we summarize our results; in Appendix A we present our simulation results for a simple square well potential.

II. THE ISOTROPIC SOFT-CORE ATTRACTIVE POTENTIAL
For attractive potentials with a sufficiently broad interaction distance, the phase diagram has a first-order gas-liquid transition ending in a gas-liquid critical point, and a first-order liquid-solid phase transition [40]. When the attractive range is small, the liquid phase and the gas-liquid critical point are metastable with respect to the solid phase [41,42,43,44,45].
For a strictly repulsive soft-core potential, simulations show a phase diagram with a first-order gas-solid phase transition and a first-order phase transition between two solids of different densities, but with the same structural symmetry, ending in a solid-solid critical point [28,29,30,31]. Recent theoretical work has suggested that systems with a broad soft-core potential have a fluid-fluid phase transition and liquid anomalies [46], or give rise to stripe phases in two dimensions [47].
We have shown in Ref. [38] that the combination of a repulsive soft core with an attractive well is sufficient to give rise to a phase diagram with two liquid phases. This simple isotropic model potential is similar to those used in the seminal work of Stell and Hemmer [48], who studied a soft-core potential in one dimension (1D). Similar potentials were studied in 2D and 3D showing phase diagrams with a possible liquid-liquid critical point [49,50].
The 3D isotropic potential we consider (Fig. 1) has a hard core (infinite repulsion) at distance a, a repulsive soft core of width w R and energy U R > 0, and an attractive square well of width w A and energy −U A < 0 [38,39]. The potential has three parameters: w R /a, w A /a, and U R /U A , where a and U A have been chosen as units of length and energy, respectively. Though this potential is discontinuous, it is similar to model potentials for complex fluids, such as colloids, protein solutions, star polymers [44,51,52,53,54,55,56], and resembles pair potentials proposed for water [52], or that have been related to liquid metals under specific conditions [57,58,59].
This potential with parameters w R /a = 1.0, w A /a = 0.2, and U R /U A = 0.5 has a phase diagram with gas-LDL and gas-HDL first-order phase transitions, each ending in a critical point in the supercooled fluid region [38]. Both liquid phases are metastable with respect to a single crystal phase and no density anomaly is observed [39].
In this paper we present systematical MD studies of the phase diagrams for this potential (Fig. 1). By varying the parameters of the potential, w A /a, w R /a, and U R /U A , we relate the attractive and repulsive components of the potential to the appearance and stability of the liquid-liquid phase transition and critical points.

III. MOLECULAR DYNAMICS SIMULATIONS
We perform MD simulations of N = 850 particles of unit mass m at constant volume V , and constant temperature T , interacting via the potential described above (Fig. 1). The details of the event-driven MD we use are presented in Refs. [38,39]. We measure time in units of aU −1/2 A m 1/2 and record potential energy and pressure every ∆t = 100 time units. To understand the effect of each parameter on the phase diagram of our system, we simulate sixteen sets of potential parameters (Table I). After a preliminary screening, we choose to study the region of parameter space where the low-density, gas-liquid critical point C 1 always has a critical temperature above that of the high-density critical point C 2 . Therefore, while in Refs. [38,39] C 2 is a gas-HDL critical point, here C 2 is a LDL-HDL critical point. As shown in in Refs. [38,39], C 2 can lie in the supercooled metastable phase, close to the line of homogeneous nucleation, as in water or silica [14,19,25]. We make certain that all our calculations are performed before the onset of crystallization, as discussed in Ref. [39]. The description of the crystal phases goes beyond the goal of this work. To optimize our analysis we use two different MD methods.

A. Isothermic method
The first method is a straightforward calculation of the phase diagram's state points. For each state point with given ρ = N/V and T , we perform typically ten independent simulations of t ≈ 2 × 10 3 time units. We estimate the error in pressure measurements from the standard deviation of the ten averaged values computed for each independent simulation. The state points along the isotherms are approximated by a two-variable polynomial P (ρ, T ) = iκ a iκ ρ i T κ obtained by the least squared fit of all the state points in the vicinity of the critical point. This fitting implies mean field critical exponents [60] and may produce incorrect results in the close vicinity of the critical point. However, this method helps us fit the state points, known with statistical errors, by approximate polynomial isotherms and thus obtain the approximate position of the critical point.
The coexistence curves are calculated using Maxwell's equal area construction and spinodal line is estimated by locating the maxima and minima of the isotherms. After calculating the state points, isotherms, coexistence curves, and spinodal lines, we estimate the critical pressure, temperature, and density for C 1 and C 2 (P C 1 , T C 1 , ρ C 1 , P C 2 , T C 2 , and ρ C 2 , respectively) as the point where coexistence and spinodal curves meet, coinciding at their maxima. We apply this method to six sets of potential parameters (ii, iii, iv, xi, xii and xiv in Table I).
The results are presented in Figs. 2-7 in the pressure-density (P -ρ) phase diagrams. The estimates of the critical points are presented in Table II.

B. Isochoric method
The isothermic method gives us fairly complete information about the details of the phase diagrams, but requires much computation to calculate enough state points for accurate isotherms. Thus, in order to find the positions of critical points for a wide range of potential parameters, we adopt a faster but less accurate MD method. For sets of parameters close to the sets of parameters studied with the isothermic method, we estimate the location of the spinodal line by evaluating the intersections of isochores in the P − T plane. We first equilibrate several configurations at a high initial temperature k B T I /U A = 2.0 for several values of density above and below the densities where we expect to find ρ C 1 and ρ C 2 . At constant density, the system is slowly cooled down from T I to a final temperature k B T F /U A = 0.1 during a simulation time of 10 4 time units [61].
The average values of T and P are recorded each 100 time units, which is comparable to the equilibration time of the system for k B T /U A > 0.5. As the temperature decreases, the equilibration time increases and the method becomes less reliable. Thus, we use this method to estimate pressure and potential energy for k B T /U A > 0.5.
The error bars of each measurement are of the order of the non-monotonic jumps of the isochores (see Fig. 8 inset). The intersection is determined by fitting isochores with smooth polynomial fits. The best results can be achieved by quadratic fits in the temperature range including the region of possible isochore crossing extending from 0.9T C to 1.5T C , so that the tentative critical temperature T C is inside this interval.
Since at the spinodal line (∂P/∂ρ) T = 0, two isochores with two close values of density must intersect in the vicinity of the spinodal line. By definition, the critical point corresponds to the maximum temperature on the spinodal. Therefore, the critical pressure and temperature can be evaluated by estimating the pressure corresponding to the maximum temperature at which isochores intersect. The critical density can be estimated as (ρ 1 + ρ 2 )/2, where ρ 1 and ρ 2 are the densities of the two isochores intersecting at the highest temperature (Fig. 8). The critical point values estimated with this method are presented in Table III.
This approximate method is allowed as long as we use it to estimate the critical points of potentials with sets of parameters close to those for which we have done a detailed study using the isothermic method. We apply the isochoric method to 16 sets of the potential parameters. The comparison of the two methods (Tables II-III) shows that the resulting estimates of critical P , T , and ρ of C 1 and C 2 are consistent.

IV. PHASE DIAGRAM RESULTS
Our results in Figs. 2-7 clearly show that the phase diagram strongly depends on the potential parameters. For example, phase diagrams in Figs. 2-4 have fluid phases (gas, LDL, and HDL) at positive pressures, while for the phase diagrams in Figs. 5-7, the high-density critical point appears at negative pressures, i.e. in the region of stretched fluid.
To investigate how the position of critical points depends on the potential parameters, we vary one of the three parameters w A /a, w R /a, U R /a at a time, keeping the other two constant. The behavior of T , P , and ρ for C 1 and C 2 ( Fig. 9 and Table IV) are presented in the following.
A. Effect of the square-well width w A By keeping w R /a = 0.5 and U R /U A = 2.0 constant, we find ( Fig. 9a-c, Fig. 10) that by increasing well width w A , ρ C 1 is almost unaffected, while ρ C 2 decreases, T C 1 and T C 2 increase, P C 1 increases, while P C 2 decreases. For w A /a > 0.7 the LDL-HDL critical point C 2 occurs at negative pressures, as in Fig. 5. Hence, C 2 lies in the stretched fluid region and, therefore, it is metastable. In order to have a stable LDL-HDL critical point, the attractive distance w A /a must be sufficiently narrow, so that C 2 occurs at positive pressures. A too narrow well, however, enhances crystallization [39,41,42,43,44,45] so that the highdensity critical point shifts below the line of spontaneous crystallization, becoming difficult to observe. Thus the liquid-liquid critical point is observable in our MD simulations only for intermediate values of w A /a.

B. Effect of the shoulder width w R
Increasing the width of the repulsive interaction w R , while keeping w A /a = 0.7 and U R /U A = 2.0 constant, we find ( Fig. 9d-f, Fig. 11) that both ρ C 1 and ρ C 2 decrease, T C 1 decreases, while T C 2 increases, and both P C 1 and P C 2 decrease. For w R /a < 0.4 the dynamics of the system in the vicinity of the expected high density critical temperature become too slow and the equilibration time becomes too long, with respect to our simulation time, to measure the equilibrium state points with sufficient accuracy. Furthermore, as expected for decreasing w R , T C 2 approaches T = 0 ( Fig. 9e), suggesting that C 2 disappears for w R /a = 0. At w R /a > 1.0 the system spontaneously crystallizes at high density without showing a second critical point C 2 . Hence, the width of the shoulder w R /a must be of an intermediate value for C 2 to be observed above the lines of spontaneous crystallization and outside the region of very slow dynamics, at least for our choice of w A and U R .

C. Effect of the shoulder height U R
For w R /a = 0.5 and w A /a = 0.9, we increase the repulsive energy U R and find ( Fig. 9gi, Fig. 12) that for increasing U R , ρ C 1 decreases, while ρ C 2 is almost unaffected, both T C 1 and T C 2 decrease, P C 1 decreases, while P C 2 rapidly increases. For U R /U A < 2.0 the highdensity phase transition occurs at very low negative pressures and the fluid phases are highly metastable. For U R /U A > 4.0 the diffusion in the system in the vicinity of the high density critical point becomes markedly slow, due to the soft core becoming less penetrable and assuming the role of an effective hard core. Therefore, an intermediate repulsive energy is needed to observe C 2 in our MD simulations.

V. MODIFIED VAN DER WAALS EQUATION
To rationalize the dependence of the temperature, pressure and density of the two critical points on the potential's parameters, we develope a simple mean field theory that gives rise to a modified van der Waals equation (MVDWE), that has the same form of the standard van der Waals equation (see Appendix A), but with an excluded volume B(ρ, T ) depending on the density and temperature of the state point and increasing with w R /a, and with a strength of attraction A that increases with w A /a and decreases with U R /U A . It should be pointed out that a different modification of the van der Waals equation [62] also gives rise to the high density critical point. In contrast with our work, Ref. [62] is particularly suitable for density dependent potentials since it assumes a constant excluded volume B and a density dependent attractive term A(ρ).
For a system with a hard core and a soft core, one can assume that the effective excluded volume B(ρ, T ) changes with temperature and density [63]. Indeed, at low densities and low temperatures, particle cannot penetrate into the soft core so B(ρ, T ) ≈ B 2 where B 2 = 2π(a + w R ) 3 /3 is the excluded volume associated with the soft core. In contrast, for high densities and high temperatures, particles easily penetrate into the soft core and B(ρ, T ) ≈ B 1 , where B 1 = 2πa 3 /3 is the excluded volume associated with the hard core. More specifically, B(ρ, T ) must be an analytical function of its parameters such that and lim T →0 from which it follows that B(ρ, T ) < 1/ρ for any ρ and T > 0. Since in any case van der Waals equation can give us only qualitative agreement with reality, we can select any model function B(ρ, T ) which satisfies the above conditions. Never the less, it is desirable to select B(ρ, T ) in such a way that it will describe the behavior of some physical system for which the analytical solution can be found. One dimensional system of particles with a pair potential provides such a solution. Applying the Takahashi method [50,64], we obtain the Gibbs potential where N 1 is the number of particles, T is temperature, P 1 is pressure of the one-dimensional system and Accordingly V 1 = ∂G/∂P 1 and S 1 = −∂G/∂T are the volume and entropy of the onedimensional system, and U 1 = G − P 1 V 1 + T S 1 is the potential energy for the one dimensional system. The fraction of the soft cores f (ρ, T ) penetrated by the particles is is the fraction of the soft cores penetrated by the particles in the high-temperature limit in which soft cores play no role. It can be computed assuming a Poisson distribution of interparticle distances: . The probability that the soft core does not reflect the neighboring particle is equal to the fraction of these two quantities f /f ∞ < 1. In this case, the excluded volume is equal to B 1 . In the opposite case with probability 1 − f /f ∞ , the excluded volume is equal to B 2 . Hence, the effective excluded volume where, and P 1 must be found from the equation Figure 13a illustrates the behavior of B(ρ, T ) for a particular set of parameters. It is clear that B(ρ, T ) satisfies all the physical conditions we impose on the effective excluded volume. The modified van der Waals equation (1) has two critical points: one for low density ρ << 1/B 2 and another for high density ρ ≈ 1/B 2 , whose positions on the phase diagram of the dimensionless variablesT = k B T /U R ,P = B 1 P/U R , andρ = B 1 ρ depend on the dimensionless parameters of the MDVWE: B 2 /B 1 and A/(U R B 1 ). Figure 13b shows a P − T diagram with two critical points C 1 , C 2 , for a particular set of parameters, for which the position of the critical points are similar to the positions found in our simulations, i.e.
Now we can relate the parameters of the Eq. 1 to the potential parameters used in our simulations. The parameters B 1 and B 2 are increasing functions of the hard core diameter a and the shoulder width w R , respectively. The parameter U R has an identical meaning in MVDWE and in simulations. The strength of attraction A is an increasing function of w A and a decreasing function of U R . Indeed, according to the formula of the second virial coefficient v 2 for our potential, we have [65] and v R are positive quantities with the dimension of a volume depending on a, w R , and can be rewritten in the form of the van der Waals equation From the equations above we can derive the functional relation for v A and v R in the By using these relations it is possible to see that in general A is an increasing function of w A and a decreasing function of U R . The derivative ∂A/∂w R may have different sign depending on other parameters. Although at finite T these relations could be valid only to the leading order, it is reasonable to assume that A increases with w A and decreases with U R at any T .
However, to simplify our qualitative study of the MVDWE, we assume the parameters A, B 2 and U R are independent. By varying these parameters one at a time and by relating B 2 to w R , and A only to w A , we found that the MVDWE predicts that the derivatives of the low-density critical point values, T C 1 , P C 1 , ρ C 1 , with respect to each of the parameters of MVDWE, have the same sign and this sign is negative for B 2 (w R ) and U R , which increase repulsion, and this sign is positive for A(w A ), which increases attraction. These results are consistent with the MD results for the low-density critical point. For the high-density critical point values the MVDWE predicts for some parameters non-monotonic behaviors and we find that there are regions of parameters where the critical values as a function of A, B and U R have the same qualitative behaviors as those found in the simulations as a function of w A , w R and U R , respectively (Fig. 14). Specifically, we observe monotonic behaviors of ρ C 2 (B 2 ), T C 2 (A), T C 2 (B 2 ), and P C 2 (U R ) which qualitatively coincide with the corresponding behaviors in simulations (Figs 9 and 14 d,b,e, and i). We find non-monotonic behaviors of ρ C 2 (A), ρ C 2 (U R ), T C 2 (U R ), P C 2 (A), and P C 2 (B 2 ), which qualitatively coincide with the corresponding behaviors in simulations in certain range of parameters (Figs 9 and 14 a,g,h,c, and f).
These observations indicate that the behavior of the critical points in simulation may also become non-monotonic in the range of parameters that we do not explore. For example, T C 2 (U R ) may start to increase for large U R /U A > 4 and small w R /a < 0.5. Another interesting prediction of the MVDWE is that for large 2) = 1.7 the high-density critical temperature becomes larger than the low-density critical temperature as in simulations of Refs. [38,39], for which the repulsive shoulder w R /a = 1 was much wider than the attractive well w A /a = 0.2. Also, MVDWE predicts the existence of the third, very high density critical point for large B 2 /B 1 and large AB 1 /U R , which was recently observed in simulations with a wide soft core [66].

VI. ROLE OF POTENTIAL PARAMETERS
In the following we will present the comparison between the MD results and MVDWE predictions.
A. The low density critical point First we note that at low densities, corresponding to the critical point C 1 , and at sufficiently low temperatures, particles do not penetrate into the repulsive region, r < a + w R . Therefore, we can assume that, at low enough temperatures and densities, the system is interacting via an effective potential given by a simple square-well with hard core a + w R , an attractive well of relative width w A /(a + w R ) and attractive energy U A .
Indeed, for increasing width of the attractive well w A , ρ C 1 is roughly constant (Fig. 9a) and T C 1 and P C 1 increase (Fig. 9b and c). This behavior is consistent with the predictions of the standard van der Waals theory for the gas-liquid critical point for a square-well potential (see Appendix A), that yields Eqs. (A7-A9). This result supports the idea that the effect of the soft core is negligible at low densities. The MVDWE also predicts strong increase of P C 1 and T C 1 with the strength of attraction A, which increases with w A . For ρ C 1 , the MVDWE predicts a weak increase, which can be observed in Fig. 9a for larger w A .
For increasing width of the repulsive shoulder w R , ρ C 1 decreases and saturates (Fig. 9d) and T C 1 and P C 1 decrease (Fig. 9e and f). As a consequence of the above considerations, the increase of w R corresponds to a decrease of this effective attractive parameter. Accordingly, T C 1 and P C 1 display the behavior predicted for the decrease of the attractive width in van der Waals theory for the square-well potential (Eqs. A8-A9). Moreover, the behavior of ρ C 1 is consistent with Eq. (A7), which predicts ρ C 1 ∼ 1/a 3 . Indeed, in our case, the hard core is replaced by the effective hard core, hence a 3 ρ C 1 ∼ a 3 /(a + w R ) 3 , which is a decreasing function of w R . The calculations using MVDWE completely confirms these predictions by showing that T C 1 , P C 1 , and ρ C 1 all decrease with increasing w R (or B 2 = 2π(a + w R ) 3 /3 as in Sec.V).
For increasing repulsive energy U R , the behaviors of ρ C 1 , T C 1 and P C 1 are the same as those observed for increasing w R (Fig. 9g-i). This can be understood by considering that the increase of U R effectively decreases the penetrability of interparticle distances r < a + w R . The soft core becomes an effective hard core as described above. In particular, the saturation of ρ C 1 , already observed in Fig. 9d, is now more evident and an analogous behavior is now also seen for T C 1 and P C 1 . This result shows that for a high-enough repulsive energy U R and low-enough T , the soft-core potential is equivalent to a hard-core potential, for which there is no dependence of the critical point on U R . Again, the MVDWE agrees with these predictions by showing that T C 1 , P C 1 , and ρ C 1 decrease with increasing U R .

B. The high density critical point
For increasing w A , the critical point density ρ C 2 decreases, the temperature T C 2 increases, and the pressure P C 2 decreases (Fig. 9a-c). This finding is in agreement with MVDWE predictions for a wide range of parameters. (Fig. 14a-c). The behavior of T C 2 is consistent with the idea that the increase of the attractive distance increases the overall attractive strength of the potential, allowing more particles to fit within the attractive interaction range. As a consequence, the system enters the low-energy and high-density fluid phase at a higher temperature, i.e. T C 2 increases. Hence, the increase of w A increases the average kinetic energy of particles at C 2 , favoring the overcoming of the soft-core shoulder at low pressures, and we can expect that P C 2 decreases. Moreover, the increase of w A decreases the number of elastic interparticle collisions at the soft-core distance, hence decreases their contribution to the virial expression for the pressure [67] (see Eq. (11) in Ref. [39]), decreasing the critical pressure P C 2 . Note that the behavior of P C 2 in this case is the opposite of the behavior for P C 1 (Fig. 9c). With the decrease of pressure, the density must also decrease, so we can conclude that ρ C 2 must decrease with increasing w A , in agreement with our simulation results.
For increasing w R , T C 2 increases and saturates, and both P C 2 and ρ C 2 decrease with a tendency toward saturation (Fig. 9d-f), in agreement with predictions of MVDWE for a wide range of parameters (Fig. 14d-f). This happens because the transition from the LDL to the HDL is characterized by the penetration of particles into the repulsive soft cores of their neighbors. Therefore the repulsive soft-core distance a + w R characterizes the typical distance between the particles at C 2 . Hence the increase of w R reduces the critical density ρ C 2 . The behavior of pressure follows the behavior of density, as in the case of w A , while the derivatives of pressure and temperature must have the opposite sign due to the same arguments as above.
For increasing U R , P C 2 increases, ρ C 2 slowly increases, and T C 2 decreases (Fig. 9g-i). The predictions of MVDWE coincide with the behavior of P C 2 and ρ C 2 in a wide range of parameters (Fig. 14g,i) . However, the theory apparently predicts an increasing T C 2 with U R , except for very small U R B 1 /A < 0.4 and large B 2 /B 1 > 1.5 (Fig. 14h). This discrepancy arises from the fact that, although it is physically clear that the attractive strength A is a decreasing function of U R , we find the explicit dependence of A on U R only in the limiting case T → ∞, while for finite T we assume them to be independent. Hence, we ignore that an increase of U R decreases A which induces, as shown in Sec.V, a decrease in T C 2 .
The behavior of P C 2 is easier to understand. Indeed, the pressure at which the repulsive shoulder can be overcome increases with U R , which is consistent with the increase of P C 2 with U R . This effect is expected to be more evident at high U R and to saturate for decreasing U R , which is consistent with our results. The critical density ρ C 2 increases with U R , for small values of U R /U A , as a consequence of the increase of P C 2 , and is practically independent of U R when the soft core plays the role of an effective hard core, i. e. for large enough U R /U A . The decrease of T C 2 with the increase of U R is more difficult to explain. Nevertheless, the same argument as in the case of w A and w R which predicts that the derivatives T C 2 and P C 2 must have the opposite signs may apply in this case as well. Finally, we note that increasing with U R and decreasing with w R , the behavior of P C 2 in three dimensions is consistent with its behavior in the one dimensional case, for which P C 2 /U A = (U R /U A + 1)/w R [68].

VII. DISCUSSION AND CONCLUSIONS
We have studied an isotropic attractive soft-core square potential in three dimensions that has a phase diagram with a gas-liquid critical point C 1 and a liquid-liquid critical point C 2 , separating a HDL and a LDL phases. We have investigated, with molecular dynamics simulations, how the critical density, temperature and pressure of the two critical points vary as a function of the three parameters of the potential, that are the repulsive energy U R /U A in units of the attractive energy U A , the repulsive width w R /a in units of the hard core a, and the attractive width w A /a in units of a.
Table IV and Fig. 9 show our results for the ρ, T , and P of C 1 and C 2 for varying parameters of the potential. To summarize, the behavior of C 1 is consistent with that of a system interacting via an effective square-well potential, with a hard core a + w R , a relative attractive well w A /(a + w R ) and attractive energy U A . The increase of U R /U A or of w R /a decreases the effective attractive strength and this effect saturates for large values of U R /U A .
This behavior is perfectly predicted by the simple mean field MVDWE. In MVDWE, as in the standard van der Waals equation, the relevant physical parameters are the excluded volume B and the stength of attraction A, but now we assume that B = B(ρ, T ) increases with w R /a and depends on the state point, while A increases with w A /a and decreases with The MVDWE predictions are consistent also with our results on C 2 . In general, this approach rationalizes why the increase of U R /U A and the decrease of w A /a have the same qualitative effect on the critical points, since they have the same effect of the attractive strength. Moreover, it rationalizes the effect of increasing w R via the increase of the excluded volume, hence the decrease of the critical densities. This decrease induces a decrease of the critical pressures P C 1 and P C 2 , as a consequence of the mechanical stability of the fluid phases.
While for the low-density critical point C 1 , the decrease of P C 1 occurs with the decrease of T C 1 , i.e. their derivative with respect to the potential parameters always have the same sign, for the high-density critical point C 2 the critical pressure and temperatures always have derivatives with the opposite sign. This behavior can be understood in terms of the number of elastic collisions with the soft core, which decreases as T C 2 increases, reducing the virial contribution to the pressure. At the same time, the increase of T C 2 reduces the pressure P C 2 necessary to overcome the repulsive soft core and enter into the HDL phase.
As a consequence, the high density critical point C 2 exists at positive pressure only in a finite region in the parameter space. Indeed, when the attraction is too strong, i.e. w A /a is too large or U R /U A is too small, the pressure P C 2 becomes negative. On the other hand, when the strength of attraction is too week, C 2 occurs in the deeply supercooled liquid phase, becoming difficult to observe as in the experimental situation of water or silica [14,19,25].
In conclusion, the behavior of both low-density and high-density critical points qualitatively obeys the mean field predictions of the modified van der Waals theory based on effective excluded volume which varies between the hard-core value for high temperature and low density and the soft-core value for low temperature and low density. The quantitative theory based on the thermodynamic perturbation approximations or various integral equation closures [69] is yet to be developed.
Confirming the results presented in Refs. [38,39] we do not find density anomaly. Our simulations show that density anomaly is unlikely to exist for the discontinuous double-step potentials shown in Fig. 1, in contrast to ramp potentials [49] and a Gaussian soft core potential [53].
Our results may be relevant for experiments on systems that can be described by an isotropic soft-core attractive potential and have no density anomaly, such as colloids, protein solutions or liquid metals. Indeed, our results show that in these systems the possibility of the existence a liquid-liquid phase transition will depend on the relative ratio between the attractive and the repulsive part. Here we recall the case of a square well attractive potential, where the only parameter is the width of the attractive well w A /a, since the hard core distance a and the attractive energy U A can be taken as units of distance and energy, respectively. In particular, we show how the gas-liquid critical point density ρ c , temperature T c , and pressure P c depend on w A /a.
Even in this simple case, the phase diagram has no exact analytical solution and one must rely on various approximations and numerical simulations [40,69,70,71,72]. Using MD simulations of N = 850 particles, we verify that the behavior of ρ c , T c , and P c , are approximately linear for a wide range of w A /a (Fig. 15, Table VI) [40]. The values of a 3 ρ c decreases for increasing w A /a and a 3 P c /U A and k B T c /U A increase with w A /a.
Except for density, these results are in agreement with the van der Waals theory (see, e.g., Ref. [65]). The equation of state in the van der Waals theory is given by where B = 2 3 a 3 π has the meaning of excluded volume per particle and A is a quantity, with the dimension of the product of energy and volume, characterizing the strength of attraction between particles. Therefore, A can be related to the product of U A and the volume of the attractive well, which is proportional to w A a 2 for small w A /a.
The position of a critical point must satisfy the equations For the van der Waals equation of state Eq. (A1), Eqs. (A2) and (A3) yield the coordinates of the critical point: Hence: which predict that k B T c /U A and a 3 P c /U A increase with w A /a, while ρ c a 3 does not depend on it.
TABLE I: Sets of parameters for the generic soft-core potential (Fig. 1) considered in this paper: w R /a and w A /a are the soft-core width and the attractive width, respectively, both in units of the hard-core distance, and U R /U A is the repulsive energy in units of the attractive energy. Sets i-vi have same w A and U R ; sets ii, vii-xii have same w R and U R ; sets xii-xvi have same w R and w A .  Tables II and  III).   TABLE II: Temperatures T C 1 and T C 2 , pressures P C 1 and P C 2 , and densities ρ C 1 and ρ C 2 , for the critical points C 1 and C 2 , respectively, computed by the isothermic method.
Set  TABLE III: Temperatures T C 1 and T C 2 , pressures P C 1 and P C 2 , and densities ρ C 1 and ρ C 2 for the critical points C 1 and C 2 , respectively, estimated by cooling the system at constant ρ (isochoric method) for the potential with the set of parameters in Table I. Set  IV: Summary of the effects on ρ C 1 , T C 1 , P C 1 , and ρ C 2 , T C 2 , P C 2 , from variation of parameters w A /a, w R /a, and U R , one at the time. The symbols ↑, ↓ and ≈ represent, respectively, an increase, a decrease and a small variation of a thermodynamic quantity as a consequence of the increase of the potential parameter.   Fig. 15 for k B T c /U A = p + qw A /a, and analogous linear expressions for a 3 P c /U A and a 3 ρ c , for the gas-liquid critical point.           Where not shown, errors are smaller than the symbol size. Lines are the linear fits T c , P c , and ρ c as functions of w A , with the parameters in Table VI.