Unsuccesful TEBD for Time Independent Fermi-Hubbard Model

+1 vote
edited

I have been trying to recreate the dynamics of a driven Fermi-Hubbard Hamiltonian in Julia:

$$\hat{H}=-t_0 \sum_{j,\sigma}\left( e^{-i\Phi(t)}\hat{c}^{\dag}_{j,\sigma}\hat{c}{j+1,\sigma} + h.c. \right)+\sum_{j} \hat{n}_{j,\uparrow}\hat{n}{j,\downarrow}$$
(don't know why this isn't formatting this is the Hamiltonian given by equation 12 in https://journals.aps.org/pra/abstract/10.1103/PhysRevA.101.053408 )

where @@ \sigma @@ is spin, but my first step is attempting this with @@\Phi(t)@@ set to 0, which results in an undriven, time independent model. I have successfully performed DMRG to obtain the ground state (which I know because I have compared energies to exact calculations in the QuSpin package), but I have been unsuccessful in using TEBD (and RK2 and RK4) to evolve the ground state using the time independent Hamiltonian. I know applying the time evolution operator as an exponential of the TI Hamiltonian to the ground state should preserve the ground state, but instead the energy of the state resulting from the application of the time evolution operator increases throughout the evolution. Though I have not tested dependence rigorously, the energy increase seems to be mostly independent of the time step even though it should be @@\mathcal{O}(\delta t)@@ according to https://www.sciencedirect.com/science/article/pii/S0003491619302532. I have implemented a dynamic time step as well, and tested a couple different maximum bond dimensions in DMRG, as well as decreasing the cutoff for both DMRG and the apply function. Any help is greatly appreciated, thanks in advance! My code is below:

"""
Use second order TEBD to evolve psi(t) to psi (t + dt)
"""
function TEBD(psi, dt, time, params, independent, cutoff)

phi = phi_tl(time + dt / 2, params, independent)

# odd gates (1,2),(3,4),(5,6),...
ogates = ITensor[]
# even gates (2,3),(4,5),(6,7),...
egates = ITensor[]
for j=1:params.nsites
s1 = params.space[j]
#periodic BC
if j == params.nsites
s2 = params.space[1]
else
s2 = params.space[j+1]
end
# we have to define the two site operator so contributions
# of each site will be counted twice
ul = ur = params.U / 2
# exp(i phi(t)) and exp(-i phi(t))
eiphi = exp(1.0im * phi)
eiphiconj = conj(eiphi)
# create operator (only scaled parameters are passed in so t is always 1)
hj = -eiphiconj * op("Cdagup",s1) * op("Cup",s2) +
-eiphiconj * op("Cdagdn",s1) * op("Cdn",s2) +
-eiphi * op("Cdagup",s2) * op("Cup",s1) +
-eiphi * op("Cdagdn",s2) * op("Cdn",s1) +
ul * op("Nupdn", s1) * op("Id",s2) +
ur * op("Id",s1) * op("Nupdn", s2)
# odd gate
if j % 2 == 1
Gj = exp(-1.0im * dt * hj)
push!(ogates, Gj)
# even gate
else
Gj = exp(-1.0im * dt / 2 * hj)
push!(egates, Gj)
end
end

gates = ITensor[]
append!(gates, egates)
append!(gates, ogates)
append!(gates, egates)

return apply(gates, psi; cutoff=cutoff)

end

"""
Propogates ground from 0 to tf using method to propogate by an adaptive timestep
Parameters:
ground - ground state (MPS)
tf - final time
method - the function used to propogate a wavefunction over single dt
dti - initial guess for timestep
epsilon - total change allowed in a single time step
independent - boolean that is True if we are evolving with a time independent
Hamiltonian, and False otherwise
"""
function propogation(ground::MPS, params, tf, method, dti, epsilon, independent, cutoff)
# error measures (en is for current step, en1 is for previous step)
en = en1 = 0.0
time = 0.0
# copy states
psi = deepcopy(ground)
# set dt and previous dt
dt = pdt = dti

# initialize vectors for saving times and expectation values
times = [0]
energies = [inner(ground, get_ham(0, params, independent), ground)]
currents = [inner(ground, get_current(0, params, independent), ground)]

# run for the entire time interval
while time < tf

# get the MPS at time + dt
next_psi = method(psi, dt, time, params, independent, cutoff)

# calculate difference between current and next psi
en = difference(psi, next_psi)

# run propogation while the difference is greater than acceptable error
while en > epsilon
dt *= epsilon / en
# get the next MPS and calculate difference
next_psi = method(psi, dt, time, params, independent, cutoff)
en = difference(psi, next_psi)
end

# incremement time and add it to the array
time += dt
times = vcat(times, [time])

# accept wavefunction
psi = deepcopy(next_psi)

# calculate expectations
energies = vcat(energies, [inner(psi, get_ham(time, params, independent), psi)])
currents = vcat(currents, [inner(psi, get_current(time, params, independent), psi)])

en1 = (en1 > 0) ? en1 : en

# adjust for next time step
# https://www.sciencedirect.com/science/article/pii/S0377042705001123
ndt = dt * ((epsilon^2 * dt)/ (en * en1 * pdt))^(1/12)

# update values for next iteration e_{n-1} -> e_n, dt_{n-1} = dt_n,
# dt_n -> dt_{n+1}
en1 = en
pdt = dt
dt = ndt
end
return times, energies, currents
end

commented by (70.1k points)
Thanks for the question. I have two questions back:
1. how quickly does the energy change in time? I don't think there's any guarantee that the TEBD method will preserve the energy of the ground state due to various technical issues about the approximations it makes, but it should preserve it moderately well

2. what is the reason for the factor dt/2 in some of your gates but only dt without the division by 2 in others? I think the correct Trotter formula should always have the /2 here
commented by (130 points)
1. My evolution should run from 0 to about 241, but half the initial energy is lost within .5. My dynamic time step converges on about .005. I compared to the Heisenberg model example (https://itensor.github.io/ITensors.jl/stable/tutorials/MPSTimeEvolution.html) where the only change I made was to start evolution from the ground state as given by DMRG. This model loses about less than a hundredth of a percent of its initial energy over the entire evolution with a timestep of .1.

2. I am using the formula for TEBD2 given by https://www.sciencedirect.com/science/article/pii/S0003491619302532. I have also tried using the same gates used in the example referenced above which made no difference.
commented by (70.1k points)
Are you normalizing your state after each call to apply?
commented by (130 points)
I just tried that, it did not make a difference.
commented by (70.1k points)
Hi, so I would like to help with this, but to do so it would be really useful to have a fully working, minimal code that reproduces the issue. For example, from what I can see, the code above requires a ground state MPS as input, but I don't have the code you used to compute that ground state. Also the code calls a function phi_tl in it that isn't defined above.

Overall, it's hard to know a priori why you are finding that the energy isn't conserved. It could be a very subtle bug in the code somewhere, or perhaps a bug in an ITensor function. I think it's a case where the only way to find out is to run the code on the simplest case that reproduces the issue and then come up with various checks or test cases that draw out the bug.

Could you please provide the minimal working code? Thanks!