Topological nature of the liquid–liquid phase transition in tetrahedral liquids

The first-order phase transition between two tetrahedral networks of different density—introduced as a hypothesis to account for the anomalous behaviour of certain thermodynamic properties of deeply supercooled water—has received strong support from a growing body of work in recent years. Here we show that this liquid–liquid phase transition in tetrahedral networks can be described as a transition between an unentangled, low-density liquid and an entangled, high-density liquid, the latter containing an ensemble of topologically complex motifs. We first reveal this distinction in a rationally designed colloidal analogue of water. We show that this colloidal water model displays the well-known water thermodynamic anomalies as well as a liquid–liquid critical point. We then investigate water, employing two widely used molecular models, to demonstrate that there is also a clear topological distinction between its two supercooled liquid networks, thereby establishing the generality of this observation, which might have far-reaching implications for understanding liquid–liquid phase transitions in tetrahedral liquids. Supercooled water undergoes a liquid–liquid phase transition. The authors show that the two phases have distinct hydrogen-bond networks, differing in their degree of entanglement, and thus the transition can be described by the topological changes of the network.

A salient feature of liquid water is the anomalous behaviour of its thermodynamic response functions upon cooling, the most famous being the density maximum at ambient pressure [1][2][3][4][5][6] . In an effort to explain the origin of water's thermodynamic anomalies, enhanced in supercooled states, Poole et al. suggested the presence of a first-order liquid-liquid phase transition (LLPT) line in the supercooled region of the pressure-temperature (P-T) phase diagram 7 . In this scenario, this transition line separates two liquid phases formed of transient hydrogen-bond networks-the low-density liquid (LDL) and high-density liquid (HDL)-and terminates at a liquid-liquid critical point (LLCP) 7 . The anomalous behaviour of liquid water, along with water polyamorphism, was then suggested to be a manifestation of the presence of this LLPT. Recent numerical studies 8,9 , as well as experiments [10][11][12] (despite the difficulties induced by the rapid formation of ice around the predicted critical temperature and pressure conditions), have provided strong support to this fascinating hypothesis.
Several order parameters, based on geometric or energetic criteria, have been proposed to characterize the transition [13][14][15][16][17][18][19][20] . In particular, the local structure index 13 and the r 5 parameter 14 -two order parameters that effectively measure the distance between the first and second coordination shells-have been routinely applied. However, along the years, it has become increasingly apparent that a proper description of the LLPT requires the development of order parameters which include information beyond that of the single-molecule neighbourhood. This led to the introduction of other order parameters, such as ζ (ref. 21 ), V 4 (ref. 15 ) and the total node communicability 20 , which include some information on the connectivity properties of the hydrogen-bond network 19 . It has also become clear that network interpenetration 22-24 plays a key role in the LLPT in tetrahedral liquids 23,25-27 , with the LDL being composed of a single random tetrahedral network and the HDL consisting of two (or more 23 ) disordered but locally interpenetrating networks. In this context, a structure-based order parameter that fundamentally distinguishes the two liquids at the microscopic level will be important to develop our understanding of the mechanism for achieving a difference in density and, hence, the physical origin of the LLPT.
A colloidal model of water. Colloidal patchy particles, due to their capacity to form network liquids via encoded self-assembly information 25 , combined with their synthetic availability 28-32 , are ideally suited to investigating universal aspects of the LLPT, as well as potentially facilitating detailed experimental investigations at single-particle resolution. Therefore, guided by recent advances in programming their hierarchical self-assembly 33 and in understanding the conditions necessary for facilitating a LLPT-bond flexibility 25,26 and network interpenetration 23,26 -we rationally designed a colloidal analogue of water. Figure 1a schematically illustrates this 'colloidal water' model. We considered designer triblock patchy particles 28,30,34 with two distinct attractive patches that encode tetrahedrality-key to the uniqueness of water 5 . The tetrahedrality emerges as the triblock patchy particles undergo two-stage assembly upon cooling. The energetics and geometry of the patches, labelled A and B, encode the desired staged-assembly information into the particles 33,35,36 . The energetics are chosen to facilitate the formation of discrete tetrahedral clusters, initially driven by strong A-A interactions, and then a tetrahedral network fluid via weaker B-B interactions ( Fig. 1a and Supplementary Fig. S1). The width of the A patches is judiciously chosen so as to favour the formation of a tetrahedral cluster fluid upon completion of the first stage of assembly 33,35,36 . The width of each B patch is instead chosen to favour its interaction with only one other B patch. The discrete tetrahedral clusters then behave as secondary building blocks, similar to tetrahedral patchy particles, that can interact with one another via B patches at lower temperatures to drive a second stage of assembly, leading to the formation of a liquid network. In this case, the patch widths are chosen not only to enforce the desired valency constraints but also to maximize the flexibility of the bonds to hinder the onset of crystallization 25 . As a result, the tetrahedral clusters formed by the triblock patchy particles can deviate from the ideal, tetrahedral symmetry. Note that the volume occupied by a tetrahedral cluster is less than that of a spherical patchy particle of equivalent size (Fig. 1a), thereby easing the formation of locally interpenetrating networks.

