Liquid-Liquid Phase Transition for an Attractive Isotropic Potential with Wide Repulsive Range

Recent experimental and theoretical results have shown the existence of a liquid-liquid phase transition in isotropic systems, such as biological solutions and colloids, whose interaction can be represented via an effective potential with a repulsive soft-core and an attractive part. We investigate how the phase diagram of a schematic general isotropic system, interacting via a soft-core squared attractive potential, changes by varying the parameters of the potential. It has been shown that this potential has a phase diagram with a liquid-liquid phase transition in addition to the standard gas-liquid phase transition and that, for a short-range soft-core, the phase diagram resulting from molecular dynamics simulations can be interpreted through a modified van der Waals equation. Here we consider the case of soft-core ranges comparable with or larger than the hard-core diameter. Because an analysis using molecular dynamics simulations of such systems or potentials is too time-demanding, we adopt an integral equation approach in the hypernetted-chain approximation. Thus we can estimate how the temperature and density of both critical points depend on the potential's parameters for large soft-core ranges. The present results confirm and extend our previous analysis, showing that this potential has two fluid-fluid critical points that are well separated in temperature and in density only if there is a balance between the attractive and repulsive part of the potential. We find that for large soft-core ranges our results satisfy a simple relation between the potential's parameters.

molecular dynamics simulations of such systems or potentials is too timedemanding, we adopt an integral equation approach in the hypernetted-chain approximation. Thus we can estimate how the temperature and density of both critical points depend on the potential's parameters for large soft-core ranges. The present results confirm and extend our previous analysis, showing that this potential has two fluid-fluid critical points that are well separated in temperature and in density only if there is a balance between the attractive and repulsive part of the potential. We find that for large soft-core ranges our results satisfy a simple relation between the potential's parameters. The phase diagram of a typical monatomic substance is comprised of solid and fluid phases, with the fluid phase separating below the critical point into gas and liquid phases.
The prototype of such substances are simple (i.e., argon-like) fluids. Interparticle interactions in these systems can be appropriately described by the well-known Lennard-Jones potential. Other simple models-such as those described by the hard-sphere square-well potential or by the hard-sphere-Yukawa potential-exhibit similar phase diagrams. All these potentials consist in a short-range harshly repulsive core plus a longer-ranged attraction.
New insights into the relationship between phase diagrams and interparticle interaction emerged recently from the finding that when the range of the attractive component is sufficiently small, the liquid phase and the gas-liquid critical point become metastable with respect to crystallization [1][2][3][4][5] . Shouldered potentials, i.e., potentials with a hard core and a finite repulsive shoulder, exhibit even more exotic phase diagrams. Simulations and theories showed that such potentials may give origin to non-trivial phase behaviors, such as isostructural solid-solid transitions and liquid-liquid transitions [6][7][8][9][10][11][12][13][14][15] . The key to this complex phase behavior resides in the peculiar penetrability of the repulsive core, a feature that gives rise to a density-dependent effective interaction.
The possible existence of a liquid-liquid phase transition for single-component systems in particular has received considerable attention in recent years. A direct evidence of this phenomenon has been observed on the experimental side in liquid phosphorous 16,17 . Experimental data consistent with a liquid-liquid phase transition have been presented for other single-component systems such as water [18][19][20] , silica 21,22 , carbon 23 , selenium 24 , and cobalt 25 , among others. A recent theory has shown that a phase transition occurring in a solution of DNA-coated colloids can be understood in terms of the liquid-liquid phase transition 26 .
We have recently shown through molecular dynamics (MD) simulations 10 that a system of particles interacting through an isotropic potential with an attractive well and a repulsive component consisting of a hard core plus a finite shoulder may possess a high-density liquid phase and a low-density liquid phase. Potentials with such characteristics were used to model interactions in a variety of systems including liquid metals, metallic mixtures, electrolytes, and colloids, as well as anomalous liquids, like water and silica [35][36][37][38][39][40][41][42][43][44] .
In spite of the simplicity of the model, the physical mechanism that causes the liquidliquid transition for a potential with a hard core plus a repulsive shoulder and an attractive well is not easy to assess since it arises from an interplay of the different components of the interparticle interaction. To disentangle the role of each component it is necessary to investigate the dependence of the phase stability of the fluid on the potential parameters.
This task was undertaken in Ref. 15 , where the results of MD calculations performed for several sets of parameters were presented. The resulting behavior of the critical points was interpreted through a modified van der Waals equation (MVDWE) 15 , a mean field approach assuming that the effect of the repulsive shoulder at different densities ρ and temperatures T can be taken into account by an effective excluded volume depending on both ρ and T .
However, the analysis was limited to cases where both the soft-core range and the attractive range are smaller than the hard-core range a and the total interaction range does not exceed 2.6a. Indeed, MD becomes quite time-demanding for larger interaction ranges. Nevertheless, there are cases such as biological solutions and colloids where the soft-core range could be as large as the hard-core or even larger 45 . For this reason and to gain a better understanding of the role played by each component of a soft-core's attractive potentials, we will now explore how the phase diagram changes when the soft-core range exceeds the hard-core diameter.
For this we use the hypernetted-chain (HNC) integral equation for the radial distribution function 46 . This approach can be considered in many respects as an intermediate between MD and MVDWE. In fact, the HNC equation is far less time-consuming than a simulation, but is by no means as accurate. On the other hand, being a microscopic theory, it is not based on the assumption typical among mean field approaches such as the MVDWE that particles experience a uniform attractive potential. Hence it is intrinsically more accurate.
Since the HNC equation provides a fast way of estimating the position of the critical points, we perform an extensive investigation in the space of the potential parameters, considering an extremely ample number of combinations. Thus we are able to frame previous results into a wider perspective and obtain a better understanding of the physical mechanism leading to a liquid-liquid transition in one-component fluids.
Our results show that the high-density critical point can be found only when there is a balance between the attractive part and the repulsive part of the potential. In Ref. 15 this balance was expressed through the mean field strength of attraction, a parameter proportional to the attractive range w A /a ( Fig. 1) and inversely proportional to the repulsive energy U R , for fixed attractive energy U A . Here we find an approximated relation between U R /U A and w A /w R that is well-verified for large soft-core ranges and quantifies the ideal balance between the repulsive and the attractive components of the potential more effectively. Our results show that the liquid-liquid phase transition could be found in systems with small repulsion if the attraction is small as well, with U R /U A ∝ w A /w R , and in systems with wide repulsion, with U R /U A ∝ 3w A /w R . Typical systems with these characteristics could be colloids, where the effective repulsion and attraction can be regulated 26,47 .

