Monte Carlo solution of the Schrödinger equation

Monte Carlo methods can be used to solve the time-independent Schrödinger equation by iteratively solving an associated integral equation. For the method used in this program we begin by considering the time-dependent Schrödinger equation, i.e.,

$$- {\partial {\phi}( {\bf x},t) \over { i \partial t} } = ( \hbox{{\curly H}}- E_T ) {\phi}( {\bf x},t) .$$ (1)

$E_T$ is an arbitrary energy shift and $\hbar$ = 1 in atomic units. Our motivation is to exploit the formal similarity between Eq. 1 (in imaginary time) and the classical diffusion equation. This will be discussed at greater length in Sec. 2.2. The formal solution of Eq. 1 may be expressed in terms of the eigenfunctions, $\Phi_k$, and eigenvalues, $E_k$, of $\hbox{{\curly H}}$,

$${\phi}( {\bf x}, t ) = \sum_{{k} = 0}^{\infty} C_k \Phi_k ({\bf x}) e^{-i( E_k - E_T ) t} .$$ (2)

The coefficients $C_k$ depend only on the initial conditions

$$C_k = \langle \Phi_k \vert {\phi}( t = 0 ) \rangle .$$ (3)

The oscillatory time behavior of Eq. 2 is not of interest here; rather we wish to extract the ground state. This extraction can be achieved by substituting an imaginary time $\tau$ = $it$ in Eqs. 1 and 2 so that the oscillatory time behavior becomes exponential, i.e.

$$- {\partial {\phi}( {\bf x},\tau) \over \partial\tau} = ( \hbox{{\curly H}}- E_T ) {\phi}( {\bf x},\tau),$$ (4)

which has the solution

$${\phi}( {\bf x}, \tau ) = \sum_{k=0}^{\infty} C_k \Phi_k ({\bf x}) e^{-(E_k - E_T ) \tau } .$$ (5)

For sufficiently large $\tau$, only one eigenfunction contributes to ${\phi}$, namely the one with the most negative eigenvalue -- the ground state. Time-independence can be achieved by choosing $E_T$ to be ${E_0}$, yielding

$${\phi}({\bf x},\tau ) = C_0 \Phi_0 + \sum_{k=1}^{\infty} C_k \Phi_k ({\bf x}) e^{-(E_k - {E_0}) \tau } ~.$$ (6)

Thus, as $\tau \to \infty$, ${\phi}\to C_0 \Phi_0$. Note that even though we shall refer to the motion of particles, no real-time dynamics can be obtained from Eq. 4.

The use of Monte Carlo methods requires an integral equation which can be solved iteratively. The integral form of Eq. 1 is

$${\phi}( {\bf y}, \tau + {\delta \tau}) = \int\limits _{-\inf... ...bf x}; {\delta \tau}) {\phi}( {\bf x}, \tau ) {\rm d} {\bf x}.$$ (7)

The Green's function $G( {\bf y}, {\bf x}; {\delta \tau})$ is a solution of the imaginary-time-dependent Schrödinger equation, i.e.

$$- { \partial G( {\bf y}, {\bf x}; {\delta \tau}) \over \part... ...\hbox{{\curly H}}- E_T ) G( {\bf y}, {\bf x}; {\delta \tau}) .$$ (8)

with the boundary condition,

$$G( {\bf y}, {\bf x}; {\delta \tau}=0 ) = \delta ( {\bf x}- {\bf y}) .$$ (9)

Formally solving for $G( {\bf y}, {\bf x}; {\delta \tau})$ we find that it is the position space representation of the imaginary-time propagator, namely,

$$G( {\bf y}, {\bf x}; {\delta \tau}) = \langle {\bf y}\vert ... ...(\hbox{{\smcurly H}}-E_T) {\delta \tau}}\vert {\bf x}\rangle .$$ (10)

The convergence of the series produced by iterating Eq. 7 can be studied by expanding $G( {\bf y}, {\bf x}; {\delta \tau})$ in $\Phi_i$, the eigenfunctions of $\hbox{{\curly H}}$. This expansion is achieved by inserting two complete sets of states into Eq. 10, and yields

$$G( {\bf y}, {\bf x}; {\delta \tau})= \sum_{i=0}^{\infty} e^{-(E_i - E_T ) {\delta \tau}} \Phi_i ({\bf y}) \Phi_i ({\bf x}) .$$ (11)

Now substitute this into Eq. 7 while also expanding ${\phi}({\bf x},0)$ = $\sum_k C_k \Phi_k$, to obtain the first iteration,

$${\phi}( {\bf y}, {\delta \tau}) = \int\limits _{-\infty}^{\infty} \sum_{i=0}^{\infty} e^{-(E_i - E_... ...f y}) \Phi_i ({\bf x}) \sum_{k=0}^{\infty} C_k \Phi_k ({\bf x}) {\rm d} {\bf x}$$

$$= \sum_{k=0}^{\infty} C_k \Phi_k ({\bf y}) e^{ - ( E_k - E_T ) {\delta \tau}} .$$ (12)

This shows that the excited states decay exponentially fast. After $n$ iterations we have,

$${\phi}( {\bf y}, n {\delta \tau}) = \sum_{k=0}^{\infty} C_k \Phi_k ({\bf y}) e^{ - ( E_k - E_T ) n {\delta \tau}} .$$ (13)

At large $n$, the lowest energy solution with a non-zero coefficient $C_k$ will dominate the sum. Generally this will be the ground state, unless ${\phi}({\bf x},0)$ is specifically chosen so that it is orthogonal to the ground state.

The ``physical'' interpretation of the time-dependent Green's function is that $G({\bf y}, {\bf x}; \tau)$ represents the probability that a particle moves from ${\bf x}$ to ${\bf y}$ in an imaginary time $\tau$. Hence, starting from an initial ensemble of random walkers and propagating them iteratively in imaginary time according to the probabilities $G({\bf y}, {\bf x}; \tau)$, the population density after a large number of generations will represent the ground state wave function.