Using purification method to prepare a finite temperature state of Bose-Hubbard model.

ITensor Julia Questions
1

Hello everyone!

I want to use purification method or ancillas method to get a finite temperature state with an characteristic energy E_{e} and with half-filling bosons in a Bosonic system.

Its Hamiltonian is a simple Bose-Hubbard model.

\hat{H}=\sum_i t(\hat{a}^{\dagger}_{i}\hat{a}_{i+1}+\hat{a}_{i}^{\dagger}\hat{a}_{i+1})+\frac{U}{2}\hat{n}_i(\hat{n}_i-1)

I start imaginary time-evolution from an infinite temperature state which is an Identity operator \rho_o=\frac{1}{d^L}\hat{I} according to the paper.

My question is how to find an specific finite temperature state which is characterized by its energy E_e and particle numbers N=L/2.

Another particularly trivial question is that how can one obtain local observables from density operator \hat{\rho}(\beta).
I have tried tr(rho[i]*dag(op("N",sites[i]))) to get particle numbers on the site i, but it doesn’t work.

tr(rho1[2]*dag(op("N",sites[2])))
ERROR: MethodError: no method matching outer()
Closest candidates are:
  outer(::NDTensors.DenseTensor{ElT1}, ::NDTensors.DenseTensor{ElT2}) where {ElT1, ElT2} at C:\Users\黄旺\.julia\packages\NDTensors\oXCSC\src\dense\tensoralgebra\outer.jl:31
  outer(::NDTensors.DiagTensor{ElT1, N1}, ::NDTensors.DiagTensor{ElT2, N2}) where {ElT1, ElT2, N1, N2} at C:\Users\黄旺\.julia\packages\NDTensors\oXCSC\src\diag\tensoralgebra\outer.jl:25
  outer(::NDTensors.DiagBlockSparseTensor{ElT1, N1}, ::NDTensors.DiagBlockSparseTensor{ElT2, N2}) where {ElT1, ElT2, N1, N2} at C:\Users\黄旺\.julia\packages\NDTensors\oXCSC\src\blocksparse\diagblocksparse.jl:425
  ...
