Nuclear quantum effects on the thermodynamic response functions of a polymorphic waterlike monatomic liquid

Yizhi Liu ,1,* Gang Sun ,2,* Ali Eltareb ,3,5,* Gustavo E. Lopez,4,5,† Nicolas Giovambattista ,3,5,† and Limei Xu1,6,† 1International Center for Quantum Materials, School of Physics, Peking University, Beijing 100871, China 2School of Chemistry, University of Sydney, Sydney, New South Wales 2006, Australia 3Department of Physics, Brooklyn College of the City University of New York, Brooklyn, New York 11210, USA 4Department of Chemistry, Lehman College of the City University of New York, Bronx, New York 10468, USA 5The Graduate Center of the City University of New York, New York, New York 10016, USA 6Collaborative Innovation Center of Quantum Matter, Beijing, China


I. INTRODUCTION
The anomalous behavior of water has been noticed for at least more than 350 years [1] and, at present, more than 70 anomalies have been identified [2]. For example, upon cooling at low pressures (approximately P < 200 MPa) water exhibits an anomalous increase in thermodynamic response functions, such as isobaric heat capacity C P (T ) and isothermal compressibility κ T (T ), and a density maximum at ≈277 K. In the 1970s, Angell [3][4][5] and collaborators performed a series of studies showing that water anomalies became more pronounced in the supercooled liquid state of water (metastable liquid water at T < 273 K) and at 0 < P < 200 MPa. In these remarkable experiments, Angell et al. showed that thermodynamic response functions such as κ T (T ) and C P (T ) seemed to diverge at T ≈ 228 K [3,6]. A few scenarios have been proposed to explain water's anomalous behavior in the liquid state; excellent reviews are available (e.g., Refs. [7][8][9]). These scenarios include (i) the stability-limit conjecture [10], (ii) the liquid-liquid critical point (LLCP) scenario [11], (iii) the critical-point-free scenario [12], and (iv) the singularity-free scenario [13]. Scenarios (i) and (iii) predict a divergence in κ T (T ) and C P (T ) upon isobaric cooling at low pressures, while in the case of scenarios (ii) and (iv), κ T (T ) and C P (T ) exhibit a maximum [14,15]. Recent experiments [16][17][18] indicate that upon isobaric cooling, κ T (T ) does not diverge but exhibits a maximum, implying that scenarios (i) and (iii) are not valid. We note that water is also anomalous in the glass state (see, e.g., Refs. [8,[19][20][21][22]), existing in two amorphous solid states at approximately P < 0.8 MPa. Moreover, the two glassy states of water can be interconverted by, e.g., isothermal compression, via an apparent first-order-like phase transition. In this regard, it is important to keep in mind that only the LLCP scenario can satisfactorily explain the sharpness of such a glass-glass first-order-like phase transition [22]. In addition, the LLCP scenario provides a natural framework to understand the anomalous behavior of binary aqueous solutions [23]. We also note that the LLCP scenario is strongly supported by computer simulations of many classical water models where a LLCP is found at low temperatures [9,[24][25][26][27][28][29][30].
In this work, we perform path-integral Monte Carlo (PIMC) simulations of a monatomic waterlike model liquid that exhibits a LLCP and the associated κ T (T )-, C P (T )-, and ρ-maxima lines in the supercritical region of the P-T plane. By tuning the value of the Planck's constant h from zero to positive values we are able to explore the nuclear quantum effects on both the LLCP and supercritical lines as the liquid evolves from classical (h = 0) to quantum (h > 0). This allows us to identify which supercritical lines should (and should not) be considered to predict the shift in the LLCP pressure and temperature due to isotope substitution. Indeed, our results indicate that as the liquid becomes more quantum, the location of the LLCP and the κ T (T )-and C P (T )maxima lines shift towards lower temperatures and higher pressures, while, instead, the location of the ρ-maxima line shifts towards higher temperatures and covers a wider range of pressures.
This work is organized as follows. In Sec. II we present the computer simulation details. The results are included in Sec. III, where we discuss the nuclear quantum effects on the location of the liquid-liquid phase transition (LLPT) and LLCP as well as supercritical maxima lines in the P-T plane. A summary is included in Sec. IV.

