Critical behavior of a water monolayer under hydrophobic confinement

The properties of water can have a strong dependence on the confinement. Here, we consider a water monolayer nanoconfined between hydrophobic parallel walls under conditions that prevent its crystallization. We investigate, by simulations of a many-body coarse-grained water model, how the properties of the liquid are affected by the confinement. We show, by studying the response functions and the correlation length and by performing finite-size scaling of the appropriate order parameter, that at low temperature the monolayer undergoes a liquid-liquid phase transition ending in a critical point in the universality class of the two-dimensional (2D) Ising model. Surprisingly, by reducing the linear size L of the walls, keeping the walls separation h constant, we find a 2D-3D crossover for the universality class of the liquid-liquid critical point for , i.e. for a monolayer thickness that is small compared to its extension. This result is drastically different from what is reported for simple liquids, where the crossover occurs for , and is consistent with experimental results and atomistic simulations. We shed light on these findings showing that they are a consequence of the strong cooperativity and the low coordination number of the hydrogen bond network that characterizes water.

and the temperature of minimum density (TminD). The TMD locus merges the TminD line as in experiments 52 and other models 53 .
At low T a discontinuous change in r is observed for 1wPv 0 = 4E ð Þ §0:5, where the parameters v 0 and E are explained in the Methods section, as it would be expected in correspondence of the HDL-LDL phase transition. At very high pressures (Pv 0 = 4E ð Þw1) the system behaves as a normal liquid, with monotonically increase of r upon decrease of T.
We estimate the liquid-to-gas (LG) spinodal at Pv 0 = 4E ð Þv0 for low T (Fig. 2) as the temperature along an isobar at which we find a discontinuous jump of r to zero value by heating the system. The LG spinodal identifies the locus of the stability limit of liquid phase with respect to the gas phase: at pressures below the LG spinodal in the P 2 T plane is no longer possible to equilibrate the system in the liquid phase. The LG spinodal continues at positive pressures ending in the LG critical point (data not shown). We observe that the TMD line approaches the LG spinodal, without touching it (Fig. 2). We recover the LG spinodal also as envelope of isochores (Fig. 2b).
We find a second envelope of isochores at lower T and higher P, pointing out the liquid-to-liquid (LL) spinodal. Indeed, the two spinodals associated to the LLPT, i.e. the HDL-to-LDL spinodal and the LDL-to-HDL spinodal, collapse one on top of the other and are indistinguishable within our numerical resolution. Nevertheless, we clearly see that isochores are gathering around the points (Tk B = 4E ð Þ*0:06, Pv 0 = 4E ð Þ*0:5) and (Tk B = 4E ð Þ~0, Pv 0 = 4E ð Þ~1), where k B is the Boltzmann constant, marking two possible critical regions (Fig. 2b).
We calculate the isothermal compressibility by its definition K T ; 2(1/AEVae) (hAEVae/hP) T and by the fluctuation-dissipation theorem K T 5 AEDV 2 ae/k B TV along isobars, K T (T), and along isotherms, K T (P) (Fig. 3), where AEVae ; V is the average volume and AEDV 2 ae the volume fluctuations. We find two loci of extrema for each quantity K T (T) and K T (P): one of strong maxima and one of weak maxima. The loci of strong maxima in K T (T) and K T (P), respectively K sMax We find also loci of weak maxima, K wMax T T ð Þ and K wMax T P ð Þ and minima K min T T ð Þ and K min T P ð Þ. The loci of weak extrema and minima of K T (T) and K T (P) do not coincide in the T 2 P plane. The locus of weak maxima along isotherms K wMax T P ð Þ merges with the locus of minima K min T P ð Þ at the point where the slope of both loci is hP/hT R '. Furthermore, both loci approach to the LL spinodal at high P. The locus of weak maxima along isobars K wMax T T ð Þ approaches the LL spinodal where K T exhibits the strongest maxima, and merges with the locus of minima K min T T ð Þ where the slope of both loci is hP/hT R 0 (data at high P and T not shown in Fig 3). This locus intersects the TMD at its turning point. Indeed, as reported in Ref. 39 and in the Methods section, the temperature derivative of isobaric K T along the TMD line is related to the slope of TMD line where all the quantities are calculated along the TMD line. Hence the locus of extrema in K T (T), where (hK T /hT) P 5 0, crosses the TMD line where the slope (hP/hT) TMD is infinite. We observe also that the weak maxima of K T (T) and K T (P) increase as they approach the LL spinodal. All loci of extrema in K T are summarized in Fig. 3. Next we calculate the isobaric specific heat C P ; (hAEHae/hT) P 5 AEDH 2 ae/k B T along isotherms and isobars, where is the Hamiltonian as defined in the Methods section, AEDH 2 ae is the enthalpy fluctuations (Fig. 4). We find two maxima at low P separated by a minimum. At high-T the maxima are broader and weaker than those at low-T. As discussed in Ref. 49, the maxima at high T are related to maxima in fluctuations of the HB number N HB , while the maxima at low T are a consequence of maxima in fluctuations of the number N coop of cooperative HBs. The lines of strong maxima at constant P and constant T, respectively C sMax P T ð Þ and C sMax P P ð Þ, overlap for all the considered pressures, and both maxima are more pronounced in the range Pv 0 = 4E ð Þ[ 0:5, 0:6 ½ and Tk B = 4E ð Þ[ 0:06, 0:07 ½ . The weak maxima C wMax P P ð Þ and C wMax P T ð Þ increase approaching the LL spinodal and have their larger maxima at the state point where they converge to the strong maxima, consistent with the occurrence of a critical point for a finite system (Fig. 4). The lines of weak maxima overlap for all positive pressures, branching off at negative pressures. At negative pressures, the locus C wMax P P ð Þ bends  . P increases from 20.5 (bottom curve) to 1.5 ð4EÞ=v 0 (top curve). Along each isobar we locate the maximum r (green squares at high T) and the minimum r (green small circles at low T) and the liquid-gas spinodal (open large circles at low P). (b) Loci of TMD, TminD, liquid-gas spinodal and liquid-liquid spinodal in T 2 P plane. Dashed lines with labels represent the isochores of the system from rv 0 5 0.43 (bottom) to rv 0 5 0.80 (top). Dashed lines without labels represent intermediate isochores. TMD and TminD correspond to the loci of minima and maxima, respectively, along isochores in the T 2 P plane. We estimate the critical isochore at rv 0 , 0.47 (red circles). All the isochores with 0.47 , rv 0 , 0.76 intersect with the critical isochore for Pv 0 = 4E ð Þ §0:5 along the LL spinodal (tick turquoise) line.

