Intra-molecular coupling as a mechanism for a liquid-liquid phase transition

We study a model for water with a tunable intra-molecular interaction $J_\sigma$, using mean field theory and off-lattice Monte Carlo simulations. For all $J_\sigma\geq 0$, the model displays a temperature of maximum density.For a finite intra-molecular interaction $J_\sigma>0$,our calculations support the presence of a liquid-liquid phase transition with a possible liquid-liquid critical point for water, likely pre-empted by inevitable freezing. For J=0 the liquid-liquid critical point disappears at T=0.


I. INTRODUCTION
Water has an anomalous density decrease for isobaric cooling below a temperature of maximum density (TMD) [1]. Other thermodynamic anomalies-such as the rapid increase of the response functions-can be fit by power laws with apparent singularity well below the freezing temperature [1]. Several interpretations for this behavior have been proposed, but it is unclear which describes water or if any describes other anomalous liquids including, among the others, S, Se, Te, Cs, Si, Ge, I, C, P, SiO 2 , BeF 2 [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18].
One of the interpretations, the stability-limit conjecture [19], assumes that the limits of stability of the superheated, supercooled and stretched liquid form a single retracing spinodal line in the pressure-temperature (P -T ) plane. This scenario predicts a divergence of the response functions at the supercooled liquid-to-liquid spinodal [20].
The singularity-free interpretation [21,22] envisages the experimental data represent apparent singularities, due to anticorrelated fluctuations of volume and entropy. In this scenario, these fluctuations are responsible for the TMD line.
The liquid-liquid phase transition hypothesis [23] proposes the presence of a first order line of phase transitions separating two liquid phases differing in density, the high density liquid (HDL) and the low density liquid (LDL). In this scenario the HDL-LDL phase transition, possibly ending in a liquid-liquid critical point, is responsible for the anomalies.
Although Refs. [8,9], by tuning parameters of the corresponding models, predict smooth transitions from the different scenarios, to help elucidate which is the most reasonable description for water, we consider a model fluid with inter-molecular and intra-molecular interactions. This model, in the particular case of a zero intra-molecular interaction, recovers the model introduced by Sastry et al. [22] that predicts the singularity-free scenario. Our aim is to understand how the presence of an intra-molecular interaction changes this prediction.
We perform analytic calculations in a mean field approximation and an off-lattice Monte Carlo (MC) simulation. Our results show that a non-zero intra-molecular interaction gives rise to a HDL-LDL phase transition, with a possible critical point [23], with the liquid-liquid critical temperature decreasing to zero and vanishing with the intra-molecular interaction. Therefore, at least for this model, the singularity-free scenario is obtained only in the particular case of a zero intra-molecular interaction, while for a finite intra-molecular interaction the HDL-LDL phase transition is predicted. General considerations suggest that the liquid-liquid phase transition for water could occur below the glass temperature, i.e. outside the accessible experimental range.
The paper is organized as follows: in Sec. II we define the model defined on a lattice. In Sec. III we describe the equation of state approach in the mean field approximation and we present the mean field results. In Sec. IV we introduce the off-lattice model, we describe the MC approach and show the simulation results. In Sec. V we discuss our results and we give the conclusions.