Stacktrace:
 [1] combiner(inds::Vector{Index{Vector{Pair{QN, Int64}}}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ITensors C:\Users\Wang\.julia\packages\ITensors\4aoLl\src\qn\qnitensor.jl:405
 [2] combiner(inds::Vector{Index{Vector{Pair{QN, Int64}}}})
   @ ITensors C:\Users\Wang\.julia\packages\ITensors\4aoLl\src\qn\qnitensor.jl:401
 [3] _tr(T::ITensor; plev::Pair{Int64, Int64}, tags::Pair{TagSet, TagSet})
   @ ITensors C:\Users\Wang\.julia\packages\ITensors\4aoLl\src\tensor_operations\matrix_algebra.jl:4
 [4] _tr
   @ C:\Users\Wang\.julia\packages\ITensors\4aoLl\src\tensor_operations\matrix_algebra.jl:2 [inlined]
 [5] #tr#256
   @ C:\Users\Wang\.julia\packages\ITensors\4aoLl\src\tensor_operations\matrix_algebra.jl:20 [inlined]
 [6] tr(T::ITensor)
   @ ITensors C:\Users\Wang\.julia\packages\ITensors\4aoLl\src\tensor_operations\matrix_algebra.jl:19
 [7] top-level scope
   @ REPL[105]:1

Here is my code.

sites = siteinds("Boson",L;conserve_qns = true,conserve_number = true,dim = 5)
#Halmitonian whcih only acts on physical sites
hal=OpSum()
    for i=1:2:L
            hal += -j1,"Adag",i,"A",i+2
            hal += -j1,"Adag",i+2,"A",i
          hal += U/2,"N",i,"N",i
          hal += -U/2,"N",i
    end
    hal += U/2,"N",L-1,"N",L-1
    hal += -U/2,"N",L-1
 H = toMPO(hal,sites)

#troter gates
gates=ITensor[]
    for i=1:2:L
        s1=sites[i]
        s2=sites[i+2]
        hj = U/2*op("N * N",s1)*op("Id",s2)
        hj += -U/2*op("N",s1)*op("Id",s2)
        hj += -j1*op("Adag",s1)*op("A",s2)
        hj += -j1*op("Adag",s2)*op("A",s1)
        Gj = exp(-tau / 2 * hj)
        push!(gates,Gj)
    end
    hj = U/2*op("N * N",sites[L-1])*op("Id",sites[L])
    hj += -U/2*op("N",sites[L-1])*op("Id",sites[L])
    Gj = exp(-tau / 2 * hj)
    push!(gates,Gj)
    append!(gates,reverse(gates))

#infinite temperature mixed state 
rho = MPO(sites, "Id")/(5^L)


#imaginary time-evolution
for beta in 0:tau:beta_max
        println("beta=",beta)
        energy_rho=real(inner(rho,H))
        println("energy_rho=$energy_rho")
        density_rho=[]
        for i=1:2:L
            density_i=tr(rho[i]*op("N",sites[i]))
            append!(density_rho,density_i)
        end
        println("density_psi=$density_rho")
       rho=apply(gates,rho)
end
2

Hi Wang,
Regarding your first question (which is a physics question, not specific to ITensor), when working with the purification method one is preparing a thermal state e^{-\beta H} (actually it is customary in many texts to make e^{-\beta H/2} then use two copies of this i.e. squaring it, when computing observables). So that is a state whose energy depends on \beta. So if you want the average energy E_e to be a certain value, you have to increase the \beta and continue computing the energy until it reaches the value you want. It’s practically impossible to know the value of \beta you will need in advance. Similarly for particle number, you will need to introduce a chemical potential and adjust it until you get the average particle number you want. In principle there could be a way to formulate the purification method to directly target a fixed particle number sector, but this technique hasn’t been developed very much at this time.

Regarding your other question, what kind of object is rho? Is it an MPO? If so, then trying to trace rho[i] times an operator acting on site i cannot give a scalar, because rho[i] will have virtual MPO indices connecting to the other MPO tensors. Please instead review the definition of the property that you are trying to compute and include the approximation of rho as an MPO. Hope that helps.

3

Hi Miles,
Thanks for your reply, it helps me a lot.
I have tried delta() and scalar() to get local obervables from MPO of density matrix at site i rho[i], it returns exactly the quantities it shoud be (at least for product state).

sites=siteinds("Boson",4;conserve_qns=true,dim=5)
states=["1","2","3","4"]
psi=MPS(sites,states)
rho=outer(psi',psi)

julia> den_2_p=rho[2]*dag(den_2)
ITensor ord=2
(dim=1|id=148|"Link,l=2") <Out>
 1: QN("Number",0) => 1
(dim=1|id=382|"Link,l=1") <In>
 1: QN("Number",0) => 1
NDTensors.BlockSparse{Float64, Vector{Float64}, 2}
julia> den_2_p*delta(dag(inds(rho[2])[3]),dag(inds(rho[2])[4]))
ITensor ord=0
NDTensors.BlockSparse{Float64, Vector{Float64}, 0}

julia> scalar(den_2_p*delta(dag(inds(rho[2])[3]),dag(inds(rho[2])[4])))
2.0

julia> den_3_p=rho[3]*dag(op("N",sites[3]))
ITensor ord=2
(dim=1|id=880|"Link,l=3") <Out>
 1: QN("Number",0) => 1
(dim=1|id=148|"Link,l=2") <In>
 1: QN("Number",0) => 1
NDTensors.BlockSparse{Float64, Vector{Float64}, 2}

julia> scalar(den_3_p*delta(dag(inds(rho[3])[3]),dag(inds(rho[3])[4])))
3.0

Further, I want to ask a how can one obtain correlation functions from MPO of density matrix?