MPO for Anyon Hubbard model

[Ashinxinl]

Hello everyone,
I am trying to reproduce the results from a paper on the Anyon Hubbard model.

(ImageSource:DOI: 10.1103/PhysRevB.97.115126)

After defining the operators, I encountered an issue: when summing terms in the Hamiltonian using the following expressions:
os += -J,"adag",j,"a",j+1,"ExpIθn",j
os += -J,"adag",j+1,"a",j,"ExpIθndag",j
my program identifies the Hamiltonian as a non-Hermitian Hamiltonian. However, in theory, the Hamiltonian should be Hermitian.

Due to my limited programming experience, I couldn’t find a solution on the forum. Additionally, I am not sure if my ExpIθn function is implemented correctly.
Any suggestions or guidance would be greatly appreciated! Thank you all for your help.
this is a code demo:

using ITensors, ITensorMPS
let
function ITensors.op(::OpName"ExpIθn", ::SiteType"Qudit", d::Int)
    mat = zeros(ComplexF64, d, d)
    for k in 1:d
        mat[k, k] = exp(im *1*pi * (k - 1))
    end
    return mat
end
function ITensors.op(::OpName"ExpIθndag", ::SiteType"Qudit", d::Int)
  mat = zeros(ComplexF64, d, d)
  for k in 1:d
      mat[k, k] = exp(im *1*pi * (k - 1))
  end
  return mat'
end

#----------------------------parameter---------------------
N=60
L_list = collect(Int64, 1:1:N)
U=10
J=1.0
V=10
alpha=1/3
delta=2*pi/3
#----------------------------parameter dmrg----------------
    maxdim = [10,10,20,20,30,30,40,40,80,80,200,200,400,400,600,600]
    nsweeps=40
    cutoff = [1E-9]
    noise=[1E-6,1E-6,1E-6,1E-7,1E-7,1E-7,1E-7,1E-7,1E-7,1E-8,1E-8,1E-8,1E-9,1E-9,1E-9,1E-9,1E-9,1E-9,1E-10,1E-10,1E-10,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,0.0]
    etol=1E-7
#---------------------------- dmrg----------------
    state = [ n % 3 == 0 ? "1" : "0" for n=1:N]
    sites = siteinds("Boson",N,conserve_qns=true,dim=5)
    psi0 =  random_mps(sites,state)
    @show flux(psi0)
    os = OpSum()
    for j=1:N-1
        os += -J,"adag",j,"a",j+1,"ExpIθn",j
        os += -J,"adag",j+1,"a",j,"ExpIθndag",j
    end
    for j=1:N
        os += U/2,"n",j, "n",j
        os -= U/2,"n",j  
    end
    for j=1:N
        os += V*cos(2*pi*alpha*j+delta),"n",j
    end
#pbc
    os += -J,"adag",N,"a",1,"ExpIθn",N
    os += -J,"ExpIθndag",N,"adag",1,"a",N
    H = MPO(os,sites)
#----------------------------strat calculate---------------
    energy,psi = dmrg(H,psi0;nsweeps,noise,maxdim,cutoff)
end

[VinceNeede]

Hi @Ashinxinl
On the top of my head, it looks like there is an error on how you are writing the hermitian conjugate part, in particular that should be

\text{e}^{-i\theta n_j} b_{j+1}^\dagger b_j

The term \text{e}^{-i\theta n_j} commutes with b_{j+1}^\dagger since they act on different sites, but it does not commute with b_j. So try to correct the line

        os += -J,"adag",j+1,"a",j,"ExpIθndag",j

to

        os += -J,"ExpIθndag",j,"adag",j+1,"a",j

and let us know if the problem persists.

EDIT: I found some time to test it, strangely enough I didn’t get any warning, but by explicitly testing if the MPO was hermitian, it was not. By changing the line as said previously, the mpo is now hermitian so your code should now work

[VinceNeede]

I’m a bit concerned about the absence of the warning in my case. Can you please check what versions do you have of ITensorMPS and of KrylovKit?