II. THE LATTICE MODEL
The fluid is represented by partitioning the system into N cells of equal size. A variable n i is associated with each cell i = 1, . . . , N , with n i = 1 if the cell is occupied by a molecule, n i = 0 otherwise.
The inter-molecular interaction [22], has a first term describing the van der Waals attraction between molecules, where ǫ > 0 is the energy gain for two nearest neighbor (nn) occupied cells and the sum is over all the possible nn cells. The second term in Eq. (1) accounts for the dynamic network of hydrogen bonds (HBs) formed by liquid wa-ter, with each molecule typically bonded to four other molecules at low T [1], with an energy-gain J > 0 per HB. We consider cells with size of a water molecule and with four arms, one per possible HB. For the molecule in the cell i, the orientation of the arm facing the cell j is represented by a Potts variable σ ij = 1, . . . , q, with a finite number q of possible orientations. Two molecules in nn cells form a HB only if they are correctly oriented [1], i.e., by assumption [22], if δ σij ,σji = 1 (δ a,b = 1 if a = b and δ a,b = 0 otherwise).
The experimental oxygen-oxygen correlation function shows that a HB is formed if and only if the intermolecular distance is within a characteristic range [1]. Hence, we assume that the formation of a HB leads to a local volume expansion [22], where V 0 is the volume of the liquid with no HBs, is the total number of HBs in the system, and v HB is the specific volume per HB. Experiments show that the relative orientations of the arms of a water molecule are correlated, with the average H-O-H angle equal to 104.45 o in an isolated molecule, 104.474 o in the gas and 106 o in the high-T liquid [24], suggesting an intra-molecular interaction between the arms. This interaction must be finite, because the angle changes with T , consistent with ab-initio calculations [25] and molecular dynamics simulations [26]. Hence, we introduce the intra-molecular term [27,28] where for each of the 4 C 2 = 6 different pairs (k, l) i of the arms of a molecule i, with the appropriate orientation (δ σ ik ,σ il = 1), there is an energy gain J σ > 0. For J σ = 0 we recover the model of Ref. [22] that predicts the singularity-free scenario, and where the HBs are uncorrelated, inhibiting the orientational long-range order. We study the general case with finite J σ , by using (i) a mean field approximation and (ii) MC simulations.

III. THE EQUATION OF STATE APPROACH
The equation of state of our system is implicitly given by where is the total internal energy and µ is the chemical potential. From Eqs. (1)-(4), we rewrite the equation of state as Here is the effective attraction energy, depending on P and the local arm configuration, is the effective HB interaction energy due to the HB volume-increase, and is the effective local chemical potential depending on the local arm configuration.

A. The mean field approximation
The mean field approximation consists in assuming a linear relation between the number density of liquid cells, and the density order parameter m ∈ [−1, 1], and between the number density of arms in the appropriate state for a HB, and the orientational order parameter m σ ∈ [0, 1], i.e.
Hence, the molar density ρ ≡ nN/V is Here v 0 ≡ V 0 /N , and p HB ≡ n 2 p σ is the probability of forming a HB between two nn molecules, where n 2 is the probability of finding two nn molecules and is the probability of having the facing arms of the two molecules in the appropriate orientational state for a HB.
For T → ∞ we expect m σ → 0, hence p σ → 1/q. For T → 0 the finite values of ǫ, J and J σ allow us to assume a cooperative effect and an orientational long-range order in a preferred state, with m σ → 1, hence p σ → 1.

B. The cooperative effect
To include the cooperative effect, we consider that each arm σ ij interacts with a mean field h generated by the other three arms on the same molecule, in addition to the effective interaction J ′ (P ) with the facing arm σ ji on a nn molecule. Since the energy is minimized when the arms are in the same orientational state, the system breaks the symmetry ordering in the preferred Potts state. Hence, we choose h proportional to n σ , to J σ and to the number of arms generating h, i.e.
Our results do not depend on this choice for h, and are recovered using higher order approximations for h. We assume i. e. that p σ [Eq. (15)], for two nn molecules interacting with the surrounding, coincides with the probability δ σij ,σji h , for the facing arms (σ ij , σ ji ) of two isolated nn molecules, of being in the same orientational state under the action of the field h. By definition is where the right-most-side is the explicit calculation of the left-side with partition function Z h . Here the sum is over all the configurations of the two variables σ ij , σ ji , and the symbols are with the Boltzmann constant k B chosen as unitary hereafter. As expected for p σ , also δ σij ,σji h → 1/q for T → ∞ and δ σij ,σji h → 1 for T → 0. Numerically we find that the solution of the Eq. (17) is m * σ (T, P ) = 0 for P > P max (T ), and m * σ (T, P ) > 0 for P ≤ P max (T ), where m * σ (T, P ) = 0 corresponds to lack of orientational order and m * σ (T, P ) > 0 corresponds to the orientational long-range order. In this mean field approximation, P max (T ) turns out to be well described by a decreasing linear function of T .
C. The Gibbs free energy Next, we write a mean field expression for the molar Gibbs free energy as a function of the two order parameters m and m * σ , where is the molar energy as derived by Eqs. (1)-(4) and Eq. (6), with the mean field approximations i,j v ≡ 1/ρ is the molar volume derived by Eq. (14); is the molar entropy, are the standard mean field expressions for the entropy of N variables n i and for the entropy of 4nN q-states Potts variables for the arms, respectively, and n * σ is the number density n σ of HBonded arms [Eq. (13)] evaluated in m * σ (T, P ).

