Nuclear quantum effects on the dynamics and glass behavior of a monatomic liquid with two liquid states

Glasses are out-of-equilibrium, amorphous solid states. They exhibit no evident long range order, like liquids, and they possess no molecular translational motion, like solids.¹ Most glass-former liquids are classical in nature. This is because common glass former liquids are usually composed of heavy atoms or molecules and exhibit a relatively high glass transition temperature; e.g., T_g ≈ 1000 K in the case of silicon.² At such conditions, quantum mechanics plays no role. However, in the case of light elements, such as H or He, or small molecules that contain H, such as water, nuclear quantum effects (NQE; i.e., quantum fluctuations due to atoms delocalization) can play a relevant role.³ The NQE on low-T glass forming liquids composed of light elements, and on the corresponding glasses, are not well understood. For example, intuitively, one would expect that the inclusion of NQE should lower the glass transition temperature of a liquid. Yet, computer simulations and theory indicate that quantum fluctuations can promote or inhibit the glass formation.^4,5 NQE seem to also play a relevant role in He glasses formed by isothermal compression.⁶ 

Understanding the role of NQE on low-T liquids becomes more complex for the case of polymorphic liquids. This is because polymorphic liquids, such as water and hydrogen, exhibit two liquid states at low temperatures, low-density, and high-density liquid (LDL and HDL). In these substances,^7–9 LDL and HDL are separated by a first-order phase transition line (in the P–T plane) that ends at a liquid–liquid critical point (LLCP) at P > 0. For example, in the case of water, the estimated location of the LLCP is T_c ≈ 220 K and P_c ≈ 50−100 MPa; in the case of hydrogen, T_c ≈ 2000 K and P_c ≈ 100−200 GPa.^10–14 Interestingly, in the case of water, vitrification of LDL and HDL leads to two different glassy states, low-density and high-density amorphous ices (LDA and HDA).^15–22 To the best of our knowledge, LDA and HDA forms have not yet been identified in the case of hydrogen.

Experiments and computer simulations of water indicate that water thermodynamic and dynamic properties are indeed affected by NQE. For example, the diffusion coefficient of H₂O, at a given temperature, is larger than for D₂O, and the estimated LLCP temperature of H₂O is lower than the corresponding temperature of D₂O.^15,23–27 In recent studies,^28,29 we performed path-integral Monte Carlo (PIMC) simulations of the Fermi–Jagla (FJ) liquid, a water-like monatomic model liquid that exhibits a liquid–liquid phase transition (LLPT) and LLCP and studied the NQE on its thermodynamic properties. It was shown that the LLCP of the FJ liquid also shifts toward lower temperatures as NQE become more pronounced, in agreement with experiments and computer simulations of water.^15,24,27,30 Moreover, all the water-like anomalous thermodynamic properties of the FJ liquid, such as the compressibility and density maxima lines, shift toward lower temperatures and change shape with the inclusion of NQE.²⁸ Motivated by these results and the reported NQE effects in water (and, potentially, hydrogen) at conditions close to the corresponding LLPT/LLCP, in this work, we extend our previous studies^28,29 and investigate the role of NQE on the dynamics of the FJ liquid at low-temperatures, close to the corresponding glass state. In the case of the classical FJ liquid,^31,32 as well as water,³³ the system has two glass transition temperatures, one for LDL $[TgLDL(P)]$ and one for HDL $[TgHDL(P)]$⁠.³⁴ These two temperatures exhibit a very different dependence on pressure. One of the goals of this work is (i) to use path integral molecular dynamics (PIMD) simulations and the ring-polymer MD (RPMD) method to study how the introduction of NQE affects both $TgLDL(P)$ and $TgHDL(P)$⁠.

The behavior of $TgLDL(P)$ and $TgHDL(P)$ in water and FJ liquids are intimately related to the presence of LDA and HDA. Hence, in this work, (ii) we also explore the glass behavior of the ring-polymer systems associated with the classical (h = 0) and quantum (h > 0) FJ liquids [via the path-integral (PI) formulation of statistical mechanics]. We note that, strictly speaking, the PI formalism provides an isomorphism between the canonical partition function of the quantum liquid and a classical ring-polymer system composed of beads connected by temperature-dependent springs.^35,36 Accordingly, the PI formalism is limited to equilibrium conditions. Yet, characterizing the behavior of the ring-polymer systems with h ≥ 0 in the glass state is useful to understand the low-temperature behavior of the quantum liquids. Indeed, a few studies have extended the isomorphism between quantum systems and the corresponding ring-polymer counterpart to out-of-equilibrium conditions.^6,37

Our PIMD/RPMD simulations show that including NQE shifts both $TgLDL(P)$ and $TgHDL(P)$ toward lower temperatures, particularly, $TgHDL(P)$ (high-pressure liquid). For all values of h studied, the FJ liquids vitrify into LDA and HDA forms. In addition, the ring-polymers simulations in the glass domain indicate that LDA and HDA forms can be interconverted by isothermal compression/decompression, as expected for polyamorphic substances. Accordingly, our results suggest that the common liquid/glass phase diagram of polyamorphic substances is robust and that glass polymorphism is not removed by the inclusion of NQE. However, microscopically, the LDA–HDA transformation is fundamentally different for the classical-like (h ≈ 0) and quantum-like FJ liquids (h > 0). In particular, for large values of h, the ring-polymers become increasingly expanded during both the LDA-to-HDA and HDA-to-LDA transformations (in the case of small h, the ring-polymers contract along these transitions). These results are important in the context of polyamorphism in light elements, such as hydrogen, for which different studies indicate the existence of a LLPT/LLCP and for the understanding of isotope substitution effects in supercooled and glassy water.