LP LT
in case of intersection between the locus of extrema (hC P /hP) T 5 0 and the TMD line, it results that (hP/hT) TMD 5 0. Note that, as we explain in the Methods section, the relation (2) does not imply any change in the slope of the TminD line at the intersection with the locus of (hC P /hP) T 5 0. We calculate also the thermal expansivity a P ; (1/AEVae) (hAEVae/hT) P along isotherms and isobars (Fig. 5). As for the other response functions, we find two loci of strong extrema, minima in this case, a smin P P ð Þ and a smin P T ð Þ, along isotherms and isobars, respectively showing a divergent behavior in the same region where we find the strong maxima of K T and C P . From this region two loci of weaker minima depart. We find that the locus of weak minima along isobars a wmin P T ð Þ bends toward the turning point of the TMD. Although our calculations for a P do not allow us to observe the crossing with the TMD line, based on the relation (see Methods) that holds at the TMD line, we can conclude that a wmin P T ð Þ should have zero T-derivative if it crosses the point where the TMD turns into the TminD line, because in this point the TMD slope approaches zero.
The locus of weaker minima along isotherms a wmin P P ð Þ, merges with the locus of maxima a Max P P ð Þ at the state point where the slope of both loci is hP/hT R ' (not shown in Fig. 5). According to the thermodynamic relation, discussed in Methods section, we find that the locus of extrema in thermal expansivity along isotherms coincides, within the error bars, with the locus of extrema of isothermal compressibility along isobars (Fig. 5c). All the loci of extrema of response functions that converge toward the same region A in Fig. 3, 4 and 5 increase in their absolute values. Because the increase of response functions is related to the increase of fluctuations and this is, in turn, related to the increase of correlation length j, to estimate j we calculate the spatial correlation function wherer i is the position of the molecule i,r i {r l j j~r the distance between molecule i and molecule l and AE?ae the thermodynamic average. The states of the water molecule, as well as the density r, the energy E and the entropy S of the system, are completely described by the bonding variables s ij . Therefore, the function G(r) accounts for the fluctuations in r, E and S and allows us to evaluate the correlation length because the order parameter of the LLPT, as we discuss in the following, is related to a linear combination of r and E. Note that, instead, the density-density correlation function would give only an approximate estimate of j.
We observe an exponential decay of G(r) , e 2r/j at high temperatures in a broad range of pressures. Approaching the region A, the correlation function can be written as G(r) , e 2r/j /r d221g where d is the dimension of the system and g a (critical) positive exponent. When j is of the order of the system size, the exponential factor approaches a constant leaving the power-law as the dominant contribution for the decay.
At P below the region A, we find that j has a maximum, j Max , along isobars and that j Max increases approaching A (Fig. 6). The j Max locus coincides with the locus of strong extrema of C P , K T and a P (Fig. 6b). We observe that this common locus converges to A and that all the extrema increase approaching A. This behavior is consistent with the identification of A with the critical region of the LLCP. Furthermore, we identify the common locus with the Widom line that, by definition, is the j Max locus departing from the LLCP in the one-phase region 54,55 . Our calculations allow us to locate the Widom line at any P down to the liquid-to-gas spinodal.
At P above the region A, we find the continuation of the j Max line, but with maxima that decrease for increasing P, as expected at the LL spinodal that ends in the LLCP (Fig. 6). Therefore, we identify the high-P part of the j Max locus with the LL spinodal. Along this line the density, the energy and the entropy of the liquid are discontinuous, as discussed in previous works 31,40,[44][45][46][47][48][49] .
To better locate and characterize the LLCP in A we need to define the correct order parameter (o.p.) describing the LLPT. According to mixed-field finite-size scaling theory 56 , a density-driven fluid-fluid phase transition is described by an o.p. M ; r* 1 su*, where r* 5 rv 0 represents the leading term (number density), u:E= EN ð Þ is the energy density (both quantities are dimensionless) and s is the mixed-field parameter. Such linear combination is necessary in order to get the right symmetry of the o.p. distribution Q N (M) at the critical is the critical exponent that governs M, n is the critical exponent that governs j, with n and b defined by the universality class, a M is a nonuniversal system-dependent parameter andp d is an universal function characteristic of the Ising fixed-point in d dimensions. We adjust B and M c so that Q N (M) has zero mean and unit variance.
We combine, using the multiple histogram reweighting method 57 described in the Methods section, a set of 3 3 10 4 MC independent configurations for , 300 state points with 0:040ƒTk B = 4E ð Þƒ0:065 and 0:40ƒPv 0 = 4E ð Þƒ0:75. We verify, by tuning s, T and P, that there is a point within the region A where the calculated Q N (x) has a symmetric shape with respect to x 5 0 (Fig. 7). We find s 5 0.25 6 0.03 for our range of N. The resulting critical parameters T c (N), P c (N) and the normalization factor B(N) follow the expected finite-size behaviors with 2D Ising critical exponents 56 . From the finite-size analysis we extract the asymptotic values T c k B = 4E ð Þ~0:0597+0:0001 and P c v 0 = 4E ð Þ~0:555+0:002. The presence of a first order phase transition ending in a critical point, associated to the o.p. M, is confirmed by the finite size analysis of the Challa-Landau-Binder parameter 58  where the symbol AE?ae N refers to the thermodynamic average for a system with N water molecules. U M quantifies the bimodality in Q N (M). The isobaric value of U M shows a minimum at the temperature where Q N (M) mostly deviates with respect to a symmetric distribution (Fig. 8). Minimum of U M converges to 2/3 in the thermodynamic limit away from a first order phase transition, while it approaches to a value ,2/3 where the bimodality of Q N (M) indicates the presence of phase coexistence.
These results are consistent with the behavior of the Gibbs free energy G calculated with the histogram reweighting method (Fig. 9). In particular, we calculate G along isotherms, for P crossing the LLPT and the loci of weak maxima in K T (T) and C P (P). We find that the behavior of G for T , T c is consistent with the occurrence of a discontinuity in volume V 5 hG/hP, in the thermodynamic limit, with a decrease of V corresponding to the transition from LDL to HDL for increasing P. Crossing the loci K T ðTÞ wMax and C P ðPÞ wMax the volume decreases with pressure without any discontinuity as expected in the one-phase region.
The distribution Q N (N) adjust well to the data only for large N. We, therefore, perform a more systematic analysis. For each N, we quantify the deviation of the calculatedp N ð Þfrom the expectedp 2 for the 2D Ising. Furthermore, due to the behavior of data for small N (Fig. 7a), we calculate the deviation from the 3D Isingp 3 56 . We estimate the Kullback-Leibler divergence 51,59 , of the probability distributionp i N ð Þ of x i from the theoretical valuẽ p d,i of x i (i 5 1, …, n) in d dimensions (Fig. 10a), and the Liu et al. deviation 51 , withp d,peak {p d,x~0 difference between the distribution peak and its value at x 5 0 (Fig. 10b). We confirm s^0:25 forp 2 and find s 5 0.10 6 0.02 forp 3 for our range of N. For both D KL d and W d , with d 5 2 and d 5 3, we find minima at T c k B = 4E ð Þ^0:06 and P c v 0 = 4E ð Þ^0:55 that become stronger for increasing N. We find that D KL 2 and W 2 decrease with increasing N, vanishing for N R ' (Fig. 10). Therefore, for an infinite monolayer between hydrophobic walls separated by h < 0.5 nm, the system has a LLCP that belongs to the 2D Ising universality class, as expected from our representation of the system as the 2D projection of the monolayer.
However, by increasing the confinement, i.e. reducing N and L at constant r, D KL 2 and W 2 become larger than D KL 3 and W 3 , respectively. Therefore, the calculatedp N ð Þdeviates from the 3D probability distribution less than from the 2D probability distribution. For N 5 2500 we find that both D KL 3 and W 3 have values approximately equal to those for D KL 2 and W 2 calculated for a system ten times larger. In particular we find D KL 3^0 for N 5 2500. Hence, by increasing the confinement of the monolayer at constant r, the LLCP has a behavior that approximates better the bulk [25][26][27][28][29][30]38 , with a crossover between 2D and 3D-behavior occurring at Þof configurations of N water molecules with energy and volume V at the LLCP. This quantity is expected to scale as DG!N d{1 d . We find that our data can be fitted as N 2 3 for small sizes and as N 1 2 for large sizes with a crossover around N 5 10 4 (Fig. 10c). Considering the value of the estimated r c in real units (^1g cm 3 ) 45 , the corresponding crossover wall-size is L^25 nm.

Discussion
Our rationale for this dimensional crossover at fixed h is that, when L/h decreases toward 1, the characteristic way the critical fluctuations spread over the system, i.e. the universality class of the LLCP, resembles closely the bulk because the asymmetry among the three spatial dimensions is reduced. A similar result was found recently by Liu et al. for the gas-liquid critical point of a Lennard-Jones (LJ) system confined between walls by fixing L and varying h 51 . However, in the case considered by Liu et al. the crossover was expected because the number of layers of particles was increased from one to several, making the system more similar to the isotropic 3D case. Here, instead, we consider always one single layer, changing the proportion L/h by varying L. Therefore, it could be expected that the system belongs to the 2D universality class for any L. Furthermore, the extrapolation of the results for the LJ liquid to our case of a monolayer with h=r 0^1 :7, where r 0 is the water van der Waals diameter, would predict a dimensional crossover at L=h^5 51 . Here, instead, we find the crossover at L=h^50, i.e. one order of magnitude larger than the LJ case. We ascribe this enhancement of the crossover to (i) the presence of a cooperative HB network and (ii) the low coordination number that water has in both the monolayer and the bulk. These are the main differences between water and a LJ fluid. The cooperativity intensifies drastically the spreading of the critical fluctuations along a network, contributing to the effective dimensionality increase of the confined monolayer. Moreover, the HB network has in 3D a coordination number (z 5 4) as low as in 2D, making the first coordination shell similar in both dimensions.
Our findings are consistent with recent atomistic simulations of water nanoconfined between surfaces. 60 60 . This is confirmed by our study, where the water-interface interaction is purely due to excluded volume. Following the authors of Ref. 60, this observation allows us to relate our finding for rigid surfaces to experimental results for water hydrating membranes 63 , reporting new types of water dynamics in thin interfacial layers, and water nanoconfined in different types of reverse micelles 64 , showing that the water dynamics is governed by the presence of the interface rather than the details (e.g., the presence charged groups) of the interface.
In conclusion, we analyze the low-T phase diagram of a water monolayer confined between hydrophobic parallel walls of size L separated by h < 0.5 nm. We study water fluctuations associated to the thermodynamic response functions and their relations to the loci of TMD, TminD. For each response function we find two loci of extrema, one stronger at lower-T and one weaker and broader at higher-T. These loci converge toward a critical region where the fluctuations diverge in the thermodynamic limit, defining the LLCP. We calculate the Widom line departing from the LLCP based on its definition as the locus of maxima of j and show that it coincides with the locus of strong maxima of the response functions. We find that the LLCP belongs to the 2D Ising universality class for L R ', with strong finite-size effects for small L. Surprisingly, the finitesize effects induce the LLCP universality class to converge toward the bulk case (3D Ising universality class) already for a system with a very pronounced plane asymmetry, i.e. a water monolayer of height h < 0.5 nm and L/h < 50. For normal liquid, instead, this is expected only for much smaller relative values of L (L/h # 5). We rationalize this result as a consequence of two properties of the HB network: (i) its high cooperativity, that enhances the fluctuations, and (ii) its low coordination number, that makes the first coordination shell for the monolayer and the bulk similar.

Methods
The model. We consider a monolayer formed by N water molecules confined in a volume V ; hL 2 between two hydrophobic flat surfaces separated by a distance h, with V=N §v 0^4 2 A 03 , where v 0 is the water excluded volume. Each water molecule has four next-neighbours 7 . We partition the volume into N equivalent cells of height h^0:5nm and square section with size r: ffiffiffiffiffiffiffiffiffiffiffi L 2 =N p , equal to the average distance between water molecules. By coarse-graining the molecules distance from the surfaces, we reduce our monolayer representation to a 2D system. We use periodic boundary conditions parallel to the walls to reduce finite-size effects. We simulate constant N, P, T, allowing V(T, P) to change, with each cell i 5 1, …, N having number density r i :r T,P ð Þ:N=Vƒr 0 :1=v 0 corresponding to a mass density^1 g cm 3 . To each cell we associate a variable n i 5 0 (n i 5 1) depending if the cell i has r i /r 0 # 0.5 (r i /r 0 . 0.5). Hence, n i is a discretized density field replacing the water translational degrees of freedom. The water-water interaction is given by The first term, summed over all the water molecules i and j at O-O distance r ij , has U(r) ; ' for rvr 0 : ffiffiffiffiffiffiffiffiffi v 0 =h p~2 :9 A 0 (water van der Waals diameter), U r ð Þ:4E r 0 =r ð Þ 12 { r 0 =r ð Þ 6 Â Ã for r $ r 0 with E:5:8 kJ=mol, and U(r) ; 0 for r . r c ; 25r 0 (cutoff).
The  The third term of the Eq.(9) accounts for the HB cooperativity due to the quantum many-body interaction 65 , with J s =4E:0:05 and N coop : X i n i X l,k ð Þ i d sik ,sil , where (l, k) i indicates each of the six different pairs of the four indices s ij of a molecule i. The value J s =J is chosen in such a way to guarantee an asymmetry between the two components of the HB interaction. To the cooperative term is due the O-O-O correlation that locally leads the molecules toward an ordered configuration. In bulk water this term would lead to a tetrahedral structure at low P up to the second shell, as observed in the experiments 66 . An increase of T or P partially disrupts the HB network and induces a more compact local structure, with smaller average volume per molecule. Therefore, for each HB we include an enthalpy increase Pv HB , where v HB /v 0 5 0.5 is the average volume increase between high-r ices VI and VIII and low-r (tetrahedral) ice Ih. Hence, the total volume is V ; V 0 1 N HB v HB , where V 0 $ Nv 0 is a stochastic continuous variable changing with Monte Carlo (MC) acceptance rule 46 . Because the HBs do not affect the n.n. distance 66 , we ignore their negligible effect on the U(r) term. Finally, we model the water-wall interaction by excluded volume.
The observables. The LLCP is identified by the mixed-field order parameter M and not by the magnetization of the Potts variables s i,j as in normal Potts model. M is related to the configuration of the system by the relation where v ; V 0 /N and s is the mixed-field parameter. M is therefore a linear combination of density and energy.
Thermodynamic response functions are calculated from and C P : The Monte Carlo method. The system is equilibrated via Monte Carlo simulation with Wolff algorithm 46 , following an annealing procedure: starting with random initial condition at high T, the temperature is slowly decreased and the system is reequilibrated and sampled with 10 4 4 10 5 independent configurations for each state point. The thermodynamic equilibrium is checked probing that the fluctuationdissipation relations, Eq. (11) and (12), hold within the error bar.
The histogram reweighting method. The probability Q N (M) is calculated in a continuous range of T and P across the j Max line. We consider an initial set of m g [10520] independent simulations within a temperature range DTk B = 4E ð Þ*10 {4 and a pressure range DPv 0 = 4E ð Þ*10 {3 . For each simulation i 5 1, …, m we calculate the histograms h i (u, r) in the energy density-density plane. The histograms h i (u, r) provide an estimate of the equilibrium probability distribution for u and r; this estimate becomes correct in the thermodynamic limit. For the NPT ensemble, the new histogram h(u, r, P9, b9) for new values of b9 5 1/k B T9 and P9 close the simulated ones, is given by the relation 57 where N i is the number of independent configurations of the run i. The constants C i , related to the Gibbs free energy value at T i and P i , are self-consistently calculated from the equation 57 e Ci~X u X r h u,r, We choose as initial set of parameters C i 5 0. The parameters C i are recursively calculated by means of Eq. (13) and (14) until the difference between the values at iteration k and k 1 1 is less then the desired numerical resolution (10 23 in our calculations). Once the new histogram is calculated, Q N (M) at T i and P i is calculated integrating h(u, r, P i , b i ) along a direction perpendicular to the line r 1 su.
Thermodynamic relations. We report here the calculations for the thermodynamic relations in Eq. (1), (2), (3) and (4) 39 . To verify the relation (4) we calculate the derivative of K T along isobars and the derivative of a P along isotherms Following 39,67 the line of extrema in density (TMD and TminD lines) is characterized by a P 5 0, hence, da P 5 0 along the TMD line. Therefore, where the index ''ED'' denotes that the derivatives are taken along the locus of extrema in density. So, the slope hP/hT of TMD is given by from which, using Eq. (15) with a P 5 0, we get Eq. (1). The Eq. (18) holds as long as both (ha P /hP) T and (ha P /hT) P do not vanish contemporary, as it occurs along the Widom line, where the loci of strong minima of a P overlap. For this reason the intersection between the Widom line and TminD line does not imply any change in the slope (hP/hT) TminD . To calculate Eq. (2) we start from C P and a P written in terms of Gibbs free energy C P T~{ from which results