II. THE ATTRACTIVE SOFT-CORE POTENTIAL
A soft-core potential with an attractive interaction at large distances was first proposed by Hemmer and Stell 35 to understand the possibility of the solid-solid critical point in materials such as Ce and Cs and was studied through an exact analysis in 1D. Other soft-core potentials with an attractive well were proposed and studied with approximate methods or with numerical simulations in 2D to rationalize the properties of liquid metals, alloys, electrolytes, colloids, and the water anomalies [36][37][38][39][40][41][42][43]10,11,44,12,14 .
The peculiarity of such potentials is the presence of two repulsive length scales. This feature is typical of systems with core-corona architecture such as, for example, star polymers. However, isotropic soft-core potentials have also been proposed as effective potentials resulting from an average over the angular degrees of freedom for systems where the distance of the minimum approach between particles depends on their relative orientation. Thus, in some respects, they have been considered 41,43,44 as simplified models of complex anisotropic interactions, such as those resulting from the hydrogen bonding between water molecules.
The model potential considered in this paper is similar to that investigated in Refs. 10,11,14,15 and is an isotropic pair potential with two characteristic short-range repulsive distances: one associated with the hard-core exclusion between two particles and the second with a weak repulsion (soft-core), which can be overcome at large pressure. More precisely, our pair potential U(r) (Fig. 1) consists of a hard-core of radius a, a repulsive square shoulder of height U R extending from r = a to r = b, and an attractive component having the form of a square well of energy −U A < 0 extending from r = b to r = c (here r is the interparticle distance). Choosing a and U A as length and energy units respectively, this potential depends on three free parameters: the width of the soft-core w R /a ≡ (b − a)/a, the width of the attractive well w A /a ≡ (c − b)/a, and the soft-core energy U R /U A .
Our aim is to understand how the position of the critical points in the thermodynamic plane changes upon varying the parameter values. In Ref. 15 we investigated a number of cases with w R < a and w A < a through MD simulations and presented a mean field approach with an MVDWE to rationalize the results. However, the MD analysis of potentials with w R ≥ a would require very large computation times, therefore we study this case with integral equations in the HNC approximation, which represents a compromise between accuracy and economy.

