Metastable liquid-liquid phase transition in a single-component system with only one crystal phase and no density anomaly

We investigate the phase behavior of a single-component system in 3 dimensions with spherically-symmetric, pairwise-additive, soft-core interactions with an attractive well at a long distance, a repulsive soft-core shoulder at an intermediate distance, and a hard-core repulsion at a short distance, similar to potentials used to describe liquid systems such as colloids, protein solutions, or liquid metals. We showed [Nature {\bf 409}, 692 (2001)] that, even with no evidences of the density anomaly, the phase diagram has two first-order fluid-fluid phase transitions, one ending in a gas--low-density liquid (LDL) critical point, and the other in a gas--high-density liquid (HDL) critical point, with a LDL-HDL phase transition at low temperatures. Here we use integral equation calculations to explore the 3-parameter space of the soft-core potential and we perform molecular dynamics simulations in the interesting region of parameters. For the equilibrium phase diagram we analyze the structure of the crystal phase and find that, within the considered range of densities, the structure is independent of the density. Then, we analyze in detail the fluid metastable phases and, by explicit thermodynamic calculation in the supercooled phase, we show the absence of the density anomaly. We suggest that this absence is related to the presence of only one stable crystal structure.


I. INTRODUCTION
Soft-core potentials have been widely used to study a variety of systems such as liquid metals, metallic mixtures, electrolytes and colloids, as well as anomalous liquids, like water and silica [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]. In these models, the specific structural characteristic at the molecular (or atomic) level are neglected and the molecules (or atoms) are represented by simple spheres. Quantum effects (such as the quantum nature of chemical interactions) and classical effects (such as the Coulomb interaction) are modeled through a phenomenological isotropic pair potential. The advantages of this approach are that while these potentials are simple enough to be treated analytically [16], they still allow a qualitative comparison with the experiments. Moreover, they can be studied by means of numerical simulations less time-consuming than those of realistic models [17].
We consider an off-lattice model in three dimensions (3D) [11] related to the soft-core potentials studied by Hemmer and Stell [2] for solid-solid critical points. Our model shows a phase diagram with two fluid-fluid phase transitions, a feature recently seen in experiments on phosphorus [18] and confirmed by specific simulations [19]. In Ref. [11] we showed that both first-order fluidfluid phase transitions end in critical points, a low-density * Present address: SMC, Dipartimento di Fisica, Università "La Sapienza", P.le A. Moro 2, I-00185 Roma, Italy; Electronic address: franzese@na.infn.it critical point C 1 , and a high-density critical point C 2 . For the considered potential, both transitions occur in the supercooled phase with respect to the crystal phase.
The hypothesis of a second critical point has been proposed [20] as a way to rationalize the density anomalyi.e., the expansion upon isobaric cooling-in networkforming fluids, such as water [1,20,21], carbon [22], silica [23] and silicon [24]. Consequently, both the experimental [25,26,27,28,29] and the theoretical [7,8,9,10,13,14,15,30,31] investigations about the possibility of a second critical point have been focused on systems with the density anomaly (Fig. 1a). However, the results of Ref. [11] have shown that the presence of the critical point C 2 does not necessarily induce the density anomaly, indicating that the quest for simple liquids with two critical points is not restricted to systems with densities exhibiting anomalous behavior (Fig. 1b). In this paper we push forward the analysis, by studying the equilibrium phase diagram and showing that the system introduced in Ref. [11] has one single crystal phase, in the range of considered densities, suggesting that the absence of density anomalies is related to the presence of only one stable crystal structure.
To reach this goal we organize the work in the following way. After the definition of our soft-core potentials (Section II), (i) we describe in detail the integral equations in the hyper-netted chain (HNC) approximation (Section III.A) and (ii) we generate new calculations using different assigned parameters for the pair potential, thus establishing a bridge between the potential studied in Ref. [11] and the potential investigated in Ref. [10,15] and rationalizing how-in addition to the critical point Phase diagram with critical point C between the gas and the uniform liquid phase at high T and low P and critical point C ′ between the low-density liquid (LDL) and the high-density liquid (HDL) at low T and high P . This phase diagram has been proposed for water and has been shown to be consistent with the density anomaly. (b) Phase diagram with critical point C1 between the gas and the LDL at low T and low P and the critical point C2 between the gas and the HDL at high T and high P . Increasing P at constant temperature Ta below C1 (dashed line at Ta), the system undergoes a first order phase transition between the gas and the LDL phase, followed by a first order phase transition between the LDL and the HDL phases. Increasing the pressure at constant temperature T b above C1 but below C2 (dashed line at T b ), the system undergoes only a first order phase transition between the gas and the HDL. The square represents the gas-LDL-HDL triple point. This phase diagram has been found in Ref. [11] with no evidence of the density anomaly.
C 1 -the critical point C 2 arises as a function of the parameters of the pair potential (Section III.B). Then, (iii) we describe in detail the molecular dynamics (MD) simulations, studying the equilibrium phase diagram and finding that the only stable crystal structure, in the range of simulated densities, has two characteristic lattice distances (Section IV.A). (iv) By MD simulations we analyze also the supercooled liquid phase and the metastable fluid-fluid phase transition (Section IV.B); (v) in order to construct an accurate phase diagram and to study the finite size effect, we perform new MD simulations in addition to the calculations presented in Ref. [11] (Section IV.C). Hence, (vi) we study the radial distribution function, comparing the HNC predictions with the MD results and analyzing the composition of the system within the fluid-fluid coexisting regions (Section V). Finally, (vii) we address the density anomaly issue by presenting the explicit thermodynamic calculation, based on the MD phase diagram, that allows us to exclude the presence of the density anomaly (Section VI). We give our conclusions in Section VII.

