Schrodinger’s equations states that \begin{align} \frac{d}{dt}|\Psi_t \rangle = \frac{-i \cdot \hat H}{\hbar} \cdot |\Psi_t \rangle \end{align} A formal solution to the equation is \begin{align} |\Psi_t \rangle = e^{\frac{-i \cdot \hat Ht}{\hbar}} \cdot |\Psi_0 \rangle \end{align}
Using \begin{align} 1 = \int{|x_i\rangle \langle x_i| dx_i} \end{align} and insert \(\langle x_f|\) to the left, we obtain \begin{align} \langle x_f| \Psi_t \rangle = \int \langle x_f | e^{\frac{-i \cdot \hat Ht}{\hbar}} | x_i \rangle \langle x_i | \Psi_0 \rangle dx_i \end{align}
Define \(U(x_f, t; x_i, 0) = \langle x_f | e^{\frac{-i \cdot \hat Ht}{\hbar}} | x_i \rangle\), then we have \begin{align} \langle x_f| \Psi_t \rangle = \int U(x_f, t; x_i, 0) \cdot \langle x_i |\Psi_0 \rangle dx_i \end{align}
How do we calculate \(U(x_f, t; x_i, 0)\)? Feynman came up with the following solution. Let’s descretize the path from \(x_i\) to \(x_f\). We have points along the path, \(x_1, x_2, ..., x_{N - 1}\), and \(x_0 = x_i , x_N = x_f\). For each point we have
\begin{align}
1=\int{|x_{1}\rangle \langle x_{1}| dx_{1}} \
1=\int{|x_{2}\rangle \langle x_{2}| dx_{2}} \
….
\end{align}
We insert all these identity in \(U(x_f, t; x_i, 0)\) and obtain \begin{align} U(x_f, t; x_i, 0)=\int dx_{N-1} \int dx_{N-2} \int dx_{N-3}…\int dx_1 \langle x_N | e^{\frac{-i \cdot \hat Ht}{N\hbar}} | x_{N-1} \rangle \cdot \langle x_{N-1} | e^{\frac{-i \cdot \hat Ht}{N\hbar}} | x_{N-2} \rangle \cdot \langle x_{N-2} | e^{\frac{-i \cdot \hat Ht}{N\hbar}} | x_{N-3} \rangle \cdot … \cdot \langle x_1 | e^{\frac{-i \cdot \hat Ht}{N\hbar}} | x_{0} \rangle \end{align}
To evaluate \(\langle x_N | e^{\frac{-i \cdot \hat Ht}{N\hbar}} | x_{N-1} \rangle\), we insert the following identity
\begin{align}
1 = \int |p\rangle \langle p | dp \
1 = \int |p’ \rangle \langle p’ | dp’
\end{align}
So,
\begin{align}
\langle x_N | e^{\frac{-i \cdot \hat Ht}{N\hbar}} | x_{N-1} \rangle = \int dp \int dp’ \langle x_N | p’ \rangle \langle p’ |e^{\frac{-i \cdot \hat Ht}{N\hbar}}|p \rangle \langle p|x_{N-1} \rangle
\end{align}
For larger enough N, in interval \([x_{N-1}, x_N]\), \(V(x)\) can be approximated as \(V(x_{N-1})\), so \begin{align} \langle p’ |e^{\frac{-i \cdot \hat Ht}{N\hbar}}|p \rangle = \langle p’ |e^{\frac{-i \cdot t }{N\hbar}(\frac{\hat p^2}{2 \cdot m} + V(x_{N-1}))}|p \rangle = e^{\frac{-i \cdot t }{N\hbar}(\frac{\hat p^2}{2 \cdot m} + V(x_{N-1}))} \cdot \delta (p - p’) \end{align}
We also have
\begin{align}
\langle x_N | p’ \rangle = \frac{1}{\sqrt{2 \cdot \pi \hbar}} \cdot e^{\frac{i \cdot p’ \cdot x_N}{\hbar}} \
\langle p | x_{N-1} \rangle = \frac{1}{\sqrt{2 \cdot \pi \hbar}} \cdot e^{-\frac{i \cdot p \cdot x_{N- 1}}{\hbar}}
\end{align}
Therefore, we have
\begin{equation}
\langle x_N | e^{\frac{-i \cdot \hat Ht}{N\hbar}} | x_{N-1} \rangle = \frac{1}{2\pi\hbar}\int dp e^{\frac{-i \cdot t }{N\hbar}(\frac{\hat p^2}{2 \cdot m} + V(x_{N-1}))} \cdot e^{\frac{i}{\hbar} \cdot p \cdot (x_N - x_{N-1})} = \frac{1}{2\pi\hbar} \cdot e^{\frac{-i\Delta t}{\hbar} V(x_{N-1})} \cdot e^{\frac{i\Delta t}{2m\hbar}(\frac{(x_N - x_{N - 1})m}{\Delta t})^2} (\frac{2m\pi\hbar}{i\Delta t})^{1/2} \
\langle x_N | e^{\frac{-i \cdot \hat Ht}{N\hbar}} | x_{N-1} \rangle = (\frac{m}{2\pi\hbar i \Delta t})^{1/2} \cdot e^{\frac{i \Delta t}{\hbar} \cdot (\frac{1}{2} m (\frac{x_N - x_{N - 1}}{\Delta t})^2 - V(x_{N - 1}))} = (\frac{m}{2\pi\hbar i \Delta t})^{1/2} \cdot e^{\frac{i \Delta t}{\hbar} \cdot L(x_{N - 1})}
\label{eq:pathint}
\tag{1}
\end{equation}
where \(\Delta t = \frac{t}{N}\) and \(L\) is the Lagrangian
Now we insert the above formula for \(j = 0...N - 1\) into the expression of \(U(x_f, t; x_i, 0)\) and obtain,
\begin{align}
U(x_f, t; x_i, 0)=\int dx_{N-1} \int dx_{N-2} \int dx_{N-3}…\int dx_1 (\frac{m}{2\pi\hbar i \Delta t})^{N/2} \cdot e^{\frac{i\Delta t}{\hbar}\sum_{j=0}^{N - 1}L(x_{j})} \
U(x_f, t; x_i, 0)=\int dx_{N-1} \int dx_{N-2} \int dx_{N-3}…\int dx_1 (\frac{m}{2\pi\hbar i \Delta t})^{N/2} \cdot e^{\frac{i}{\hbar} \int_0^t L d \tau} \
U(x_f, t; x_i, 0)=\int dx_{N-1} \int dx_{N-2} \int dx_{N-3}…\int dx_1 (\frac{m}{2\pi\hbar i \Delta t})^{N/2} \cdot e^{\frac{i}{\hbar} S(t)}
\end{align}
Where \(S(t)\) is the Action Integral
Insert the above expression for \(U\) into the expression for \(\langle x_f| \Psi_t \rangle\), we obtain \begin{align} \langle x_f| \Psi_t \rangle = \int dx_0 \int dx_1 …\int dx_{N - 1} (\frac{m}{2\pi\hbar i \Delta t})^{N/2} \cdot e^{\frac{i}{\hbar} S(t)} \cdot \langle x_0|\Psi_0 \rangle \end{align}
This equation is Exact and it is called Feynman’s Path Integral Formulation of Quantum Mechanics
But why we spend so much time derive the Path Integral Formulation? How it is related to Ring Polymer Molecular Dynamics or Statistical Mechanics?
To see the connection with statistic mechanics, we observe that the quantum statistical partition function \begin{align} Q(\beta) = Tr(e^{-\beta \hat H}) \end{align} where \(\beta = \frac{1}{k_B T}\) and \(Tr\) is the trace of the matrix. If we expand Q in the basis of coordinate states, we obtain \begin{align} Q(\beta) = \int dx \langle x | e^{-\beta \hat H} | x \rangle \end{align}
We notice that \begin{align} \rho_{xx} = \langle x | e^{-\beta \hat H} | x \rangle \end{align} is very close to \begin{align} U(x_f, t; x_i, 0) = \langle x_f | e^{\frac{-i \cdot \hat Ht}{\hbar}} | x_i \rangle \end{align} in the form with \(x_f = x_i\) and \(t = -i\beta\hbar\)
So, doing equilibrium statistical mechanics is equivalent to doing quantum mechanics with imaginary time.
Similarly, we can use Feynman’s Path Integral Formulation to calculate \(\rho_{xx}\) too. In order to do that, we just simply replace \(x_f\) with \(x_0\)and t with \(-i\beta\hbar\). Substitue t in Eq. \eqref{eq:pathint}.
So, \begin{align} \langle x_N | e^{\frac{-\beta \cdot \hat H}{N}} | x_{N-1} \rangle = (\frac{m}{2\pi\beta_{N}\hbar^2})^{1/2} \cdot e^{-\beta_N V(x_{N - 1})} \cdot e^{-\frac{\beta_N}{2}m\omega_N^2 \cdot (x_N - x_{N - 1})^2} \end{align} Where \(\beta_N = \frac{\beta}{N}\) and \(\omega_N = \frac{1}{\beta_N \hbar}\).
Inserting above equation to \(\rho_{xx}\) , we obtain \begin{align} \rho_{xx} = \int dx_1 \int dx_2 \int dx_3 … \int dx_{N - 1} (\frac{m}{2\pi\beta_{N}\hbar^2})^{N/2} \cdot e^{-\beta_N \cdot (\sum_0^{N-1} \frac{m}{2} \omega_N^2(x_{j+1} - x_j)^2 + V(x_j))} \end{align} where \(x_N = x_0\)
Therefore, \begin{align} Q(\beta) = \int dx \langle x | e^{-\beta \hat H} | x \rangle = \int dx_0 \int dx_1 … \int dx_{N-1} (\frac{m}{2\pi\beta_{N}\hbar^2})^{N/2} \cdot e^{-\beta_N \cdot (\sum_0^{N-1} \frac{m}{2} \omega_N^2(x_{j+1} - x_j)^2 + V(x_j))} \end{align}
Notice for N distinguishable classical bead with interacting potential \(V_{eff} = (\sum_0^{N-1} \frac{m}{2} \omega_N^2(x_{j+1} - x_j)^2 + V(x_j))\) (Only adjacent beads are connected through a spring with force constant of \(\omega_N\). Its partition function is \begin{align} Q_{classical}(\beta) = (\frac{m}{2\pi\beta\hbar^2})^{N/2} \int dx_0 \int dx_1 … \int dx_{N - 1} e^{-\beta V_{eff}} \end{align}
The form of \(Q_{classical}(\beta)\) is almost the same as \(Q(\beta)\), except the temperature. If we evaluate \(Q_{classical}\) at a different temperature \(T' = NT\), we would obtain \begin{align} Q(\beta) = Q_{classcial}(\beta_N) \end{align}
In summary, the quantum mechanical partition function \(Q(\beta)\) can be obtained by classcial mechanics simulation of N distinguishable beads with effective Hamiltonian as at temperature NT \begin{align} H_{eff} = \sum_0^{N - 1} (\frac{p_j^2}{2m} + \frac{m}{2} \omega_N^2(x_{j+1} - x_j)^2 + V(x_j)) \end{align} Where \(x_N = x_0\).