Path integral Monte Carlo ground state approach: Formalism, implementation, and applications

Path integral Monte Carlo ground state approach: Formalism, implementation, and applications

Yangqian Yan and D. Blume Department of Physics, Indiana University Purdue University Indianapolis (IUPUI), Indianapolis, Indiana 46202, USA Department of Physics and Astronomy, Washington State University, Pullman, Washington 99164-2814, USA Homer L. Dodge Department of Physics and Astronomy, The University of Oklahoma, 440 W. Brooks Street, Norman, OK 73019, USA
July 12, 2019

Monte Carlo techniques have played an important role in understanding strongly-correlated systems across many areas of physics, covering a wide range of energy and length scales. Among the many Monte Carlo methods applicable to quantum mechanical systems, the path integral Monte Carlo approach with its variants has been employed widely. Since semi-classical or classical approaches will not be discussed in this review, path integral based approaches can for our purposes be divided into two categories: approaches applicable to quantum mechanical systems at zero temperature and approaches applicable to quantum mechanical systems at finite temperature. While these two approaches are related to each other, the underlying formulation and aspects of the algorithm differ. This paper reviews the path integral Monte Carlo ground state (PIGS) approach, which solves the time-independent Schrödinger equation. Specifically, the PIGS approach allows for the determination of expectation values with respect to eigen states of the few- or many-body Schrödinger equation provided the system Hamiltonian is known. The theoretical framework behind the PIGS algorithm, implementation details, and sample applications for fermionic systems are presented.

1 Introduction

Monte Carlo techniques have found many applications, ranging from the modeling of the stock market to the simulation of classical and quantum spin models [1]. This review introduces the path integral Monte Carlo ground state (PIGS) method [2, 3, 4, 5], which allows for the treatment of quantum mechanical systems with continuous spatial degrees of freedom at zero temperature. The PIGS method is a variant of the finite temperature path integral Monte Carlo method [6]. The key quantity in the finite temperature path integral Monte Carlo approach is the density matrix , , where denotes the quantum mechanical system Hamiltonian, the Boltzmann constant, and the temperature. Making the formal replacement , where is the time, the density matrix turns into the time evolution operator. Introducing the imaginary time , , and repeatedly acting with the “imaginary time evolution operator” , (assuming is, using a metric to be defined later, small) onto an initial state, the ground state wave function, or more precisely the lowest energy state that has finite overlap with the initial state, is projected out. This projection idea is the key concept behind the PIGS approach as well as many other imaginary time propagation schemes [7, 8]. Unlike grid based or basis set expansion based approaches, the PIGS approach is applicable to systems with varying degrees of freedom, i.e., few- and many-body systems. This versatility of the PIGS approach stems from the fact that the action of the imaginary time evolution operator on the initial (or propagated) state is evaluated stochastically, i.e., by means of Monte Carlo Metropolis sampling.

While the finite temperature path integral Monte Carlo algorithm, which—as has already been aluded to—has many features in common with the PIGS algorithm, has been reviewed quite extensively [6, 9], the PIGS algorithm has not, despite its generality, been reviewed in detail in the literature. The present paper thus serves three main purposes: (i) It develops, starting from equations that should be familiar to an advanced undergraduate student, the theoretical concepts behind the PIGS algorithm. (ii) It details how the relevant equations can be evaluated numerically, provides a good number of implementation details, and discusses various aspects regarding the algorithm performance. (iii) It presents applications of the PIGS algorithm to fermionic systems.

The PIGS approach allows one to solve the time-independent non-relativistic Schrödinger equation. Since the PIGS algorithm does not provide the full wave function in numerical or analytical form, the type of expectation values that one would like to determine needs to be specified a priori rather than a posteriori. In particular, an estimator has to be derived and implemented for each observable. The PIGS algorithm works, as already mentioned, through imaginary time propagation. It is imperative to clarify upfront that the imaginary time propagation is a numerical tool that facilitates projecting out unwanted excited state contributions. An extension to the real time evolution is, in general, not feasible, at least not for systems with a large number of degrees of freedom (for small systems, grid based real time propagation schemes do, of course, exist). A key ingredient of the PIGS algorithm is the stochastic evaluation of high-dimensional integrals, which arise from the discretization of the imaginary time and from the intrinsic degrees of freedom (particle coordinates) of the system under study. The stochastic Monte Carlo based approach to evaluating these integrals makes the PIGS method applicable to large systems containing as many as hundreds of particles. However, as in many Monte Carlo techniques, the treatment of identical fermions leads to the infamous Fermi sign problem. This tutorial applies the PIGS algorithm to small fermionic systems with zero-range interactions. It is shown that the sign problem can be “postponed” but not be avoided, i.e., application of the PIGS algorithm to larger fermionic systems with zero-range interactions will necessarily fail.

