Path Integral Monte Carlo Simulations Research Articles

Neutron diffraction determination of the microscopic structure of solid deuterium close to melting

The static structure factor of solid normal deuterium, close to its triple point, has been measured using time-of-flight neutron diffraction. Experimental results have been compared to path-integral Monte Carlo simulations as well as to a simple uncorrelated lattice calculation including the mean-squared displacement as an adjustable parameter. We found a very good agreement between experimental data and quantum simulations in all the measured range, while the uncorrelated model, as expected, turned out to be satisfactory only at long distances and for a mean-squared displacement value almost twice as large as the low-temperature one. This fact can be ascribed to the large anharmonic character of the ${\text{D}}_{2}$ quantum crystal close to its melting.

Ordered structure formation in 2D mass asymmetric electron–hole plasmas

We study strong Coulomb correlations in dense two-dimensional electron–hole plasmas by means of direct path integral Monte Carlo simulations. In particular, the formation and dissociation of bound states, such as excitons, bi-excitons and many particle clusters, is analyzed and the density-temperature regions of their occurrence are identified. At high density, the Mott transition to the fully ionized state (electron–hole hexatic liquid) is detected. Particular attention is paid to the influence of the hole to electron mass ratio M on the properties of the plasma. For high enough values of M we observed the formation of Coulomb hole crystal-like structures.

Critical Temperature of Interacting Bose Gases in Two and Three Dimensions

We calculate the superfluid transition temperature of homogeneous interacting Bose gases in three and two spatial dimensions using large-scale path integral Monte Carlo simulations (with up to N=10;{5} particles). In 3D we investigate the limits of the universal critical behavior in terms of the scattering length alone by using different models for the interatomic potential. We find that this type of universality sets in at small values of the gas parameter na3 < or approximately 10(-4). This value is different from the estimate na3 < or approximately 10(-6) for the validity of the asymptotic expansion in the limit of vanishing na3. In 2D we study the Berezinskii-Kosterlitz-Thouless transition of a gas with hard-core interactions. For this system we find good agreement with the classical lattice |psi|4 model up to very large densities. We also explain the origin of the existing discrepancy between previous studies of the same problem.

Survey of classical density functionals for modelling hydrogen physisorption at77K

This work surveys techniques based on classical density functionals for modeling the quantum dispersion of physisorbed hydrogen at 77K. Two such techniques are examined in detail. The first is based on the "open ring approximation" (ORA) of Broukhno et al., and it is compared with a technique based on the semiclassical approximation of Feynman and Hibbs (FH). For both techniques, a standard classical density functional is used to model hydrogen molecule-hydrogen molecule (i.e., excess) interactions. The three-dimensional (3D) quantum harmonic oscillator (QHO) system and a model of molecular hydrogen adsorption into a graphitic slit pore at 77K are used as benchmarks. Density functional results are compared with path-integral Monte Carlo simulations and with exact solutions for the 3D QHO system. It is found that neither of the density functional treatments are entirely satisfactory. However, for hydrogen physisorption studies at 77K the ORA based technique is generally superior to the FH based technique due to a fortunate cancellation of errors in the density functionals used. But, if more accurate excess functionals are used, the FH technique would be superior.

Anharmonic Effects in Neutron Scattering Studies of Lattice Excitations in BCC 4He

We report high resolution measurements of lattice excitations of bcc solid 4He using inelastic neutron scattering. The resulting dispersion relations of the phonons agree well with Path Integral Monte Carlo Simulations which we carried out in parallel with the experimental work. In addition, we studied the recently discovered two “optic-like” excitation branches. One of these branches is weakly dispersive, while the second one seems dispersionless. The dispersive optical branch interacts with the phonons. We suggest that this interaction may be responsible for some of the anharmonic effects reported in the past, such as phonon interference. We argue that the agreement of experimental results with the “Self Consistent Phonon Theory” may be improved if the coupling between the phonons and one of the “optic-like” excitation branches is included in the calculations.

Electron wave packets: Quantum statistics in path integral representation

Fundamental relations between ensemble-averaged kinetic energy, pressure, and virial are obtained for quantum many-particle systems at finite temperatures. A path-integral method is developed for calculating kinetic energy for mixed (quantum-classical) systems. The method eliminates the problem of divergent variance. Path-integral Monte Carlo simulations are performed to validate the method as applied to single-electron and electron-pair wave packets, with exact treatment of exchange symmetry and spin states. Quantitative results are compared with corresponding characteristics calculated for point charges.

