Computing correlation function

Question

Hi there! I am using Julia to compute the one-particle spectral function for the Hubbard Holstein system. While I have some degree of confidence in my DMRG calculation (I get the same results as ED) and time evolution (the entropy of the unmodified ground state function, evolved in time, remains very flat), I am getting correlation results that are exactly sign flipped relative to the ED calculation. I've attached a plot of the results for the case where @@N=8, t=1, U=8, \Delta t=0.01, T=1@@. The correlation function I (believe I am) plotting is @@\langle c_j^\dagger (t) c_i(0)\rangle@@ for @@i=j=4@@, the midpoint.

As you can see, there is an exact sign flip between the ED analysis and the results I've obtained by running the algorithm described here on my DMRG ground state. Here is a plot showing the corresponding entropy of both wavefunctions:

Because I think the issue might be localized to how I am computing the correlation function, I tried a few different ways of doing so, but ultimately was unable to fix the problem. Here is a snippet of the code I am using to time-evolve the wavefunction and compute the correlation function:

function apply_onesite_operator(ϕ::MPS, opname::String, sites, siteidx::Int)
    ϕ = copy(ϕ) 

    ## Account for fermion sign using Jordan-Wigner strings ##
    if opname == "Cup" || opname == "Cdn"
        ϕ = apply_op(ϕ, opname, sites, siteidx)
        for i in reverse(1:(siteidx-1)) # Don't act with string on-site
            ϕ = apply_op(ϕ, "F", sites, i)
        end
        return ϕ
    elseif opname == "Cdagup" || opname == "Cdagdn"
        for i in 1:(siteidx-1) # Don't act with string on-site
            ϕ = apply_op(ϕ, "F", sites, i)
        end
        ϕ = apply_op(ϕ, opname, sites, siteidx)
        return ϕ
    end

    # Otherwise, just apply the operator as usual
    return apply_op(ϕ, opname, sites, siteidx)
end

function apply_op(ϕ::MPS, opname::String, sites, siteidx::Int)
    ϕ = copy(ϕ) # Make a copy of the original state

    orthogonalize!(ϕ, siteidx)
    new_ϕj = op(opname,sites[siteidx]) * ϕ[siteidx] # Apply the local operator
    noprime!(new_ϕj) 
    ϕ[siteidx] = new_ϕj
    return ϕ
end

And now for the compute_correlation functions itself:

function compute_correlations(dmrg_results::DMRGResults, A_t0::String, A_t::String, HH::HubbardHolsteinModel, p::Parameters)
    # Results 
    corrs = []

    # The wavefunction being acted upon at t=0, |ψ⟩ = A_t0|ϕ⟩
    ϕ = copy(dmrg_results.ground_state)
    ψ = copy(ϕ)

    # Apply A_t0 to middle site
    ψ = apply_onesite_operator(ψ, A_t0, HH.sites, p.mid)

    nsteps = floor(p.T/p.τ) # Number of time steps for time evolution
    t = 0.0
    for step in 1:nsteps
        ϕ = apply(HH.gates, ϕ; maxdim=p.TEBD_maxdim, cutoff=p.TEBD_cutoff)
        ψ = apply(HH.gates, ψ; maxdim=p.TEBD_maxdim, cutoff=p.TEBD_cutoff) # evolve forward

        t += p.τ 

        function measure_corr(j::Int)
            A_tψ = apply_onesite_operator(ψ, A_t, HH.sites, j)
            return inner(ϕ,A_tψ)
        end

        # Measure the correlation fcn 
        push!(corrs,measure_corr.(collect(1:p.N)))
    end
    hcat(corrs...)
end

Due to how I've set up the code, I pass in Cdagup for At and Cup for At0. Here is the Hamiltonian AMPO, which I use during DMRG:

ampo = OpSum()
for j=1:N-1
    ampo += -t,"Cdagup",j,"Cup",j+1
    ampo += -t,"Cdagup",j+1,"Cup",j
    ampo += -t,"Cdagdn",j,"Cdn",j+1
    ampo += -t,"Cdagdn",j+1,"Cdn",j

    ampo += U,"Nupdn",j,"I",j

    ampo += ω,"Nb",j

    ampo += g0,"Ntot(Bd+B)",j

    ampo += g1,"Ntot",j,"Bdag+B",j+1 
    ampo += g1,"Ntot",j+1,"Bdag+B",j
end
# Edge site
ampo += U,"Nupdn",N
ampo += ω,"Nb",N
ampo += g0,"Ntot(Bd+B)",N
H = MPO(ampo,sites)

And the trotter gates for time evolution:

gates = ITensor[]
for j=1:N-1
    s1 = sites[j] # site j
    s2 = sites[j+1] # site j+1

    hj_twosite = -t*(op("Cdagup*F",s1) * op("Cup",s2)  # t * (c^†_jσ c_{j+1}σ + h.c.)
             -op("Cup*F",s1) * op("Cdagup",s2) 
             +op("Cdagdn*F",s1) * op("Cdn",s2) 
             -op("Cdn*F",s1) * op("Cdagdn",s2)) 
            + g1*(op("Ntot",s1) * op("Bdag+B",s2))
            + g1*(op("Bdag+B",s1) * op("Ntot",s2))

    hj_onesite = U*(op("Nupdn",s1) * op("I",s2))    
                + ω*(op("Nb",s1) * op("I",s2))   
                + g0*(op("Ntot(Bd+B)",s1) * op("I",s2))

    Gj_twosite = exp(-1.0im * τ/2 * hj_twosite)
    Gj_onesite = exp(-1.0im * τ/2 * hj_onesite)
    push!(gates,Gj_twosite)
    push!(gates,Gj_onesite)
end
# End site 
hn = U*op("Nupdn",sites[N]) 
    + ω*op("Nb",sites[N]) 
    + g0*op("Ntot(Bd+B)",sites[N])
Gn = exp(-1.0im * τ/2 * hn)
push!(gates,Gn)
append!(gates,reverse(gates))

Thanks so much for the help!

Comments

To follow up on this, for those reading, a key step which was done correctly here is in the first box of code, adding code like:

for i in reverse(1:(siteidx-1)) # Don't act with string on-site
   ϕ = apply_op(ϕ, "F", sites, i)
end

where note that apply_op is a custom function defined by this user but the key is the idea of putting Jordan-Wigner string to the left of any fermionic operator like "C" or "Cdag" when acting them one at a time onto an MPS, with a time evolution operator in between them.

Answer 1

Wow okay this is so embarrassing but the bug was in my plotting function all along. The code as described above works perfectly fine! Thank you to everyone who tried to help, and I am so sorry for wasting your time like this

Comments

No worries, thanks for posting the question because this kind of this is still challenging to do correctly in ITensors.jl as it is (though in fairness would be in most other tensor libraries too). We are working to add some features to make fermions even easier to use in the future. It's not an obvious thing to get right!