Calculating finite temperature two time correlation function for fixed magnetisation sector

ITensor Julia Questions
1

Hi,

I am working with the following system:

H = H_{XXZ} + h \sum_{j}jS_{j}^{z} + W_{z}\sum_{j}(-1)^{j}S_{j}^{z}S_{j+1}^{z}.

My objective is to compute \text{Tr}(e^{-\beta H}D(t)D(0)), where D = \sum_{j}jS_{j}^{z}. To meet my end, I first create the the finite \beta density matrix via ‘purification’ in the same spirit as demonstrated in ITensors.jl/examples/finite_temperature/purification.jl at main · ITensor/ITensors.jl · GitHub. i.e., first create an infinite temp. density matrix and then ‘cool’ it down to inverse temperature \beta.
Then, I do further an imaginary time evolution of the operator MPO and contract the correct MPOs with that of density matrix and finally take trace. My code gives correct results as I can crosscheck with ED. However, I am interested in a scenario in which I confine my calculation to only Sz=0 sector of H, which I can easily do in ED. However, at this moment I am not sure how to get it done in ITensor in the set up I am working. I shall be grateful if you can help.

Best,
Sourav.

2

The standard way to do that would be to define your site indices with QN conservation enabled, for example:

n = 10 # Number of sites
s = siteinds("S=1/2", n; conserve_qns=true)
3

Hi Matt,

Thanks for your reply. Indeed I tried that as also written in the ITensor tutorial. However, I get some errors which I do not understand. For convenience , let me paste down my code first and then the error

#---------------Define Basic Params------------------------#

N, Jz =8, 1.0
hstark = 1.0
beta_max = 0.2
t_max = 2.
Wz = 0.1
hz = zeros(N)
for i=1:N
    hz[i] = i*hstark
end

sites = siteinds("S=1/2", N; conserve_qns=true)


#------------Initial Density Matrix---------------------------------------#

rho = MPO(sites, "Id") ./ √2 #---------Initial product state with infinite temp.--------------------#

#------time steps----------------#
τ = beta_max/200

#----------Gates for TEBD evo-----------------------------#
G = ITensor[]
for j=1:N-1
    s1 = sites[j]
    s2 = sites[j+1]
    #print(s2)
    hj = Jz* op("Sz",s1) * op("Sz",s2) + 1.0* op("Sy",s1) * op("Sy",s2)+ 1.0* op("Sx",s1) * op("Sx",s2) + hz[j] * op("Sz",s1)*op("Id", s2) + Wz*(-1)^(j)*op("Sz",s1) * op("Sz",s2)
    Gj = exp(-τ*0.5*hj)
    push!(G,Gj)
end
s1 = sites[N-1]
s2 = sites[N]
hL = hz[N] * op("Id",s1)*op("Sz", s2)
GL = exp(-τ*0.5* hL)
push!(G,GL)
append!(G, reverse(G))

#---------------Propagation---------------------------#
for β in 0:τ:beta_max-τ
    rho = apply(G, rho; cutoff=1E-10)
    rho = rho/tr(rho)
end

The error I am getting is the following : (I am just pasting it).

In `setindex!`, the element (1, 2) of ITensor: 
Dim 1: (dim=2|id=382|"S=1/2,Site,n=1")' <Out>
 1: QN("Sz",1) => 1
 2: QN("Sz",-1) => 1
Dim 2: (dim=2|id=382|"S=1/2,Site,n=1") <In>
 1: QN("Sz",1) => 1
 2: QN("Sz",-1) => 1
NDTensors.BlockSparse{ComplexF64, Vector{ComplexF64}, 2}
 2×2
