0 votes
asked by (250 points)
edited by

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.

Calculating the one-particle correlation function for the midsite, comparing my results to ED

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:

Entropy for time evolution

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)
        return ϕ
    elseif opname == "Cdagup" || opname == "Cdagdn"
        for i in 1:(siteidx-1) # Don't act with string on-site
            ϕ = apply_op(ϕ, "F", sites, i)
        ϕ = apply_op(ϕ, opname, sites, siteidx)
        return ϕ

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

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
    ϕ[siteidx] = new_ϕj
    return ϕ

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ψ)

        # Measure the correlation fcn 

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
# 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)
# 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)

Thanks so much for the help!

commented by (70.1k points)
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)

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.

1 Answer

+1 vote
answered by (250 points)

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

commented by (70.1k points)
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!
Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.