Time Evolution of Fermions in Quantum Spin Hall state w/ Coulomb Repulsion

ITensor Julia Questions
1

Hi,

I am interested in comparing the time evolution of two colliding fermions in the Quantum Spin Hall state (QSH) to that of a trivial model. QSH edge states are a pair of helical 1D edge modes, one spin-up fermion moving to the right and one spin-down fermion moving to the left. They are protected from back-scattering and have linear dispersion. I’m wondering if this would lead to any qualitatively different behavior, or to reduced scattering via the Klein paradox (though I don’t think a \frac{1}{r} Coulomb potential is a steep enough change for Klein tunneling to occur).

I’d like to explore these two scenarios in a one-dimensional system, as I vary the relative strengths of Coulomb repulsion, nearest neighbor hopping, and the scattering rate out of the states. The trivial model is easy to set up and I’ve found some examples in the forums. As for the QSH state, I derived a hopping matrix T that satisfies spin-up states only moving right and spin-down states only moving left:

T \mid \psi_\uparrow \rangle = t \mid \psi_\uparrow \rangle \\ T \mid \psi_\downarrow \rangle =0 \\ T^\dagger \mid \psi_\uparrow \rangle =0 \\ T^\dagger \mid \psi_\downarrow \rangle = t \mid \psi_\downarrow \rangle, \\

for some scalar hopping energy t:

T = \frac{t}{\sqrt{2}}\begin{bmatrix} 1 & -i \\ i & 1 \end{bmatrix}.

Time evolution looks good for this, and is qualitatively different vs. the trivial system, however I have to set ishermitian=false for it to work. I’m assuming this is because I’ve only included the edge mode in the Hamiltonian, rather than the full two-dimensional QSH phase.

I am using ITensorTDVP.jl to run this, from what I’ve read on this post TDVP is preferred to Trotter gates when the Hamiltonian includes long-range interactions, and I want to include the effects of Coulomb repulsion at a small length scale across several next nearest neighbors. My full code is pasted at bottom.

I have a few questions, hopefully requiring short answers, to sanity check I’m going down the right path here. I see that @miles has worked on topological condensed matter in the past (https://arxiv.org/pdf/1104.5493.pdf), so I figured I’d ask here first, but please let me know if anything is better suited for Physics StackExchange or other forum.

  1. Is this an appropriate model to study fermions in QSH edge states? I’m wondering if instead I need to move to the full 2D Hamiltonian, which would add further complexity and there are fewer examples for it. If so, are there any good examples for TDVP in a 2D lattice?
  2. I am most interested in the interplay between hopping and Coulomb repulsion. Is TDVP a suitable technique when 1/r interactions are the key interest? If I’m interested in really small length scales, I wonder if I have to move to a completely different approach like path integral Monte Carlo.
  3. I would also like to add scattering out of the 1D model, which involve a single destruction operator, but I hit the following error when I include that term:
    LoadError: Parity-odd fermionic terms not yet supported by AutoMPO
    I have seen this error in some other posts, is there a work-around to include parity-odd terms any other way?

Full Code is below:

using ITensors: ITensors, MPO, OpSum, inner, randomMPS, siteinds, MPS, expect
using ITensorTDVP: ITensorTDVP, tdvp
using Observers: observer

function ITensors.state(::ITensors.StateName"+", ::ITensors.SiteType"Electron")
  return [0, 1 / sqrt(2), im / sqrt(2), 0]
end
function ITensors.state(::ITensors.StateName"-", ::ITensors.SiteType"Electron")
  return [0, 1 / sqrt(2), -im / sqrt(2), 0]
end

#Number of sites
n_sites = 64
#left state will start at x_init, right state will start at (N-(x_init-1))
x_init = 24
cutoff = 1E-8
#time step
t_step = 0.02
#stop at this time
t_total = 6.0
#nearest neighbor hopping energy
t_hop = -3.0
#scattering energy
t_scatt = 1.0
#on-site coulomb potential ("Ntot")
U = 10
use_topological = false

s = siteinds("Electron", n_sites)

os = OpSum()
#on-site coulomb interaction
if U != 0.0
  for j in 1:n_sites
    for i in j:n_sites
      global os += U / (1 + abs(j - i)), "Ntot", j, "Ntot", i
    end
  end
end
#nearest-neighbor hopping
if use_topological
  for j in 1:(n_sites - 1)
    #attempt #1 with "Up" and "Dn" states
    # global os += (-im * t_hop / 2), "Cdagup", j + 1, "Cup", j
    # global os += (-im * t_hop / 2), "Cdagdn", j, "Cdn", j + 1

    #attempt #2 with "+" and "-" states
    global os += (t_hop / (2 * sqrt(2))), "Cdagup", j, "Cup", j + 1
    global os += (-im * t_hop / (2 * sqrt(2))), "Cdagup", j, "Cdn", j + 1
    global os += (im * t_hop / (2 * sqrt(2))), "Cdagdn", j, "Cup", j + 1
    global os += (t_hop / (2 * sqrt(2))), "Cdagdn", j, "Cdn", j + 1

    global os += (t_hop / (2 * sqrt(2))), "Cdagup", j + 1, "Cup", j
    global os += (im * t_hop / (2 * sqrt(2))), "Cdagup", j + 1, "Cdn", j
    global os += (-im * t_hop / (2 * sqrt(2))), "Cdagdn", j + 1, "Cup", j
    global os += (t_hop / (2 * sqrt(2))), "Cdagdn", j + 1, "Cdn", j
  end
else
  for j in 1:(n_sites - 1)
    global os += t_hop / 2, "Cdagup", j, "Cup", j + 1
    global os += t_hop / 2, "Cdagdn", j, "Cdn", j + 1
    global os += t_hop / 2, "Cdagup", j + 1, "Cup", j
    global os += t_hop / 2, "Cdagdn", j + 1, "Cdn", j
  end
end
#scattering into bulk states
if t_scatt != 0.0
  #TODO AutoMPO doesn't support odd-parity operators
  for j in 1:(n_sites - 1)
    global os += -t_scatt, "Cup", j
    global os += -t_scatt, "Cdn", j
  end
end

H = MPO(os, s)
ψ = MPS(ComplexF64, s, n -> begin
  if n == x_init
    "-" #"Up"
  elseif n == n_sites - (x_init - 1)
    "+" #"Dn"
  else
    "Emp"
  end
end)

function measure_nupdn(; psi, bond, half_sweep)
  if bond == 1 && half_sweep == 2
    return expect(psi, "Nupdn")
  end
  return nothing
end

function measure_ntot(; psi, bond, half_sweep)
  if bond == 1 && half_sweep == 2
    return expect(psi, "Ntot")
  end
  return nothing
end

obs = observer("nupdn" => measure_nupdn, "ntot" => measure_ntot)

ϕ = tdvp(
  H,
  -im * t_total,
  ψ;
  #TODO topological Hamiltonian is not hermitian?
  ishermitian=false,
  time_step=-im * t_step,
  normalize=true,
  maxdim=16,
  cutoff=cutoff,
  outputlevel=1,
  (observer!)=obs,
)

#reduce(vcat,transpose.(x)) converts Vector{Vector{Float64}} to Matrix{Float64}
nupdn = transpose(reduce(vcat, transpose.(obs.nupdn)))
ntot = transpose(reduce(vcat, transpose.(obs.ntot)))

using Plots
heatmap(ntot)

Sorry for the really long post. We use ITensors for machine learning-based things at my work, and my research background is in condensed matter physics, but I’ve never used techniques like TDVP before, so any support is much appreciated! Thank you for building and maintaining such a cool ecosystem.

-Brian