Hydrogen bond correlated percolation in a supercooled water monolayer as a hallmark of the critical region

TL;DR

This study models supercooled water monolayers, revealing that hydrogen bond percolation correlates with the liquid-liquid phase transition and critical point, providing insights into the thermodynamics of supercooled water.

Contribution

It introduces a theoretical model mapping supercooled water monolayers onto a percolation problem, linking hydrogen bond clusters to phase transition phenomena.

Findings

Percolation line coincides with the liquid-liquid phase transition.

Percolation corresponds to maxima of thermodynamic response functions.

A positive slope percolation transition is observed at higher pressures.

Abstract

Numerical simulations for a number of water models have supported the possibility of a metastable liquid-liquid critical point (LLCP) in the deep super-cooled region. Here we consider a theoretical model for a supercooled liquid water monolayer and its mathematical mapping onto a percolation problem. The mapping allows us to identify the finite-size clusters at any state-point, and the infinite cluster at the critical point, with the regions of correlated hydrogen bonds (HBs). We show that the percolation line coincides with the first-order liquid-liquid phase transition ending at the LLCP. At pressures below the LLCP, the percolation line corresponds to the strong maxima of the thermodynamic response functions and to the locus of maximum correlation length (Widom line). At higher pressures, we find a percolation transition with a positive slope and we discuss its possible relation with the thermodynamics.

Equations

$$V_{\mathrm{tot}}\equiv V+N_{\mathrm{HB}}{v}_{\mathrm{HB}},$$

$$\mathcal{H}\equiv \sum _{ij}U(r_{ij})-J{N}_{\mathrm{HB}}-J_{\sigma }{N}_{\mathrm{coop}}.$$