This work is organized as follows: The computer simulation details are provided in Sec. II. In Sec. III, we discuss the NQE on the dynamical properties of the FJ liquid. In Secs. IV and V, we discuss the role of NQE on the glass polymorphism and glass transition of the FJ liquid. A summary is included in Sec. VI.

We perform equilibrium and out-of equilibrium PIMD simulations of a system of particles with isotropic pair interactions given by the Fermi–Jagla (FJ) potential;^29,32 see Fig. S1 of the supplementary material. Following Ref. 32, the FJ potential is truncated at the cutoff distance r_c = 4.0 (all quantities are reported in reduced units³²). The FJ liquid exhibits many of water anomalous properties, such as the increase of the isothermal compressibility and a density maximum upon isobaric cooling and a diffusivity maximum upon isothermal compression. The FJ liquid also exhibits a LLPT and an associated LLCP,³² consistent with computational and experimental studies of water, sulfur, and phosphorus, among others.^38–52 In the glass state, the FJ model shows glass polymorphism with transformations between a low-density and a high-density amorphous solid (LDA and HDA).^31,34 The transformation between LDA and HDA is reminiscent of first-order phase transitions in equilibrium systems, and it is similar to the LDA–HDA transformation observed in water.¹⁵ 

The PIMD simulations technique is based on the path integral formulation of statistical mechanics.^35,36 It allows one to calculate static properties in the presence of NQE, including thermodynamic (e.g., energy and compressibility) and structural properties [e.g., radial distribution functions (RDFs)]. In the path integral formulation, the canonical partition function of the quantum liquid is shown to be isomorphic to the canonical partition function of a classical liquid composed of ring polymers. Each ring-polymer is composed of n_b beads, where the beads are connected by springs with spring constant k_sp = 2πn_b/(βλ²), where $λ=h/2πmkBT$ is de Broglie’s thermal wavelength, β = 1/k_BT, and m is the mass of the particle. The interactions between beads of different polymers are subtle; bead n (n = 1, 2, …, n_b) of ring-polymer A interacts only with bead n of ring polymer B (in our case, the corresponding bead–bead interactions are given by the FJ pair potential re-scaled by 1/n_b).

PIMD simulations cannot be used to calculate dynamical properties, such as the diffusion coefficient of the atoms in the liquid. Hence, for dynamical properties, we use the RPMD method.⁵³ The RPMD approach uses the classical evolution of the ring-polymer system in constant-(N, V, E) PIMD simulations (with the mass of the ring-polymer beads equal to the corresponding atomic mass) to calculate, approximately, the Kubo-transformed correlation functions of the quantum system. From the Kubo-transformed correlation functions, one can obtain dynamical properties of the quantum system. For example, the RPMD approach allows one to calculate the diffusion coefficient of the atoms in the liquid, D, from the Kubo-transformed velocity autocorrelation function evaluated from constant-(N, V, E) PIMD simulations. The use of constant-(N, V, E) PIMD simulations is not strictly necessary. For example, in the thermostatted RPMD (T-RPMD)⁵⁴ method, the Kubo-transformed velocity autocorrelation functions are calculated from constant-(N, V, T) PIMD simulations where the temperature is controlled using a stochastic (local) path integral Langevin equation (PILE) thermostat.⁵⁴ 

In this work, the diffusion coefficient D of the quantum liquid is calculated from constant-(N, V, T) PIMD simulations using a local PILE thermostat, as in the T-RPMD method.⁵⁴ However, in T-RPMD using the local PILE thermostat, one sets the friction coefficient of the zero-frequency mode to zero (γ = 0),⁵⁴ while in our case, γ = 0.1 ps⁻¹. As shown in Ref. 23, the diffusion coefficients calculated in this way are consistent with true T-RPMD as long as the value of γ remains small. In addition, we calculate D from the mean-square displacement (MSD) of the ring-polymer centroids, as implemented in Ref. 55. As shown in Ref. 23, the diffusion coefficient obtained by the method followed here gives values of D that are consistent with those obtained by the original RPMD method, based on the Kubo-transformed velocity autocorrelation functions and calculated from constant-(N, V, E) PIMD simulations.

In order to study the NQE on the properties of the FJ liquid, we study a family of FJ liquids, each liquid differing by the value of Planck’s constant, h, in the expression for k_sp. As the value of h increases, k_sp decreases, allowing the beads to spread. We note that the classical FJ liquid corresponds to the case h = 0 for which k_sp → ∞ and the ring-polymers collapse. In this work, we focus on three values of h: h₀ = 0, h₁ = 0.2474, and h₃ = 0.7948. We note that 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.^28,29 The thermodynamic properties and phase behavior of the quantum FJ liquids with h = h₁, h₃ are discussed in Ref. 29. Equilibrium PIMD simulations are run for n_eq = 5 × 10⁵ to 1.5 × 10⁷ simulation steps during equilibration, depending on temperature and volume. Production runs extend for an additional n_prod = 5 × 10⁵ to 1.5 × 10⁷ number of steps. We confirm that n_eq, n_prod are large enough so the mean-square displacement of the particles is $>10$ (i.e., particles diffuse for at least one nearest-neighbors average separation).