Solvation Structure and Rotational Dynamics of LiH in 4He Clusters

We present results of path integral Monte Carlo simulations of LiH solvated in superfluid 4He clusters of size up to N = 100. Despite the light mass of LiH and the strongly anisotropic LiH-He potential with a large repulsion at the hydrogen end, LiH is solvated inside the cluster for sufficiently large N. Using path integral correlation function analysis, we have determined the dipole (J = 1) rotational excitations of the cluster and a corresponding effective rotational constant Beff of the solvated LiH. We predict that Beff is greatly reduced with respect to the gas-phase rotational constant B, to a value of only about 6% of B. This exceptionally large reduction of the rotational constant is due to the highly anisotropic 4He solvation structure around LiH. It does not follow the previously established trend of a relatively small B reduction for light molecules, showing the strongest reduction of all molecules in 4He to date. Comparison of the calculated rotational spectra of LiH in helium obeying Bose and Boltzmann statistics, respectively, demonstrates that the Bose statistics of helium is an essential requirement for obtaining well-defined molecule rotational spectra in helium-4.

Activity expansion calculation of shock-compressed helium: The liquid Hugoniot

The Hugoniot of shock-compressed liquid helium was calculated using the activity expansion dense plasma program. The predicted maximum compression is at $100\phantom{\rule{0.3em}{0ex}}\mathrm{GPa}$ and $58\phantom{\rule{0.2em}{0ex}}000\phantom{\rule{0.3em}{0ex}}\mathrm{K}$, at a density nearly six-fold that of the initial state. Comparisons are made with recent path-integral Monte Carlo simulations, which predict a smaller maximum compression of 5.24-fold, near $360\phantom{\rule{0.3em}{0ex}}\mathrm{GPa}$ and $150\phantom{\rule{0.2em}{0ex}}000\phantom{\rule{0.3em}{0ex}}\mathrm{K}$.

Rovibrational Spectra for the HCCCN·HCN and HCN·HCCCN Binary Complexes in 4He Droplets

Rovibrational spectra are measured for the HCCCN*HCN and HCN*HCCCN binary complexes in helium droplets at low temperature. Though no Q-branch is observed in the infrared spectrum of the linear HCN*HCCCN dimer, which is consistent with previous experimental results obtained for other linear molecules, a prominent Q-branch is found in the corresponding infrared spectrum of the HCCCN*HCN complex. This Q-branch, which is reminiscent of the spectrum of a parallel band of a prolate symmetric top, implies that some component of the total angular momentum is parallel to the molecular axis. The appearance of this particular spectroscopic feature is analyzed here in terms of a nonsuperfluid helium density induced by the molecular interactions. Finite temperature path integral Monte Carlo simulations are performed using potential energy surfaces calculated with second-order Möller-Plesset perturbation theory, to investigate the structural and superfluid properties of both HCCCN*HCN(4He)N and HCN*HCCCN(4He)N clusters with N < or = 200. Explicit calculation of local and global nonsuperfluid densities demonstrates that this difference in the rovibrational spectra of the HCCCN*HCN and HCN*HCCCN binary complexes in helium can be accounted for by local differences in the superfluid response to rotations about the molecular axis, i.e., different parallel nonsuperfluid densities. The parallel and perpendicular nonsuperfluid densities are found to be correlated with the locations and strengths of extrema in the dimer interaction potentials with helium, differences between which derive from the variable extent of polarization of the CN bond in cyanoacetylene and the hydrogen-bonded CH unit in the two isomers. Calculation of the corresponding helium moments of inertia and effective rotational constants of the binary complexes yields overall good agreement with the experimental values.

Interacting electrons in one dimension: a path integral Monte Carlo study

Path integral Monte Carlo simulations have been carried out for a system of interacting fermionic particles in one dimension. Due to the fact that the sign problem can be completely eliminated for one-dimensional systems, such simulations allow one to obtain accurate energies and particle densities in arbitrary external potentials. Two cases were considered: electrons in the field of a uniform neutralizing background and in an effective field of atoms on a lattice. In the latter case, a clear asymmetry of the charge density distribution of conducting electrons and holes has been observed.

USING DISCRETE VARIABLE REPRESENTATION PATH INTEGRAL MONTE CARLO WITH METROPOLIS SAMPLING TO COMPUTE GROUND STATE WAVEFUNCTIONS

