Hello again,
I am working on the one dimensional extended Hubbard model for spinless fermions and a light matter coupling that is introduced via the quantized Peierls phase. I want to implement this for open and periodic boundary conditions. This is the corresponding full Hamiltonian for open boundary conditions:
I am trying to implement this using OpSum(), where I compute the matrix form of Peierls phase on a truncated Hilbert space.
function build_peierls_phase(g::Real, dim_ph::Int)::Matrix{ComplexF64}
# zeros on diagonal
d = zeros(Float64, dim_ph)
# off-diagonal entries: sqrt(1), …, sqrt(dim_ph-1)
e = sqrt.(collect(1:dim_ph-1)) # i.e. from 1 to N_ph-1
# diagonalize a + a^\dagger =: A which is tridiagonal in the number basis
A = SymTridiagonal(d, e)
eigenvals, eigenvecs = eigen(A)
# write A = V D V^\dagger with D diagonal matrix of eigenvals and V matrix of eigenvecs
phases = exp.(1im * g .* eigenvals) # .* element-wise multiplication
# U = V * diag(phases) * V†
U = eigenvecs * Diagonal(phases) * eigenvecs'
return ComplexF64.(U) # ensure complex type
end
function xxz_cavity(
sites::Vector{<:Index},
t::Real=1.0,
U::Real=1.0,
g::Real=1.0,
omega::Real=1.0
)::MPO
# unpack sites with convention that the last one is the bosonic one
f_sites = sites[1:end-1]
b_site = sites[end]
L = length(f_sites) # number of fermionic sites
L ≥ 2 || throw(ArgumentError("Need at least two lattice sites"))
P = build_peierls_phase(g, dim(b_site))
os = OpSum()
b = L + 1 # boson site index
for j in 1:(L-1)
# dressed hopping
os += -t, P, b, "c†", j, "c", j+1
os += -t, P', b, "c†", j+1, "c", j
# interaction term
os += U, "n", j, "n", j+1
end
# add boson energy term
os += omega, "N", b
return MPO(os, sites)
end
which can be used for example like this
function main()
L = 4
N_ph = 5
conserve_qns = false
f_sites = siteinds("Fermion", L; conserve_qns=conserve_qns)
b_site = siteind("Boson"; dim=N_ph+1, conserve_qns=conserve_qns)
sites = vcat(f_sites, [b_site])
H = xxz_cavity(sites)
f_state = [isodd(n) ? "1" : "0" for n=1:L]
b_state = "0"
state = vcat(f_state, [b_state])
psi0 = MPS(sites,state)
@show inner(psi0', H, psi0)
end
main() # yields inner(psi0', H, psi0) = 0.0 + 0.0im
The problem is that I want the number of fermions to be conserved, but not the number of bosons (as can be seen in the Hamiltonian). Setting conserve_qns = true yields “ERROR: LoadError: Fluxes not all equal”, since the quantum number (QN) structure of the fermion and boson site are not meant for this.
I tried defining a custom site “CavityMode“ that is in all accounts the same as the “Boson” site except for its QN structure:
using ITensors
using ITensorMPS
using LinearAlgebra
alias(::SiteType"CavityMode") = SiteType"Boson"()
"""
space(::SiteType"CavityMode";
dim = 2,
conserve_qns = false,
conserve_number = false,
qnname_number = "N_f")
Create the Hilbert space for a site of type "CavityMode".
Optionally specify the conserved symmetries and their quantum number labels.
"""
function ITensors.space(
::SiteType"CavityMode";
dim = 2,
conserve_qns = false,
conserve_number = conserve_qns,
qnname_number = "N_f",
)
if conserve_number
return [QN(qnname_number, 0) => dim] # one large QN block with QN=0
end
return dim
end
ITensors.val(vn::ValName, st::SiteType"CavityMode") = val(vn, alias(st))
function ITensors.state(sn::StateName, st::SiteType"CavityMode", s::Index; kwargs...)
return state(sn, alias(st), s; kwargs...)
end
function ITensors.op(on::OpName, st::SiteType"CavityMode", ds::Int...; kwargs...)
return op(on, alias(st), ds...; kwargs...)
end
function ITensors.op(on::OpName, st::SiteType"CavityMode", s1::Index, s_tail::Index...; kwargs...)
rs = reverse((s1, s_tail...))
ds = dim.(rs)
opmat = op(on, st, ds...; kwargs...)
return itensor(opmat, prime.(rs)..., dag.(rs)...)
end
Now, however, I can’t figure out how to implement the Peierls phase. Everything I try leads to errors (defining it again as a matrix fails when conserve_qns = true, trying to implement it as an operator with an extra argument g does not work with OpSum, defining it as an operator where the argument g is inferred from its name does not work when I try to implement it and manually defining it will probably at the latest run into problems when going to periodic boundary conditions).
I found another topic that deals with this exact problem ( Customized Operator Error: "Fluxes not all equal" - #2 by VinceNeede ), but sadly, the final working code is not provided, and I can’t figure it out myself so I am happy for every support that you can give me!