We perform PIMD/RPMD simulations of a system composed of N = 1000 particles in a cubic box with periodic boundary conditions. Most of the PIMD/RPMD simulations are performed with n_b = 10 beads per ring-polymer. However, at very low temperatures, we also run simulations using n_b = 50, 100, 150, 250 to test for the convergence of our results with n_b. Equilibrium PIMD simulations are performed at constant (N, V, T) and extend over a wide range of temperatures and volumes. Out-of-equilibrium simulations, in the low-temperature liquid and glass domains, involve isobaric cooling/heating and isothermal compression/decompression runs. During isobaric cooling/heating, the thermostat temperature changes linearly with time, T(t) = T₀ ± q_T t, where q_T = 0.001 is the cooling/heating rate and t is the time. Similarly, during isothermal compression/decompression runs, the barostat pressure changes linearly with time, P(t) = P₀ ± q_P t, where q_P = 0.0001 is the compression/decompression rate. In all cases, the temperature is controlled using a local PILE thermostat,⁵⁶ while the pressure is controlled using a Monte Carlo (MC) barostat.^57,58 The thermostat collision frequency γ and the frequency of the MC barostat Δ used in the equilibrium and non-equilibrium simulations are important; see Fig. S2 of the supplementary material. In this work, we use γ = 0.215 05 (γ = 0.1 ps⁻¹ in OpenMM default units) and Δ = 1 (which corresponds to performing one MC volume change trial every time step dt = 0.000 232). As shown in Fig. S2 of the supplementary material, for these particular values of γ and Δ, we can reproduce the equilibrium and non-equilibrium phase diagram of the classical FJ liquid (e.g., Fig. 2 of Ref. 31 and Fig. 1 of Ref. 34).

For a given quantum liquid, characterized by a value of h > 0, we study the pressure-induced LDA–HDA transition of the associate ring-polymer system at temperature T₁(h). The starting LDA used for the compression runs are obtained by cooling the liquid at constant pressure P₁(h) from a high temperature, where the liquid is in equilibrium, to the compression/decompression temperature T₁(h). LDA is then compressed from P = P₁(h) up to pressure P₂(h) [at T = T₁(h)] to produce HDA; the so produced HDA is then decompressed to very low pressures [at T = T₁(h)]. The values of (T₁, P₁, P₂) vary with h because changing h shifts the location of the LLCP in the P–T plane.²⁹ We set T₁(h) = 0.0555 T_c(h), where T_c(h) is the LLCP temperature of the FJ liquid characterized by a specific value of h; see Table I. To define P₁(h) and P₂(h), we take advantage of the fact that the LLCP volume is v_c(h) = 2.9 for all the FJ liquids studied (i.e., v_c is independent of h). The value of P₁(h) is chosen so the volume of the starting liquid configuration for all values of h is v ≈ 4.5 > v_c; this corresponds to cooling the liquids (LDL) isobarically at P₁(h₀) = −0.10, P₁(h₁) = −0.10, and P₁(h₃) = −0.15. P₂(h) is chosen so as to keep the same ratio for −P₁/P₂. For the classical FJ glass, P₂(h₀) = 1.5 and, hence, −P₁(h₀)/P₂(h₀) = 15. Accordingly, P₂(h₁) = 1.5 and P₂(h₃) = 2.25. For a given cooling/heating/compression/decompression simulation, we perform ten independent runs, following the same protocol explained above. All of the PIMD simulations are performed using the OpenMM (version 7.4.0) software package,⁵⁹ modified to include the FJ pair potential interactions. In our computer simulations using the OpenMM software package, the atoms mass is m = 39.948 amu, a = 3.40 Å, and ϵ₀ = 0.9978 kJ/mol (but all quantities are reported in reduced units^29,32,34).

In this section, we focus on the diffusion coefficient D of the different FJ liquids studied. From the RPMD simulations, we first calculate the mean-square displacement (MSD) of the ring-polymer centroids, MSD$(t)=|r⃗(t0+t)−r⃗(t0)|2$⁠, where $⋯$ indicates average over particles/ring-polymers and initial times t₀. The diffusion coefficient is obtained from the long term behavior of the MSD, MSD(t) → 6Dt.^23,55

To determine the NQE on D as the FJ liquid becomes increasingly quantum (increasing h), we calculate the iso-diffusivity line in the P–T plane for the case D = 0.019. The (D = 0.019)-line in the P–T plane for the quantum liquid with h = h₃ is shown in Fig. 1(a); as shown in Fig. S3 of the supplementary material, similar results hold for the cases h = h₀, h₁. Also included in Fig. 1(a) are the LLCP reported in Ref. 29 as well as the lines of maxima diffusivity and maxima density (D_max- and ρ_max-lines). These maxima lines are well-known anomalous properties found in water and other anomalous liquids.^60,61 The D_max-line indicates the pressures at which D reaches a maximum upon isothermal compression and the ρ_max-line indicates the temperature at which ρ reaches a maximum upon isobaric cooling.