D. The mean field results
By numerically minimizing g(T, P ) with respect to m and m * σ with the constraint that m * σ is solution of Eq. (17), we find the equilibrium values of m(T, P ) and m * σ (T, P ). By using Eq. (14), we find of ρ(T, P ) at equilibrium ( Fig. 1).
At high P the mean field theory predicts that ρ(T ) increases when T decreases (Fig. 1). At low P , for decreasing T , the theory predicts (i) a discontinuity in ρ(T ), corresponding to the liquid-gas first order phase transition ending in the liquid-gas critical point C (Fig. 2), (ii) the TMD decreasing for increasing P , (iii) a discontinuity in ρ(T ) at low T , disappearing at lower P .
The first two predictions are consistent with either the singularity-free scenario or the liquid-liquid phase transition hypothesis, while the third prediction is consistent only with the HDL-LDL first order phase transition P v0/ǫ 0.25, by decreasing T , ρ has a discontinuity at high T , from a low value (in the gas phase) to a high value (in the liquid phase), then ρ has a maximum followed by a smooth saturation to the finite value ρHB = 0.5/v0 corresponding to the full-HBonded liquid. For 0.3 P v0/ǫ 0.4, by decreasing T , there is no discontinuity in ρ, but there is a maximum in ρ and the saturation to ρHB. For 0.5 P v0/ǫ 1.2, by decreasing T , ρ has a maximum and then a discontinuity to ρHB. For 1.25 P v0/ǫ 1.6, ρ has only a maximum, and, for higher P , ρ regularly increases, by decreasing T . (b) Blowup of the low-T region. Both discontinuities reported show a first order phase transition, each ending in a critical point (Fig. 2). hypothesized in the latter scenario. In particular, the smooth disappearing of the discontinuity at lower P is consistent with a phase transition line, with a negative slope in the P -T phase diagram, ending in a HDL-LDL critical point C ′ (Fig. 2).

IV. THE OFF-LATTICE MODEL
To show that our mean field predictions are robust, we now use a completely different approach based on an offlattice (OL) model representing a system with a homogeneous distribution of molecules in the available volume V , that we divide in N equivalent cells of volume V /N .
As a consequence of the homogeneity of the system, for each cell the N degrees of occupancy freedom n i are set to n i = 1. In analogy with Eq. (2), the total volume V is defined as Here N HB and v HB are defined as in Eqs.  Also, following the lattice model, the molecules have four arms described by four q-state variables σ ij , with the HB interaction defined by the second term in Eq. (1) and the intra-molecular interaction by Eq. (4). These interactions are both independent of the distance among first-neighbor molecules and depend only on the arms orientation σ ij .
In this off-lattice model, the average distance between two molecules r is a continuous variable, so we replace the first term in Eq. (1) with where R 0 ≡ √ v 0 is the hard-core diameter of each molecule. In analogy with the Eq. (1), we consider this off-lattice van der Waals energy independent of the HB expansion, so in U W (r) we use

