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 packages.

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;

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);

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)

Here, 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 an arbitrary N-particle starting address to build the Hamiltonian.

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

Then, use the 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.9999999999999991
 -0.9999999999999982
 -0.9999999999999973
  1.0000000000000002
  1.0000000000000027
  1.0000000000000027
  3.0
  3.0000000000000004
vectors:
10×10 Matrix{ComplexF64}:
       0.0+0.0im       -5.55112e-17-0.353553im         0.534676-0.217841im      2.77556e-17-0.204124im             0.0+0.0im           0.534676+0.217841im      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.378073+0.154037im              0.0-0.288675im             0.0+0.0im          -0.378073-0.154037im     -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+0.0im           -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      5.55112e-17+6.93889e-17im           0.0-0.612372im             0.0+0.0im        5.55112e-17-1.38778e-16im           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   5.55112e-17-2.77556e-17im  -2.46519e-32-2.77556e-17im          0.0+0.612372im  -1.11022e-16-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.154037-0.378073im         0.288675-1.95105e-16im          0.0+0.0im          -0.154037+0.378073im         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.217841+0.534676im         0.204124-1.1389e-16im   8.32667e-17+0.0im           0.217841-0.534676im         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}:
 -1.6810327175633306e-17
  3.469966302164114e-17
  0.15915494309189535
  0.1591549430918954
  0.1591549430918954

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.

This page was generated using Literate.jl.