Path Integral Monte Carlo Calculations Research Articles

Energy response and spatial alignment of the perturbed electron gas.

We study the linear energy response of the uniform electron gas to an external harmonic perturbation with a focus on resolving different contributions to the total energy. This has been achieved by carrying out highly accurate abinitio path integral Monte Carlo (PIMC) calculations for a variety of densities and temperatures. We report a number of physical insights into effects such as screening and the relative importance of kinetic and potential energies for different wave numbers. A particularly interesting finding is obtained from the observed non-monotonic behavior of the induced change in the interaction energy, which becomes negative for intermediate wave numbers. This effect is strongly dependent on the coupling strength and constitutes further direct evidence for the spatial alignment of electrons introduced in earlier works [T. Dornheim et al., Commun. Phys. 5, 304 (2022)]. The observed quadratic dependence on the perturbation amplitude in the limit of weak perturbations and the quartic dependence of perturbation amplitude corrections are consistent with linear and nonlinear versions of the density stiffness theorem. All PIMC simulation results are freely available online and can be used to benchmark new methods or as input for other calculations.

Observation of Multiple Ordered Solvation Shells in Doped Helium Droplets: The Case of HeNCa2.

In this Letter, we report the experimental detection of likely the largest ordered structure of helium atoms surrounding a monatomic impurity observed to date using a recently developed technique. The mass spectrometry investigation of HeNCa2+ clusters, formed in multiply charged helium nanodroplets, reveals magic numbers at N = 12, 32, 44, and 74. Classical optimization and path integral Monte Carlo calculations suggest the existence of up to four shells surrounding the calcium dication which are closed with well-ordered Mozartkugel-like structures: He12Ca2+ with an icosahedron, the second at He32Ca2+ with a dodecahedron, the third at He44Ca2+ with a larger icosahedron, and finally for He74Ca2+, we find that the outermost He atoms form an icosidodecahedron which contains the other inner shells. We analyze the reasons for the formation of such ordered shells in order to guide the selection of possible candidates to exhibit a similar behavior.

Vacancy-induced supersolidity of the second He4 layer on a biphenylene carbon sheet

We have performed path-integral Monte Carlo calculations to investigate various quantum phases of $^{4}\mathrm{He}$ layers adsorbed on a biphenylene sheet, a recently synthesized two-dimensional network of carbon atoms. Three different commensurate solid phases are identified in the first $^{4}\mathrm{He}$ layer, which is found to be completed to a ${\mathrm{C}}_{2/2}$ commensurate solid at an areal density of 0.119 ${\AA{}}^{\ensuremath{-}2}$. We have also found that the commensurate structure of the completed first layer allows second-layer $^{4}\mathrm{He}$ atoms to constitute a stable commensurate ${\mathrm{C}}_{2/3}$ structure at a total helium coverage of 0.199 ${\AA{}}^{\ensuremath{-}2}$. Furthermore, this second-layer commensurate structure is observed to be stable even with the formation of vacancies that are mobile to instigate the superfluid response at low temperatures. This unequivocally shows that vacancy-induced supersolidity can be realized in the second $^{4}\mathrm{He}$ layer on biphenylene.

Collision-induced three-body polarizability of helium.

We present the first-principles determination of the three-body polarizability and the third dielectric virial coefficient of helium. Coupled-cluster and full configuration interaction methods were used to perform electronic structure calculations. The mean absolute relative uncertainty of the trace of the polarizability tensor, resulting from the incompleteness of the orbital basis set, was found to be 4.7%. Additional uncertainty due to the approximate treatment of triple and the neglect of higher excitations was estimated at 5.7%. An analytic function was developed to describe the short-range behavior of the polarizability and its asymptotics in all fragmentation channels. We calculated the third dielectric virial coefficient and its uncertainty using the classical and semiclassical Feynman-Hibbs approaches. The results of our calculations were compared with experimental data and with recent Path-Integral Monte Carlo (PIMC) calculations [Garberoglio et al., J. Chem. Phys. 155, 234103 (2021)] employing the so-called superposition approximation of the three-body polarizability. For temperatures above 200K, we observed a significant discrepancy between the classical results obtained using superposition approximation and the ab initio computed polarizability. For temperatures from 10K up to 200K, the differences between PIMC and semiclassical calculations are several times smaller than the uncertainties of our results. Except at low temperatures, our results agree very well with the available experimental data but have much smaller uncertainties. The data reported in this work eliminate the main accuracy bottleneck in the optical pressure standard [Gaiser et al., Ann. Phys. 534, 2200336 (2022)] and facilitate further progress in the field of quantum metrology.

