Early prediction of spinodal-like relaxation events in supercooled liquid water

Bulk liquid water presents many anomalies, such as the apparent divergence of the specific heat (C_p) and the isothermal compressibility factor (κ_T) as the temperature decreases. Many hypotheses have been suggested in order to justify these phenomena, and among them, the most credited one is the existence of a first-order liquid–liquid phase transition (LLPT).^1,2 This scenario relies on the existence of a liquid–liquid critical point (LLCP) localized at high pressure in the deeply supercooled region of the pressure–temperature diagram (the so-called No Man’s Land).^3–5 This point marks the end of the coexistence line between two different phases known as low-density liquid (LDL) and high-density liquid (HDL).² Besides differing in density, the LDL and HDL states also have dissimilar structural, thermodynamic, and dynamic properties.^6–13 The short-range LDL structure is comparable to that of ice (I_h),¹⁴ with four nearest neighbors bonded via hydrogen bond (HB) in a regular tetrahedral arrangement; however, differently from ice I_h, it lacks long-range order. The HDL structure is characterized by a distorted HB pattern with a fifth non-hydrogen bonded nearest neighbor (the so-called interstitial water molecule). Despite these distinctions, the existence of the LLCP in bulk water has not yet been proven experimentally due to the fast crystallization of water upon cooling, and thus, the LLPT scenario is mainly supported by computational studies and experiments on confined water and amorphous ice.

Many computational works provided compelling evidence of the existence of the LLPT, and associated LLCP, in several classical and ab initio-based water models and thoroughly studied the thermodynamic, dynamic, and structural properties of the two liquid phases differing in density.^15–20 In general, studies conducted on different water models in different conditions point to the result that, microscopically, the LLPT is associated with a rearrangement of the HB network.^3,15 This evidence is further supported by the fact that this behavior has been observed experimentally in other systems that exhibit a liquid–liquid transition, such as cerium²¹ and melts of yttrium-alumina,²² in which a rearrangement of the bond network is also observed.

Despite the vast amount of work on the topic, the microscopic interconversion mechanism between HDL and LDL is still unclear. Given the fact that both liquid phases are metastable with respect to ice-Ih, but stable with respect to each other in the supercooled region near the coexistence line, it is plausible that the mechanism that may take place in the LLPT is the nucleation, and subsequent growth, of the LDL from the HDL phase and vice versa.³ 

According to Classical Nucleation Theory (CNT),²³ nucleation processes can be homogeneous or heterogeneous and take place in classical or non-classical pathways. In homogeneous nucleation, a group of molecules comes together to form a critical nucleus. This is mainly due to density or thermal fluctuations that allow the overcoming of the free energy formation barrier. In contrast, heterogeneous nucleation processes involve the formation of a new phase due to the presence of nucleation sites or impurities. In computer simulations of pure water, which use periodic boundary conditions (i.e., no walls are present), only homogeneous nucleation can take place. Furthermore, the nucleation time of a stable state from a metastable one is inversely proportional to the volume of the system. Given that the experimental volumes are much larger when compared to computer simulation volumes, the latter exhibit longer nucleation times. This allows the time scales of the LLPT and ice nucleation to be better separated and characterized.³ On the other hand, this also implies that the study of the transition can be very computationally demanding, especially when equilibrium states are involved. Moreover, the transition kinetics widely vary depending on the T and P simulation conditions. As of now, the nucleation process of LDL into HDL and vice versa has not been fully characterized due to these contributing factors.

Conversely, a transition process is found to be much faster when there is no free energy barrier to overcome; as in the case of a spinodal decomposition, a mechanism in which an unstable phase spontaneously separates into two phases, ultimately leading to the formation of the stable phase.^24,25 Unlike nucleation, in which the stable phase nuclei form in a discrete number of points in the volume, spinodal decomposition occurs throughout the whole volume of the system, and the stable phase is organized in fluctuating microdomains that immediately start growing. The growth of the domains takes place via diffusion and coalescence until they become large enough to become stable.^26,27 Therefore, from a transition kinetics point of view, both spinodal decomposition and nucleation mechanisms are followed by a similar growth phase.

