Example 5: Degenerate perturbation theory in a harmonic oscillator basis

Rimu can also handle non-lattice systems. This example looks at weakly-interacting bosonic particles in a harmonic oscillator external potential using a basis of (Cartesian product) single-particle eigenstates of the harmonic oscillator potential. Blocks of degenerate non-interacting states are coupled by a contact interaction in first order degenerate perturbation theory. This example shows how to generate these blocks and find the energy and angular momentum eigenstates.

A runnable script for this example is located here. Run it with julia HO-example.jl.

First, load all needed modules.

using Rimu
using DataFrames
using LinearAlgebra

Define the system size for $N=2$ particles in a 2D harmonic oscillator allowing $M=4$ levels in each dimension, including the groundstate.

N = 2
M = 4

4

Use a tuple S to define the range of harmonic oscillator states in a Cartesian basis, in this isotropic case $n_x,n_y=0,1,\ldots,M-1$.

S = (M, M)

(4, 4)

In Rimu the $N$-particle states are still stored as Fock states

P = prod(S)
addr = BoseFS(P, M => N)

BoseFS{2,16}(0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

where the numbering of the modes folds in the two spatial dimensions. Use the utility function fock_to_cart to convert a Fock address to human-readable Cartesian quantum numbers for inspection.

fock_to_cart(addr, S)

2-element StaticArraysCore.SVector{2, Tuple{Int64, Int64}} with indices SOneTo(2):
 (3, 0)
 (3, 0)

The output shows that all $N$ particles are in single-particle state $n_x=M-1, n_y=0$.

The harmonic oscillator Hamiltonian HOCartesianContactInteractions handles contact interactions with first-order perturbation theory, so the matrix representation will block according to the non-interacting energy of the basis states. The first task is to find all blocks of basis states with the same energy. The strength of the interaction is not relevant at this point, just that it is non-zero. Use a dummy groundstate address to build the Hamiltonian

H = HOCartesianContactInteractions(BoseFS(P, 1 => N); S)

HOCartesianContactInteractions(BoseFS{2,16}(2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); S=(4, 4), η=(1.0, 1.0), g=1.0, interaction_only=false)

and then a utility function get_all_blocks to find all blocks. The blocks are found by looping over all possible states with N particles in Cartesian states defined by S. Note that this will only work for total energy up to the maximum accessible by a single particle. The $N$-particle groundstate energy for a 2D harmonic oscillator is $E_0 = N \hbar \omega$ and the maximum single-particle energy is $E = (E_0 + M - 1) \hbar \omega$.

block_df = get_all_blocks(H; max_energy = N + M - 1)

This outputs a list of blocks in H indexed by the noninteracting energy of all states in the block, and a single address that can be used to rebuild the block for further analysis.

addr1 = block_df[7,:addr]
E = block_df[7,:block_E0]

5.0

First, notice that all basis states have the same energy, defined by the block

basis1 = build_basis(H, addr1)
map(b -> Hamiltonians.noninteracting_energy(H, b), basis1)

5-element Vector{Float64}:
 5.0
 5.0
 5.0
 5.0
 5.0

There are two blocks at each energy level (except the groundstate), which are different due to parity conservation, which is the only other symmetry in the Cartesian harmonic oscillator. The basis of this other block is different

addr2 = block_df[4,:addr]
basis2 = build_basis(H, addr2);
basis1 ≠ basis2

true

but its basis elements have the same energy

map(b -> Hamiltonians.noninteracting_energy(H, b), basis2)

5-element Vector{Float64}:
 5.0
 5.0
 5.0
 5.0
 5.0

However, since this system is an isotropic harmonic oscillator, it is possible to build simultaneous eigenstates of the angular momentum operator $L_z$, implemented with AxialAngularMomentumHO

Lz = AxialAngularMomentumHO(S)

AxialAngularMomentumHO((4, 4); z_dim = 3, addr = BoseFS{0,16}(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))

$L_z$ does not conserve parity so both blocks are required. First combine the bases of each block and convert to DVecs

dvs = map(b -> DVec(b => 1.0), vcat(basis1, basis2));

and then compute overlaps for the matrix elements of $L_z$

Lz_mat = [dot(v, Lz, w) for v in dvs, w in dvs]

10×10 Matrix{ComplexF64}:
 0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+1.0im      0.0+0.0im      0.0-1.41421im  0.0+0.0im      0.0+0.0im
 0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+1.41421im  0.0+0.0im      0.0-1.0im      0.0-1.41421im  0.0+0.0im
 0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+1.73205im  0.0+0.0im      0.0+0.0im      0.0-2.0im
 0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+1.41421im  0.0+1.0im      0.0+0.0im
 0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+1.73205im
 0.0-1.0im      0.0-1.41421im  0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im
 0.0+0.0im      0.0+0.0im      0.0-1.73205im  0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im
 0.0+1.41421im  0.0+1.0im      0.0+0.0im      0.0-1.41421im  0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im
 0.0+0.0im      0.0+1.41421im  0.0+0.0im      0.0-1.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im
 0.0+0.0im      0.0+0.0im      0.0+2.0im      0.0+0.0im      0.0-1.73205im  0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im      0.0+0.0im

By diagonalising this matrix the eigenstate have energy E and well-defined angular momentum

Diagonalise this matrix to obtain the eigenstates of $L_z$. The eigenvectors provide the linear combinations of basis states with well-defined angular momentum, within the subspace of energy $E$.

Lz_vals, Lz_vecs = eigen(Lz_mat)

Eigen{ComplexF64, Float64, Matrix{ComplexF64}, Vector{Float64}}
values:
10-element Vector{Float64}:
 -2.9999999999999956
 -2.999999999999991
 -0.9999999999999989
 -0.9999999999999982
 -0.9999999999999973
  1.0
  1.0000000000000027
  1.0000000000000027
  3.0
  3.0000000000000004
vectors:
10×10 Matrix{ComplexF64}:
       0.0+0.0im       -5.55112e-17-0.353553im         0.526707-0.23646im       2.77556e-17-0.204124im             0.0+0.0im          0.526707+0.23646im       2.77556e-17+0.204124im             0.0+0.0im       -5.55112e-17+0.353553im           0.0+0.0im
       0.0+0.0im        5.55112e-17-0.5im             -0.372438+0.167202im              0.0-0.288675im             0.0+0.0im         -0.372438-0.167202im     -5.55112e-17+0.288675im             0.0+0.0im                0.0+0.5im                0.0+0.0im
       0.0+0.612372im           0.0+0.0im          -5.55112e-17+2.77556e-17im   -1.2326e-32+2.77556e-17im          0.0+0.353553im          0.0+0.0im                   0.0+2.77556e-17im          0.0+0.353553im           0.0+0.0im                0.0+0.612372im
       0.0+0.0im                0.0+0.353553im      1.11022e-16+1.66533e-16im           0.0-0.612372im             0.0+0.0im       5.55112e-17+2.77556e-17im           0.0+0.612372im             0.0+0.0im                0.0-0.353553im           0.0+0.0im
       0.0-0.353553im    1.2326e-32+1.11022e-16im   1.11022e-16-2.77556e-17im  -2.46519e-32-2.77556e-17im          0.0+0.612372im          0.0-2.77556e-17im   4.93038e-32-2.77556e-17im          0.0+0.612372im   -1.2326e-32+8.32667e-17im        0.0-0.353553im
       0.0+0.0im           0.353553+0.0im                   0.0+0.0im              0.612372+0.0im                  0.0+0.0im               0.0+0.0im              0.612372+0.0im                  0.0+0.0im           0.353553+0.0im                0.0+0.0im
 -0.353553+0.0im                0.0+0.0im                   0.0+0.0im                   0.0+0.0im            -0.612372+0.0im               0.0+0.0im                   0.0+0.0im             0.612372+0.0im                0.0+0.0im           0.353553+0.0im
       0.0+0.0im               -0.5+7.02973e-17im     -0.167202-0.372438im         0.288675-1.95105e-16im          0.0+0.0im         -0.167202+0.372438im         0.288675+2.29062e-16im          0.0+0.0im               -0.5-7.48398e-17im        0.0+0.0im
       0.0+0.0im          -0.353553+8.84171e-17im       0.23646+0.526707im         0.204124-1.1389e-16im   8.32667e-17+0.0im           0.23646-0.526707im         0.204124+2.4427e-16im   8.32667e-17+0.0im          -0.353553-8.31718e-17im        0.0+0.0im
  0.612372+0.0im                0.0+0.0im                   0.0+0.0im                   0.0+0.0im            -0.353553-0.0im               0.0+0.0im                   0.0+0.0im             0.353553+0.0im                0.0+0.0im          -0.612372-0.0im

Finally, consider the effect of interactions by looking at how states in a single block are perturbed. Only the energy shift due to the interaction is relevant so now rebuild the Hamiltonian without the non-interacting energy

Hint = HOCartesianContactInteractions(addr1; S, interaction_only = true)
ΔE = eigvals(Matrix(Hint, addr1))

5-element Vector{Float64}:
 7.155734338400866e-18
 2.7755575615628883e-17
 0.15915494309189537
 0.15915494309189543
 0.15915494309189548

Two eigenstates in this block are unaffected by the interaction and three have a non-zero energy shift.

The default strength of the interaction is g = 1.0. Other interactions strengths can be obtained by using keyword argument g in HOCartesianContactInteractions or by rescaling ΔE since the interactions are handled with first-order perturbation theory.

Rimu also contains HOCartesianEnergyConservedPerDim which is a similar Hamiltonian but with the stricter condition that the contact interaction only connects states that have the same total energy in each dimension, rather than conserving the overall total energy. Both Hamiltonians can handle anisotropic systems by passing a tuple S whose elements are not all the same. This will alter which states are connected by the interaction, but assumes that the harmonic trapping frequencies in each dimension are commensurate

Finished!

This page was generated using Literate.jl.