The PIGS algorithm has been applied to systems relevant to physics and chemistry. For example, our understanding of pristine and doped bosonic helium clusters of varying size has been informed by PIGS calculations [10]. Unlike alternative zero-temperature methods such as the variational Monte Carlo method [11, 7] and the diffusion quantum Monte Carlo method with mixed estimators [12, 13], the PIGS approach is known to yield unbiased results for structural properties like the radial density and pair distribution function. Here the term “unbiased” refers to the fact that the resulting structural properties are independent of the initial state, provided the initial state has finite overlap with the state of interest and provided the state of interest is the lowest energy state with a particular symmetry. Moreover, the condensate fraction [6, 14] and Renyi entropy [15, 16] are observables that can be calculated relatively straightforwardly within the PIGS algorithm (or at least more straightforwardly than within a number of other approaches). The PIGS algorithm has also been, among others, applied to bulk helium in varying spatial dimensions [14, 17, 18, 19, 20, 21, 22, 23, 24], liquid helium in nano-pores [25, 26], molecular para-hydrogen in nano-pores [27, 28], molecular para-hydrogen clusters [29, 30, 31], hardsphere bosons [32, 33], dipolar systems [34, 35], and cold atoms loaded into optical lattices [36].

The applications presented in this review deal with cold atom systems with infinitely large two-body -wave scattering length  [37, 38, 39], which are—like helium droplets—strongly interacting. However, the average interparticle spacing in cold atom systems tends to be significantly larger than that in helium droplets. This implies that the sample applications presented in this review deal with Hamiltonian that are characterized by vastly different length scales. To describe systems for which the average interparticle distance is large compared to the two-body interaction range, we employ two-body zero-range interactions. The use of two-body zero-range interactions removes the two-body range from the problem. If is send to infinity, as in the applications presented in this review, then two-component fermions are characterized by the same number of length scales as the corresponding non-interacting system [37, 38, 40, 41].

For bosons, in contrast, a three-body parameter, which can be defined in terms of the size of one of the extremely weakly-bound Efimov trimers, sets a length scale of the interacting system even if the range of the two-body interactions is zero. Since the use of two-body zero-range interactions in continuum Monte Carlo calculations is a fairly novel development [42, 43, 44, 45, 46, 47, 48], the associated implementation details are discussed in detail. Simulation results are presented for fermionic systems. Application of the algorithm to bosons requires only a few changes in the code; however, due to the existence of a three-body parameter, the number of time slices, e.g., is much larger than for fermionic systems with two-body zero-range interactions.

The remainder of this article is organized as follows. Section 2 introduces, starting from the non-relativistic Schrödinger equation, the key quantum mechanical equations behind the PIGS algorithm. Section 3 discusses a number of theoretical concepts that are needed to reformulate the basic quantum mechanical equations in a form amenable to computer simulations; many considerations in this section do not only apply to the PIGS algorithm but also to other Monte Carlo algorithms. Section 3.1 introduces some basic ideas. Sections 3.2 and 3.3 discuss two different approaches for approximating the short-time propagator, namely a Trotter formula based approach and an approach that utilizes the so-called pair product approximation; these two approaches are compared in Sec. 3.4. The use of two-body zero-range interactions within the pair product approximation is discussed in Sec. 3.5. Section 4 “translates” the formalism introduced in Secs. 2 and 3 into an algorithm. Section 4.1 introduces the basics of Monte Carlo sampling of high-dimensional integrals while Sec. 4.2 reviews formal aspects of the Monte Carlo Metropolis sampling. Section 4.3 discusses the generation of new configurations, i.e., the moves employed in the PIGS algorithm; as applicable, differences to the path integral Monte Carlo algorithm are pointed out. Sections 4.4 and 4.5 discuss the determination of expectation values and associated error bars, respectively. Last, Sec. 4.6 discusses how to treat permutations in the PIGS algorithm; this discussion is particularly relevant if the system contains two or more identical fermions.

Section 5 presents a number of applications to harmonically trapped equal-mass two-component Fermi gases. The simulation results are discussed from two different angles. On the one hand, “technical aspects” such as convergence with respect to the propagation time and the time step are discussed. On the other hand, the physical relevance of the simulation results presented is highlighted. Spin-balanced systems with up to particles and a non-interacting Fermi gas with a single impurity with up to particles are considered. In both cases, interspecies two-body zero-range interactions with infinitely large -wave scattering length are employed. The construction of different types of trial functions is discussed and the dependence of the simulation results on is elucidated. PIGS results for the energy, pair distribution function, and contact are presented and compared to literature results where available. Last, Sec. 6 provides a summary and an outlook.

2 Quantum mechanical foundation

We consider non-relativistic particles described by the time-independent Hamiltonian at zero temperature. The Hamiltonian may contain single-particle potentials, two-body potentials, and higher-body potentials. We work in position space, where the potentials are local, i.e., we consider potentials that only depend on the position vectors and not on the momentum vectors as would be the case if, e.g., spin-orbit coupling terms were present [49, 50, 51]. The position vector for the -th particle with mass is denoted by and we collectively denote the position vectors of all the particles by , . The stationary eigen states and corresponding eigen energies are denoted by and , where . The form a complete set and we are, throughout this article, interested in systems that support at least one -body bound state. The treatment of scattering states by means of quantum Monte Carlo approaches is, in general, a challenging task [52, 53, 54, 55] that is beyond the scope of this paper. The PIGS algorithm allows one to calculate a subset of the bound state energies as well as expectation values such as the pair distribution functions associated with the corresponding eigen states.

The PIGS algorithm is rooted in imaginary time propagation, a concept that is used widely to find the ground state or selected excited states of linear and non-linear Schrödinger equations [8, 56]. The concept of imaginary time propagation is also used to solve non-quantum mechanical wave equations. In what follows, we restrict ourselves, for concreteness, to the linear Schrödinger equation. To illustrate the key idea behind imaginary time propagation algorithms, we assume that the ground state is non-degenerate, i.e., that for . We consider an initial trial function , which does not have to be normalized, that has finite overlap with the ground state wave function . To analyze what happens when this trial function is propagated in imaginary time, we decompose the trial function into the eigen states of the Hamiltonian ,


