Bose Hubbard model with cavity interaction

Hello everyone I trying to simulate bose hubbard model with cavity interaction for 2 label atom. I want to see the phase diagram using DMRG, for that my Hamiltonian is,

Screenshot from 2025-09-16 11-38-34

where,
Screenshot from 2025-09-16 11-39-59



using ITensors, ITensorMPS, LinearAlgebra

ITensors.space(::SiteType"polaritons") = 10

n = [0 0 0 0 0;
     0 1 0 0 0;
     0 0 2 0 0;
     0 0 0 3 0;
     0 0 0 0 4]

a = [0 1 0 0 0;
     0 0 sqrt(2) 0 0;
     0 0 0 sqrt(3) 0;
     0 0 0 0 sqrt(4);
     0 0 0 0 0]

a_t = [0 0 0 0 0;
      1 0 0 0 0;
      0 sqrt(2) 0 0 0;
      0 0 sqrt(3) 0 0;
      0 0 0 sqrt(4) 0]

I_atom = Diagonal(ones(2))
I_boson = Diagonal(ones(5))

Sz = [0.5 0;
     0 -0.5]
sigma_x = [0 1;
          1 0]
sigma_y = [0 -1.0*im;
          1.0*im 0]
sigma_up = sigma_x + 1.0*im*sigma_y
sigma_dn = sigma_x - 1.0*im*sigma_y

ITensors.op(::OpName"A",::SiteType"polaritons") =
kron(a, I_atom)

ITensors.op(::OpName"Adag",::SiteType"polaritons") =
kron(a_t, I_atom)

ITensors.op(::OpName"Num",::SiteType"polaritons") =
kron(n, I_atom)

ITensors.op(::OpName"Sz",::SiteType"polaritons") =
kron(I_boson, Sz)

ITensors.op(::OpName"sig_up",::SiteType"polaritons") =
kron(I_boson, sigma_up)

ITensors.op(::OpName"sig_dn",::SiteType"polaritons") =
kron(I_boson, sigma_dn)

ITensors.op(::OpName"Id",::SiteType"polaritons") =
kron(I_boson, I_atom)

N = 10
s = Index(N, "polaritons")
a = op("A", s)
adag = op("Adag", s)
N = op("Num", s)
Sz = op("Sz", s)
sig_up = op("sig_up", s)
sig_dn = op("sig_dn", s)
Id = op("Id", s)

let
        N  = 10
        t = 0.1
        mu = 0.2
        omega = 0.3
        beta = 0.1
        epsilon = 0.2

        s = Index(N, "polaritons")
        ampo = AutoMPO()

        Sup = AutoMPO()
        Sdn = AutoMPO()
        S_z = AutoMPO()
        for i in 1:N
                Sup += (1.0, "sig_up", i)
        end
        for i in 1:N
                Sdn += (1.0, "sig_dn", i)
        end
        for i in 1:N
                Sdn += (1.0, "Sz", i)
        end

for i in 1:N-1
                ampo += -t, "adag", i, "a", i+1
                ampo += -t, "a", i, "adag", i
        end
        for i in 1:N
                ampo += -mu, "N", i
        end
        for i in 1:N
                ampo += beta, "S_up", i, "a", i
                ampo += beta, "S_dn", i, "adag", i
        end
        for i in 1:N
                ampo += omega, "adag", i, "a", i
                ampo += epsilon, "S_z", i
                ampo += epsilon, 0.5*N, "Id", i
        end

        H = MPO(ampo, s)
        return
end
ERROR: LoadError: MethodError: no method matching MPO(::Sum{Scaled{ComplexF64, Prod{Op}}}, ::Index{Int32})
The type `MPO` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  MPO(::Any, ::Any, ::Any)
   @ ITensors ~/.julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/mpo.jl:13
  MPO(::Vector{<:Index}, ::Any)
   @ ITensors ~/.julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/mpo.jl:97
  MPO(::ITensors.LazyApply.Applied{typeof(sum), Tuple{Array{ITensors.LazyApply.Applied{typeof(*), Tuple{C, Prod{Op}}, @NamedTuple{}}, 1}}, @NamedTuple{}} where C, ::Vector{<:Index}; splitblocks, kwargs...)
   @ ITensors ~/.julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/opsum_to_mpo/opsum_to_mpo_generic.jl:289
  ...