Path-integral quantum Monte Carlo calculations of light nuclei

We describe a path-integral ground-state quantum Monte Carlo method for light nuclei in continuous space. We show how to efficiently update and sample the paths with spin-isospin dependent and spin-orbit interactions. We apply the method to the triton and alpha particle using both local chiral interactions with next-to-next-to-leading-order %(N$^2$LO) and the Argonne interactions. For operators, like the total energy, that commute with the Hamiltonian, our results agree with Green's function Monte Carlo and auxiliary field diffusion Monte Carlo calculations. For operators that do not commute with the Hamiltonian and for Euclidean response functions, the path-integral formulation allows straightforward calculation without forward walking or the increased variance typical of diffusion methods. We demonstrate this by calculating density distributions, root mean square radii, and Euclidean response functions for single-nucleon couplings.

Path-integral Monte Carlo simulations on the thermodynamic properties of single-layer hexagonal boron nitride

We investigate the quantum effects and anharmonicity on the thermodynamical properties of single-layer hexagonal boron nitride (SL-hBN) using quantum path-integral Monte Carlo (PIMC) and classical Monte Carlo (CMC) simulations in the isothermal–isobaric ensemble at zero pressure. The quantum effects are extracted from the differences between the results obtained from the PIMC and CMC calculations. Our numerical simulations show that the out-of-plane motion of the atoms plays a decisive role in the negative thermal expansion of the SL-hBN crystal and the quantum zero-point vibrational energy greatly enhances the negative thermal expansion coefficients. We demonstrate for the first time the anisotropic thermal expansions in SL-hBN crystal manifested in the different linear thermal expansion coefficients in the zig-zag and armchair directions. Such an anisotropic effect can be significant at high temperatures. From an analysis of the atomic kinetic and potential energy contributions in the thermal vibration energy of the crystal lattice, we find that the anharmonicity of the SL-hBN crystal lies on the ground of the asymmetric thermal vibrations of the boron and nitrogen atoms and their coupling in the same unit cell of the hexagonal lattice.

The relevance of electronic perturbations in the warm dense electron gas.

Warm dense matter (WDM) has emerged as one of the frontiers of both experimental physics and theoretical physics and is a challenging traditional concept of plasma, atomic, and condensed-matter physics. While it has become common practice to model correlated electrons in WDM within the framework of Kohn-Sham density functional theory, quantitative benchmarks of exchange-correlation (XC) functionals under WDM conditions are yet incomplete. Here, we present the first assessment of common XC functionals against exact path-integral Monte Carlo calculations of the harmonically perturbed thermal electron gas. This system is directly related to the numerical modeling of x-ray scattering experiments on warm dense samples. Our assessment yields the parameter space where common XC functionals are applicable. More importantly, we pinpointwhere the tested XC functionals fail when perturbations on the electronic structure are imposed. We indicate the lack of XC functionals that take into account the needs of WDM physics in terms of perturbed electronic structures.

Experimental and theoretical determinations of hydrogen isotopic equilibrium in the system CH4[sbnd]H2[sbnd]H2O from 3 to 200 °C