II. SOFT-CORE POTENTIALS
Among the isotropic potentials, much attention has been devoted to soft-core potentials, which have a finite (soft-core) repulsion added to the infinite (hard-core) repulsion. The infinite repulsion is due to the impenetrability of the spheres. The finite repulsion represents the combination of all the quantum and classical repulsive effects averaged over the angular part. It has been shown [4] that a weak effective repulsion can be derived by a first principle calculation for liquid metals [3]. To understand the possibility of the solid-solid critical point in such material as Ce and Cs, Hemmer and Stell [2] proposed a soft-core potential with an attractive interaction at large distances, performing an exact analysis in 1D. Over the past thirty years, several other soft-core potentials with one or more attractive wells were proposed and studied with approximate methods, or numerical simulations in 2D, to rationalize the properties of liquid metals, alloys, electrolytes, colloids and the water anomalies [1,2,3,4,5,6,7,8,10,11,12,13,14,15]. It has been recently shown [11], for the first time in 3D with MD, that an appropriate soft-core potential with an attractive well is able to show two supercooled liquids of different densities, with two critical points. Similar results have been reported for a soft-core potentials with a linear repulsive ramp by Monte Carlo simulations [14]. Following Ref. [11], we define the isotropic pair potential U (r), as a function of the pair distance r (inset in Fig.2): where a is the hard-core distance and b is the soft-core distance. For a ≤ r < b, the spheres interact with a finite (soft-core) repulsive energy U R . For b ≤ r < c, the pair's interaction is attractive with energy −U A < 0. The distance c is the cut-off radius beyond which the pair's interaction is considered negligible. For sake of comparison with Ref. [10,15], we will consider both U R > 0 and U R < 0. The step-like shape in Eq. (1) has the advantage of being defined by only three parameters: the width of the soft-core in units of the hard-core distance w R /a ≡ (b − a)/a, the width of the attractive well in the same units w A /a ≡ (c−b)/a and the soft-core energy in units of the attractive well energy U R /U A . To explore the phase diagram of the model as a function of these three parameters in an approximate, yet fast way, we use the integral equations in the HNC approximation.

