# Problem with TEBD on Fermi-Hubbard model with Electron indices

**URL:** <https://itensor.discourse.group/t/problem-with-tebd-on-fermi-hubbard-model-with-electron-indices/2570>\
**Category:** ITensor Julia Questions\
**Created:** [January 28, 2026, 11:16am UTC](https://itensor.discourse.group/t/problem-with-tebd-on-fermi-hubbard-model-with-electron-indices/2570 "2026-01-28T11:16:34Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![MatteDAnna](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/mattedanna/32/861_2.png) [@MatteDAnna](https://itensor.discourse.group/u/MatteDAnna)\
**Post date:** [January 28, 2026, 11:16am UTC](https://itensor.discourse.group/t/problem-with-tebd-on-fermi-hubbard-model-with-electron-indices/2570/1 "2026-01-28T11:16:34Z")

</div>

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.

```auto
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](https://github.com/ITensor/ITensorMPS.jl/issues/159) on [ITensors.jl](https://github.com/ITensor/ITensorMPS.jl)]

---

<div class="post-metadata">

**Author:** ![corbett5](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/corbett5/32/886_2.png) [@corbett5](https://itensor.discourse.group/u/corbett5)\
**Post date:** [January 28, 2026, 6:38pm UTC](https://itensor.discourse.group/t/problem-with-tebd-on-fermi-hubbard-model-with-electron-indices/2570/2 "2026-01-28T18:38:43Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![MatteDAnna](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/mattedanna/32/861_2.png) [@MatteDAnna](https://itensor.discourse.group/u/MatteDAnna)\
**Post date:** [February 16, 2026, 2:10pm UTC](https://itensor.discourse.group/t/problem-with-tebd-on-fermi-hubbard-model-with-electron-indices/2570/3 "2026-02-16T14:10:54Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![corbett5](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/corbett5/32/886_2.png) [@corbett5](https://itensor.discourse.group/u/corbett5)\
**Post date:** [February 16, 2026, 7:41pm UTC](https://itensor.discourse.group/t/problem-with-tebd-on-fermi-hubbard-model-with-electron-indices/2570/4 "2026-02-16T19:41:11Z")

</div>

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.

```julia
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)

```