In normal liquids, the diffusion coefficient at a given temperature decreases monotonically with increasing pressure. Hence, one would expect to find an iso-diffusivity line with a positive slope in Fig. 1(a). Instead, for the classical and quantum FJ liquids considered, the (D = 0.019)-line is negatively sloped at low pressures, in the LDL-like domain, and positively sloped at large pressures, in the HDL-like domain. Indeed, the LDL and HDL branches of the (D = 0.019)-line seem to intersect at the D_max-line. The non-monotonic behavior of the (D = 0.019)-line shown in Fig. 1(a) and its relationship with the D_max-line are fully consistent with previous classical MD simulations of water; see Ref. 33. Specifically, in the case of water, the glass transition temperature of LDL is negatively sloped in the P vs T plane, while the glass transition temperature of HDL is positively sloped. In addition, the glass transition temperature of water reaches a minimum at the intersection with the D_max-line.

As shown in Fig. 1(b), the LLCP and the D_max-line both shift toward lower temperatures with increasing h, i.e., as the liquid becomes more quantum. Interestingly, including NQE also alters the shape of the D_max-line. A similar shift in temperature and change in shape were found in the lines of maxima in thermodynamic properties, such as the ρ_max-line specific heat and compressibility maxima lines;²⁸ see Fig. S3 of the supplementary material.

In Fig. S4 of the supplementary material, it is shown that, in the classical liquid case, the (D = 0.019)-line expands above and below T_c. Instead, in the quantum liquids, the (D = 0.019)-line is always located at T > T_c. Moreover, the separation in temperature between the LLCP and the (D = 0.019)-line increases with increasing h. This implies that while introducing NQE shifts both the iso-diffusivity line and the LLCP toward lower temperatures, NQE are more pronounced on the LLCP. Accordingly, it is possible that at very large h, the quantum liquid still exhibits an anomalous D_max-line, with non-monotonic iso-diffusivity lines, without an accessible LLPT/LLCP.

The LLCP found in the quantum liquids with h = h₁, h₃ and the anomalous shape of the iso-diffusivity line in Fig. 1 are consistent with the possibility that the quantum systems also exhibit glass polymorphism, such as the classical FJ liquid. That glass polymorphism may exist in the quantum case is not evident. For example, crystallization could intervene before the glass state is reached,²⁹ and NQE and/or quantum tunneling could affect the transformation between different glass states at very low temperatures.⁶² Here, we extend the PIMD simulations of the quantum FJ model liquids to the low-temperature, out-of-equilibrium domain. We explore the glass behavior of the ring-polymer systems associated with the quantum FJ liquids studied. Again, while the PI formalism does not guarantee that the ring-polymer system in the glass state and its glass quantum counterpart are isomorphic, studying the glass behavior of the ring-polymer systems can shed light onto the unusual dynamics of the quantum liquids at low temperatures (see Sec. V). Briefly, we show that for all values of h considered, the ring-polymer systems also exhibit glass polymorphism with sharp, first-order-like phase transitions between LDA and HDA. However, differences in the LDA–HDA transformations in the liquids with different values of h exist at the microscopic level, which have implications in the reversibility of the LDA–HDA transformation.

In this section, we focus on the quantum systems with h = h₁, h₃. The glass behavior of the classical FJ model system is described in Refs. 31 and 34. Following Ref. 31, we prepare LDA by cooling equilibrium configurations of LDL at low pressure P₁ = −0.10 (h = h₁) and P₁ = −0.15 (h = h₃); q_T = 0.001. Figure 2 shows the potential energy (PE) and volume of the system during cooling for h = h₃. Depending on the value of n_b, some of the cooling runs exhibit crystallization; those trajectories that crystallized are shown in the insets of the figure. The main point of Fig. 2 is that the properties obtained upon cooling are very sensitive to the number of beads employed, particularly, at low temperatures. Specifically, for n_b = 50, the volume and PE of the system decrease abruptly as T → 0. This is unexpected since in the classical FJ liquid, the PE and volume are linear functions of temperature at low temperatures, within the glass domain. The unusual sharp decrease in PE and volume shown in Fig. 2 are artifacts that vanish as the number of beads per ring-polymer increases. For example, at approximately T ≥ 0.03, the PE and volume of the system converge (i.e., they are n_b-independent) for n_b ≥ 50. Instead, for approximately T ≥ 0.01, the PIMD simulations converge for n_b ≥ 150. That an increasing number of beads n_b is needed to obtain converged results with decreasing temperatures is not unexpected.^63,64 However, it is interesting that as n_b increases, the behaviors of the volume and PE with temperature become normal at all temperatures, i.e., as in the classical systems. Specifically, Fig. 2 indicates that, at low temperatures (i.e., in the glass state), both the PE and volume depend linearly with temperature for sufficiently large n_b.