Stacktrace:
 [1] top-level scope



Please help

Two problems I can point out:
1.

This should be OpSum() everywhere. AutoMPO() is deprecated code.

  1. The error is probably coming from here

You can change this to ampo += epsilon*0.5*N, "Id", i

I hope this helps!

Still facing the same issue . The MPO operator showing the the error.

ERROR: LoadError: MethodError: no method matching MPO(::Sum{Scaled{ComplexF64, Prod{Op}}}, ::Index{Int32})
The type `MPO` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  MPO(::Any, ::Any, ::Any)
   @ ITensors ~/.julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/mpo.jl:13
  MPO(::Vector{<:Index}, ::Any)
   @ ITensors ~/.julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/mpo.jl:97
  MPO(::ITensors.LazyApply.Applied{typeof(sum), Tuple{Array{ITensors.LazyApply.Applied{typeof(*), Tuple{C, Prod{Op}}, @NamedTuple{}}, 1}}, @NamedTuple{}} where C, ::Vector{<:Index}; splitblocks, kwargs...)
   @ ITensors ~/.julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/opsum_to_mpo/opsum_to_mpo_generic.jl:289
  ...

Stacktrace:
 [1] top-level scope

I see you are constructing the MPO like this: H = MPO(ampo, s), however you define s as: s = Index(N, "polaritons") which is single Index with a dimension N. Maybe you want s = [Index(10, “polaritons”) for _ in 1:N] to make a chain of N polaritons (each of dimension 10).

Thank you for replying but It still showing error.

ERROR: LoadError: ArgumentError: Overload of "op" or "op!" functions not found for operator name "N" and Index tags: ("polaritons",).
Stacktrace:
  [1] op(name::String, s::Index{Int64}; adjoint::Bool, kwargs::@Kwargs{})
    @ ITensors.SiteTypes ~/.julia/packages/ITensors/AbUYR/src/lib/SiteTypes/src/sitetype.jl:402
  [2] op
    @ ~/.julia/packages/ITensors/AbUYR/src/lib/SiteTypes/src/sitetype.jl:236 [inlined]
  [3] op(s::Index{Int64}, opname::String)
    @ ITensors.SiteTypes ~/.julia/packages/ITensors/AbUYR/src/lib/SiteTypes/src/sitetype.jl:449
  [4] computeSiteProd(sites::Vector{Index{Int64}}, ops::Prod{Op})
    @ ITensors.ITensorMPS ~/.julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/opsum_to_mpo/opsum_to_mpo_generic.jl:136
  [5] svdMPO(ValType::Type{Float64}, os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index{Int64}}; mindim::Int64, maxdim::Int64, cutoff::Float64)
    @ ITensors.ITensorMPS ~/.julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/opsum_to_mpo/opsum_to_mpo.jl:118
  [6] svdMPO
    @ ~/.julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/opsum_to_mpo/opsum_to_mpo.jl:5 [inlined]
  [7] #svdMPO#541
    @ ~/.julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/opsum_to_mpo/opsum_to_mpo.jl:147 [inlined]
  [8] svdMPO
    @ ~/.julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/opsum_to_mpo/opsum_to_mpo.jl:144 [inlined]
  [9] MPO(os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index{Int64}}; splitblocks::Bool, kwargs::@Kwargs{})
    @ ITensors.ITensorMPS ~/.julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/opsum_to_mpo/opsum_to_mpo_generic.jl:299
 [10] MPO(os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index{Int64}})
    @ ITensors.ITensorMPS ~/.julia/packages/ITensors/AbUYR/src/lib/ITensorMPS/src/opsum_to_mpo/opsum_to_mpo_generic.jl:289
 [11] top-level scope

