Adding a field to Hamiltonian for TEBD gates

Hey everyone, I just have a question about creating the correct gets when the Hamiltonian has an added field. Specifically,

H = \sum_{i=1}^{N-1}J (S_{i}^{x}S_{i+1}^{x} + S_{i}^{y}S_{i+1}^{y}) + \sum_{i=1}^{N}V(i)S_{i}^{z}

Here is what I tried to do but it does not seem to match my small system ED code.

using ITensors

function Gates(N::Int, dt::Float64, J, λ, q, phi, sindex)
    gates = ITensor[]
    for j=1:N-1
        s1 = sindex[j]
        s2 = sindex[j+1]
        if j != (N-1)
            Wj = 2*λ*cos(2*pi*q*j + phi)
            hj = J[j]*0.5*(op("S+", s1)*op("S-", s2) + op("S-", s1)*op("S+", s2)) + Wj*op("Sz", s1)*op("Id", s2)
            Gj = exp(-1.0im*0.5*dt*hj)
        elseif j==(N-1)
            WjN = 2*λ*cos(2*pi*q*N + phi)
            Wj = 2*λ*cos(2*pi*q*j + phi)
            hj = J[j]*0.5*(op("S+", s1)*op("S-", s2) + op("S-", s1)*op("S+", s2)) + WjN*op("Id", s1)*op("Sz", s2) + Wj*op("Sz", s1)*op("Id", s2)
            Gj = exp(-1.0im*0.5*dt*hj)
        end
        push!(gates, Gj)
    end
    append!(gates, reverse(gates))
end

I really appreciate any feedback on if this the correct approach or if there is something better I can do.

Your approach of just including the field term on the left side of each bond, then also including the right-side field term on the last bond looks very reasonable to me.

If you aren’t getting the results you expect, I’d just suggest checking over all the details very carefull and trying limits like N=2 sites. It’s really good that you are comparing to ED.

When you say the results don’t match ED, are they very different or just a little bit different?

Thank you for the quick response. First, if I do not include the field i.e \lambda =0 then my ED and TEBD code match exactly. For short times they match well but then start to deviate away from one another.

Most likely you just have a bug somewhere, either in your ITensor or ED code (i.e. they aren’t really doing the same calculation). Your code looks reasonable to me, but that does not mean it doesn’t have any bugs. Please use the usual approach to finding bugs, such as testing various limits (e.g. \lambda=0 as you already did, and then other limits) or making your same terms but using them to compute the ground state, etc. to try to find what is causing the discrepancy between your two codes.