Given the wide range and the complexity of the possible pathways that can take place, many models have been developed in order to explain the kinetics of growth processes that occur in nature. Among them, we can find the sigmoidal models used to describe “cooperative” phenomena.^28–32 In most cases, the sigmoidal curve is made up by a lag phase (often reported as induction period), a growth phase, and a final plateau phase. All of these models require a deep knowledge of the properties of the species involved in the growth event. In the LLPT scenario, this implies an unambiguous distinction between both phases. However, even from a computational point of view, it is not a trivial task to determine if a molecule belongs to the LDL or HDL phase, and, as a result, many structural order parameters have been developed with the purpose of characterizing simulated water at different temperatures and pressure conditions. Among them, we can find the local structure index (LSI),³³ d₅,³⁴ V₄,^35, ψ,^36, ζ,³⁷ and the node total communicability (NTC) we recently proposed.^38,39 In our previous works, we showed that the NTC is able to identify the two liquid states differing in density both along an isobar that crosses the liquid–liquid coexistence line³⁸ and at ambient pressure where both states coexist.³⁹ We also showed that the HDL-like state is characterized by the presence of high connectivity patches, i.e., highly connected networks composed of nodes with an above-average connectivity and an increased number of interstitial water molecules. We observed the presence of small highly connected patches also at low temperatures, where the prevailing state is LDL-like, and we speculated that these small patches might function as initial sites for the formation of extended HDL regions. In this work, we expand on this hypothesis by investigating the liquid–liquid phase transition using the NTC to distinguish the two phases along the transition. Using as a starting point LDL water, we perform molecular dynamics (MD) simulations at a temperature at which the equilibrium state is HDL-like, thus investigating how the HDL-like phase is formed from the LDL phase. We propose a “spinodal-like” scenario and investigate the role of highly connected HDL-like patches in the formation of the growing phase.

We use a 300 ns-long MD simulation of 710 water molecules at 1950 bar at 170 K to extract five LDL conformations. These five structures are then used as starting points for five independent 15 ns-long simulations, performed at the same pressure (1950 bar) and at 200 K with the TIP4P/2005 water model⁴ that exhibits a metastable liquid–liquid critical point under deeply supercooled conditions.⁹ In the discussion, we will refer to the five trajectories as T1-2-3-4-5. We then use the same protocol for the investigation of a four-times bigger box. We perform a 250 ns-long MD simulation of 2840 water molecules at 1950 bar at 170 K to extract two LDL conformations. These two structures are then used as starting points for two independent 15 ns-long simulations, performed at the same pressure (1950 bar) and at 200 K. In the discussion, we will refer to these two trajectories as T1B and T2B.

All MD simulations are performed in the NPT ensemble with the 5.1.2 version of the GROMACS software⁴⁰ using a cubic simulation box. Temperature and pressure are kept constant by using the Nosé-Hoover temperature coupling⁴¹ and the Parrinello–Rahman barostat with 2 ps relaxation times.⁴² Periodic boundary conditions are used, long-range electrostatic interactions are treated with the particle mesh Ewald method⁴³ with a real space cutoff of 0.9 nm, and for short-range interactions, a cut-off radius of 0.9 nm is employed. All bonds are constrained using the LINCS algorithm⁴⁴ along with a 2 fs time step.

In this section, we recall some basic concepts of graph theory,⁴⁵ we explain how we construct the graphs given an MD trajectory, and we give a definition of patches. We refer to our previous works^38,39,46 for a more detailed description.

From each frame of an MD trajectory, we construct the associated graph in the following way: we consider the oxygen atoms as the nodes of our graphs, and two oxygen atoms are connected if their distance is $≤0.35$ nm; see the work of Faccio et al.³⁸ for explanations on the choice of this threshold. To calculate the distance between the oxygen atoms, we use the minimum image convention. In this way, we introduce the periodic boundary conditions also in our graphs.

As explained in the previous works,^38,39,46 we define a molecule in the LDL/LDL-like phase if its NTC value is $≤95$⁠; otherwise, it is in the HDL/HDL-like phase. In this work, we also define two types of patches: HDL-like patches and high connectivity (HC) patches. HDL-like patches are connected graphs composed of nodes in the HDL-like phase (i.e., the NTC of which is $>95$⁠). HC patches are connected graphs composed of nodes that, besides being in the HDL-like phase, have a number of nearest neighbors (within 0.35 nm) greater than 5. Since 5 is the typical number of nearest neighbors in the HDL-like states (i.e., one interstitial water molecule), these patches are composed of molecules with an increased number of interstitial water molecules.