$$U(r)\equiv \{\begin{array}{l}\infty \text{ for }r<r_0\\ \\ 4ϵ\left[{\left(\frac{r_0}{r}\right)}^{12}-{\left(\frac{r_0}{r}\right)}^6\right]\text{ for }{r}_0\le r\le 25r_0\\ \\ 0\text{ for }r>25r_0\end{array},$$

$$N_{\mathrm{HB}}\equiv \sum _{⟨ij⟩}{n}_i{n}_j{\delta }_{{\sigma }_{ij},{\sigma }_{ji}},$$

$$N_{\mathrm{coop}}\equiv \sum _i{n}_i\sum _{(l,k{)}_i}{\delta }_{{\sigma }_{ik},{\sigma }_{il}},$$

$$H\equiv U(v)-J_{\mathrm{eff}}\sum _{⟨i,j⟩}{n}_i{n}_j{\delta }_{{\sigma }_{ij},{\sigma }_{ji}}-J_{\sigma }\sum _i{n}_i\sum _{(k,l{)}_i}{\delta }_{{\sigma }_{ik},{\sigma }_{il}}+Pv,$$

$$Z(T,P)\equiv \sum _{\{\sigma \}\{v\}}{e}^{-H/k_{B}T},$$

$$Z(T,P)=\sum _{\{v\}}{e}^{-[U(v)+Pv]/k_{B}T}\times Z_{\{\sigma \}}$$

$$\begin{array}[]{llll}Z_{\{\sigma\}}&\equiv&\sum_{\{\sigma\}}&e^{(J_{\rm eff}/k_{B}T)\sum_{\langle i,j\rangle}n_{i}n_{j}\delta_{\sigma_{ij},\sigma_{ji}}}~{}\times~{}e^{(J_{\sigma}/k_{B}T)\sum_{i}n_{i}\sum_{(k,l)_{i}}\delta_{\sigma_{ik},\sigma_{il}}}\\ \\ &=&\sum_{\{\sigma\}}&\prod_{\langle i,j\rangle}\left[1+\left(e^{(J_{\rm eff}/k_{B}T)}-1\right)n_{i}n_{j}\delta_{\sigma_{ij},\sigma_{ji}}\right]\\ \\ &{}&{}&\times\prod_{i=1}^{N}\prod_{(k,l)_{i}}\left[1+\left(e^{(J_{\sigma}/k_{B}T)}-1\right)n_{i}\delta_{\sigma_{ik},\sigma_{il}}\right]\end{array}$$

$$\begin{array}[]{llll}Z_{\{\sigma\}}&=&\sum_{\{\sigma\}}&\prod_{\langle i,j\rangle}\left[1+\left(e^{(J_{\rm eff}/k_{B}T)}-1\right)\delta_{\sigma_{ij},\sigma_{ji}}\right]\\ \\ &&&\times\prod_{i=1}^{N}\prod_{(k,l)_{i}}\left[1+\left(e^{(J_{\sigma}/k_{B}T)}-1\right)\delta_{\sigma_{ik},\sigma_{il}}\right].\end{array}$$

$$p_{J_{\sigma }}\equiv 1-e^{-J_{\sigma }/k_{B}T}$$

$$|C|+|B|=\mathcal{N}-|D|$$

$$\sum _C{p}_{J_{\sigma }}^{|C|}(1-p_{J_{\sigma }}{)}^{|B|}=1,$$

$$\begin{array}[]{lcl}\sum_{C}\!\!\!\!\!\!\!\!\!\!\raisebox{-9.0pt}{${\scriptstyle fixed\ |D|}$}~{}~{}p_{J_{\sigma}}^{|C|}(1-p_{J_{\sigma}})^{|B|}&\equiv&\sum_{|C|=0}^{{\cal N}-|D|}\left(\raisebox{8.0pt}{${\scriptstyle|C|+|B|}$}\!\!\!\!\!\!\!\!\!\!\raisebox{-8.0pt}{${\scriptstyle|C|}$}~{}~{}~{}\right)p_{J_{\sigma}}^{|C|}(1-p_{J_{\sigma}})^{|B|}=\\ &=&[p_{J_{\sigma}}+(1-p_{J_{\sigma}})]^{{\cal N}-|D|}\end{array}$$

$$\begin{array}[]{cl}Z_{\{\sigma\},J_{\sigma}}&\equiv\sum_{\{\sigma\}}\prod_{i=1}^{N}\prod_{(k,l)_{i}}\left[1+\left(e^{(J_{\sigma}/k_{B}T)}-1\right)\delta_{\sigma_{ik},\sigma_{il}}\right]\sum_{C}p_{J_{\sigma}}^{|C|}(1-p_{J_{\sigma}})^{|B|}\\ &=\sum_{\{\sigma\}}\sum_{C}W(\{\sigma\},C)\end{array}$$

$$W(\{\sigma \},C)=\prod _{i=1}^N\prod _{(k,l{)}_i}\left[1+\left(e^{(J_{\sigma }/k_{B}T)}-1\right){\delta }_{{\sigma }_{ik},{\sigma }_{il}}\right]p_{J_{\sigma }}^{|C|}(1-p_{J_{\sigma }}{)}^{|B|}$$

$$W(\{\sigma \},C)=0$$

$$Z_{\{\sigma \},J_{\sigma }}=\sum _C\sum _{\{\sigma \}}W(\{\sigma \},C),$$

$$1+\left(e^{(J_{\sigma}/k_{B}T)}-1\right)\delta_{\sigma_{ik},\sigma_{il}}=\left\{\begin{array}[]{ll}1&\mbox{if }\sigma_{ik}=\sigma_{il}\\ 1-p_{J_{\sigma}}&\mbox{otherwise}\end{array}\right.$$

$$W(\{\sigma \},C)=p_{J_{\sigma }}^{|C|}(1-p_{J_{\sigma }}{)}^{|A|}{\delta }_{\{\sigma \},C}$$

$${\delta }_{\{\sigma \},C}=\prod _{⟨k,l⟩\in C}{\delta }_{{\sigma }_{ik},{\sigma }_{il}}$$

$$Z_{\{\sigma \},J_{\sigma }}=\sum _C{p}_{J_{\sigma }}^{|C|}(1-p_{J_{\sigma }}{)}^{\mathcal{N}-|C|}{q}^{N(C)}.$$

$$p_J\equiv 1-e^{-J/k_{B}T}$$

$$Z_{\{\sigma \}}=\sum _C{p}_J^{|C_J|}(1-p_J{)}^{{\mathcal{N}}_J-|C_J|}\times p_{J_{\sigma }}^{|C|}(1-p_{J_{\sigma }}{)}^{\mathcal{N}-|C|}{q}^{N(C)}$$

$${\mathcal{P}}_{\infty }\equiv 1-{\mathcal{P}}_{<\infty }.$$

$$h_3(h_1,h_2)\approx |h_2{|}^{2-\alpha }{f}^{\pm }\left(\frac{h_1}{|h_2{|}^{\gamma +\beta }}\right)$$

$${\varphi }_1\equiv {\left(\frac{\partial h_3}{\partial h_1}\right)}_{h_2},{\varphi }_2\equiv {\left(\frac{\partial h_3}{\partial h_2}\right)}_{h_1}.$$

$${\chi }_1\equiv {\left(\frac{\partial {\varphi }_1}{\partial h_1}\right)}_{h_2},{\chi }_2\equiv {\left(\frac{\partial {\varphi }_2}{\partial h_2}\right)}_{h_1},$$

$${\chi }_{12}\equiv {\left(\frac{\partial {\varphi }_1}{\partial h_2}\right)}_{h_1}={\left(\frac{\partial {\varphi }_2}{\partial h_1}\right)}_{h_2}$$

$$h_1=h_2=h_3=0$$

$$h_1=0.$$

Full text

††journal: Journal of Molecular Liquids

1 Introduction

A fundamental feature of water molecules is their capability to form hydrogen bonds (HBs) [1]. In the liquid phase the HBs organize in a quasi-tetrahedral network that is continuously built up and broken on timescales of picoseconds [2]. Although there is scientific consensus about the relevance of the HBs for the anomalous properties of water, the consequences of the HB peculiar properties on the supercooled water phase diagram are under debate since decades [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. In bulk water near freezing each molecule minimize the free energy by favoring a hydration shell made of four tetrahedral H-bonded molecules. For decreasing temperature $T$, this mechanism leads to increasing fluctuations in volume $\langle V^{2}\rangle$ and in energy $\langle E^{2}\rangle$–associated to the isothermal compressibility $K_{T}$ and the isobaric specific heat $C_{P}$, respectively– and to increasing cross-fluctuations in volume and energy, related to the isobaric thermal expansion coefficient $\alpha_{P}$. The thermodynamics response functions $K_{T}$ and $C_{P}$ increase below 46 and 35 oC, respectively, while $\alpha_{P}$ becomes negative below 4 oC and increases in absolute value [31, 32, 33, 34]. However, their largest increase occur in the supercooled region, below the melting line, where water can be liquid down to temperatures between -46 and -42 oC, approximately [35, 36, 37]. Similar results have been observed at negative pressure and moderate supercooling [38].

The rapid increase of the fluctuations has suggested the possible existence of a divergency at approximately -45 oC [33], below or near the experimental limit at which the present technology allows us to probe the metastable liquid [35, 36, 37]. This possibility has been supported by the seminal work of Poole et al. in 1992, in which the existence of a liquid-liquid critical point (LLCP) in the deep supercooled regime was suggested based on numerical simulations of the ST2 model [6], one of the many atomistic models for water. The LLCP hypothesis since then has been corroborated by a series of theoretical and experimental results [39, 40, 13, 41, 42, 43, 44, 45, 46, 18, 19, 47, 48, 26], although different scenarios are still under debate [49, 50, 10, 13, 51]. One of them [51, 16] recently was invalidated due to conceptual errors at its origin [29, 30], while another [49] was found in a colloidal model [24], but does not seem to be compatible with recent experimental data for water [36, 38].

The LLCP scenario hypothesizes the existence of a first order phase transition line separating two metastable liquid phases, the high-density liquid (HDL) at higher pressures and temperatures, and the low-density liquid (LDL) at lower pressures and temperatures. The presence of a critical point (CP) would imply the existence of the Widom line [8, 9, 52, 53], that emanates from the CP and is the locus in the $T$–$P$ plane where the statistical correlation length $\xi$ has a maximum. The Widom line is calculated along the constant thermal field [52, 53] and not necessarily coincides with the maxima of $K_{T}$, $C_{P}$ or $\alpha_{P}$ [54], neither represents a unique thermodynamic separatrix of the supercritical region [55]. A detailed definition of the Widom line is presented in the Appendix. The $\xi$ maximum along the Widom line increases approaching the CP and diverges at the CP.

Nonetheless, no experiment so far have found a direct evidence of the LLCP in supercooled water. The reason is that the LLCP, if present in bulk water, should occur where the lifetime of the metastable liquid is too short for being probed with the present experimental techniques. Therefore, different strategies have been developed to perform experiments that could help us in establishing if the LLCP occurs in bulk water. Among them water solutions [56] or mixtures [57] have been used, as well as water under confinement [58].

However, confinement-dependent effects on the LLCP are difficult to predict without detailed studies. In the following we study the case of a supercooled single layer of water confined in two-dimensional geometry. To this goal, we investigate a coarse-grain water model by means of a percolation mapping that allow us to identify the critical region of the system.

Such a description is based on a mathematical mapping developed by Kasteleyn and Fortuin [59], and Coniglio and Klein [60]–namely the site-bond correlated percolation–between physical and geometric properties [59, 61, 60, 62, 63, 64, 65, 66, 67, 68] which guarantees that the correlation length $\xi$ statistically coincides with the connectivity length $\xi_{c}$, i.e. the average size of the cluster.

A few works studying the percolation region in the super-critical liquid water have been published so far [69, 70, 71, 72]. They are based on a definition of a cluster as a contiguous region of H-bonded molecules. As we will describe in the next sections, this definition provides an overestimation of the region of statistically correlated molecules and cannot be stricly related to the thermodynamic behavior of water. Indeed, the $T$-$P$ locus of the percolation transition and its nature change according to the definition of cluster of molecules, which is in principle arbitrary. For example, it is possible to define a cluster as the contiguous region of molecules whose distance is below a given threshold, or as the contiguous region of molecules bonded via HB, etc. Each definition is rational and identifies a specific feature of the system, but none of them is in principle related to the size $\xi$ over which the statistical fluctuations of the degrees of freedom spread. On the other hand, the Kasteleyn-Fortuin/Coniglio-Klein approach (described in the section 3.2), allows to define a probability $p$ that two degrees of freedom belong to the same cluster, in such a way that the critical thermodynamic behavior is preserved in the percolation description or, in other words, that the clusters have a characteristic length-scale (connectivity length) that coincides with $\xi$.

Here, hence, we characterize the formation of the HB network in relation with the clusters of correlated water molecules. According to the $T$ and $P$ conditions, we compute the mean cluster-size and cluster-size distribution, identifying the critical region where the clusters span all over the system and mark the onset of a percolation transition. We show how the percolation transition is associated to the building up of the ordered HB network, although with different mechanisms at higher and lower pressures. Finally, we discuss the low-$P$ percolation transition in relation to the Widom line in supercooled water.

2 The FS model for a water monolayer

We adopt the many-body model for a water monolayer introduced in 2002 by Franzese and Stanley (FS) [73, 74] and that has been shown to reproduce, at least qualitatively, many of the water properties [75, 13, 76, 77, 78, 79, 80, 81, 82, 14, 83, 44, 84, 85, 86]. When the many-body interaction is zero the FS model coincides with the model introduced in Ref.[50] by Sastry, Debenedetti, Sciortino and Stanley. The monolayer is kept between parallel walls at sub-nanometer distance, $h=0.5$ nm, with purely repulsive (excluded volume) hydrophobic interaction, allowing deep supercooling without crystallization, consistent with simulations and experiments [87, 88].

We consider the model in the Gibbs ensemble, keeping the number $N$ of molecules, the pressure $P$ and the temperature $T$ fixed. Therefore, the volume $V$ changes according to the equation of state and fluctuates when we perform our Monte Carlo simulations. Assuming that the density is homogenous, this implies that the (average) distance $r$ between the molecules changes with $P$ and $T$. We partition the system into cells whose volume coincides with the (average) volume occupied by each molecules. Hence, by construction, each cell has the size equal to the inter-molecule distance $r$ and the cell’s size changes with $P$ and $T$. Because we consider here a monolayer, we adopt a square partition, whose coordination number coincides with the number of HBs that each molecule forms at low temperature. In three dimensions other partitions are more appropriate [89].

To reduce the number of degrees of freedom of the system, in accordance with FS, we coarse grain the molecules translational coordinates as follows. By construction each cell includes always a water molecule and its volume $v$ is at least equal to the van der Waals volume $v_{0}$ of the molecule. Hence, $v_{0}/v$ is the cell density in van der Waals units. We associate an index $n_{i}=0$ if $v_{0}/v\leq 0.5$. For our choice of geometry and parameters ($v_{0}\equiv hr_{0}^{2}$ with $r_{0}\simeq 2.9\AA$ van der Waals diameter), it is $n_{i}=0$ when $r\geq\sqrt{2}~{}2.9\AA\geq 3.7\AA\equiv r_{\rm max}$ maximum O–O elongation of a straight HB as calculated by ab initio molecular dynamics simulations from the proton-transfer coordinate [90], assuming a covalent distance O–H $\simeq 1\AA$. Therefore, a cell with $n_{i}=0$ cannot form HBs, although it includes a molecule. On the other hand, we associate to the cell an index $n_{i}=1$ if $v_{0}/v>0.5$, i.e. $r\leq r_{\rm max}$, and the molecule in it can form HBs. Hence, the index $n_{i}$ is a discretized density field that marks if the molecule $i$ can, or cannot, form HBs. If we consider that all the HBs are formed in the liquid phase, then $n_{i}=1$ is associated to a liquid-like density, and $n_{i}=0$ to a gas-like density, recalling a lattice-gas model for argon-like systems. Hence, in the liquid phase, all the indices $n_{i}$ are equal to 1, while they are all 0 in the gas phase. Therefore, the average $\langle n_{i}\rangle=n_{i}$ does not play the role of an order parameter as in the lattice-gas model. Instead, here the order parameter for the liquid-gas transition is associated to the total density, that, as the total volume, is a continuous function of $P$ and $T$.

Furthermore, water unlike argon and other usual liquids, has HBs. Following FS, we assume that all the possible heterogeneities of water are due to the HBs. When water forms HBs within its hydration shell, its coordination number reduces to 4, in an almost perfectly tetrahedral configuration. In this fully-bonded configuration the volume per molecule is larger than in configurations with larger coordination number, because the volume occupied by the hydration shell as a whole is not changing in a sensible way, while the number of molecules in it does change [91]. This can be represented, approximately, by associating a proper volume $v_{\rm HB}$ to each new HB – the larger the number of HBs in the system, the better the approximation. We choose $v_{\rm HB}/v_{0}=0.5$, equal to the average volume increase between high-$\rho$ ices VI and VIII and low-$\rho$ (tetrahedral) ice Ih.

As a consequence, since the number $N$ of molecules is fixed, the total volume $V_{\rm tot}$ occupied by the system increases linearly with the number of HBs, $N_{\rm HB}$, i.e.,

[TABLE]

where $V\equiv Nv$ is the volume without HBs. Despite $V$ is homogeneously distributed among the $N$ water molecules, the number of HBs changes among molecules, giving rise to local heterogeneities in the density field.

The Hamiltonian of the FS model is given by

[TABLE]

The first term, summed over all the water molecules $i$ and $j$ at O–O distance $r_{ij}$, is given by a truncated Lennard-Jones (LJ) model

[TABLE]

where $\epsilon\equiv 5.8$ kJ/mol, close to the estimate based on isoelectronic molecules at optimal separation $\simeq 5.5$ kJ/mol [92]. The potential $U(r)$ takes into account the van der Walls (dispersive) attraction and hard core (electron) repulsion between water molecules.

As already noted, the formation of HBs, does not affect the volume occupied by the full-bonded hydration shell, but only the number of molecules included in it. Therefore, the distance $r$ between the molecules is not modified by the HBs [91]. Hence, the van der Waals interaction is not affected by the HBs 111We acknowledge discussion with the late Prof. David Chandler for noticing that this point is crucial to state that the FS model is not mean field.. Because in the FS model $r$ is a continuous variable, it is appropriate to represent the van der Waals interaction with a LJ model. We truncate the LJ potential at large distances, as usually done in continuous models for numerical efficiency, and also at short distance. In particular, we replace the repulsive power 12 222While the power 6 of the LJ model can be derived from first principke calculations, the power 12 is arbitrary. with a hard-core at the van der Waals diameter. Our previous analysis show that both truncations do not affect the results and simplify the implementation of the model [93]. The truncated LJ potential drives the liquid-gas phase transition. It plays a fundamental role in the vicinity of the liquid-gas spinodal and is the only relevant interaction for temperatures above the spinodal temperature. Without this term the model would not reproduce correctly the fluid phases of water.

The second term in Eq.(2) represents the short-range, directional, covalent [94] component of the HB, with $J/4\epsilon=0.5$, i.e. $J\simeq 11$ kJ/mol, close to the estimate of this energy constant that can be derived from the optimal HB energy and a HB cluster analysis [13]. Here

[TABLE]

the number of HBs, is the sum over nearest neighbors (n.n.) pairs of water molecules $i$ and $j$ at a distance $r\leq r_{\it max}$ [95, 96] (i.e., with $n_{i}n_{j}=1$) and in the same bonding state ($\delta_{\sigma_{ij},\sigma_{ji}}=1$), where

$\sigma_{ij}=1,\dots,q$ is the bonding variable of molecule $i$ facing the n.n. molecule $j$, with $\delta_{ab}=1$ if $a=b$, 0 otherwise.

The bonding variables are introduced to account correctly for the variation of entropy and energy associated to the formation or breaking of a HB, as explained in the following. If two n.n. molecules $i$ and $j$ form a HB, the system energy decreases by $-J$ and the system entropy decreases by $-k_{B}\ln q$ ($k_{B}$ is the Boltzmann constant) because both molecules have bonding variables in the same state, $\sigma_{ij}=\sigma_{ji}$.

Each molecule has 4 bonding variables, one for each possible HB, and $q^{4}$ possible bonding configurations. The HB between two oxygens is broken if ${\widehat{\rm OOH}}>30^{\circ}$, or ${\widehat{\rm OOH}}<-30^{\circ}$, as estimated from Debye-Waller factors [97, 96]. Therefore, only $1/6$ of the entire range of possible orientations $[0,360^{\circ}]$ in the OH—O plane is associated to a bonded state, and each HB formation leads to an entropy decrease equal to $-k_{B}\ln 6$. Therefore, by choosing $q=6$ we guarantee the correct HB definition and entropy loss.

The third term in Eq.(2) accounts for the HB cooperativity due to O–O–O correlation, that in bulk leads the molecules toward an ordered tetrahedral configuration [91]. Such an effect originates from quantum many-body interactions of the HB [98, 99]. The number of cooperative pairs of bonding variables in the system is

[TABLE]

where, for each molecule $i$, $(l,k)_{i}$ indicates each of the six different pairs of the four indexes $\sigma_{ij}$ of the molecule.

Therefore, the enthalpy of the system can be written as

[TABLE]

where $J_{\rm eff}\equiv J-Pv_{HB}$ is the effective interaction between $\sigma$-variables of n.n. molecules that depends on $P$ and $U(v)\equiv\sum_{i,j}U(r_{ij})$, where $r_{ij}$ is a function of $r=\sqrt{Nv/h}$: $r_{ij}=r$ if $i$ and $j$ are n.n., $r_{ij}=\sqrt{2}r$ if $i$ and $j$ are next n.n., etc.

In the $NPT$ ensemble the partition function of the system is

[TABLE]

where the sum is over all the possible configurations of bonding variables $\{\sigma\}$ and cell volumes $\{v\}$. The Eq. (7) can be rewritten as

[TABLE]

where

[TABLE]

where $\prod_{\langle i,j\rangle}$ runs over all the n.n. molecules $j$ of the molecule $i$, $\prod_{i=1}^{N}$ runs over all molecules and $\prod_{(k,l)_{i}}$ extends over all the six pairs of bonding variables of a specific molecule $i$.

In the following we set $J_{\sigma}/4\epsilon=0.05$, i.e. $J_{\sigma}\simeq 1$ kJ/mol, in such a way that, for each molecule whose bonding configuration changes from three bonding variables in the same state to four (adding a cooperative HB), the energy decreases by $-3J_{\sigma}\simeq-3$ kJ/mol. This value is consistent with the energy decrease

of -3 kJ/mol [100] of each HB observed in ice Ih with respect to liquid water, in the reasonable hypothesis that this energy change is entirely attributed to the HB cooperativity [101, 13]. Our choice of the model parameters, with $J_{\sigma}\ll J$, guarantees an asymmetry between the two components of the HB interaction, such that bonding variables can organize cooperatively only if each single HB is formed. From here on we express $T$ and $P$ in internal units, $4\epsilon/k_{B}$ and $4\epsilon/v_{0}$ respectively.

3 Site-bond correlated percolation

3.1 Percolation mapping

The percolation theory studies the connectivity properties of clusters defined on a $d$-dimensional lattice. We talk of site or bond percolation depending on whether the clusters are made of n.n. vertices or links of the lattice, respectively, that have been randomly occupied with a probability $p$.

For each lattice size $L$, the smallest value of $p$ at which there is 50% probability of finding at least a cluster spanning the entire lattice is indicated as $p^{*}$. The thermodynamic limit $p^{*}(L\rightarrow\infty)\equiv p_{c}$ is called the percolation threshold.

The theory shows that the percolation on large lattices displays a second-order phase transition where $p$ is the control field and the order parameter is the probability $\mathscr{P}_{\infty}$ for an arbitrary lattice site to belong to the percolating cluster [102, 103]. The percolation transition is characterized by critical exponents for $p\rightarrow p_{c}$, in analogy with thermodynamic second-order phase transitions where $p$ is replaced by the thermal field [104, 102, 103].

In order to recover the behavior of a fluid near the liquid-gas (LG) critical point, Fisher listed [105] the properties that the “right” clusters (droplet) must have in order to represent thermodynamically correlated regions [102, 106]: i) an infinite droplets is formed only at the critical point; ii) the order parameter of the phase transition is related to the size of the infinite droplet; iii) the compressibility (or susceptibility for magnetic systems) of the order parameter is proportional to the mean cluster-size; iv) the thermodynamic correlation length $\xi$, quantifying the spatial extent of thermal fluctuations, is proportional to the average radius of the clusters, i.e. the connectivity length $\xi_{c}$.