A. The Monte Carlo simulation
We perform MC simulations, in two dimensions [29], at constant N , P , T and variable V (N P T ensemble) with N ∈ [10 2 , 10 4 ]. The MC dynamics consists in updating the variables σ ij and the variable V OL 0 , accepting the new state with probability exp[−(∆U W +P ∆V /k B T )] if ∆U + P ∆V > 0, or with probability 1 if ∆U + P ∆V < 0. Here ∆U ≡ ∆(U W + H IM ) and ∆V are the changes of total internal energy and total volume Eq.(25), respectively, after the update. Our results for the average density ρ MC ≡ N/V , averaged on 6×10 5 MC steps after 1.2×10 5 MC steps of thermalization at each T , are qualitatively consistent with the mean field prediction (Fig. 3).
By MC simulations, we find (i) for J σ = 0, no liquidliquid phase transition and the TMD line, consistent with mean field in Ref. [22]; (ii) for J σ ≫ J, for any P in the liquid phase, at T below the TMD line, a discontinuity in density, suggesting a first order phase transition along a line with a negative slope in the P -T phase diagram, consistent with mean field for J σ → ∞ [27,28]; (iii) for 0 < J σ < J, a phase transition line ending in a critical point C ′ (Fig. 4), and occurring at increasing P and decreasing T for decreasing J σ .
To verify that the jumps found in the MC density are marking a first order phase transition, instead of a narrow continuous phase transition, we study the isothermal compressibility,

Its maximum
increases linearly with the number of particles N at a first order phase transition [30], while at a second order phase transition K max T is proportional to both N and the fluctuation of the ordering parameter, scaling as a power of N [31]. Therefore, the finite size scaling analysis on K max T (N ) allows us to discriminate between a continuous and a discontinuous phase transition.
To estimate the maximum with a great precision, we use a continuous T algorithm, the histogram reweighting method [32]. By checking the minimum T and the maximum P where the behavior of K max T (N ) fails to be linear, we estimate the critical point at T C ′ /ǫ = 0.045 ± 0.005 and P C ′ v 0 /ǫ = 0.841 ± 0.042 for N → ∞ (inset Fig.4).
Next, we obtain the coexisting lines by extrapolating to N → ∞ the values P (T, N ) corresponding to K max T (N ) at fixed T , both for T ≤ T C and T ≤ T C ′ . Furthermore, our results are consistent with the necessary condition that the K max T (T ) line emanates from the critical point [22,33], both for T > T C ′ and for T > T C (Fig.4).

