 # Numerical solution of PDE:s, Part 5: Schrödinger equation in imaginary time

In the last post, I mentioned that the solution of the time dependent 1D Schrödinger equation can be written by expanding the initial state $\Psi (x,0)$ in the basis of the solutions of time-independent Schrödinger equation and multiplying each term in the expansion by a complex-valued time dependent phase factor: Now, assuming that the energies $E_n$ are all positive and are ordered so that $E_0$ is the smallest of them, we get an interesting result by replacing the time variable $t$ with an imaginary time variable $s = it$. The function $\Psi (x,s)$ is then which is a sum of exponentially decaying terms, of which the first one, $c_0 \exp(-E_0 s)\psi_0 (x)$ decays slowest. So, in the limit of large $s$, the wave function is i.e. after normalizing it, it’s approximately the same as the ground state $\psi_0$. Also, if $s$ is a large number, we have or So, here we have a way to use the TDSE to numerically calculate the ground state wavefunction and corresponding energy eigenvalue for any potential energy function $V(x)$. This is very useful, as the ground state can’t usually be solved analytically, except for some very simple potential energy functions such as the harmonic oscillator potential.

To test this imaginary time method, I will approximately calculate the ground state of an anharmonic oscillator, described by a Schrödinger equation as an initial state, I will choose a Gaussian function and the computational domain is defined by 0 < x < 6, 0 < t < 3, $\Delta x = 0.05$, $\Delta t = 0.01$.

An R-Code that performs this calculation is shown below:

library(graphics) #load the graphics library needed for plotting

lx <- 6.0 #length of the computational domain
lt <- 3.0 #length of the simulation time interval
nx <- 120 #number of discrete lattice points
nt <- 300 #number of timesteps
dx <- lx/nx #length of one discrete lattice cell
dt <- (-1i)*lt/nt #length of imaginary timestep

V = c(1:nx) #potential energies at discrete points

for(j in c(1:nx)) {
V[j] = as.complex(0.25*(j*dx-3)*(j*dx-3)+0.15*(j*dx-3)*(j*dx-3)*(j*dx-3)*(j*dx-3)) #anharmonic potential
}

kappa1 = (1i)*dt/(2*dx*dx) #an element needed for the matrices
kappa2 <- c(1:nx) #another element

for(j in c(1:nx)) {
kappa2[j] <- as.complex(kappa1*2*dx*dx*V[j])
}

psi = as.complex(c(1:nx)) #array for the wave function values

for(j in c(1:nx)) {
psi[j] = as.complex(exp(-2*(j*dx-3)*(j*dx-3))) #Gasussian initial wavefunction
}

xaxis <- c(1:nx)*dx #the x values corresponding to the discrete lattice points

A = matrix(nrow=nx,ncol=nx) #matrix for forward time evolution
B = matrix(nrow=nx,ncol=nx) #matrix for backward time evolution

for(j in c(1:nx)) {
for(k in c(1:nx)) {
A[j,k]=0
B[j,k]=0
if(j==k) {
A[j,k] = 1 + 2*kappa1 + kappa2[j]
B[j,k] = 1 - 2*kappa1 - kappa2[j]
}
if((j==k+1) || (j==k-1)) {
A[j,k] = -kappa1
B[j,k] = kappa1
}
}
}

for (k in c(1:nt)) { #main time stepping loop

sol <- solve(A,B%*%psi) #solve the system of equations

for (l in c(1:nx)) {
psi[l] <- sol[l]
}

nrm = 0
for (l in c(1:nx)) nrm <- nrm + dx*abs(psi[l])*abs(psi[l])

if(k %% 3 == 1) { #make plots of psi(x) on every third timestep
jpeg(file = paste("plot_",k,".jpg",sep=""))
plot(xaxis, abs(psi)^2,xlab="position (x)", ylab="Abs(Psi)^2")
title(paste("|psi(x,t)|^2 at t =",k*dt))
lines(xaxis, abs(psi)^2)
lines(xaxis, V*max(abs(psi)^2))
dev.off()
}
}

And a set of plots of $|\Psi (x,t)|^2$ for several values of $t$ look like this:

Here the shape of the anharmonic potential has been plotted to the same images.

The problem with this method for finding a ground state is that if the system has more degrees of freedom than a single coordinate x, the matrices in the linear systems needed in the time stepping quickly become very large and the calculation becomes prohibitively slow. To make this worse, the Crank-Nicolson method of time stepping can’t be parallelized for multiple processors. However, there is another way to compute the evolution of a wave function in imaginary time, which is called Diffusion Monte Carlo, and that is easily parallelizable. DMC is one practical way for calculating ground state wave functions for multi-particle systems such as a helium or a lithium atom.