The stable isotopic composition of methane (CH4) is commonly used to fingerprint natural gas origins. Over the past 50 years, there have been numerous proposals that both microbial and thermogenic CH4 can form in or later attain hydrogen isotopic equilibrium with water (H2O) and carbon isotopic equilibrium with carbon dioxide (CO2). Evaluation of such proposals requires knowledge of the equilibrium fractionation factors between CH4 and H2O or CO2 at the temperatures where microbial and thermogenic CH4 form in or are found in the environment, which is generally less than 200 °C. Experimental determinations of these fractionation factors are only available above 200 °C, requiring extrapolation of these results beyond the calibrated range or the use of theoretical calculations at lower temperatures. Here, we provide a calibration of the equilibrium hydrogen isotopic fractionation factor for CH4 and hydrogen gas (H2) (DαCH4(g)–H2(g)) based on experiments using γ-Al2O3 and Ni catalysts from 3 to 200 °C. Results were regressed as a 2nd order polynomial of 1000 × lnDαCH4(g)–H2(g) vs. 1/T (K−1) yielding:1000×lnDαCH4(g)-H2(g)=3.5317×107T2+2.7749×105T-179.48We combine this calibration with previous experimental determinations of hydrogen isotope equilibrium between H2, H2O(g), and H2O(l) and we provide an interpolatable experimental calibration of 1000 × lnDαCH4(g)–H2O(l) from 3 to 200 °C. Our resulting 4th order polynomial is the following equation:1000×lnDαCH4(g)-H2Ol=-7.9443×1012T4+8.7772×1010T3-3.4973×108T2+5.4398×105T-382.05At 3 °C, the value from our calibration differs by 93‰ relative to what would be calculated based on the extrapolation of the only experimental calibration currently available to temperatures below its calibrated range (lowest temperature of 200 °C; Horibe and Craig, 1995). We additionally provide new theoretical estimates of hydrogen isotopic equilibrium between CH4(g), H2(g), and H2O(g) and carbon isotopic equilibrium between CH4(g) and CO2(g) using Path Integral Monte Carlo (PIMC) calculations. Our PIMC calculations for hydrogen isotopic equilibrium between CH4 and H2 agree 1:1 with our experiments. Finally, we compile carbon and hydrogen isotopic measurements of CH4, CO2, and H2O from various environmental systems and compare observed differences between carbon and hydrogen isotopes to those expected based on isotopic equilibrium. We find that isotopic compositions of some microbial gases from marine sedimentary, coalbed, and shale environments are consistent with those expected for CH4H2O(l) hydrogen and CH4CO2 carbon isotopic equilibrium. In contrast, microbial terrestrial and pure culture gases are not consistent with both CH4H2O(l) hydrogen and CH4CO2 carbon isotopic equilibrium. These results are explained qualitatively using previously developed conceptual models that link free energy gradients available to microorganisms to the degree that their enzymes can promote isotope-exchange reactions between CH4, CO2, and H2O.

Path-integral calculation of the fourth virial coefficient of helium isotopes.

We use the path-integral Monte Carlo (PIMC) method and state-of-the-art two-body and three-body potentials to calculate the fourth virial coefficients D(T) of 4He and 3He as functions of temperature from 2.6 K to 2000 K. We derive expressions for the contributions of exchange effects due to the bosonic or fermionic nature of the helium isotope; these effects have been omitted from previous calculations. The exchange effects are relatively insignificant for 4He at the temperatures considered, but for 3He, they are necessary for quantitative accuracy below about 4 K. Our results are consistent with previous theoretical work (also with some of the limited and scattered experimental data) for 4He; for 3He, there are no experimental values, and this work provides the first values of D(T) calculated at this level. The uncertainty of the results depends on the statistical uncertainty of the PIMC calculation, the estimated effect of omitting four-body terms in the potential energy, and the uncertainty contribution propagated from the uncertainty of the potentials. At low temperatures, the uncertainty is dominated by the statistical uncertainty of the PIMC calculations, while at high temperatures, the uncertainties related to the three-body potential and omitted higher-order contributions become dominant.

Nonideal mixing effects in warm dense matter studied with first-principles computer simulations.