[Ashinxinl]

Hi VinceNeede
Apologies for the delayed response. Your suggestion was excellent, and my code is now working correctly. However, I don’t quite understand why the order of operators in the opsum() function matters. As far as I know, such an issue usually only arises in the case of fermions. Could this be due to the Jordan-Wigner transformation?

For reference, my package versions are:

• ITensorMPS: v0.2.6
• KrylovKit: v0.8.3

In my original code, warnings started appearing after 16 sweeps.

[VinceNeede]

@Ashinxinl It isn’t related to Jorgan Wigner transformations at all, it’s just how creation and annihilation operators are defined for bosons, i.e. such that

[a,a^\dagger]=1

So the order matters, cause a^\dagger a \neq a a^\dagger but a a^\dagger=1+a^\dagger a.
The difference in the ordering that you are mentioning is related to operators that acts on different sites, that is because the commutation rule for bosons on a lattice is

[a_i,a^\dagger_j]=\delta_{i,\,j}

i.e. it is one if they act on the same site, otherwise it is zero, that means you can write either a_j a_i^\dagger or a_i^\dagger a_j as long as i\neq j. For fermions it is different cause on different sites they anticommute, therefore you would have f_i^\dagger f_j=-f_j f_i^\dagger for i\neq j, but if they are acting on the same site, order matters both for fermions and for bosons cause the commutator is not zero.

[VinceNeede]

Anyway, I would recommend you to upgrade julia and the other packages to the latest versions.
You can download the last julia release using juliaup if you already have it installed, otherwise you can download it using the instructions on the official Julia page Download Julia .
For the packages it should be easier, once upgraded julia, from the REPL just digit the following commands:

julia>]
pkg> up

[Ashinxinl]

@ VinceNeedeThank you for your patience and detailed response. I have already updated my version.
But When attempting to reproduce the image from Figure 2 in the paper,

(ImageSource:DOI: 10.1103/PhysRevB.97.115126)
I found that the results I obtained for ⟨ψₙ₊₁|nj|ψₙ₊₁⟩ - ⟨ψₙ|nj|ψₙ⟩ did not match those in the figure when setting θ to π or any other angle.
⟨ψₙ₊₁|nj|ψₙ₊₁⟩

⟨ψₙ|nj|ψₙ⟩

Δnj

My conditions are consistent with those in the paper, so I suspect that the issue might still lie in the construction of the MPO.

[VinceNeede]

The graph you’re showing is obtained using Open Boundary Condition (OBC), while you are using PBC, could that be the problem? Otherwise, can you show how you are taking the expectation values?

[Ashinxinl]

Yes, I am using OBC conditions, and I obtained the expectation values using the following code.

[Ashinxinl]

@ VinceNeedeThis is my whole code:

using ITensors, ITensorMPS
using CairoMakie

let

function ITensors.op(::OpName"ExpIθn", ::SiteType"Qudit", d::Int)
    mat = zeros(ComplexF64, d, d)
    for k in 1:d
        mat[k, k] = exp(im *pi * (k - 1))
    end
    return mat
end
function ITensors.op(::OpName"ExpIθndag", ::SiteType"Qudit", d::Int)
  mat = zeros(ComplexF64, d, d)
  for k in 1:d
      mat[k, k] = exp(im *pi * (k - 1))
  end
  return mat'
end

#----------------------------parameter---------------------
N=90
L_list = collect(Int64, 1:1:N)
f=30
U=10
J=1
V=10
alpha=1/3
delta=2*pi/3
#----------------------------parameter dmrg----------------
    maxdim = [10,10,20,20,30,30,40,40,80,80,200,200,400,400,600,600]
    nsweeps=60
    cutoff = [1E-10]
    noise=[1E-6,1E-6,1E-6,1E-7,1E-7,1E-7,1E-7,1E-7,1E-7,1E-8,1E-8,1E-8,1E-9,1E-9,1E-9,1E-9,1E-9,1E-9,1E-10,1E-10,1E-10,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,1E-11,0.0]
    etol=1E-7
