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

**URL:** <https://itensor.discourse.group/t/using-purification-method-to-prepare-a-finite-temperature-state-of-bose-hubbard-model/971>\
**Category:** ITensor Julia Questions\
**Tags:** mps\
**Created:** [June 24, 2023, 3:52pm UTC](https://itensor.discourse.group/t/using-purification-method-to-prepare-a-finite-temperature-state-of-bose-hubbard-model/971 "2023-06-24T15:52:11Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Wang](https://avatars.discourse-cdn.com/v4/letter/w/edb3f5/32.png) [@Wang](https://itensor.discourse.group/u/Wang)\
**Post date:** [June 24, 2023, 3:52pm UTC](https://itensor.discourse.group/t/using-purification-method-to-prepare-a-finite-temperature-state-of-bose-hubbard-model/971/1 "2023-06-24T15:52:11Z")

</div>

Hello everyone!

I want to use [purification method or ancillas method](http://dx.doi.org/10.1103/PhysRevB.72.220401) 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.

```auto
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.

```auto
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

```

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [June 26, 2023, 8:11pm UTC](https://itensor.discourse.group/t/using-purification-method-to-prepare-a-finite-temperature-state-of-bose-hubbard-model/971/2 "2023-06-26T20:11:32Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Wang](https://avatars.discourse-cdn.com/v4/letter/w/edb3f5/32.png) [@Wang](https://itensor.discourse.group/u/Wang)\
**Post date:** [June 28, 2023, 4:33pm UTC](https://itensor.discourse.group/t/using-purification-method-to-prepare-a-finite-temperature-state-of-bose-hubbard-model/971/3 "2023-06-28T16:33:01Z")

</div>

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).

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

```

```auto
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?