The Ising CFT spectrum

This tutorial is meant to show the finite size CFT spectrum for the quantum Ising model. We do this by first employing an exact diagonalization technique, and then extending the analysis to larger system sizes through the use of MPS techniques.

using MPSKit, MPSKitModels, TensorKit, Plots, KrylovKit
using LinearAlgebra: eigvals, diagm, Hermitian

The hamiltonian is defined on a finite lattice with periodic boundary conditions, which can be implemented as follows:

L = 12
H = periodic_boundary_conditions(transverse_field_ising(), L)

12-site FiniteMPOHamiltonian{MPSKit.JordanMPOTensor{ComplexF64, TensorKit.ComplexSpace, Union{TensorKit.BraidingTensor{ComplexF64, TensorKit.ComplexSpace}, TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 2, Vector{ComplexF64}}}, TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}, TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 2, Vector{ComplexF64}}, TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}}}:
┬ W[12]: 6×1×1×1 JordanMPOTensor(((ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1) ⊗ ⊕(ℂ^2)) ← (⊕(ℂ^2) ⊗ ⊕(ℂ^1)))
┼ W[11]: 6×1×1×6 JordanMPOTensor(((ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1) ⊗ ⊕(ℂ^2)) ← (⊕(ℂ^2) ⊗ (ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1)))
┼ W[10]: 6×1×1×6 JordanMPOTensor(((ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1) ⊗ ⊕(ℂ^2)) ← (⊕(ℂ^2) ⊗ (ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1)))
┼ W[9]:  6×1×1×6 JordanMPOTensor(((ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1) ⊗ ⊕(ℂ^2)) ← (⊕(ℂ^2) ⊗ (ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1)))
┼ W[8]:  6×1×1×6 JordanMPOTensor(((ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1) ⊗ ⊕(ℂ^2)) ← (⊕(ℂ^2) ⊗ (ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1)))
┼ W[7]:  6×1×1×6 JordanMPOTensor(((ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1) ⊗ ⊕(ℂ^2)) ← (⊕(ℂ^2) ⊗ (ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1)))
┼ W[6]:  6×1×1×6 JordanMPOTensor(((ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1) ⊗ ⊕(ℂ^2)) ← (⊕(ℂ^2) ⊗ (ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1)))
┼ W[5]:  6×1×1×6 JordanMPOTensor(((ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1) ⊗ ⊕(ℂ^2)) ← (⊕(ℂ^2) ⊗ (ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1)))
┼ W[4]:  6×1×1×6 JordanMPOTensor(((ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1) ⊗ ⊕(ℂ^2)) ← (⊕(ℂ^2) ⊗ (ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1)))
┼ W[3]:  6×1×1×6 JordanMPOTensor(((ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1) ⊗ ⊕(ℂ^2)) ← (⊕(ℂ^2) ⊗ (ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1)))
┼ W[2]:  6×1×1×6 JordanMPOTensor(((ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1) ⊗ ⊕(ℂ^2)) ← (⊕(ℂ^2) ⊗ (ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1)))
┴ W[1]:  1×1×1×6 JordanMPOTensor((⊕(ℂ^1) ⊗ ⊕(ℂ^2)) ← (⊕(ℂ^2) ⊗ (ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1 ⊕ ℂ^1)))

Exact diagonalisation

In MPSKit, there is support for exact diagonalisation by leveraging the fact that applying the hamiltonian to an untruncated MPS will result in an effective hamiltonian on the center site which implements the action of the entire hamiltonian. Thus, optimizing the middle tensor is equivalent to optimixing a state in the entire Hilbert space, as all other tensors are just unitary matrices that mix the basis.

energies, states = exact_diagonalization(H; num=18, alg=Lanczos(; krylovdim=200));
plot(real.(energies);
     seriestype=:scatter,
     legend=false,
     ylabel="energy",
     xlabel="#eigenvalue")

Extracting momentum

Given a state, it is possible to assign a momentum label through the use of the translation operator. This operator can be defined in MPO language either diagramatically as

or in the code as:

function O_shift(L)
    τ = BraidingTensor{ComplexF64}(ℂ^2, ℂ^2)
    O = TensorMap(τ)
    return periodic_boundary_conditions(InfiniteMPO([O]), L)
end

O_shift (generic function with 1 method)

We can then calculate the momentum of the groundstate as the expectation value of this operator. However, there is a subtlety because of the degeneracies in the energy eigenvalues. The eigensolver will find an orthonormal basis within each energy subspace, but this basis is not necessarily a basis of eigenstates of the translation operator. In order to fix this, we diagonalize the translation operator within each energy subspace. The resulting energy levels have one-to-one correspondence to the operators in CFT, where the momentum is related to their conformal spin as $P_n = \frac{2\pi}{L}S_n$.

function fix_degeneracies(basis)
    L = length(basis[1])
    M = zeros(ComplexF64, length(basis), length(basis))
    T = O_shift(L)
    for j in eachindex(basis), i in eachindex(basis)
        M[i, j] = dot(basis[i], T, basis[j])
    end

    vals = eigvals(M)
    return angle.(vals)
end

momenta = Float64[]
append!(momenta, fix_degeneracies(states[1:1]))
append!(momenta, fix_degeneracies(states[2:2]))
append!(momenta, fix_degeneracies(states[3:3]))
append!(momenta, fix_degeneracies(states[4:5]))
append!(momenta, fix_degeneracies(states[6:9]))
append!(momenta, fix_degeneracies(states[10:11]))
append!(momenta, fix_degeneracies(states[12:12]))
append!(momenta, fix_degeneracies(states[13:16]))
append!(momenta, fix_degeneracies(states[17:18]))