Topological nature of the liquid-liquid phase transition in tetrahedral liquids
Andreas Neophytou 1 ✉ , Dwaipayan Chakrabarti 1 ✉ and Francesco Sciortino 2 ✉ The first-order phase transition between two tetrahedral networks of different density-introduced as a hypothesis to account for the anomalous behaviour of certain thermodynamic properties of deeply supercooled water-has received strong support from a growing body of work in recent years. Here we show that this liquid-liquid phase transition in tetrahedral networks can be described as a transition between an unentangled, low-density liquid and an entangled, high-density liquid, the latter containing an ensemble of topologically complex motifs. We first reveal this distinction in a rationally designed colloidal analogue of water. We show that this colloidal water model displays the well-known water thermodynamic anomalies as well as a liquidliquid critical point. We then investigate water, employing two widely used molecular models, to demonstrate that there is also a clear topological distinction between its two supercooled liquid networks, thereby establishing the generality of this observation, which might have far-reaching implications for understanding liquid-liquid phase transitions in tetrahedral liquids.
Thermodynamic anomalies of colloidal water. To establish that the chosen system of triblock patchy particles can indeed be considered as a colloidal analogue of water, we first demonstrate that it captures water's thermodynamic anomalies. Figure 1b shows the reduced temperature (T ⋆ ≡ k B T/ϵ BB , where k B is the Boltzmann constant and ϵ BB is the well depth of patch B-patch B interactions) dependence of the reduced density (ρ ⋆ ≡ ρσ 3 , where σ is the hard-sphere diameter of a triblock patchy particle) along different isobars as determined by Monte Carlo simulations (Methods). Below a certain reduced pressure (P ⋆ ≡ Pσ 3 /ϵ BB ), the colloidal water model displays a density maximum upon cooling. Upon further cooling, ρ ⋆ approaches either a low or a high density value, depending on the pressure. Figure 1c shows the temperature dependence of the isobaric thermal expansion coefficient (α ⋆ P ), isothermal compressibility (κ ⋆ T ) and isobaric heat capacity (c ⋆ P ) at a select pressure, highlighting that, upon cooling, colloidal water shows non-monotonic behaviour in its thermodynamic response functions. Figure 1d shows the dependence of ρ ⋆ and the fraction of bonds f b (defined as the fraction of bonded B patches, a quantity directly reflective of the potential energy of the system) on P ⋆ . We find that, below a certain T ⋆ , there is a discontinuous jump in ρ ⋆ and f b as P ⋆ changes, signalling a first-order phase transition. In contrast, ρ ⋆ and f b change continuously with P ⋆ above this critical temperature, indicating that there is a LLCP. To locate the critical point, we followed a procedure similar to that applied recently in the study of a molecular model of water 9 to calculate the order parameter distribution P(M), where M = ρ ⋆ − sv ⋆ is a combination of the density, the potential energy per particle (v ⋆ = V/(Nε BB ), where V is the total potential energy) and a field-mixing parameter s (ref. 37 ). Upon performing a non-linear fit to identify the critical temperature (T c ) and critical pressure (P c ), we find that P(M) closely matches the distribution of the magnetization at the critical point of the three-dimensional (3D) Ising model (Fig. 1e), confirming the presence of an LLCP in colloidal water.
Topological nature of the LLPT in colloidal water. Topological concepts have become central in the description of physical 38-40 , biological 41,42 , social 43 and financial 44,45 networks. The topological properties of the interaction matrix, 39 owing to their unexpected connections with hidden conservation laws, have been shown to affect the response of the system to external perturbation and to play a key role in phase transitions. In this context, a pertinent question is whether the LLPT also has a topological character. Therefore, to identify any topological features, we inspected configurations of the LDL and HDL, finding that linked and knotted ring structures were noticeably present in the networks of the HDL phase but seemingly absent from the LDL phase. Figure 2a,b shows representative snapshots of the two liquid phases, where examples of these topologically complex motifs in the HDL networks are highlighted: a trefoil knot and a Hopf link 46 (Supplementary Video 1). In addition to knots and links, we found a large number of knotted theta curves 47,48 in the HDL networks, which can be viewed as two entangled fused rings (that is, rings that share at least two vertices). Examples of these knotted theta curves are shown in Fig. 2c and Supplementary Videos 2-4.
To investigate the importance of such topological features in the LLPT, we quantified the degree of 'entanglement' in the LDL and HDL networks using the concept of helicity. The helicity has previously been used to provide a measure of the topological properties of vector fields, including magnetic fields 49 and vortex fields in fluids 50 . This scalar quantity measures the amount of entwining in the set of closed curves contained within a system. For the particle-based networks considered here, these closed curves correspond to the set of rings formed by bonded particles. Each pair of disjoint (and therefore distinct) rings R i and R j contributes to the helicity with a term proportional, in essence, to how many times the two rings wind around one another (alternatively, this can be thought of as the number of times ring R i pierces the surface bounded by ring R j ). Also, each individual ring contributes with a term proportional to how many times it loops around itself. Formally, the helicity is defined as a double sum over the total number of rings N R of the double line integral, known as the Gauss linking integral: form stronger bonds than the B patches (blue) so as to encode two-stage assembly upon cooling. b, the evolution of the reduced density ρ ⋆ as a function of the reduced temperature T ⋆ for different reduced pressures P ⋆ , where P ⋆ × 10 3 = 5, 6, 7, 7.5, 8.5, 9, 10, 12, 14 and 16. the arrow indicates the direction of increasing P ⋆ . the inset highlights the density maximum for P ⋆ × 10 3 = 5, 6, 7 and 7.5. c, the evolution of the reduced thermal expansion coefficient (α ⋆ P ), isothermal compressibility (κ ⋆ T ) and isobaric heat capacity (c ⋆ P ) as functions of T ⋆ at P ⋆ = 0.0085 (close to the critical pressure). error bars represent the standard deviation over five sets of Monte Carlo trajectories, each of 1 × 10 8 cycles. d, the dependence of ρ ⋆ and the fraction of BB bonds formed (f b ) on P ⋆ at T ⋆ = 0.105 and T ⋆ = 0.112. e, the distribution of the order parameter M for colloidal water (blue symbols), calculated using histogram reweighing, with T ⋆ ≈ 0.1075, P ⋆ ≈ 0.0082 and s ≈ 0.627, compared with the corresponding 3D Ising universal distribution (solid red line).
where r i and r j are points on the two rings R i and R j , respectively. In equation (1), dr i and dr j are infinitesimal vectors tangential to the two rings R i and R j located at points r i and r j , respectively. In the case where i ≠ j, the Gauss linking integral gives the linking number (Lk ij ) of the two curves and takes on only integer values, whereas for i = j, it corresponds to the writhe (Wr i ) of a single curve and can take on any real value. Equation (1) can therefore be rewritten as the sum of the linking number and the writhe as A visual representation of the Gauss linking integral for a set of closed loops, along with their corresponding writhe or linking number, is shown in Fig. 2d.
To establish the prevalence of topologically complex motifs in a liquid network, we adapted the concept of the helicity to define two quantities: the network linking number 39 (Ln) and the network writhe (Wn). To emphasize the degree of entanglement in a network, Ln and Wn are computed by summing up the absolute values of Lk ij and Wr i , respectively, to ignore the intrinsic orientation of the closed curves (in other words, we ignore the chirality of the topological objects). Additionally, the writhe is only computed for the set N k of knotted objects that includes knotted rings and pairs of fused rings that contain a knotted path. Including the latter in N k ensures that the contributions of theta curves (such as those shown in Fig. 2c) are captured by Wn. We formally define Ln and Wn as We also note that, in some instances, the helicity can include a twist term when the structures under consideration are constructed from bundles of closed curves (as is the case for supercoiled DNA). However, since this is not the case with the bond rings considered here, the twist can practically be ignored.
We computed Ln and Wn by first identifying all the rings in the system up to a maximum size of lmax (which defines the maximum number of particles forming a ring) and then calculating either (1) the linking number for all disjoint rings or (2) the writhe for each ring that is knotted and for entangled fused rings (such as the knotted theta curves shown in Fig. 2c). This then becomes akin to  (1)). e, the dependence of the average network linking number (L n) and the network writhe (Wn) divided by the total number of rings in the network ( N R ) as a function of P ⋆ at T ⋆ = 0.105 (T < T c ) and T ⋆ = 0.112 (T > T c ). Ln and Wn were computed using rings of size up to lmax = 13. error bars represent the standard deviation over two Monte Carlo trajectories, each of 4 × 10 7 cycles. f, Frequency distribution for the writhe (Wr), separated for different unique knotted structures observed in the HDL network, at T ⋆ = 0.105 and P ⋆ = 0.016 with lmax = 15. the idealized knot associated with each distribution is colour coded and shown to the right of the plot. the dashed vertical lines indicate the writhe values of the idealized knots 57 , and the solid lines through the data points are a guide to the eye.
identifying linked and knotted paths in the spatial embedding of a graph 48,51 (full details on how Ln and Wn are computed in practice are given in Methods). Because we only identified links between rings of sizes l ≤ lmax and knotted structures of sizes l ≤ (2lmax − 2), the computed values of Ln and Wn depend on the chosen lmax. The LDL and HDL networks were found to show asymptotic limiting behaviour in their topological properties as a function of lmax, above a critical value (Supplementary Figs. S2 and S3). Therefore, we computed Ln and Wn using the smallest lmax value that was found to be in the asymptotic regions of both the LDL and HDL phases (lmax ≥ 13). Figure 2e shows the average values of Ln and Wn as a function of pressure for T < T c and T > T c with lmax = 13. Both Ln and Wn display a discontinuous change with pressure for sub-critical temperatures (T < T c ), concomitant with the discontinuity in the density (Fig. 1d), while they change continuously for super-critical temperatures.
The transition between the LDL and HDL phases is clearly captured by Ln when lmax ≥ 5 (Extended Data Fig. 1). However, as knots require longer chains of particles to form, they require a larger lmax to be identified, and so, Wn only shows a clear discontinuity when lmax ≥ 10 (Extended Data Fig. 1). Additionally, this means that there are fewer knots than links in the network, which in turn, means that Wn is smaller than Ln for a given lmax. Yet, there are a variety of knotted objects in the HDL networks that contribute to the nonzero value of Wn. Figure 2f shows the frequency distribution of the writhe for different knotted objects found in the HDL network with lmax = 15. The majority of these structures contain trefoil knots, but we also find figure-of-eight, cinquefoil, three-twist and Stevedore knots 46 . Similar to randomly generated knotted rings, the ensemble average of the writhe for each of these knots correlates well with the writhe of the corresponding idealized configuration (an idealized knot configuration is one which allows the volume-to-surface area ratio to be maximized 52 ). hydrogen-bond networks of molecular water, simulated using the tIP4P/Ice potential, at T = 188K and P = 1,100 bar and P = 2,000 bar, respectively. this temperature is slightly below the critical temperature (T c = 188.8 K (ref. 9 )), and the pressures are below and above the critical pressure (P c = 1,725 bar (ref. 9 )). In the LDL, only a few links are present (one is highlighted in a). In the HDL, both linked and knotted motifs are present (highlighted in b). c, the frequency distribution of the writhe (Wr) of the knotted paths found in the HDL network (essentially all being trefoil knots) at T = 188 K and P = 2,000 bar. the inset shows the idealized configuration of a trefoil knot, and the dotted vertical line indicates the value of the writhe for the idealized trefoil knot 57 . d, Dependence of Ln and Wn, computed using rings of sizes up to lmax = 13, divided by N R , as a function of P at T = 188 K for N = 1,000 molecules. error bars represent the standard deviation as calculated along the corresponding molecular dynamics trajectory. e, the fluctuations in the density (ρ) and Ln (computed using lmax = 13) with time (t) along an isobaric-isothermal molecular dynamics trajectory for N = 300 molecules of tIP4P/Ice water at a temperature of T = 188 K for P = 1,100 bar, P = 1,725 bar and P = 2,000 bar.
Topological nature of the LLPT in molecular water. It is clear that both Ln and Wn can serve as a topological order parameter to describe the LLPT in colloidal water, providing strong evidence that the LLPT is a transition from an 'unentangled' to an 'entangled' network. This raises the question of whether entanglement also underpins the LLPT in supercooled water. Recently, it was shown that both the TIP4P/Ice 53 and TIP4P/2005 54 molecular models of water, which do not explicitly enforce tetrahedral coordination, exhibit a first-order LLPT and corresponding LLCP 9 . We therefore performed a similar topological analysis for TIP4P/Ice and TIP4P/2005. Figure  3a,b shows representative snapshots of the hydrogen-bond networks for the LDL and HDL of TIP4P/Ice. We find that the LDL contains a small number of links and practically no knotted structures (with only one or two knots occasionally appearing in a frame), while the HDL contains trefoil knots and theta curves in addition to links (which are also present in more abundance than in the LDL). Figure  3c shows that the ensemble average of the writhe for these trefoil motifs in the HDL of TIP4P/Ice also correlates well with the writhe of the idealized trefoil knot, as was the case for colloidal water. Additionally, as for colloidal water, Fig. 3d shows that, for TIP4P/ Ice, both Ln and Wn jump at the critical pressure, suggesting that entanglement also underpins the LLPT in molecular water. A universal property of systems close to a second-order critical point (when T ≈ T c and P ≈ P c ) is wide oscillations in relevant order parameters, as configurations belonging to both phases are sampled along the trajectory 37 . To check whether the topological order parameter Ln does indeed properly distinguish between the two liquid phases of the TIP4P/Ice and TIP4P/2005 models, we calculated Ln along 40 μs molecular dynamics trajectories of supercooled water for different values of P close to T c . Figure 3e shows the evolution of ρ and Ln along the molecular dynamics trajectory for a system of N = 300 TIP4P/Ice molecules at three different pressures. Above and below the critical pressure, Ln displays quite different values, indicating a notable topological difference between the two liquids. At the critical pressure, wide fluctuations in Ln are observed, closely correlating with those of ρ. Extended Data Fig. 2 shows that the TIP4P/2005 model also displays these wide fluctuations in Ln and ρ when P ≈ P c , confirming that the LLPT can indeed be interpreted as a topological transition.