The average radius of gyration of the ring-polymers R_g during the cooling runs for h = h₃ is shown in Fig. 3 (in the equilibrium domain, R_g is identical to the quantum uncertainty in the position of the quantum liquid atoms³⁵). It is expected that R_g increases upon cooling since k_sp ∝ T², i.e., the springs of the ring-polymers become weaker with decreasing temperature. However, for small values of n_b, R_g decreases anomalously at very low temperatures. Again, this anomalous ring-polymer contraction, analogous to the behavior of the PE and volume shown in Fig. 2, is an artifact that vanishes with increasing n_b. Indeed, Fig. 3 shows that for n_b = 250, R_g(T) increases upon cooling. We note that, qualitatively, similar results to those shown in Figs. 2 and 3 hold for the case of h = h₁. In this case, however, the volume, PE, and R_g(T) already converge for n_b = 50; see Fig. S5 of the supplementary material.

In Figs. 2 and 3, we distinguished among those trajectories that crystallized upon cooling from those that lead to an amorphous solid state (LDA). To do this, we follow the same criteria used in Ref. 31 to identify crystallization. Specifically, we evaluate the (global) order parameter Q₆ introduced by Steinhardt et al. in Ref. 65 (using the centroids coordinates). A configuration with Q₆ ≥ 0.20 is classified as crystalline; otherwise, it corresponds to an amorphous structure. The behavior of Q₆(T) for the cooling runs performed for h = h₃ and n_b = 250 is shown in Fig. 4(a). From these ten runs, only two trajectories show crystallization. The radial distribution functions (RDFs) of the system at the final temperature T = 0.004 are shown in Fig. 4(b). When crystallization is absent, the RDFs of the system show smooth and wide extrema that become weaker at large separations, consistent with the absence of long-range order. Instead, the RDF of the system that shows (partial) crystallization (Q₆ ≥ 0.20) exhibits distinct peaks (the corresponding RDF maxima become larger, the more crystal-like the system becomes).

In this section, we describe the pressure-induced LDA–HDA transformations at T₁(h). As it will be shown, the results are sensitive to the quantumness of the system, as quantified by Planck’s constant h. Hence, we will focus on the following two cases: (a) h = h₁ and n_b = 50 and (b) h = h₃ and n_b = 250.

(a) h = h₁. Figure 5(a) shows the volume as a function of pressure during the compression of the ring-polymer system with h = h₁ (n_b = 50) at T₁ = 0.008. The starting LDA configurations are prepared by isobaric cooling of the liquid at P₁ = −0.10, as explained in Sec. IV A. The sudden decrease in volume at P ≈ 0.75 indicates the transformation of LDA to HDA. From the nine compression runs performed, only six of them show no sign of crystallization. Specifically, for these trajectories, Q₆(P) < 0.20 for all pressures considered; see Fig. 5(b).

Included in Figs. 5(a) and 5(b) are the v(P) and Q₆(P) during the decompression-induced HDA-to-LDA transformation (red lines). The sudden increase in v(P) at P ≈ 0.1 indicates the conversion of HDA back to an LDA-like state. Further decompression of the system at P < 0.1 leads to the fracture of the LDA [at these conditions, the ring-polymer interactions cannot withhold the negative pressures (tension) applied to the system and cavitation occurs]. The results shown in Figs. 5(a) and 5(b) are qualitatively similar to those obtained in Ref. 31 for the classical FJ glass. This is consistent with the fact that h₁ is small and, hence, only mild NQE are expected. Indeed, as shown in Fig. 5(c), the ring-polymers for h = h₁ are rather localized at all pressures considered. Specifically, R_g(P) < 0.13, implying that the polymers expand for approximately $<26%$ of the hard-core radius of the classical FJ pair potential (the hard-core radius of the FJ pair potential is r = 1 in reduced units).

During compression, R_g(P) decreases smoothly, except at the LDA-to-HDA transition where R_g(P) shows a small maximum followed by a sudden decrease. During decompression, R_g(P) increases rather smoothly until it reaches the values characteristic of the starting LDA state (at P < 0.10); at the HDA-to-LDA transition, R_g(P) also shows a small maximum. Overall, for the case of h = h₁, the changes in R_g(P) during the compression and decompression runs are rather small, ∼0.08 < R_g(P) < 0.13, indicating that the spread of the ring-polymers is mildly affected during the LDA–HDA transformations (h = h₁). Figure 6 shows snapshots of ring-polymers in the LDA and HDA states for h = h₁.

Figure 7 shows the RDF of the ring-polymer system at selected volumes during the compression and decompression runs. The structural changes during the LDA–HDA transformations are similar to those observed in the classical case (h = 0). Accordingly, we refer the reader to Refs. 29, 31, and 32 for details. Briefly, the HDA state is characterized by a first peak in the RDF at r ≈ 1, corresponding to the hard-core distance of the FJ pair potential. Instead, in the LDA state, the ring-polymer centroids/beads are further apart by a distance r ≈ 1.7, close to the distance at which the FJ potential has its minimum. Not surprisingly, since R_g(P) is small, we find no differences in the RDFs of the ring-polymer centroids and beads (that belong to the same replica). Interestingly, this implies that the collapse of the LDA-to-HDA transformation in the quantum system at the h = h₁ level occurs within each replica as well. Specifically, during the LDA-to-HDA transformations, all ring-polymer beads with the same number n (n = 1, 2, ..., n_b) move closer to one another (from r = 1.7 to r = 1.0), i.e., the beads in replica n collapse as do the atoms in the classical FJ liquid.³¹ In addition, this implies that when one LDA replica transforms to HDA, all of the remaining replicas transform to HDA as well. The results of Figs. 5–7 show that in the low-quantum regime (h = h₁), the LDA before the compression and the LDA form after the decompression are practically identical. Specifically, they both have approximately the same volume, R_g, and replica RDF. Accordingly, the LDA–HDA transformation in the low-quantum regime (h = h₁) can be considered to be reversible, as found previously for the classical (h = 0) case.

