TEBD on 2D fermionic system

ITensor Julia Questions
1

Hello everyone,

My code below is meant to perform Trotterized imaginary-time evolution on an initial Neel state, for a 2 x 2 Hubbard square lattice at half filling (open boundary conditions). There is nothing wrong with the construction of the Hamiltonian MPO, etc. as confirmed by the DMRG run at the end that yield the correct energy.

using ITensors

function hubbard_TEBD_gates(NN,sites,t,U,τ)
    # NN: list of nearest-neighbor sites
    # sites: siteinds of MPS
    # t, U: Hubbard parameters
    # τ: imaginary-time step
    
    gates = ITensor[]
    
    for edge in NN
        edge[2] < edge[1] && reverse!(edge)
        
        s1, s2 = sites[edge[1]], sites[edge[2]]
        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)
        Gj = exp(-1.0 * τ/2 * hj)
        push!(gates,Gj)
    end

    for i in 1:length(sites)
        si = sites[i]
        hi = U*op("Nupdn",si)
        Gi = exp(-1.0 * τ/2 * hi)
        push!(gates,Gi)
    end
    
    append!(gates,reverse(gates))
    
    return gates
end

function lattice_hubbard_MPO(NN,sites,nx,ny,t,U)
    N = nx*ny
    ampo = OpSum()
    for edge in NN
        ampo += (-t,"Cdagup",edge[1],"Cup",edge[2])
        ampo += (-t,"Cdagup",edge[2],"Cup",edge[1])
        ampo += (-t,"Cdagdn",edge[1],"Cdn",edge[2])
        ampo += (-t,"Cdagdn",edge[2],"Cdn",edge[1])
    end
    for i in 1:N
        ampo += U,"Nupdn",i
    end
    H = MPO(ampo,sites)
    return H
end

let
    nx = 2 #number of rows in lattice
    ny = 2 #number of columns
    t, U = 1, 6 #parameters in Hubbard model
    τ = 0.1 #imaginary-time evolution step-size
    NN = [[1,2],[1,4],[2,3],[3,4]] #nearest-neighbor sites
    
    sites = siteinds("Electron",nx*ny;conserve_qns=true)
    gates = hubbard_TEBD_gates(NN,sites,t,U,τ)
    
    H = lattice_hubbard_MPO(NN,sites,nx,ny,t,U)
    
    ψTEBD = productMPS(sites,n -> isodd(n) ? "Up" : "Dn")
    
    for i in 1:10
        ψTEBD = apply(gates, ψTEBD)
        normalize!(ψTEBD)
        E = inner(ψTEBD',H,ψTEBD)
        
        println()
        @show i
        @show E
    end
    
    ψDMRG = productMPS(sites,n -> isodd(n) ? "Up" : "Dn")
    dims = [4,8,12,16]
    sweeps = Sweeps(length(dims))
    setmaxdim!(sweeps,dims...)
    EDMRG, ψDMRG = dmrg(H,ψDMRG,sweeps)
end

Instead, the problem is coming from construction of my gates in the function hubbard_TEBD_gates( ). I have not performed swaps to make sites adjacent to each other before applying gates, because my understanding is that the apply( ) function should handle all of this for me in ITensor. I apologize if I am asking a version of a question you have already answered, but I genuinely unsure of the best approach for applying these gates for cases with long-range interactions.

2

Hi Max, to make sure I understand your question, are you finding that when you use your approach you are getting incorrect results? Or are you just in the planning stage and asking if you need to perform any extra steps when making fermionic evolution gates that are not nearest-neighbor?

3

Just wanted to quickly add that if you are doing imaginary time evolution you should use

Gj = exp(-1.0im * τ/2 * hj)

so the operators are properly complex

4

Hi Miles, actually I am really curious about the latter question you proposed. When the interactions are beyond nearest-neighbor (NN), a usual workaround is to use SWAP gate to make the interacting sites being NN (Minimally entangled typical thermal state algorithms - IOPscience). For fermions, SWAP gates will cause additional signs. I would like to ask: (1) if the apply function in Itensor.jl utilizes the SWAP gates method to deal with long-range interaction? (2) should I worry about the fermionic cases?

Thank you!
Fo-Hong