The average lifetime of the biggest patch in a time interval [t, t + T] is computed in this way: let B be the biggest patch at time t; at time t + 1, we evaluate the intersection of B with all the existing patches; if at least one of these intersections has cardinality equal or greater than a given value μ (here, we choose μ as 15% of the cardinality of B), the lifetime of B increases by 1; the algorithm then iterates the procedure intersecting B with the patches formed at subsequent time frames until, at a given frame, all the intersections between B and the patches have cardinality $<μ$⁠; when this occurs, the biggest patch is considered vanished. To start with a new experiment, we jump 75 ps to have statistically uncorrelated experiments and restart the procedure. At the end, we consider the average of the obtained lifetimes.

We note that $d2f(t)dt2=0$ ⇔ $(r2e−r(t−τc)−r2)=0$⁠, i.e., t = τ_c (we suppose 0 ≤ a < 1 and r ≠ 0). Furthermore if t < τ_c ⇒ $e−r(t−τc)>1⇒d2f(t)dt2>0$⁠; otherwise, if $t>τc⇒0<e−r(t−τc)<1⇒d2f(t)dt2<0$⁠. Accordingly, τ_c is again an inflection point and it is the time of maximum rate, and $f(τc)=1+a2$⁠.

Starting from five configurations extracted from an MD simulation of 710 TIP4P/2005 water molecules at 170 K and 1950 bar, we simulate a temperature jump to 200 K maintaining the same pressure. The second critical point for TIP4P/2005 has been located at different pressures and temperatures, with a recent accurate estimate at 1860 bar and 172 K.¹⁸ Therefore, the initial simulation conditions (i.e., 170 K and 1950 bar) are very close to the coexistence line, as recently shown.⁵¹ This proximity, combined with the well-known issues related to long relaxation times at low temperatures, does not allow us to state that the extracted structures are fully equilibrated. However, for the scope of this work, we require the starting points for our temperature jump simulations to be representative of LDL structures. To ensure this, we benchmark the density of the extracted structures (Table I) against that of a previous MD simulation under the same conditions,⁵¹ in which LDL configurations are retained for several microseconds (977 kg/m³). Thus, we start our non-equilibrium MD simulations from an LDL-water state (extracted at 170 K) and induce a transition to an HDL-like state (at 200 K) (see Fig. 1). It has to be noted that our simulation protocol that uses a temperature jump determines an abrupt transition to a region past the LDL spinodal line obtained from the Potential Energy Landscape of the TIP4P/2005 water model.^51,52 Therefore, we expect along our simulations a behavior compatible with a spinodal-like decomposition scenario in which the instability of the LDL phase leads to an HDL-like composition.

Along the five trajectories, we calculate the NTC value for each molecule (see Sec. II B) and use it to calculate the fraction of HDL-like molecules. As done in our previous works, for NTC ≤95, a water molecule is assigned to the LDL state, and for NTC $>$ 95 to the HDL-like state.^38,39 In the five starting configurations, the HDL fraction is ≈0.1 (Table I), matching the average HDL fraction calculated in the 170 K MD simulation.³⁸ Consistently, the NTC of the initial configurations averaged over the molecules of the box matches the average NTC of the MD simulation at 170 K, i.e., 71. Along time, we observe an increase in density, showing the transition from the LDL to the HDL-like state, and a corresponding increase in the average NTC [Figs. 2(a) and 2(b)]. As already observed along a simulation of near-critical supercooled water,³⁶ the time evolution of NTC and density are highly correlated, showing that the NTC well catches the LDL to HDL-like transition. Consistently, the HDL-like fraction increases along time with an approximately sigmoidal trend [Fig. 2(c), black].

To better describe the growth of the HDL-like population, we apply a logistic model to fit the HDL-like fraction [Fig. 2(c), red]. Logistic curves are commonly used to describe restrained population growth, i.e., reaching an asymptotic position. As outlined in Sec. II C, the logistic fit provides as an output parameter the time τ_c of the inflection point of the curve, corresponding to the maximum growth rate. In addition, the calculation of the second derivative of the curve provides the times τ₁ and τ₂ corresponding to the maximum and minimum growth accelerations, respectively. For the five trajectories, τ_c, τ₁, and τ₂ are reported in Table I. More details on the fits performed can be found in the supplementary material in Table S1 and Figs. S1–S5.

To investigate the formation of the HDL-like state, we are mainly interested in the events occurring in the first stage of the process, i.e., before τ_c. Within this time interval, two distinct kinetic regimes can be identified, with τ₁ marking the switch between them. In what follows, we investigate the behavior of the system in these two kinetic regimes, with a special focus on the early stages of the process (before τ₁).