The “right” cluster definition for Ising-like systems was proposed in the ’70s by Kasteleyn and Fortuin (KF). They shown that it is possible to map a ferromagnetic Potts model onto a corresponding percolating model [59, 61]. Later, Coniglio and Klein (CK) [60, 62] introduced the “random site-correlated bond” percolation (site-bond correlated percolation, for short) to prove

the equivalence between the thermodynamic and percolation critical behavior. The CK approach is at the base of the “Cluster Monte Carlo” dynamics introduced by Swendsen and Wang (SW) [107], representing a major improvement in our ability to generate equilibrium configurations near a critical point and at lower temperatures [108].

3.2 Site-bond correlated percolation for the FS model

Following the KF and CK approaches, we map the physical system onto a site-bond correlated percolation model where bonds are set between n.n. variables $\sigma$ with a probability that depends on their state, their interaction, $P$ and $T$, defining clusters of thermodynamically correlated degrees of freedom. Because in the FS model the $\sigma$-variables can interact with two coupling constants—$J$, between n.n. molecules, and $J_{\sigma}$, within the same molecule—we define two probabilities functions—$p_{J}$ and $p_{J_{\sigma}}$, respectively—to set bonds between, and within, molecules.

In order to find the functional forms for such probabilities in the condensed phase, we observe that the $\sigma$-part of the partition function, Eq. (9,) in the liquid phase, where is $n_{i}=1$ $\forall i$, reduces to