Conclusions.
We have designed a colloidal analogue of water-a tetrahedral network liquid self-assembled from designer triblock patchy colloidal particles via tetrahedral clusters. This colloidal water model captures the anomalous thermodynamic behaviour of supercooled water, including the well-known density maximum, and displays two structurally different network liquid phases. With the synthetic feasibility of the triblock patchy particles under consideration 30 , colloidal water presents itself as an ideal model system, amenable to experimental investigation at the single-particle level, for gleaning fundamental understanding of LLPTs in tetrahedral liquids.
By analysing configurations from this colloidal model and from two widely used molecular models 53 for water, we have established that the LDL and HDL networks associated with the LLPT in all three cases are topologically distinguishable. The LDL in all these cases is practically an 'unentangled network' , consisting primarily of unknots, while the HDL is an 'entangled network' containing an ensemble of topologically non-trivial motifs. For colloidal as well as molecular models, Ln and Wn-metrics able to capture the degree of entanglement-act as suitable order parameters, highlighting that the LLPT can also be interpreted as a topological phase transition.
Uncovering the topological distinction between the LDL and HDL networks allows us to understand the physical mechanism underpinning the LLPT from a new perspective that sheds light on its microscopic origin. The LLPT occurs at low temperatures at which the system tends to form a liquid network with the maximum number of bonds possible. At low pressures, it is well established that this is achieved by forming an open random tetrahedral network. However, as the pressure increases, the system tends also to minimize its volume. In tetrahedral liquids, associated with this drive to maximize the number of bonds, there is an increase in the number of rings in the network. Due to the directionality of the interactions, the system cannot minimize its volume by simply deforming the rings while preserving the network connectivity because there is a limit to the flexibility of the bonds. Instead, the system is able to simultaneously minimize its volume and maximize the number of bonds in the network by forming knots and links. Additionally, we note that, unlike knots and links formed by covalently bonded chains or networks 40 , the LDL and HDL are transient networks (that is, bonds are constantly breaking and forming) and so their topological state is not constant. Despite this evolution, the topological signatures for each phase are properly defined.
The helicity has been used previously to study knotted vortices in classical fluids 55 and superfluids 56 , where these knots were found to untie via a universal topological mechanism. It would be of interest to use the metrics for the degree of entanglement used here to follow the nucleation of the LDL out of the HDL, and vice versa, to see whether this universal topological mechanism is also relevant in the context of the LLPT. Finally, owing to their topological complexity, these entangled liquids should possess physical properties 38,40 distinct from simple liquids, which will merit exploration.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41567-022-01698-6.

