next up previous
Next: 7.3 Wave Packet Dynamics Up: 7. Quantum Mechanical Simulation Previous: 7.1 Diffusion Monte Carlo


7.2 Path Integral Monte Carlo (PIMC)

So far: only ground state, i.e. $T=0$. $\Longrightarrow$How to treat quantum systems at finite temperature?

Let a quantum system be in contact with a heat bath of temperature $kT>0$; it is then in a mixed state made up of the energy eigenstates:
\begin{displaymath}
\Psi=\sum_{n}c_{n}\Psi_{n}\,,\;\;\;\;{\rm where}\;\;
\mbox{\rm H} \Psi_{n}=E_{n}\,\Psi_{n}
\end{displaymath} (7.17)

In classical statistical mechanics the probability density of states is given by the Boltzmann factor. In quantum mechanics the density matrix serves the same purpose:
$\displaystyle \rho(\mbox{$\bf r$},\mbox{$\bf r$}';kT)$ $\textstyle \equiv$ $\displaystyle \sum_{n}\Psi_{n}^{*}(\mbox{$\bf r$})\,
e^{\textstyle - \mbox{\rm H}/kT}\,\Psi_{n}(\mbox{$\bf r'$})$  
  $\textstyle =$ $\displaystyle \sum_{n}\Psi_{n}^{*}(\mbox{$\bf r$})\,e^{\textstyle -E_{n}/kT}\,\Psi_{n}(\mbox{$\bf r'$})$ (7.18)

The average of some observable $a(\mbox{$\bf r$})$ is given by
\begin{displaymath}
\langle a \rangle =\int a(\mbox{$\bf r$})\, \rho(\mbox{$\bf ...
...};\beta)\, d\mbox{$\bf r$} \equiv Sp[a \rho]
\right/ Sp[\rho]
\end{displaymath} (7.19)

where $\beta \equiv 1/kT$.

Problem: The explicit form of $\rho(\mbox{$\bf r$},\mbox{$\bf r$};\beta)$ is usually complex or unknown.
Solution: Reduce the above expression to one that contains only known and tractable density matrices.

Back to school: What are the explicit forms of the density matrix for very simple models?

Strategy of PIMC: Express the density matrix of any given system in terms of the free particle density 7.24. This can be done because:
$\displaystyle \rho(x,x';\beta)$ $\textstyle =$ $\displaystyle \sum_{n}\Psi_{n}^{*}(x)\,e^{-\beta \mbox{\rm H}} \Psi_{n}(x')=$  
  $\textstyle =$ $\displaystyle \sum_{n}\Psi_{n}^{*}(x)\,e^{-\beta \mbox{\rm H}/2}\,
e^{-\beta \mbox{\rm H}/2} \Psi_{n}(x')=$  
  $\textstyle =$ $\displaystyle \sum_{n}\Psi_{n}^{*}(x)\,e^{-\beta \mbox{\rm H}/2} \int dx'' \delta(x'-x'')\,
e^{-\beta \mbox{\rm H}/2} \Psi_{n}(x'')=$  
  $\textstyle =$ $\displaystyle \sum_{n}\Psi_{n}^{*}(x)\,e^{-\beta \mbox{\rm H}/2} \int dx'' \sum_{m}
\Psi_{m}(x') \, \Psi_{m}^{*}(x'')\,e^{-\beta \mbox{\rm H}/2} \Psi_{n}(x'')=$  
  $\textstyle =$ $\displaystyle \int dx'' \left[ \sum_{n}\Psi_{n}^{*}(x)\,e^{-\beta \mbox{\rm H}/...
...ft[ \sum_{m}\Psi_{m}^{*}(x'') \, e^{-\beta \mbox{\rm H}/2}
\Psi_{m}(x') \right]$  

The density matrix may therefore be written as
\begin{displaymath}\fbox{$ \displaystyle
\rho(x,x';\beta)=\int dx'' \rho(x,x''; \frac{\beta}{2})\,
\rho(x'',x';\frac{\beta}{2})
$}
\end{displaymath} (7.29)

The integrand in the path integral on the r.h.s. consists of two density matrices at $\beta/2$, i.e. double the original temperature.

Now we iterate, writing $x_{0},x_{1},x_{2}\dots$ instead of $x,x',x'',\dots$:
\begin{displaymath}
\rho(x_{0},x_{P};\beta)=\int \dots \int dx_{1}\, dx_{2}\dots...
..._{1};\frac{\beta}{P})\dots
\rho(x_{P-1},x_{P};\frac{\beta}{P})
\end{displaymath} (7.30)



The number $P$ of intermediate steps is called the Trotter number. If $P$ is large ($5-100$) the high temperature $PT$ will act to damp the effect of the potential $U(x)$:
$\displaystyle \rho(x_{p}, x_{p+1}; \frac{\beta}{P})$ $\textstyle \longrightarrow$ $\displaystyle \sum_{n}\Psi_{n}^{*}(x_{p})\,
e^{\textstyle -(\beta/P)\left[\mbox{\rm H}_{free}+U(x) \right]}\,
\Psi_{n}(x_{p+1})$  
    $\displaystyle \approx \sum_{n}\Psi_{n}^{*}(x_{p})\,
e^{\textstyle -(\beta/P)\mb...
...} \Psi_{n}(x_{p+1})
e^{\textstyle -(\beta/2P)\left[U(x_{p})+U(x_{p+1}) \right]}$  
    $\displaystyle = \rho_{0}(x_{p},x_{p+1};\frac{\beta}{P})\,
e^{\textstyle -(\beta/2P)\left[U(x_{p})+U(x_{p+1}) \right]}$ (7.31)

with diagonal element
\begin{displaymath}
\rho(x_{0},x_{0};\beta)= A^{P/2} \int \dots \int dx_{1} \dots dx_{P-1}\,
e^{\textstyle -\beta (U_{int}+U_{ext})}
\end{displaymath} (7.32)

where
\begin{displaymath}
A \equiv \frac{mP}{2\pi \beta \hbar^{2}}\,,\;\;\;\; U_{ext} \equiv
\sum_{p=0}^{P-1} U(x_{p})/P
\end{displaymath} (7.33)

and
\begin{displaymath}
U_{int} \equiv \frac{A \pi}{\beta} \sum_{p=0}^{P-1} \left(x_{p}-x_{p+1}
\right)^{2}
\end{displaymath} (7.34)



Three dimensions:

\begin{displaymath}
\rho(\mbox{$\bf r$}_{0},\mbox{$\bf r$}_{0};\beta)= A^{3P/2} ...
...mbox{$\bf r$}_{P-1}\,
e^{\textstyle -\beta (U_{int}+U_{ext})}
\end{displaymath} (7.35)

with the same $A$ as above, and
\begin{displaymath}
U_{ext} \equiv \sum_{p=0}^{P-1} U(\mbox{$\bf r$}_{p})/P\,,\;...
...ft\vert\mbox{$\bf r$}_{p}-\mbox{$\bf r$}_{p+1} \right\vert^{2}
\end{displaymath} (7.36)





Visualization / Classical isomorphism: The expression 7.35 looks like a classical Boltzmann factor for a ring polymer of $P$ links [CHANDLER 81]:
Each link is under the influence of an external potential
\begin{displaymath}
u_{ext}(x_{p})=U(x_{p})/P
\end{displaymath} (7.37)

while successive links are coupled by a harmonic bond potential
\begin{displaymath}
u_{int}(x_{p},x_{p+1})=\frac{A \pi}{\beta}\,
\left(x_{p}-x_{p+1} \right)^{2}
\end{displaymath} (7.38)

with $x_{P}=x_{0}$.
Figure 7.2: Classical isomorphism for one particle
\begin{figure}\includegraphics[width=240pt]{figures/f7pimc1.ps}\end{figure}

Therefore: To solve the quantum problem, apply Monte Carlo method to the classical ring polymer.

Figure 7.3: PIMC for one particle
\begin{figure}{\bf Path integral Monte Carlo for one particle:}
\\ [12pt]
A part...
...box{$\bf r$}$\ remain unchanged.
\item Return to 1.
\end{enumerate}
\end{figure}


Note:
- The oscillator strength $k=Pm/\beta^{2} \hbar^{2}$ increases with larger Trotter numbers $P$, requiring smaller MC steps
- The external potential $u_{ext}(\mbox{$\bf r$}_{p})=U(\mbox{$\bf r$}_{p})/P$ is damped by larger $P$, allowing for larger steps
$\Longrightarrow$Move the entire ring polymer by a large random step, then displace the individual elements by a small amount.
Or:
Displace the center of mass of the chain by a wide random step, then construct the entire ring polymer anew, sampling the single element positions from the narrow multivariate Gauss distribution (see Chapter ``Stochastics''):
\begin{displaymath}
p(\mbox{$\bf r$}_{0},\dots \mbox{$\bf r$}_{P-1}) \propto exp...
...=0}^{P-1}(\mbox{$\bf r$}_{p}-\mbox{$\bf r$}_{p+1})^{2}
\right]
\end{displaymath} (7.39)



$N$ interacting quantum objects:
- Assume a pair potential $u(\vert\mbox{$\bf r$}_{j}-\mbox{$\bf r$}_{i}\vert)$.
- Each quantum particle is represented by a classical $P$-element ring chain.
- Let $\mbox{$\bf r$}_{i,p}$ be the position of element $p\,$ in chain (=particle) $i\,$
- The diagonal element of the total density matrix is then
$\displaystyle \rho(\mbox{$\bf r$}_{0},\mbox{$\bf r$}_{0};\beta)$ $\textstyle =$ $\displaystyle A^{3NP/2} \int \dots \int d\mbox{$\bf r$}_{1,1} \dots d\mbox{$\bf...
...}^{P-1}
\left( \mbox{$\bf r$}_{i,p}- \mbox{$\bf r$}_{i,p+1} \right)^{2} \right]$  
    $\displaystyle \cdot exp \left[- \frac{\beta}{P}\sum_{i=1}^{N}\sum_{j>i}\sum_{p=0}^{P-1}
u(\vert\mbox{$\bf r$}_{j,p}-\mbox{$\bf r$}_{i,p}\vert) \right]$ (7.40)

with $\mbox{$\bf r$}_{0}\equiv \{ \mbox{$\bf r$}_{1,0} \dots \mbox{$\bf r$}_{N,0} \}$. The pair potential acts only between respective links ($p$) of different chains.

EXERCISE: Write a PIMC program treating the case of one particle of mass $m$ in a two-dimensional oscillator potential $U(\mbox{$\bf r$})=kr^{2}/2$. Let the Trotter number vary between $2$ and $10$. Determine the positional probability $p(\mbox{$\bf r$})$ of the particle from the relative frequency of residence at $r$, averaged over all chain links. Noting that
\begin{displaymath}
p(\mbox{$\bf r$}) \equiv \rho(\mbox{$\bf r$}, \mbox{$\bf r$};\beta)
\end{displaymath} (7.41)

we would expect for the two-dimensional harmonic oscillator (with $\omega_{0}^{2}=k/m$)
\begin{displaymath}
p(r)=2\pi r \left[ \frac{A}{\pi}\right] \, e^{-Ar^{2}}\,,\;\...
... \omega_{0}^{2}}{\hbar}\,
tanh\frac{\beta \hbar \omega_{0}}{2}
\end{displaymath} (7.42)

(For convenience, put $\hbar =1$.) Draw several configurations of the ring polymer that occur in the course of the simulation.


Sample applications:
next up previous
Next: 7.3 Wave Packet Dynamics Up: 7. Quantum Mechanical Simulation Previous: 7.1 Diffusion Monte Carlo
Franz J. Vesely Oct 2005
See also:
"Computational Physics - An Introduction," Kluwer-Plenum 2001