Estimated values of the kinetic energy for liquid $^3$He

TL;DR

This paper estimates the ground-state kinetic energy of liquid $^3$He using a variational path integral quantum Monte Carlo method, resolving previous discrepancies between experimental and theoretical values.

Contribution

It applies an extended fermionic path integral Monte Carlo approach with a fixed-node approximation to accurately compute kinetic energy and other properties of liquid $^3$He.

Findings

Estimated kinetic energy: 10.16±0.05 K/atom.

Results agree with most experimental data.

Method effectively reduces approximation errors.

Abstract

The kinetic energy is estimated for the ground-state of liquid $^3$He at equilibrium density. The obtained value for this quantity, $10.16\pm0.05$ K/atom at density $0.0163~\mbox{\AA}$, is in agreement with most of the experimental data found in the literature. This result resolves a long-standing controversy between experimental and theoretical values of this quantity. The variational path integral method, an "exact" quantum Monte Carlo method extended for fermionic systems, is applied in the calculations. The results obtained are subjected only to the restrictions imposed by a chosen nodal structure without any further approximation, even for quantities that do not commute with the Hamiltonian. The required fixed-node approximation entails an implementation that allows a more effective estimation of the quantities of interest. Total and potential energies together with the radial distribution function are also computed.

Tables

Table 1. Table 1: Kinetic and potential energies in units of K/atom evaluated with the direct ( 𝒟 ) 𝒟 (\cal D) and thermodynamic ( 𝒯 ) 𝒯 (\cal T) estimators using configurations projected from the JS+BF+T wave function.

Estimator	Kinetic Energy	Potential Energy
$𝒟$	$10.19\pm 0.07$	$-12.75\pm 0.01$
$𝒯$	$10.14\pm 0.07$	$-12.75\pm 0.01$

Equations

$$\rho (R,R^{\prime },\beta )=$$

$$\rho (R^{\prime \prime },R^{\prime \prime \prime },\tau )\approx {\rho }_0(R^{\prime \prime },R^{\prime \prime \prime },\tau )e^{-\frac{\tau }{2}\left[V(R^{\prime \prime })+V(R^{\prime \prime \prime })\right]},$$

$$\rho (R_k,R_l,\tau )\approx {\rho }_0(R_k,R_l,\tau )\prod _{i<j}{e}^{-U(r_{ij}^{(k)},r_{ij}^{(l)})},$$

$$U(r^{(k)},r^{(l)})=\frac{\tau }{3}\left[2v(r^{(k)})+v(r^{(l)})\right]+$$

$$\frac{{\tau }^3\lambda }{9}{\left[\frac{\partial v}{\partial r}(r^{(k)})\right]}^2,$$

$$\rho (R,R^{\prime },\beta ){\Psi }_T(R^{\prime })=$$

$$(-1{)}^{n_p}\rho (R,\mathcal{P}{R}^{\prime },\beta ){\Psi }_T(R^{\prime })=$$

$$\rho (R,\mathcal{P}{R}^{\prime },\beta ){\Psi }_T(\mathcal{P}{R}^{\prime }),$$

$$O(\beta )$$

$$=⟨{\Psi }_T|\rho (\beta )\mathcal{O}\rho (\beta )|{\Psi }_T⟩,$$

$$O(\beta )=\int d{R}_1\dots d{R}_{2M+1}P(R_1,\dots ,R_{2M+1})O_L^{\mathcal{X}},$$

$$P(R_1,\dots ,R_{2M+1})$$

$$\times \dots \rho (R_{2M},R_{2M+1},\tau ){\Psi }_T(R_{2M+1}).$$

$$\mathcal{H}=\frac{1}{2m}\sum _{i=1}^N{\mathbf{p}}_i^2+\sum _{i<j}^{N}v(r_{ij}),$$

$${\Psi }_T(R)=e^{-\frac{1}{2}{\sum }_{i<j}u(r_{ij})}\times$$

$${\det }_{\uparrow }\left(e^{i{\mathbf{k}}_l\cdot {\mathbf{r}}_m}\right){\det }_{\downarrow }\left(e^{i{\mathbf{k}}_l\cdot {\mathbf{r}}_n}\right),$$

$$\eta (r)={\lambda }_B{e}^{-{\left(\frac{r-s_B}{w_B}\right)}^2}+\frac{{\lambda }_B^{\prime }}{r^3}.$$

$$\exp \left[-\frac{1}{2}\sum _{i<j}\tilde{u}(|{\mathbf{r}}_i-{\mathbf{r}}_j|)-\frac{{\lambda }_T}{4}\sum _l\mathbf{G}(l)\cdot \mathbf{G}(l)\right],$$

$$\xi (r)=e^{-{\left(\frac{r-s_T}{w_T}\right)}^2}.$$

Full text

