Path Integrals

In previous sections we have discussed how to compute thermal correlation functions of classical systems by propagating the usual equations of motion (Hamilton's equations) according to the Velocity-Verlet algorithm. Coordinates and momenta $q(t)$ and $p(t)$ are propagated for a sufficiently long trajectory and correlation functions are obtained as follows:

$$C(t) = \langle A(0) B(t) \rangle = \frac{1}{\tau} \int_0^{\tau} dt' A(q(t'),p(t')) B(q(t'+t),p(t'+t)),$$ (505)

where $A(0)$ and $B(t)$ represent the quantities of interest at time 0 and t, respectively ¹. The goal of this section is to show how to compute $C(t)$ for systems where quantum mechanical effects are important. The quantum mechanical expression of $C(t)$ is,

$$C(t) = Tr[\hat{\rho} \hat{A} \hat{B}(t)],$$ (506)

where $\hat{\rho}=Z^{-1} exp(-\beta \hat{H})$ is the density operator and the operators $\hat{A}$ and $\hat{B}(t)$ are defined so that $A(0) = \langle \Psi_0 \vert \hat{A} \vert \Psi_0 \rangle$ is the expectation value of $A$ at $t=0$ and

$$B(t)=\langle \Psi_0 \vert \hat{B}(t) \vert \Psi_0 \rangle = \lang... ... e^{(i/\hbar) \hat{H} t} \hat{B} e^{-(i/\hbar) \hat{H} t} \vert \Psi_0 \rangle,$$ (507)

is the expectation value of $B$ at time t when the system is initially prepared in state $\vert \Psi_0 \rangle$ and evolves according to the Hamiltonian,

$$\hat{H} = \hat{p}^2/(2m) + \hat{V},$$ (508)

as follows: $\vert \Psi_t \rangle = e^{-(i/\hbar) \hat{H} t} \vert \Psi_0 \rangle$. Note that $\hat{B}(t) = e^{(i/\hbar) \hat{H} t} \hat{B} e^{-(i/\hbar) \hat{H} t}$ is the Heisenberg operator associated with quantity $B$. Thermal correlation functions can therefore be expressed as,

$$C(t) = Z^{-1} Tr[e^{-\beta \hat{H}} \hat{A} e^{(i/\hbar) \hat{H} t} \hat{B} e^{-(i/\hbar) \hat{H} t}],$$ (509)

an expression that can be re-written in coordinate representation as follows:

$$\begin{array}{ll} C(t) = Z^{-1} \int dx \int dx' \int dx'' \i... ...''' \vert e^{-(i/\hbar) \hat{H} t} \vert x \rangle. \end{array}$$ (510)

Note that in order to compute $C(t)$ it is necessary to obtain expressions for the Boltzmann operator matrix elements $\langle x \vert e^{-\beta \hat{H}} \vert x' \rangle$ as well as for the forward and backward time-evolution operator matrix elements $\langle x \vert e^{-(i/\hbar) \hat{H} t} \vert x' \rangle$ and $\langle x \vert e^{(i/\hbar) \hat{H} t} \vert x' \rangle$, respectively. In order to obtain an expression of the matrix elements of the Boltzmann operator, we express the exponential operator as a product of a large number $n$ of exponential operators,

$$\langle x_0 \vert e^{-\beta \hat{H}} \vert x_n \rangle = \langle ... ...lon \hat{H}} e^{-\epsilon \hat{H}} ... e^{-\epsilon \hat{H}} \vert x_n \rangle,$$ (511)

where $\epsilon \equiv \beta/n << 1$. Inserting the closure relation in between exponential operators we obtain,

$$\langle x_0 \vert e^{-\beta \hat{H}} \vert x_n \rangle = \int dx_... ...j=1}^n \int dx_j \langle x_{j-1} \vert e^{-\epsilon \hat{H}} \vert x_j \rangle.$$ (512)

The high-temperature Boltzmann operator $e^{-\epsilon \hat{H}}$ can be written in the form of the Trotter expansion,

$$e^{-\epsilon \hat{H}} \approx e^{-\epsilon \hat{V}/2} e^{-\epsilon \hat{p}^2/(2 m)} e^{-\epsilon \hat{V}/2},$$ (513)

to second order accuracy.


Exercise: Prove that the Trotter expansion, introduced by Eq. (513), is accurate to second-order in $\epsilon$ by substituting the exponential operators according to their expansions in powers of $\epsilon$.


Therefore matrix elements of the Boltzmann operator at high-temperature can be obtained as follows:

$$\langle x_0 \vert e^{-\epsilon \hat{H}} \vert x_1 \rangle = \int ... ...le p \vert x \rangle \langle x \vert e^{-\epsilon \hat{V}/2} \vert x_1 \rangle,$$ (514)

where

$$\langle x \vert p \rangle = \frac{1}{\sqrt{2 \pi \hbar}} e^{\frac{i}{\hbar} x p},$$ (515)

since

$$-i \hbar \frac{\partial}{\partial x} \langle x \vert p \rangle = p \langle x \vert p \rangle.$$ (516)

Furthermore,

$$\langle x \vert e^{-\epsilon \hat{V}/2} \vert x' \rangle = e^{-\epsilon V(x)/2} \delta(x-x').$$ (517)

Therefore,