I am still confused about how to put the Hamiltonian in MPO form what I posted earlier. It would be really helpful for me if anyone from this discussion group could help.

I see that you’re using an operator "N" when you are constructing your OpSum but you don’t define an operator "N" for your "polaritons" SiteType, maybe you’re missing that definition?

Right now my code is working but the local density plot not showing good results. It shows that the system reaches superfluid state much earlier. I want to reproduce the PRL 99, 186401 (2007) results. Here I have shared the code please check where I am making mistakes.

using ITensors, ITensorMPS, LinearAlgebra

ITensors.space(::SiteType"polaritons") = 10

n = [0 0 0 0 0;
     0 1 0 0 0;
     0 0 2 0 0;
     0 0 0 3 0;
     0 0 0 0 4]

a = [0 1 0 0 0;
     0 0 sqrt(2) 0 0;
     0 0 0 sqrt(3) 0;
     0 0 0 0 sqrt(4);
     0 0 0 0 0]

a_t = [0 0 0 0 0;
      1 0 0 0 0;
      0 sqrt(2) 0 0 0;
      0 0 sqrt(3) 0 0;
      0 0 0 sqrt(4) 0]

I_atom = Diagonal(ones(2))
I_boson = Diagonal(ones(5))
sz = [0.5 0;
     0 -0.5]
sigma_x = 0.5*[0 1;
          1 0]
sigma_y = 0.5*[0 -1.0*im;
          1.0*im 0]
sigma_up = sigma_x + 1.0*im*sigma_y
sigma_dn = sigma_x - 1.0*im*sigma_y

ITensors.op(::OpName"A",::SiteType"polaritons") =
kron(I_atom, a)

ITensors.op(::OpName"Adag",::SiteType"polaritons") =
kron(I_atom, a_t)

ITensors.op(::OpName"Num",::SiteType"polaritons") =
kron(I_atom, n) + kron(sigma_up*sigma_dn, I_boson)

ITensors.op(::OpName"Sz",::SiteType"polaritons") =
kron(sz, I_boson)

ITensors.op(::OpName"sig_up",::SiteType"polaritons") =
kron(sigma_up, I_boson)
ITensors.op(::OpName"sig_dn",::SiteType"polaritons") =
kron(sigma_dn, I_boson)

ITensors.op(::OpName"Id",::SiteType"polaritons") =
kron(I_atom, I_boson)

s = siteind("polaritons")
A = op("A", s)
Adag = op("Adag", s)
Num = op("Num", s)
Sz = op("Sz", s)
sig_up = op("sig_up", s)
sig_dn = op("sig_dn", s)
Id = op("Id", s)
let
        L  = 5
        #t = 0.24
        mu = 0.5
        omega = 0.5
        beta = 0.1
        epsilon = 0.5
        Energy = Float64[]
        values = Float64[]

        s = siteinds("polaritons", L)

        for t in 0.0:0.012:0.24
                println("t = ",t)

                state = fill(2, L)
                state[L] = 1
                psi0 = productMPS(s,state)

                ampo = OpSum()
                for i in 1:L-1
                        ampo += -t, "Adag", i, "A", i+1
                        ampo += -t, "A", i, "Adag", i+1
                end
                for i in 1:L
                        ampo += -mu, "Num", i
                end

                for i in 1:L
                        ampo += beta, "sig_up", i, "A", i
                        ampo += beta, "sig_dn", i, "Adag", i
                end
                for i in 1:L
                        ampo += epsilon, "Sz", i
                        ampo += epsilon*0.5, "Id", i
                end
                for i in 1:L
                        ampo += omega, "Adag", i, "A", i
                end

                H = MPO(ampo, s)
                nsweeps = 100
                maxdim = [10,20,40,60,100,200]
                cutoff = [10E-10]
                energy,psi = dmrg(H, psi0; nsweeps, maxdim, cutoff)
                 println(exp_val)
                push!(values, exp_val[2])
        end
        println(values)
     
        return
end


Please help me as early as possible.

Thank you