Hi,
Following my previous post (Reproducing results of the purification code snippet for Heisenberg chain), I think I am now able to fully understand the purification code for the Heisenberg chain. I give some details below of further steps I’ve taken with this model. Now the issue I have is when trying to adapt it to tackle my system of interest, Fermi-Hubbard chains.
1. Further steps with Heisenberg model
I have adapted it to track both the energy density e(\beta)=E(\beta)/N and the (von Neumann) entropy density s(\beta)=S(\beta)/N = -\mathrm{Tr}(\rho(\beta) \log_2 \rho(\beta)) with β. I’ve done that by keeping track of the normalisation factors for \rho(\beta) (note that compared to the code snippet provided in ITensor documentation, I have changed the \sqrt{2} factor in the definition of the initial state in favour of 2 (the local dimension) as I think \sqrt{2} was erroneous.)
# Initial state is infinite-temperature mixed state
rho = MPO(s, "Id") ./ 2
energies = []
entropies = []
normalization_factor = tr(rho)
# Do the time evolution by applying the gates
# for Nsteps steps
for β in 0:δτ:beta_max
energy = inner(rho, H)
@printf("β = %.2f energy = %.8f\n", β, energy)
Z_β = (2^N)*normalization_factor
entropy = β*energy + log(Z_β)
append!(energies, energy/N)
append!(entropies, entropy/N/log(2))
rho = apply(gates, rho; cutoff)
normalization_factor *= tr(rho)
rho = rho / tr(rho)
end
When I plot the points (s, e), as expected I can observe s=1 for β=0 and as β goes to infinity (I stopped at β=10), the entropy density saturates at 1/3 which reflects the two-fold degeneracy of the ground state (I am considering a chain of length N=3 in this example and consider log in base 2).
2. Adaption to Fermi-Hubbard
Being back on my feet with the Heisenberg model, I have tried adapting the code to tackle my system of interest, Fermi-Hubbard chains, namely H=U\sum_{i=1}^N n_{i\uparrow}n_{i\downarrow} - t\sum_{i=1}^{N-1} \sum_{\sigma} (c^{\dagger}_{i \sigma}c_{i+1 \sigma} + \mathrm{h.c.}). As of now, this does not work, because as I increase β I observe increasing entropies.
Here is my code:
using ITensors
using ITensorMPS
using Plots
# I define the operators intervening in the trotterized imaginary time evolution
function ITensors.op(::OpName"expτcc", ::SiteType"Electron", s1::Index, s2::Index; τ, t)
h = -t*op( "Cdagup", s1) * op( "Cup", s2) +
-t*op( "Cdagup", s2) * op( "Cup", s1) +
-t*op( "Cdagdn", s1) * op( "Cdn", s2) +
-t*op( "Cdagdn", s2) * op( "Cdn", s1)
return exp(τ * h)
end
function ITensors.op(::OpName"expτnn", ::SiteType"Electron", s::Index; τ, U)
h = U*op("Nupdn", s)
return exp(τ * h)
end
function main(; N=3, cutoff=1E-8, δτ=0.1, beta_max=2.0, U=2, t=1)
# Make an array of 'site' indices
s = siteinds("Electron", N; conserve_qns=false)
# Hopping gates
gates = ITensors.ops([("expτcc", (n, n + 1), (τ=-δτ / 2, t=t)) for n in 1:(N - 1)], s)
# Onsite interaction gates
Ugates = ITensors.ops([("expτnn", (n), (τ=-δτ / 2, U=U)) for n in 1:(N - 1)], s)
append!(gates, Ugates)
# Include gates in reverse order to complete Trotter formula
append!(gates, reverse(gates))
# Initial state is infinite-temperature mixed state
rho = MPO(s, "Id") ./ 4 # local dimension is four (empty, sp up, sp dn, doubly-occ)
hubbard_MPO = AutoMPO()
for i=1:N-1
hubbard_MPO .+= -t, "Cdagup", i, "Cup", i+1
hubbard_MPO .+= -t, "Cdagup", i+1, "Cup", i
hubbard_MPO .+= -t, "Cdagdn", i, "Cdn", i+1
hubbard_MPO .+= -t, "Cdagdn", i+1, "Cdn", i
end
for i=1:N
hubbard_MPO .+= U, "Nupdn", i
end
H = MPO(hubbard_MPO, s)
energies = []
entropies = []
normalization_factor = tr(rho)
# Do the time evolution by applying the gates for Nsteps steps
for β in 0:δτ:beta_max
energy = inner(rho, H)
Z_β = (4^N)*normalization_factor
entropy = β*energy + log(Z_β)
append!(energies, energy/N)
append!(entropies, entropy/N/log(4))
rho = apply(gates, rho; cutoff)
normalization_factor *= tr(rho)
rho = rho / tr(rho)
end
return energies, entropies
end
energies, entropies = main()
plot(entropies, energies, seriestype=:scatter, xlabel = "entropy density", ylabel = "energy density")
Calling this function and plotting the results I get the following scatter plot:
Is there an obvious issue in my code? I can imagine problems could arise from:
- the anti commutation relations, but using the Electron type for the sites and “Cdn” , etc operators, in my understanding, this should automatically be taken care of
- the fact that depending on the initial state we may be restricted to a certain particle number sector (which I don’t want). Here I have considered the initial state to be the richest mixture in that we consider the equal-weight superposition of all basis states, so I shouldn’t be restricted to a given particle number sector.
Best regards,
SpSn