We study nonideal mixing effects in the regime of warm dense matter (WDM) by computing the shock Hugoniot curves of BN, MgO, and MgSiO3. First, we derive these curves from the equations of state (EOS) of the fully interacting systems, which were obtained using a combination of path integral Monte Carlo calculations at high temperature and density functional molecular dynamics simulations at lower temperatures. We then use the ideal mixing approximation at constant pressure and temperature to rederive these Hugoniot curves from the EOS tables of the individual elements. We find that the linear mixing approximation works remarkably well at temperatures above ∼2 × 105 K, where the shock compression ratio exceeds ∼3.2. The shape of the Hugoniot curve of each compound is well reproduced. Regions of increased shock compression, which emerge because of the ionization of L and K shell electrons, are well represented, and the maximum compression ratio of the Hugoniot curves is reproduced with high precision. Some deviations are seen near the onset of the L shell ionization regime, where ionization equilibrium in the fully interacting system cannot be well reproduced by the ideal mixing approximation. This approximation also breaks down at lower temperatures, where chemical bonds play an increasingly important role. However, the results imply that the equilibrium properties of binary and ternary mixtures in the regime of WDM can be derived from the EOS tables of the individual elements. This significantly simplifies the characterization of binary and ternary mixtures in the WDM and plasma phases, which otherwise requires large numbers of more computationally expensive first-principles computer simulations.

Two-dimensional superfluidity in He4 clusters intercalated into graphite

Motivated by recent technical developments on selective intercalation of inert gases just beneath the surface graphene of graphite, we perform path-integral Monte Carlo calculations to study structural and superfluid properties of $^{4}\mathrm{He}$ clusters encapsulated between a bulged surface graphene layer and the graphite (0001) substrate. The intercalated $^{4}\mathrm{He}$ atoms initially decorate the inner surface of the bulged graphene to form a shell structure, leaving a void inside it. Then, the additional $^{4}\mathrm{He}$ atoms form a disk-shaped liquid platelet in the void, which shows finite superfluid response at temperatures below 1 K. The temperature-dependent superfluid fractions of the platelet follow modified Kosterlitz-Thouless recursion relations, indicating that two-dimensional superfluidity can be realized in an intercalated $^{4}\mathrm{He}$ system.

Solving fermion problems without solving the sign problem: Symmetry-breaking wave functions from similarity-transformed propagators for solving two-dimensional quantum dots.

It is well known that the use of the primitive second-order propagator in path-integral Monte Carlo calculations of many-fermion systems leads to the sign problem. This work will show that by using the similarity-transformed Fokker-Planck propagator, it is possible to solve for the ground state of a large quantum dot, with up to 100 polarized electrons, without solving the sign problem. These similarity-transformed propagators naturally produce rotational symmetry-breaking ground-state wave functions previously used in the study of quantum dots and quantum Hall effects. However, instead of localizing the electrons at positions that minimize the potential energy, this derivation shows that they should be located at positions that maximize the bosonic ground-state wave function. Further improvements in the energy can be obtained by using these as initial wave functions in a ground-state path-integral Monte Carlo calculation with second- and fourth-order propagators.

Quantum Monte Carlo study of strongly interacting bosonic one-dimensional systems in periodic potentials

We present diffusion Monte Carlo (DMC) and path-integral Monte Carlo (PIMC) calculations of a one-dimensional Bose system with realistic interparticle interactions in a periodic external potential. Our main aim is to test the predictions of the Luttinger liquid (LL) theory, in particular with respect to the superfluid-Mott insulator transition at both zero and finite temperatures, in the predicted robust and fragile superfluid regimes. For that purpose, we present our results of the superfluid fraction $\rho_s/\rho_0$, the one-body density matrix, the two-body correlation functions, and the static structure factor. The DMC and PIMC results in the limit of very low temperature for $\rho_s/\rho_0$ agree, but the LL model for scaling $\rho_s/\rho_0$ does not fit the data well. The critical depth of the periodic potential is close to the values obtained for ultracold gases with different models of interaction, but with the same value of the bare LL parameter, demonstrating the universality of LL description. Algebraic decay of correlation functions is observed in the superfluid regime and exponential decay in the Mott-insulator one, as well as in all regimes at finite temperature for distances larger than a characteristic length.

Equation of state and first principles prediction of the vibrational matrix shift of solid parahydrogen.