V. DISCUSSION AND CONCLUSIONS
The mean field and the MC results for our water-like fluid model for a finite intra-molecular interaction J σ show qualitatively the same phase diagram. In both approaches the TMD line decreases with increasing P , consistent with the experiments [4]. In particular, both approaches predict a first order phase transition in the liquid phase, occurring at low T and at a pressure P max (T ) decreasing for increasing T . For J σ < J we find that the first order phase transition line ends in a critical point, separating, by necessity, two phases with the same symmetry, in this case two liquids, HDL and LDL, being the critical point in the liquid phase. Three considerations are in order here.
The first consideration is related to the comparison with the result of Ref. [22], here recovered for J σ = 0. The Sastry et al. model [22], upon HB formation, accounts for (i) inter-molecular orientational correlation [Eq.(1)]; (ii) local expansion with lowering temperature [Eq.(2)]; (iii) anticorrelation between V and S, because the formation of HBs decreases the number of possible orientational configurations for the system, hence the entropy S decreases for increasing N HB , i.e. for increasing V . This is expected in a system with a density anomaly, because (∂V /∂T ) P < 0 implies (∂S/∂V ) T < 0. Finally, the Sastry et al. model assumes the arms of a molecule completely independent (J σ = 0), and predicts the singularity-free scenario.
We tested, by preliminary MC calculations, that for J σ → 0 the HDL-LDL critical point C ′ moves to lower T and to P max (T = 0). In particular, the phase transition disappears at T = 0 for J σ = 0, while the effect of the decreasing J σ on the location of the TMD line is weak. Although these results require a longer analysis, beyond the scope of the present work, they show that the predictions of Ref. [22] are recovered in the limit J σ = 0, confirming the validity of our MC approach [34]. We, therefore, conclude that in this model the presence of a finite J σ is responsible for the appearance of the first order phase transition, with a possible HDL-LDL critical point C ′ .
The second consideration is related to the possibility that this low-T HDL-LDL phase transition is pre-empted by inevitable freezing in real water. Recent analysis of the realistic model for water ST2 [15] suggests that the HDL-LDL critical point may occur above the glass temperature T g , though as yet still outside the easily accessible experimental range.
We compare our MC results with data for real water. From the location of the liquid-gas critical point (T C , P C ) and the TMD line at ambient pressure (T * , P * ) (Table I), we find the ratios P * /P C and T * /T C in real water. By assuming that the same P * /P C holds in our MC case, we calculate the corresponding P * ∼ 2.76 10 −4 ǫ/v 0 in our model and then we estimate the T * corresponding to P * from the TMD line in the MC phase diagram (Fig. 4). In this way we find a ratio T * /T C from the MC results that is consistent with the real water data, suggesting the validity of our assumption on P * /P C . Therefore, we use the same kind of assumption also to estimate the glass temperature T g for our phase diagram. In particular, from real water data we obtain the ratio T g /T C at ambient pressure and, assuming that it holds also for our model, we estimate T g /ǫ ∼ 0.42 at P v 0 /ǫ ∼ 2.76 10 −4 for the MC phase diagram. Hence, for our model with the parameters chosen in this paper, is T g > T C ′ , i.e. the HDL-LDL critical temperature at P C ′ ≃ 0.841ǫ/v 0 is below the glass temperature at P * ∼ 2.76 10 −4 ǫ/v 0 . From the study of the phase diagram of real water [4] is reasonable to expect that T g (P ) decreases for increasing P , therefore our analysis does not exclude that T C ′ is above T g (P C ′ ). However, by considering a very large J σ , such that the HDL-LDL critical pressure is ∼ 2.76 10 −4 ǫ/v 0 , we can compare T C ′ and T g at the same pressure P * . Our preliminary results show that the HDL-LDL critical temperature is in this case very close to T g .
As a consequence of this analysis, our model supports the possibility that the HDL-LDL critical point is located deep into the supercooled region, below or close to the glass temperature, depending on the value of J σ . Therefore, the liquid-liquid phase transition could be preempted, if our model is representative of the thermodynamic properties of the real system. This results is analogous to what has been proposed for silica [15], another liquid with density anomaly, suggesting that the present model could provide a general theoretical framework for anomalous molecular liquids.
The third consideration is about the role of the tetrahedrality in determining the properties of anomalous liquids. For the present model we do not consider a tetrahedral geometry in the two-dimensional MC approach, and the geometry is not explicitly defined in the mean field approach. Nevertheless, our results are consistent with the experimentally accessible phase diagram of real water, suggesting that the tetrahedral network is not an essential feature for the anomalous behavior of water-like liquids.
This conclusion is consistent with what has been observed by Angell in Ref. [7] and is well described by a general cooperative model [6] with a generic drive to phase separate the excitations into distinct regions of space (clustering). In our model the drive is given by the intra-molecular interaction, that mimics the geometrical drive in tetrahedral liquids even if is not necessarily limited to the tetrahedral case.
In conclusion, we studied the effect of an intramolecular interaction J σ in a model for anomalous molecular liquids with a mean field approach, valid for J σ > 0, and with an off-lattice MC simulation. For J σ > 0 we found a HDL-LDL phase transition while our MC results confirm that for J σ = 0 the singularity-free scenario holds [22]. Hence, the two interpretations originate from the same mechanism with a different hypothesis on the intra-molecular interaction; the latter is strictly valid I: Characteristic temperatures and pressures for real water and for the present model, for the gas-liquid critical point (TC, PC ), the TMD at ambient pressure (T * , P * ) and the glass temperature Tg at ambient pressure. The ratios P * /PC and Tg/TC are not available for the model. We assume that the corresponding H2O values are valid also for the model, and we use these values to estimate P * and Tg, respectively, for the model. Temperatures are measured in K for H2O and in ǫ for the model. Pressures are measured in MPa for H2O and in ǫ/v0 for the model. only for J σ = 0. Within this framework, the most reasonable scenario for water includes a HDL-LDL phase transition, probably hindered by inevitable freezing. Our results suggest also that the tetrahedrality is not essential to understand the anomalous behavior in water-like liquids.