Skip to content

Commit

Permalink
interaction strength for HO hamiltonians (#207)
Browse files Browse the repository at this point in the history 
* init

* parse show output tests

* Apply suggestions from code review

Co-authored-by: Joachim Brand <j.brand@massey.ac.nz>

---------

Co-authored-by: Joachim Brand <j.brand@massey.ac.nz>
  • Loading branch information
@christofbradly @joachimbrand
christofbradly and joachimbrand authored Jul 24, 2023
1 parent 2eb3917 commit 9a61ea4
Show file tree
Hide file tree
Showing 3 changed files with 40 additions and 21 deletions.
29 changes: 19 additions & 10 deletions src/Hamiltonians/HOCartesianEnergyConserved.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,21 +100,24 @@ end
"""
HOCartesianEnergyConserved(addr; S, η, g = 1.0, interaction_only = false)
Implements a bosonic harmonic oscillator in Cartesian basis with contact interactions.
Implements a bosonic harmonic oscillator in Cartesian basis with contact interactions
```math
\\hat{H} = \\sum_{i} \\epsilon_i n_i + \\frac{g}{2}\\sum_{ijkl} V_{ijkl} a^†_i a^†_j a_k a_l
\\hat{H} = \\sum_{i} \\epsilon_i n_i + \\frac{g}{2}\\sum_{ijkl} V_{ijkl} a^†_i a^†_j a_k a_l ,
```
Indices ``\\mathbf{i}, \\ldots`` are ``D``-tuples for a ``D``-dimensional harmonic oscillator.
The energy scale is defined by the first dimension i.e. ``\\hbar \\omega_x`` so that
single particle energies are
with the additional restriction that the interactions only couple states with the same
non-interacting energy. See [`HOCartesianEnergyConservedPerDim`](@ref) for a model
that conserves energy separately per spatial dimension.
For a ``D``-dimensional harmonic oscillator indices ``\\mathbf{i}, \\mathbf{j}, \\ldots``
are ``D``-tuples. The energy scale is defined by the first dimension i.e. ``\\hbar \\omega_x``
so that single particle energies are
```math
\\frac{\\epsilon_\\mathbf{i}}{\\hbar \\omega_x} = (i_x + 1/2) + \\eta_y (i_y+1/2) + \\ldots``.
\\frac{\\epsilon_\\mathbf{i}}{\\hbar \\omega_x} = (i_x + 1/2) + \\eta_y (i_y+1/2) + \\ldots.
```
The factors ``\\eta_y, \\ldots`` allow for anisotropic trapping geometries and are assumed to
be greater than `1` so that ``\\omega_x`` is the smallest trapping frequency.
Matrix elements ``V_{\\mathbf{ijkl}}`` are for a contact interaction calculated in this basis using
first-order degenerate perturbation theory.
Matrix elements ``V_{\\mathbf{ijkl}}`` are for a contact interaction calculated in this
basis using first-order degenerate perturbation theory.
```math
V_{\\mathbf{ijkl}} = \\delta_{\\mathbf{i} + \\mathbf{j}}^{\\mathbf{k} + \\mathbf{l}}
\\prod_{d \\in x, y,\\ldots} \\mathcal{I}(i_d,j_d,k_d,l_d),
Expand Down Expand Up @@ -198,7 +201,8 @@ function HOCartesianEnergyConserved(
energies = reshape(map(x -> dot(aspect1, Tuple(x) .- 1/2), states), P)
end

u = sqrt(prod(aspect1)) * g / 2
# the aspect ratio appears from a change of variable when calculating the interaction integrals
u = g / (2 * sqrt(prod(aspect1)))

max_level = S_sort[1] - 1
r = 0:max_level
Expand All @@ -208,9 +212,14 @@ function HOCartesianEnergyConserved(
end

function Base.show(io::IO, h::HOCartesianEnergyConserved)
print(io, "HOCartesianEnergyConserved($(h.addr); S=$(h.S), η=$(h.aspect1), u=$(h.u))")
flag = iszero(h.energies)
# invert the scaling of u parameter
g = h.u * 2 * sqrt(prod(h.aspect1))
print(io, "HOCartesianEnergyConserved($(h.addr); S=$(h.S), η=$(h.aspect1), g=$g, interaction_only=$flag)")
end

Base.:(==)(H::HOCartesianEnergyConserved, G::HOCartesianEnergyConserved) = all(map(p -> getproperty(H, p) == getproperty(G, p), propertynames(H)))

function starting_address(h::HOCartesianEnergyConserved)
return h.addr
end
Expand Down
24 changes: 15 additions & 9 deletions src/Hamiltonians/HOCartesianEnergyConservedPerDim.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ function four_oscillator_integral_1D(i, j, k, l; max_level = typemax(Int))
p = sqrt(gamma(k + 1) * gamma(l + 1)) / sqrt(2 * gamma(i + 1) * gamma(j + 1)) / pi^2
s = sum(gamma(t + 1/2) * gamma(i - t + 1/2) * gamma(j - t + 1/2) / (gamma(t + 1) * gamma(k - t + 1) * gamma(l - t + 1)) for t in 0 : min_kl)

return p * s / 2
return p * s
end

"""
Expand Down Expand Up @@ -79,17 +79,17 @@ end
Implements a bosonic harmonic oscillator in Cartesian basis with contact interactions
```math
\\hat{H} = \\sum_{i} ϵ_i n_i + \\frac{g}{2}\\sum_{ijkl} V_{ijkl} a^†_i a^†_j a_k a_l
\\hat{H} = \\sum_{i} ϵ_i n_i + \\frac{g}{2}\\sum_{ijkl} V_{ijkl} a^†_i a^†_j a_k a_l,
```
with the adiitional restriction that the interactions only couple states with the same
with the additional restriction that the interactions only couple states with the same
energy in each dimension separately. See [`HOCartesianEnergyConserved`](@ref) for a model that
conserves total energy.
Indices ``\\mathbf{i}, \\ldots`` are ``D``-tuples for a ``D``-dimensional harmonic oscillator.
The energy scale is defined by the first dimension i.e. ``\\hbar \\omega_x`` so that
single particle energies are
For a ``D``-dimensional harmonic oscillator indices ``\\mathbf{i}, \\mathbf{j}, \\ldots``
are ``D``-tuples. The energy scale is defined by the first dimension i.e. ``\\hbar \\omega_x``
so that single particle energies are
```math
\\frac{\\epsilon_\\mathbf{i}}{\\hbar \\omega_x} = (i_x + 1/2) + \\eta_y (i_y+1/2) + \\ldots``.
\\frac{\\epsilon_\\mathbf{i}}{\\hbar \\omega_x} = (i_x + 1/2) + \\eta_y (i_y+1/2) + \\ldots.
```
The factors ``\\eta_y, \\ldots`` allow for anisotropic trapping geometries and are assumed to
be greater than `1` so that ``\\omega_x`` is the smallest trapping frequency.
Expand Down Expand Up @@ -162,7 +162,8 @@ function HOCartesianEnergyConservedPerDim(
energies = reshape(map(x -> dot(aspect1, Tuple(x) .- 1/2), states), P)
end

u = sqrt(prod(aspect1)) * g / 2
# the aspect ratio appears from a change of variable when calculating the interaction integrals
u = g / (2 * sqrt(prod(aspect1)))

M = maximum(S)
bigrange = 0:M-1
Expand All @@ -176,9 +177,14 @@ function HOCartesianEnergyConservedPerDim(
end

function Base.show(io::IO, h::HOCartesianEnergyConservedPerDim)
print(io, "HOCartesianEnergyConservedPerDim($(h.addr); S=$(h.S), η=$(h.aspect1), u=$(h.u))")
flag = iszero(h.energies)
# invert the scaling of u parameter
g = h.u * 2 * sqrt(prod(h.aspect1))
print(io, "HOCartesianEnergyConservedPerDim($(h.addr); S=$(h.S), η=$(h.aspect1), g=$g, interaction_only=$flag)")
end

Base.:(==)(H::HOCartesianEnergyConservedPerDim, G::HOCartesianEnergyConservedPerDim) = all(map(p -> getproperty(H, p) == getproperty(G, p), propertynames(H)))

function starting_address(h::HOCartesianEnergyConservedPerDim)
return h.addr
end
Expand Down
8 changes: 6 additions & 2 deletions test/Hamiltonians.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1340,6 +1340,8 @@ end
b2 = Hamiltonians.find_Ebounds(3, 2, S, H.aspect)
@test b1 == b2
@test !(b1 === b2)

@test eval(Meta.parse(repr(H))) == H
end

@testset "HOCartesianEnergyConservedPerDim" begin
Expand All @@ -1364,7 +1366,7 @@ end

# interaction matrix elements
@test count(H.vtable .== 0) == 70
@test sum(H.vtable) ≈ 3.3630246382916664
@test sum(H.vtable) ≈ 2 * 3.3630246382916664

# aspect ratio
S = (4,2,2)
Expand All @@ -1374,7 +1376,9 @@ end
H = HOCartesianEnergyConservedPerDim(addr; S, η = (1,2,3))
@test H.aspect1 == (1.0,2.0,3.0)
H = HOCartesianEnergyConservedPerDim(addr; S, η = 2)
@test H.aspect1 == (1.0,2.0,2.0)
@test H.aspect1 == (1.0,2.0,2.0)

@test eval(Meta.parse(repr(H))) == H
end

@testset "Angular momentum" begin
Expand Down

0 comments on commit 9a61ea4

Please sign in to comment.