# Numerical solution of PDE:s, Part 6: Adiabatic approximation for quantum dynamics

Posted by

Having solved the time-dependent Schrödinger equation both in real and imaginary time, we can move forward to investigate systems where the potential energy function $V$ has an explicit time dependence in it:

In this kind of systems, the expectation value of the Hamiltonian operator doesn’t have to stay constant.

Time-dependent perturbation theory is one method for finding approximate solutions for this kind of problems, but here I will handle a simpler example, which is called adiabatic approximation.

Suppose that the potential energy function $V(x,t)$ is known. Now, let’s say that we also know the solutions of the time-independent Schrödinger equation

for any value of $t$. I denote the solutions as $\psi_n (x;t)$, where it is understood that $x$ is a variable and $t$ is a parameter. Now, if the function $V(x,t)$ changes very slowly as a function of time, i.e. its partial derivative with respect to $t$ is small at all points of the domain, we can use the adiabatic approximation, which says that if the initial state $\Psi (x,0)$ is the ground state for the potential $V(x,0)$, then the state at time $t$ is the ground state for the potential $V(x,t)$.

So, we can change a ground state of one potential $V_1 (x)$ into the ground state of another potential $V_2 (x)&bg=ffffff&fg=000000$ by making a continuous change from $V_1 (x)$ to $V_2 (x)$ slowly enough.

Let’s test this by chooosing a function $V$ as

i.e. a Hookean potential that moves to the positive x-direction with constant speed. If we set the wavefunction at $t=0$ to

which is the ground state corresponding to $V(x,0)$, the time depelopment of the wavepacket should be like

which means that is moves with the same constant speed as the bottom of the potential $V$. This can be calculated with the R-Code below:

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

lx <- 6.0 #length of the computational domain
lt <- 15.0 #length of the simulation time interval
nx <- 100 #number of discrete lattice points
nt <- 300 #number of timesteps
dx <- lx/nx #length of one discrete lattice cell
dt <- lt/nt #length of timestep

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

for(j in c(1:nx)) {
V[j] = as.complex(2*(j*dx-3)*(j*dx-3)) #harmonic 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(-(j*dx-3)*(j*dx-3))) #Gaussian 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 #off-diagonal elements
B[j,k] = kappa1
}
}
}

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

for(j in c(1:nx)) {
V[j] = as.complex(2*(j*dx-3-k*dt*0.05)*(j*dx-3-k*dt*0.05)) #time dependent potential
}

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

for(l in c(1:nx)) {
for(m in c(1:nx)) {
A[l,m]=0
B[l,m]=0
if(l==m) {
A[l,m] = 1 + 2*kappa1 + kappa2[m]
B[l,m] = 1 - 2*kappa1 - kappa2[m]
}
if((l==m+1) || (l==m-1)) {
A[l,m] = -kappa1
B[l,m] = kappa1
}
}
}
sol <- solve(A,B%*%psi) #solve the system of equations

for (l in c(1:nx)) {
psi[l] <- sol[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",ylim=c(0,2))
title(paste("|psi(x,t)|^2 at t =",k*dt))
lines(xaxis, abs(psi)^2)
lines(xaxis, V)
dev.off()
}
}

and in the following sequences of images you see that the approximation is quite good for $v = 0.05$

As an animation, this process looks like shown below:

By doing the same calculation again, but this time with $v = 3$, the image sequence looks like this:

where it is obvious that the approximation doesn’t work anymore.