Problem with TEBD on Fermi-Hubbard model with Electron indices

ITensor Julia Questions
1

I am trying to prepare the ground state of a simple 1d Fermi-Hubbard model and then I want to time-evolve it with TEBD to benchmark some other method where the dynamics is given in terms of a Trotterized circuit.

Ideally I would prepare the ground state of FH at a certain repulsion U_0, and then quench it with the evolution given by the Hamiltonian at different repulsion U. However, at the moment I am trying to do it for the simple case of same U_0 where the ground state should remain unchanged by Hamiltonian dynamics (up to small Trotter errors). So I expect that both energy and overlap with the original state remain constant over time.
However, if I use the dagger operators with the Jordan-Wigner strings attached (like Cdagup , Cdn etc), I don’t get this conservation of energy and fidelity wrt. to the initial state.
Changing this to the JW-less operators A... fixes the issue, but this leaves me very suspicious. What am I doing wrong in my implementation?
I am attaching a MWE of my issue, where changing ladder_op_type = "C" to ladder_op_type = "A" reproduces what I described.

From a recent chat with Miles and from what I could read on the forum my understanding is that I should use some OpSum to build my gates, it takes into account the JW string automagically, but I could not find informations on how to build my gates \exp(-i \theta~ \mathrm{op\_sum}) starting from the OpSum object.

using ITensors, ITensorMPS

function hopping_gate(sites, i, j, spin, theta)
    ladder_op_type = "C"
    hj = op("$(ladder_op_type)dag$spin", sites[i])* op("$(ladder_op_type)$spin", sites[j]) + op("$(ladder_op_type)dag$spin", sites[j])* op("$(ladder_op_type)$spin", sites[i])
    return exp(-im * theta * hj)
end

let 

#number of sites, interactions etc 
N_spinful_sites = 8
U = 1.
t = 1.
dτ = 0.04

#choose between conserving quanties or not 
conserve = true
if conserve
    sites = siteinds("Electron", N_spinful_sites; conserve_qns=true)
else
    sites = siteinds("Electron", N_spinful_sites; conserve_qns=false)
end


#Build Hamiltonian and Trotter layer 

ladder_op_type = "C"

os = OpSum()
trotter_layer::Vector{ITensor} = []

#up hopping 
for j in 1:(N_spinful_sites - 1)
    os += -t, "$(ladder_op_type)dagup", j, "$(ladder_op_type)up", j + 1
    os += -t, "$(ladder_op_type)dagup", j + 1, "$(ladder_op_type)up", j
    push!(trotter_layer, hopping_gate(sites, j, j + 1, "up", -dτ*t/2))
end

#down hopping
for j in 1:(N_spinful_sites - 1)
    os += -t, "$(ladder_op_type)dagdn", j, "$(ladder_op_type)dn", j + 1
    os += -t, "$(ladder_op_type)dagdn", j + 1, "$(ladder_op_type)dn", j
    push!(trotter_layer, hopping_gate(sites, j, j + 1, "dn", -dτ*t/2))
end

#onsite repulsion 
for j in 1:(N_spinful_sites)
    os += U, "Nupdn", j 
    push!(trotter_layer, exp(-im * dτ * U * op("Nupdn", sites[j])))
end

#down hopping
for j in 1:(N_spinful_sites - 1)
    push!(trotter_layer, hopping_gate(sites, j, j + 1, "dn", -dτ*t/2))
end 
#up hopping
for j in 1:(N_spinful_sites - 1)
    push!(trotter_layer, hopping_gate(sites, j, j + 1, "up", -dτ*t/2))
end

#build Hamiltonian 
H = MPO(os, sites)

#build initial state for DMRG
if conserve 
    Nup = 3
    Ndown = 3
    state = fill("Emp", N_spinful_sites)
    @assert Nup + Ndown <= N_spinful_sites
    for j=1:Nup
        state[j] = "Up"
    end
    for j=(Nup+1):(Nup+Ndown)
        state[j] = "Dn"
    end

    psi0 = random_mps(sites, state)
else
    psi0 = random_mps(sites)
end
@show sum(expect(psi0, "Nup")), sum(expect(psi0, "Ndn"))

#DMRG
maxdim = [20, 40, 80, 160, 500]
cutoff=1.e-12
nsweeps = 2*length(maxdim)
energy, psi_gs = dmrg(H, psi0; nsweeps, maxdim, cutoff)
@show sum(expect(psi_gs, "Nup")), sum(expect(psi_gs, "Ndn"))


#dynamics 
χ = 2000
nlayers = 10
psi = deepcopy(psi_gs)
@show typeof(psi)
for k=1:nlayers
    psi = ITensorMPS.apply(trotter_layer, psi; maxdim=χ, cutoff=1.e-18)
    energy_val = real(inner(psi, Apply(H, psi)))
    overlap_with_gs = abs(inner(psi_gs, psi))^2
    println("<psi|H|psi> = $energy_val, |<psi_gs|psi>|^2 = $overlap_with_gs, maxlinkdim(psi) = $(maxlinkdim(psi))")
end


end

Thank you for your help,
Matteo

[this is a duplicate of Issue#159 on ITensors.jl]

Problem with exp and sparse matrices
2

The answer is that the Jordan-Wigner string isn’t really handled by the op function. Looking at the matrix for Cdagdn we see that the JW string is handled on-site because of -1 entry, but in reality with the JW mapping Cdagdn is a non-local operator that acts on an extensive number of sites.

julia> matrix(op("Cdagdn", siteind("Electron")))
4×4 Matrix{Float64}:
 0.0   0.0  0.0  0.0
 0.0   0.0  0.0  0.0
 1.0   0.0  0.0  0.0
 0.0  -1.0  0.0  0.0

julia> matrix(op("Cdagup", siteind("Electron")))
4×4 Matrix{Float64}:
 0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0

The JW strings are handled correctly in OpSum so you should always use that for things that require a full MPO (i.e. DMRG and TDVP). My guess is to why using A appears to work is that you are constructing the MPO with A as well.

3

Thank you very much for your reply!
What about the case I want to do the dynamics using TEBD (say, a simple Trotterization) instead of TDVP? How can I create the individual exponentials while taking in a systematic way that always handles the JW strings correctly?

4

The safest way would be to construct the MPO corresponding to the gate using OpSum, then contract the gate into a single tensor and exponentiate it.

sites = siteinds("Electron", 2; conserve_qns=true)

os = OpSum()
os .+= -1, "Cdagup", 1, "Cup", 2
os .+= -1, "Cdagup", 2, "Cup", 1
H = contract(MPO(os, sites))

dt = 0.01
gate = exp(-im * dt * H)
1 Like