where is non-zero by assumption. Using Eq. (1), ,


can be written as


Since is, by assumption, greater than , the excited states contained in decay out during the imaginary time propagation. In the limit, approaches, except for an overall factor, the eigen state . Correspondingly, the energy ,


approaches the exact ground state energy exponentially in the limit. For finite , provides an upper bound to the exact eigen energy. This suggests that one can obtain a reliable estimate of by extrapolating the for various finite to the limit. Expectation values of an arbitrary operator can be written analogously,


where denotes the -dependent expectation value. The convergence of toward the exact expectation value with respect to may not be simply exponential and needs to be analyzed carefully for each operator (see Sec. 5 for examples).

Equations (2), (4), and (5) constitute the starting point of the PIGS algorithm (see Sec. 4). Based on these equations, two ingredients or components of the PIGS algorithm can already be identified. (i) An initial trial function needs to be supplied by the “simulator”. From Eq. (3) it is clear that the efficiency of the PIGS algorithm depends on the overlap between and : If all with vanish, then the imaginary time propagation is not needed at all. If the for the states that lie energetically close to vanish, then the decay of the excited states is fast, i.e., small should suffice. The construction of is, of course, strongly dependent on the Hamiltonian under study. Examples are discussed in Sec. 5. (ii) Given an initial trial function , the action of onto needs to be evaluated. The PIGS algorithm as well as many other algorithms accomplish this by dividing into multiple smaller imaginary time steps . While non-Monte Carlo based approaches are, typically, restricted to relatively small system sizes, the PIGS algorithm as well as other Monte Carlo algorithms are designed to treat systems for which can be a high-dimensional vector.

The discussion thus far focused on determining the absolute ground state of the system. The outlined formalism can be readily adopted to the determination of the energetically lowest-lying state with a given symmetry. For concreteness, let us assume that the total angular momentum and the total parity are good quantum numbers and that the absolute ground state has vanishing angular momentum () and positive parity (). If is chosen to have a symmetry other than , say symmetry, then the imaginary time propagation projects out the eigen state with symmetry that has the lowest energy. Said differently, the imaginary time propagation preserves the symmetry of and acts in a subspace of the full Hilbert space.

It is instructive to compare the PIGS formalism discussed above with another imaginary time propagation based Monte Carlo technique, namely the diffusion Monte Carlo technique (for the purpose of the discussion that follows, the Green’s function Monte Carlo technique behaves identically) [7]. The diffusion Monte Carlo approach utilizes, in addition to a trial function, a reference energy that is adjusted continually during the simulation. If Eq. (3) is multiplied by , then the right hand side is, except for an overall -independent factor, independent of for sufficiently large and . This is the key idea behind the diffusion Monte Carlo approach. The accumulation of expectation values is started after the excited state contributions have decayed out and after has been adjusted such that . While the diffusion Monte Carlo and PIGS approaches share, as just discussed, similarities, the treatment of particle permutations differs notably if the system contains two or more identical fermions. The diffusion Monte Carlo algorithm does not explicitly apply sequences of two-particle permutation operators; the identical particle characteristics (bosons and/or fermions) of the -particle wave function are instead encoded via the trial function, combined with the fixed- or released-node approach in the case of identical fermions [57, 58]. The PIGS algorithm, in contrast, explicitly anti-symmetrizes the paths if the system contains identical fermions. If the system contains identical bosons, explicit symmetrization operations are not needed provided the ground state of the system where the bosons are replaced by “Boltzmann particles” is the same as that of the system containing bosons.

3 PIGS algorithm: General considerations

3.1 Basic concepts

This section rewrites Eq. (4) in a form amenable to evaluation by Monte Carlo techniques. The actual Monte Carlo implementation is discussed in Sec. 4. Using and , Eq. (4) reads


The denominator is commonly denoted by ,


To obtain a prescription for evaluating the operators in the integrands, we project and onto the position basis. Formally, this amounts to inserting the identity


where denotes the unit operator, multiple times into Eq. (6),




We refer to ,


as the imaginary time evolution operator projected onto the position basis or, in short, as the imaginary time evolution operator or propagator. Using Eq. (11), we obtain


The normalization factor , Eq. (7), plays a key role in the simulations. For example, if is known, instead of evaluating Eq. (12), one can calculate the energy expectation value directly,


Equations (12) and (13) generate two distinct energy estimators (see Sec. 4.4 for details).

In the zero propagation time limit, i.e., for , becomes the identity operator. This implies that the propagator is simply a -function in position space,


To propagate to finite imaginary time, one can solve the Bloch equation [6]


which is obtained by taking the partial derivative of the propagator with respect to . Equation (15) can be interpreted as a diffusion equation in the imaginary time . For the remainder of this section, we write the Hamiltonian as a sum of the kinetic energy operator and the potential energy operator . Moreover, we assume that all particles have the same mass ; this assumption, which can be readily relaxed, simplifies the notation. If the kinetic energy operator is zero, the propagator can be readily solved for. Similarly, if the potential energy operator is zero, the propagator can also be solved for. In this case, the solution corresponds to free particles diffusing in space (the subscript “0” is used to indicate that the corresponding Hamiltonian contains only kinetic energy terms), i.e., is a product of single-particle Gaussians,