III. THE INTEGRAL EQUATIONS IN THE HYPER-NETTED CHAIN APPROXIMATION
In this section we present details of the integral equations and the HNC approximation adopted in Ref. [11] and in this work. The radial distribution function g(r) plays a central role in the physics of fluids [16]. This quantity is proportional to the probability of finding a particle at a distance r from a reference one and is the ratio of local to bulk density at distance r, with where n(r) is the number of particles at a distance between r and r+dr from the reference particle and ρ is the number density, assumed to be independent of r (uniform system). The radial distribution function goes to 1 for large r and is always 1 for a random spatial distribution of particles. To represent deviations from randomness, the total pair correlation function h(r) ≡ g(r) − 1 is introduced. These functions are relevant because they are directly measurable by radiation scattering experiments and are related to the thermodynamic properties of the fluid. A fundamental relation between structure and thermodynamics is given by where K T ≡ (∂ρ/∂P ) T /ρ is the isothermal compressibility, P is the pressure, and k B is the Boltzmann constant.
Provided that the particles interact through pairwiseadditive forces, other thermodynamic properties of the fluid-such as the internal energy-can be expressed using integrals over the pair correlation function. The function h(r) is the result of the interaction of all the particles in the system. Formally, h(r) can be decomposed into (i) the contribution coming from the direct interaction between two particles at distance r, called c(r), and (ii) the contribution due to the indirect interaction propagated through any other particle in the system. This second contribution is written in turn as an integral convolution of direct correlations and total pair correlation.
This decomposition, for uniform systems, is expressed by the Ornstein-Zernike relation Equation (4) is also the formal definition of c(r). Both h(r) and c(r) in Eq. (4) are unknown functions, thus to solve this equation, one needs another relation (closure) between these two functions. This relation is provided by the diagrammatic expansion of g(r) [16] which, after formal summation, yields the functional relation where U (r) is the inter-particle potential, β ≡ 1/(k B T ) and d(r) is the sum over a specific class of diagrams (bridge diagrams) [16]. Since d(r) cannot be calculated exactly, one resorts to approximate expressions. The simplest approximation assumes d(r) = 0 (HNC closure) [32]. One expects this approximation to work better at lower ρ, where the direct correlation function c(r) is more relevant then the correlation propagated through the other particles. However, our results (see Section V) will show that this intuitive observation is not straightforward, at least for soft-core potentials.
A. The iterative procedure The solution of the integral Eqs. (4,5) with the HNC closure is obtained through a numerical iterative procedure whose essential scheme is the following. Under the assumption d(r) = 0, one can write g(r) as where the function θ(r) ≡ h(r) − c(r) has the remarkable property of being a continuous function of r, even for discontinuous potentials (as in this paper). From the definitions of h(r), θ(r) and Eq. (6), one can derive the equation By using the Fourier transformf ( q) ≡ f ( r) exp(i q· r) d r defined for a generic function f ( r), from Eq. (4) we obtain h(q) =ĉ(q) + ρĉ(q)ĥ(q) .
Or, using the definition of θ(r), we havê The numerical iteration is based on Eqs. (7) and (9). We start by choosing an initial guess for θ(r). A reasonable input, at least at high temperatures, is the θ(r) of a fluid of hard spheres with diameter a. In fact, at high temperatures our potential can be approximated with a simple hard-core repulsion. We can calculate the corresponding θ(r) by making use of the Percus-Yevick integral equation [16] which for hard spheres can be solved analytically. Next, at constant ρ, we decrease the temperature of δT and we perform calculations at fixed ρ and T by using as input the θ(r) obtained as solution at ρ and T + δT .
From the chosen guess of θ(r) we calculate c(r) by using Eq. (7). Its Fourier transformĉ(q) is used in Eq. (9) to calculateθ(q). Its inverse Fourier transform provides a new θ(r) that is used as a new input for the next cycle. We evaluate the functions on M = 2048 discrete points r m = mδr, with m = 1, . . . , M and δr = 0.01a. Successive iterations of the elementary cycle defines a succession θ (k) (r), where k = 1, 2, . . . is the number of the iteration. If the difference between two consecutive elements of this succession, decreases for increasing k, the succession converges towards a θ * (r) that is solution of our integral equations. The iteration process is stopped when ∆ ≤ 10 −7 .
Based on this iterative procedure, different algorithms can be used to improve the accuracy and rapidity of convergence of the numerical solution of HNC equations. However, independent of the algorithm used, there exists 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, defining an instability line in the ρ-T plane.

B. The HNC instability line
The nature of the locus of instabilities of the HNC equation and its relationship with the spinodal line of the fluid was investigated for a hard-core potential plus an attractive Yukawa tail in a number of papers [33,34]. These studies showed that the isothermal compressibility does not diverge as the temperature is lowered and the instability region is approached from above. This conclusion was definitively assessed through extensive numerical calculations [35] both for the hard-core Yukawa fluid and other model potentials, showing that this behavior is directly correlated to the existence of multiple HNC solutions. The analysis developed was based on a careful treatment of the low-k behavior of the Fourier transforms of the correlation functions required by the iterative procedure. A further theoretical support to these results was given by an analysis [36] on models for an ionic fluid and a monoatomic Lennard-Jones fluid.
In the light of the above mentioned studies, an identification of the instability line of the HNC equation with the spinodal line of the fluid, which is characterized by a diverging compressibility, is not possible. Keeping in mind this limitation, one can nevertheless observe that for a large number of simple fluid pair potentials the shape of the instability line qualitatively resembles the region of spinodal decomposition of the fluid. Also for our potential the comparison of the HNC calculations with the MD results (Section V) shows that the HNC instability line is qualitatively consistent with the spinodal line [37]. Thus, studying the modifications of the instability line as the potential parameters are varied can yield some approximate, yet useful, informations on the phase behavior of the fluid. For UR/UA = −0.5, the pair potential recovers the one studied in 1D and 2D in Ref. [10,15]. The symbols represent the calculations and the lines are guides for the eyes. Inset: the isotropic pairwise-additive potential U (r) as a function of the inter-particle distance r; a is the hard-core distance, b is the soft-core distance, c is the maximum attractive distance, −UA < 0 is the attractive energy and UR is the repulsive energy.

