{'fraction_non_alphanumeric': 0.03801533903701858, 'fraction_numerical': 0.008998532547000027, 'mean_word_length': 4.737569499602859, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cluster-variation approximation for a network-forming lattice-fluid model

We consider a 3-dimensional lattice model of a network-forming fluid, which has been recently investigated by Girardi and coworkers by means of Monte Carlo simulations [J. Chem. Phys. \textbf{126}, 064503 (2007)], with the aim of describing water anomalies. We develop an approximate semi-analytical calculation, based on a cluster-variation technique, which turns out to reproduce almost quantitatively different thermodynamic properties and phase transitions determined by the Monte Carlo method. Nevertheless, our calculation points out the existence of two different phases characterized by long-range orientational order, and of critical transitions between them and to a high-temperature orientationally-disordered phase. Also, the existence of such critical lines allows us to explain certain ``kinks'' in the isotherms and isobars determined by the Monte Carlo analysis. The picture of the phase diagram becomes much more complex and richer, though unfortunately less suitable to describe real water.


I. INTRODUCTION
Water attracts great interest from physicists, due to its different anomalous properties [1,2,3]. Just to mention a few of them, it is well known that, at ordinary pressures, the solid phase (ice) is less dense than the corresponding liquid phase, and the latter displays a temperature of maximum density, just above the freezing transition. Furthermore, the heat capacity of liquid water is unusually large, whereas both isothermal compressibility and isobaric heat capacity display a minimum as a function of temperature. Although a full prediction of such and other anomalies from first principles has not been given yet, it is widely believed that it should be related to the ability of water molecules to form a network of hydrogen bonds.
In the early model proposed by Bell [10,11,12], molecules can point their arms only to nearest neighbor sites. An attractive energy is assigned for every pair of occupied nearest neighbors, with an extra contribution if a hydrogen bond is formed, i.e., if a donor arm points to an acceptor arm. Furthermore, a repulsive energy is assigned for certain triplets of occupied sites, in order to account for the difficulty of forming H bonds by closely packed water molecules. Minor variations of the Bell model have been proposed by Meijer and coworkers [13,14], which replaced the three-body interaction with a simple next-nearest-neighbor repulsion, and by Lavis and Southern [15], who defined a simplified model with no distinction between donors and acceptors, first pointing out that such a distinction is scarcely important for a qualitative description of water thermodynamics. All these studies are mainly focused on the equilibrium phase diagram of water, and predict a liquid-vapor coexistence and two different low-temperature phases characterized by long-range orientational order and different densities. At zero temperature, the low-density phase is an ideal diamond network made up of H bonded water molecules, with half the lattice sites left empty, resembling the structure of ice Ic (cubic ice). At the melting temperature, this phase correctly displays a lower density than the liquid phase. Conversely, the high-density phase at zero temperature is made up of two interpenetrating diamond structures, with all sites occupied, resembling the structure of ice VII (a high pressure form of ice).
More recently, a similar model has been proposed by Roberts and Debenedetti [18,19,20], with the aim of exploring the phase diagram of water deep into the supercooled liquid region, and, more generally, the possibility of liquid-liquid immiscibility in a network-forming fluid. There are two differences with respect to the original Bell model. First, water molecules have an additional number of configurations, in which they cannot form H-bonds. Such a number is an adjustable parameter, related to the entropy of breaking H bonds. Second, the energy penalty for the triplets occurs only when two molecules in a triplet form a H bond. At zero temperature, this model still predicts the two ordered H-bond networks of the Bell model, but these phases have never been investigated at finite temperature, due to the main attention devoted by the cited papers to (stable or metastable) liquid phases. In particular, these studies have given a contribution in the framework of a debate between two different conjectures put forth respectively by Speedy and Angell [29] and Poole and coworkers [30] to explain apparent power-law divergences in the response functions of supercooled water [29]. In a few words, Speedy and Angell postulated a reentrance of the stability limit (spinodal) of liquid water, whereas Poole and coworkers postulated the existence of two different metastable liquid phases (in analogy with the two different amorphous ices observed in experiments), whose coexistence line could possibly end in a (metastable) critical point. The model by Roberts and Debenedetti turned out to be compatible with the second critical point conjecture, which has also collected several other indirect evidences, both in experiments [31] and computer simulations [32]. Two authors of the present paper have studied a simplified version of the Roberts-Debenedetti model, without the donor-acceptor asymmetry [21,22], pointing out that the model is indeed compatible with both the aforementioned conjectures (depending on parameter values), and is also able to describe the anomalous properties of water as a solvent for apolar molecules (hydrophobic effect). A different variation of the Bell model has been considered by Besseling and Lyklema [16,17], who have taken into account only nearest-neighbor interactions, namely, an attractive term for H-bonded molecules and a repulsive term for nonbonded molecules. Even in this case, the two different ordered phases described above are stable at zero temperature, but they have not been investigated at finite temperature. The cited studies were indeed focused on the liquid-vapor interface properties [16] and on hydrophobic hydration thermodynamics [17], for which the authors obtained good agreement with experiments. These results were obtained in the so-called quasichemical or Bethe approximation, i.e., a first-order approximation which takes into account correlations over clusters made up of two nearest-neighbor sites.
Very recently, Girardi and coworkers [23] have performed extensive Monte Carlo simulations for the previously described model, in the simplified version obtained by removing the donor-acceptor distinction. This work suggests the existence of two liquid phases of different densities, which can coexist in equilibrium, in agreement with the conjecture of Poole and coworkers [30]. The coexistence line terminates in a critical point, whereas the low-density phase exhibits a temperature of maximum density (TMD), depending on pressure. A subsequent work [24] suggests a relationship between the density anomaly and the anomalous behavior of the diffusion coefficient.
In this paper we analyze the above model by a generalized first-order approximation, based on a 4-site tetrahedral cluster. Such a work is motivated by the fact that the reliability of cluster-variational approximations, though widely employed for this kind of models with orientationdependent interactions, has been scarcely tested against Monte Carlo simulations [33]. Let us anticipate that our calculation reproduces most physical properties obtained by simulations (namely, isotherms, isobars, binodals, TMD locus) with remarkable accuracy. Nevertheless, the picture of the phase diagram turns out to be quite different, as the two different condensed phases exhibit long-range orientational order, which lead us to identify them as a continuous temperature evolution of the two zero-temperature network structures introduced above. As a consequence, the two critical points are indeed tricritical, and two critical lines appear, indicating two different kinds of symmetry breaking, respectively from an orientationally-disordered phase to the highdensity ordered phase and from the latter to the lowdensity one. In order to check the existence of critical lines and symmetry breaking we have performed some new Monte Carlo simulations in specific regions of the phase diagram.
The paper is organized as follows. In Section II we describe the model hamiltonian and analyze the groundstate. In Sec. III we introduce the cluster-variational technique employed for the calculation. Sec. IV describes our results, discussing them in comparison with those obtained from numerical simulations in Ref. [23]. Sec. V reports on the new Monte Carlo results, while Sec. VI is devoted to some concluding remarks.

II. MODEL AND GROUND STATE
As mentioned in the Introduction, the model is defined on a body-centered cubic lattice (see Fig. 1). Each site can be empty or occupied by a molecule. The model molecule possesses the tetrahedral structure introduced above, with 4 equivalent bonding arms, which can point towards 4 out of 8 nearest neighbors of a given site. A hydrogen bond is formed, yielding an attractive energy −η < 0, whenever two nearest-neighbor molecules point an arm to each other, with no distinction between donor and acceptor (see Fig. 2). Moreover, a repulsive interaction ǫ > 0 is assigned to any pair of nearest-neighbor sites occupied by water molecules. This choice of the interaction parameters implies an energetic penalty for nearestneighbor molecules not forming H bonds. Finally, as we find it convenient to work in the grand-canonical ensemble, a chemical potential contribution −µ is taken into account for every occupied site.
In principle, the hamiltonian of the system can be written as a sum of coupling terms between nearest-neighbor sites. Nevertheless, with a view to subsequent analytical development, let us write the hamiltonian as a sum over irregular tetrahedral clusters, whose vertices lie on 4 different face-centered cubic sublattices [see Figs dra sharing a given site, but it is sufficient to choose only a subset of 4, in order to cover all the nearest-neighbor pairs [ Fig. 3(b)]. There are indeed 6 different possibilities to choose a proper subset, but the hamiltonian turns out to be independent of this choice. We thus finally write where the elementary contribution H ijkl will be denoted as tetrahedron hamiltonian. The subscripts i α , i β , i γ , i δ denote site configurations for the 4 vertices α, β, γ, δ of each tetrahedron. Possible site configurations are as follows: i = 0 denotes a vacancy (empty site), whereas i = 1, 2 denote a molecule in its 2 possible orientations (see Fig. 2). Due to the presence of only nearest-neighbor interactions, and assuming that (i, j), (j, k), (k, l), and (l, i) just refer to nearest-neighbor configurations, the tetrahedron hamiltonian can be written as where In the latter equation, n i is an occupation variable, defined as n i = 0 if i = 0 (empty site), n i = 1 other-wise (occupied site), whereas h ij is a bond variable, defined as h ij = 1 if the pair configuration (i, j) represents a H bond, and h ij = 0 otherwise. Let us also assume that i, j, k, l (in this order) denote configurations of sites placed on, say, A, B, C, D sublattices respectively. If A, B, C, D are defined as in Fig. 1, we can define h ij = 1 if i = 1 and j = 2, and h ij = 0 otherwise. The 1/4 coefficient is meant to avoid multiple countings of the chemical potential terms.
Let us now define the tetrahedron configuration probability by p ijkl , with the same convention about the subscript order, and assume that the probability distribution is equal for every tetrahedron. Such an assumption seems to be reasonable, as it allows one to characterize the two ordered structures present in the ground state. In both cases, the distribution is a delta function concentrated on a fixed tetrahedron configuration (i, j, k, l) (i.e., p ijkl = 1 for that configuration, and 0 otherwise), so that there is no residual entropy. As far as the low-density (LD) structure is concerned, it indeed turns out to be fourfold degenerate, since it can be represented by four alternative tetrahedron configurations, namely, (i, j, k, l) = (1, 2, 0, 0), (0, 1, 2, 0), (0, 0, 1, 2), (2, 0, 0, 1) (obtained from one another by a circular permutation). In each structure, a pair of sublattices (respectively, AB, BC, CD, and DA) is occupied by fully bonded molecules (forming a diamond structure), while the other two sublattices are empty. Conversely, the high-density (HD) structure turns out to be two-fold degenerate, being represented by the two alternative configurations (i, j, k, l) = (1, 2, 1, 2), (2, 1, 2, 1). Both these configurations have all lattice sites occupied, but they differ in the pairs of fully bonded sublattices, that is, AB and CD in the former case and BC and DA in the latter.
The grand-canonical energy per site in the thermodynamic limit can be written as which takes into account that there is one tetrahedron per site. By grand-canonical energy, we mean where u is the internal energy per site and ρ is the density, i.e., the average occupation probability We are now in a position to evaluate the energy of the two ground-state structures described above. For the LD phase, we have to replace into Eq. (4), obtaining while Eq. (6) correctly yields the density ρ LD = 1/2. Moreover, for the HD phase, we consider which, replaced into Eqs. (4) and (6), yields respectively and ρ HD = 1. Finally, we have also to consider a gas (G) phase, in which all but a finite number of sites are empty, so that both density and energy are 0 in the thermodynamic limit. It is easy to show that the G phase is thermodynamically favored (0 < w LD and 0 < w HD ) for the LD phase is favored (w LD < 0 and w LD < w HD ) for and the HD phase is favored (w HD < 0 and w HD < w LD ) for µ > µ LD-HD . Let us note that the LD phase has always a stability region, as µ G-LD < µ LD-HD , because of the repulsive term ǫ > 0.

III. CLUSTER-VARIATION APPROXIMATION
We have performed the finite temperature analysis by means of a cluster-variational approximation, previously applied by two of the authors for investigating a similar water-like lattice model [21,22]. The clustervariation method is an improved mean-field theory, which in principle can take into account correlations at arbitrarily large, though finite, distances. In Kikuchi's original work [34], an approximate entropy expression was obtained by heuristic arguments, while, in more recent and rigorous formulations [35], the approximation is shown to be equivalent to a truncation of a cluster cumulant expansion of the entropy. The approximation is expected to work, because of a rapid decreasing of the cumulant magnitude, upon increasing the cluster size, namely when the latter becomes larger than the correlation length of the system [36]. A particular approximation is defined by the largest clusters left in the truncated expansion, usually denoted as basic clusters. One obtains a free energy functional in the cluster probability distributions, to be minimized, according to the variational principle of statistical mechanics.
For our model we have chosen, as basic clusters, the irregular tetrahedra which we have previously used to express the model hamiltonian Eq. (1). Indeed, this choice turns out to coincide with the (generalized) firstorder approximation (on the tetrahedron cluster), which is also equivalent to the exact calculation for a Husimi lattice [37], whose (tetrahedral) building blocks are just arranged as in Fig. 3(b). The main advantage of this approximation is its simplicity, as the only clusters retained in the expansion are basic clusters and single sites (it is sometimes denoted as cluster-site approximation [38]), together with a substantial improvement over the ordinary mean-field theory, which has been recognized for different models, even with orientation dependent interactions [9]. In particular, let us note that the internal energy is treated exactly, since the range of interactions is smaller than the basic cluster size. The grand canonical free energy per site ω = w − T s (T being the temperature, expressed in energy units, and s the entropy per site, in natural units) can be written as where p X i is the probability of the i configuration for a site on the X sublattice (X = A, B, C, D). These probabilities can be obtained as marginals of the tetrahedron probability distribution as Accordingly, the variational free energy in Eq. (13) turns out to be a function of only the tetrahedron distribution p ijkl . Let us note that the free energy expression contains two different logarithmic terms, respectively corresponding to the cluster and site entropies. The former takes into account correlations among four configuration variables on a tetrahedron, while the latter can be viewed as a correction term ensuring that, if the tetrahedron distribution factorizes into a product of single site probabilities, the mean field (Bragg-Williams) entropy approximation is recovered. By the way, such a heuristic argument was first used by Guggenheim [39] to derive an equivalent expression for the case of a two-site cluster, also equivalent to the well-known Bethe approximation [40]. Minimization of ω with respect to the set of tetrahedron probabilities {p ijkl }, with the normalization constraint 2 i,j,k,l=0 can be performed by the Lagrange multiplier method, yielding the equation where ξ is related to the Lagrange multiplier, and can be obtained imposing the constraint Eq. (15) Eq. (16) is in a fixed point form, and can be solved numerically by simple iteration (natural iteration method [41]). For the cluster-site approximation, the numerical procedure can be proved to reduce the free energy at each iteration [37,41], and therefore to converge to local minima. From the solution of Eq. (16) one obtains a tetrahedron distribution p ijkl , whence one can compute the thermal average of every observable; density by Eq. (6), internal energy by Eqs. (4) and (6), and free energy by Eq. (13). The latter can be also related to the normalization constant as [37] ω = −T ln ξ.
Finally, assuming the volume per site equal to unity, pressure can be determined, in energy units, as P = −ω.
In the presence of multiple solutions, i.e., of competing phases, the thermodynamically stable one is selected by the lowest free energy (highest pressure) value. At zero temperature, we have P = −w, so that we can easily determine the transition pressures as follows. Replacing Eq. (11) into (8), we obtain the G-LD transition pressure while, replacing Eq. (12) into (10), we obtain the LD-HD transition pressure

IV. RESULTS AND DISCUSSION
In this section we present and discuss our results, including a detailed comparison with the Monte Carlo simulations of Girardi and coworkers [23]. With this purpose, we consider the same ratio of interaction parameters as in Ref. 23, namely, η/ǫ = 2. Let us note that the ground-state transition pressures computed above are independent of this parameter.

A. Phase diagram
In Figs. 4 and 5 we report two projections of the phase diagram, respectively in the temperature-pressure and density-temperature planes. We can observe three different phases, which we denote as G, LD, and HD, as they can be respectively identified as continuous evolutions of the three ground-state phases discussed above. The G phase is homogenous, in that the single-site probability distribution is independent of the sublattice, whereas the LD and HD phases exhibit a symmetry-breaking, as different sublattices behave differently. In the homogeneous phase, molecules can locally form H bonds, but a long-range ordered H-bond network does not appear. Conversely, the other two phases exhibit long-range order, and one can identify the same types of symmetrybreaking observed in the zero-temperature LD and HD phases, respectively. More precisely, upon decreasing temperature, the tetrahedron distributions tend to be concentrated on the fixed configurations, corresponding respectively to the different (LD or HD) H-bond networks described above. For instance, in the LD phase, H bonds are preferentially formed through AB (or BC, CD, DA) sublattices, whereas, in the HD phase, through AB and CD (or BC and DA). Let us have a closer look at the phase diagram topol-ogy. In the pressure vs temperature diagram (Fig. 4), we observe three different first-order transition lines. The first one separates the G and LD phases, whereas the other two occur between the LD and HD phases. All these lines are mapped onto coexistence regions in the temperature vs density diagram (Fig. 5). Both the LD-HD first-order transition lines terminates in tricritical points (denoted as TP and TP'), which turn out to be connected by a second-order transition line (i.e., a line of critical points), which encloses the LD phase. Another critical line separates the G phase from the HD phase, which correctly prevents any continuous transformation between these two phases, the latter having a lower symmetry. This critical line terminates in a critical end-point (CEP). Let us note, by the way, that, in the LD-HD coexistence region between CEP and TP, the LD phase has in fact a higher density than HD, so that the two phases are identified on the basis of their different symmetries. Let us also note that the LD phase exhibits a density anomaly, namely, a temperature of maximum density (TMD), at constant pressure.
Figs. 4 and 5 also display some data point obtained by the Monte Carlo analysis of Ref. [23]. We find a remarkably good agreement both for first-order transitions (coexistence regions in Fig. 5 and transition lines in Fig. 4), and for the TMD locus in the LD phase. Deviations are slightly larger only in the vicinity of critical lines, where the correlation length of the system diverges, so that we expect a break-down of our approximation performances. In spite of such an impressive agreement between Monte Carlo results and our approximate calculation, which suggests that the tetrahedron approximation captures the most relevant correlations present in the system, an important difference emerges between the two studies. Indeed, Girardi and coworkers [23] recognize the two long-range ordered network structures only in the ground-state analysis. Such a long-range order seems to be lost at finite temperature, so that they denote the two denser phases as low density liquid (LDL) and high density liquid (HDL). In this scenario, they only observe two first-order transitions lines (G-LDL and LDL-HDL), terminating in two different critical points.
Although in principle such a discrepancy might be an artifact of our approximate method, different evidences lead us to conjecture that this is not the case, and that in Ref. 23 the authors have simply not explored the possibility of a symmetry-breaking. A first evidence is of course the striking quantitative agreement shown above. Moreover, as previously mentioned, our symmetry-broken phases at T > 0 are continuous evolutions of those predicted at T = 0, whereas in Ref. 23 it is not clear how the transitions from the zero-temperature ordered phases to the finite temperature liquid phases occur. Looking for further evidences, we have analyzed isotherms and isobars.

B. Isotherms and isobars
In Fig. 6 we report isotherms in a pressure vs density diagram, together with corresponding Monte Carlo data, available from Ref. 23. In particular, we consider three different temperatures, showing qualitatively different behaviors. The first one is T /ǫ = 0.8 (Fig. 6a), where the isotherm displays two large plateaus, resulting from the two different (G-LD and LD-HD) first-order transitions discussed above. These plateaus can be clearly observed even from the Monte Carlo data, which are very well fitted by our curve. The second temperature is T /ǫ = 1.2 (Fig. 6b), lying between the two tricritical points TP and TP' (see Figs. 4 and 5). In this case, upon increasing pressure, the system first undergoes a G-LD first-order transition, and then a LD-HD secondorder transition. The former results in a plateau as well, whereas the latter gives rise to a simple "kink" in the isotherm. The Monte Carlo data exhibit a similar plateau and, noticeably, they also seem to be compatible with the possibility of a critical-line crossing, i.e., with the existence of a kink in the thermodynamic limit. Finally, we consider T /ǫ = 1.6 ( Fig. 6b), lying above both tricritical points. Here the system directly undergoes a (secondorder) G-HD transition, so that the isotherm displays just a slight kink. Even in this case, the slope change observed in the Monte Carlo isotherm seems to be well explained by the critical-line crossing.
Let us note that, looking at Figs. 4 and 5, we indeed expect two more qualitatively different regimes, respectively corresponding to the narrow temperature intervals between CEP and TP and between TP and the maximum temperature reached along the TP-TP' critical line. In the former case, upon increasing pressure, we find a G-HD second-order transition, followed by a HD-LD firstorder transition, and finally by a LD-HD second-order transition. In the latter case, we find the same sequence of transitions, all second-order. Corresponding isotherms are not displayed in Fig. 6. Indeed, we suppose that a comparison with Monte Carlo simulations may be very difficult in these temperature region, where maximum quantitative discrepancies are observed in the phase diagram (see Fig. 5). To be more precise, it is likely that similar behaviors might also appear from Monte Carlo data, but at a slightly lower temperature than observed from the approximate calculation.
In Fig. 7 we report isobars in a density vs temperature diagram (more precisely, Fig. 7a displays only results from our calculation, whereas in Fig. 7b some Monte Carlo results are superimposed). We have considered different pressure values, below the TP' pressure, pointing out four different regimes. The lowest pressure isobars lie below the CEP pressure P CEP /ǫ ≈ 1.4, so that the system only undergoes a first-order G-LD transition, upon decreasing temperature (see Fig. 4). At the lowest pressures, the density maximum turns out to be extremely shallow, but it becomes more pronounced upon increasing pressure (the TMD locus correctly connects all den-sity maxima). At pressure values lying between P CEP and the TP pressure P TP /ǫ ≈ 1.8, the system phase behavior becomes more complicated. Upon decreasing temperature, we first observe a second-order G-HD transition followed by a weak first-order HD-LD transition. The HD-LD transition also becomes continuous above P TP . Finally at P/ǫ > 2 [i.e., above the ground-state LD-HD transition pressure, see Eq. (20)], a further LD-HD firstorder transition appears at low temperature. A comparison of these results with Monte Carlo isobars drawn from Ref. 23 generally confirms the good agreement obtained for phase diagrams and isotherms. In particular, let us note that the G-HD critical line provides a clear interpretation of the slope changes observed in the two highest pressure isobars in Fig. 7b.
Looking at Fig. 4, we expect two more regimes, for pressure values above TP'. More precisely, between TP' and the maximum pressure value reached along the TP-TP' critical line, we observe a sequence of three continuous transitions, namely, G-HD, HD-LD, and LD-HD, upon decreasing temperature. For higher pressure, we finally observe a direct G-HD continuous transition. Even in these cases, corresponding isobars have not been reported in Fig. 7.

V. A MONTE CARLO TEST
As pointed out in the previous section, different evidences support the existence of long-range ordered symmetry-broken phases, and related critical transitions. Nevertheless, due to the approximate nature of our results, we have found it necessary to perform some specific checks by Monte Carlo simulations. We have only investigated certain regions around which the transitions are expected to appear.
In order to verify the G-HD transition, we consider a constant temperature T /ǫ = 1.7 (see Figs. 5 and 4), and work out a grand-canonical simulation for varying chemical potential. The lattice is made up of L × L × L bcc conventional cells (see Fig. 1), with L ranging from 10 up to 30, and periodic boundary conditions. The simulation is of the Metropolis type, starting from a random lattice configuration. We randomly select a lattice site, and flip its configuration according to the Metropolis rule. We define the Monte Carlo step as 2L 3 iterations of this process. Measurements are performed every 30-100 Monte Carlo steps, depending on the phase (the Monte Carlo dynamics is faster in the disordered G phase than in the LD and HD phases). We collect up to 10 5 measurements, waiting a thermalization time of order 10 3 measurements.
In Fig. 8 we report the fourth-order cumulant [42] of the grand-canonical energy V w as a function of the chemical potential, for different values of the lattice size L. We can observe that V w fluctuates around zero, except in a narrow interval (which shrinks upon increasing L) where V w exhibits a minimum followed by a maximum. This is a typical signature of a second-order transition [43], which in our case we expect to be the G-HD transition. Actually, in both phases far from the transition, the energy distribution can be described by a single peak, which can be approximated by a gaussian, tending to a Dirac delta in the L → ∞ limit. In this case V w should vanish independently of L. Upon approaching the critical line, the symmetric gaussian peak description is no longer valid, so that V w = 0 [43]. Though we do not report the details here, we have worked out an equivalent analysis for the LD-HD critical line, obtaining similar results.
The results displayed in Fig. 8 seem to indicate a transition around 1.3 µ/ǫ 1.4, where the simulation predicts a density 0.59 ρ 0.61. The latter values turn out to be slightly higher than the transition density predicted by the cluster-variation approximation ρ ≈ 0.56, at the corresponding transition chemical potential µ/ǫ ≈ 1.00. We then expect that the actual critical line is shifted towards higher density (i.e., higher pressure) with respect to the approximate line shown in Figs. 4 and 5. This result is in agreement with the known behavior of cluster-variational approximations, which, due to their mean-field nature, usually overestimate the ordered phase region.
In order to better characterize the transition, we define the following order parameter where p XY 1 is the probability of the i = 1 configuration (see Fig. 2) on the X and Y sublattices. A totally equivalent order parameter φ 2 could be obtained considering the i = 2 configuration. Let us note that, as previously mentioned, the HD phase is two-fold degenerate. Therefore, we expect that our simulated system can go into one of the two possible states, depending on initial conditions. In one case, preferential H bonding occurs either through AB and CD sublattices, with preferred i = 1 configurations on A and C. In the other case, preferential bonding occurs through BC and DA, with preferred i = 1 configurations on B and D. In the two cases we respectively obtain positive and negative values of the order parameter.
In Fig. 9 we report the absolute value of φ 1 , computed both by simulations (for different lattice sizes L) and by the cluster-variation approximation (in the thermodynamic limit), as a function of the chemical potential. In the latter case, we compute the two-sublattice proba-bilities as where the one-sublattice probabilities are derived from Eqs. (14). It is easy to see that φ 1 behaves like a good order parameter for the G-HD transition. Indeed, for low µ values, where we expect to observe the disordered G phase, φ 1 tends to zero upon increasing L. Conversely, for higher µ values, where we expect to observe the ordered HD phase, φ 1 remains nonzero, with a smaller effect of L. In the intermediate region, the slope of the curves increases upon increasing L, revealing the phase transition. Let us note, by the way, that for L = 10 (the size used in Ref. 23) finite size effects are definitely not negligible in the vicinity of the transition, where the correlation length of the system is expected to diverge. Accordingly, also the cluster-variational result fails in this region, although it progressively improves upon moving away from the transition.

VI. CONCLUSIONS
In this paper, we have considered a lattice model with orientation-dependent interactions, meant to describe a so-called network-forming fluid. As mentioned in the Introduction, this model is an instance of quite a large class of similar models (defined on the body-centered cubic lattice, with tetrahedral model molecules), originally conceived to describe water and investigate its various anomalies. Several studies have been carried out using a variety of approximate techniques, whose reliability has not always been checked. Motivated by a recent work by Girardi and coworkers [23], performing extensive Monte Carlo simulations for one of such models, we have worked out a corresponding analysis by means of a cluster-variational technique, with the aim of comparing the results and evaluating the reliability of the approximate method. The latter, based on a tetrahedral cluster, had been already employed by two authors of the present paper for investigating a different water-like lattice model [21], and actually turns out to be a generalization of the original first-order approximation worked out by Bell [10]. The generalization consists in taking into account the possibility of a symmetry-breaking at the level of four face-centered cubic sublattices of the original body-centered cubic lattice. Such a symmetry breaking is indeed known to occur in the model groundstate.
As already pointed out in the text, we obtain a remarkably good agreement with Monte Carlo results for most thermodynamic quantities. Nevertheless, our results clearly suggest that the the two denser phases, identified as liquid in Ref. 23, are in fact endowed with a long range ordered structure, closely resembling the hydrogen bond networks observed at zero temperature. Accordingly, the emerging picture of the phase diagram turns out to be much more complex and richer. In fact, due to symmetry-breaking, the two critical points observed in Ref. 23 turns out to be tricritical, and two new critical lines emerge, separating the two symmetry-broken phases from each other and from the disordered phase. Such critical lines are hardly detectable by observing macroscopically averaged quantities (such as the density), since they manifest themselves only as kinks (not plateaus) in isotherms and isobars, which are further smoothed out by finite size effects. We have indeed performed a Monte Carlo simulation, in order to check the qualitative correctness of our results, around the critical lines determined by the approximate method. The density exhibits in fact a smooth behavior, but the analysis of a Binder cumulant [43], and of a suitable order parameter, shows clear evidence of an order-disorder critical transition.
As previously mentioned, these results suggest that the tetrahedral cluster approximation is able to capture the most relevant correlations present in the system. On the other hand, unfortunately, the observed phase diagram seems to make the model less relevant for investigating water. Concerning this issue, the main problem is of course that the two ordered phases, which indeed may well represent two different ice forms [11], undergo second-order transitions, which have never been observed experimentally. It is likely that further variations of the model might improve the similarity of its phase diagram to that of real water, but a discussion of this point is beyond the scope of the present paper and is left for future work.