Methods
Colloidal water model. We employed the Kern-Frenkel pair potential 58,59 , which has been used extensively to study patchy particles. Each particle, with a hard-core diameter σ, has (at opposite poles) two distinct, circular attractive patches, labelled A and B, having half-angles of θ A and θ B , respectively. The effective potential for a pair of patchy particles v ij is vij(rij, Ωi, Ωj) = v hs where r ij = |r ij | is the centre-to-centre distance between particles i and j and Ω i and Ω j describe the orientations of particles i and j, respectively. In equation (4), v hs ij is the hard-sphere pair potential, given by v hs and v sw ij is a square-well potential, given by v sw where ε αβ and δ αβ control the depth and range, respectively, of patch α-patch β attraction. The factor f(rij,n α i ,n β j ) controls the angular dependence of the interaction between two patches and is given by where n α i is a normalized vector from the centre of particle i in the direction of the centre of patch α on its surface, thus depending on Ω i . Similarly, n β j is a normalized vector from the centre of particle j in the direction of the centre of patch β, thus depending on Ω j . The total potential energy of the system was calculated by summing over the contributions from all distinct pairs (V = ∑ i̸ =j vij). We employed ε AA = 5ε BB along with ε AB = 0 to encode a two-stage self-assembly process upon cooling 33 . We selected θ A = 50° and θ B = 26° because this choice, combined with the interaction strengths, ensures that an A patch interacts with three other A patches (to yield discrete tetrahedral clusters at relatively high temperatures) while a B patch interacts with only one other B patch (to bring the tetrahedral clusters together to form a tetrahedral network liquid at lower temperatures). The choice of ε AB = 0, which is realizable via DNA-mediated interactions, facilitates staged assembly and allows for the use of a wider A patch without compromising the formation of self-limiting tetrahedral clusters at intermediate temperatures. Note that the use of an A patch with θ A = 50° makes the tetrahedral clusters flexible enough to hinder crystallization from the tetrahedral liquids formed at lower temperatures. The ranges of the interactions are set by δ AA = 0.05 and δ BB = 0.2 (the value of δ AB being irrelevant since ε AB = 0). The longer range of B-B interactions facilitates the formation of a tetrahedral network liquid at low temperatures. In the present study, we used reduced units: the length in units of σ, the energy in units of ε BB and the temperature in units of ε BB /k B , with the Boltzmann constant set to k B = 1.
Monte Carlo simulations of colloidal water. All the Monte Carlo simulations were carried out with systems of N = 1,000 triblock patchy particles contained in a cubic box under periodic boundary conditions, using the minimum image convention. Each triblock patchy particle was treated as a rigid body whose orientational degrees of freedom were represented by quaternions. The potential energy was calculated using a spherical cutoff of σ(1 + δ BB ), and a cell list was used for efficiency 60 . To enhance sampling in systems composed of clusters, we employed a combination of translational and rotational single-particle and cluster moves, where the acceptance probability for a translational or rotational cluster move was taken as 60 where β = 1/k B T, V(o) and V(n) correspond to the potential energy of the system before and after the cluster move, respectively, and P(i, j) represents the probability that particles i and j share a bond, that is, that the patches interact. We defined a cycle as the number of attempts, equal to the number of particles N, to move (rotate or translate) a particle or a cluster (with equal probability). Additionally, in isothermal-isobaric (NPT) simulations, we included one volume move, on average, per cycle. We performed cluster volume moves in the logarithmic space of the volume (that is, ln(Vn) = ln(Vo) + ΔV ), where each cluster's centre of mass was scaled isotropically. Volume moves were accepted with probability 61 where N clstr is the number of discrete clusters identified in the system. In the context of an anticipated two-stage assembly of the patchy colloidal particles via a discrete set of clusters, formed by A-A bonds, we constructed a simple rule to move the discrete clusters. We took P(i, j) = 1 if two particles share an A-A bond, and P(i, j) = 0 otherwise. Therefore, to ensure detailed balance is satisfied, any proposed cluster or volume moves that created a new A-A bond between particles i and j were rejected. At low temperatures (T ⋆ ≤ 0.12), two sets of simulations were carried out over the same range of state points, starting from a tetrahedral cluster fluid of triblock patchy particles, where the initial density was set to ρ ⋆ = 0.2 for one set of simulations and to ρ ⋆ = 0.4 for the other set. For each of these simulations, 3 × 10 8 Monte Carlo cycles were performed to ensure that equilibrium was attained. At higher temperatures, 1 × 10 8 Monte Carlo cycles were performed, starting with a density of ρ ⋆ = 0.4.
From the NPT simulations, we calculated the reduced isobaric thermal expansion coefficient (α ⋆ P ), isothermal compressibility (κ ⋆ T ) and isobaric heat capacity (c ⋆ P ) by computing the covariance and variance in the volume and enthalpy (H) as These thermodynamic quantities, as well as the average density, were computed from five independent NPT simulations, consisting of 1 × 10 8 Monte Carlo cycles, following equilibration at each temperature and pressure considered.
Molecular dynamics simulations of TIP4P/Ice and TIP4P/2005. Molecular dynamics simulations were performed for the TIP4P/Ice 53 and TIP4P/2005 54 water models in the NPT ensemble, using GROMACS. 62 The integration time step was 2 fs. The Nosé-Hoover thermostat and the isotropic Parrinello-Rahman barostat, both with characteristic times of a few picoseconds, were used. Rigid constraints for the molecular model were implemented with a sixth-order linear constraint solver. The long-range electrostatic interactions were dealt with by using the particle-mesh Ewald method at fourth order. The cut-off distances for both the van der Waals and the real-space electrostatic interactions were fixed at 0.9 nm. For the TIP4P/Ice model, we focused on the T = 188 K isotherm for systems of N = 300 and 1,000 molecules, with pressures ranging from 1 bar to 4,000 bar. For the TIP4P/2005 model, we performed simulations at T = 177 K and P = 1,750 bar for a system of N = 300 molecules. The analysed trajectories were longer than 40 μs for both models, and equilibrated data from ref. 9 were used as starting configurations.
Cluster identification in the colloidal water model. We identified the tetrahedral clusters formed by the triblock patchy particles using the local order parameter, q Td (which equals 1 for a perfect tetrahedron), defined as 63 where N c = 4 is the number of particles in the cluster under consideration and ψ ij is the angle subtended at the centre of the cluster by the two vectors joining the centre of particle i and particle j. We defined the set of neighbours for each particle i as those with which it shares a patch A-patch A bond. q Td was then calculated for all distinct combinations of four particles including i from this set. All unique sets of four particles satisfying the condition q Td > 0.99 were considered to form a tetrahedron.