In a spinodal-like decomposition scenario, the early formation of delocalized HDL-like fluctuating domains is expected. Along time, the size of these fluctuating domains grows until they are large enough to become stable. Subsequent coarsening processes occur at a slower rate eventually leading to the pure HDL-like phase. In this view, we inspect the presence of HDL-like domains, and their behavior, at the earliest time of our MD simulations. Within our graph-based framework, we define as domains/patches connected networks of HDL-like molecules, i.e., with NTC $>$ 95 (see Sec. II B). Given that, according to the NTC, the HDL-like fraction along the 170 K MD simulation is ≈0.1, we started by checking the existence of any HDL-like patch at low temperature. We obtain that almost all (83%) HDL-like molecules present in the 170 K simulation belong to patches made up of at least 5 HDL-like molecules. On average, 3.1 patches made up of ≈19 molecules are present at each MD frame. The most probable dimension of the biggest patch is ≈15–25 molecules and that of the second biggest patch is ≈10 molecules, as shown in Fig. S6 in the supplementary material.

The HDL-like patches are then investigated along the 5 MD simulations after the temperature jump. Figures 3(a) and 3(b) show that, along time, there is an enlargement of the HDL-like patches and a corresponding reduction in the number of patches. Both trends show that after τ₂, there is a single patch that almost coincides with the entire simulation box, i.e., almost the whole system is HDL-like. In addition, Fig. 3(c) shows that these patches are highly fluctuating: by monitoring the residence time of the molecules in the HDL-like phase, it can be observed that, on average, an initially HDL-like molecule becomes LDL-like in a few ps. Figure 3(c) (black curve) is obtained as follows: at t = 0, HDL-like molecules are identified using the NTC; the fraction of these molecules that are still HDL-like is monitored from t = 1 to t = 25 ps; the same calculation is then repeated taking as initial reference time multiples of 100 ps (i.e., 75 ps are skipped between two consecutive time intervals) until t = τ₁; the trends obtained are then averaged to produce the final curve. The same procedure is repeated for the time intervals τ₁ − τ_c and τ_c − τ₂ (red and green curves). In Fig. 3, the data obtained for T4 are reported; those for T1–T3 and T5 are reported in Figs. S7–S10 in the supplementary material. The results show that there is a ps time scale exchange of molecules between the HDL-like and LDL-like phases, which starts to slow down after τ_c. For comparison, in Fig. S11 of the supplementary material, we also report the HDL-like fraction, patch characterization, and residence time of the molecules in the HDL-like phase along the 170 K MD. The plots show a slightly faster exchange of molecules coupled with the presence of small HDL-like patches. Therefore, without any temperature jump, a similar residence time is not associated with patch growth.

The data in Fig. 3, taken together, suggest a scenario in which transient fluctuating HDL-like domains emerge in the LDL-like unstable phase; these domains appear to grow in size through coalescence. All these features are consistent with the expected spinodal-like decomposition scenario, within which it is anticipated that, at a certain point, the fluctuating domains become large enough to also become stable, initiating the growth phase. Therefore, we focus on the behavior of connected patches including at least 25 molecules (i.e., the same order of magnitude of the biggest patch in the 170 K MD simulation that corresponds, at 200 K and 1950 bar, to approximately two hydration shells). By monitoring the time evolution of the ratio between the number of HDL-like molecules inside patches of at least 25 molecules and the total number of HDL-like molecules [Fig. 4(a)], it can be observed that this fraction increases until it reaches a plateau. It can also be observed that t = τ₁ corresponds, for all trajectories, to a fraction above 0.95, suggesting that the kinetic crossover occurs when almost all HDL-like molecules are organized in sufficiently big patches. In addition, the fraction of HDL-like molecules is ≈0.30–0.35 at t = τ₁ for all trajectories [Fig. 4(b)]. In Fig. 4, T3 is not fully in line with the above conclusions, reaching a fraction of 0.95 before τ₁ [panel (a)] and showing a larger HDL-like fraction at t = τ₁ [panel (b)]. We believe that these slight discrepancies can be due to the lower quality of the fit used to estimate τ₁ in T3, see Fig. S3 in the supplementary material.

