Numerous studies have attempted to provide evidence of a second critical point in water2,7,8,9, silicon10 and other tetrahedral liquids. Unfortunately, in all cases its expected location is well within the region where these liquids spontaneously crystallize and so no conclusive evidence has ever been reported. To observe critical fluctuations, the liquid must survive in its metastable supercooled state long enough to equilibrate, and the experimental observation time needs to be smaller than the spontaneous nucleation time. In principle, these conditions are easier to meet in simulation studies, in which heterogeneous nucleation is absent by construction and homogeneous nucleation is slowed down by the small size of the investigated sample. Despite these advantages, numerical evidence for the existence of a liquid–liquid critical point (LLCP) has been reported only for a small set of models1,11,12,13,14,15,16, and in some of the cases, not in a definitive way. It has even been speculated17 (but not confirmed by studies repeated by independent groups15,16,18) that the LLCP is a misinterpreted liquid–crystal transition, casting doubt on the relevance and appropriateness of this elegant hypothesis for interpreting the behaviour of water and silica, two of the most common substances on Earth.

If the liquid–liquid (LL) scenario has universal aspects, it should be possible to isolate its key ingredients in terms of the interparticle interactions to pinpoint the physical origin of the phenomenon. It should also be possible to clarify the contradictory results available in the literature and identify conditions in which the LLCP can be unambiguously observed in experiments. This requires a thorough investigation, not only of the liquid(s), but also of the relative stability of the disordered and ordered phases. Here, our investigation is inspired by two distinct recent theoretical discoveries: that network interpenetration can facilitate a LL transition19, and that bond directionality is crucial in controlling crystallization20. According to the first, particles with long arms and strong directional interactions can form fully bonded interpenetrating network structures separated by a first-order transition line. According to the second, when bonds are highly flexible, the liquid phase can remain stable in an intermediate range of densities down to very low temperatures. If the bond flexibility destabilizes the crystal state more rapidly than the LL transition, it might be possible to observe a genuine phase separation into two thermodynamically stable liquid phases, with no interference of crystallization at any temperature T.

Guided by the previously described discoveries, and mimicking the DNA tetravalent nanostars that have recently been experimentally investigated21,22,23, we develop a simple limited-valence model in which bond length and bond flexibility can be continuously varied. We attach to a central spherical core of diameter σ four identical spheres, generating four arms arranged in a tetrahedral geometry as depicted in Fig. 1. The centres of the four arms are located a distance L from the centre of the core sphere. The position of each of the arms in a tetramer is not fixed, but rather can deviate from its ‘ideal’ tetrahedral position by a maximum angle ϕ at no energetic cost, thus simply modelling arm flexibility. An attractive patch is located on each of the four arms, pointing directly away from the centre of the core. The patches interact through the commonly used angular square-well (of depth ε) Kern–Frenkel potential24 (Methods), with parameters (the angular width cosθ and the interaction range δ) satisfying the single bond per patch condition. When L = 0 and ϕ = 0, the model reduces to a sphere of diameter σ with four tetrahedrally oriented patches; for this case, previous studies have shown that the liquid crystallizes on cooling, without any evidence of a LLCP (ref. 25). In summary, the angle ϕ increases the bond flexibility beyond θ and L controls the effective ‘softness’ of the interaction. In this study we investigate more than 20 models, generated by different combinations of L and cosϕ. For each of these, we examine (using successive umbrella sampling (SUS) grand-canonical simulations26) approximately five to ten different T and two hundred different densities to provide an accurate description of the thermodynamic behaviour (Methods).

Figure 1: Tetramer model. Each tetramer consists of five hard spheres of diameter σ: a core at the centre and four ‘arms’ oriented along a tetrahedral geometry. The centres of the arms are located a distance L from the centre of the core. The vector connecting the centre of the core and each arm is allowed to fluctuate within a variable angle ϕ from its ideal tetrahedral position. In addition, there is an attractive patch on the surface of each arm characterized by an angular width cosθ = 0.95 and a bond range δ = 0.251σ, pointing directly away from the core. a, A two-dimensional schematic diagram of the model, showing a central sphere with three (out of four) arms, indicating the relevant angles and distances. b, The average configuration (solid) and typical fluctuations (partly transparent) for a tetramer with cosϕ = 0.9. Full size image