C. The results
First, we calculate the instability line of the HNC equations for the potential investigated in Refs. [10,15]. The corresponding parameters are w R /a = 0.4, w A /a = 0.3 and U R /U A = −0.5. In this case, the soft-core is given by two attractive wells with different depths. Calculations in 1D and 2D [10,15] have shown a water-like density anomaly. Therefore, it is interesting to analyze the phase diagram in 3D. However, the instability line for this case (Fig. 2) is similar to the spinodal line usually exhibited by a simple fluid, e.g. interacting via a Lennard-Jones potential with the maximum of the spinodal line corresponding to the liquid-gas critical point. Upon increasing U R /U A to 0.5 ( Fig. 2), the only evident change of the instability line is a shift toward a lower T as a result of the overall decrease of the inter-particle attraction, with no hints of a second critical point. A small shift to lower ρ is also seen. This behavior is more evident for larger w R .
Next, we consider a potential with a larger w R (w R /a = 0.7, w A /a = 0.3). The instability line is calculated for several values of U R /U A (Fig. 3). Upon increasing U R /U A , we now find not only the shift to a lower T , but also an evident shift of the maximum of the line (i.e. the critical point, assuming that the instability line represents the behavior of the spinodal line) to a lower ρ. This result can be rationalized by observing that, passing from U R < 0 to U R > 0, the soft-core becomes more and more difficult to penetrate and the system passes from a potential with a hard-core a and an effective attractive range w A + w R to a potential with an effective hard-core b, for U R /U A large enough, and an attractive range w A . As a consequence of the increase of the effective hardcore, the critical density decreases and, as a consequence of the decrease of the effective overall attraction, the critical temperature decreases. Comparing Fig. 2 and Fig. 3, we notice an important difference. In the case of larger w R (Fig. 3), as U R /U A increases, the temperature of the instability line does not decrease for increasing ρ, but becomes rather flat. This result suggests that the instability line might develop a second maximum at larger values of ρ for even larger w R .
We thus consider a potential with w R /a = 1.0 and w A /a = 0.2. The results (Fig. 4) show that for 0.4 ≤ U R /U A ≤ 0.6, the instability line has two well-distinct local maxima, suggesting the possibility of two critical points in the phase diagram for the fluid phases [38]. For 7, the instability line shows just one maximum, similar to the typical spinodal line of a fluid of hard spheres with diameter a or b, respectively, attracting via a square well of width w A . As a consequence of this analysis, we choose w R /a = 1, w A /a = 0.2 and U R /U A = 0.5 as the set of parameters for the potential used in the MD calculations in 3D [11].

IV. THE MOLECULAR DYNAMICS APPROACH
In this Section we give extensive details on the MD method and we extend the analysis performed in Ref. [11], including new calculations for the crystal phase, the crystal nucleation process and the metastable phases. We perform MD simulations at constant number of particles N of unit mass m, at constant volume V , with periodic boundary conditions and at constant average temperature T . We present the results for N = 490, 720 and N = 1728. The average temperature is set by coupling the system to a thermal bath at the assigned T , with a thermal exchange coefficient per particle between the system and the bath equal to k = 0.015 (U A /m) 1/2 k B /a. We use a standard collision event list algorithm [39] to evolve the system and a modified Berendsen method to achieve the desired T [40].
The pressure is calculated by using the virial expression for a step potential [17] with ′ N i,j sum over the particles pairs (i, j) undergoing a collision in the time interval ∆t ≡ (10 5 ma 2 /U A ) 1/2 , hereafter used as unit of time, and with ∆v i are the velocities of the particle i at position r i before and after the collision with particle j at position r j .

