Path Integral Molecular Dynamics for Bosons

TL;DR

This paper introduces a recursive path integral molecular dynamics method for simulating Bosons, efficiently handling exchange symmetry and accurately modeling Bosons in a 2D trap.

Contribution

A novel recursive approach for including exchange symmetry in path integral molecular dynamics for Bosons, avoiding explicit permutation enumeration.

Findings

Method scales cubically with system size

Results agree with exact solutions for 2D trap Bosons

Analysis of exchange effects at different temperatures

Abstract

Trapped Bosons exhibit fundamental physical phenomena and are potentially useful for quantum technologies. We present a method for simulating Bosons using path integral molecular dynamics. A main challenge for simulations is including all permutations due to exchange symmetry. We show that evaluation of the potential can be done recursively, avoiding explicit enumeration of permutations, and scales cubically with system size. The method is applied to Bosons in a 2D trap and agrees with essentially exact results. An analysis of the role of exchange with decreasing temperature is also presented.

Equations

$$Z_D^{(N)}\sim \int e^{-\beta \left(V_D^{(N)}+U\right)}d{R}_1...d{R}_N,$$

$$V_D^{(N)}(R_1,...,R_N)=\sum _{l=1}^N{E}_o(R_l)$$

$$\equiv \frac{1}{2}m{\omega }_P^2\sum _{l=1}^N\sum _{j=1}^P{\left({\mathbf{r}}_l^{j+1}-{\mathbf{r}}_l^j\right)}^2.$$

$$Z_B^{(3)}$$

$$e^{-\beta V_B^{(N=3)}}=\frac{1}{6}\left(e^{-\beta E_{ooo}}+3e^{-\beta E_{Oo}}+2e^{-\beta E_{◯}}\right).$$

$$E_N^{(k)}(R_{N-k+1},...,R_N)\equiv \frac{1}{2}m{\omega }_P^2\sum _{l=N-k+1}^N\sum _{j=1}^P{\left({\mathbf{r}}_l^{j+1}-{\mathbf{r}}_l^j\right)}^2.$$

$$e^{-\beta V_B^{(N)}}=\frac{1}{N}\sum _{k=1}^N{e}^{-\beta \left(E_N^{(k)}+V_B^{(N-k)}\right)},$$

$$Z_B^{(N)}\sim \int e^{-\beta \left(V_B^{(N)}+U\right)}d{R}_1...d{R}_N.$$

$$\rho (\mathbf{r})=⟨\frac{1}{P}\sum _{j=1}^P\sum _{l=1}^N\delta ({\mathbf{r}}_l^j-\mathbf{r})⟩.$$

$$⟨E⟩=\frac{PdN}{2\beta }+⟨U⟩+⟨V_B^{(N)}+\beta \frac{\partial V_B^{(N)}}{\partial \beta }⟩,$$

$$V(|{\mathbf{r}}_l-{\mathbf{r}}_m|)=\frac{g}{\pi s^2}{e}^{-\frac{({\mathbf{r}}_l-{\mathbf{r}}_m{)}^2}{s^2}}.$$

$$P_c=\frac{w_c\cdot e^{-\beta E_c}}{\frac{1}{N}{\sum }_{k=1}^N{e}^{-\beta \left(E_N^{(k)}+V_B^{(N-k)}\right)}},$$

Full text

††preprint: APS/123-QED

Trapped cold atoms are fascinating systems that allow studying fundamental physical phenomena. These range from the observations of Bose-Einstein condensation in macroscopic systems Anderson et al. (1995); Davis et al. (1995) to the study of exotic supersolids Tanzi et al. (2019) and the formation of ultra-cold diatomic molecules Regal et al. (2003); Ni et al. (2008). Experiments can control many properties of such systems. They can be studied in optical lattices with single-atom resolution Bakr et al. (2010); Sherson et al. (2010), at varying interaction strengths Paredes et al. (2004); Kinoshita et al. (2004); Zürn et al. (2012), under the influence of external fields Büchler et al. (2007) and more. Developing numerical simulations with microscopic resolution of such systems is thus highly desirable. In this Letter, we report a new method for performing simulations of Bosons using path integral molecular dynamics (PIMD) Parrinello and Rahman (1984).

