The Haldane gap
In this tutorial we will calculate the Haldane gap (the energy gap in the $S = 1$ Heisenberg model) in two different ways. To follow the tutorial you need the following packages:
using MPSKit, MPSKitModels, TensorKit, Plots, Polynomials
import TensorOperations; TensorOperations.disable_cache(); # hideThe Heisenberg model is defined by the following hamiltonian:
\[H = -J∑_{⟨i,j⟩} (X_i X_j + Y_i Y_j + Z_i Z_j)\]
This hamiltonian has an SU(2) symmetry, which we can enforce by using SU(2)-symmetric tensors:
H = heisenberg_XXX(ComplexF64, SU2Irrep; spin=1, J=-1);Finite size extrapolation
We can start the analysis using finite-size methods. The groundstate of this model can be approximated using finite MPS through the use of DMRG.
The typical way to find excited states is to minimize the energy while adding an error term $λ \left|gs\right> \left< gs\right|$ Here we will instead use the quasiparticle ansatz.
In Steven White's original DMRG paper it was remarked that the $S = 1$ excitations correspond to edge states, and that one should define the Haldane gap as the difference in energy between the $S = 2$ and $S = 1$ states. This can be done as follows.
L = 10
physical_space = Rep[SU₂](1 => 1)
virtual_space = Rep[SU₂](0 => 15, 1 => 15, 2 => 10, 3 => 5)
ψ₀ = FiniteMPS(rand, ComplexF64, L, physical_space, virtual_space)
ψ, envs, delta = find_groundstate(ψ₀, H, DMRG(; verbose=false))
E₀ = real.(expectation_value(ψ, H))
En_1, st_1 = excitations(H, QuasiparticleAnsatz(), ψ, envs; sector=SU₂(1))
En_2, st_2 = excitations(H, QuasiparticleAnsatz(), ψ, envs; sector=SU₂(2))
ΔE_finite = real(En_2[1] - En_1[1])-0.2897302384697755We can go even further and doublecheck the claim that $S = 1$ is an edge excitation, by plotting the energy density.
p_density = plot(; xaxis="position", yaxis="energy density")
excited_1 = convert(FiniteMPS, st_1[1])
excited_2 = convert(FiniteMPS, st_2[1])
plot!(p_density, real.(expectation_value(excited_1, H) .- E₀) ./ L; label="S = 1")
plot!(p_density, real.(expectation_value(excited_2, H) .- E₀) ./ L; label="S = 2")Finally, we can obtain a value for the Haldane gap by extrapolating our results for different system sizes.
Ls = 5:5:50
ΔEs = map(Ls) do L
ψ₀ = FiniteMPS(rand, ComplexF64, L, physical_space, virtual_space)
ψ, envs, delta = find_groundstate(ψ₀, H, DMRG(; verbose=false))
En_1, st_1 = excitations(H, QuasiparticleAnsatz(), ψ, envs; sector=SU₂(1))
En_2, st_2 = excitations(H, QuasiparticleAnsatz(), ψ, envs; sector=SU₂(2))
return real(En_2[1] - En_1[1])
end
f = fit(Ls .^ (-2), ΔEs, 1)
ΔE_extrapolated = f.coeffs[1]-0.2155532516993712p_size_extrapolation = plot(; xaxis="L^(-2)", yaxis="ΔE",
xlims=(0, 0.015))
plot!(p_size_extrapolation, Ls .^ (-2), ΔEs;
seriestype=:scatter, label="numerical")
plot!(p_size_extrapolation, x -> f(x); label="fit")Thermodynamic limit
A much nicer way of obtaining the Haldane gap is by working directly in the thermodynamic limit. As was already hinted at by the edge modes, this model is in a non-trivial SPT phase. Thus, care must be taken when selecting the symmetry sectors. The groundstate has half-integer edge modes, thus the virtual spaces must also all carry half-integer charges.
In contrast with the finite size case, we now should specify a momentum label to the excitations. This way, it is possible to scan the dispersion relation over the entire momentum space.
virtual_space_inf = Rep[SU₂](1 // 2 => 20, 3 // 2 => 20, 5 // 2 => 10, 7 // 2 => 10,
9 // 2 => 5)
ψ₀_inf = InfiniteMPS([physical_space], [virtual_space_inf])
ψ_inf, envs_inf, delta_inf = find_groundstate(ψ₀_inf, H)
kspace = range(0, π, 16)
Es, _ = excitations(H, QuasiparticleAnsatz(), kspace, ψ_inf, envs_inf; sector=SU₂(1))
ΔE, idx = findmin(real.(Es))
println("minimum @k = $(kspace[idx]):\t ΔE = $(ΔE)")[ Info: vumps @iteration 1 galerkin = 0.6648889732866761
[ Info: vumps @iteration 2 galerkin = 0.06359011437557995
[ Info: vumps @iteration 3 galerkin = 0.0164265217276693
[ Info: vumps @iteration 4 galerkin = 0.002373916187780843
[ Info: vumps @iteration 5 galerkin = 0.0006289175273869992
[ Info: vumps @iteration 6 galerkin = 0.00022983228453123588
[ Info: vumps @iteration 7 galerkin = 0.00010467216817002606
[ Info: vumps @iteration 8 galerkin = 5.438736021221384e-5
[ Info: CG: initializing with f = -0.936207316315, ‖∇f‖ = 7.6838e-05
[ Info: CG: iter 1: f = -0.936207316802, ‖∇f‖ = 3.3956e-05, α = 1.77e-01, β = 0.00e+00, nfg = 5
[ Info: CG: iter 2: f = -0.936207316979, ‖∇f‖ = 2.4905e-05, α = 3.25e-01, β = 1.65e-01, nfg = 2
[ Info: CG: iter 3: f = -0.936207317061, ‖∇f‖ = 1.3668e-05, α = 2.55e-01, β = 6.41e-01, nfg = 2
[ Info: CG: iter 4: f = -0.936207317085, ‖∇f‖ = 8.2456e-06, α = 2.51e-01, β = 3.10e-01, nfg = 2
[ Info: CG: iter 5: f = -0.936207317094, ‖∇f‖ = 5.4141e-06, α = 2.88e-01, β = 3.53e-01, nfg = 2
[ Info: CG: iter 6: f = -0.936207317098, ‖∇f‖ = 3.5886e-06, α = 2.58e-01, β = 4.29e-01, nfg = 2
[ Info: CG: iter 7: f = -0.936207317100, ‖∇f‖ = 2.3617e-06, α = 2.28e-01, β = 4.39e-01, nfg = 2
[ Info: CG: iter 8: f = -0.936207317100, ‖∇f‖ = 1.5488e-06, α = 2.07e-01, β = 4.38e-01, nfg = 2
[ Info: CG: iter 9: f = -0.936207317101, ‖∇f‖ = 1.0240e-06, α = 2.29e-01, β = 4.16e-01, nfg = 2
[ Info: CG: iter 10: f = -0.936207317101, ‖∇f‖ = 6.7229e-07, α = 2.44e-01, β = 4.52e-01, nfg = 2
[ Info: CG: iter 11: f = -0.936207317101, ‖∇f‖ = 4.5048e-07, α = 2.40e-01, β = 4.23e-01, nfg = 2
[ Info: CG: iter 12: f = -0.936207317101, ‖∇f‖ = 2.8542e-07, α = 2.10e-01, β = 4.53e-01, nfg = 2
[ Info: CG: iter 13: f = -0.936207317101, ‖∇f‖ = 1.7638e-07, α = 2.10e-01, β = 3.99e-01, nfg = 2
[ Info: CG: iter 14: f = -0.936207317101, ‖∇f‖ = 1.1830e-07, α = 2.33e-01, β = 3.84e-01, nfg = 2
[ Info: CG: iter 15: f = -0.936207317101, ‖∇f‖ = 8.0831e-08, α = 2.09e-01, β = 4.50e-01, nfg = 2
[ Info: CG: iter 16: f = -0.936207317101, ‖∇f‖ = 5.8991e-08, α = 2.06e-01, β = 4.67e-01, nfg = 2
[ Info: CG: iter 17: f = -0.936207317101, ‖∇f‖ = 4.3695e-08, α = 2.09e-01, β = 5.30e-01, nfg = 2
[ Info: CG: iter 18: f = -0.936207317101, ‖∇f‖ = 3.2093e-08, α = 1.97e-01, β = 5.52e-01, nfg = 2
[ Info: CG: iter 19: f = -0.936207317101, ‖∇f‖ = 2.3767e-08, α = 2.23e-01, β = 5.39e-01, nfg = 2
[ Info: CG: iter 20: f = -0.936207317101, ‖∇f‖ = 1.7986e-08, α = 2.39e-01, β = 5.49e-01, nfg = 2
[ Info: CG: iter 21: f = -0.936207317101, ‖∇f‖ = 1.3820e-08, α = 2.37e-01, β = 5.73e-01, nfg = 2
[ Info: CG: iter 22: f = -0.936207317101, ‖∇f‖ = 1.0492e-08, α = 2.20e-01, β = 5.91e-01, nfg = 2
[ Info: CG: iter 23: f = -0.936207317101, ‖∇f‖ = 8.5940e-09, α = 2.39e-01, β = 5.76e-01, nfg = 2
[ Info: CG: iter 24: f = -0.936207317101, ‖∇f‖ = 7.4331e-09, α = 2.67e-01, β = 6.72e-01, nfg = 2
[ Info: CG: iter 25: f = -0.936207317101, ‖∇f‖ = 6.0167e-09, α = 2.45e-01, β = 7.48e-01, nfg = 2
[ Info: CG: iter 26: f = -0.936207317101, ‖∇f‖ = 4.5352e-09, α = 2.12e-01, β = 6.55e-01, nfg = 2
[ Info: CG: iter 27: f = -0.936207317101, ‖∇f‖ = 2.9325e-09, α = 1.91e-01, β = 5.68e-01, nfg = 2
[ Info: CG: iter 28: f = -0.936207317101, ‖∇f‖ = 1.7500e-09, α = 1.76e-01, β = 4.18e-01, nfg = 2
[ Info: CG: iter 29: f = -0.936207317101, ‖∇f‖ = 1.0053e-09, α = 1.85e-01, β = 3.56e-01, nfg = 2
[ Info: CG: iter 30: f = -0.936207317101, ‖∇f‖ = 6.0296e-10, α = 2.02e-01, β = 3.30e-01, nfg = 2
[ Info: CG: iter 31: f = -0.936207317101, ‖∇f‖ = 3.9119e-10, α = 2.11e-01, β = 3.60e-01, nfg = 2
[ Info: CG: iter 32: f = -0.936207317101, ‖∇f‖ = 2.8042e-10, α = 2.19e-01, β = 4.21e-01, nfg = 2
[ Info: CG: iter 33: f = -0.936207317101, ‖∇f‖ = 2.1988e-10, α = 2.21e-01, β = 5.14e-01, nfg = 2
[ Info: CG: iter 34: f = -0.936207317101, ‖∇f‖ = 1.9308e-10, α = 2.56e-01, β = 6.15e-01, nfg = 2
[ Info: CG: iter 35: f = -0.936207317101, ‖∇f‖ = 1.6850e-10, α = 2.54e-01, β = 7.71e-01, nfg = 2
[ Info: CG: iter 36: f = -0.936207317101, ‖∇f‖ = 1.4140e-10, α = 2.25e-01, β = 7.61e-01, nfg = 2
[ Info: CG: iter 37: f = -0.936207317101, ‖∇f‖ = 1.1451e-10, α = 2.29e-01, β = 7.04e-01, nfg = 2
[ Info: CG: iter 38: f = -0.936207317101, ‖∇f‖ = 9.1346e-11, α = 2.33e-01, β = 6.55e-01, nfg = 2
[ Info: CG: iter 39: f = -0.936207317101, ‖∇f‖ = 7.0586e-11, α = 2.35e-01, β = 6.37e-01, nfg = 2
[ Info: CG: iter 40: f = -0.936207317101, ‖∇f‖ = 5.8120e-11, α = 2.60e-01, β = 5.97e-01, nfg = 2
[ Info: CG: iter 41: f = -0.936207317101, ‖∇f‖ = 4.7058e-11, α = 2.59e-01, β = 6.78e-01, nfg = 2
[ Info: CG: iter 42: f = -0.936207317101, ‖∇f‖ = 3.2408e-11, α = 2.24e-01, β = 6.55e-01, nfg = 2
[ Info: CG: iter 43: f = -0.936207317101, ‖∇f‖ = 2.1412e-11, α = 2.07e-01, β = 4.75e-01, nfg = 2
[ Info: CG: iter 44: f = -0.936207317101, ‖∇f‖ = 1.3598e-11, α = 1.86e-01, β = 4.35e-01, nfg = 2
[ Info: CG: iter 45: f = -0.936207317101, ‖∇f‖ = 8.9005e-12, α = 1.76e-01, β = 4.05e-01, nfg = 2
[ Info: CG: iter 46: f = -0.936207317101, ‖∇f‖ = 5.5988e-12, α = 1.85e-01, β = 4.27e-01, nfg = 2
[ Info: CG: iter 47: f = -0.936207317101, ‖∇f‖ = 3.6664e-12, α = 2.10e-01, β = 3.97e-01, nfg = 2
[ Info: CG: iter 48: f = -0.936207317101, ‖∇f‖ = 2.5881e-12, α = 2.24e-01, β = 4.29e-01, nfg = 2
[ Info: CG: iter 49: f = -0.936207317101, ‖∇f‖ = 1.9556e-12, α = 2.36e-01, β = 4.95e-01, nfg = 2
[ Info: CG: iter 50: f = -0.936207317101, ‖∇f‖ = 1.5149e-12, α = 2.25e-01, β = 5.73e-01, nfg = 2
[ Info: CG: iter 51: f = -0.936207317101, ‖∇f‖ = 1.2152e-12, α = 2.33e-01, β = 6.02e-01, nfg = 2
[ Info: CG: converged after 52 iterations: f = -0.936207317101, ‖∇f‖ = 9.9968e-13
[ Info: Found excitations for p = 2.9321531433504737
[ Info: Found excitations for p = 3.141592653589793
[ Info: Found excitations for p = 2.722713633111154
[ Info: Found excitations for p = 2.5132741228718345
[ Info: Found excitations for p = 2.303834612632515
[ Info: Found excitations for p = 0.0
[ Info: Found excitations for p = 2.0943951023931953
[ Info: Found excitations for p = 1.8849555921538759
[ Info: Found excitations for p = 1.6755160819145563
[ Info: Found excitations for p = 1.4660765716752369
[ Info: Found excitations for p = 0.8377580409572781
[ Info: Found excitations for p = 1.0471975511965976
[ Info: Found excitations for p = 1.2566370614359172
[ Info: Found excitations for p = 0.6283185307179586
[ Info: Found excitations for p = 0.41887902047863906
[ Info: Found excitations for p = 0.20943951023931953
minimum @k = 0.0: ΔE = -0.2595611997968351
plot(kspace, real.(Es); xaxis="momentum", yaxis="ΔE", label="S = 1")This page was generated using Literate.jl.