where is equal to . Equation (16) shows that the off-diagonal terms (terms for which ) of , expressed in the position basis, are non-zero. This shows explicitly that the kinetic energy operator is non-local in position space. If and are both non-zero, then the propagator at finite is known only for a few selected problems such as non-interacting particles in a harmonic trap [59] and two particles with zero-range interactions [60, 61, 62, 63, 45]. In general, the -particle imaginary time evolution operator or propagator is unknown. If it was known, the problem would be “trivial”.

The PIGS algorithm is based on the idea of writing the imaginary time evolution operator for large as a product over imaginary time evolution operators for small imaginary time steps,


Using Eq. (17) in Eq. (11) and inserting the unit operator [Eq. (8)] times, we obtain




The problem of evaluating the propagator at the desired imaginary time has been converted to evaluating propagators at and integrating over (potentially high-dimensional) auxiliary coordinates . The key points are that one can typically find a fairly accurate but approximate short-time propagator for finite (see Secs. 3.2-3.5) and that the associated “auxiliary” integrations can be performed efficiently by Monte Carlo techniques (see Sec. 4).

To simplify the notation, the product in the denominator of Eq. (12) is rewritten as . Each set of coordinates inserted in Eq. (18) is referred to as a “time slice”. There are three “special” time slices: the initial time slice , the middle time slice , and the final time slice . The initial and middle time slices are connected by the propagator and the middle and final time slices are connected by the propagator . Since both propagators are rewritten by inserting auxiliary time slices, the “expanded” partition function contains a total of time slices. The position vector of the -th particle in the set of coordinates is referred to as a “bead”. Thus, a single particle is represented by beads. The propagator that “connects” two consecutive time slices is referred to as a “link”. The inverse temperature corresponding to this link is , where . The propagator that “connects” two consecutive beads is referred to as a “single-particle link”. In addition, the set of all time slices is referred to as a configuration. The definitions are summarized in Table 1.

bead a single coordinate of the -th particle
at the -th imaginary time index
time slice a set of beads at the -th
imaginary time index
configuration the set of all time slices
link the propagator connecting two
consecutive time slices
single-particle link the propagator connecting two
consecutive beads
Table 1: PIGS terminology used in this article. Columns 1-3 show the term, symbol, and explanation, respectively.

Figure 1

Figure 1: Schematic world-line representation for a single particle with mass in a one-dimensional harmonic trap with angular trapping frequency . Panels (a), (b), (c), and (d) show paths for , 5, 9, and 17 time slices, respectively.

shows the world-line representation of a single particle in a one-dimensional harmonic trap. World lines move in position space (-axis) and imaginary time (-axis). Figures 1(a), 1(b), 1(c), and 1(d) show paths for 3, 5, 9, and 17 beads, respectively. As increases, the path is resolved in more detail (each link is more accurate) and observables calculated based on the sampled paths become more accurate.

Figure 2 depicts a single particle (the Hamiltonian only contains the kinetic energy term) in two-dimensional space [64]. Two consecutive beads (circles in Fig. 2) are connected by a single-particle link (wiggly line in Fig. 2). The kinetic energy is “carried” by the propagators represented by the links. The expression for the propagator in free space reads [Eq. (16) for a single particle with position vector ]


The action  [6],


of the single-particle link that connects the beads labeled and reads


It can be seen that the action has the same form as that of a “spring potential” for two classical particles with position vectors and connected via Hooke’s law. The propagator can thus be interpreted as being proportional to the Boltzmann factor of a classical system of springs. Note that and in Eqs. (20) and (22) correspond to position vectors of consecutive beads for one single particle while and in the classical isomorphism correspond to position vectors of two different particles.

Figure 2: Isomorphism between path integrals for a single free particle in two-dimensional space and classical particles connected by springs. In the path integral interpretation, the circles and wiggly lines depict the beads and (single-particle) links of a single particle, respectively. In the classical mechanics formulation, the circles and wiggly lines depict particles and springs, respectively. The position vector of the particle is denoted by , where indicates the imaginary time index (, and ).

In addition to the propagators and , the partition function contains the trial functions (or “weights”) and . For the single two-dimensional particle in free space, and can be interpreted as potentials that are felt by the first and last particle of the chain of classical particles. Thus, we can interpret the PIGS simulation of a single particle as a simulation of a chain of classical particles connected by springs (or a polymer with nearest neighbor interactions) and two external forces that act on the particles at the ends of the chain. The stiffness of the springs is determined by , i.e., by the inverse of the imaginary time step associated with the links.

The following two sections introduce two different approaches for approximating the short-time propagator, namely the Trotter formula and the pair product approximation.

3.2 Trotter formula

One way to approximate the short-time propagator is to use the Trotter formula [65]. For sufficiently small time steps , the kinetic energy contribution and the potential energy contribution to the propagator can be split,


where the notation indicates that the leading-order error scales, in general, as to the power of 2. More specifically, by Taylor expanding the exponentials, one can prove that the leading-order error is , where is the commutator between and , . In the limit (this corresponds to the insertion of infinitely many time slices, i.e., the limit), the Trotter formula becomes exact. Since cannot be infinitely large in practice, one performs calculations for different and extrapolates the observables of interest to the infinite limit.