II. COMPUTER SIMULATION DETAILS
We perform PIMC simulations [52] of an isotropic pairinteraction potential given by a modified version of the Fermi-Jagla (FJ) potential defined in Ref. [53]. Specifically, the modified potential is defined by the same function that defines the FJ potential, U (r) = 0 a n r n + but with the set of parameters given in Table I. The modified Fermi-Jagla (mFJ) potential is compared to the original FJ potential in Fig. 1. We also include a Lennard-Jones potential parametrized so its minimum coincides with the minimum of the FJ potential. As for the FJ potential, the mFJ model is a core-softened potential characterized by a hard-core radius r = a and an attractive minimum at r ≈ 1.7a. Moreover, both model liquids exhibit many of water's anomalous properties [8,9,18], including the increase in compressibility and the development of a density maximum upon isobaric cooling. In particular, both model liquids also exhibit a LLPT and LLCP, consistent with simulations of several full-atomistic water models [11,28,[54][55][56][57][58][59][60][61][62]. However, the mFJ liquid (but not the FJ liquid [53]) also exhibits maxima in the isobaric heat capacity C P (T ) upon isobaric cooling. This is because the slope of the LLPT line dP/dT | LLPT of the mFJ is positive, while in the FJ model, dP/dT | LLPT ≈ 0. It has been shown that when dP/dT | LLPT = 0, there is no maxima in C P (T ) above the LLCP [63,64]. In the case of water models that exhibit a LLCP, C P (T ) exhibits a maximum in the supercritical region of the P-T plane. It is for this reason that in this work we consider the mFJ model as opposite to the FJ model. We note, however, that in the computer simulations of water models that exhibit a LLPT line, dP/dT | LLPT < 0, i.e., opposite to the case of the mFJ liquid. Our computational techniques are identical to those employed in Ref. [65], where we focused on the FJ model and the nuclear quantum effects on the corresponding LLCP (the thermodynamic response functions in the supercritical region were not explored). Accordingly, we refer the reader to Ref. [65] for details. Briefly, PIMC simulations are performed at constant (N, V, T ) for a system of N = 1000 particles, each composed of n = 10 beads. Following Ref. [65], we consider a family of mFJ liquids by considering different values of the Planck's constant, h. Specifically, we perform simulations for the classical mFJ liquid, corresponding to h = 0, as well as three quantum versions of the mFJ liquid corresponding to h = h 1 = 0.2474, h 2 = 0.515, and h 3 = 0.7948 (these values of h are not negligible and are relevant in the context of real substances such as hydrogen and water which have been suggested to exhibit a LLCP). It follows that the mFJ liquids are more quantum with increasing h. At a given (N, V, T ), PIMC simulations are performed for 10 5 Monte Carlo (MC) steps during equilibration. This is followed by additional 10 5 − 3 × 10 5 MC steps, depending on the working conditions. For data analysis, we only considered state points for which the mean-square displacement of the ring-polymer centroids is larger than twice the mFJ hard-core radius a. In one MC step, we first move all the 10 000 beads of the system and then the centroids of all N ring-polymers are displaced. We use reduced units where the particle mass is m = 1 and the Boltzmann constant k B = 1; the energy and distance units are given by the mFJ potential parameters 0 and a, respectively. It follows that the units of h are a( 0 m) 1/2 [66]. We note that our work addresses the effects of quantum mechanics on a given classical model (i.e., by tuning the Planck's constant from h = 0 to h > 0). Accordingly, the underlying interaction potential between particles (mFJ pair potential) is not altered as the liquid becomes more quantumlike. In a real liquid, increasing h could also lead to additional changes in the effective interacting potential. Hence, in the context of this work, these additional effects should be considered to be negligible or secondary. Such effects may include changes at the electron distribution level in the atoms/molecules of the system and/or changes in the "geometry" of the molecules considered. At present, it is unclear how relevant such effects may be and how much they depend on the system studied. The evaluation of the quantum effects on the underlying interacting potential considered falls beyond the scope of this work and requires theoretical and computational techniques other than PIMC simulations.
From the PIMC simulations we extract the thermodynamic properties of the ring-polymer system, such as the P(v) isotherms, from which κ T can be obtained. The case of C P is subtle. From the PIMC simulations we obtain the heat capacity of the ring-polymer system C P (T ) RP = (∂H/∂T ) P . The heat capacity of the corresponding quantum liquid is given by C P (T ) = C P (T ) RP + C P , where C P is given by is the potential energy of the ring-polymers' springs at the given working conditions (N, V, T ); see, e.g., Ref. [67] for details.

III. RESULTS
For a given value of h, we perform PIMC simulations for volumes v = V/N expanding within the range 1.0 < v < 3.2 and at temperatures T = 0.10, 0.15, 0.20 . . . 0.95. From the PIMC simulations data, we construct the corresponding P(v) isotherms. As an example, we include in Fig. 2(a) all the P(v) isotherms obtained for h = h 2 . As shown in the Supplementary Material [68], similar results are obtained at all the values of h studied. The main point of Fig. 2(a) is the development of van der Waals loops in the P(v) isotherms   upon cooling, which indicates the existence of a LLPT and associated LLCP at low temperature. The corresponding P-T phase diagram (h = h 2 ) showing the LLCP and LLPT line is shown in Fig. 2(b). P(v) isotherms similar to those shown in Fig. 2(a) and a P-T phase diagram analogous to Fig. 2(b) were also found for the classical and quantum FJ liquids in Refs. [53,65]. Accordingly, we refer the reader to these works for a detailed discussion of the phenomenology associated to the LLPT and LLCP. As shown in Refs. [53,65], the two liquids at T < T c differ in their local structure (and hence, in thermodynamic properties). As shown in Fig. 3, in the high-density liquid (HDL) the nearest neighbors of a given Centroid-centroid RDF for the quantum mFJ liquid at h = h 2 and T = T c , and for different volumes. As the liquid evolves from the LDL-like (v > v c ≈ 2.1) to the HDL-like (v < v c ) state, the peak of g cc (r) at r/a ≈ 1.7 decreases while an extra peak develops at r/a ≈ 1, consistent with structural changes found in the FJ liquid in the LDL and HDL states [53,65]. The inset shows the probability distribution P(r bc ) to find a bead of a given ring polymer (atom) at a distance r bc from the corresponding centroid. P(r bc ) varies weakly with v and is nonzero up to r bc /a ≈ 0.15, indicating that atoms for the mFJ liquid at h = h 2 are weakly delocalized. See Supplemental Material [68] for results obtained at h = h 1 , h 3 . atom are at a distance r nn ≈ a, while in the low-density liquid (LDL) r nn ≈ 1.7a, which corresponds approximately to the minimum of the mFJ pair potential (Fig. 1).
The LLCP temperature T c for the different quantum mFJ liquids is shown in Fig. 4(a). As for the original FJ liquids, the LLCP of the mFJ liquids move towards lower temperatures as the liquid becomes more quantum (i.e., as h increases). In addition, we find the same nuclear effects on critical pressure and critical volume of the FJ and mFJ. As shown in Figs. 4(b) and 4(c), P c increases while v c remains constant as h increases. We note that our results indicate that the nuclear quantum effects on (T c , P c , v c ) are not related to the slope of the LLPT line in the P-T plane, i.e., whether the slope of the LLPT line is zero (FJ liquid) or positive (mFJ liquid). This is important since the slope of the LLPT line determines the behavior of the maxima lines (in the P-T plane) of thermodynamic response functions, such as the heat capacity and compressibility maxima lines.
Interestingly, the slope of the LLPT line can be affected by the nuclear quantum effects. As shown in Fig. 4(d), the slope of the LLPT line, dP/dT | LLCP , decreases with increasing h. The Clausius-Clapeyron relationship [69] states that dP/dT | LLCP = S/ V , where S = S HDL − S LDL is the difference in entropy of HDL and LDL and V = V HDL − V LDL is the corresponding difference in volume. Hence, it follows from Fig. 4(d) that as the liquid becomes more quantum, the ratio S/ V decreases. For all values of h studied, we find that −1.0 < V < −0.1, depending on the specific values of T < T c and h. It follows that since V < 0 and dP/dT | LLCP > 0 for all h, we also have S < 0. Hence, HDL is always more ordered (lower entropy) than LDL [15]. Interestingly, we find that the entropy difference between HDL and LDL becomes less pronounced as the liquid becomes more quantum (i.e., | S| decreases with increasing h).
Next we discuss how the lines of maxima density ρ, isobaric heat capacity C P , and isothermal compressibility κ T change location in the P-T plane as the mFJ liquid evolves from classical (h = 0) to quantum (h > 0). The results are summarized in Fig. 5. Figure 5(a) shows that the C P -maxima line shifts towards lower temperatures and higher pressures as h increases, consistent with the corresponding shift in the location of the LLCP. Similarly, increasing h shifts the κ Tmaxima lines towards lower temperatures and higher pressures; see Fig. 5(b). The effect of increasing h on the ρmaxima line is shown in Fig. 5(c). Surprisingly, the shift in the ρ-maxima line is different from the shift in the location of the LLCP. Specifically, as the liquid becomes more quantum, the LLCP shifts towards lower temperatures and larger pressures, while the ρ-maximum line extends towards higher temperatures and covers a wider range of pressures. This means that, for example, in the case of D 2 O and H 2 O, the relative difference in pressure and temperature of the corresponding ρ-maxima lines should not be used to estimate the relative location of the corresponding LLCPs.
It follows from Figs. 5(a) and 5(c) that nuclear quantum effects not only shift the location of the ρand C P -maxima lines but can also change their shape. This implies that the ρand C P -maxima lines among isotopes, e.g., between H 2 O and D 2 O, cannot be used to predict the relative location of the corresponding LLCPs. The case of the κ T -maxima lines is more subtle. The lower κ T -maxima line in Fig. 5(b) is sensitive to h while the upper κ T -maxima line is not. Computer simulations of water models show only one κ T maxima, and The insets show the quantities in the main panels shifted in pressure and temperature relative to the LLCP values. As the liquid becomes more quantum (i.e., h increases), the C P -and κ T -maxima lines shift towards lower temperatures and higher pressures, as for the case of the LLCP. This is not the case for the ρ-maxima line. this line joins the C P -maxima line in the proximity of the LLCP, defining the so-called Widom line (see, e.g., Ref. [9]). In other words, the κ T -maxima line relevant to the case of water (and experiments) corresponds to the upper κ T -maxima line in Fig. 5(b). In the case of the mFJ liquids, quantum nuclear effects shift this κ T -maxima line in the P-T plane but do not affect its shape. If this is also the case of water, then one could indeed determine the relative locations of the LLCP of H 2 O and D 2 O by experimentally measuring the relative shift in pressure and temperature of the corresponding κ T -maxima line. At present, the κ T -maxima line for H 2 O and D 2 O has been measured at P = 0.1 MPa, with maxima observed at 229 and 233 K, respectively [16]. This shift in temperature is not inconsistent with the estimated LLCP temperatures of H 2 O and D 2 O, respectively, T c ≈ 220-223 K (P c ≈ 50-100 MPa) [22,70], and T c = 230 ± 5 K (P c ≈ 50 ± 20 MPa) [51].

IV. SUMMARY AND DISCUSSION
We performed PIMC simulations of the mFJ liquid over a wide range of volumes and temperatures using different values of Planck's constant. The mFJ model is a waterlike, monatomic model liquid that exhibits a LLPT and LLCP, as well as lines of maxima in κ T , C P , and ρ in the supercritical region of the P-T plane. The same anomalous properties are predicted by different classical water models.
We find that adding nuclear quantum effects to the classical mFJ liquid shifts the location of the LLPT and LLCP towards higher pressures and lower temperatures. Similar effects were found in the case of the original FJ liquid [65]. The FJ and mFJ liquids have LLPT lines in the P-T plane with different slopes. In the FJ liquid, dP LLCP /dT ≈ 0, while dP LLCP /dT > 0 in the case of the mFJ liquid. This suggests that the reported nuclear quantum effects on the LLCP/LLPT are rather general, independent of dP LLCP /dT . This is important since, in the case of water, computer simulations indicate that dP LLCP /dT < 0.
We also find that the nuclear quantum effects shift the location of the κ T -and C P -maxima lines towards higher pressures and lower temperatures. It follows that these maxima lines can be used to estimate the nuclear quantum effects on the LLCP/LLPT. Indeed, experiments in H 2 O and D 2 O show that the maximum in κ T at P = 1 bar shift by −4 K in H 2 O relative to D 2 O [16], i.e., as nuclear quantum effects become more pronounced. This shift in temperature is also consistent with the estimated location of the LLCP in H 2 O (T c ≈ 220 − 223 K) and D 2 O (≈230 ± 5 K) [50,51]. Instead, we find that the nuclear quantum effects on the ρ-maxima line are different to those observed on the LLCP/LLPT. Specifically, the ρ-maxima line shifts towards higher temperatures and covers a wider pressure range as the mFJ becomes more quantum. It follows that the nuclear quantum effect on the LLCP/LLPT should not be estimated by the corresponding effect on the ρ line.
Our results may be important for experiments in water and hydrogen searching for evidence of a LLCP/LLPT. In both cases, computer simulations predict the existence of a LLCP/LLPT. However, accessing the LLCP/LLPT has been challenging in experiments. In the case of water, crystallization occurs rapidly in the proximity of the hypothesized LLCP. In the case of hydrogen, the LLCP is predicted to exist at very high pressures and temperatures. In both cases, supercritical lines in thermodynamic response functions are predicted by computer simulations, and these maxima lines are easier to access experimentally than the LLCP itself. In this regard, our results provide an additional test of the LLCP hypothesis in experiments. The present PIMC simulations suggest how isotope substitution should shift the supercritical lines (and LLCP/LLPT).