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 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 are all positive and are ordered so that is the smallest of them, we get an interesting result by replacing the time variable with an imaginary time variable . The function is then

which is a sum of exponentially decaying terms, of which the first one, decays slowest. So, in the limit of large , the wave function is

i.e. after normalizing it, it’s approximately the same as the ground state . Also, if 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 . 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, , .

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 for several values of 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.