Symmetry under permutation of indistinguishable particles is a fundamental property of many-body quantum systems. Therefore, methods which include exchange effects are an important goal in physics and chemistry. The path integral formulation of quantum mechanics Feynman and Hibbs (2005) is a powerful tool for studying many-body systems Markland and Ceriotti (2018). Applications of path integral methods often neglect exchange effects, treating the particles as distinguishable. This approximation leads to the well-known isomorphism between the partition function of a quantum-mechanical system and that of a classical system, in which each particle is represented by a ring-polymer composed of $P$ beads connected through harmonic springs Chandler and Wolynes (1981). The partition function of the classical system can be sampled using molecular dynamics (MD) Parrinello and Rahman (1984) or Monte Carlo (MC) Pollock and Ceperley (1984).

In the case of indistinguishable particles, one must include all possible $N!$ permutations  Ceperley (1995). This leads to ring-polymer configurations in which permuted particles are connected sequentially into longer rings (see Figure 1). For Bosons the weights of all configurations are positive while for Fermions they can be negative, giving rise to the infamous sign problem DuBois et al. (2014); Runeson et al. (2018). Here we focus on Bosons, the treatment of Fermions being beyond the scope of this Letter. Enumerating all Bosonic ring-polymer configurations is impractical in all but the smallest systems. To alleviate this difficulty, one can address it as a sampling problem and develop Monte Carlo algorithms that introduce permutations. The most notable examples are the pioneering applications of path integral Monte Carlo (PIMC) to superfluid Helium Ceperley (1995). Significant advancement in efficiency was achieved with the introduction of the worm algorithm Boninsegni et al. (2006).

Far less progress has been made in the use of MD to sample the Bosonic partition function. Miura and Okazaki have suggested a method for performing PIMD simulations of Bosons Miura and Okazaki (2000). However, their approach involved calculating a permanent. Since this is equivalent to enumerating all possible permutations, the method was only applied to very small systems. In this Letter, we show how to evaluate the potential energy and forces acting on the particles avoiding calculating the permanent. This results in an algorithm that scales cubically with system size, allowing larger systems to be studied. The resulting approach is simpler than available rather-complex sampling algorithms. Developing alternatives to MC sampling can extend the applicability of path integral methods.

In the following, we first present the method and benchmark it against analytical results for up to 64 non-interacting Bosons and numerical diagonalization of the Hamiltonian for small interacting systems Mujal et al. (2017). We also apply the method to 32 interacting particles in a 2D trap, for which the number of permutations to be considered is already prohibitive. Finally, we present an approach to analyze the importance of exchange effects from the different ring-polymer configurations.

We consider a system of $N$ identical particles of mass $m$ at inverse temperature $\beta=(k_{B}T)^{-1}$. The interaction between the particles is denoted as $V(\bm{r}_{1},...,\bm{r}_{N})$.

Neglecting effects due to quantum symmetry, the path integral representation of the partition function is given Tuckerman (2010) by

[TABLE]

where $R_{l}$ represents collectively the coordinates $(\bm{r}_{l}^{1},...,\bm{r}_{l}^{P})$ of the $P$ beads composing the ring polymer corresponding to particle $l$ and $dR_{l}=d\bm{r}_{l}^{1}...d\bm{r}_{l}^{P}$. The exact result is obtained in the limit $P\rightarrow\infty$. In Equation 1, the interaction potential between beads $j$ of different particles is scaled by the number of beads and $U(R_{1},...,R_{N})=\frac{1}{P}\sum_{j=1}^{P}V(\bm{r}_{1}^{j},...,\bm{r}_{N}^{j})$. The potential $V_{D}^{(N)}(R_{1},...,R_{N})$ is the sum of ring polymer spring energies $E_{o}(R_{l})$ for each particle,

[TABLE]

The bead index in Equation 2 is cyclic, $\bm{r}_{l}^{P+1}=\bm{r}_{l}^{1}$, and $\omega_{P}=\sqrt{P}/\beta\hbar$ is the frequency of the springs connecting the beads in the ring polymer.

To include exchange effects one must consider all permutations. This results in a sum over $N!$ terms in the partition function Miura and Okazaki (2000). Each term in the sum corresponds to a configuration of the ring polymers in which some of the particles are connected into longer rings Runeson et al. (2018). Figure 1 shows all configurations for $N=3$. However, when integrated over all particles, some of the configurations contribute equally to the partition function. Therefore, one can sum over a much smaller number of terms arising from the number of unique configurations for $N$ particles Lyubartsev and Vorontsov-Velyaminov (1993). Each term is weighted by the number of different permutations which lead to the same configuration. This observation has been used recently in attempts to overcome the Fermion sign problem in PIMC simulations Voznesenskiy et al. (2009); DuBois et al. (2014). For Bosons, all configurations give a positive contribution and no such sign problem arises.