First, we investigate the role of L at fixed arm flexibility. At each state point, we record whether a LL phase separation occurs, and whether the system crystallizes spontaneously during the simulations. Typically, we very clearly observe a LLCP at sufficiently low T. The critical number density ρ c (L) is found to be approximately the same if measured in units of the typical distance between the centres of two bonded particles σ + L + δ/2, that is ρ c ∗ = ρ c (L)(σ + L + δ/2)3 ≈ 0.97. Thus, the free volume at ρ c increases as a function of L, affecting both the bonding probability at a given T and the ability of the system to form an interpenetrating network. As a result, the corresponding LL critical temperature T c depends significantly on the arm length. The results, summarized in Fig. 2, show a non-monotonous dependence of T c (solid line) on the arm length. For L → 0, the high packing fraction of the system near the LLCP hinders phase separation, as the hard cores of the particles make it difficult to compress part of the system into a denser liquid phase. In fact, it seems that in our model, a small amount of softness (L > 0) is necessary to observe a LL transition at finite T. At L ≈ 0.5σ, there is a maximum in T c beyond which the LL transition is shifted to smaller T. This is a trivial consequence of the decrease in ρ c σ3 as a function of L; that is, the bonds required for phase separation at larger L are formed at lower T because increasing the volume per particle disfavours bond formation.

Figure 2: Conditions for liquid–liquid (LL) phase separation: arm length. Diagram showing the region in which LL phase separation and crystallization occur as a function of arm length L and temperture T in units of the bonding energy (ε) for cosϕ = 0.9. The state points where LL phase separation was observed are denoted by red crosses; the blue dots indicate points where this was not the case. The solid line indicates the LL critical temperature T c , as determined from the successive umbrella sampling simulations. The dashed line shows the equilibrium crystallization temperature T x at ρ c , obtained from free-energy calculations performed for L/σ = 0.5,0.625,0.75,0.8375 and 1. The grey circles indicate spontaneous crystallization into a body-centred cubic crystal in the density region sampled by the LL critical fluctuations. Full size image

Figure 2 also shows the equilibrium crystallization temperature T x (L) along the critical isochore ρ c . This line intersects the LL critical-point line T c (L) around L = 0.82σ. The two lines separate the parameter space into four regions with different phase behaviour: the standard liquid region, present for all L at large T (denoted by I in Fig. 2); the region where the system crystallizes and there is no LL separation (II); the region where a LL phase separation exists, but is metastable with respect to crystallization (III); and the region where the LL is truly thermodynamically stable (IV).

We note that for small L we also observe spontaneous crystallization into a body-centred cubic (bcc) crystal near the LLCP. This is an effect of the large packing fraction sampled by the critical fluctuations near the LL phase transition, which favours the formation of the closely packed bcc crystal. As expected for this bond flexibility (cosϕ = 0.9), we never observe spontaneous crystallization of the lower density diamond crystal20,25. For comparison, we note that the typical (low-pressure) densities of diamond and bcc, when measured in the same reduced units as the critical density, are and respectively.

Next, we turn our attention to the effect of arm flexibility on the LL phase transition. Figure 3 shows that increasing the flexibility of the arms (that is, decreasing cosϕ) lowers T c (cosϕ). Stiff bonds thus also favour LL phase separation in addition to crystal stability20,25, explaining the difficulty in observing a stable LLCP in atomic and molecular network-forming systems, where interactions (for example, hydrogen bonds in water or sp3 electronic orbitals in silicon) are highly directional27. As in Fig. 2, the equilibrium crystallization temperature T x intersects T c , around cosϕ = 0.85, dividing the diagram into four regions. Again, we find a region (III) where the LL is metastable with respect to crystallization, as proposed for the case of water-like systems, and a region (IV) where the LL is truly thermodynamically stable. In this region, the LL critical point can be experimentally accessed, without any fear of crystallization.

Figure 3: Conditions for liquid–liquid (LL) phase separation: bond flexibility. Diagram showing the phase behaviour at ρ c (the LL critical density) in the bond flexibility (cosϕ) versus temperature (T) plane for L = 0.5σ. The indicated regions are the same as those shown in Fig. 2. The red crosses and blue dots indicate points where LL phase separation was or was not detected in the successive umbrella sampling simulations; the grey circles indicate spontaneous crystallization. The dashed and solid black lines are based on free-energy calculations performed at cosϕ = 0.8, 0.825, 0.8375, 0.85, 0.875, 0.9 and 1. The four vertical dotted lines denote the bond flexibility values for which the phase diagrams in Fig. 4 are drawn. ε is the bonding energy. Full size image

