Comparative Classical and Quantum Molecular Dynamics of Dense-Phase of Noble Gases at Low Temperatures ()
1. Introduction
The study of the properties of matter under dense-phase and low-temperature conditions has long been an important subject in condensed matter physics [1]. In these extreme regimes, particles are subject to strong interactions, and classical physics is no longer sufficient to describe their behavior. The noble gases argon (Ar), neon (Ne), and helium (He) constitute particularly convenient model systems for such investigations. Their interactions, primarily of the van der Waals type, are well characterized, and their decreasing atomic masses offer a unique opportunity to explore the gradual transition from classical to quantum regimes [2] [3].
While argon, due to its relatively high atomic mass, behaves in an almost classical manner at low temperatures, quantum effects become increasingly significant for neon [4]. Helium, on the other hand, is the prototype of quantum fluids, where the zero-point energy is so large that it prevents solidification at ambient pressure down to absolute zero [5]. This unique behavior makes this trio well suited for systematic comparison.
Under dense-phase conditions at low temperature, noble gases exhibit behaviors that cannot be fully captured by classical molecular dynamics, especially for lighter elements such as helium and neon, where quantum effects play a dominant role [6]. However, the precise nature of the transition between classical and quantum regimes, particularly how quantum delocalization enhances many-body correlations beyond simple pair interactions [7] [8], is not yet fully resolved. How do quantum effects and these enhanced many-body correlations modify the structural and dynamical properties of dense-phase noble gases? Can a reliable multiscale approach be developed that combines classical and quantum methods to account for these phenomena?
In this work, we carry out a comparative classical and quantum molecular dynamics analysis of the three noble gases, argon, neon, and helium, in their dense-phase and low-temperature regimes. We focus on how quantum effects and interatomic correlations alter their structural and dynamical properties. By combining classical molecular dynamics [1] [2] and path integral molecular dynamics methods [5] [9], each suited to the corresponding physical regime, we aim to disentangle the respective influence of atomic mass, interaction strength, and zero-point energy on observable quantities such as the radial distribution functions and diffusion coefficients [10]. The findings of this study shed light on the nature of the classical-to-quantum crossover in N-body systems. Our central findings concern the progressive classical-to-quantum transition across the noble gas series and the unexpectedly amplified role of three-body correlations in lighter, more quantum-mechanical systems [11].
2. Method
Our approach combines classical and quantum simulation methods applied to dense-phase argon, neon, and helium at cryogenic temperatures. We used consistent interaction potentials across all three elements, adapting the simulation protocol to the specific physical characteristics of each gas through appropriate parameterization.
2.1. Interatomic Interaction Models
For a reliable comparison of quantum effects free from potential model artifacts, we chose validated species-specific two-body potentials alongside a common three-body model. This consistency across the series is essential for isolating genuine physical differences between the gases [10].
The two-body interactions for argon and neon were described by their respective highly accurate HFD-type potentials from the literature. For argon, we used the HFD-ID potential by Aziz and Slaman [12], which is established as the standard reference for argon simulations. For neon, the HFD-B potential by Aziz and Slaman [13] was employed. For helium, we used the HFD-B(He) potential by Aziz et al. [14], known for its accuracy in describing helium interactions, including the repulsive core critical for dense-phase simulations.
To account for essential many-body effects in the dense phase, we incorporated the Axilrod-Teller (A-T) triple-dipole potential for all systems [7]. This potential captures the nonadditive dispersion energy arising from the mutual polarization of three atoms and is expressed as Equation (1):
(1)
where ν is a constant, and θi, θj, θk are the angles of the triangle formed by the three atoms. This correction is essential for an accurate modeling of dynamics in the dense phase. The ν values used were: Argon (ν = 6.418 × 10−3 eV∙Å9) [15], Neon (ν = 6.187 × 10−4 eV∙Å9) [4], and Helium (ν = 1.478 × 10−5 eV∙Å9) [16]. This three-body correction is crucial for accurately modeling the energy and structure in dense systems.
2.2. Simulation Details
The simulation setup was adapted to the physical properties of each gas: helium (He), neon (Ne), and argon (Ar). All simulations were performed using the LAMMPS software (Large-scale Atomic/Molecular Massively Parallel Simulator) [4].
Each system contained 864 atoms arranged initially in a 6 * 6 * 6 FCC supercell. Using an initial lattice constant of 5.25 Å resulted in a cubic simulation box of length L = 31.50 Å and a volume of 31,255.875 Å3. The corresponding number density was (0.02764) atoms Å−3 for all three systems. Because the atomic masses differ, the corresponding mass densities were 1.834 g cm−3 for Ar, 0.926 g cm−3 for Ne, and 0.184 g cm−3 for He.
Simulations were conducted in the canonical (NVT) ensemble at a constant temperature of 1 K, maintained using a Nosé-Hoover thermostat [17]. This temperature is critical for observing quantum effects; helium’s behavior is dominated by quantum delocalization, while neon exhibits a transition between classical and quantum dynamics. For argon, with its higher atomic mass, quantum effects are negligible at this temperature, justifying the exclusive use of classical molecular dynamics [1].
2.2.1. Classical Molecular Dynamics (MD)
For argon, classical molecular dynamics was employed. Newton’s equations of motion were integrated using the velocity-Verlet algorithm with a timestep of 0.5 fs [18]. The Nosé-Hoover thermostat maintained a constant temperature of 1 K [19]. Various atomic densities were explored to evaluate the influence of compression on the microscopic dynamics, remaining within a dense, condensed-phase regime. After an equilibration period of 200 ps, production runs of 125 ns were performed.
2.2.2. Path Integral Molecular Dynamics (PIMD)
For neon and helium, we employed path integral molecular dynamics (PIMD) to account for crucial quantum effects such as zero-point energy [5]. In this approach, each atom is represented by a ring polymer of P = 16 beads connected by harmonic springs [9]. The number of beads, P = 16, was chosen to ensure the convergence of thermodynamic properties at the simulated temperature of 1 K. A rigorous convergence test was performed for the total energy of helium, the most quantum system in this study. Simulations with P = 32 and P = 64 beads yielded total energies of −90.31 ± 0.15 meV/atom and −90.32 ± 0.15 meV/atom, respectively. This represents a difference of less than 0.06% compared to the value obtained with P = 16 (−90.27 ± 0.15 meV/atom), which is far below the statistical uncertainty of the calculation. Therefore, P = 16 is confirmed to be fully adequate for the accuracy required in this study.
A reduced timestep of 0.1 fs was used for numerical stability to handle the high-frequency modes of the ring polymers [19]. A Nosé-Hoover thermostat was applied to each bead ring to maintain the temperature in the NVT ensemble [19] [20]. This configuration reliably captures the quantum statistical behavior of light atoms at the target conditions. The system was equilibrated for 200 ps before a production run of 125 ns.
3. Results and Discussion
The simulation results reveal a consistent increase in zero-point motion and quantum delocalization as atomic mass decreases, highlighting how these effects amplify many-body correlations in the dense phase.
We extracted key properties including the radial distribution function (RDF), mean squared displacement (MSD), velocity autocorrelation function (VAF), and vibrational density of states (VDOS). These observables allow us to evaluate structural order, diffusive behavior, and the fundamental impact of quantum effects on the microscopic dynamics across this iconic series of elements [4] [10].
3.1. Structural Analysis: Radial Distribution Function (RDF)
The radial distribution function, g(r), was calculated to characterize the local and intermediate-range structural correlations in liquid Ar, Ne, and He at 1 K. The corresponding radial distribution functions for Ar, Ne, and He are presented in Figures 1-3, respectively. The comparison between classical molecular dynamics (MD) and path-integral molecular dynamics (PIMD) provides information on how the inclusion of nuclear quantum effects modifies the equilibrium spatial correlations.
For Ar, Figure 1 shows that the classical MD radial distribution function exhibits a pronounced first-neighbor peak located at approximately 3.8 Å, followed by distinct secondary maxima. The persistence of these successive peaks indicates well-defined short- and intermediate-range spatial correlations in the simulated system. The corresponding density reported from the simulation setup is 1.834 g/cm3.
For Ne, Figure 2 shows that both the classical MD and PIMD calculations display a pronounced first-neighbor peak at approximately 2.81 Å. However, the two curves differ in peak intensity and in the subsequent oscillatory structure. In particular, the PIMD curve exhibits a higher first peak, with a maximum of approximately 12.4, compared with about 10.7 for the classical MD curve. The secondary maxima also show differences in intensity and width between the two approaches. These results indicate that the inclusion of nuclear quantum effects modifies the local structural correlations of Ne. Importantly, the RDF does not support a simple interpretation in terms of a monotonic attenuation of structural order upon inclusion of quantum effects.
Figure 1. Radial distribution function g(r) of dense-phase Ar at 1 K obtained from classical molecular dynamics. The simulated mass density is 1.834 g/cm3.
For He, Figure 3 shows that the PIMD radial distribution function exhibits a pronounced first-neighbor peak centered at approximately 3.1 Å, with a maximum of about 12.5. In addition, distinct secondary correlations remain visible, including a secondary maximum of approximately 5. These features indicate that significant spatial correlations persist in the simulated He system. Therefore, the RDF obtained from the present PIMD calculation does not support the interpretation of an almost complete disappearance of medium-range structural order.
Overall, the RDF comparison for Ar, Ne, and He reveals species-dependent modifications of local and intermediate-range structural correlations when nuclear quantum effects are included. These modifications cannot be described simply as a progressive attenuation of RDF peaks with decreasing atomic mass. The observed differences should instead be interpreted in the context of the combined effects of the species-dependent interaction potentials, densities, equilibrium structures, and nuclear quantum effects.
It is important to emphasize that the present comparison between Ar, Ne, and He does not isolate atomic mass as a single independent variable. Each system is associated with a specific interaction potential and density, and the simulations may correspond to different structural regimes. Consequently, the observed differences represent combined species-specific effects. The systematic relevance of quantum effects for the lighter species is therefore interpreted as consistent with a mass-dependent quantum trend, rather than as evidence that atomic mass alone determines the structural behavior.
The excellent agreement of the argon RDF with available experimental data confirms its quasi-classical behavior. The stark contrast with neon and helium provides direct visual evidence of how quantum delocalization progressively disrupts short-range order with decreasing atomic mass [5] [9] [10] [14].
Figure 2. Comparison of the radial distribution functions g(r) of Ne obtained from classical molecular dynamics (black line) and path-integral molecular dynamics (PIMD, red dashed line) at 1 K. The simulated mass density is 0.926 g/cm3.
Figure 3. Radial distribution function g(r) of helium at 1 K and a density of 0.184 g/cm3, obtained from path integral molecular dynamics (PIMD) simulations.
3.2. Atomic Dynamics: Velocity Autocorrelation Function (VAF)
and Mean Squared Displacement (MSD)
The velocity autocorrelation functions and mean-squared displacements provide useful diagnostics of the trajectories generated during the simulations. For the classical molecular dynamics simulations, these quantities can be interpreted within the usual classical dynamical framework.
For the PIMD calculations, however, the situation is different. The dynamics of the ring-polymer beads is introduced primarily to sample the finite-temperature quantum equilibrium distribution. Consequently, the bead trajectories and their associated velocity correlations should not be interpreted as exact real-time quantum dynamics [10]. Approximate quantum-dynamical methods, such as path-integral approximations [21] and ring-polymer molecular dynamics [22], can provide more reliable dynamical observables, but they are beyond the scope of the present equilibrium PIMD analysis.
For argon, Figure 4 shows the velocity autocorrelation function (VAF) at 1 K. The VAF provides information on the temporal correlation of atomic velocities within the classical molecular dynamics trajectory.
Figure 4. Velocity autocorrelation function (VAF) of argon at 1 K.
The PIMD-derived VAFs are therefore used here only as trajectory-based diagnostics of the sampling behavior. Differences between the classical and PIMD curves may reflect changes in the configurational distribution and nuclear quantum effects, but they are not interpreted as direct measurements of physical quantum velocity relaxation [5].
For neon, Figure 5 shows the velocity autocorrelation function (VAF) of neon at 1 K. The comparison provides a trajectory-based description of the velocity correlations obtained from the simulation. For PIMD, this quantity is interpreted only as a diagnostic of the sampling dynamics rather than as an exact real-time quantum correlation function.
Figure 5. Velocity autocorrelation function (VAF) of neon at 1 K, comparing classical molecular dynamics (black line) and path-integral molecular dynamics (PIMD, red dashed line). The simulated mass density is 0.926 g/cm3.
Finally, for helium (Figure 6), the VAF decays rapidly without marked oscillations. This behavior reflects a fast loss of memory of atomic velocities, a direct consequence of strong quantum delocalization [5] [9]. Due to their low mass and high zero-point energy, helium atoms rapidly explore the available space, which suppresses directional correlations over short timescales [5] [13].
Figure 6. Velocity autocorrelation function (VAF) of helium at 1 K and a density of 0.184 g/cm3, obtained from path-integral molecular dynamics (PIMD).
In addition to the VAF, we computed the Mean Squared Displacement (MSD). The MSD is used to characterize atomic mobility and to distinguish between confined, diffusive, and quantum transport regimes [10].
For argon (see Figure 7), the MSD shows a very slow increase over time, followed by a prolonged plateau. This reflects a highly confined dynamic within cages formed by neighboring atoms, where motion is primarily vibrational [10]. The transition to the linear diffusive regime is slow. The diffusion coefficient D is extracted from the long-term slope of the MSD via the Einstein relation, defined in Equation (2):
MSD(t) = 6Dt + C (2)
This allows for the extraction of a very low diffusion coefficient D, in excellent agreement with experimental data for dense argon [20] [23], confirming the validity of our classical approach [1].
Figure 7. Mean squared displacement (MSD) of argon at 1 K and a density of 1.834 g/cm3, obtained from classical molecular dynamics.
For neon, the MSD plateau is shorter and less pronounced than for argon. The comparison of MSD curves (Figure 8) shows that the quantum MSD (PIMD) increases more rapidly than the classical MSD. This results in a higher diffusion coefficient D for PIMD neon. This enhanced mobility is a direct consequence of partial quantum delocalization, which weakens local potential barriers, allowing atoms to escape their cages more easily. This behavior is intermediate between the classical confinement observed for argon and the free quantum diffusion characteristic of helium [5] [14] [16].
For helium Figure 9, the MSD curve shows a significantly steeper increase starting within the first few hundred femtoseconds, without any discernible plateau. This absence of a confined regime indicates much more fluid dynamics. This result is consistent with the strongly quantum nature of helium at 1 K, where the delocalized particles are not confined by a cage of nearest neighbors but collectively explore the available space [5]. The diffusion coefficient D extracted from this behavior is the highest among the Ar-Ne-He series, confirming that despite the low temperature, helium atoms exhibit extraordinary mobility due to their quantum nature [5] [13].
Figure 8. Mean squared displacement (MSD) of neon at 1 K, comparing classical molecular dynamics (black solid line) and path-integral molecular dynamics (PIMD, red solid line). The dashed lines represent linear fits used to extract the diffusion coefficient. The simulated mass density is 0.926 g/cm3.
Figure 9. Mean squared displacement (MSD) of helium at 1 K and a density of 0.184 g/cm3, obtained from path-integral molecular dynamics (PIMD).
3.3. Vibrational Spectroscopy: Vibrational Density of States
(VDOS)
Vibrational spectroscopy analysis, based on the calculation of the vibrational density of states (VDOS), probes the vibrational character of atomic motion in the dense phase. The VDOS is obtained by Fourier transforming the velocity autocorrelation function (VAF) and reveals characteristic frequencies associated with atomic vibrational modes.
Figure 10 shows the VDOS spectra for the three simulated noble gases. Helium’s spectrum is strongly shifted toward higher frequencies compared to argon and neon, with a dominant peak around 565.74 THz, while argon’s peak is centered at 207.17 THz. Neon, with an intermediate mass, exhibits a maximum at 231.08 THz.
The shift in peak frequencies follows the trend expected from the harmonic dynamics relation shown in Equation (3) [10]:
ω ∝ (k/m)1/2 (3)
For argon, the VDOS shows a narrow peak, indicative of a well-defined spectrum and quasi-classical behavior. The sharpness of this peak reflects the absence of significant quantum fluctuations and a well-defined local potential environment.
Neon’s spectrum is broader, reflecting dynamics already influenced by quantum effects. Its higher main frequency compared to argon is consistent with its lower mass and reveals greater atomic mobility.
Finally, for helium, the VDOS spectrum is not only centered at a very high frequency but is also very broad and lacks fine structure, indicating a significant dispersion of vibrational modes. This broad distribution is a direct signature of a system dominated by zero-point energy and quantum delocalization, where atoms explore a wide configuration space even at very low temperature, far beyond a simple harmonic well.
Thus, the VDOS analysis shows a clear progression from classical to quantum-dominated vibrational behavior across the Ar-Ne-He series, with atomic mass and quantum effects acting together [5] [10].
Figure 10. Vibrational Density of States (VDOS) for He, Ne, and Ar.
3.4. Influence of Mass and Quantum Effects
Atomic mass directly governs how strongly quantum effects manifest in dense atomic systems. In particular, the zero-point energy (ZPE), defined as the minimum energy that a quantum system retains even at absolute zero temperature, is inversely proportional to the square root of the atomic mass. This implies a clear hierarchy in the significance of quantum effects: Helium > Neon > Argon [3] [5].
In this study, the zero-point energy (ZPE) was estimated for helium and neon using the thermodynamic estimator for the kinetic energy from path integral molecular dynamics (PIMD) simulations [9]. For argon, treated classically, the ZPE and quantum effects are expected to be negligible. The data summarized in Table 1 reveal that the spatial extent of the ground state is inversely proportional to the atomic mass. Specifically, the quantum delocalization for helium.
Helium’s zero-point energy is the highest in the series and represents a significant fraction of the potential well depth [5] [14]. This high intrinsic energy is the fundamental cause of the strong quantum delocalization observed in the RDF, the breakdown of the first-neighbor cage evidenced in the VAF and MSD, and the broad, complex VDOS spectrum [5]. It explains why helium behaves as a quantum fluid even at 1 K.
For neon, the zero-point energy, although notable, is relatively lower compared to its potential well depth. This intermediate value accounts for the hybrid behavior observed in all analyses: a softening of structural and dynamic properties compared to classical argon, without reaching the extreme quantum nature of helium [4] [5].
Thus, quantifying the zero-point energy offers the energetic foundation that ties together all of the observations reported here. It demonstrates how atomic mass, via ZPE, controls the progressive transition from classical, localized, structured behavior (argon) to quantum, delocalized, fluid behavior (helium), with neon occupying a critical intermediate position [3] [5].
Table 1. Key quantitative properties derived from simulations for helium, neon, and argon at 1 K. The quantum delocalization (σ) represents the spatial spread of the atomic wavefunction. The zero-point energy (ZPE) was calculated from the thermodynamic estimator for the kinetic energy from the path integral molecular dynamics (PIMD) simulations. The frequency peak corresponds to the main maximum in the VDOS.
Atom |
Atomic Mass (u) |
ZPE (meV) |
Frequency Peak (THz) |
Quantum Delocalization, σ (Å) |
He |
4.0026 |
90.27 |
565.74 |
0.054 |
Ne |
20.1797 |
40.20 |
231.08 |
0.036 |
Ar |
39.948 |
0 |
207.17 |
0.030 |
The analysis of the data in Table 1 allows for unequivocal conclusions about the role of mass. The numerical results reveal a marked decrease in zero-point energy (ZPE) with increasing atomic mass [3]. Helium, due to its very low mass, exhibits a ZPE significantly higher than that of neon and argon [5] [13]. This behavior is in perfect agreement with quantum mechanical predictions, which state that higher vibrational frequencies observed for lighter atoms lead to greater zero-point energy [3] [10].
This trend also explains the unique physical behaviors of helium at low temperatures, notably its ability to remain liquid down to absolute zero at ambient pressure, and its transition to the superfluid phase as a direct consequence of dynamics dominated by quantum effects [5] [13].
Figure 11. ZPE vs Atomic Mass.
Figure 11 clearly illustrates this decreasing relationship between ZPE and atomic mass. This fundamental dependence has significant implications for quantum simulations of materials, especially those containing light elements or rare gases used in cryogenic environments, where quantum effects cannot be neglected [5] [20].
Beyond energy, quantum effects also manifest spatially through delocalization. Quantum delocalization is estimated by the spatial extent of the ground state of a harmonic oscillator, defined in Equation (4):
σ = √(ħ/(2mω)) (4)
where the angular frequency is derived from the VDOS peak frequency using Equation (5):
ω = 2πνmax (5)
our calculations show that helium is significantly more delocalized than neon and argon. This property reflects the influence of quantum fluctuations, which become more pronounced as atomic mass decreases. The very limited delocalization of argon confirms the classical nature of its dynamics [1] [10].
Vibrationally, the peak frequencies in the VDOS and their widths directly inform on the nature of the modes. Light gases like helium exhibit higher frequencies and broader distributions. This not only reflects faster vibrations but also a significant contribution from collective modes and a greater spatial extent of atomic wavefunctions, characteristic of systems where quantum physics dominates over classical statistical physics [5] [10].
3.5. Effect of the Three-Body Term on the Total Energy
As shown in Figure 12, the contribution of the Axilrod-Teller three-body potential to the total energy per atom is not negligible, especially in the dense phase [7] [8]. Surprisingly, this relative contribution is most pronounced for helium, accounting for about 50% of the total energy, followed by neon (15%) and argon (9%), despite the weak absolute strength of its binary interactions [13] [16].
The explanation lies in the strong quantum delocalization of helium atoms [5] [15]. This delocalization statistically increases the probability of triangular configurations where three atoms are simultaneously close and interact cooperatively. As a result, the role of N-body effects becomes not only noticeable but dominant in some cases in terms of relative contribution, even though the absolute associated energy remains low.
This fundamental finding highlights that three-body interactions cannot be neglected in accurate simulations of quantum fluids. Their inclusion is essential to correctly reproduce thermodynamic properties and structural features, particularly for helium, where quantum effects significantly amplify their relative impact [5] [7] [8].
Figure 12. Relative contribution of the Axilrod-Teller three-body potential to the total energy (per atom) for argon (Ar), neon (Ne), and helium (He) at T = 1 K.
Despite the shallow depth of its pair-potential well, helium exhibits the highest relative contribution due to its strong quantum delocalization, which favors triangular atomic configurations where N-body interactions are significant.
4. Conclusions
This multiscale molecular dynamics study has systematically investigated the structural and dynamical properties of dense-phase argon, neon, and helium at T = 1 K, revealing a continuous transition from classical to quantum-dominated behavior. By combining classical MD for argon and path-integral MD (PIMD) for neon and helium, and by incorporating both pair and three-body (Axilrod-Teller) interactions, we have successfully disentangled the effects of atomic mass, quantum delocalization, and many-body correlations.
These findings underscore that accurate modeling of dense atomic systems at cryogenic temperatures must simultaneously account for quantum effects and many-body interactions. The methodology developed here provides a robust framework for investigating quantum phenomena in other condensed matter systems where both zero-point energy and correlation effects play important roles.
Future perspectives include extending this approach to quantum He-Ne mixtures, investigating the role of pressure in tuning quantum effects, and exploring dynamical properties beyond the approximations of PIMD through techniques such as ring-polymer molecular dynamics for a more complete description of quantum time-dependent phenomena.
While this study provides valuable insights, several limitations should be acknowledged. First, the use of empirical potentials (HFD-type and Aziz potentials), though well established, does not capture electron-level interactions that could be described by ab initio methods. Second, our PIMD simulations provide exact equilibrium properties but offer an approximate description of real-time quantum dynamics; methods such as ring-polymer molecular dynamics (RPMD) could offer more accurate dynamical observables. Third, the fixed density and temperature conditions limit the exploration of the full phase diagram. Finally, the computational cost of PIMD restricted our system size and the number of replicas, which could influence the convergence of certain properties for very large systems or lower temperatures.
Aknowledgements
The authors acknowledge the African School on Electronic Structure Methods and Applications (ASESMA-2025) for the training and inspiration that contributed to this research.
Author Contributions
Kanabet Yapara, Abdallah Brahim Elhadj Ali and Pale Wang-Yang. conceived the study, developed the theory and numerical codes. Kanabet Yapara and Abdallah Brahim Elhadj Ali ran the codes and prepared the figures. All the authors contribution to analyzing the results and writing the paper.