User:Statmech2020/sandbox

Path integral Monte Carlo (PIMC) is a quantum Monte Carlo method in the path integral formulation of quantum statistical mechanics. In this approach, each quantum-mechanical particle is approximated by a path of particles ("beads") obeying an effective classical Hamiltonian. In the limit of infinitely many beads, the distribution of classical beads approximates the true quantum-mechanical distribution.

The equations are often applied assuming that quantum exchange does not matter (often, the particles are assumed to be Boltzmann particles, not the physically realistic fermion and boson particles). The theory usually is applied to calculate thermodynamic properties such as the internal energy, heat capacity, or free energy. As with all Monte Carlo method based approaches, a large number of points must be calculated. As more "replicas" are used to integrate the path integral, the more quantum and the less classical the result is. But, the answer might become less accurate initially as more beads are added, until a point where the method starts to converge to the correct quantum answer. Because it is a statistical sampling method, PIMC takes into account all the anharmonicity, and because it is quantum, it takes into account all quantum effects (with the exception of the exchange interaction usually). An early application was to the study of liquid helium. It has been extended to include the grand canonical ensemble and the microcanonical ensemble.

With agent-based PIMC the perimeter and sum borderlines of objects can be calculated.

Formulation
Following the treatment presented in, we briefly describe the theoretical formulation of path integral Monte Carlo. For simplicity, we consider a single 1D particle in a potential $$V(x)$$.

Recall the definition of the density matrix:

$$\rho = e^{-\beta H}$$

where $$\beta = 1/k_B T$$ is the thermodynamic beta (inverse temperature).

We can expand the density matrix in an eigenbasis of the Hamiltonian

$$\rho = \sum_i |\phi_i \rangle e^{-\beta E_i} \langle \phi_i |.$$

This allows us to compute the average of an observable of interest using the trace of the density matrix and the operator:

$$\langle A \rangle = \frac{1}{Z} \text{Tr}\rho A= \frac{\sum_i e^{-\beta E_i}\langle \phi_i| A|\phi_i\rangle}{\sum_i e^{-\beta E_i}} $$

where $$Z= \text{Tr}\,\rho$$ is the partition function. Similar to the how the time evolution is factorized in the real-time path integral, path integral Monte Carlo relies on dividing the density matrix into many factors. In particular, the following relationship holds exactly:

$$\rho(\beta_1+\beta_2) = e^{-(\beta_1+\beta_2)H}=e^{-\beta_1 H} e^{-\beta_2 H} = \rho(\beta_1) \rho(\beta_2)$$

This relates the density matrix at high $$\beta$$ (low temperature) to products of density matrices at low $$\beta$$ (high temperature). The partition function can be written (expanding the trace in the position basis):

$$Z=\text{Tr}\,\rho = \int dx\, \langle x | e^{-\beta H} |x\rangle.$$

To evaluate this integral, we assume that the Hamiltonian is separable $$H=T(p) + V(x)$$, where $$p$$ and $$x$$ follow the usual canonical commutation relation $$[x, p] = i\hbar$$. Because $$T $$ and $$V $$

do not commute, $$\rho=e^{-\beta(T + V)} \neq e^{-\beta T} e^{-\beta V}$$. However, we can use the Trotter product formula to factorize the density matrix:

$$e^{-\beta(T + V)} = \lim_{M\rightarrow \infty} \left(e^{-\beta T/M} e^{-\beta V/M}\right)^M$$

Now, using this fact, we factorize the density matrix into $$M$$ terms.

$$\rho(\beta) = \left[\rho\left(\frac{\beta}{M}\right)\right]^M = \left[e^{-\beta(T + V)/M}\right]^M \approx \left[e^{-\beta T/M}e^{-\beta V/M}\right]^M$$

where the final approximation is based on the Trotter product formula (assuming $$M$$ is large). Now we insert resolutions of the identity and sum over other positions $$x_2,.\ldots,x_M$$:

$$\rho(\beta) \approx \int dx_2\ldots dx_M \, e^{-\beta T/M} e^{-\beta V / M} |x_2\rangle\langle x_2|e^{-\beta T/M} e^{-\beta V / M} |x_3\rangle\langle x_3|e^{-\beta T/M} e^{-\beta V / M}|x_4\rangle\langle x_4|\cdots |x_M\rangle\langle x_M|e^{-\beta T/M} e^{-\beta V / M}$$

Because $$V$$ is diagonal in the position basis, this gives $$\rho(\beta) \approx \int dx_2\ldots dx_M \, e^{-\frac{\beta}{M}\left[V(x_2)+\cdots+V(x_M)\right]}e^{-\beta T/M} |x_2\rangle\langle x_2|e^{-\beta T/M} |x_3\rangle\langle x_3|e^{-\beta T/M}|x_4\rangle\langle x_4|\cdots |x_M\rangle\langle x_M|e^{-\beta T/M} e^{-\beta V / M}$$

Factors like $$\langle x'| e^{-\beta T / M} |x''\rangle$$ correspond to a free particle propagator, which can be shown to take the form

$$\langle x'|e^{-\beta T / M}|x''\rangle = \left(\frac{M}{4\pi \lambda \beta}\right)^{1/2} \exp{\left(-\frac{(x' - x'')^2}{4\lambda \beta / M}\right)}$$

When computing expectation values of operators, we care about traces of products of the density matrix and observables. For instance, the simplest example is the partition function:

$$Z = \text{Tr}\,\rho = \int dx_1 \langle x_1 | \rho(\beta) |x_1 \rangle =

\int dx_1 dx_2 \ldots dx_M\,

\exp\left(-\sum_{m=1}^M \left[\frac{(x_{m-1} - x_m)^2}{4\lambda \beta/M} + \frac{\beta}{M} V(x_m)\right]\right)$$

where by notation we take $$x_0=x_M$$. Thus, the term $$(x_0 - x_1)^2$$ connects the initial position of the path back to the start (i.e., $$(x_M - x_1)^2$$); the path is a "ring."

The integrand can be interpreted as a Boltzmann weight $$e^{-\beta H_\text{eff}}$$ for an effective classical Hamiltonian of the form

$$H_\text{eff} =

\sum_{m=1}^M \left[\frac{(x_{m-1} - x_m)^2}{4\lambda \beta^2/M} + \frac{1}{M} V(x_m)\right] $$

This is like the Hamiltonian for a classical ring polymer consisting of point particles joined by harmonic springs in a potential $$V(x)$$. Thus, in the limit of large $$M$$, we can approximate the low-temperature quantum distribution by the distribution of the classical ring polymer.

Because the Boltzmann weight is positive, it corresponds to a valid un-normalized probability measure. To estimate quantities of interest, we can sample from this density using techniques like Metropolis-Hastings sampling. This allows us to compute estimated expected values of relevant observables like the energy (Figure 1).