Importantly, the Trotter formula can be extended to higher orders. We can readily adopt a scheme by further splitting the kinetic energy term or the potential energy term,




In practice, Eq. (25), which is accurate up to second order [the error is ], is easier to use than Eq. (24). In position space, Eq. (25) reads


where [see Eq. (16)] is the propagator that accounts for the kinetic energy term.

One can reach successively higher accuracy by the repeated use of the Baker-Campell-Hausdorff formula (see, e.g., Ref. [66])




Using Eqs. (27) and (28) twice, we obtain [67]




Applying Eqs. (29) and (30) twice to


we can check that the fourth-order factorization [67]


where is given by , holds. The term , in position space, can be simplified to , where the gradient in three spatial dimensions is given by


with , , and denoting unit vectors that point in the , , and directions, respectively. The term corresponds to the square of the force on the -th particle. Care needs to be taken in evaluating the derivatives, since usually contains a double sum over two-body potentials or even a triple sum over three-body potentials. In most cases, the evaluation of the force terms cannot be simplified analytically, implying that the evaluation of double commutators involves double or triple sums over the total number of particles. This makes the numerical evaluation comparatively expensive. Note that the exponentials in Eq. (32) that contain the potential energy can have different numerical factors. In addition, Eq. (32) is not the only form of the fourth-order factorization [67, 68].

Using the Trotter formula, the isomorphism between a single quantum particle in free space and the classical spring system can be extended to multiple quantum particles with interactions. Figure 3 depicts two interacting particles in two-dimensional space (it is assumed that the particles do not feel a single-particle potential).

Figure 3: Isomorphism between path integrals for two interacting particles in two-dimensional space and classical particles connected by springs. In the path integral interpretation, the circles and wiggly lines depict the beads and links of a single particle, respectively, and the dotted lines depict the two-body interaction between beads with the same imaginary time index. In the classical mechanics interpretation, the circles and wiggly lines depict particles and springs, respectively, and the dotted lines depict the two-body interaction between selected classical particles. The position vector of the -th particle at the -th imaginary time index is denoted by .

In this case, the propagator for the link that connects the beads labeled by , , , and reads [see Eq. (25)]


where denotes the two-body interaction potential between particles 1 and 2 and the single-particle propagator of the -th particle [see Eq. (20)]. As in Fig. 2, two consecutive beads for the same particle (e.g., the circles labeled by and in Fig. 3) are connected by single-particle links (wiggly lines in Fig. 3) that represent the propagators . Since the two-body interaction potential (dotted lines in Fig. 3) is diagonal in position space [see the exponentials on the right hand side of Eq. (34)], it connects beads of different particles with the same index, i.e., it connects and , and , and and (or, in general, it connects and ). Each particle in the PIGS simulation corresponds to classical particles connected by springs. Classical particles associated with different chains interact only if they have the same bead index.

3.3 Pair product approximation

To introduce the pair product approximation, we assume for simplicity that the potential energy operator can be written as a sum of two-body terms,


i.e., we assume for now that single-particle and three- and higher-body forces are absent. Under these assumptions, the short-time propagator can be evaluated using the pair product approximation [6]. It is convenient to define the two-body kinetic energy operator for the -th and -th particle in position space,


The relative non-interacting and interacting two-body Hamiltonian are and , respectively. The pair product approximation considers two-body correlations explicitly, but not higher-body correlations, and writes the many-body propagator as a product over single-particle propagators and two-body propagators,


where ,


is the reduced pair propagator. The denominator of the reduced pair propagator coincides with the known relative non-interacting two-body propagator. Thus the only “non-trivial input” is the relative propagator of the interacting two-body system. One can readily see that the pair product approximation is exact for two particles for any propagation time because the center-of-mass and relative degrees of freedom separate in this case. In some cases such as for the two-body zero-range interaction potential, the exact reduced pair propagator is known analytically [60, 62, 61, 63]. In other cases such as for the two-body hardcore potential, an approximate reduced pair propagator is known analytically in closed form [69, 70]. If the reduced propagator is not known analytically, one can perform a partial wave decomposition and obtain a numerical representation of the reduced two-body propagator [6].

In dilute gases or weakly-bound droplets, the interparticle spacing is typically so large that two-body collisions dominate over three- and higher-body collisions. The leading-order error of the pair product approximation is determined by the importance of three-body correlations. For two-component equal-mass Fermi gases with two-body zero-range interactions, three-body correlations are suppressed by the Pauli exclusion principle. For this system, we found that the pair product approximation provides an extremely good description of the propagator. Specifically, we obtain accurate simulation results for a small number of beads (see Sec. 5 for details). For bosons, in contrast, three-body correlations can be significant. As a consequence, the pair product approximation is not as efficient as for two-component fermions and simulations typically employ a large number of beads (“large” in this context means about two orders of magnitude more number of beads as in the simulations for fermions [45]).

To illustrate the pair product approximation, we cannot use the classical isomorphism because the kinetic and potential energy contributions are mixed. One needs to evaluate the single-particle propagator, which can be represented by springs as in Figs. 2 and 3. However, one also needs to evaluate the reduced two-body propagator, which connects two consecutive beads of one particle’s path with the same consecutive beads of another particle’s path. These “four-bead connections” do not have a simple classical analog.

3.4 Comparison of the two approximations

This section discusses the advantages and disadvantages of approximating the short-time propagator with the help of the Trotter formula and the pair product approximation.

