This notebook can be found on github

Particle in harmonic trap potential

A particle in a harmonic trap is described by a Hamiltonian of the form

\[ H = \frac{\hat{p}^2}{2m} + \frac{1}{2}m \omega^2 \hat{x}^2 \]

For numerical simulations we are forced to work in a basis. For particles there are two common choices. We can either work in real space or alternatively in momentum space.

using QuantumOptics
using PyPlot

# System Parameters
m = 1.
ω = 0.5 # Strength of trapping potential

0.5

The position basis is simply defined by the range of positions it spans by defining a minimal position ($x_\min$) and a maximal one ($x_\max$), respectively. The amount of points in position we account for, which ultimately governs the dimension of the basis, is given by Npoints.

# Position Basis
xmin = -5
xmax = 5
Npoints = 100
b_position = PositionBasis(xmin, xmax, Npoints)

# Hamiltonian in real space basis
p = momentum(b_position) # Dense operator
x = position(b_position) # Sparse operator

H = p^2/2m + 1/2*m*ω^2*dense(x^2)

Of course we could also choose to work in momentum space:

From a PositionBasis QuantumOptics.jl can automatically infer the corresponding MomentumBasis by calculating $p_\mathrm{min} = -\pi/dx$ and $p_\mathrm{max} = \pi/dx$ where $dx = (x_\mathrm{max} - x_\mathrm{min})/N$

b_momentum = MomentumBasis(b_position);

# Hamiltonian
p = momentum(b_momentum) # Sparse operator
x = position(b_momentum) # Dense operator

H = dense(p^2)/2m + 1/2*m*ω^2*x^2

However, both choices are not optimal since in real space the position operator is diagonal while the momentum operator is a completely dense matrix and vice versa for the momentum space. Therefore, the calculation will scale with $N^2$ where $N$ is the dimension of the Hilbert space. A commonly used trick is to utilize fast Fourier transformation to convert the state of the system between real and momentum space. This allows us to always use the diagonal form of the operators which all in all speeds up the calculations to $N \log N$.

This idea is implemented by the FFTOperator which performs a fast Fourier transformation on the multiplied state.

# Transforms a state multiplied from the right side from real space
# to momentum space.
T_px = transform(b_momentum, b_position)

To use this operator in a Hamiltonian we additionally need the concept of lazy operators which allow us to delay certain operations to a later point in the simulation. E.g. the LazyProduct allows us to do $A*(B*x)$ instead of $(A*B)*x$ which means for our case that the matrix-matrix product never has to be calculated directly but only two matrix-vector multiplications instead.

T_xp = dagger(T_px)

x = position(b_position)
p = momentum(b_momentum)

H_kin = LazyProduct(T_xp, p^2/2m, T_px)
V = ω*x^2
H = LazySum(H_kin, V)

Finally we can simulate the time evolution according to a Schroedinger equation.

# Initial state
x0 = 1.5
p0 = 0
sigma0 = 0.6
Ψ0 = gaussianstate(b_position, x0, p0, sigma0);

# Time evolution
T = [0:0.1:3;]
tout, Ψt = timeevolution.schroedinger(T, Ψ0, H);

# Plot dynamics of particle density
x_points = samplepoints(b_position)

n = abs.(Ψ0.data).^2
V = ω*x_points.^2
C = maximum(V)/maximum(n)

figure(figsize=(6,3))
xlabel(L"x")
ylabel(L"| \Psi(t) |^2")
plot(x_points, (V.-3)./C, "k--")

for i=1:length(T)
    Ψ = Ψt[i]
    n .= abs.(Ψ.data).^2
    plot(x_points, n, "C0", alpha=0.9*(float(i)/length(T))^8+0.1)
end