#----------------------------parameter dmrg----------------
    state = [ n <= f ? "1" : "0" for n=1:N]
   # state[Int64(N/3+1)] = "1" 
    sites = siteinds("Boson",N,conserve_qns=true,dim=5)
    psi0 =  random_mps(sites,state)
    @show flux(psi0)
    os = OpSum()
    for j=1:N-1
        os += -J,"adag",j,"a",j+1,"ExpIθn",j
        os += -J,"ExpIθndag",j,"adag",j+1,"a",j
    end
    for j=1:N
        os += U/2,"n",j, "n",j
        os += -U/2,"n",j  
    end

    for j=1:N
        os += V*cos(2*pi*alpha*j+delta),"n",j
    end
    #pbc
   # os += -J,"adag",N,"a",1,"ExpIθn",N
   # os += -J,"ExpIθndag",N,"adag",1,"a",N
    H = MPO(os,sites)


    H = MPO(os,sites)
#----------------------------strat calculate---------------
    energy,psi = dmrg(H,psi0;nsweeps,eigsolve_krylovdim=6,noise,maxdim,cutoff)



density= expect(psi, "n")
print(density)

f1 = Figure(backgroundcolor = :white) 
ax = Axis(f1[1,1],#xticks = [1,L/2,L],#yticks = [-0.25,0,0.25],
#title = "L=($L)N=($(L/2))h=($h)beta=($beta)V=($V)U=($U)t=($t)",
title = "alpha=(1/3)delta=(2*pi/3)L=($N)N=($(f))V=($V)U=($U)_θ=pi",
xlabel = "L",
)
#ylims!(ax,0.0, 0.7)
resize_to_layout!(f1)
scatterlines!(L_list,density,label="ni")
axislegend(position = :rb)
display(f1)


  end

[VinceNeede]

@Ashinxinl have you actually checked whether the dmrg converged? i.e. incrementing the number of sweeps or the maxdim does not change the result? Cause this is the result I’ve obtained:

The mode is localized on the last site instead of the first one but it shouldn’t really matter, what is important is that it is localized on the extrema.

Anyway, I took the freedom to rewrite your code in a more readable way. Please take a look at (Style Guide · The Julia Language), in particular use functions when possible, add logging messages, and use variable names that actually represent what the variable means (i.e. try to avoid naming variables that are used for the whole program with a single letter, like N, L or f)

Rewriting of the code

using ITensors, ITensorMPS
using LinearAlgebra
using CairoMakie

