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

モデル

    H=\sum_j{S}_j\cdot{S}_{j+1}

where

    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)
    .

結合定数は仕事をしないので今は無視します。

実装

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

using LinearAlgebra
using Plots
using Arpack

const ⊗ = kron

はじめにHamiltonianを用意。
テンソル積は安直にKronecker積で実装可能です。
VS Code では \otimes [tab] で ⊗ 記号が出てきました。

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秒近くかかりました。
厳密対角化の限界がわかります。

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

また、実行してみると実際に出てくる固有値の数が想定よりも小さいようです。
hamiltonian(L) の次元はしっかり $2^L$ で正しいので、Arpack.eigs の仕様と齟齬があるのだと思います。
取り出す固有値の数はアジャストしてくれるので、(少なくとも見た目は) 結果に影響ないようです。

描画してみます。
L=10 に注意してください。

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)

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

ギャップの描画

先の画像からは、あまりgaplessさが読み取れません。
Spin-1/2 反強磁性 Heisenberg モデルの励起エネルギーが $\mathcal{O}(L^{-1})$ であることは Lieb-Schultz-Mattis の定理 [1] として物性業界では比較的よく知られています。
確認のためサイト数 L と第1励起エネルギーをプロットしてみます。

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)

手元の計算環境では14サイトで150秒近くかかったので、これ以上の探索は断念します。

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

Reference

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