[TABLE]

We now introduce the cluster definition and show how the partition function can be equally expressed as a sum over all possible clusters associated to a given configuration $\sigma$ [109]. Following CK, we define a cluster $C$ as a region of variables connected through fictitious bonds. The bonds do not affect the interaction energy but are introduced to define the connectivity between the variables [110].

In our coarse-grain model the degrees of freedom that are relevant for the HB formation are the variables $\sigma$. We therefore, perform the percolation mapping on the part of the partition function that depends on the bonding variables $\sigma$, $Z_{\{\sigma\}}$.

For sake of simplicity, let’s assume that we have only $J_{\sigma}$ interactions. Then we will generalize to the case including also $J$ interactions. We put fictitious bonds on $J_{\sigma}$ interactions with probability

[TABLE]

between two $\sigma$-variables in the same state and define a cluster as maximal set of $\sigma$-variables connected by bonds. We call these bonds $b_{J_{\sigma}}$. We will show now that the definition in Eq.(11) allows us to map the thermodynamic model into the percolation model.

First we observe that for a given configuration $\{\sigma\}$ we have several bond configurations that are “compatible” with $\{\sigma\}$, i.e. such that the fictitious bonds are set only between $\sigma$-variables in the same state. Let’s call $C$ one of these bond configurations, $|C|$ the number of bonds in the configuration $C$, $|B|$ the number of “missing” bonds, i.e. the number of couples of $\sigma$-variables in the same state without a bond, $|D|$ the number of couples of n.n. $\sigma$-variables in different states (Fig. 1).