Interestingly, a fraction of ≈0.2–0.3 of the growing phase was previously identified, using a series of equilibrium MD simulations at different temperature and pressure conditions, as a critical concentration inducing changes in the dynamical and thermodynamic behavior of supercooled water.⁵³ At this concentration, the fragmentation of the dominating phase was observed along with the formation of connected patches of the growing phase. Our data support this observation and suggest that this critical concentration also corresponds to a kinetic crossover along the phase transition. Furthermore, at the same kinetic crossover, the fluctuating patches that form at the initial stages of the simulation have become stable. As a matter of fact, for t $<τ1$⁠, the patches’ mean lifetime is ≈60 ps, while the biggest patch present at t = τ₁ is stable until the end of the simulation (see Sec. II B for details on the patches mean lifetime calculation). Therefore, after τ₁, the HDL-like fraction growth is governed by the coalescence and coarsening among stable HDL-like domains.

While the above data clearly show that HDL-like patches are highly fluctuating before τ₁ [Fig. 3(c)], visual inspection of the trajectories shows that these domains tend to appear preferentially in specific regions of the simulation box in all trajectories. This observation suggests that the formation of HDL-like domains in the earliest stages of the process could have an effect on the dynamic evolution of the transition. To substantiate this hypothesis, we focus on the behavior of the system in the first 50 ps of the five MD trajectories. In particular, we monitor the time evolution of the ten molecules that have the highest NTC value in the first 50 ps of each trajectory. As shown in Fig. 5, the ten molecules that exhibit the highest NTC in the first 50 ps of the simulation display a higher NTC with respect to ten randomly extracted molecules also at much longer times. This points to the existence of regions in which the formation of fluctuating HDL-like domains is favored.

In our previous works, we observed the presence of small highly connected HDL-like domains in the LDL-like phase.³⁹ These highly connected (HC) patches differ from the above defined HDL-like patches in their increased connectivity with respect to typical HDL-like networks. In other words, these patches are composed by molecules that are not only in the HDL-like phase as defined by their NTC value, but that also feature an increased number of interstitial water molecules (see Sec. II B). By checking the formation of HC patches at the very earliest times of our MD simulations, we note that the origin of the regions in which the formation of HDL-like domains is favored seems to be related to these HC patches. As a matter of fact, we observe that, on average, 75% of the ten molecules that show the highest NTC in the first 50 ps and that show an above average NTC along the simulation belong to 1–2 HC patches that form at 5–15 ps. As shown in the representative snapshots in Fig. 5, HDL-like patches preferentially form around these HC clusters. In line with what is observed in a binary fluid, we suggest that early forming HC domains can act as “defects” around which the formation of the growing phase is favored.⁵⁴ This can be appreciated in a more quantitative fashion in Fig. 6. In that figure, panel (b), molecules are colored differently according to their LDL- or HDL-like state and ordered according to their distance from the molecule displaying the highest NTC in the first 50 ps of T4 (similar plots for T1–T3 and are reported in Figs. S12–S15 of the supplementary material). It can be observed that as the HDL-like fraction increases [panel (a)], the region around the molecule with the early highest NTC stays persistently in the HDL-like phase and that the HDL-like region around such molecule becomes larger, although the formation of less stable HDL-like regions can also be observed. This is also shown in panel (c), in which the probability of each molecule to be in the HDL-like state is reported. The figure shows that the region surrounding the highest NTC molecule has an increased probability of populating the HDL-like state. It was previously shown that regions characterized by “structural defectiveness,” i.e., by an HDL-like distorted HB network, can act as early time predictors of intermittent glassy relaxation events in supercooled water.^55,56

To investigate if and how the obtained results are sensitive to the box dimension, we repeat the above analyses on a four times larger box (2840 water molecules). We extracted two LDL conformations from an MD simulation of 2840 water molecules at 1950 bar at 170 K and used these configurations as starting points for two independent 15 ns-long simulations, performed at the same pressure (1950 bar) and at 200 K. The results show that the main features observed in the smaller box (710 molecules) are similar to those observed in the bigger box (2840 molecules). More specifically, we observe that the kinetic crossover as obtained from the logistic fit (Figs. S16 and S17 in the supplementary material) occurs at a similar HDL fraction (0.35–0.40) and that, as in the smaller box, at the kinetic crossover all HDL molecules are organized in patches with a dimension equal or bigger than the average HDL patches dimension in the MD simulation at 170 K (Fig. 7). In addition, we also observe that in the early stages of the temperature jump simulations, the mean residence time of a molecule in the HDL-like state is low, supporting the spinodal-like decomposition scenario observed in the smaller box (Figs. S18 and S19 in the supplementary material). Most importantly, we observe, in the bigger box, the preferential growth around highly connected patches that are identified with the NTC at the earliest stages of the temperature jump simulations (Figs. S20 and S21 in the supplementary material). The parameters of the logistic fit and details on the starting structures are reported in Tables S2 and S3 of the supplementary material. The above results show that early forming HC patches, as identified through the NTC, foster the development of HDL-like domains along the transition. Notably, these outcomes also confer the NTC the capability of being a predictive indicator for later relaxation events.

