Spin 1/2 反強磁性Heisenbergスペクトルを厳密対角化で計算 in Julia

Juliaの勉強がてら数値計算をしてみます。

## モデル

```math
    H=\sum_j{S}_j\cdot{S}_{j+1}
```
where
```math
    S_j^x
    =
    \left(\begin{matrix}
        & 1
        \\
        1 &
    \end{matrix}\right),
    \qquad
    S_j^y
    =
    \left(\begin{matrix}
        & -i
        \\
        i &
    \end{matrix}\right),
    \qquad
    S_j^z
    =
    \left(\begin{matrix}
        1 & 
        \\
        & -1
    \end{matrix}\right)
    .
```
結合定数は仕事をしないので今は無視します。


## 実装

### エネルギースペクトルの描画

```julia
using LinearAlgebra
using Plots
using Arpack

const ⊗ = kron
```
はじめにHamiltonianを用意。
テンソル積は安直にKronecker積で実装可能です。
VS Code では \otimes [tab] で ⊗ 記号が出てきました。
```julia
function pauli()
    Sx = [0 1; 1 0]
    Sy = [0 -im; im 0]
    Sz = [1 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, 2^L,2^L)
    for i in 1:(L-1)
        H += I(2^(i-1)) ⊗ heisenbergint(Sx,Sy,Sz) ⊗ I(2^(L-i-1))
    end
    return H
end
```

`Arpack.eigs`を使って実行してみます。
サイト数 `L=9` まではストレスなく計算できますが、`L=10` にて30秒近くかかりました。
厳密対角化の限界がわかります。
```julia
L = 10
@time eigenvalues = eigs(hamiltonian(L), nev=2^L, which=:SR)[1]
@show eigenvalues
```
また、実行してみると実際に出てくる固有値の数が想定よりも小さいようです。
`hamiltonian(L)` の次元はしっかり $2^L$ で正しいので、`Arpack.eigs` の仕様と齟齬があるのだと思います。
取り出す固有値の数はアジャストしてくれるので、(少なくとも見た目は) 結果に影響ないようです。

描画してみます。
`L=10` に注意してください。
```julia
plt = plot(yshowaxis=false)
xlabel!("energy")
title!("Spin-1/2 AF Heisenberg chain energy spectrum")

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

display(plt)
```
![spin-halfHeisenberg-spectrum.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3875030/6ecb1ab5-6123-c3f7-fd50-0356f4f16f62.png)

2024/09/11追記:
コメントの指摘に基づき、forループの中身にて`plot!(fill(abs(...`を`plot!(fill(real(...`へ変更しました。

### ギャップの描画

先の画像からは、あまりgaplessさが読み取れません。
Spin-1/2 反強磁性 Heisenberg モデルの励起エネルギーが $\mathcal{O}(L^{-1})$ であることは Lieb-Schultz-Mattis の定理 [[1]](#reference) として物性業界では比較的よく知られています。
確認のためサイト数 `L` と第1励起エネルギーをプロットしてみます。
```julia
gaps = []
range = 3:14
for L in range
    @time eigvals = eigs(hamiltonian(L), nev=2, which=:SR)[1]
    gaps = append!(gaps, eigvals[2] - eigvals[1])
end

plt = plot(range, real(gaps), ylim=(0,2.7), label="")
xlabel!("sites")
ylabel!("E1 - E0")

display(plt)
```
![Shalf_gap.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3875030/1ed363b2-6767-a630-c326-3814c02e88ab.png)
手元の計算環境では14サイトで150秒近くかかったので、これ以上の探索は断念します。

奇数サイトにてギャップが 0 (基底状態が縮退) になっていますが、[[2]](#reference) の結果から、基底状態は合計スピンが $S=1/2$ となって厳密に二重縮退することがわかります。
一方、偶数サイトであれば基底状態は unique になることが [[2], [3]](#reference) の結果から知られています。
この辺りの話は田崎さんの本 [[4]](#reference) が非常にわかりやすいです。




## Reference

1. E. Lieb, T. Schultz, D. Mattis, _Two soluble models of an antiferromagnetic chain_. Annals of Physics **16**, 407 (1961)
2. E. Lieb, D. Mattis, *Ordering energy levels in interacting spin chains*. J. Math. Phys. **3**, 749 (1962)
3. W. Marshall, *Antiferromagnetism*. Proc. R. Soc. A **232**, 48 (1955)
4. H. Tasaki *Physics and Mathematics of Quantum Many-Body Systems*. Springer Cham, 2020. ISBN: 978-3-030-41264-7.