18-element Vector{Float64}:
  3.1857542530442574e-30
  1.4051172205135508e-30
 -6.311400512248983e-31
 -1.110223024625158e-16
 -1.3602461870053588e-18
  1.0408340855860866e-16
 -4.1470733097570545e-18
  3.371727499927069e-18
  1.5959455978986586e-16
  7.055770862123108e-18
  5.5511151231257815e-17
 -2.973502681954427e-31
 -6.331740687315362e-17
  5.359415723744804e-18
  1.301042606982604e-18
 -9.920449878242352e-17
  0.0
 -6.247725059346555e-19

Calculating scaling dimensions from the energy gap

v = 2.0
Δ = real.(energies[1:18] .- energies[1]) ./ (2π * v / L)
p = plot(momenta, real.(Δ);
         seriestype=:scatter, xlabel="momentum", ylabel="energy", legend=false)
vline!(p, [2π / L * i for i in -3:3]; color="gray", linestyle=:dash)
hline!(p, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color="gray", linestyle=:dash)
p

Finite bond dimension

If we limit the maximum bond dimension of the MPS, we get an approximate solution, but we can reach higher system sizes.

L_mps = 20
H_mps = periodic_boundary_conditions(transverse_field_ising(), L_mps)
D = 64
ψ, envs, δ = find_groundstate(FiniteMPS(L_mps, ℂ^2, ℂ^D), H_mps, DMRG());

[ Info: DMRG init:	obj = -1.971321475690e+01	err = 7.4003e-02
[ Info: DMRG   1:	obj = -2.549098964429e+01	err = 6.9440494815e-03	time = 2.67 sec
[ Info: DMRG   2:	obj = -2.549098968636e+01	err = 9.0032288881e-07	time = 0.66 sec
[ Info: DMRG   3:	obj = -2.549098968636e+01	err = 1.3353746489e-07	time = 0.57 sec
[ Info: DMRG   4:	obj = -2.549098968636e+01	err = 1.0102919353e-08	time = 0.35 sec
[ Info: DMRG   5:	obj = -2.549098968636e+01	err = 3.6147044323e-09	time = 0.40 sec
[ Info: DMRG   6:	obj = -2.549098968636e+01	err = 2.3487622576e-09	time = 0.39 sec
[ Info: DMRG   7:	obj = -2.549098968636e+01	err = 1.6417769963e-09	time = 0.47 sec
[ Info: DMRG   8:	obj = -2.549098968636e+01	err = 1.1602119470e-09	time = 0.64 sec
[ Info: DMRG   9:	obj = -2.549098968636e+01	err = 8.1567961851e-10	time = 0.39 sec
[ Info: DMRG  10:	obj = -2.549098968636e+01	err = 5.7049959097e-10	time = 0.39 sec
[ Info: DMRG  11:	obj = -2.549098968636e+01	err = 4.6429052737e-10	time = 0.40 sec
[ Info: DMRG  12:	obj = -2.549098968636e+01	err = 3.8662080912e-10	time = 0.40 sec
[ Info: DMRG  13:	obj = -2.549098968636e+01	err = 3.2276986394e-10	time = 0.40 sec
[ Info: DMRG  14:	obj = -2.549098968636e+01	err = 2.7010926832e-10	time = 0.40 sec
[ Info: DMRG  15:	obj = -2.549098968636e+01	err = 2.2653819206e-10	time = 0.38 sec
[ Info: DMRG  16:	obj = -2.549098968636e+01	err = 1.9037561223e-10	time = 0.40 sec
[ Info: DMRG  17:	obj = -2.549098968636e+01	err = 1.6027409710e-10	time = 0.49 sec
[ Info: DMRG  18:	obj = -2.549098968636e+01	err = 1.3515030264e-10	time = 0.54 sec
[ Info: DMRG  19:	obj = -2.549098968636e+01	err = 1.1412981266e-10	time = 0.39 sec
[ Info: DMRG conv 20:	obj = -2.549098968636e+01	err = 9.6503774976e-11	time = 11.14 sec

Excitations on top of the groundstate can be found through the use of the quasiparticle ansatz. This returns quasiparticle states, which can be converted to regular FiniteMPS objects.

E_ex, qps = excitations(H_mps, QuasiparticleAnsatz(), ψ, envs; num=16)
states_mps = vcat(ψ, map(qp -> convert(FiniteMPS, qp), qps))
E_mps = map(x -> expectation_value(x, H_mps), states_mps)

momenta_mps = Float64[]
append!(momenta_mps, fix_degeneracies(states[1:1]))
append!(momenta_mps, fix_degeneracies(states[2:2]))
append!(momenta_mps, fix_degeneracies(states[3:3]))
append!(momenta_mps, fix_degeneracies(states[4:5]))
append!(momenta_mps, fix_degeneracies(states[6:9]))
append!(momenta_mps, fix_degeneracies(states[10:11]))
append!(momenta_mps, fix_degeneracies(states[12:12]))
append!(momenta_mps, fix_degeneracies(states[13:16]))

v = 2.0
Δ = real.(energies[1:18] .- energies[1]) ./ (2π * v / L)
plot(momenta_mps, Δ;
     seriestype=:scatter, xlabel="momentum", ylabel="energy", legend=false)
vline!(p, [2π / L * i for i in -3:3]; color="gray", linestyle=:dash)
hline!(p, [0, 1 / 8, 1, 9 / 8, 2, 17 / 8]; color="gray", linestyle=:dash)
p

This page was generated using Literate.jl.