We investigated properties of normal condensed 3He at the equilibrium density and compared to experimental values. Neither experimental or theoretical quantities of this system are easily obtained. Direct experimental information about single-particle dynamical properties such as the mean kinetic energy $\langle E_{K}\rangle$ Bryan et al. (2016); Dimeo et al. (1998); Azuah et al. (1995); Mook (1985); Sokol et al. (1985) of this strongly interacting Fermion system can be obtained by deep inelastic neutron scattering. These are challenging experiments since the absorption cross-section for thermal neutrons is about three orders of magnitude higher than for inelastic scattering. On the other hand, theories using quantum Monte Carlo many-body methods must avoid the Fermion sign problem that so far has resisted an entirely satisfactory answer. Most of the experiments report kinetic energies in the range of 8 to 11 K/atom Bryan et al. (2016); Dimeo et al. (1998); Azuah et al. (1995); Mook (1985); Sokol et al. (1985), whereas theory predicted values between 12 and 13 K/atom Mazzanti et al. (2004); Moroni et al. (1997); Manousakis et al. (1983). This is a small, but a significant difference for an “exact” method.

In calculations made at zero temperature, we employed the variational path integral (VPI) method introduced by Ceperley Ceperley (1995), who computed the total energy of 4He at equilibrium density. This is a well established method, also known as path-integral ground-state (PIGS), employed in the recent investigation of a variety of bosonic systems, see for instance references Rossoti et al., 2017; Abolins et al., 2018; Bertaina et al., 2016. We extended the method to deal with fermionic systems, in order to estimate properties of liquid 3He. In this approach, a projection to the ground-state of the system is made from a given initial state using ideas of path-integrals over imaginary time Ceperley (1995). The employed projector and how it is used in the VPI method is reminiscent of how particles are treated in a path-integral Monte Carlo calculation. The necklace describing a particle can be thought of as having been cut and the coordinates at the extremities are assumed to be those of a trial function. This is what we refer to as an open path or polymer. Configurations associated to monomers at the middle of long enough polymers allow one to estimate any quantity, regardless whether their expected values are associated to operators that do or do not commute with the Hamiltonian. “Exact” values are always obtained without the need for any extrapolation. However, since we are dealing with a fermionic system, the usual fixed-node approximation needs to be used. In our context, configurations associated to the trial function at each end of the polymers need to be considered. Results obtained for all quantities of interest are only subjected to the restrictions imposed by a chosen nodal structure.

Our main aim is the investigation of properties of liquid 3He associated with operators that do not commute with the Hamiltonian. We especially want to study the kinetic energy of these systems, since there are controversies between experimental and theoretical results that continue up until the present Bryan et al. (2016). We show that the VPI method gives estimates that are in agreement with most of the experimental results.

Ground-state properties estimated by the VPI method are made by applying the imaginary time evolution operator, $\rho(\beta)=\exp\left(-\beta{\cal H}\right)$, with ${\cal H}$ being the system Hamiltonian, in an initial state $\ket{\Psi_{T}}$ to project out the ground-state $\ket{\phi_{0}}$. The state $\ket{\phi(\beta)}=\rho(\beta)\ket{\Psi_{T}}$ converges exponentially to $\ket{\phi_{0}}$ as $\beta$ increases.

The matrix element $\rho(R,R^{\prime},\beta)=\bra{R}\rho(\beta)\ket{R^{\prime}}$, propagates configuration $R$ to $R^{\prime}$ in a “time” $\beta$ Ceperley (1995), where $R$ stands for all particle coordinates. It is written as the exponential of the action integrated over all paths. The integration can be made by factorizing $\rho(\beta)$ into the product of $M$ projectors $\rho(\tau)$, $\tau=\beta/M$, and using the convolution property

[TABLE]

The intermediary configurations or beads $R_{n}$, $n=1,\ldots,M-1$, can be seen as the set of atomic coordinates at “time” $t=n\tau$. The beads stand for a sort of discretization of the path from $R$ to $R^{\prime}$ in a “time” $\beta$. Therefore the integration of Eq.(Estimated values of the kinetic energy for liquid 3He) converges to the integration over all paths if $\tau$ is small enough. In this case, it is possible to employ the primitive approximation,

[TABLE]

where $V(R)$ is the potential energy of configuration $R$, and $\rho_{0}(R^{\prime\prime},R^{\prime\prime\prime},\tau)$ is the projector of non-interacting atoms, $\rho_{0}(R^{\prime\prime},R^{\prime\prime\prime},\tau)\propto\exp[-(R^{\prime\prime}-R^{\prime\prime\prime})^{2}/4\lambda\tau]$, where $\lambda$ is $\hbar^{2}/2m$. The primitive approximation is accurate to the second order in $\tau$. We also implemented calculations with the Suzuki pair approximation Cuervo et al. (2005); Rossi et al. (2009), which is a fourth order in $\tau$ approximation,