We generate the equation of state (EOS) of solid parahydrogen (para-H2) using a path-integral Monte Carlo (PIMC) simulation based on a highly accurate first-principles adiabatic hindered rotor potential energy curve for the para-H2 dimer. The EOS curves for the fcc and hcp structures of solid para-H2 near the equilibrium density show that the hcp structure is the more stable of the two, in agreement with experiment. To accurately reproduce the structural and energy properties of solid para-H2, we eliminated by extrapolation the systematic errors associated with the choice of simulation parameters used in the PIMC calculation. We also investigate the temperature dependence of the EOS curves, and the invariance of the equilibrium density with temperature is satisfyingly reproduced. The pressure as a function of density and the compressibility as a function of pressure are both calculated using the obtained EOS and are compared with previous simulation results and experiments. We also report the first ever a priori prediction of a vibrational matrix shift from first-principles two-body potential functions, and its result for the equilibrium state agrees well with experiment.

A combined experimental and theoretical investigation of Cs+ ions solvated in HeN clusters.

Solvation of Cs+ ions inside helium droplets has been investigated both experimentally and theoretically. On the one hand, mass spectra of doped helium clusters ionized with a crossed electron beam, HeNCs+, have been recorded for sizes up to N = 60. The analysis of the ratio between the observed peaks for each size N reveals evidences of the closure of the first solvation shell when 17 He atoms surround the alkali ion. On the other hand, we have obtained energies and geometrical structures of the title clusters by means of basin-hopping, diffusion Monte Carlo (DMC), and path integral Monte Carlo (PIMC) methods. The analytical He-Cs+ interaction potential employed in our calculations is represented by the improved Lennard-Jones expression optimized on high level ab initio energies. The weakness of the existing interaction between helium and Cs+ in comparison with some other alkali ions such as Li+ is found to play a crucial role. Our theoretical findings confirm that the first solvation layer is completed at N = 17 and both evaporation and second difference energies obtained with the PIMC calculation seem to reproduce a feature observed at N = 12 for the experimental ion abundance. The analysis of the DMC probability distributions reveals the important contribution from the icosahedral structure to the overall configuration for He12Cs+.

Symmetry-changing commensurate-incommensurate solid transition in the He4 monolayer on 6,6,12-graphyne

Path-integral Monte Carlo calculations have been carried out to investigate physical properties of a $^{4}\mathrm{He}$ monolayer adsorbed on a single 6,6,12-graphyne sheet, which is one of the graphyne families possessing a rectangular symmetry. To characterize elusive quantum phases of an adsorbed $^{4}\mathrm{He}$ monolayer on 6,6,12-graphyne, we model the $^{4}\mathrm{He}$-graphyne interaction by the pairwise sum of empirical $^{4}\mathrm{He}\text{\ensuremath{-}}\mathrm{C}$ interatomic potentials. At partially filled $^{4}\mathrm{He}$ coverages, we identify three commensurate solids of the ${C}_{3/4},{C}_{4/4}$, and ${C}_{6/4}$ structures from the two-dimensional density distribution. These solids show the rectangular symmetry inherited from the symmetry of 6,6,12-graphyne, which were confirmed with the analysis of their static structure factors. At high helium coverages near its completion, the $^{4}\mathrm{He}$ monolayer is predicted to exhibit a transition from a rectangular commensurate structure to a triangular incommensurate structure, after going through inhomogeneous structures mixed with domains of triangular and rectangular orderings. This symmetry-changing transition has not been observed in $^{4}\mathrm{He}$ monolayers adsorbed on other carbon substrates.

Path integral Monte Carlo calculation of the momentum distribution of the homogeneous electron gas at finite temperature

Path integral Monte Carlo (PIMC) simulations are employed to calculate the momentum distribution of the homogeneous electron gas at finite temperature. Open paths are introduced to sample off-diagonal elements of the real-space density matrix. It is demonstrated how the restricted PIMC method can be extended to incorporate open paths in order to allow for the simulations in fermionic systems where a sign problem is present. The computed momentum distribution shows significant deviations from behavior of free fermions when strong correlations are present but agrees with predictions from variational methods.

Path integral Monte Carlo simulations of warm dense aluminum

