Unexpected oscillations in long-time TEBD results

ITensor team,

I’m trying to time-evolve a Hubbard Hamiltonian, but, no matter the timescale, timestep, or maximum bond dimension, the expectation values of any given observable oscillate up and down without converging to a constant value, which disagrees with intuition and many papers on time evolution using TEBD. Here’s a minimal example for SvN, but the same thing happens with spin correlations, etc.

using ITensors, ITensorMPS

let
  N = 200
  Npart = 100
  maxdim = 150
  tau = 0.025
  ttotal = 1.25

  s = siteinds("Electron", N; conserve_qns=true)
  U = 4.0
  t = 1

  gates = ITensor[]
  for j in 1:(N - 1)
    s1 = s[j]
    s2 = s[j+1]
    hj = 
         -t*op("Cdagup", s1) * op("Cup", s2) +
         -t*op("Cdagup", s2) * op("Cup", s1) +
         -t*op("Cdagdn", s1) * op("Cdn", s2) +
         -t*op("Cdagdn", s2) *op("Cdn", s1) +
          U/2*op("Nupdn", s1) * op("Id", s2) +
          U/2*op("Id", s1) * op("Nupdn", s2)
    Gj = exp(-im * tau/2 * hj)
    push!(gates, Gj)
  end

  append!(gates, reverse(gates))

  state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
  psi = MPS(s, state)


  for t in 0.0:tau:ttotal
    vonneuman_entropy = 0.0
    b = 50
    psi = orthogonalize!(psi, b)
    U,S,V = svd(psi[b], (linkinds(psi, b-1)..., siteinds(psi, b)...))

    for n in 1:dim(S, 1)
        addition = S[n, n]^2
        vonneuman_entropy -= addition * log2(addition)
    end
    @show vonneuman_entropy, t, maxlinkdim(psi)
    
    t≈ttotal && break

    psi = apply(gates, psi; maxdim)
    normalize!(psi)
  end
end

This code is so simple that I figured my problem is the result of a broad misunderstanding of TEBD in ITensor, so I figured this might be a good place to seek help.

Hm, sorry to hear that but it’s hard to say just by looking at your code, which looks reasonable.

Here are a few initial things to try:

  1. you could try replacing all of your gates with identity operators, just to test that the system stays in the initial state and that all of your measurements correspond to what you’d expect your initial state to be like. This would just be a way of testing your measurement code.

  2. how are you measuring the other local observables you mentioned? Are you using e.g. the expect function or a different way?

Thanks a lot for the reply! I am using expect for other local observables, for example, I would compute the average energy cost of double occupancy with something like

sum  =  0.0
for c in 1:N
        sum += 1/N*expect(psi, "Nupdn"; sites = c)
end

Your first question actually lead to something strange: the SvN of a lattice with perfect Neel order evolved under identity matrices stayed very close to 0, only going as far as around -10e-15, but was much more erratic at the beginning of the simulation before decaying a bit. I’ve attached a picture if you’re curious.

That said, my intention is certainly not to get the developers of ITensor to debug my code for me. If there are no mistakes in my code indicative of a pressing misunderstanding, I’ll keep looking for the issue with some classic trial and error. Thanks so much for your help!

For your code snippet, a quick comment - I would recommend not using sum as a variable since that is a function in Julia, and you can do it more efficiently with

total = 1/N*sum(expect(psi, "Nupdn"; sites = 1:N))
1 Like

Hi,

Just to point out that just because your SvN and other observables are always oscillatory and never settling down to some value doesn’t by itself mean that something is going wrong in your code or in TEBD. This happens due to physics reasons, because not every system relaxes (“converges”) to some equilibrium state under time evolution, in which case entanglement entropies, expectation values of certain observables etc will generally show oscillatory behaviour for any simulatable times, or show a very slow convergence to some equilibrium state meaning that these oscillations decay very slowly over a very long period of time. It can happen in one parameter regime and not elsewhere, with some particular initial state and not others, with some kind of observables and not others, etc. Underlyingly these things are typically due to existence of confining/bound excitations, scarred states in the spectrum, etc but all that is physics and inappropriate to discuss at length here but I just wished to point this out (because nothing is wrong per se in your code). So in case you have not already done so, I’d suggest you to benchmark with respect to some particular case explicitly demonstrated in some paper where relaxation has been demonstrated for specified initial states and parameter values over an explicitly shown time-scale. If in such a case your code is not reproducing their results, then perhaps we can conclude that something is off in your code.

1 Like

Thanks for the thoughtful response!

This is a direct recreation of a graph from a paper, where the expectation values converge on a constant value rather quickly. Something to point out is that the SvN isn’t supposed to “converge” but to increase with constant velocity, with the same slope as a line of best fit for my values, which increase similarly but oscillate while doing so. I can’t imagine what could possibly be messing my results up, considering how extremely simple the objective of the code is. That said, possible source could be that the paper used iTBED, although I’m not sure of the underlying physics.

Either way, thanks very much for all of your comments, really helpful stuff!

So I just ran your code (for total time=2.5 secs). With the Neel initial state, at U/t=4 (very strong coupling regime), I am not surprised at the oscillatory increase of SvN. When I do U/t=1, it shows a much more linear increase of SvN. Both conform to physical expectations applicable to the two respective regimes. Substantially increasing the allowed maxdim and ttotal, I am pretty sure SvN will attain a saturation (eventually) as expected from theory.

So I am quite sure there’s nothing off with your code : )

Thanks so much, thats a relief. This makes me question if there’s a problem in the calculation of my ‘double occupancy’ observable (equation 3 in https://arxiv.org/pdf/1503.02010) , which differs by only a few lines from my original code. The problem I’m noticing is that the amplitudes of the oscillations towards a constant value are about 3 times larger than they should be. Please excuse me if it’s not appropriate to post this here, and thank you very much for your help.

using ITensors, ITensorMPS

let
  N = 200
  Npart = 100
  maxdim = 150
  tau = 0.025
  ttotal = 1.25

  s = siteinds("Electron", N; conserve_qns=true)
  U = 4.0
  t = 1

  gates = ITensor[]
  for j in 1:(N - 1)
    s1 = s[j]
    s2 = s[j+1]
    hj = 
         -t*op("Cdagup", s1) * op("Cup", s2) +
         -t*op("Cdagup", s2) * op("Cup", s1) +
         -t*op("Cdagdn", s1) * op("Cdn", s2) +
         -t*op("Cdagdn", s2) *op("Cdn", s1) +
          U/2*op("Nupdn", s1) * op("Id", s2) +
          U/2*op("Id", s1) * op("Nupdn", s2)
    Gj = exp(-im * tau/2 * hj)
    push!(gates, Gj)
  end

  append!(gates, reverse(gates))

  state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
  psi = MPS(s, state)


  for t in 0.0:tau:ttotal
    @show 1/N*sum((expect(psi, "Nupdn"; sites = 1:N)))
    
    t≈ttotal && break

    psi = apply(gates, psi; maxdim)
    normalize!(psi)
  end
end