Multiple histogram reweighting.
To identify the location of the LLCP in the colloidal water system, we merged information collected from the Monte Carlo simulations to generate an estimate of the density of states, Ω (E, V). In the NPT ensemble, the partition function at a given pressure P and temperature T can be written as 9,64 where E is the energy of the system, V is the volume and G is the Gibbs free energy. We estimated the density of states by performing R simulations at different temperatures and pressures, where the total set of state points considered was given by {T i , P i }, i = 1, …, R. For each simulation, a histogram (H i ) of the energies and volumes sampled was generated using data from C i configurations spaced equally along the trajectory. We defined an estimate for the density of states at discrete values of energy and volume, using data from a single simulation, as However, each of the R simulations only provides information of the density of states for a finite portion of the phase space. Hence, a better estimate was obtained by combining all the information gathered from the R simulations. This 'best guess' for the density of states was written as a linear combination of these individual estimates, combined with a weight factor ω i (E, V), aŝ Conventionally, it is assumed that the 'random' variable H i (E, V) has Poisson statistics. It then becomes possible to optimize equation (16) by minimizing the variance of Ω b (E, V) using the method of Lagrange multipliers, under the constraint that and When the number of samples used to generate each of the R histograms is the same (that is, ∀i: C i = C), the above equations can be simplified by writing and ) . (20) Evaluation of the exponential terms can be numerically unstable, but rewriting these equations in terms of the entropy allows us to make use of ln(e x ) = a + ln(e x−a ) to improve the numerical precision. Equations (19) and (20) can then be solved self-consistently via an iterative procedure, taking G i = 0 for all i as an initial guess. We assume that the values of the Gibbs free energies have converged once the sum squared difference between consecutive estimates is less than a tolerance value of 10 −7 . The converged solutions provided an estimate for the values of G(P, β) up to some additive constant c. From these estimates, we obtained values for the density of states. It was, therefore, possible to predict the probability of observing a state with energy E and volume V at any relevant state point (P, β) using Bond definition. In the case of colloidal water, the network is created by the patch B-patch B interactions, forming bonds that join discrete tetrahedral clusters.
Owing to the square-well nature of the patch-patch interaction, searching for the existing bonds becomes equivalent to searching for all pairs of particles with pair interaction energy equal to −ϵ BB . For the TIP4P/Ice and TIP4P/2005 models, a bond between two water molecules was deemed to exist when a hydrogen bond was present. This was defined on the basis of a geometric criterion 65 by the conditions r < 0.35 nm and θ < 30°. Here, r is the intermolecular oxygen-oxygen distance and θ is the smallest among the four angles between the intramolecular O-H and the intermolecular O-O lines. As shown previously 19 , no ambiguity in the identification of hydrogen bonds remains if the Chandler-Luzar geometric criterion is applied using inherent structure configurations. These are the local potential energy minimum (21) configurations reached via a steepest descent path that quenches the vibrational degrees of freedom. To evaluate the inherent structures, we used the steepest descent algorithm in GROMACS (compiled in double precision) with a force tolerance of 1 J mol −1 nm −1 and a maximum step size of 5 × 10 −4 nm.