For example, the partition function for $N=3$ Bosons can be written Tuckerman (2010) as

[TABLE]

where the potential $V_{B}^{(N=3)}(R_{1},R_{2},R_{3})$ is given by

[TABLE]

Here we define $E_{ooo}\equiv V_{D}^{(3)}$, $E_{Oo}\equiv E_{o}(R_{1})+E^{(2)}(R_{2},R_{3})$ and $E_{\bigcirc}\equiv E^{(3)}(R_{1},R_{2},R_{3})$ and denote the total spring energy of a ring polymer constructed by connecting all of the beads of $k$ particles sequentially as

[TABLE]

In Equation 5 it is implied that $\bm{r}_{N}^{P+1}=\bm{r}_{1}^{1}$ or otherwise $\bm{r}_{l}^{P+1}=\bm{r}_{l+1}^{1}$. We note that $U(R_{1},...,R_{N})$ is identical to the case of distinguishable particles and that, using the Trotter decomposition, only the terms arising from the kinetic energy operator are affected by the permutations.

While the computational cost is significantly reduced, as compared to summing over all permutations, the number of all the unique configurations still scales exponentially with system size. However, taking $N=3$ as an example, one finds that all configurations can be generated by the following procedure (see Figure 2): Starting from a single ring polymer for particle 1, we generate all configurations for $N=2$ by adding a ring polymer representing particle 2, to the configuration of $N=1$ and then adding another configuration composed of particles $1-2$ connected sequentially to a long ring. We generate the configurations for $N=3$ by adding a ring of double-length, composed of particles $2-3$ connected sequentially, to the configuration of $N=1$ and a ring polymer representing particle 3 to the configurations of $N=2$. Lastly, we add another configuration composed of all 3 particles connected sequentially in a long ring.

This procedure can be generalized to $N$ particles using the recurrence relation

[TABLE]

where $V_{B}^{(0)}=0$. The partition function is then given by

[TABLE]

The forces required for PIMD simulations, can also be written using a recurrence relation that is given in the SI. We note that in practical evaluation of Equation 6 it is useful to rewrite it so that only small energy differences appear in the exponent. Further details of the implementation are given in the SI.

We note that for non-interacting particles, the integration in Equation 7 can be done analytically. Using Equation 6, this results in a known Borrmann and Franke (1993); Schmidt and Schnack (2002); Krauth (2006) recurrence relation for the partition function of $N$ non-interacting Bosons, $Z_{B}^{(N)}=\sum_{k=1}^{N}z_{k}Z_{B}^{(N-k)}$, where $z_{k}$ is the partition function of a single particle evaluated at inverse temperature $k\beta$.

A key observation is that Equation 6 can be used to evaluate the potential energy for $N$ Bosons, at a particular point $(R^{\prime}_{1},...,R^{\prime}_{N})$ without enumerating explicitly all ring-polymer configurations. This leads to an algorithm that scales cubically with system size (see SI).

To evaluate thermal expectation values for quantum-mechanical operators using PIMD simulations, expressions for estimators must be derived. Estimators for operators which are local in coordinate representation are obtained as in PIMD simulations for distinguishable particles. For example, the estimator for the density is

[TABLE]

For the energy, we derive

[TABLE]

where $V_{B}^{(N)}+\beta\frac{\partial V_{B}^{(N)}}{\partial{\beta}}$ is also calculated using a recurrence relation, given in the SI. It is well known that, for distinguishable particles, the variance of the thermodynamic estimator grows with the number of beads $P$ Herman et al. (1982). Instead, the estimator above does is well-behaved due to the fact that all terms carry a positive sign.

To test the method, we first applied it to systems of up to $N=64$ non-interacting Bosons. The agreement for the ground-state energy was excellent $(\sim 0.5\%)$ as shown in the SI. To verify that the correct Bose-Einstein statistics are obtained we also checked that the energy as a function of temperature agrees with analytical results (see SI). To further validate the method, we performed simulations of small systems of interacting Bosons. For such small systems very accurate numerical results, obtained from diagonalizing the Hamiltonian, are available Mujal et al. (2017). Following Mujal et al., we simulate $N=2-4$ Bosons confined in an isotropic 2D harmonic trap and interacting through the repulsive Gaussian pair-potential:

[TABLE]

