スピン1 AF Heisenbergのギャップを厳密対角化で計算in Julia

モデル

    H
    =
    \sum_{i}{\tau}_i\cdot{\tau}_{i+1}

where

    \tau^x
    =
    \frac{1}{\sqrt{2}}
    \left(\begin{matrix}
        & 1 &
        \\
        1 & & 1
        \\
        & 1 &
    \end{matrix}\right),
    \qquad
    \tau^y
    =
    \frac{1}{\sqrt{2}}
    \left(\begin{matrix}
        & -i &
        \\
        i & & -i 
        \\
        & i &
    \end{matrix}\right),
    \qquad
    \tau^z
    =
    \left(\begin{matrix}
        1& & 
        \\
        & 0 & 
        \\
        & & -1
    \end{matrix}\right)

実装

Hamiltonianを組み立てるところまでは以下と同じ。
https://qiita.com/saitotakuma/items/0d12895b9978d78bf0f2

using LinearAlgebra
using Arpack
using Plots

const ⊗ = kron

function pauli()
    Sx = 1/sqrt(2)*[0 1 0; 1 0 1; 0 1 0]
    Sy = 1/sqrt(2)*[0 -im 0; im 0 -im; 0 im 0]
    Sz = [1 0 0; 0 0 0; 0 0 1]
    return Sx, Sy, Sz
end

function heisenbergint(Sx,Sy,Sz)
    return Sx ⊗ Sx + Sy ⊗ Sy + Sz ⊗ Sz 
end

function hamiltonian(L::Int)
    Sx, Sy, Sz = pauli()
    H=zeros(ComplexF64, 3^L,3^L)
    for i in 1:(L-1)
        H += I(3^(i-1)) ⊗ heisenbergint(Sx,Sy,Sz) ⊗ I(3^(L-i-1))
    end
    return H
end

エネルギースペクトル

厳密対角化でスペクトル計算をしてみます。
L=6で20 sかかりました。
S=1/2のときよりも計算時間は重めです。

L=9
@time eigenvalues=eigs(hamiltonian(L), nev=3^L,which=:SR)[1]
@show eigenvalues

# energy spectrum
plt=plot(yshowaxis=false,xlim=(0,2))
xlabel!("energy")
title!("Spin-1 AF Heisenberg chain energy spectrum")

for j in 1:3^L-5
    plot!(fill(abs(eigenvalues[j]),2),[0,.5],label="", lc=:black)
end

display(plt)

基底状態とその上にしっかりギャップがあるのが見えます。
これだけではgapless系の有限サイズ効果なのか、はたまた無限系でもgappedになるのかは判別がつきません。

Haldane gap

S=1 反強磁性 Heisenberg模型は無限系でもギャップ (Haldane gap) を持つことが知られています。
以下、サイト数を増やしてギャップが残るかを確認します。

gaps = []
range=3:8
for L in range
    gap=eigs(hamiltonian(L),nev=2,which=:SR)[1]
    @time gaps = append!(gaps, gap[2]-gap[1])
end
@show gaps

プロットしてみます。

# gaps
plot(range, real(gaps), yscale=:log10)
xlabel!("length")
ylabel!("Haldane gap")
title!("Haldane gaps of Spin-1 AF Heisenberg chain")

L=9までやるとこんな感じです。

めちゃくちゃgaplessへ近づきそうな見た目をしています。
実際サイト数をさらに増やしても全然有限値に収束するようには見えないと言われています。

Shanks変換というテクニックがあります。
深入りすると地獄を見るので、ここでは指数収束の場合だけを紹介します。
[1] から引っ張ってきました。
指数収束する

    \begin{cases}
        A_{n-1}=a+be^{c(n-1)}
        \\
        A_{n}=a+be^{cn}
        \\
        A_{n+1}=a+be^{c(n+1)}
    \end{cases}

から直ちに漸化式

    A_n-A_{n-1}
    =
    be^{c(n-1)}(e^c-1),
    \qquad
    A_{n+1}-A_n
    =
    be^{cn}(e^c-1)

を建てることで

    e^c
    =
    \frac{A_{n+1}-A_n}{A_{n}-A_{n-1}}

    \therefore
    \quad
    a
    =
    A_{n-1}
    -
    \frac{A_n-A_{n-1}}{e^c-1}
    =
    A_{n-1}
    -
    \frac{(A_n-A_{n-1})^2}{A_{n+1}-2A_n+A_{n-1}}
    =
    \frac{A_{n-1}A_{n+1}-A_n^2}{A_{n+1}-2A_n+A_{n-1}}

がわかります。
この $a$ をプロットしてギャップを求めてみます。

shanks = [(gaps[i - 1] * gaps[i + 1] - gaps[i]^2) / (gaps[i + 1] - 2gaps[i] + gaps[i - 1]) for i in 2:length(gaps) - 1]

plot(range[2:end - 1], real(shanks))

これでも結局gaplessへ行くのか有限値に留まるのかよくわかりませんね。

Reference

1. Stack Exchange Mathematics Simplified Explanation of Shanks Transformation (https://math.stackexchange.com/questions/3871781/simplified-explanation-of-shanks-transformation)