In the Trotter formula based scheme, the kinetic and the potential energy terms are treated separately. Inserting the identity , Eq. (8), multiple times into Eqs. (25) or (32), it can be seen that the potential energy is diagonal in position space. This means that one can directly evaluate the potential energy term at each time slice. The kinetic energy term contains off-diagonal terms and needs to be evaluated at each link instead of at each time slice. Nevertheless, since the kinetic energy term corresponds to a simple Gaussian, the sampling of the kinetic energy piece of the propagator can, in general, be performed efficiently (see Sec. 4.3.2 for details).

Even though the Trotter formula can formally be generalized to expressions that are accurate to order  [71, 72], many of these expressions are of limited use in practice because they contain either commutators that involve rather complicated expressions or terms that correspond to negative imaginary time, which are not normalizable. There exists a multi-product expansion for the propagator [72]; however, applications thereof are still rare [73]. Thus present-day algorithms mostly employ Trotter formula based decompositions that are accurate to order .

In the pair product approximation, the two-body reduced propagator contains kinetic energy and potential energy contributions. This means that the reduced two-body propagator has to be evaluated at each link. Because the reduced two-body propagator is, in general, not a simple Gaussian, the sampling is typically less efficient than in the case where the Trotter formula is used. Furthermore, as discussed in Sec. 4.6, the evaluation of the permutations is computationally more involved.

Our discussion of the pair product approximation assumed that the potential energy can be written as a sum over two-body terms. If the Hamiltonian contains one-, three-, or higher-body potential energy terms, one can include them by combining the Trotter formula and the pair product approximation. To this end, one first splits the propagator into two terms using the Trotter formula. The first term contains the kinetic energy operator and the two-body interactions while the second term contains all other potential terms. One then applies the pair product approximation to the first term. In this approach, it is most convenient to use the second-order Trotter formula for two reasons. First, if a higher-order Trotter formula was used, one would need to evaluate the commutator between the one-, three- and higher-body potential terms and the two-body potential terms. For the two-body zero-range interactions considered in Sec. 5 this is a rather challenging task. Second, both the second-order Trotter formula and the pair product approximation yield errors for the energy that scale quadratically with the time step . While the error in the pair product approximation tends to be smaller than that associated with the pair product approximation, ultimately it is the scaling with that determines the accuracy and use of the fourth-order Trotter formula typically leads only to a small overall improvement.

From our perspective, the pair product approximation has one key advantage: It can deal with a class of two-body potentials that the Trotter formula based scheme cannot deal with (at least no such treatment is known to us). For example, the two-body hardcore and zero-range potentials contain infinities and can thus not be treated by the Trotter formula based scheme. However, the infinities can, as discussed in the next section exemplarily for the two-body zero-range potential, be dealt with analytically in the pair product approximation.

3.5 Propagator for two-body zero-range interactions

As alluded to in the previous section, two-body zero-range interactions can be incorporated into continuum Monte Carlo simulations through the pair product approximation [45, 42, 44, 46], which employs the relative propagator [see Eqs. (37) and (38)]. In what follows, we limit our discussion to three spatial dimensions. To determine , one considers the relative Hamiltonian for two particles interacting through the regularized zero-range Fermi-Huang pseudopotential  [74] in free space,




Here, denotes the two-body reduced mass, the interparticle distance vector, and the -wave scattering length. The regularization operator in Eq. (40) ensures that the Hamiltonian is well-behaved. Without this operator, the two-body coupling constant would have to be renormalized. With the regulator, however, the coupling strength is uniquely determined and given by .

The reduced (or normalized) relative propagator corresponding to the Hamiltonian given in Eq. (39) reads [61, 63]


where and . For , the length scale drops out of the expression for the propagator and Eq. (41) simplifies to


Importantly, the reduced relative propagator diverges when or approach zero. These divergencies have implications for the Monte Carlo sampling of the paths. As discussed in detail in Sec. 4.3.3, moves have to be designed carefully such that detailed balance is fulfilled. For example, while diverges for and , the probability to find two particles at vanishing interparticle distance does not diverge. The treatment of systems with two-body hardcore interactions is similar in spirit to that detailed here for two-body zero-range interactions.

Adding the spherically symmetric harmonic confining potential for the relative degrees of freedom to the Hamiltonian given in Eq. (39) and assuming that is infinitely large, the reduced relative propagator reads


where . In the limit that the angular trapping frequency goes to zero, Eq. (43) reduces to Eq. (42). Expression (43) is used in Sec. 5 to treat harmonically trapped two-component Fermi gases with two-body zero-range interactions at unitarity using the pair product approximation.

4 Monte Carlo techniques and the PIGS algorithm

Throughout this section we assume that the trial function is given and that its value can be determined for any set of coordinates . The functional form of depends sensitively on the system under study. The choice of and the dependence of the PIGS results on will be discussed in Sec. 5 for harmonically trapped two-component Fermi gases.

4.1 General sampling scheme: Importance sampling

Equation (19) writes the long-time propagator as a high-dimensional integral over a product of short-time propagators. This implies that the evaluation of the normalization factor is a high-dimensional integral. This section discusses the Monte Carlo sampling of this high-dimensional integral over (there are independent coordinates if we are considering three spatial dimensions). To proceed, we write explicitly in terms of the short-time propagator,