A. The crystal
First, to locate the equilibrium crystal line we simulate a crystal seed surrounded by the gas. We prepare a crystal seed by cooling at T = 0.45U A /k B a gas configuration with density ρ = 0.018/a 3 .
The crystal (Fig. 5) is the effect of the competition between the hard-core repulsion at distance r = a and the attraction at distance r = b. The resulting structure is reminiscent of the close packing of hard spheres with diameter a or b, but the competition gives rise to new symmetries (Fig. 5). The minimum in the interparticle interaction potential at b ≤ r < c would induce a face-centered-cubic crystal with lattice space ranging from b to c and a characteristic six-fold symmetry on one projection plane with seven particles to form a triangular lattice. However, due to the soft-core, the system can allocate particles at the hard-core distance a = b/2. This induces a twelve-fold symmetry, placing an average number of twelve, almost on-plane, nearest neighbor (nn) particles at a distance 1 ≤ r/a 1.2 to form a dodecagon around two particles. These two nn particles are next nearest neighbors to the dodecagon, at a distance 2 ≤ r/a ≤ 2.2, and are placed on a line almost perpendicular to the plane individuated by the dodecagon (Fig. 5b). This structure is distorted in such a way to form non-closed chains of nn particles, that wrap along another axis to give rise to an eight-fold symmetry (Fig. 5c,d).
By analyzing the crystal structure obtained from the MD simulations, we conclude that the position of the particles in the crystal can be described by r = i · a +  Table II), whose long sides form two triangles with two particles on the same plane ( rm with m = 5, 6 in Table II), and the short sides forms two tetrahedra, each with two more particles ( rm with m = 7, . . . 10 in Table II).
(b) The crystal configuration rotated by π/4 around a central horizontal axis shows the 8-fold symmetry seen in Fig. 5a and d. (c) A further rotation of π/4 around the same axis shows the dodecagons seen in Fig. 5b and c. (d) A rotation of π/2 around a central vertical axis shows again the octagons seen in Fig. 5a and d.  j · b + k · c + r m , where r m , for m = 1, . . . 10, are the coordinates, with respect to the center of the cell, of the 10 particles forming a crystal cell, a, b and c are the lattice vectors describing the position of the center of the cell and i, j and k are integers such that i + j + k is even. We estimate the lattice vectors (Table I) and the coordinates of the particles forming the cell (Table II) after an equilibration time 10 2 ∆t at T = 0.03U A /k B for of an artificial crystal placed in the vacuum (Fig. 6). The resulting density of the crystal is a 3 ρ ≃ 0.39. Surface effects could be responsible for the tilt that can be seen in Fig. 6d. In a system with N = 720, this tilt II: The coordinates rm = (xm, ym, zm) for m = 1, . . . 10 of the 10 particles forming a crystal cell, with respect to the center of the cell, as obtained after an equilibration time 10 2 ∆t at T = 0.03UA/kB of the artificial crystal (Fig. 6). The characteristic distances, with an error on the last decimal digit, are: l1/a = 0.53, l2/a = 0.59, l3/a = 0.07, l4/a = 1.43, l5/a = 0.05, l6/a = 0.03, l7/a = 1.40, l8/a = 0.52. The errors decrease as the square root of the time of averaging. For each particle m of the cell, we denote with n sc m the number of particles in the crystal at a distance r ≤ b (in the soft-core) and with n aw m the number of particles at a distance b < r ≤ c (in the attractive well). disappears when the sample is equilibrated at higher T (Fig. 7). We compare the g(r) (Fig. 8) of the MD crystal in Fig. 5 and of the artificial crystal in Fig. 6, both equilibrated at T = 0.48U A /k B . Both functions show peaks located at the same distances, with two large peaks at r/a = 1 and r/a = 2, consistent with the presence of the two characteristic distances a and b in the potential. The comparison confirms that the proposed crystal is a good representation of the crystal structure generated by the MD simulation. The slightly different intensities of the peaks of the g(r) of the two systems are probably due to the defects of the MD crystal.
The validity of the artificial crystal as a good description of the real crystal structure is confirmed also by the evolution of the potential energy per particle (inset in Fig. 8) when the MD crystal and the artificial crystal are heated, from the configuration equilibrated at T = 0.48U A /k B , to T = 0.60U A /k B . Both samples equilibrate to the same energy. The starting potential energy is, as expected, in both cases greater than the ground state energy U 0 /N = −8.45U A , calculated from the number of particles at distance a ≤ r < b and at distance b ≤ r < c (Table II), due to surface effects. We find analogous results for the evolution of the kinetic energy.
To test if the system has more than one crystal structure as function of the density, we cool at T = 0.6U A /k B a fluid configuration equilibrated at T = 0.8U A /k B and ρ = 0.267/a 3 , and compare the resulting g(r) with the case of the crystal seeds at T = 0.45U A /k B and ρ = 0.018/a 3 , finding no relevant differences ( Fig. 9). At the same time, the attempt of finding alternative artificial crystal structures has revealed, after an appropriate equilibration, that the only stable structure is the one presented in Fig. 6. We therefore assume that the system, at least for this choice of the potential's parameters, has one single crystal structure independent of ρ, within the considered range of densities.
Starting from a configuration with the crystal seed described above, we equilibrate the system at different densities and temperatures. We define the system to be in the solid phase if, after a time 10 6 (ma 2 /U A ) 1/2 ≃ 3 × 10 3 ∆t, the crystal seed is growing, or we consider it in a fluid phase if the seed is melting. The cases in which the trend is not clear within the simulation time are considered as belonging to the first-order transition region [41]. The crystallization pressure rapidly increases with ρ and T , giving a first-order transition line (in the thermodynamic limit) that separates the equilibrium P - ρ phase diagram in a high-T fluid and a low-T crystal (Fig. 10) [11].

