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 QuantumOpticsusing PyPlot

# System Parametersm = 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 Basisxmin = -5xmax = 5Npoints = 100b_position = PositionBasis(xmin, xmax, Npoints)# Hamiltonian in real space basisp = momentum(b_position) # Dense operatorx = position(b_position) # Sparse operatorH = 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);# Hamiltonianp = momentum(b_momentum) # Sparse operatorx = position(b_momentum) # Dense operatorH = 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^2H = LazySum(H_kin, V)

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

# Initial statex0 = 1.5p0 = 0sigma0 = 0.6Ψ0 = gaussianstate(b_position, x0, p0, sigma0);# Time evolutionT = [0:0.1:3;]tout, Ψt = timeevolution.schroedinger(T, Ψ0, H);# Plot dynamics of particle densityx_points = samplepoints(b_position)n = abs.(Ψ0.data).^2V = ω*x_points.^2C = 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