III. THE HYPERNETTED CHAIN INTEGRAL EQUATION APPROACH
The spatial distribution of a system of particles may be conveniently described by the radial distribution function g(r) 46 , a quantity directly measurable by scattering experiments and related to the thermodynamic properties of the fluid. One of the theoretical approaches most used to calculate this function is represented by integral equations. These are based on the so-called Ornstein-Zernike (OZ) relation which relates the total pair correlation function h(r) ≡ g(r) − 1 to the direct correlation function c(r), which describes the contribution coming from the direct interaction between two particles at distance r. Both h(r) and c(r) are unknown functions, thus to solve the OZ relation, one needs another relation (closure) between these two functions. This is by necessity an approximate relation and it is obtained in the HNC equation 48 by neglecting the sum over a specific class of diagrams in the diagrammatic expansion of g(r) 46 . In principle, this approximation is expected to work better at lower ρ, where the direct correlation function c(r) is more relevant than the correlation propagated through the other particles.
The solution of the system formed by the OZ relation plus the HNC closure is obtained through a numerical iterative procedure which is stopped when the difference between two consecutive elements of the succession generated (in the space of the distribution functions) is smaller than some given value. Based on the standard iterative procedure, different algorithms can be used to improve the accuracy and rapidity of convergence of the HNC equation's numerical solution. However, independent of the algorithm used, there is a region in the ρ-T plane where no solution can be found, i.e., for any ρ, there is a T below which the numerical algorithm does not converge. This defines an instability line (IL) of the theory in the ρ-T plane. The nature of the locus of instabilities of the HNC equation and its relationship with the spinodal line of the fluid was the object of extensive investigation 49 which showed that though the isothermal compressibility increases considerably when the instability line is approached from above, there is no real divergence of this quantity. Thus, identifying the IL of the HNC equation with the spinodal line of the fluid, which is characterized by a diverging compressibility, is not possible. Another well-known deficiency of the HNC is its thermodynamical inconsistency. This can at least be partially removed through suitable modifications of the equation which, however, make the method considerably less rapid.
In spite of the above limitations, knowledge of the instability region of the HNC equation may allow us to estimate the topology of the region of spinodal decomposition of the fluid.
In particular, for the potential defined in Sect. II, with parameters w R /a = 1, w A /a = 0.2, and U R /U A = 0.5, it was found that the IL is qualitatively similar to the spinodal line calculated through MD calculations 14 . More precisely, the density and temperature of the low-density critical point estimated through the HNC equation are in satisfactory agreement with simulation results, while the density of the second-critical point is overestimated by the theory 14 . This is not surprising since the theory is an approximate one and becomes progressively less accurate as the density increases. However, the ability of the HNC equation When the shoulder height increases, the potential gradually changes from potential A to potential B starting from U R /U A = −1. In any intermediate configuration, the potential has a penetrable finite repulsive shoulder. It is also possible to go from potential A to B following a different "path", i.e., by increasing the hard core radius from a to b. In this case, however, the potential always consists of an impenetrable hard-core (of varying radius) plus an attractive square well.
The IL is shown in Fig. 2 at fixed shoulder and well widths (w R /a = 1, w A /a = 0.2) for several values of the shoulder height U R . In the two limit cases corresponding to potentials A and B, the IL exhibits a single maximum corresponding to a phase diagram with a single liquid-gas critical point, a well-known behavior for a fluid of hard spheres with an attractive well. The position of the critical points in the ρ, T plane is considerably different in the two cases. The critical point corresponding to potential B is at a lower temperature than that corresponding to potential A, due to the weaker attraction, i.e., shorter attractive range c−b of potential B with respect to the largest attractive range c − a of potential A. Furthermore, the critical density for potential B is smaller than that for potential A and rescales as the hard-core volume (a/b) 3 of the two potentials. We observe that this rescaling overshadows the shift of the critical point toward higher densities due to the decrease of the attraction range (e.g., see Appendix A in Ref. 15 ), unless b/a ≃ 1.
As U R increases, starting from U R /U A = −1 (potential A), the IL moves toward lower temperatures as a consequence of the overall reduction of the interparticle potential's attractive component. At the same time, the IL undergoes a topological change which eventually yields a line with two maxima (Fig. 2b). This peculiar topology of the IL becomes most evident for intermediate values of U R (0.4 ≤ U R /U A ≤ 0.6). As U R increases further, the second maximum disappears and again the shape of the IL becomes more and more similar to the shape typical of the hard-core square-well potential. Thus, when w R and w A are fixed, two maxima are observed in the IL only for a finite range U max R ≤ U R ≤ U min R . In this range of values, as U R increases, the density ρ 1 of the low-density maximum becomes smaller, while that of the other maximum ρ 2 slightly increases. The critical temperatures T 1 and T 2 respectively corresponding to these two maxima, decrease-this behavior being more evident for the second maximum. These results agree with the behavior found with MD simulations for the two critical points (reported in Fig. 9g, 9h of Ref. 15 ). Thus for increasing U R , the two maxima move away from each other both in density and temperature.
We now consider a fixed shoulder width (w R /a = 1) and several values of the well width (w A /a = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6). We calculate the IL for each of them, letting the height U R of the shoulder vary (Figs. 3a, 2, 3b, 4a, 4b and 5a). The values of the shoulder height U R for which two maxima are observed increase with w A . By increasing w A , the second maximum is less and less evident and for large values of w A /a (Figs. 4b and 5a) the second maximum is not observed for any U R . For small values of w A , the decrease of the attraction flattens the curve and the second maximum becomes difficult to observe (Fig. 3a).
For IL's with two maxima and the same w R /a = 1 and U R but different w A , both maxima move toward higher temperatures for increasing w A , due to the increased attraction.
Moreover, by increasing w A , ρ 2 becomes smaller while ρ 1 does not vary significantly. This behavior agrees with that predicted by MD simulations for the two critical points (shown in Fig. 9a,9b of Ref. 15 ).
In Fig. 5b we show the behavior of the IL for w R /a = 0.8 and w A /a = 0.2. Comparing these results with those shown in Fig. 2, we observe that for fixed w A and U R , as w R increases, ρ 1 and T 1 are almost constant and at variance with MD results, but ρ 2 decreases and T 2 increases in agreement with the MD simulations (see Fig. 9d,9e of Ref. 15 ).
We next consider a potential with a wider repulsive shoulder (w R /a = 1.5) and several values of the well width (w A /a = 0.5, 1.0, 1.5, 2.5). The behavior of the IL (Figs. 6 and 7) for varying U R is quite similar to that observed in the previous cases. For fixed w R and w A , the IL shows only two maxima in a finite range of values of U R ; these values increase with w A and, for large values of w A , the two-maxima topology is not observed regardless of the value of U R (Fig. 7b). The range of values of w A in which we observe two maxima is larger with respect to the case w R /a = 1.
For one particular set of parameters (w R /a = 1.5, w A /a = 0.5 and U R /U A = 0.8), it is possible to compare the results obtained using the HNC equation with the phase diagram calculated through a theoretical approach based on a thermodynamically consistent integral equation 11 . Once again it appears evident that the main flaw of the HNC equation is to overestimate the critical density of the second critical point.
However, a direct comparison of HNC results with those obtained through MD simula-tions can be disappointing. For some of the parameter sets investigated in Fig. 9 of Ref. 15 the IL shows only one maximum, while for others the two-maxima topology is barely observable. As an example, we show the IL corresponding to the parameters w R /a = 0.5, w A /a = 0.5 with 1.0 ≤ U R /U A ≤ 1.7 (Fig. 8). It was not possible to directly analyze the value U R /U A = 2 (considered in Ref. 15 ) since, in this case, the HNC cannot be solved at high densities before any considerable increase of the compressibility can be observed (in general, this occurs when the finite repulsion is considerably stronger than the attraction). The results obtained at slightly smaller values of U R show, however, a non-monotonic behavior of the IL, suggesting the presence of a liquid-liquid critical point.