B. The supercooled liquids
At equilibrium, there is no (stable) liquid phase. A phase diagram without equilibrium liquid phase is expected [42,43] for systems with an inter-particles potential with a narrow attractive part, such as the one we are considering here. However, the liquid is present as a metastable (supercooled) phase with respect to the crystal phase [11]. To study the metastable phase diagram, we equilibrate the system for each ρ from a gas configuration at T = 0.70U A /k B and then rapidly cool it to the desired T ≥ 0.57U A /k B , calculating P , g(r) and the total potential energy U ≡ N i<j U (|r i − r j |). We find that the supercooled fluid phase has a lifetime longer than 3 × 10 3 ∆t (the standard length of our simulations) for ρ 0.20/a 3 at T ≈ 0.57U A /k B , for ρ 0.27/a 3 at T ≈ 0.65U A /k B and for ρ 0.34/a 3 at T ≈ 0.70U A /k B . The system is equilibrated in the fluid phase for t ≥ 20∆t, after which we average P , g(r) and U over the time. We calculate each state point by averaging the configurations for 3 × 10 2 ∆t ≤ t ≤ 3 × 10 3 ∆t. We estimate the errors by dividing the configurations in 90 non-overlapping intervals of 30∆t, that we assume independent.
For larger ρ, the system spontaneously crystallizes (homogeneous nucleation process). Thus, we only average over configurations that occur before nucleation. To be certain that our estimates are carried out in the fluid phase, we study where r j (t) is the position of particle j at time t and q is the wavevector. At equilibrium, the average of S( q, t), over the time and the wavevectors with the same module, is the structure factor S(q), describing the spatial correlation in the system. Therefore, S( q, t) describes the timeevolution of the spatial correlation along the wavevector q. In particular, for a crystal-like configuration, with a long-range order, there is at least one wavevector such that S( q, t) ∼ O(N ), while for a fluid-like configuration is S( q, t) ∼ O(1) for all q.
The time-evolution of S( q, t) for a typical simulation inside the nucleation region is presented in Fig. 11.  (1) for any q; in the time interval "B" with 700∆t ≤ t ≤ 1800∆t, for six wavevectors there is an increase in S(q = 1/a, t); in the interval "C" with 2100∆t ≤ t ≤ 3000∆t, for the same six wavevectors there is a larger increase. (b) As in (a) but for q ≈ 12/a ≈ 4π/a; in this case there is a large increase in S(q ≈ 12/a, t), more then one order of magnitude, only in the time interval "C" for several wavevectors, revealing the formation of a crystal seed. (c) The structure factor S(q), given by the average over the dimensionless wavevectors aq with the same modulus and the average over the time intervals "A", "B" and "C" of S( q, t). The curves for "B" and "C" are offset by 1.5 and 3, respectively. All the curves go to 1 for large q. In the interval "A", S(q) is liquid-like. In the interval "B", S(q) is still liquid-like but with an increase for q → 0, while in the interval "C" is solid-like, with two large peaks at q ≈ 2π/(b/2) = 2π/a and q ≈ 2π/(a/2), corresponding to the soft-core radius and the hard-core radius, respectively, and a large value for q → 0.
To limit the computational effort, we consider 9 × 10 4 wavevectors with modulus q ≤ 100/a, that is much larger than the wavevectors of the largest peak of the crystal structure factor, q ≈ 2π/(a/2), corresponding to the hard-core radius (Fig. 11c). Three different regimes can be distinguished in the example in Fig. 11.
(i) A short-time regime "A", in which S( q, t) ∼ O(1) for any q (q = 1/a and q ≈ 12/a ≈ 4π/a are shown in Fig. 11a,b). Averaged on this interval (curve A in Fig. 11c) the S(q) is fluid-like.
(ii) An intermediate-time regime "B", in which S( q, t) for q = 1/a has an increase, but has no increase for q ≈ 4π/a. Averaged on this interval (curve B in Fig. 11c) the S(q) is fluid-like, but with an increase for q → 0. This increase indicates an increase of K T , according to the equation The number of pairs of particles at a relative distance ri < r ≤ ri+1 for the MD configuration in the inset, with r0 = 0, r1/a = 1 and ri+1 − ri = a/10 for i > 1. The histogram shows a large maximum corresponding to the attractive range 2 ≤ r/a < 2.2, a broad maximum around r/a = 1.2 and a small number of pairs at the hard-core distance r = a. Therefore, the preferred relative distance for pairs of particles within the soft-core range 1 ≤ r/a < 2 is, for this configuration, r/a ≃ 1.2.
where we use the Eq. (3) and the definition S( q) ≡ 1 + ρĥ( q). The increase of K T is associated with the phase separation into two fluids with different densities.
To help visualize the phase separation, in Fig. 12 we show the three planar projections of the 3D configuration corresponding to the largest peak in the time interval "B" for q = 1/a (Fig. 11a). By dividing the box into two equal parts, the histograms of the number of particles in each part (Fig. 12) show a separation in density approximately at half boxlength, corresponding to q = 4π(ρ/N ) 1/3 ≈ 1/a for ρ = 0.27/a 3 and N = 490, in agreement with the peak at q = 1/a for the curve B in Fig. 11c. In each projection, it is possible to see regions of high density and low density (Fig. 12).
To better quantify the phase separation occurring in the configuration in Fig. 12, we present (main panel in Fig. 12)  pairs of particles at a relative distance r i < r ≤ r i+1 , where r i+1 − r i = a/10. The histogram has a broad maximum around the distance r/a = 1.2 in the soft-core range, showing that there exists a subset of pairs of particles that are at a preferred distance 1.1 < r/a ≤ 1.2. This subset is the HDL that has a non-uniform distribution over space (Fig. 13b), consistent with the phase separation.
(iii) A long-time regime "C", in which is S( q, t) ∼ O(N ) for q ≈ 4π/a, revealing the crystal nucleation process. The S(q) averaged over this time interval (curve C in Fig. 11c) is solid-like. In the same interval, S( q, t) for q = 1/a has a large increase, corresponding to the increase of K T (S(q) increases for q → 0), which is consistent with the phase separation between the fluid and the crystal. As an example, in Fig. 14 we show the last configuration of the time series in Fig. 11, where the crystal structure, already observed in Fig. 5, is clearly seen.
The example in Fig. 11b shows the formation of a high density fluid phase within the time interval "B", followed by the nucleation of the crystal phase. The onset of the nucleation is marked by a large increase of S( q, t) for all FIG. 14: The last MD configuration in the time series in Fig. 11. A crystal nucleus surrounded by gas is clearly seen. Bonds connect particles at distance r/a ≤ 1.1. The 3D perspective is given as in Fig. 5. The radius of particles is not in scale with the distances.
the wavevectors corresponding to the peaks in the crystal S(q) and by a large step-like decrease of energy.
C. The phase diagram and the finite size effect By repeating the analysis described above for all the simulations inside the region with nucleation-and discarding the data corresponding to the formation of the nucleus-it is possible to calculate the state points corresponding to the metastable fluid phase. The phase diagram in Fig. 10 is based on averages over a total of 10 5 -10 6 configurations in the fluid phase, accumulated in independent runs.
For completeness we recall here the main features of the phase diagram in Fig. 10 and presented in Ref. [11]. The (mechanically unstable) region at high ρ for T 0.67U A /k B , where P decreases for increasing ρ, denotes the coexistence of gas and HDL. The unstable region at low ρ for T 0.603U A /k B (inset in Fig. 10) denotes the coexistence of gas and LDL. The coexistence lines are obtained by using the Maxwell construction of the equal areas [1], suggesting the presence of a gas-LDL-HDL triple point.
By definition, the spinodal lines (limit of stability of each phase with respect to the coexisting phase) meet the coexistence lines in a critical point. Therefore, by interpolation we estimate the gas-LDL critical point C 1 at k B T 1 /U A = 0.603±0.003, a 3 ρ 1 = 0.10±0.01, a 3 P 1 /U A = 0.0171 ± 0.0005 and, the gas-HDL critical point C 2 at k B T 2 /U A = 0.665 ± 0.005, a 3 ρ 2 = 0.306 ± 0.020, a 3 P 2 /U A = 0.10 ± 0.01. These values are consistent with the linear interpolations of the MD isotherms (Fig. 10). The phase diagram resulting from the MD calculations is, as expected, in agreement with the time-dependent analysis of the structure factor presented above. For example, the case presented in Fig. 11 corresponds to a state point inside the gas-HDL coexistence region at a density higher then the crystal nucleation density for T = 0.62U A /k B . The nucleation of the (metastable) HDL phase is thus followed by the crystal nucleation.
To estimate the finite size effect in our calculations, we compare the results for N = 490 and N = 1728 for an isotherm below both critical points (Fig. 15). The calculations do not show any relevant finite size effect, suggesting that the MD results for N = 490 are reliable.

