in our model , each particle consists of a core sphere with four tetrahedrally oriented arms ( see fig . the central sphere and the arms are modeled as interpenetrating hard spheres with diameter  , with the centre of each arm located a distance l from the centre of the core . on the surface of each arm the patch is modeled , via a standard kern - frenkel potential , as a square - well type attraction of depth  : two particles are either bound or unbound , with no intermediate energy levels . the opening angle of the patch is given by cos  = 0.95 , and the maximum interaction range  = 0.251. this choice ensures that each patch can form a bond with only a single other patch . furthermore , in order to model arm flexibility , each arm can freely rotate around the core sphere within a maximum angle  from its ideal ( tetrahedral ) position . note that there are no hard - core interactions within a particle : spheres belonging to the same tetramer can overlap . each mc simulation included translation moves , rotation moves rotating a full tetramer , rotation moves rotating only the core sphere , and rotation moves that rotate a single arm around the core . additionally , we added volume moves for isobaric simulations , and particle insertion and deletion moves for grand - canonical successive umbrella sampling ( sus )   simulations . to detect ll phase separation , we used sus simulations , generating the relative free energy as a function of the density . to calculate the phase diagrams in fig . 4 , and the crystallization lines in figs . 2 and 3 for the liquid phases , we used a high temperature fluid as a reference state , where we obtained the helmholtz free energy f by combining the chemical potential  ( evaluated from a sus simulation ) and the pressure p at the same density  ( taken from an npt mc simulation ) , using : 
 ( 1)ffn=p , 
 where n is the number of particles . we then used thermodynamic integration by integrating the temperature along an isochore : 
 ( 2)2f(,t2)=1f(,t1)+12du(,t ) , 
 with u( , t ) the average potential energy per particle ,  = 1/kbt and kb boltzmann s constant . typically , at sufficiently low temperatures , the energy as a function of temperature takes the form 
 ( 3)u=2+cexp(2 ) , 
 where 2 is the ground - state energy ( where each particle has four bonds ) and c is a ( density - dependent ) constant . this allows for straightforward extrapolation of the free energies to arbitrarily low temperatures after the regime in eq . this regime is accessible in our simulations for all cases where the fluid does not crystallize at low t ( i.e. cos   0.9 ) . system sizes for the sus simulations were chosen such that near the ll transition the system contained n  150 particles . other free - energy calculations were performed with similar or larger n. the investigated crystal structures were diamond and bcc ( corresponding to ice ic and ice vii in water , respectively ) . both are fully bonded structures and as such have the ( same ) lowest possible potential energy . depending on cos  and l , there will be a variety of other fully bonded crystal structures . however , the structures commensurate with an ideal tetrahedral geometry will always have the highest entropy and hence will be more stable at the densities where they are accessible . thus , the only crystals which might compete with the ll phase separation are diamond and bcc . at high density , i.e. beyond the bcc density , the stable crystalline phases will be controlled by a competition between packing , energy , and entropy , but as these phases will not change the stability of the ll phase transition , they are outside of the scope of this letter . similarly , we do not expect these crystal phases to affect the phase diagrams in fig . 4 , with the possible exception of the ( dashed ) high - density liquid - bcc lines in panels ( c ) and ( d ) . to determine free energies for the crystal phases , we used an einstein crystal approach , for each t where a phase coexistence was calculated . subsequently , using npt mc simulations , we employed thermodynamic integration along the density at constant temperature : 
 ( 4)f(,t)=f(ref , t)+refd()p2 , 
 where ref is the density where the reference free energy was calculated . s4 in the si for examples ) . to determine whether spontaneous crystallization occured in our sus simulations , we used a bond - order parameter to find the largest crystalline cluster in the system . first , we determine the neighbors of each particle using a solid - angle based nearest - neighbor method . we then calculate for each particle the complex vector ql , the expansion of their environment in terms of spherical harmonics ylm , with l the order of the symmetry of the crystal , and l  m  l : 
 ( 5)qlm(i)=1nb(i)j=1nbylm(r^ij ) . 
 here , nb is the number of neighbors of particle i , and r^ij is the normalized vector connecting particle i to particle j. for the purpose of finding crystalline clusters , two particles are considered bonded when their environments are sufficiently similar . to determine this , we calculate : 
 ( 6)dl(i , j)=re(ql(i)ql(j)ql(i)ql(j) ) . 
 to detect bcc clusters , we choose l = 6 , consider two particles bonded if dl(i , j ) > 0.6 , and label any particle that has at least 5 bonds as a crystalline particle . for diamond , two particles are considered bonded if either dl(i , j ) < 0.87 or 0.3 < dl(i , j ) < 0.1 , and only particles with 4 bonds are considered crystalline . systems were considered to have crystallized spontaneously when more than 10% of the system was crystalline . note that although the q6 based order parameter is sensitive to multiple crystal structures , visual inspection of the detected clusters consistently showed a bcc structure in all cases .
one of the most controversial hypotheses for explaining the origin of the thermodynamic anomalies characterizing liquid water postulates the presence of a metastable second - order liquid - liquid critical point [ 1 ] located in the  no - man s land  [ 2 ] . in this scenario , 
 two liquids with distinct local structure emerge near the critical temperature . unfortunately , since spontaneous crystallization is rapid in this region , experimental support for this hypothesis relies on significant extrapolations , either from the metastable liquid or from amorphous solid water [ 3 , 4 ] . 
 although the liquid - liquid transition is expected to feature in many tetrahedrally coordinated liquids , including silicon [ 5 ] , carbon [ 6 ] and silica , even numerical studies of atomic and molecular models have been unable to conclusively prove the existence of this transition . 
 here we provide such evidence for a model in which it is possible to continuously tune the softness of the interparticle interaction and the flexibility of the bonds , the key ingredients controlling the existence of the critical point . 
 we show that conditions exist where the full coexistence is thermodynamically stable with respect to crystallization . 
 our work offers a basis for designing colloidal analogues of water exhibiting liquid - liquid transitions in equilibrium , opening the way for experimental confirmation of the original hypothesis .