Time dependent simulations

I will focus on the pseudo-spectral Fourier split operator (or split step) method.

Firsly, we need to quantised position space by making a grid over a specific range of position values. Assume $-10$ $\mu$m to $+10$ $\mu$m, giving a grid divided into $xDim=2^8$ equispaced elements. Let us call this grid $\mathbf{x}$, and the maximum value $x_\textrm{max}$. The following Julia code will carry out the above:

xDim = 2^8; # The resolution of your grid.
xMax = 4e-4;
x = linspace(-xMax, xMax, xDim);
dx = abs(x[2]-x[1]);

Great, so now that we have our grid in position space, we need to also consider generating another grid, this time in momentum space (also known as reciprocal space, $k$ space). Position and momentum space are related by the Fourier transform.

kxMax = (2*pi/xMax)*(xDim/2);
kx = circshift(linspace(-kxMax,kxMax,xDim),xDim/2);
dkx = abs(kx[2] - kx[1]);

However, as you may know, knowledge of position can affect your knowledge of momentum, and vice-versa. In this case, our fundamental limit is defined by the size of our grid dimension, $d$. Just as the wavenumber is given by the relation $k = \frac{2\pi}{\lambda}$, we can obtain our smallest $k$ value using our largest position space value, $x_\textrm{max}$.

#Define properties of the system
omegaX = 2*pi*1 # Harmonic trapping potential frequency
hbar = 1.0545718e-34;
m = 1.4431607e-25; #Rb87 mass
a_s = 4.76e-9; #S-wave scattering length
dt = 1e-4; #Timestep

#Make an initial guess
wfc = pdf(Normal(0,xMax/100),x) + 0*im;

#Define the trapping potential and kinetic energy operators
V = 0.5*m*x.^2.*omegaX^2;
K = (0.5*hbar^2/m).*(kx.^2);

GV = exp(-0.5*V*dt/hbar);
GK = exp(-K*dt/hbar);

wfc /= sqrt(sum(abs(wfc).^2)*dx);
for t=1:10000
        wfc = wfc.*GV;
        wfc = ifft(fft(wfc).*GK);
        wfc = wfc.*GV;
        wfc /= sqrt(sum(abs(wfc).^2)*dx);
end

function energy(wfc,dx,V,K)
        EV = V.*wfc;
        EK = ifft(K.*fft(wfc));
        E = conj(wfc).*(EV + EK)
        return real(sum(E)*dx)/hbar
end

e_0 = 1e100 #Arbitrarily large unrealistic value
e_1 = energy(wfc,dx,V,K);

while e_0-e_1 > 1e-7
        push!(E,e_1)
        wfc = wfc.*GV;
        wfc = ifft(fft(wfc).*GK);
        wfc = wfc.*GV;
        wfc /= sqrt(sum(abs(wfc).^2)*dx);
        e_0 = e_1;
        e_1 = energy(wfc,dx,V,K)
        println(e_1)
end

$$ k_\textrm{max} = \frac{2\pi}{x_\textrm{max}}d $$ n A