Path Integral Monte Carlo

The Path Integral Monte Carlo (PIMC) has become very popular recently. This method is based on Feynamn's formulation of quantum mechanics in terms of path integrals [53].

At finite temperature, equilibrium properties of quantum systems can be calculated using the thermal density matrix, $\hat \rho$.

$$\hat \rho = e^{-\beta \hat H},$$ (47)

where $\beta = 1/k_b T$ is the inverse temperature. The expectation value of an observable $\hat O$ at equilibrium:

$$<\hat O> = Z^{-1} Tr \left( \hat \rho \hat O \right) = Z^{-1} Tr \left( \hat O e^{-\beta \hat H} \right),$$ (48)

where $Z = \rm Tr( e^{-\beta \hat H})$ is the partition function.

The PIMC method is usually implemented in a position-space representation $\vert R> = \vert\vec r_1,\vec r_2,\vec r_3,\vec r_4, ... >$, where the density matrix is given by:

$$\rho (\bR,\bRs;\beta ) = <\bR\vert e^{-\beta \hat H}\vert\bRs>,$$ (49)

where the elements of $\rho (\bR,\bRs;\beta )$ are positive and can be interpreted as probabilities. The $\rho (\bR,\bRs;\beta )$ contains the exponential of the Hamiltonian $\hat H = \hat K + \hat V$, where $V$ is the potential energy and $K$ kinetic energy. These operators are non-commutative, which makes them difficult to evaluate density_matrix. For small $\beta$ (high temperature) one can approximate

$$e^{-\beta (\hat K + \hat V)} \simeq e^{-\beta \hat K}e^{-\beta \hat V} +O(\beta^2).$$ (50)

If $\beta$ is large (low temperature), then one divide $\beta$ into many small segments, $M$, i.e. $\tau = \beta/M$, which can be treated as a case of small $\beta$. The partition function, $Z$,

$$Z = \int d\bR_0 <\bR_0\vert-\beta \hat H \vert\bR_0 >$$ (51)

can written in the following form:

$$Z = \int d\bR_0 d\bR_1 d\bR_2 ... d\bR_{M-1} <\bR_0\vert-\t... ...t H \vert\bR_2 >...<\bR_{M-1}\vert-\tau \hat H \vert\bR_{0} >,$$ (52)

where the $M-1$ units operators $\int \bR_i <\bR_i\vert\bR_i>$ are inserted, and the first $\vert\bR_0 >$ and the last $\vert\bR_M >$ elements are identified as required by the trace operations. In the short-time approximation the expression for $\rho (\bR, \bRs; \tau )$ is given by:

$$\rho (\bR,\bRs;\tau ) = <\bRs \vert e^{- \tau (\hat K + \hat... ... -\bRs\vert^2}{4\tau } \right) \exp\left(-\tau V(\bR) \right),$$ (53)

which is the same form as Green_full. Substituting the result into expas one obtains

$$Z = 1/(2\pi\tau)^{3MN/2}\int d\bR_0 ...d\bR_{M-1} \exp \left... ...frac{\bR_{m+1} -\bR_{m}}{\tau} \right)^2 + V(\bR_{m}) \right).$$ (54)

The integral PIMC looks like a partition function of a classical ring polymer, where the nearest neighbors (beads) are connected with springs, while the $V(\bR)$ is an external potential. Moving a quantum particle is equivalent to evolving this polymer. The system of quantum particles can be represented by an ensemble of these ring polymers, as shown in Fig. 5.1.

Figure 5.1: The trace of the closed paths of 54 helium atoms in the bcc lattice at a temperature of 1.5 K with M = 25 beads. Straight lines connect successive beads of the polymers. The PIMC (upicode 9.0) is used for calculations.

Figure 5.2: Cartoon of paths for 3 particles (A, B, C) with M=4, and $\tau = \beta /4$. The kinetic terms are represented by zig-zag lines and potential by dotted lines, which connect only beads of the polymers at the same time, from Bernu and Ceperley [20]

The only difference between the classical and quantum systems is that the potential interaction between these polymers occurs only at the same time slices (beads), while in a real polymer all monomers interact with each other.(See Fig. 5.2).

The expectation values are calculated by sampling the paths with the MC method. The simplest move is the translation of a whole chain. This move changes the potential energy of the system. If the relative position of the monomers is altered, then one changes the kinetic energy. Several algorithms are used to sample the paths efficiently [22].

The Bose and Fermi symmetry of the quantum particles increases significantly the complexity of the PIMC method. In the case of bosons, the density matrix is written as:

$$\rho (\bR,\bRs;\beta ) = \frac{1}{N!}\sum_P\rho (P\bR,\bRs;\beta ),$$ (55)

where $P$ is the permutation operator, for example a pair permutation $P \psi(1,2) =\psi(2,1)$. The action of the permutation operator is equivalent to opening the ring polymers involved in the permutation and making a single polymer out of them.

In PIMC (for bosons) no a priori is needed, one can use only the Hamiltonian and let the code run for a long time. It is possible that the paths themselves make the transition to a new, unexpected state. There is no problems with scalability in the method: increase of the number of polymers does not lead to instability, the computational effort increases, the correlation times become larger, but the system remains stable.