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

## モデル

```math
    H
    =
    \sum_{i}{\tau}_i\cdot{\tau}_{i+1}
```
where
```math
    \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

```julia
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のときよりも計算時間は重めです。
```julia
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)
```
![spin1Heisenberg-spectrum.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3875030/e541babc-5e18-e2c9-972b-2e6cec979198.png)
基底状態とその上にしっかりギャップがあるのが見えます。
これだけではgapless系の有限サイズ効果なのか、はたまた無限系でもgappedになるのかは判別がつきません。


### Haldane gap

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

プロットしてみます。
```julia
# gaps
plot(range, real(gaps), yscale=:log10)
xlabel!("length")
ylabel!("Haldane gap")
title!("Haldane gaps of Spin-1 AF Heisenberg chain")
```
`L=9`までやるとこんな感じです。
![haldaneS1gap.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3875030/baf4f431-9e66-85f7-026b-d4bef3ac5d92.png)
めちゃくちゃgaplessへ近づきそうな見た目をしています。
実際サイト数をさらに増やしても全然有限値に収束するようには見えないと言われています。

Shanks変換というテクニックがあります。
深入りすると地獄を見るので、ここでは指数収束の場合だけを紹介します。
[[1]](#reference) から引っ張ってきました。
指数収束する
```math
    \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}
```
から直ちに漸化式
```math
    A_n-A_{n-1}
    =
    be^{c(n-1)}(e^c-1),
    \qquad
    A_{n+1}-A_n
    =
    be^{cn}(e^c-1)
```
を建てることで
```math
    e^c
    =
    \frac{A_{n+1}-A_n}{A_{n}-A_{n-1}}
```
```math
    \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$ をプロットしてギャップを求めてみます。
```julia
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))
```
![gapS1Shanks.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3875030/dabbdab1-d2e5-de15-aea0-eb588271a4f5.png)
これでも結局gaplessへ行くのか有限値に留まるのかよくわかりませんね。


## Reference

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