V. DISCUSSION
The overall behavior of the IL's is synthesized in Fig. 9 which shows, for different values of the shoulder width (w R /a = 0.6, 0.8, 1.0, 1.5), the points in the w A , U R plane where two maxima are found by using the HNC approximation. We observe that both ranges of w A and U R where two maxima are observed increase with w R .
The general behavior of U R as a function of w A at constant w R can be rationalized by using the modified van der Waals approach (MVDWE) presented in Ref. 15 . First we approximate the interval of values of U min R ≤ U R ≤ U max R for each w R and w A in Fig. 9 with its middle point U * R = (U max R + U min R )/2 (Fig. 10). Next, we recall from Ref. 15 the relation between the potential's parameters and the strength of attraction A, a parameter increasing with w A /a and decreasing with U R /U A . In particular, as T → ∞ it is with and The relation can be rewritten as where the coefficients are the volumes, and the ratios of the surfaces and radii of the hard core (HC) and the soft core (SC), respectively, and all depend only on the parameter w R /a. Hence, at a fixed value of w R , the function U R (w A ) in Eq. (4) only has A as an unknown parameter.
This MVDWE prediction can be verified by using the HNC results. In Fig. 10, Eq. (4) is used to fit the values of U * R (w A ) resulting from HNC calculations for different values of w R , with A as the only fitting parameter. As expected from Eq. (4), when w A /a < 1, the leading order in U * R (w A ) is linear, while when w A /a > 1 (corresponding to larger w R ), the non-linear behavior is evident. Fig. 10 also shows that, by increasing w R , the coefficients of the third-degree polynomial in w A decrease as predicted by Eqs. (4) and (5).
Moreover, the fitting parameter A in Fig. 10 shows a non-monotonic behavior with w R . This is consistent with the MVDWE prediction in Ref. 15 that ∂A/∂w R may have different signs, depending on the other parameters. Therefore, Eqs. (4) and (5) give us a fair description of how the three parameters U R , w R , and w A are related to each other when the phase diagram has two critical points at positive pressure and finite temperature. However, Eqs. (4) and (5) do not help us understand why the phase diagram has an accessible liquidliquid critical point only for limited ranges of w A and U R , given a value of w R .
To gain some insight into this point, we observe that if we plot Eq.(4) with A = 0 and no fitting parameters (Fig. 11), we get a rough approximation of the calculated U * that becomes fair for the largest w R . This suggests that as a first approximation we can assume that A = 0 at least for w R /a > 1, which is consistent with the conclusion of Ref. 15 that in order to have two accessible critical points in the fluid phase, the attractive and repulsive part of the potential must compensate, i.e., Hence from Eqs. (2) and (3) we get the approximation First we observe that to get an accessible liquid-liquid critical point, U R /U A ∼ O(1) is the relevant case. Indeed, the case with U R /U A ≫ 1 at high-enough T and small-enough P corresponds to an effective attractive potential with no repulsive shoulder and a hard core at a distance of a + w R with no liquid-liquid phase transition, or with a liquid-liquid phase transition at vanishing T and very high P (see MD results in Fig.s 9h, 9i in Ref. 15 ). On the other hand, for U R /U A ≃ 0, Eq. (6) gives w A ≃ 0, leading to a simple hard-core potential with no attractive well. Hence we consider the case with U R /U A ∼ O(1).
Next, we observe that for increasing w R (w R → ∞), Eq. (6) goes to from which the condition U R /U A ≃ 3w A /w R follows for small w A /w R . This relation is reasonably satisfied by HNC data for w R /a = 1 (w A /a ≃ 0.4 for U R /U A ≃ 0.9) and is better approximated for w R /a = 1.5 (w A /a ≃ 0.6 for U R /U A ≃ 0.9). We can deduce from these considerations that to get a phase diagram with two accessible critical points in the fluid region, the three parameters of the potential should be related by the approximate relation in Eq. (7) for w R /a ≫ 1, which reduces to w A ≃ w R /3 for U R /U A ≃ 1.