(b) h = h₃. v(P), Q₆(P), and R_g(P) for the case h = h₃ (n_b = 250), during the LDA–HDA transformations at T₁ = 0.004, are included in Fig. 8. At this level of quantumness, only two of the eight trajectories show no crystallization during the compression–decompression cycle [Fig. 8(b)]. Overall, Figs. 5(a), 5(b), 8(a), and 8(b) are qualitatively similar with a sudden change in v(P) and Q₆(P) accompanying the LDA–HDA transformations. However, increasing Planck’s constant from h₁ to h₃ reduced considerably the hysteresis ΔP during the LDA–HDA transformations: ΔP ≈ 0.7, 0.2 for h = h₁ and h₃, respectively. We note that, as shown in Fig. S6 of the supplementary material, the properties we report during the LDA–HDA transformation for the case h = h₃ vary with n_b and require a large number of beads to converge (n_b > 150).

The most relevant NQE on the LDA–HDA transformations occur at the microscopic level. Specifically, as shown in Fig. 8(c), R_g(P) suddenly increases during the LDA-to-HDA transition for h = h₃ while, instead, R_g(P) suddenly decreases for h = h₁ [Fig. 5(c)]. It follows that for the more quantum system, the ring-polymers are not able to contract once the HDA state is formed. Indeed, a large increase in the ring-polymer expansion occurs during the LDA-to-HDA transformation, with R_g ≈ 0.25 in the LDA state and R_g ≈ 0.42 in the HDA state. Remarkably, during the decompression, the ring-polymers continue to expand considerably. In the final LDA state, after decompression, R_g ≈ 0.65, implying that the ring-polymers expand by ∼160% during the LDA–HDA–LDA sequence of transformations.

An extrapolation of the results obtained for the ring-polymer system to the case of the quantum liquid suggests that, during the pressure-induced LDA-to-HDA and HDA-to-LDA cycle, the quantum particles of the system become increasingly delocalized. Yet, the changes in v(P) between the final and starting LDA state are small, $<1%$⁠. It follows that while the behavior of v(P) suggest that the LDA forms before the compression and after the decompression have the same volume, the starting LDA and final LDA states are not identical microscopically, a consequence of the quantum fluctuations due to NQE. This is evident in the snapshots shown in Fig. 9. The snapshots in Figs. 9(b) and 9(d) confirm that the microscopic structures of both LDA forms are different. We note that in Ref. 6, a similar sudden atom delocalization was found during isothermal fast compressions of ⁴He.⁶⁶ 

The centroid and replica RDFs for the system with h = h₃ are shown in Fig. 10. A comparison with Fig. 7 shows that the replica RDFs are very similar for the quantum systems with h = h₁ and h₃. For example, in both cases, the first two maxima of the replica RDF for LDA are located at r = 1.7, 3.2; similarly, the first two maxima of the replica RDFs for HDA are located at r = 1.2, 1.7. This implies that at the replica level, the LDA–HDA transformations involve the same structural changes among the ring-polymer beads, independent of h (i.e., during the LDA-to-HDA transformations, the replicas collapse and the corresponding beads move closer from r = 1.7 to r = 1). The situation is different for the centroids RDF. At h = h₁, at which NQE are mild, the ring-polymer beads are very close to corresponding centroids, and hence, the centroids and replica RDFs show only minor differences (Fig. 7). Instead, at h = h₃, the ring-polymer are delocalized and, accordingly, the replica and centroid RDFs differ; see, e.g., the centroid and replica RDFs of LDA at r = 2.2, 3.9 in Fig. 10(a) and the centroids and replica RDFs of HDA at r = 1.2 in Fig. 10(b).

While the LDA–HDA transformation is reversible in the low-quantum regime (h = h₁), the reversibility of this transformation is less evident in the high-quantum regime (h = h₃). On the one hand, the LDA forms before the compression and after the decompression have practically the same volume and replica RDFs.⁶⁷ On the other hand, they exhibit different R_g (i.e., atoms delocalization). Thus, in principle, one may conclude that both LDA forms are different glasses. However, the subtlety in this interpretation is that LDA and HDA are not analogous to crystalline states, with unique RDFs and thermodynamic properties. As discussed in Ref. 16 for the case of water, LDA and HDA represent two families of glasses; glassy forms within each family exhibit similar RDF and (out-of-equilibrium) thermodynamic behavior, including the glass transition temperature. This is because glasses are out-of-equilibrium systems with properties that vary with the specific method of preparation (including cooling and compression/decompression rates). The case of water is a nice example. In Refs. 68 and 69, it is shown that the water LDA and HDA “families” are characterized by different regions of the potential energy landscape (PEL), each region corresponding to a wide set of slightly different LDA/HDA forms. Therefore, with the available data, it remains unclear whether the LDA forms obtained before the compression and after the decompression are (i) the same LDA form, (ii) two glassy forms belonging to the same LDA family of glasses, or (iii) two glasses, say, LDA_I and LDA_II. This question will be addressed in a future study.