let
  function ITensors.op(::OpName"ExpIθn", ::SiteType"Boson", dim::Int64; kwargs...)
    return Diagonal([exp(im * pi * (k - 1)) for k in 1:dim])
  end
  function ITensors.op(::OpName"ExpIθndag", ::SiteType"Boson", dim::Int64; kwargs...)
    return Diagonal([exp(-im * pi * (k - 1)) for k in 1:dim])
  end
  #----------------------------parameter---------------------
  chain_length = 90
  Nbosons = 30
  U = 10.
  J = 1.
  V = 10.
  alpha = 1 // 3
  delta = 2 * pi / 3
  #----------------------------parameter dmrg----------------
  maxdim = [10, 10, 20, 20, 30, 30, 40, 40, 80, 80, 200, 200, 400, 400, 600, 600]
  nsweeps = 120
  cutoff = 1.e-12
  #----------------------------parameter dmrg----------------
  """
  	function constructBoseHubbardHamiltonian(J::T, U::T, V::T, alpha::T, delta::T, chain_length::Int, pbc::Bool)::OpSum where T<:Float64

  	Construct the Bose-Hubbard Hamiltonian from the parameters, returns it as a sum of operators.
  """
  function constructBoseHubbardHamiltonian(J::T, U::T, V::T, alpha::T, delta::T, chain_length::Int, pbc::Bool)::OpSum where T<:Float64
    os = OpSum()
    for j = 1:chain_length-1
      os += -J, "adag", j, "a", j + 1, "ExpIθn", j
      os += -J, "ExpIθndag", j, "adag", j + 1, "a", j
    end
    for j = 1:chain_length
      os += U / 2, "n", j, "n", j
      os += -U / 2, "n", j
    end

    for j = 1:chain_length
      os += V * cos(2 * pi * alpha * j + delta), "n", j
    end
    if pbc
      os += -J, "adag", chain_length, "a", 1, "ExpIθn", chain_length
      os += -J, "ExpIθndag", chain_length, "adag", 1, "a", chain_length
    end
    return os
  end

  """
  	function constructBoseHubbardHamiltonian(J::T, U::T, V::T, alpha::T, delta::T, chain_length::Int, pbc::Bool, sites::Vector{<:Index})::MPO where T<:Float64

  	Construct the Bose-Hubbard Hamiltonian from the parameters, returns it as an MPO acting on the sites.
  """
  function constructBoseHubbardHamiltonian(J::T, U::T, V::T, alpha::T, delta::T, chain_length::Int, pbc::Bool, sites::Vector{<:Index})::MPO where T<:Float64
    os = constructBoseHubbardHamiltonian(J, U, V, alpha, delta, chain_length, pbc)
    return MPO(os, sites)
  end

  """
  	function constructInitialState(chain_length::Int, Nbosons::Int, sites::Vector{<:Index})::Vector{String}

  	Construct an initial random state with fixed number of bosons.
  """
  function constructInitialState(chain_length::Int, Nbosons::Int, sites::Vector{<:Index})
    state = [n <= Nbosons ? "1" : "0" for n = 1:chain_length]
    psi0 = random_mps(sites, state)
    return psi0
  end

  #----------------------------strat calculate---------------
  sites = siteinds("Boson", chain_length, conserve_qns=true, dim=5)

  global densities = Dict{Int,Vector{Float64}}()
  global energies = Dict{Int,Float64}()
  if !(@isdefined ψs)
    global ψs = Dict{Int,MPS}()
  end
  for N in Nbosons-1:Nbosons+1
    tmp = get(ψs, N, nothing)
    if !isnothing(tmp)
      @info "Taking previous result"
      psi0 = tmp
      maxdim = 2 * maxlinkdim(psi0)
      sites = siteinds(psi0)
    else
      psi0 = constructInitialState(chain_length, N, sites)
    end
    @info "Constructing Hamiltonian for $N bosons"
    H = constructBoseHubbardHamiltonian(J, U, V, float(alpha), delta, chain_length, false, sites)
    @info "Starting DMRG for $N bosons"
    energy, psi = dmrg(H, psi0; nsweeps=nsweeps, maxdim=maxdim, cutoff=cutoff)
    ψs[N] = psi
		energies[N] = energy
    densities[N] = expect(psi, "n")
  end

  quasiparticle_density = densities[Nbosons+1] - densities[Nbosons]


  f1 = Figure(backgroundcolor=:white)
  ax = Axis(f1[1, 1],
	  title = "α = $alpha δ = $(rationalize(delta/π))π L = $chain_length N = $(Nbosons) V = $V U = $U ",
    xlabel = "Lattice Site i",
		ylabel = "Quasiparticle Density"
  )
  resize_to_layout!(f1)
  scatterlines!(collect(1:chain_length), quasiparticle_density, label="ni")
  axislegend(position=:rb)
  display(f1)

end

Something more complicated that I’ve added is that the eigenvector obtained is saved in a global variable, and if this variable is already defined it doesn’t get redefined. In this way when working from the REPL and calling include("file_name.jl") more than one time, the previous ground state is used as starting state, this allows to simply increment the number of sweeps and check for the convergence in a more easy way. If you have more questions on the code feel free to ask .

[Ashinxinl]

@VinceNeede
Thank you for your patient and detailed response. I believe I previously encountered the issue of DMRG getting stuck in a local minimum during my calculations. After following your advice and increasing the precision, I was able to achieve the results shown in the figure. Wishing you a happy life and success in your work!