Ring statistics.
We calculated the number of rings ( NR l ) with sizes l ∈ [4 − lmax] where each bond network was abstracted into a periodic undirected graph G. The vertices (VG) of the graph represent the positions of the objects of interest (either the tetrahedral cluster centres or the oxygen atoms of the water molecules), and each of the edges (EG) connects two bonded objects as defined above. We defined a ring as a closed path in this graph (that is, a path whose first and final vertices are the same). We excluded rings in which non-adjacent vertices are directly connected via an edge. For each vertex V i G in the graph, we determined all distinct rings, sequentially identifying all rings associated with each vertex, using a depth-first search traversal, up to the chosen lmax. To improve the efficiency of the analysis, following the extraction of all relevant paths starting from V i G (and hence rings containing V i G ), all edges connected to that vertex were removed from the graph before initiating the depth-first search traversal from V i+1 G . We also accounted for the double counting of a ring, which arises from the undirected nature of the graph.
Linking number and writhe calculations. We measured the entanglement of liquid networks using metrics derived from the helicity, which is defined in equation (2), in terms of Lk ij , the linking number between two independent closed curves R i and R j , and Wr i , the writhe of a single closed curve. Both the linking number and the writhe can be defined as double line integrals 46,66-68 : Note that equation (22) and (23) are combined together in equation (1) (the contributions i ≠ j and i = j, respectively). The closed curves from which we computed Ln and Wn are constructed from the set of rings identified in the liquid networks from both models, where r i and r j represent the positions of the objects that make up rings R i and R j , respectively. We calculated the linking number only for pairs of disjoint rings (rings not sharing any vertices), and the writhe for both knotted individual rings and for pairs of rings that share at least one edge that were also knotted (these pairs of rings can essentially be considered as composite rings with intra-ring edges).
Computing the linking number. The rings are composed of discrete line segments and so can be considered as closed polygonal paths. As a result, the linking number can be computed using the expression for polylines 69 : where λ kl is the contribution to the Lk ij from a pair of line segments k and l belonging to rings i and j, respectively, and λ kl = atan where a = r l − r k , b = r l − r k+1 , c = r l+1 − r k+1 and d = r l+1 − r k (ref. 69 ).
Computing the writhe. The writhe has previously been used to study folding in protein systems 70 and supercoiling in DNA 68 . In these studies, proteins and DNA strands are modelled as polymer chains made up of straight line segments. Similarly, we computed the writhe using 68 where ω kl is the contribution to the writhe from the segments k and l on ring i, defined as 68 ω kl = sgn { [(r l+1 − r l ) × (r k+1 − r k )] · a } ( sin −1 (n1 · n2) + sin −1 (n2 · n3) + sin −1 (n3 · n4) + sin −1 (n4 · n1) ) , where sgn(⋅) is the sign function and where a, b, c and d are defined as above.