If we fix the configuration $\{\sigma\}$, then $|D|$ is fixed and

[TABLE]

where ${\cal N}$ is the total number of n.n. pairs of $\sigma$-variables in the system. Therefore, for a given configuration $\{\sigma\}$ is

[TABLE]

because

[TABLE]

for the Newton binomial. Therefore, we can write

[TABLE]

where

[TABLE]

is the weight associated to the $\sigma$-variables configurations $\{\sigma\}$ together with the compatible bond configuration $C$.

If we set

[TABLE]

when $\{\sigma\}$ and $C$ are not compatible, then we can write

[TABLE]

where now the $\sum_{C}$ is done over all the bond configurations, independent of their compatibility with $\{\sigma\}$.

For n.n. $\sigma$-variables

[TABLE]

therefore,

[TABLE]

where $|A|\equiv|B|+|D|$ and

[TABLE]

takes into account the constraint on the compatibility of $\{\sigma\}$ and $C$.

If $N(C)$ is the number of clusters in the bond configuration $C$, the number of configurations $\{\sigma\}$ compatible with $C$ is $q^{N(C)}$, since every cluster can have only $q$ states with all the $\sigma$-variables equals. Therefore, summing over $\{\sigma\}$ in Eq.(16) and taking into account Eq.(17) and Eq.(12), we have

[TABLE]

In the case in which we consider also $J$ interactions, we set fictitious bonds, called $b_{J}$, between n.n. $\sigma$-variables interacting with $J$ coupling with probability

[TABLE]

and, generalizing the previous discussion and the notation in a straightforward way (Fig. 1), we can rewrite Eq.(10) as

[TABLE]

where the cluster configuration $C$ includes bonds of both type $J$ and $J_{\sigma}$. This expression of the partition function depends only on percolation quantities, completing the mapping between the thermodynamic system and the percolation problem.

Therefore, the mapping into the site-bond correlated percolation does not change the partition function of the model and its free energy. Furthermore, following FK/CK [59, 61, 60, 62], it can be shown that this definition satisfies the Fisher’s conditions for the “right” clusters, giving a geometric representation of the thermodynamic droplets. As already mentioned, this observation allows us to adopt the well know SW Cluster Monte Carlo dynamics [107, 111] that is able to generate equilibrium configurations in a very efficient way near phase transitions, where other Monte Carlo dynamics slow down and freeze [108].

4 Results

4.1 The thermodynamic behavior

Based on the percolation description of the system it is possible to implement a cluster Monte Carlo algorithms [111, 108] allowing faster equilibration in the deep supercooled region. In each Monte Carlo (MC) step we generate a cluster of fictitious bonds following the prescription of KF-CK-SW, with all the variables $\sigma_{ij}$ of the cluster in the same state by construction. Next, we change at random the $\sigma$ state of the entire cluster, whose average linear dimension $\xi_{c}$ is related to the correlation length $\xi$. Hence, approaching to a critical point, where $\xi$ becomes as large as the system size, the update involves a number of HBs of the order of the entire system, allowing a fast decorrelation of the new configurations and a more efficient MC sampling of independent configurations at equilibrium. This efficiency strongly contrasts with the critical slowing down that the local dynamics, such as Metropolis or Heat-Bath [108], suffer near a critical point due to the increase of the correlation length $\xi$.