To better elucidate the effect of the arm flexibility on the topology of the phase diagram, we calculate the full phase diagrams, including ordered and disordered phases, in the T–ρ plane for four different values of cosϕ at fixed L. The different panels in Fig. 4 show the progression in the topology of the phase diagram on increasing flexibility. For highly flexible bonds, (Fig. 4a, cosϕ = 0.8), we observe a liquid–gas but no LL phase separation. In addition, there is a wide region of densities where the liquid is the stable phase for all T. For cosϕ = 0.825 (Fig. 4b), we observe both gas–liquid and LL critical points, but still no stable bcc or diamond crystal phase. This is the main finding of this work: the evidence that there are cases in which the LL transition can be studied without the interference of crystallization at any T. For cosϕ = 0.8375 (Fig. 4c), we still observe both the gas–liquid and the LL critical point. However, a triple point between a low-density liquid (LDL), high-density liquid (HDL) and bcc now exists: the LL phase separation becomes metastable with respect to a LDL–bcc phase separation at k B T/ε ≈ 0.057. We note that in all cases, the gas–liquid and LL phase separations occur at distinctly separated chemical potentials: a gas–HDL coexistence never occurs. Finally, for cosϕ = 0.9 (Fig. 4d), the LL phase separation becomes fully metastable with respect to crystallization. In addition, as a result of the strong bond directionality, the diamond phase becomes stable at intermediate density and low T (ref. 20). Interestingly, for all cosϕ values, the liquid side of the liquid–gas coexistence line ‘bulges out’ to higher densities, analogous to the phase diagram of water28.

Figure 4: Effects of bond flexibility on the phase behaviour. a–d, Phase diagrams for different values of the bond flexibility cosϕ, for L = 0.5σ. Symbols indicate state points where phase coexistences were calculated, with red symbols denoting metastable phase coexistences. The green regions indicate a stable liquid–liquid coexistence, and the red shaded regions in c and d show the metastable LL coexistence region. The dotted lines in c and d delimiting the body-centred cubic (bcc) region at high T indicate the expected behaviour of the fluid–bcc coexistence. Crystal structures with a higher density than the bcc phase were not investigated, but are not expected in the density range where the LL phase separation occurs. Note that in b no crystal phases compete with the LL phase separation, at any temperature. LDL: low-density liquid. HDL: high-density liquid. DC: diamond cubic. Full size image

Several important results stem from the extensive study reported in this Letter. In addition to the identification of two important control parameters (softness and flexibility), our study demonstrates the generality of the LL phenomenon in tetrahedral particles. Beyond a minimum flexibility (cosϕ) and a minimum softness (L) all models show a genuine LL critical point. Our study also shows that both the LL transition and crystallization are favoured by bond directionality, which destabilizes the homogeneous fluid at densities around ρ c by constraining the topology of the bonding pattern. As shown in Fig. 3, for short bonds and highly directional interactions (that is, for flexibility values typical of atomic and molecular systems) the competition between these phenomena is inevitably won by crystallization. This is the region of ϕ and L parameters typical of water, silica and silicon, and explains the exquisite sensitivity of the numerical results to the model parameters29. Variations in excluded-volume effects and/or directionality can significantly promote crystallization and suppress the possibility of approaching the LL critical point under metastable conditions. As shown in Figs 3 and 4, our calculations confirm that both for L → 0 and cosϕ → 1 (the values appropriate for the case of water, see ref. 27 and Supplementary Information), spontaneous crystallization prevails, consistent with experimental results. Nonetheless, our simple model shows that, depending on two key parameters, a LL phase transition is a general feature of tetrahedral liquids that has to be considered when studying such systems. Moreover, as speculated a long time ago in the case of water, the pattern of anomalies departing from the critical point affects the behaviour of the system even in the region where the fluid is stable8. These anomalies, for example, extrema in the density, isobaric heat capacity, and isothermal compressibility, occur in the present model as well, and are discussed in the Supplementary Information.

Finally, we have shown that sufficiently flexible bonds and sufficient softness give rise to a window where the full LL coexistence is thermodynamically stable. Interestingly, this range of flexibilities is within reach in the soft-matter realm. Patchy particles where the number and position of patches is fully controllable (including tetrahedral geometry)21,22,23 are nowadays experimentally realizable. These advances open up the possibility of synthesizing colloidal particles with the finely tuned softness and bonding flexibility required to experimentally explore the LLCP in the absence of crystallization. The new generation of patchy colloids may well provide the long-sought experimental evidence that a single-component liquid can phase separate into two distinct liquids.