[TABLE]

$r^{(\cdot)}_{ij}$ is the relative distance between atoms $i$ and $j$ within configuration $R_{\cdot}$, if $k$ is even

[TABLE]

and if $k$ is odd then $U(r^{(k)},r^{(l)})=(\tau/3)[v(r^{(k)})+2v(r^{(l)})]$; $v(r)$ is the inter-atomic potential.

By substituting Eq.(2) or Eq.(3) into Eq.(Estimated values of the kinetic energy for liquid 3He) we obtain a formula for $\rho(R,R^{\prime},\beta)$. Any error introduced by one of these approximations can, in general, be made smaller than the statistical uncertainties of the Monte Carlo method. The choice of Eq.(2) or Eq.(3) did not affect our results.

In a system made of identical Fermions such as the one we are interested in, the expression for $\rho(R,R^{\prime},\beta)$ needs to be anti-symmetric under the permutation of any pair of particles in the configuration $R$Ceperley (1995, 1996). However, if $\Psi_{T}(R)$ is anti-symmetric, it is possible to incorporate the minus sign rising from odd permutations in $\rho(R,R^{\prime},\beta)$ into $\Psi_{T}(R)$ since

[TABLE]

where ${\cal P}$ changes the coordinates of $n_{p}$ particles in a given configuration. And so, after integration in $R^{\prime}$, all permutation will have the same result (more details will be given elsewhere).

Any property of the system in its ground-state can be estimated in a straightforward manner. If a given property is associated to an operator ${\cal O}$, its expected value can be written as

[TABLE]

or as

[TABLE]

in terms of the probability distribution function $P$ of a given path

[TABLE]

In Eq. (7), $O_{L}^{{\cal X}}(R_{.})$ is the local value of the operator at a given bead and the index ${\cal X}$ labels different estimators this method can allow us to use. If ${\cal O}$ commutes with the Hamiltonian, by using its coordinate representation it is possible to estimate its value for a given configuration $R$ at the end of the path through $O_{L}^{{\cal E}}(R)={\cal O}\Psi_{T}(R)/\Psi_{T}(R)$. This is the local value of $\cal O$ evaluated for configurations at the end of the path associated to $\Psi_{T}(R)$.

An estimate of the “exact” average value of ${\cal O}$, even if it does not commute with the system Hamiltonian, can be obtained through the so called direct estimator given in the coordinate representation by $O_{L}^{{\cal D}}(R_{i},R_{i+1})={\cal O}\rho(R_{i},R_{i+1},\tau)/\rho(R_{i},R_{i+1},\tau)$, applied at the polymer middle. For efficiency, the best approach is to consider the average value $O_{L}^{{\cal D}}(R_{i},R_{i+1})$ for $i=M$ and $i=M+1$.

For the total and kinetic energy estimates, we can also use the thermodynamic estimators $O^{\cal T}_{L}$ to consider configurations at the middle of the polymer. In this context, derivatives of $\rho(\beta)$ with respect to $\beta$ and the mass $m$ are associated with the total and kinetic energy respectively Ceperley (1995). For any of these estimators, care must be taken when utilizing the Suzuki pair approximation of Eq.(3), since the operators must be inserted in odd beads Rossi et al. (2009).

Since we want to investigate fermionic systems the probability density given by Eq.(8) can be negative. This is the sign problem common to most of the ground-state Monte Carlo methods for fermionic systems. Here we avoid this problem by rejecting sampled paths where $\Psi_{T}(R_{1})\Psi_{T}(R_{2M+1})<0$. This is a fixed-node approximation that has more degrees of freedom than the restriction $\Psi_{T}(R)>0$ imposed when one applies an importance function transformation to sample $\phi_{0}(R)\Psi_{T}(R)$, where $\phi_{0}(R)$ is unknown. We believe that the extra degrees of freedom we have in this instance improves the exploration of the configuration space, especially to regions where the nodal structure of $\Psi_{T}(R)$ is not identical to that of the ground-state.

The system we consider is made of $N$ atoms of 3He inside a cubic box with periodic boundary conditions applied to the faces of the box. In our model, the atoms interact through the well-tested pairwise potential $v(r)$, HFD-B3-FCI1 Aziz et al. (1995), and the Hamiltonian can be written as,

[TABLE]

where $m$ is the 3He mass, ${\mathbf{r}}_{i}$ and ${\mathbf{p}}_{i}$ are respectively the coordinates and the momentum associated to a $i$-th atom and $r_{ij}=|{\mathbf{r}}_{i}-{\mathbf{r}}_{j}|$.