In the following we consider systems with increasing sizes, from $N=2.5\times 10^{3}$ up to $N=160\times 10^{3}$ water molecules, with periodic boundary conditions. We perform annealing simulations at constant $P$ adopting the Cluster MC algorithm described above, following Ref.s [44, 79].

In particular, our previous finite-size analysis of the model [44] shown that the system is consistent with a thermodynamic LLCP for $P_{C}v_{0}/(4\epsilon)\simeq 0.555\pm 0.002$ and $T_{C}k_{B}/(4\epsilon)\simeq 0.0597\pm 0.0001$ with respect to a mixed-field order parameter $M=\rho^{*}+b^{\prime}u^{*}$ [112, 113, 114], given by a linear combination of number density $\rho^{*}\equiv\rho v_{0}$ and energy density $u^{*}\equiv E/(4\epsilon N)$, where $b^{\prime}$ is the mixing-coefficient defined in the Appendix, with $b^{\prime}=0.25\pm 0.03$ for the explored range of system sizes333This definition of the order parameter, as shown in Ref. [112, 113, 114], allows to recover the Isigin-like symmetry of the order parameter distribution, at the critical point, between the ordered and the disordered phases.. Such a critical point represents the ending point of a first order phase transition line with a negative slope in the $P$-$T$ plane, separating two metastable liquid phases, the HDL at higher $P$ and $T$, and the LDL at lower $P$ and $T$. Starting from the LLCP we found a line of maxima of $\xi$ that extends toward negative pressures until reaching the stability limit of the stretched liquid. Following the definition given in the Appendix we identified this line of $\xi$-maxima as the Widom line [44].

4.2 Percolation probability.

The size $s$ of a cluster is defined as the number of $\sigma$ variables belonging to it. In Fig. 2 we show an example of clusters configuration for a specific $\sigma$-variables configuration. The occurrence of a percolation transition is marked by the appearance of an infinite cluster spanning the entire system, i.e. a cluster with a linear size $\xi_{c}$ comparable or equal to the system’s size (Fig.3). The probability that an arbitrary $\sigma$-variable is part of a finite cluster of size $s$ is $\mathscr{P}_{s}\equiv n_{s}s$, where $4Nn_{s}$ is the average number of finite clusters of size $s$ , i.e. the fraction of $\sigma$-variables that belong to clusters of size $s$. Therefore, the probability that an arbitrary $\sigma$-variable belongs to any finite cluster is $\mathscr{P}_{<\infty}\equiv\sum_{s}\mathscr{P}_{s}$, and its probability to belong to the infinite cluster, i.e. the percolation probability for any $\sigma$-variable, is [102]

[TABLE]

The quantity $\mathscr{P}_{\infty}$ represents the order parameter for the percolation transition, with $\mathscr{P}_{\infty}=0$ when there is no percolation and $0<\mathscr{P}_{\infty}\leq 1$ otherwise (Fig. 4.a). For finite systems $\mathscr{P}_{\infty}$ increases rapidly near the percolation transition, with a larger slope for larger systems (finite size effect, Fig. 4.b,c,d.). By convention, for any finite size the percolation threshold is at the thermodynamic point where $\mathscr{P}_{\infty}=0.5$.

Our numerical findings reveal the occurrence of a percolation transition in all the range of simulated pressures. According to $\mathscr{P}_{\infty}$ we can distinguish at least three regions 444This separation in three regions is reinforced by the analysis of the cluster-size distribution presented elsewhere. . The first region corresponds to $Pv_{0}/(4\epsilon)\lesssim 0.5$, where the percolation transition is sharp, becoming sharper for larger $N$, with an evident finite-size effect for $N<25.6\times 10^{3}$ (Fig. 4.b). The second region is for $0.6<Pv_{0}/(4\epsilon)<1$, where $\mathscr{P}_{\infty}$ increases sharply at the percolation threshold and the finite-size effect is weaker than at lower $P$, being more evident only for $N<10\times 10^{3}$ (Fig.s 4.c). The third region corresponds to $Pv_{0}/(4\epsilon)>1$, where the percolation transition is smoother for any system size $N$ (Fig. 4.d).

Our analysis reveals two percolation lines in the thermodynamic plane ($T$,$P$) emanating from the state point (0,1), expressed in units $4\epsilon/k_{B}$ and $4\epsilon/v_{0}$, respectively (Fig. 5). The first line, for $Pv_{0}/(4\epsilon)<1$, has a negative slope that decreases for decreasing $P$ and increasing $T$, becoming almost $T$-independent for sufficiently low $P$. This percolation line coincides, according to the previous analysis [44], with the first-order LLPT for $0.55\lesssim Pv_{0}/(4\epsilon)<1$, while for $Pv_{0}/(4\epsilon)\lesssim 0.55$ it coincides with the locus of strong maxima of $C_{P}$, $K_{T}$ and thermal expansivity $\alpha_{P}$. Furthermore, in Ref. [44] we shown that the locus of strong maxima of the response functions is also where the thermodynamic correlation length $\xi$ of the HBs has a maximum. Hence, the percolation line at low $P$ corresponds to the Widom line.

The second percolation line is observed for $Pv_{0}/(4\epsilon)>1$ and has a positive slope that increases for increasing $P$ and $T$ (Fig. 5). It coincides approximately with the thermodynamic locus of weak maxima of $C_{P}$ along isotherms that we found at high $P$ [44]. As we have shown [44], this locus coincides also with that of weak maxima of $K_{T}$ along isobars and with that of weak minima of $\alpha_{P}$ along isotherms, consistent with similar lines found in other models with the LLCP [52, 115, 116, 117].

5 Discussion

5.1 Hydrogen bond network and the percolation line.

For our choice of parameters and $Pv_{0}/(4\epsilon)\lesssim 1$, i.e. for $P$ such that $J_{\sigma}\lesssim J_{\rm eff}$, the percolation line marks the threshold for the formation of a network of correlated HBs. In this range of pressures, the number of HBs, $N_{\rm HB}$, increases monotonically and saturates to its maximum value upon cooling [83]. The increase is smooth at low $P$ and becomes sharper approaching $Pv_{0}/(4\epsilon)=1$ (Fig. 6).

Mazza et al. shown that in this range of $P$ the number of cooperative HBs, i.e. those HBs that become correlated, sharply increases at $Tk_{B}/(4\epsilon)\simeq 0.05$, with a locus of maximum-derivative that coincides with (i) the LLPT at $P>P_{C}$, (ii) the LLCP at the critical pressure, and (iii) the locus of strong maxima of $K_{T}$ (indicated as $K_{T}^{\rm coop}$ in Ref. [83]) at $P<P_{C}$. Here, these loci coincide with the percolation threshold for $Pv_{0}/(4\epsilon)\leq 1$ (Fig.5), as we discuss next.