In this section, we test the results from Secs. III and IV B by studying the transformations in LDA and HDA induced by isobaric heating at different pressures. Specifically, we show that for the ring-polymer systems with h = h₁, h₃, the glass transition temperature of LDA, $TgLDA(P)$⁠, is anomalous, while the glass transition temperature of HDA, $TgHDA(P)$⁠, is not. This is consistent with the iso-diffusivity lines of the quantum liquids discussed in Sec. III. In addition, we show that, at least for h = h₁, isobaric heating of HDA at low pressures results in LDA before the liquid state is reached; similarly, isobaric heating of LDA at high pressures results in HDA before the liquid state is reached. These LDA–HDA transformations induced by heating are consistent with the pressure-induced LDA–HDA transformation discussed in Sec. IV B.

We follow the same procedure used in Ref. 34 to study the heating-induced transformations in LDA and HDA for the classical FJ liquid. Accordingly, we refer the reader to Ref. 34 for details. Briefly, upon isobaric heating, we detect the glass-to-liquid transition temperatures by monitoring the MSD of the ring-polymer centroids as a function of temperature. T_g is defined as the temperature at which MSD(T) = 1; see Figs. S7(a) and S7(b) and Figs. S8(a) and S8(b) of the supplementary material. The glass–glass transition temperature can be detected by monitoring the volume of the system upon heating, v(T). During the LDA-to-HDA and HDA-to-LDA transitions, v(T) exhibits a sudden change, which is used to identify the corresponding transition temperature; see Figs. S7(c) and S7(d) of the supplementary material. Overall, our results for the QFJ model with h = h₁, h₃ are consistent with those reported in Ref. 34 for h = h₀ = 0; see also Fig. S2(e) of the supplementary material. Accordingly, we refer the reader to Ref. 34 for details.

Figure 11 shows the transition temperatures induced upon heating compressed LDA (blue triangles) and decompressed HDA (red triangles) for the ring-polymer system with Planck’s constant h = h₁. The compressed LDA glasses are obtained during the compression runs indicated by the blue lines in Fig. 5(a), at approximately P < 0.7. Depending on the pressure considered, heating LDA leads to either (i) the LDA-to-LDL transformation or (ii) the LDA-to-HDA-to-HDL sequence of transformations. Case (i) occurs at low pressures; the LDA-to-LDL transformation defines $TgLDA(P)$ (blue solid triangles). Case (ii) occurs at higher pressures. The LDA-to-HDA transformation temperatures are indicated by the blue empty triangles in Fig. 11; the HDA-to-HDL transition is observed upon further heating, but it is not indicated in Fig. 11 for clarity. As discussed in Ref. 34 for the case h = h₀, the LDA-to-LDL and LDA-to-HDA transformation temperatures define a single path in the P–T plane (blue dashed line), which indicates the limit of “stability” of LDA relative to LDL and HDA upon isobaric heating.

A similar set of transformations occur upon heating decompressed HDA. The HDA forms are obtained during the decompression runs indicated by red lines in Fig. 5(a) at approximately P > 0.10. Depending on the pressure considered, heating HDA leads to either (i) the HDA-to-HDL transformation or (ii) the HDA-to-LDA-to-LDL sequence of transformations. Case (i) occurs at high pressures; the HDA-to-HDL transformation defines $TgHDA(P)$ (red solid triangles). Case (ii) occurs at lower pressures. The HDA-to-LDA transformation temperatures are indicated by the red empty triangles in Fig. 11; the LDA-to-LDL transition is observed upon further heating, but it is not indicated in Fig. 11 for clarity. As discussed in Ref. 34 for the case h = h₀, the HDA-to-HDL and HDA-to-LDA transformation temperatures define a single path in the P–T plane (red dashed line), which indicates the limit of “stability” of HDA relative to HDL and LDA upon isobaric heating.

The heating-induced transformation for the ring-polymer systems with Planck’s constant h = h₃ is indicated in Fig. 12. At this level of quantumness, we were not able to identify the LDA-to-HDA and HDA-to-LDA transformation temperatures. Accordingly, there are no empty red and blue triangles in Fig. 12, and only $TgLDA(P)$ and $TgHDA(P)$ are indicated (solid blue and red triangles, respectively). The difficulty in detecting the glass–glass transformations is due to the small hysteresis in the pressure-induced LDA–HDA transformations shown in Fig. 8(a). Such a small hysteresis implies that there is a very small pressure interval for which isobaric heating can bring the system from LDA (HDA) to HDA (LDA). In other words, for the case h = h₃, the liquid spinodals and glass–glass transformation lines are very close to one another.

We performed PIMD simulations of the FJ liquid using different values of Planck’s constant h and over a wide range of volumes and temperatures. Tuning h allowed us to study the behavior of the FJ liquid as NQE (i.e., quantum fluctuations due to atoms delocalization) are introduced, including the classical (h = 0) and quantum liquid cases (h > 0).