Using the NTC as an order parameter, we delved into the investigation of the density relaxation process of supercooled liquid TIP4P/2005 water exposed to an instantaneous temperature jump from 170 to 200 K at 1950 bar. By inspecting the evolution of the composition of the system and the HDL-like patches, we find the following:

• The relaxation process is compatible with a spinodal-like decomposition scenario schematized in Fig. 8, in which, in the initial stage, the domains are in a fluctuating regime and the lifetime of the patches is of the order of a few tens of picoseconds, while at later stages, the coarsening process occurs and the domains become stable.

• When the HDL-like fraction reaches $≈0.30−0.35$ and the majority of the HDL-like molecules are found in patches that are large enough to become stable, the system undergoes a kinetic crossover. A fraction of $≈0.2−0.3$ of the growing phase was previously identified as a “critical” concentration inducing changes in the dynamical and thermodynamic behaviors of supercooled water.⁵³ 

• There are regions where the formation of extended and stable patches is more favored and these regions can be predicted using the NTC. As a matter of fact, HDL-like patches preferentially grow around highly connected patches that form in the first stages (5–15 ps) of the relaxation process.

Since the transition kinetics can drastically vary depending on the simulation conditions (i.e., pressure and temperature),³ we perform a preliminary test on some configurations extracted from a ≈40 μs-long MD simulation of the TIP4P-2005 water model at 177 K and 1750 bar.¹⁸ We observe that, at the kinetic crossover (i.e., at τ₁), the HDL-like fraction is close to ≈0.3 (precisely, between 0.25 and 0.3) also for three of the spontaneous transitions occurring in proximity to the critical point. We also observe that along the MD simulation at HDL fractions above 0.25–0.3 all HDL molecules are arranged in connected patches including at least 25 water molecules, as observed along the temperature jump simulations (Figs. S22 and S23 in the supplementary material). We acknowledge that this preliminary test does not imply that the spinodal-like decomposition scenario we observe along our temperature jump simulations remains unaltered in the vicinity of the second critical point. Nonetheless, it is worth highlighting the efficacy of the NTC as an order parameter in unraveling the complex dynamics of the relaxation process, offering valuable insights into the behavior of supercooled liquid water.

The supplementary material includes logistic fits for trajectories T1–T5, parameters a, τc, and r obtained by curve-fitting simulated data to the logistic function, characterization of the patches in the 170 K MD simulation, residence time of the molecules in the HDL-like phase for trajectories T1–T3 and T5, patches characterization and residence time for the 170 K MD simulation of the 710 water molecules box, figures analogous to Fig. 6 for trajectories T1–T3 and T5, density of the initial structure, ρ0, logistic fits for trajectories T1B and T2B, parameters a, τc, and r obtained by curve-fitting simulated data to the logistic function, fraction of HDL-like molecules at time 0, average NTC value at time 0, τ1 and τ2 for trajectories T1B and T2B, residence time of the molecules in the HDL-like phase for trajectories T1B and T2B, figures analogous to Fig. 6 for trajectories T1B and T2B, and time evolution of the HDL fraction and patches investigation for the trajectory in proximity to the critical point.

The authors acknowledged Francesco Sciortino for providing the configurations along the MD simulation in proximity to the critical point. They acknowledged the CINECA Award under the ISCRA Initiative for the availability of high-performance computing resources and support. This work was supported by the European Union - NextGenerationEU under the Italian Ministry of University and Research (MUR), PRIN2022-P2022MC742PNRR, CUP E53D23018440001.

The authors have no conflicts to disclose.

Nico Di Fonte: Data curation (equal); Formal analysis (equal); Investigation (equal); Software (equal); Writing – original draft (equal). Chiara Faccio: Data curation (equal); Formal analysis (equal); Investigation (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Laura Zanetti-Polzi: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Isabella Daidone: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article or supplementary material. Additional data are available from the corresponding author upon reasonable request.