V. THE RADIAL DISTRIBUTION FUNCTION ANALYSIS
The interpretation of the HNC instability line is qualitatively consistent with the MD spinodal line for the corresponding set of the potential's parameters. The projection of the MD spinodal line in the T -ρ plane (not shown) has the same characteristics of the HNC instability line, with two local maxima and one local minimum. In both approaches, the high-ρ local maximum occurs at a temperature higher then the temperature of the low-ρ maximum and the presence of a triple point is suggested by the presence of the local minimum. The quantitative HNC predictions for the locations of the two critical points are, as expected, only partially consistent with the MD results. It is remarkable that the HNC estimates of the density of the low-ρ local maximum (ρ ≈ 1/a 3 ) and the temperature of the high-ρ local maximum (T ≈ 0.65U A /k B ) of the instability line are close to the corresponding MD results for C 1 and C 2 , respectively.
An estimate of the agreement between the two methods can be evaluated by comparing the calculations for g(r) within the two approaches (Fig. 16). In contrast with what could be suggested by the nature of the HNC approximation-i.e. the under-estimate of the indirect correlation-the agreement is better at intermediate ρ than at low ρ (Fig. 16). In particular, at low ρ the HNC approximation underestimates the probability of a particle penetrating the soft-core or entering the attractive well. At higher ρ, instead, the estimates of the g(r) within the two approaches are almost indistinguishable.
The g(r) of the low-ρ fluid is characterized by a large peak at r = b corresponding to the shortest attractive distance. As a consequence of the increase of the density, the peak at the hard-core distance r = a increases while the peak at r = b decreases, and additional peaks at r/a = 3, 4, ... appear. In Fig. 17 we present the calculation of the g(r) for the gas phase, the gas-HDL coexisting region and the HDL phase. In particular, by combining the radial distribution functions evaluated in each pure phase, we can estimate the composition of the mixed phase. For example, at T = 0.64U A /k B the radial distribution function calculated at ρ 0 = 0.302/a 3 is g 0 (r) ≃ X 1 g 1 (r) + X 2 g 2 (r), where g 1 (r) and g 2 (r) are the radial distribution functions at the same T and at ρ 1 = 0.223/a 3 and ρ 2 = 0.349/a 3 , respectively, with X 1 = 0.3 and X 2 = 0.7 (Fig. 18), revealing that the system is composed approximately by 30% of gas and 70% of HDL. From the Eq. (2), by using the g(r) calculated from the MD simulations, we evaluate the average number of particles N (r) = dN (r) within a sphere of radius r (Fig. 19). This analysis reveals that the number of particles ∆N within the repulsive range and within the attractive range increases linearly with ρ (inset Fig. 19) and that the increase is faster within the attractive range ( Fig. 19), for the densities we studied. In particular, the number of particles within the attractive range b ≤ r < c increases from 2.5 to 9 ≪ 10 m=1 n aw m /10 = 19.2 estimated for the artificial crystal (Table II).

