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
# adjust time step
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