Block(2, 1)
 [2:2, 1:1]
 0.0 + 0.5im
 you are trying to set is in a block with flux QN("Sz",2), which is different from the flux QN("Sz",-2) of the other blocks of the ITensor. You may be trying to create an ITensor that does not have a well defined quantum number flux.

Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] _setindex!!(::ITensors.HasQNs, ::NDTensors.BlockSparseTensor{ComplexF64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{ComplexF64, Vector{ComplexF64}, 2}}, ::ComplexF64, ::Int64, ::Int64)
    @ ITensors ~/.julia/packages/ITensors/LXBUp/src/qn/qnitensor.jl:7
  [3] _setindex!!
    @ ~/.julia/packages/ITensors/LXBUp/src/itensor.jl:1176 [inlined]
  [4] setindex!
    @ ~/.julia/packages/ITensors/LXBUp/src/itensor.jl:1218 [inlined]
  [5] setindex!(T::ITensor, x::ComplexF64, I::CartesianIndex{2})
    @ ITensors ~/.julia/packages/ITensors/LXBUp/src/itensor.jl:1222
  [6] ITensor(::NDTensors.AllowAlias, ::Type{ComplexF64}, A::Matrix{ComplexF64}, inds::Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}; tol::Int64)
    @ ITensors ~/.julia/packages/ITensors/LXBUp/src/qn/qnitensor.jl:301
  [7] ITensor(::NDTensors.AllowAlias, ::Type{ComplexF64}, A::Matrix{ComplexF64}, inds::Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}})
    @ ITensors ~/.julia/packages/ITensors/LXBUp/src/qn/qnitensor.jl:290
  [8] ITensor(::NDTensors.AllowAlias, ::Matrix{ComplexF64}, ::Index{Vector{Pair{QN, Int64}}}, ::Vararg{Index{Vector{Pair{QN, Int64}}}, N} where N; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors ~/.julia/packages/ITensors/LXBUp/src/itensor.jl:402
  [9] ITensor(::NDTensors.AllowAlias, ::Matrix{ComplexF64}, ::Index{Vector{Pair{QN, Int64}}}, ::Vararg{Index{Vector{Pair{QN, Int64}}}, N} where N)
    @ ITensors ~/.julia/packages/ITensors/LXBUp/src/itensor.jl:402
 [10] itensor(::Matrix{ComplexF64}, ::Vararg{Any, N} where N; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors ~/.julia/packages/ITensors/LXBUp/src/itensor.jl:134
 [11] itensor(::Matrix{ComplexF64}, ::Index{Vector{Pair{QN, Int64}}}, ::Vararg{Index{Vector{Pair{QN, Int64}}}, N} where N)
    @ ITensors ~/.julia/packages/ITensors/LXBUp/src/itensor.jl:134
 [12] op(name::String, s::Index{Vector{Pair{QN, Int64}}}; adjoint::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors ~/.julia/packages/ITensors/LXBUp/src/physics/sitetype.jl:324
 [13] op(name::String, s::Index{Vector{Pair{QN, Int64}}})
    @ ITensors ~/.julia/packages/ITensors/LXBUp/src/physics/sitetype.jl:231
 [14] top-level scope
    @ ./In[42]:29
 [15] eval
    @ ./boot.jl:360 [inlined]
 [16] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1116

I am not sure what’s going wrong.

Best,
S

4

Hi Sourav, try re-writing this line:

in terms of S+ and S- instead of Sx & Sy.
JR

5

Hi Jan,

Thanks. It partially solves the issue; in the sense that I do not get any error. But problem is that it does not give the correct answer. I checked with ED that it still gives the answer when one considers all spin sectors in the Hamiltonian. To put it another way, changing

sites = siteinds(“S=1/2”, N) -----> sites = siteinds(“S=1/2”, N; conserve_qns=true)

yields the same result which is puzzling. Could you please point out what is possibly going wrong?

Best,
S

6

Hi @sourav, to understand about how to conserve total Sz within the purification method (which is a physics algorithm question, not totally specific to ITensor), the most you can do is conserve the total of Sz across both the physical and ancilla sites both. Is that what you are also doing within ED? If you’d like to fully constrain the total Sz within the physical sites only (e.g. microcanonical ensemble), this is theoretically possible but the best way to do it is kind of an open research question and there aren’t any efficient ways to do it that have been published in the literature.

Assuming that you are ok with the first kind of conservation (total Sz across both physical and ancilla), then the way to do it involves setting up your initial infinite temperature state to be made of total Sz=0 Bell pairs, such as these:

\ket{\uparrow}_P \ket{\downarrow}_A + \ket{\downarrow}_P \ket{\uparrow}_A

properly normalized. As you can see, an up spin on the physical site is paired with a down spin on the ancilla site, keeping the total to be zero.

Then after setting up the infinite-T state in this way, the time evolution ought to preserve the total Sz (assuming the Hamiltonian actually conserves it), if you set up things in the way that Matt and Jan mentioned above.

7

Hi Miles,

Thanks for reply.

To answer your first question: to calculate the finite temp correlator in ED I do not resort to ‘purification’ anyway. I just rewrite the expression in terms of eigenvector and eigenstates and proceed.

For the second part: I think the idea you provided should work in principle. However, I need a further (naive) clarification as follows:
I want to translate the whole thing in the set up described in ITensors.jl/examples/finite_temperature/purification.jl at main · ITensor/ITensors.jl · GitHub that initialises infinite temp density matrix MPO as:

rho = MPO(sites, "Id") ./ √2

My question is:: for the state you suggested what should be the this initial MPO? Problem is outer product of |\psi\rangle= \frac{1}{\sqrt{2}} (|\uparrow\rangle_{A } |\downarrow\rangle_{P } + |\downarrow\rangle_{A} |\uparrow\rangle_{P }) will generate a 4\times4 matrix which in the set up of the aforementioned example code gives obvious error of DimensionMismatch. Would be nice if it is clarified what I am getting wrong.

Best,
S

8

To answer the first part of your question, no the MPO(sites,"Id") approach would not be implementing the approach I mentioned above. That code takes a simpler approach of just making the initial state to be an identity operator, so that on each site the density matrix would be like (\ket{\uparrow}_P \ket{\uparrow}_A + \ket{\downarrow}_P \ket{\downarrow}_A). As you can see that is like an identity operator because ancilla and physical states match each other.

So if you use an MPO, you would need to make a different MPO where on each site the operator is instead like a matrix that flips the spin. As a state it would be a product of “Bell pairs” or singlets that have a total quantum number of zero for each pair.

I’m not sure what DimensionMismatch error you are referring to. If you are making an MPO acting on a S=1/2 system, you would make 2x2 matrices and set each MPO tensor to be these 2x2 matrices. If you are making an MPS you would do a different but closely related thing.

9

I see. Thanks a lot :slight_smile:

1 Like