VI. ABSENCE OF A DENSITY ANOMALY
In Ref. [11] it has been noted that the possibility of a second fluid-fluid critical point is not necessarily restricted to systems with a density anomaly, at least from a theoretical point of view. Here we present the explicit thermodynamic calculations for this result.
The defining relation for the density anomaly is given by or for the Maxwell relation where S is the entropy. Since holds for a mechanically stable phase, Eq. (15) can be rewritten as From the differential expression of the thermodynamic potential at constant T , we know that where E ≡ U + K is the total energy, with U and K total potential and kinetic energy, respectively. Therefore, it is at constant T and we can rewrite the density anomaly condition in Eq. (18) as at constant T .
To calculate the left-hand side of Eq. (21), we need to evaluate (∂U/∂V ) T . In Fig. 20, we show our MD   Table III calculation for U (ρ) at constant T . All the MD points can be fitted with a third degree polynomial in ρ. The fitting parameters are given in Table III and are used to calculate the derivative (∂U/∂V ) T , shown in Fig. 21.
Our calculations show a potential energy U increasing with V (inset Fig. 21), with a derivative always positive, thus wherever P is positive, the condition in Eq. (21) is not satisfied and there is no density anomaly.
In the region where P < 0 (at low T and small V ), the derivative (∂U/∂V ) T rapidly increases in such a way that Eq. (21) is never satisfied. Particularly, in the range of volumes considered, it is always (∂U/∂V ) T − 0.025U A /a 3 > 0, where P = −0.025U A /a 3 is the minimum pressure, reached for T = 0.6U A /k B and V /N = 3.31a 3 (Fig. 10). These results suggest that the density anomaly is ruled out for this choice of parameters. At this stage is not clear if it is ruled out for any choice of parameters for our potential in 3D (see Ref. [15]).

VII. SUMMARY AND CONCLUSIONS
We analyzed the phase diagram of a soft-core potential, similar to potentials used in systems such as protein solutions, colloids, melts and in pure systems such as liquid metals. We use two independent numerical methods, integral equations in the HNC approximation and MD simulations. The comparison of the HNC results with previously-proposed soft-core potentials suggests that the system has two fluid-fluid phase transitions for an appropriate choice of parameters-energy and width-of the repulsive soft-core. We select a set of potential parameters with a narrow attractive well that gives a HNC instability line with two maxima and suggests the presence of two critical points.
The MD analysis shows, in agreement with the previous results for potentials with a short range attraction [42], a phase diagram with no stable liquid phase. We analyze the crystal structure, characterized by the competition between the attractive interaction at distance r = b and the repulsive interaction at r = a < b. We show that the crystal, with 8-fold and 12-fold symmetries, is independent on the density, within the considered range of densities.
Hence, we study the metastable liquid phase, at temperatures above and below the line of spontaneous crystal nucleation. We find two liquids in the supercooled phase, the LDL and the HDL, with two fluid-fluid transitions ending in two critical points, the gas-LDL critical point C 1 and the gas-HDL critical point C 2 , as already shown in Ref. [11]. Here we improve our estimate of the phase diagram and we verify the absence of relevant finite size effects in the MD results.
We compare these results with the HNC calculations, concluding that the HNC approximation underestimates the effect of the attractive interaction and overestimates the effect of the repulsive interaction at low ρ, and is in good agreement with the MD results at intermediate ρ.
Finally, by explicit calculation, we show that the condition for the density anomaly is never satisfied in the range of T and V considered here, as announced in Ref. [11]. Our results suggest that the density anomaly is always ruled out for this choice of potential parameters.
In conclusion, the results of this paper evoke an intriguing relation between the absence of the density anomaly and the presence of a single crystal phase, with higher density than the liquid phases, in systems with two fluid-fluid phase transitions. This relation, that deserves greater investigation, is consistent with the fact that the substances with the density anomaly, and a hypothesized second liquid-liquid critical point, have more than one crystal structure, as in the case of water or carbon or silica.