We perform first-principles path integral Monte Carlo (PIMC) and density functional theory molecular dynamics (DFT-MD) calculations to explore warm dense matter states of aluminum. Our equation of state (EOS) simulations cover a wide density-temperature range of 0.1-32.4gcm^{-3} and 10^{4}-10^{8} K. Since PIMC and DFT-MD accurately treat effects of the atomic shell structure, we find two compression maxima along the principal Hugoniot curve attributed to K-shell and L-shell ionization. The results provide a benchmark for widely used EOS tables, such as SESAME, QEOS, and models based on Thomas-Fermi and average-atom techniques. A subsequent multishock analysis provides a quantitative assessment for how much heating occurs relative to an isentrope in multishock experiments. Finally, we compute heat capacity, pair-correlation functions, the electronic density of states, and 〈Z〉 to reveal the evolution of the plasma structure and ionization behavior.

Superfluidity, Bose-Einstein condensation, and structure in one-dimensional Luttinger liquids

We report diffusion Monte Carlo (DMC) and path integral Monte Carlo (PIMC) calculations of the properties of a one-dimensional (1D) Bose quantum fluid. The equation of state, the superfluid fraction ${\ensuremath{\rho}}_{S}/{\ensuremath{\rho}}_{0}$, the one-body density matrix $n(x)$, the pair distribution function $g(x)$, and the static structure factor $S(q)$ are evaluated. The aim is to test Luttinger liquid (LL) predictions for 1D fluids over a wide range of fluid density and LL parameter $K$. The 1D Bose fluid examined is a single chain of $^{4}\mathrm{He}$ atoms confined to a line in the center of a narrow nanopore. The atoms cannot exchange positions in the nanopore, the criterion for 1D. The fluid density is varied from the spinodal density where the 1D liquid is unstable to droplet formation to the density of bulk liquid $^{4}\mathrm{He}$. In this range, $K$ varies from $Kg2$ at low density, where a robust superfluid is predicted, to $Kl0.5$, where fragile 1D superflow and solidlike peaks in $S(q)$ are predicted. For uniform pore walls, the ${\ensuremath{\rho}}_{S}/{\ensuremath{\rho}}_{0}$ scales as predicted by LL theory. The $n(x)$ and $g(x)$ show long range oscillations and decay with $x$ as predicted by LL theory. The amplitude of the oscillations is large at high density (small $K)$ and small at low density (large $K)$. The $K$ values obtained from different properties agree well verifying the internal structure of LL theory. In the presence of disorder, the ${\ensuremath{\rho}}_{S}/{\ensuremath{\rho}}_{0}$ does not scale as predicted by LL theory. A single ${v}_{J}$ parameter in the LL theory that recovers LL scaling was not found. The one body density matrix (OBDM) in disorder is well predicted by LL theory. The ``dynamical'' superfluid fraction, ${\ensuremath{\rho}}_{S}^{D}/{\ensuremath{\rho}}_{0}$, is determined. The physics of the deviation from LL theory in disorder and the ``dynamical'' ${\ensuremath{\rho}}_{S}^{D}/{\ensuremath{\rho}}_{0}$ are discussed.

4He adsorption on a single C40 molecule: Path integral Monte Carlo study

Path-integral Monte Carlo calculations are performed to study the adsorption of 4He atoms on two different isomers of C40: one with the D2 symmetry and the other with the D5d symmetry. Here, we employ the 4He-C40 interaction described by the sum of 4He-C interatomic pair potentials, which allows us to examine the 4He corrugations on the fullerene molecular surfaces. We first observe the layer-by-layer growth of 4He on both C40 isomers, where the first layers are located at a distance of ~ 5.7 A from the molecular centers. From the investigation of the angular 4He density profiles, we find that the 4He layers display different quantum states depending on the number of 4He adatoms N. Between N = 22 and N = 40, the first 4He layer on the D2 isomer shows several different commensurate structures whereas the corresponding 4He layer on D5d undergoes a commensurate-to-incommensurate solid transition like the first 4He layer on graphite. As N increases, the superfluid responses of the first 4He layers on both isomers exhibit reentrant behavior with superfluidity being completely quenched at structurally-ordered states. In addition, the superfluid fraction of the first 4He layer is found to be about three times as large on D2 as on D5d for the same number of 4He atoms. The different structural and superfluid properties of 4He adatoms depending on the isomeric type could provide an experimental way of selecting a specific C40 isomer.