To simplify the notation, we denote the configuration by and the probability distribution by . The notation of these and other quantities is summarized in Tables 1 and 2.

probability distribution Eq. (45)
probability density function Eq. (48)
weight function (observable specific) Eqs. (46)-(47); Sec. 4.4
transition probability Eq. (49)
proposal distribution (selected by simulator) around Eqs. (49)-(50)
acceptance distribution Eq. (50)
trial function (selected by simulator) Eq. (1)
Table 2: Definitions of Monte Carlo sampling terminology used in this article. Columns 1-3 show the symbol, name, and related equation number, respectively. The configuration is defined as .

The expectation value of an arbitrary observable can be written as


where the integration goes over coordinates and the weight function needs to be determined, as will be discussed in Sec. 4.4, for each observable. To see the structure of more clearly, we rewrite Eq. (46) as


where the probability density function is defined as


In contrast to the probability distribution , the probability density function is normalized; and represent the weight contributed to the observable by the configuration and the normalized probability to be in the configuration , respectively. Equation (47) provides the basis of importance sampling: Configurations are not blindly distributed uniformly in space but instead are distributed according to . The advantage of importance sampling is that most computer time is used to sample configurations that are physically relevant and little time to sample configurations that do not contribute significantly to .

The general idea of the PIGS algorithm is to generate configurations according to and to use the generated configurations to accumulate the weight functions for a set of observables. Thus, it is crucial to have correct and efficient sampling schemes that explore the full configuration space with a relatively high acceptance ratio and without getting stuck around a local minimum. Section 4.2 reviews the basics of selected Monte Carlo methods, which are then used in the subsequent sections.

4.2 Some background on Monte Carlo methods

This section discusses how to update or generate configurations using the Metropolis algorithm. A Markov process is uniquely defined by the transition probability to go from configuration to configuration . The Metropolis algorithm satisfies the detailed balance condition [7]


which states that the flow of probability from to is equal to that from to . This means that there is no net flow of probability. The Metropolis algorithm needs to ensure ergodicity of the Markov process. If the process is ergodic, the Markov chain (i) returns to any previously generated configuration after a sufficiently long simulation time and (ii) is not periodic (a Markov chain of , e.g., is periodic). The ergodicity ensures that the probability distribution gets sampled fully. For example, as discussed in Ref. [75], if we use the traditional scheme of treating the permutations [6], for a two-component Fermi gas with zero-range interactions, the Markov process ends up with a configuration in which all particles sit on top of each other and the configuration never returns to the original configuration . This means that ergodicity is violated and that the Markov process does not generate samples according to . This renders the sampled configurations meaningless. We note, however, that while the detailed balance condition together with the ergodicity guarantees that the equilibrium distribution coincides with the desired probability distribution , there exist other Monte Carlo methods that do not satisfy the detailed balance condition but yield an equilibrium distribution that coincides with the desired probability distribution  [76].

The Metropolis algorithm consists of two steps [7]: (i) the generation of a proposed configuration (or move) and (ii) the acceptance or rejection of the proposed configuration (or move). The combination of (i) and (ii) leads to a new configuration. Starting from the configuration , we propose a new configuration according to a proposal distribution and accept (the new configuration would be ) or reject the new configuration (the new configuration would be ) according to the acceptance distribution . This implies that is given by . The Metropolis algorithm chooses such that [7]


We verify that the detailed balance condition [Eq. (49)] is satisfied in the following. If is smaller than , we obtain from Eq. (50) that




Plugging the right-hand sides of Eqs. (51) and (52) into Eq. (49), we confirm that Eq. (49) holds. If is larger or equal to , it can be checked similarly that Eq. (49) holds. Thus, we have shown that detailed balance is fulfilled.

A key task is to design proposal distributions that ensure the complete and, ideally, efficient exploration of the entire configuration space. In most cases, efficient simulation schemes are achieved if more than one proposal distribution (and hence type of move) is utilized. As discussed more in the next section, the proposal distribution might be designed based on the knowledge of the non-interacting system (see, e.g., Sec. 4.3.2) or based on knowledge of certain limiting behaviors of the interacting system (see, e.g., Sec. 4.3.3).

In practice, the acceptance ratio , i.e., one minus the fraction of rejected moves, should be monitored (note, is a number and not an - and -dependent function). The acceptance ratio for Metropolis sampling is different from the acceptance ratio encountered in the rejection sampling. In the rejection sampling, a rejected configuration does not lead to a new configuration. In the Metropolis sampling, in contrast, a rejected configuration does lead to a new configuration. When a configuration is rejected, the old configuration becomes the new configuration. For most of the updates (i.e., the generation of proposed new configurations), the acceptance ratio should not be too large and not be too small. A high acceptance ratio typically implies that the deviation between the old and new configurations is, on average, small. This means that the configuration space is explored comparatively slowly. A low acceptance ratio, in contrast, means that the Markov chain contains many identical configurations; again, typically this means that the configuration space is explored comparatively slowly. Both cases can result in large correlations of the sample and should be avoided. As a rule of thumb, the acceptance ratio should lie roughly between 30% and 50% [77].

4.3 Moves

