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 A*t and Cup for A*t0. 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!