It is interesting to experiment with different trial functions at the end of the polymer. This allows us to investigate the convergence behavior towards the exact ground-state of the quantities of interest. In this way, two wave functions with different degrees of superposition with the ground-state were considered. We performed two series of independent runs, one for each of the functions used at the extremities of the polymer. The simplest function we have considered at the extremities was the Jastrow-Slater (JS) wave function,

[TABLE]

where $u(r)=(b/r)^{5}$. The nodal structure of this wave function was improved by adding backflow correlations in the Slater determinantSchmidt et al. (1981); Panoff and Carlson (1989). These correlations are introduced by a change in the particle coordinates, ${\mathbf{r}}_{.}\rightarrow{\mathbf{r}}_{.}+\sum_{j\neq\cdot}\eta(r_{\cdot j})(\mathbf{r}_{.}-\mathbf{r}_{j})$, of the Slater determinant, where

[TABLE]

and $\lambda_{B},\ s_{B},\ w_{B},\ \lambda_{B}^{\prime}$ are parameters. Three-body correlationsSchmidt et al. (1981); Panoff and Carlson (1989) were also introduced at the extremities of the open path. Its functional form is given by

[TABLE]

where $\mathbf{G}(l)=\sum_{i\neq l}\xi(r_{ij})\mathbf{r}_{ij}$,

[TABLE]

and $s_{T},\ w_{T}$ are parameters. The pseudopotential $\tilde{u}(r)=u(r)-\lambda_{T}\xi^{2}(r)r^{2}$ cancels two-body factors arising from $\mathbf{G}(l)$. We refer to this improved wave function as JS+BF+T. In order for the wave function to be periodic it is required that the correlation functions and its derivatives go smoothly to zero at half of the side of the simulation box, $L$. This can be achieved by the replacement $f(r)\rightarrow f(r)+f(L-r)-2f(L/2)$, where $f$ is either $u$, $\eta$ or $\xi$.

Our calculations were performed with $N=54$ atoms in a non-polarized system at the equilibrium density, $0.0163\ \mbox{\AA}^{-3}$. The sampling of the beads were made by the multi-Metropolis algorithm described in reference Ceperley (1995). The configurations at the extremities of the path were sampled with the usual Metropolis algorithm.

The total energy as a function of $\beta$, $H(\beta)$, was calculated using the estimator at the end of the path for the two different trial wave functions, see Fig. 1. As $\beta$ increases the energy decreases, almost exponentially, creating a sequence of upper bound values to the ground-state energy. The results show that improvements to the trial wave function due to backflow and three-body correlations accelerate the convergence to the ground-state of the system. The ground-state energy itself can only be achieved if the nodal structure of $\Psi_{T}(R)$ is identical to that of the ground-state. In this sense the improvement in the nodal structure of $\Psi_{T}(R)$ is noticeable due to the addition of backflow correlations.

The tail ($\beta\geq 1.5\times 10^{-2}$ K*-1*) of the curve in Fig. 1 associated to the JS+BF+T wave function was fitted to a constant straight line resulting in a total energy of $-2.41\pm 0.01$ K/atom, which is a very good upper bound to the experimental data $-2.47\pm 0.01$ K/atom Roberts et al. (1964). From now on, all results we report are in reference to the results obtained from the wave function above.

For the estimation of the kinetic energy we use the same procedure of considering all the converged values we have obtained for this quantity. The straight line fit to these results gave us the value we adopt for the ground-state kinetic energy, $10.16\pm 0.05$ K/atom. In Fig. 2, we plotted this value together with experimental data from the literature for liquid 3He at equilibrium density. Most of the experimental data lies in a range from $8$ to $11$ K/atom, which is in excellent agreement with our estimates, thus resolving a long-standing disagreement between experimental data and theoretical Monte Carlo calculations For completeness in Table 1, we give values of the potential and kinetic energies calculated with the direct and thermodynamic estimators.

We have also calculated the radial distribution of atoms and its spin-resolved components for atoms with parallel and anti-parallel spins. These quantities were calculated with the direct estimator, see Fig. 3. The anti-parallel spin curve has a more pronounced peak because atoms with different spins do not suffer the Pauli exclusion.

In summary, as we have shown, the VPI approach to study the ground-state of fermionic systems is robust and reliable. Any quantity, associated with operators that do or do not commute with the Hamiltonian is readily estimated without the need of extrapolations. By avoiding them, estimates can be obtained completely free from any possible bias introduced by variational calculations. Moreover, a long-standing disagreement between experimental data and theoretical calculations of the ground-state kinetic energy was resolved. To what degree the findings in this work will be reflected in other systems is still an open question. Nevertheless, this question is very important, because many results for physical properties of great interest were obtained in the literature using extrapolations. Acknowledgments: The authors acknowledge financial support from the Brazilian agencies Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq). Part of the computations were performed at the Centro Nacional de Processamento de Alto Desempenho em São Paulo (CENAPAD-SP).