For $P<P_{C}$, the HB network is gradually formed and the percolation threshold occurs when the system is highly H-bonded, i.e. $N_{\rm HB}\simeq 95\%N_{\rm HB}^{\rm max}$ (Fig. 6).

Nevertheless, at $T$ above the percolation line the clusters of H-bonded molecules are made of isolated HBs, i.e. the clusters are made, mainly, of $b_{J}$ bonds on $J$ interactions. On the other hand, at $T$ below the percolation line these clusters coalesce into clusters that span the entire system as a consequence of the setting of a macroscopic number of $b_{J_{\sigma}}$ bonds on $J_{\sigma}$ interactions. This setting implies local rearrangements of the $\sigma$-variables toward more ordered configurations under the action of the $J_{\sigma}$ interactions, i.e. a local reordering of the HB network due to the many-body interaction among the water molecules.

The local reordering of the molecules does not affect the density of the system, but involves large fluctuations in energy and entropy. As a consequence, the percolation line coincides with the loci of strong maxima in isobaric specific heat, isothermal compressibility and isobaric heat expansion coefficient previously reported in Refs. [83, 44].

The comparison of the percolation line for $Pv_{0}/(4\epsilon)<1$ (Fig. 5) with the locus of maxima of $\xi$, reported in Fig.6.a,b of Ref. [44], reveals that both loci coincides within the numerical error. Hence, this confirms that our percolation mapping is capable to capture the statistical fluctuations of the water molecules

The mapping between the percolation problem and the thermodynamic system strictly holds only at a critical transition, because the cluster definition assumes a diverging correlation length at the percolation threshold. Therefore, from the percolation mapping we can draw useful conclusions for the thermodynamic system only along the critical isobar. Nevertheless, the coincidence of the percolation line with the LLPT, above $P_{C}$, and the loci of strong fluctuations, below $P_{C}$, shows that the clusters of correlated HBs play a role in a large range of pressures, even if their spanning is not associated to critical fluctuations of thermodynamic quantities.

This observation reconciles two classes of thermodynamic scenarios that have been debated as possible explanations of the anomalous water behavior. The first class assumes the occurrence of a LLPT and includes two scenarios: the LLCP-scenario [6] and the LLCP-free scenario [10]. The second class includes the Singularity-free scenario [50] that can be interpreted assuming a percolation behavior [118, 119]. Here, as in Refs. [118, 119], the correlation length $\xi$ increases as $T$ decreases as a consequence of the increasing size of clusters made of correlated water molecules. However, the correlated percolation approach presented here is different from the one proposed in the “polychromatic correlated-site percolation” [118, 119], in which tiny “patches”, made of water molecules with all HBs formed, are of lower density than the rest of the HB network and give rise to anomalous density fluctuations. Furthermore, here we give the bond probability as an explicit function of $T$ and $P$, while in Refs. [118, 119] the authors give only qualitative relation between the percolation quantities and the thermodynamic observables, leaving the determination of the functional form to the experiments.

5.2 The percolation at extremely high $P$

We find a percolation line, for our choice of parameters, also for $Pv_{0}/(4\epsilon)>1$ (Fig. 5). At these pressures the percolation line is characterized by a smoother change in the order parameter $\mathscr{P}_{\infty}$ (Fig. 4d). The HB formation is enthalpically unfavoured, because $Pv_{0}>J$ and $J_{\rm eff}<0$. Hence, a decrease of $T$ induces HBs breaking with a decrease of $N_{\rm HB}$ (Fig. 6).

Nevertheless, $J_{\sigma}>0$, hence the four $\sigma$-variable of the same water molecule assume the same state at low $T$, maximizing the probability to set bonds $b_{J_{\sigma}}$. At the same time, $J_{\rm eff}<0$ induces different $\sigma$-states on n.n. water molecules. Under these conditions the bonds $b_{J}$ have the maximum probability when the $\sigma$-variables have different states. Hence, at $Pv_{0}/(4\epsilon)>1$ the clusters are made of anticorrelated water molecules. It can be shown555Evidences will be provided elsewhere. that under these conditions the percolation belongs to the class of the the random-bond percolation [102] and that its percolation temperature increases for $Pv_{0}/(4\epsilon)>1$ as long as $|J_{\rm eff}|<J_{\sigma}$, i.e. up to $Pv_{0}/(4\epsilon)\simeq 1.1$, consistent with our finding (Fig.5).

Interestingly, this random-bond percolation line could be identified with the “Kertész line” or “Coniglio-Klein line” [120, 121] extending from a percolation critical point at $T=0$ and $Pv_{0}/(4\epsilon)=1$ up to $P=\infty$. It is still under debate if the Kertész line corresponds to a specific thermodynamic locus or not. For example, it has been proposed that this line exists in the supercritical phase of a simple fluid [122, 123] and is related to the locus where vanishes the surface tension of droplets made of the denser fluid.

Finally, it is worth noting that the possibility of the presence of a second critical point in liquid supercooled water, at higher $P$ and lower $T$ with respect the LLCP, has been predicted by Strekalova et al. [80] in the case of water confined in a disordered matrix of hydrophobic nano-particles. Such a second CP could coincide with the critical percolation point at $T=0$ and $Pv_{0}/(4\epsilon)=1$.

6 Conclusions

In this work we present a percolation study for a supercooled water monolayer adopting a many-body water model that has the LLCP [73, 74, 44]. The geometrical description of the system is based on the KF-CK-SW mathematical mapping [59, 61, 62, 124, 107] of the correlated site-bond percolation that allows us to compute clusters of correlated degrees of freedom along the critical isobar of the LLCP. In such a way, at $P_{C}$ (i) the average size of a cluster coincides with the thermodynamically correlated regions, (ii) the cluster connectivity length $\xi_{C}$ coincides with the thermodynamical correlation length $\xi$, and (iii) the percolation set of critical exponents belongs to the same universality class of the thermodynamic critical exponents.

The percolation description can be extended to any $P$, although the mapping with the thermodynamically correlated regions at $P\neq P_{C}$ is not exact. The discrepancy comes from the fact that the percolation phenomena is critical, by construction, at any $P$, while the thermodynamic system not necessarily is. Hence at any $P$ we always reach a percolation threshold with a diverging $\xi_{C}$, even where $\xi$ is not critical. Nevertheless, the percolation approach reveals interesting results not only at $P_{C}$.