The discrete variable representation (DVR) matrix dynamics formulation of the path integral Monte Carlo (PIMC) method, implemented numerically in a way that enables Metropolis sampling to be employed, is proposed as a means of computing ground state quantum wavefunctions. A key advantage of the DVR-PIMC approach is that customized marginal potentials may be employed, leading to significantly larger PIMC time step sizes, and substantial reductions in computational (CPU) effort. An additional key advantage of the present implementation is that the DVR provides a natural set of interpolant functions that can be used for accurate interpolation and extrapolation of function and tensor quantities away from predefined grid points. The new method is applied here to compute the ground state wavefunction of a model one degree-of-freedom (1 DOF) Morse oscillator system. A one-to-two order-of-magnitude reduction in CPU effort is observed, in comparison with a conventional PIMC simulation. The generalization for many DOFs is straightforward, and expected to result in even greater performance enhancement.

QUANTUM BOLTZMANN LIQUIDS

We study homogeneous normal systems of bosons under the influence of interparticle forces with a strongly repulsive component at short relative particle-particle distances. The repulsion prevents short-ranged exchange between the bosonic constituents in the quantum fluid. Consequently, the bosons remain distinguishable at temperatures far below the classical high-temperature regime. At these low temperatures such fluids and liquids display nevertheless distinct quantum effects due to quantum-mechanical phase-phase correlations. Typical examples are liquid para-hydrogen and fluid 4 He under certain thermodynamic conditions. The study employs Correlated Density-Matrix theory and Path-Integral Monte-Carlo simulations.

Temperature Effects on Electron Correlations in Two Coupled Quantum Dots

The path-integral Monte Carlo simulation method is used to examine one and two electrons in a system of two coupled disc-like quantum dots (QD) in a zero magnetic field. With this approach we are able to evaluate the one-electron distributions and two-electron correlation functions, and finite temperature effects on both. Increase of temperature broadens the distributions as expected, the effect being smaller for correlated electrons than for single ones. The simulated one- and two-particle distributions of a single and two coupled quantum dots are also compared to those from other theoretical methods. For the one-particle distributions we find a good agreement with those from the DFT approach. The effect of the third dimension or the thickness of the almost two-dimensional disc-like QDs is small for the one-particle distributions, but it is clearly seen in the electron-electron correlation or the two-particle distribution function at low temperatures. The mutual Coulomb energy of the two electrons is found to be temperature-independent, and also, independent of the correlation effects on the dynamics. Computational capacity is found to become the limiting factor in simulations with increasing accuracy or increasing number of particles, and in case of fermions in particular. This and other aspects of PIMC and its capability for this type of calculations are also discussed.

Screening of particle exchange in quantum Boltzmann liquids

We employ correlated density-matrix theory of strongly correlated Bose fluids to analyze the structural properties of quantum Boltzmann liquids. The constituents of such a normal quantum system are distinguishable as in a classical fluid since the interparticle forces prevent any exchange of identical bosons at short relative distances. Our study focuses on this particle-screening effect and on its consequences for the properties of various correlation functions, structure functions, momentum distributions, and quasiparticle and collective excitations. The formalism of the adopted microscopic theory is applied to a detailed numerical investigation of particle-screening properties and the quantum behavior of liquid para-hydrogen close to the triple point temperature. The theoretical results are compared with numerical data of path-integral Monte Carlo simulations and with available experimental results of recent cross-section measurements by neutron scattering.

Computational study of the melting-freezing transition in the quantum hard-sphere system for intermediate densities. I. Thermodynamic results

The points where the fluid-solid (face-centered-cubic) transition takes place in the quantum hard-sphere system, for reduced densities 0.85>rhoN*>0.5 (reduced de Broglie wavelengths lambdaB*<or=0.8), have been determined via calculations of Helmholtz free energies. A number of complementary methods have been utilized, namely, path-integral Monte Carlo simulations for fixing the basic thermodynamic and structural quantities, Ornstein-Zernike computations of the fluid isothermal compressibilities using the centroid correlations, and applications of the Einstein crystal technique. Attention is paid to the evaluation of the statistical uncertainties in the isothermal compressibilities and also to the quantum implementation of the Einstein crystal technique by including explicitly the constraint of fixed center of mass. The equation of state along the fluid lambdaB* branches studied has been determined with two methods, one based on the isothermal compressibilities and the other on the usual virial estimator. Along the solid lambdaB* branches the equation of state has been fixed with the virial estimator. The results indicate that the phase transition investigated is governed by entropic effects and that the fluid-solid coexistence densities are arranged along a straight line rhoFCC*=rho(rhoF*), a behavior which at least holds even for lambdaB*<2, as revealed by completing the present analysis with data available in the literature.