VI. CONCLUSIONS
The purpose of the present investigation is to understand the role that the different components of the interparticle interaction play in the physical mechanism underlying the liquid-liquid phase transition in one-component systems. Thus we investigate the phase diagram associated with an isotropic pair potential with an attractive well and a repulsive shoulder, described by three parameters, by analyzing for which combinations of the parameters do we find a phase diagram with two critical points in the fluid phase. In a first paper 15 we used MD simulations finding limited ranges of the parameters so that the liquidliquid phase transition was accessible. We also presented a general description based on the MVDWE approach, which rationalized our MD results.
Since MD simulations are precise but time-demanding, we completed this analysis by adopting a different approach with well-known limitations, but one that was extremely fast in terms of computational time, consisting of integral equations in the HNC approximation.
It is important to stress that the drawbacks of the HNC equation are not critical for our purposes. Indeed, we found that the theory, though at best only in qualitative agreement with MD simulations, correctly reproduces the trend according to which the simulated critical points move in the ρ, T plane as the potential parameters are changed 15 . On this basis, we use the theoretical results to estimate the phase behavior of our system over a portion of the parameter space much wider than that explored by numerical simulations. For w R ≪ a it is difficult to extract a clear relation between the potential's parameters.
However, we note that Eq. (4) shows a leading linear relation between U R /U A and w A /w R for w R ≪ a and w A ≪ a, suggesting that the liquid-liquid phase transition could also be found in systems with short repulsive range, if the attractive range is short as well.
For w R ≫ a the situation is more clear. The MVDWE predictions for A = 0 compare well with the HNC results (Fig. 11), leading to the Eq. (6). This equation gives us the intuitive understanding that the repulsive and attractive components of the interaction potential compensate when the attractive volume, weighted by the attractive energy, is equal to the repulsive volume, weighted by the repulsive energy. Moreover, for large w R , Eq. (6) reduces to the simple Eq. (7), whose leading order is 3w A /w R .
In conclusion, we have studied a large number of parameter combinations of a general squared isotropic pair potential with an attraction and soft-core repulsion by using the HNC equation. We verify the general trend previously observed in MD simulations and extend the analysis to a larger number of parameter combinations, with more emphasis on cases with wide soft-core repulsion, particularly difficult to simulate by MD. By using the MVDWE approach introduced in Ref. 15 to interpret the HNC results, we find that the condition A = 0 is well verified for potentials with a large repulsive range and two maxima in the IL. We expect that this condition for a liquid-liquid phase coexistence could be generalized to a continuous isotropic attractive potential U(r) with a wide soft-core repulsion such as where a is the hard-core distance. In particular, this condition brings to Eq. We use U R /U A , w A = c − b, and w R = b − a as independent parameters. sets with w R /a = 0.8 (squares) and with w R /a = 1.0 (diamonds) for 0 ≤ U R /U A ≤ 1.2 and 0 ≤ w A /a ≤ 0.6; sets with w R /a = 1.5 (triangles) for 0 ≤ U R /U A ≤ 6 and 0 ≤ w A /a ≤ 3.
Parameters outside these regions have not been investigated.  Fig. 9 for w R /a = 0.6 (set 1 denoted by circles), w R /a = 0.8 (set 2 denoted by squares), w R /a = 1.0 (set 3 denoted by diamonds), w R /a = 1.5 (set 4 denoted by triangles). Error bars represent the interval in Fig. 9.
Lines are one-parameter fits with Eq. (4) of the sets: for set 2 (squares) the fitting parameter is Since we only have three points for set 1, to avoid a fit with a large indeterminacy on the parameters we arbitrarily chose A/(U A V SC ) = 1 to show that the data are consistent with Eq. (4).