Nonmonotonic dependence of the absolute entropy on temperature in supercooled Stillinger-Weber silicon

Using a recently developed thermodynamic integration method, we compute the precise values of the excess Gibbs free energy (G^e) of the high density liquid (HDL) phase with respect to the crystalline phase at different temperatures (T) in the supercooled region of the Stillinger-Weber (SW) silicon [F. H. Stillinger and T. A. Weber, Phys. Rev. B. 32, 5262 (1985)]. Based on the slope of G^e with respect to T, we find that the absolute entropy of the HDL phase increases as its enthalpy changes from the equilibrium value at T \ge 1065 K to the value corresponding to a non-equilibrium state at 1060 K. We find that the volume distribution in the equilibrium HDL phases become progressively broader as the temperature is reduced to 1060 K, exhibiting van-der-Waals (VDW) loop in the pressure-volume curves. Our results provides insight into the thermodynamic cause of the transition from the HDL phase to the low density phases in SW silicon, observed in earlier studies near 1060 K at zero pressure.

💡 Research Summary

In this work the authors investigate the thermodynamic behavior of supercooled silicon modeled with the Stillinger‑Weber (SW) potential, focusing on the high‑density liquid (HDL) phase and its transition to low‑density phases near 1060 K at zero pressure. Using a recently developed thermodynamic integration (TI) scheme combined with NPT Monte‑Carlo simulations, they compute the excess Gibbs free energy G⁽ᵉ⁾ = G_HDL – G_crystal for a system of 512 atoms over the temperature range 1060 K–1070 K. By differentiating G⁽ᵉ⁾ with respect to temperature they obtain the absolute entropy of the HDL phase. The key findings can be summarized as follows:

1. Equilibration of HDL – Long Monte‑Carlo trajectories (tens of millions of steps) were generated at each temperature. At 1065 K and 1070 K the trajectories remain in the HDL basin long enough to be considered equilibrated. At 1060 K, however, two distinct types of trajectories appear: a “long” trajectory that stays in the HDL region for a substantial time before crossing a shallow free‑energy barrier, and “short” trajectories that quickly descend into low‑density amorphous or crystalline basins. The authors argue that the shallow barrier prevents full equilibration at 1060 K, leading to non‑equilibrium HDL states.

2. Volume‑pressure characteristics – By histogramming the instantaneous volume per particle v = V/N and counting the number of configurations N_c(v) in a small bin, they construct the Helmholtz free‑energy profile F(v) = –k_BT ln N_c + const. The resulting pressure‑volume (p‑v) curves display a region of nearly zero slope at 1070 K, indicative of a two‑phase coexistence. At 1060 K and 1065 K the p‑v curves acquire a positive slope and exhibit a classic van‑der‑Waals (VDW) loop. Double‑tangent constructions across the loop identify two states with identical chemical potential at the same temperature and pressure, confirming the presence of metastable HDL and low‑density liquid (LDL) branches.

3. Thermodynamic integration path – The excess Gibbs free energy is obtained via a three‑stage reversible λ‑integration path: (i) linear scaling down of the SW pair potential to approach an ideal‑gas‑like state under a volume constraint, (ii) imposition of Gaussian external wells to force particles onto a crystalline lattice, and (iii) gradual removal of the external wells while restoring the full SW interaction. For each λ‑value the ensemble average ⟨∂φ/∂λ⟩ is evaluated, and the free‑energy differences ΔG_i are integrated numerically. The Bennett Acceptance Ratio (BAR) method is employed to compute ΔG between adjacent λ‑states, ensuring high overlap and negligible hysteresis.

4. Entropy behavior – From the temperature dependence of G⁽ᵉ⁾ the absolute entropy S⁽ᵉ⁾ = –∂G⁽ᵉ⁾/∂T is extracted. For T ≥ 1065 K the entropy decreases with decreasing temperature, as expected for a stable liquid. Strikingly, at 1060 K the entropy of the HDL phase rises sharply when the system is in the non‑equilibrium state identified from the short trajectories. This increase is linked to a larger enthalpy reduction (ΔU < 0) together with a positive pΔV term, reflecting a substantial change in both internal energy and volume.

5. Free‑energy expansion around spinodals – The authors fit the Helmholtz free‑energy profile near the unstable region with a third‑order Taylor expansion about the spinodal volume v_s: F(v) ≈ F_s – p_s(v – v_s) + (1/3!)F’’’_s(v – v_s)³. The fitted coefficients reproduce the simulated F(v) very accurately, suggesting that the probability distribution of volumes is governed by the two spinodal points that delimit the VDW loop. The asymmetry of the loop implies that each branch corresponds to a distinct phase with its own equation of state, terminating at its respective spinodal.

6. Implications for HDL‑LDL transition – The presence of a VDW loop, the double‑tangent construction, and the non‑monotonic entropy signal that the HDL‑LDL transition in SW silicon is not a simple first‑order transition driven solely by latent heat. Instead, it involves a delicate balance of enthalpy, volume, and configurational entropy, with the free‑energy landscape becoming highly anharmonic near 1060 K. The shallow barrier permits the system to linger in a metastable HDL basin, during which the entropy can increase despite a drop in enthalpy.

Overall, the paper provides a rigorous thermodynamic quantification of the HDL phase in supercooled silicon, demonstrates the emergence of van‑der‑Waals loops in the p‑v isotherms, and reveals a non‑monotonic entropy trend that underlies the HDL‑to‑LDL transformation. These findings deepen our understanding of liquid‑liquid transitions in tetrahedral network formers and offer a solid computational framework for future studies of metastable liquids.