In particular, we find a line of percolation transitions, with negative slope in the $P$–$T$ plane, that starts at $Pv_{0}/(4\epsilon)=1$ and $T=0$ and extends up to the liquid-gas spinodal at negative pressures and $T>0$. This line coincides with the loci of strong extrema of the specific heat $C_{P}$, compressibility $K_{T}$ and thermal expansivity $\alpha_{P}$, as resulting from previous thermodynamic analysis [83, 44]. For pressures above $P_{C}v_{0}/(4\epsilon)\simeq 0.6$ and below $Pv_{0}/(4\epsilon)=1$, the percolation transition coincides with the thermodynamic LLPT separating two liquids with different energy and density [13, 83, 44].

Near $P_{C}$ the critical behavior of the site-bond correlated clusters is consistent with our thermodynamic analysis [44] and the appearance of the spanning cluster corresponds to a cooperative reorganization of the network of correlated HBs. At lower pressures, $Pv_{0}/(4\epsilon)\lesssim 0.5$, we observe the coincidence of the percolation line with the Widom line estimated from the direct calculation of isobaric maxima of $\xi$ [44]. Hence, at low pressures the cluster size increases more than the correlated regions, because $\xi_{C}$ diverges while $\xi$ only has a maximum marking the Widom line.

Furthermore, a modified version of the model for the water monolayer, developed by Vilanova and Franzese and including crystal phases and water polymorfism [125], reveals the occurrence of a hexatic phase separated from the liquid by a critical line in the deep supercooled region. This finding has been confirmed by molecular dynamic simulations of an atomistic model for a bilayer of water [126]. It is, therefore, interesting to analyze if the percolation transition could be eventually related to the hexatic critical line at low $P$ in the extended model. We are currently investigating this possibility.

Our percolation approach leads, instead, to a different description at pressures above the LLPT. Emanating from the state point at $T=0$ and $Pv_{0}/(4\epsilon)=1$, we find a percolation line with positive slope, extending up to infinite pressure, in the $P$–$T$ plane. This line is characterized by a smooth percolation transition as in the random-bond percolation and approximately corresponds to the loci of weak extrema found for $K_{T}$, $C_{P}$ and $\alpha_{P}$ at these pressures in this model [44], consistent with similar lines found in other models with the LLCP [52, 115, 116, 117].

An intriguing hypothesis is that this percolation line at $Pv_{0}/(4\epsilon)>1$ could be identified with the Kertész line [120], emanating from the percolation point at $T=0$ and $Pv_{0}/(4\epsilon)=1$. Similar lines have been observed in the Ising model around the Curie point and in Lennard Jones systems around the LG critical point [106]. Such line, that in our model marks the decrease of the number of HBs, $N_{\rm HB}$, at very high pressure, could be related to the locus where vanishes the surface tension of droplets made of the denser fluid [122, 123]. Further investigation, beyond the scope of this work, is necessary to understand this possibility.

Acknowledgments

V. B. acknowledges the support from the Austrian Science Fund (FWF) grant No. M 2150-N36, and from the the European Commission through the Marie Skłodowska-Curie Fellowship No. 748170 ProFrost. G. F. acknowledges support from the ICREA Academia prize, funded by the Generalitat de Catalunya and from Spanish MINECO grant FIS2015-66879-C2-2-P.

Appendix

According to the scaling theory, and following Refs. [53], [52], the critical behavior of a fluid can be described in terms of two independent scaling fields: the ordering (strong) field $h_{1}$ and the thermal (weak) field $h_{2}$ and the Widom line is identified as the locus in the $T$-$P$ plane where $\xi$ has a maximum, calculated/measured along the path $h_{2}=$constant.

This can be seen by considering that a generalized expression for the free energy is given by the scaling field $h_{3}$, depending on $h_{1}$ and $h_{2}$, that close to the critical point is

[TABLE]

where $\alpha$, $\beta$ and $\gamma$, related to each other by the scaling relation $\alpha+2\beta+\gamma=2$, are the critical exponents characteristic of the universality class to which the fluid belongs. The superscript $\pm$ of the analytical scaling function $f^{\pm}$ refer to $h_{2}>0$ and $h_{2}<0$, respectively.

The critical phase transition is described by two scaling densities, associated to the fields $h_{1}$ and $h_{2}$: the ordering density (order parameter) $\phi_{1}$ and the thermal density $\phi_{2}$

[TABLE]

The variation of the scaling densities with respect to the scaling fields define the “susceptibilities” of the fluid

[TABLE]

respectively known as strong susceptibility, weak susceptibility and cross susceptibility. For a liquid system, the thermodynamic response functions experimentally accessible are i) the isothermal compressibility $K_{T}\equiv-(1/\langle V\rangle)(\partial\langle V\rangle/\partial P)_{T}$, ii) the isobaric coefficient of thermal expansion $\alpha_{P}\equiv(1/\langle V\rangle)(\partial\langle V\rangle/\partial T)_{P}$ and iii) the isobaric specific heat $C_{P}\equiv(\partial\langle H\rangle/\partial T)_{P}$. In the previous expressions the symbol $\langle...\rangle$ refers to the thermodynamic average, $V$, $P$, $T$, and $H$ refer to the volume, pressure, temperature and enthalpy of the system respectively, with $H\equiv E+PV$, where is the $E$ energy. The quantities $K_{T}$, $\alpha_{P}$ and $C_{P}$ can be expressed as linear combination of the susceptibilities $\chi_{1}$, $\chi_{2}$ and $\chi_{12}$. Note that the generalized susceptibilities $\chi_{1}$ and $\chi_{2}$, rather than the response functions $K_{T}$ and $C_{P}$, are expected to diverge approaching to the critical point with the characteristic exponents $\gamma$ and $\alpha$ respectively.

At the critical point all the scaling fields vanishes

[TABLE]

while the coexistence line is given by

[TABLE]

The equation (27) defines also a locus that emanates from the critical point into the supercritical region as analytical continuation of the coexistence line. Such locus is, by definition, the Widom line. The line $h_{1}=0$ corresponds, by construction, also to the line of maxima of the statistical correlation length of the system $\xi$. Indeed, it corresponds to the locus where the ordering field vanishes, allowing the fluctuations to spread over broader distances.

According to Holten et al. [53], in the vicinity of the liquid-liquid critical point $(T_{c},P_{c})$, the scaling fields $h_{1}$ and $h_{2}$ can be expressed as linear combination of the physical fields $P$ and $T$

[TABLE]

[TABLE]

where

[TABLE]

with $V_{c}$ the critical volume and $R$ the gas universal constant. The coefficient $a^{\prime}\equiv-d\hat{T}/d\hat{P}$ of Eq. (28) is the slope in the $T$–$P$ plane of the coexistence line (or Widom line) at the critical point, while the coefficient $b^{\prime}$ is the mixing-coefficient of the mixed-field order parameter that accounts for the lack of symmetry in the critical density distribution [112], [44].