In Equation 10, $g$ and $s$ are a measure of the interaction strength and range, respectively. We take $s=0.5$ and we vary the interaction strength in the interval $g=0-16$ (in harmonic oscillator units). These values have been shown to span a wide range of interactions from weakly interacting particles, for which a mean-field approximation Dalfovo et al. (1999) is valid, to strongly interacting particles, for which a correlated ansatz for the wave function must be used Mujal et al. (2017). Moreover, for small systems it was shown that the high values of $g$ lead to a pair correlation function which resembles that of a system of non-interacting Fermions of the same size Mujal et al. (2018).

The new method was implemented in a development version of LAMMPS Plimpton (1995). For all the results below simulations using 72 beads were performed with a time step of 1 fs. Statistics was accumulated for 14 ns following an initial equilibration run of 1 ns. The statistical error was evaluated using block averaging. If not shown, it is smaller than the symbol size. A 2D trap frequency corresponding to $\hbar\omega_{0}=3$ meV was used. The simulations were performed with $m=1$ a.u. and at a temperature corresponding to $\beta\hbar\omega_{0}=6$ using the Nosé-Hoover chains thermostat Martyna et al. (1992). The results were converged with respect to the number of beads at all range of interaction strengths (see SI for details).

Figure 3 shows the energy as a function of $g$ for $N=2-4$. Very good agreement is obtained with the results of ref. Mujal et al. (2017). The maximum absolute deviation is  2.9% and the mean absolute deviation is 1.3%. The comparison between the density obtained and exact diagonalization results is also shown. The agreement is very good and the mean absolute deviation for all range of interactions is $\sim 2.7\%$. This slight discrepancy could be due to the finite basis-set used in ref. Mujal et al. (2017).

As mentioned previously, the number of ring-polymer configurations scales exponentially with $N$. For $N=32$ the number of unique configurations is 8439 while the number of energy terms required to evaluate Equation 6 is only 528 (see SI for details). Therefore, we chose $N=32$ to demonstrate the ability of the method to handle a large number of permutations. We used $g=3$ since already for such a repulsive potential significant differences from the non-interacting case are observed. The simulations were performed with $P=36$. Other simulation parameters are identical to the ones used for $N=2-4$. Figure 4 presents the density and pair-correlation function. Due to the repulsion between the particles, the density at the center of the trap is significantly lowered and the probability to find particles away from the center of the trap is increased. The pair-correlation function shows a complementary trend. The probability of finding any two particles at the same position is lowered as the repulsive interaction is introduced and the most probable distance between the two particles is $r\sqrt{m\omega_{0}/\hbar}\approx 0.9$ for $g=3$.

For superfluid Helium, Ceperley Ceperley (1995) has shown that the probability of a particle to be incorporated into a long ring, composed of several particles, grows significantly below the critical temperature. Therefore, the relative probability of ring-polymer configurations of different lengths can serve as an indication of the importance of exchange effects. Moreover, it was shown that the probability for the longest ring tends to $1/N$ at low temperatures. This analysis is readily performed in PIMC simulations, in which permutations are individually introduced. However, it is desirable also for PIMD simulations. To this effect, we define the instantaneous probability of a given ring-polymer configuration as

[TABLE]

where $w_{c}$ is the number of permutations pertaining to configuration $c$. The average probability for $N=2$ and $32$ is presented in Figure 5. It can be seen that the probability for the ring-polymer configuration in which all $N$ particles are connected in a ring grows with decreasing temperature. Moreover, for $g=0$ it has the correct asymptotic behaviour already at $\beta\hbar\omega_{0}=6$. Introducing repulsive interactions, the increase in probability is slower, which suggests that exchange effects become more important at temperatures lower than those of a non-interacting system. The probability for the configuration representing distinguishable particles decreases with decreasing temperature (not shown). For $N=32$ it is negligible, compared to the longest ring-polymer configuration, at all temperatures studied.

In conclusion, the aim of this paper is to present a new method for performing finite-temperature simulations of Bosons. The method has been tested and applied to systems of trapped Bosons with up to $64$ non-interacting and $32$ interacting particles. This represents an order of magnitude increase in system size, compared to diagonalization of the Hamiltonian Mujal et al. (2017) and previous PIMD simulations of Bosons Miura and Okazaki (2000). Our approach is not limited to systems of this size and we have already been experimenting with much larger systems. Since the computational cost scales relatively modestly with the number of particles, we hope that it will find useful applications in the study of cold atoms. For example, one could imagine using this approach in combination with ring-polymer MD Habershon et al. (2013) to study Boson dynamics.