Computational study of the melting-freezing transition in the quantum hard-sphere system for intermediate densities. II. Structural features

The structural features of the quantum hard-sphere system in the region of the fluid-face-centered-cubic-solid transition, for reduced number densities 0.45<rhoN*<or=0.9 (reduced de Broglie wavelengths lambdaB*<or=0.8), are presented. The parameters obtained with path-integral Monte Carlo simulations for the fluid, amorphous, and solid phases are related to the distinct sorts of pair correlations that can be defined in a path-integral quantum fluid (instantaneous, continuous linear response and centroids). These parameters cover the pair radial correlation functions, the configurational structure factors, the order parameters Q4 and Q6, and the radii of gyration of the path-integral necklaces. Also, the fluid static structure factors have been computed by solving appropriate Ornstein-Zernike equations. A number of significant regularities in the above parameters involving both sides of the crystallization line are reported, and a comparison with results for Lennard-Jones quantum systems that can be found in the literature is made. On the other hand, the main amplitudes of the quantum fluid structure factors follow a complex behavior along the crystallization line, which points to difficulties in identifying a neat rule, similar to that of Hansen-Verlet for classical fluids, for these quantum amplitudes. To complete this study a further analysis of the instantaneous and centroid triplet correlations in the vicinities of the fluid-face-centered-cubic-solid phase transition of hard spheres has been performed, and some interesting differences between the classical and quantum melting-freezing transition are observed.

Structure, superfluidity, and quantum melting of hydrogen clusters

We present results of a theoretical study of para-hydrogen and ortho-deuterium clusters at low temperature (0.5 K < T < 3.5 K), based on Path Integral Monte Carlo simulations. Clusters of size up to N=21 para-hydrogen molecules are nearly entirely superfluid at T < 1 K. For 21 < N < 30, the superfluid response displays strong variations with N, reflecting structural changes that occur on adding or removing even a single molecule. Some clusters in this size range display quantum melting, going from solid- to liquid-like as T tends to 0. Melting is caused by quantum exchanges of molecules. The largest para-hydrogen cluster for which a significant superfluid response is observed comprises 27 molecules. Evidence of a finite superfluid response is presented for ortho-deuterium clusters of size up to 14 molecules. Magic numbers are observed, at which both types of clusters feature pronounced stability.

Correlation effects in partially ionized mass asymmetric electron-hole plasmas

The effects of strong Coulomb correlations in dense three-dimensional electron-hole plasmas are studied by means of unbiased direct path integral Monte Carlo simulations. The formation and dissociation of bound states, such as excitons and biexcitons, is analyzed and the density-temperature region of their appearance is identified. At high density, the Mott transition to the fully ionized metallic state (electron-hole liquid) is detected. Particular attention is paid to the influence of the hole to electron mass ratio M on the properties of the plasma. Above a critical value of about M=80 formation of a hole Coulomb crystal was recently verified [Bonitz, Phys. Rev. Lett. 95, 235006 (2005)] which is supported by additional results. Results are related to the excitonic phase diagram of intermediate valent Tm[Se,Te], where large values of M have been observed experimentally.

Influence of nuclear fluctuations on the NMR parameters of bullvalene: A Feynman path integral – Ab initio study

The isotropic magnetic shieldings and chemical shifts in bullvalene have been determined by Feynman path integral Monte Carlo simulations which have been combined with the gauge-including atomic orbital formalism. This setup renders possible the treatment of the quantum and thermal nuclear degrees of freedom in the calculation of NMR quantities. Ab initio Hartree–Fock Hamiltonians have been employed in the present ensemble averaging of NMR parameters. The vibrational motions cause a nuclear deshielding which depends on the topology of the centre studied. The consideration of these corrections leads to a better agreement between measured and calculated chemical shifts.

First Principles Calculations of Shock Compressed Fluid Helium

The properties of hot dense helium at megabar pressures are studied with two first principles computer simulation techniques: path integral Monte Carlo simulation and density functional molecular dynamics. The simulations predict that the compressibility of helium is substantially increased by electronic excitations that are present in the hot fluid at thermodynamic equilibrium. A maximum compression ratio of 5.24(4)-fold the initial density was predicted for 360 GPa and 150,000 K. This result distinguishes helium from deuterium, for which simulations predicted a maximum compression ratio of 4.3(1). Hugoniot curves for statically precompressed samples are also discussed.