Our PIMD simulations in the equilibrium liquid state, combined with the RPMD method, show that introducing NQE increases the mobility of the atoms, suppressing the glass state. Specifically, the iso-diffusivity lines move to lower temperatures as h increases and the atoms become more delocalized (Fig. 1). In addition, we found that the maximum diffusivity line, an anomalous water-like property exhibited by the FJ liquid, also shifts toward lower temperatures as NQE are introduced. Overall, the observed NQE on the dynamics of the quantum FJ liquid are consistent with the NQE on the thermodynamic properties of the FJ liquid reported in Refs. 28 and 29 based on PIMC simulations. In those works, it was shown that introducing NQE shifts the LLCP and anomalous properties, such as the density maxima and compressibility maxima lines, to lower temperatures. Interestingly, our results suggest that with increasing h, the LLCP temperature → 0 while, simultaneously, the vitrification temperature is suppressed. This means that, for large h, it may be possible to have a quantum liquid with anomalous water-like dynamic and thermodynamic properties (e.g., density and diffusivity maxima lines) and, yet, with no accessible LLCP.

We note that all the classical/quantum FJ liquids considered here exhibit an anomalous iso-diffusivity lines (Fig. 1). Specifically, at low pressures, within the LDL domain, the iso-diffusivity lines have an (anomalous) negative slope (Fig. 1). In this case, isothermal compression enhanced the dynamics of the liquid atoms. Instead, at high pressures, in the HDL domain, the iso-diffusivity lines are positively sloped meaning that, as one would expect, isothermal compression decreases the motion of atoms. Interestingly, the NQE are more pronounced at high pressure in the HDL state. Increasing h shifts considerably the iso-diffusivity line in the HDL region toward lower temperatures, while in the case of the LDL domain, the iso-diffusivity line barely shifts with temperature.

The properties of the quantum FJ liquids are intimately connected to the behavior of the ring-polymer system in the glass state. Accordingly, we extended our PIMD/RPMD simulations to the glass domain. We found that all ring-polymer systems (h ≥ 0) exhibit glass polymorphism, such as water. For all values of h studied, compression of LDA at low temperatures produces HDA; the subsequent decompression of the HDA produces the LDA-like state. In all cases studied, the LDA–HDA transformation is accompanied by sudden changes in the volume of the system, consistent with the results of Refs. 31 and 34 for the classical (h = 0) case. The effects of introducing NQE (increasing h) is to reduce the hysteresis of the LDA–HDA transformation.

Microscopically, the LDA–HDA transformations differ profoundly in the classical and quantum cases. In the classical limit (low h), the ring-polymers become more localized during the LDA-to-HDA transformation and return to their original size during the decompression-induced HDA-to-LDA transformation. Instead, at large h, the ring-polymers expand during both the LDA-to-HDA transformation and the subsequent HDA-to-LDA transformation. In the case of h = h₃, R_g increases by 160% during the LDA-to-HDA-to-LDA cycle. Remarkably, there is no net volume change during this cycle. This has important implications in the reversibility of the LDA transformation. We found that in the classical limit, the LDA–HDA transformation is indeed reversible since the starting LDA (before the compression) and final LDA (after the decompression) have the same volume, R_g, and replica RDFs. Instead, in the large-quantum regime, the difference in the localization of the atoms in the starting and final LDA form makes it unclear to conclude whether these two glassy forms belong to the same family of LDA. The reversibility of the LDA–HDA transformation at large h-values will be studied in a separate work.

In the last section of this work, we discuss the NQE on the glass transition temperatures of LDA and HDA. As for the classical case,^31,34 the glass transition of the LDA is anomalous and decreases with increasing pressure, while, instead, the glass transition temperature of HDA increases with increasing temperature. This finding is fully consistent with the iso-diffusivity lines of LDL and HDL calculated from the equilibrium RPMD simulations. The role of introducing NQE is to shift the glass transitions of LDA and HDA to lower temperatures.

The results presented in this work are important to understand polyamorphism in light elements, such as hydrogen, for which different studies indicate the existence of an LLPT/LLCP. For example, glass polymorphism in hydrogen may be accompanied by an increasing delocalization of the atoms in the system. In such cases, the atoms in the final LDA state (of the pressure-induced LDA-to-HDA-to-LDA cycle) would be much more delocalized than in the starting LDA state. Our results may also shed light on isotope substitution effects in supercooled and glassy water.

See the supplementary material for the definition of the FJ pair potential; the effects of γ, Δ, and n_b on the PIMD/RPMD simulations; and the thermodynamic and dynamical properties of the FJ liquid with Planck’s constant h = h₁.

This work was supported by the SCORE Program of the National Institutes of Health under Award No. 1SC3GM139673 and the NSF CREST Center for Interface Design and Engineered Assembly of Low Dimensional systems (IDEALS) (NSF Grant No. HRD-1547380). This work was also supported, in part, by a grant of computer time from the City University of New York High Performance Computing Center under NSF Grant Nos. CNS-0855217, CNS-0958379, and ALI-1126113.

The authors have no conflicts to disclose.

The data that support the findings of this study are available within the article and its supplementary material.