The previous section outlined the basics of the Metropolis algorithm. This section discusses the PIGS moves that are used to update the configurations. For all moves, the proposed new configuration is chosen based on the proposal distribution and accepted/rejected based on the acceptance distribution . Once is specified, follows from Eq. (50). This section discusses three different moves. The “naive move” and the “wiggle move” are “all purpose” moves, which have proven to be useful for nearly all systems. In some cases, the use of these two types of moves alone does not lead to an efficient (or even correct) exploration of the configuration space. For systems with two-body zero-range interactions, e.g., the “pair distance move” is needed. In general, the simulator decides on the frequency with which the individual moves are used. The optimal ratio can be found empirically from the performance of the simulation itself or through the implementation of some sort of learning algorithm.

The list of moves discussed below does not include a “permutation move”, i.e., the stochastic sampling of the permutations is not discussed. The reason for this is twofold. If the system contains identical bosons, the ground state wave function is typically identical to that of Boltzmann particles, eliminating the need for an explicit symmetrization. If the system contains identical fermions, we employ the on-the-fly anti-symmetrization scheme discussed in Sec. 4.6, which is particularly useful if two-body zero-range interactions are employed.

The moves employed in the PIGS algorithm have much in common with the moves employed in the finite-temperature path integral Monte Carlo algorithm. As already discussed, one difference is that the PIGS algorithm contains the trial function while the finite-temperature path integral Monte Carlo algorithm does not. Quite generally, whether a move depends on the trial function or not depends on whether or not the beginning time slice or the ending time slice are being updated. In the implementations of the moves discussed below, the wiggle move does not depend on , and the naive move and the pair move may depend on (it depends on whether or not the randomly selected time slice to be updated is the -th or -th time slice). It should be noted, though, that the moves can, in principle, be implemented in a variety of ways, i.e., slightly different proposal distributions might be used in different implementations and be referred to by the same move name.

4.3.1 Naive move

The simplest move (the naive move) consists of shifting the position vector of a single bead by , where is drawn uniformly from the interval . The basic idea behind this move is that the propagator is a smooth function of and that a small change in does not introduce a huge change in the probability distribution . The size of the interval (if we are simulating a three-dimensional system, then the interval corresponds to a cube) can be adjusted such that the acceptance ratio of the proposed new position vector is around 50%. The proposal distribution is equal to a constant if the new bead lies in the interval ; otherwise, is equal to . Following Eq. (50), the move is accepted according to


Importantly, one cannot choose an “unbalanced” interval like , where , since the detailed balance condition, Eq. (49), is not satisfied in this case. The reason is that it is possible to go from to in one move but that it is impossible to go from to in one move.

The algorithm for the naive move can be summarized as follows. (i) Let the current configuration be . Randomly select a particle index and a time slice index , where can take any value from to . Set and calculate the old probability distribution . (ii) Generate a new position , where is drawn uniformly from the interval . Define and calculate the new probability distribution . (iii) Calculate the ratio . If this ratio is larger than a random number drawn uniformly from 0 to 1, accept the move and set ; otherwise, reject the move and set (i.e., do not change ). Figure 4 illustrates the naive move for a single particle in a one-dimensional harmonic trap. It can be seen that the proposed move involves only one bead.

Although the naive move attemps to change only one bead at a time, whether the proposed move gets accepted or rejected depends, in principle, on all the beads, i.e., the coordinates of all particles at all time slices, since the acceptance/rejection depends on the ratio . However, because the probability distribution is a product over propagators and the trial function evaluated at the two ends [Eq. (45)], the only terms that contribute to the ratio are the propagators and and, if equals or , the trial function or . If one uses the second-order Trotter formula [Eq. (26)], which treats the potential and kinetic energy terms separately, the terms that contribute to and are the potential energy term and the kinetic energy terms and .

The caveat of the naive move is that the correlation length is typically large. In the best case scenario (i.e., in the case where all beads of all particles are considered exactly once and all proposed moves are accepted), moves are needed to generate a configuration in which every bead differs from the starting configuration. Thus, we calculate observables for every -th configuration, where is a constant greater than 1 that is adjusted to ensure that the observables are calculated from configurations with small correlations. In practice, we find that lies between 2 and 20 for the applications considered in this review.

Figure 4: Illustration of the naive move for a single particle in a one-dimensional harmonic trap for beads. The black circles depict the old bead positions. The red square shows the proposed bead position for the 16-th time slice index. It can be seen that only two links (namely the link involving the 15-th and 16-th beads and that involving the 16-th and 17-th beads) are changed.

4.3.2 Wiggle move

The wiggle move uses the non-interacting propagator to design the proposal distribution . Since the non-interacting propagator is a product of simple Gaussians for particles in free space [Eq. (16)] and for particles confined in a harmonic trap [59, 45], one can generate configurations efficiently with 100 % acceptance ratio using the Box-Muller transformation [78] or with finite acceptance ratio using the Marsaglia polar method (the acceptance ratio is around 80 %) [79] or the Ziggurat algorithm (the acceptance ratio is around 98 %) [80, 81]. If the difference between the propagator of the system to be simulated and the non-interacting propagator is small, the acceptance ratio for a move generated based on the propagator of the non-interacting system is high. Despite the large acceptance ratio, the correlation between consecutive configurations is, in general, small. In the non-interacting limit, the acceptance ratio is exactly 1.

Depending on the number of beads changed simultaneously, the wiggle move is a single-bead move or a multi-bead move. The single-bead version of the wiggle move randomly selects a particle index and a time slice index ( and ). Since the beads to be moved exclude the time slices and , the wiggle move does not involve the trial function . We denote the new proposed position vector by (how to choose