$$\begin{array}{ll} \langle x_0 \vert e^{-\epsilon \hat{H}} \ve... ...{\hbar} x p} \delta (x-x_1) e^{-\epsilon V(x_1)/2}, \end{array}$$ (518)

which gives,

$$\langle x_0 \vert e^{-\epsilon \hat{H}} \vert x_1 \rangle = \frac... ...[V(x_0)+V(x_1)]} \int dp e^{-\epsilon p^2/(2 m) + \frac{i}{\hbar} (x_0-x_1) p},$$ (519)

or,

$$\langle x_0 \vert e^{-\epsilon \hat{H}} \vert x_1 \rangle = \frac... ...^{-\frac{1}{2} m \left [ \frac{(x_1-x_0)}{\hbar \epsilon} \right ]^2 \epsilon}.$$ (520)

Matrix elements of the Boltzmann operator at finite-temperature can be obtained by substituting Eq. (520) into Eq. (512):

$$\langle x_0 \vert e^{-\beta \hat{H}} \vert x_n \rangle = \int dx_... ...=1}^n \frac{1}{2}[V(x_j)-V(x_{j-1})] + \frac{1}{2} m \omega^2 (x_j-x_{j-1})^2},$$ (521)

where $\omega=1/(\hbar \sqrt{\epsilon})$. Note that the r.h.s of Eq. (521) corresponds to the partition function of a chain of $n$-harmonic oscillators with cordinates $x_j$ under the influence of an external potential $V(x_j)$. Each chain of harmonic oscillators describes a path from $x_0$ to $x_n$. The multidimentional integral, introduced by Eq. (521), can be computed by importance sampling Monte Carlo by sampling sets of coordinates $x_1, ..., x_{n-1}$ with sampling functions defined by the Gaussians associated with the linked harmonic oscillators. Such a computational approach for obtaining thermal equilibrium density matrices is called Path Integral Monte Carlo.


Exercise: Compute $\langle x_0 \vert e^{-\beta \hat{H}} \vert x_n \rangle$ for the Harmonic oscillator defined by the Hamiltonian

$$\hat{H} = \frac{\hat{p}^2}{2 m} + \frac{1}{2} m \omega_0^2 \hat{x}^2,$$ (522)

by using the Path Integral Monte Carlo method, with $n=2,4,6,8$ and 10 and show that for larger values of $n$ the calculation converges to the analytic expression:

$$\langle x \vert e^{-\beta \hat{H}} \vert x' \rangle = \sqrt{\frac... ...eta \hbar)} \left [ (x^2 + x'^2) cosh(\omega \beta) - 2 x x' \right ] \right ),$$ (523)

which in the free particle limit ( $\omega \rightarrow 0$) becomes

$$\langle x \vert e^{-\beta \hat{H}} \vert x' \rangle = \sqrt{\frac... ...} exp \left ( - \frac{m}{2 \beta \hbar^2} \left [ (x - x')^2 \right ] \right ),$$ (524)

since $sinh(\beta \hbar \omega) \rightarrow \beta \hbar \omega$ and $cosh(\beta \hbar \omega) \rightarrow 1$.


Matrix elements of the time-evolution operator $e^{-\frac{i}{\hbar} \hat{H} t}$ can be obtained by following the same methodology implemented for the Boltzmann matrix $e^{-\beta \hat{H} \tau}$. We first introduce the variable substitution $\epsilon \equiv i \tau/\hbar$ in Eq. (520) and we obtain the short-time propagator as follows:

$$\langle x \vert e^{-\frac{i}{\hbar} \hat{H} \tau} \vert x' \rangl... ...eft [ \frac{(x-x')}{\tau} \right ]^2 - \frac{1}{2}[V(x)+V(x')] \right ) \tau }.$$ (525)

Then, we concatenate the short-time propagators introduced by Eq. (525) and we obtain the finite-time propagator,

$$\langle x_0 \vert e^{-\frac{i}{\hbar} \hat{H} t} \vert x_n \rangl... ...} \right ] ^2 - \frac{1}{2} \left [ V(x_j)+V(x_{j-1}) \right ] \right ) \tau },$$ (526)

which in the limit when $\tau \rightarrow 0$ and $n \rightarrow \infty$ with $t = n \tau$ becomes,

$$\langle x_0 \vert e^{-\frac{i}{\hbar} \hat{H} t} \vert x_n \rangle = \int \mathfrak{D}[x(t)] e^{\frac{i}{\hbar} S_{c}(t)},$$ (527)

where $S_c(t)$ is the classical action associated with the arbitrary trajectory $x(t)$,

$$S_c(t') \equiv \int_0^{t'} dt \hspace{.2cm} \left [ \frac{1}{2} m \left ( \frac{\partial}{\partial t}x(t) \right )^2 - V(x(t)) \right ],$$ (528)

and $\mathfrak{D}[x(t)]$ is defined as follows,

$$\int \mathfrak{D}[x(t)] f(x(t)) \equiv \int dx_1 ... \int dx_{n-1} \left ( \frac{m}{2 \pi \hbar i \tau} \right )^{n/2} f(x(t)),$$ (529)

representing the integral over all paths $x(t)$ from $x_0$ to $x_n$, with intermediate coordinates $x_1, x_2, ..., x_{n-1}$ at times $\tau , 2 \tau, ..., (n-1) \tau$, respectively.