多体電子系の厳密対角化 part2

# この記事は？
緑Julia本[「1週間で学べる！ Julia数値計算プログラミング」](https://www.kspub.co.jp/book/detail/5282823.html)の6日目に入れるために書いたのですが、あまりにも分量が多くなってしまって泣く泣くカットした部分です。発売から2年経ちましたし、供養のために公開します。多いので、二つの記事に分けます。この記事はpart2です。part1は[こちら](https://qiita.com/cometscome_phys/private/ee91276bff3876848b34)。

# 多体電子系の厳密対角化

part1は[こちら](https://qiita.com/cometscome_phys/private/ee91276bff3876848b34)。


### ハミルトニアンの行列表示
さて、Juliaで多体問題を取り扱う最低限の準備は整いました。ここでは詳細を述べませんが、生成消滅演算子を用いた場合、相互作用のある$N$サイトの1次元強束縛模型の一番単純な模型は
```math
\sum_{i} \sum_{\sigma } -t (c_{\sigma i+1}^{\dagger} c_{\sigma i}+c_{\sigma i-1}^{\dagger} c_{\sigma i}) +\sum_{i} \sum_{\sigma} E_{\sigma} c_{\sigma i}^{\dagger} c_{\sigma i}+ U \sum_i c_{\uparrow i}^{\dagger} c_{\uparrow i}c_{\downarrow i}^{\dagger} c_{\downarrow i}
```
という形で書かれます。この模型は1次元ハバード模型と呼ばれ、強く相互作用する電子集団を解析するために非常に盛んに研究されてきました。この模型の2次元版は銅酸化物高温超伝導体のミニマル模型としても有名です。
例えば2サイトの問題の場合は、模型は
```math
-t \sum_{\sigma } (c_{\sigma 2}^{\dagger} c_{\sigma 1}+c_{\sigma 1}^{\dagger} c_{\sigma 2}) +\sum_{i=1}^2 \sum_{\sigma} E_{\sigma} c_{\sigma i}^{\dagger} c_{\sigma i} + U \sum_{i=1}^2 c_{\uparrow i}^{\dagger} c_{\uparrow i}c_{\downarrow i}^{\dagger} c_{\downarrow i}
```
となります。
あとはJuliaでハミルトニアンの固有値固有ベクトルを計算してみましょう。まず、$i$番目のサイトのスピン$\sigma$を持つ消滅演算子をまとめて
```julia 
function make_Cs(nsite)
    Cs = Dict{Tuple,AbstractMatrix{Int64}}() #引数がタプル、値が行列の辞書型
    totalsite = 2*nsite
    for isite=1:totalsite
        spinindex = ifelse(isite <= nsite,"up","down")
        siteindex = ifelse(isite <= nsite,isite,isite -nsite)
        Cs[(siteindex,spinindex)] = make_C(isite,nsite)
    end
    return Cs
end
```
で定義してしまいます。ここで```Cs```は辞書型でして、```Cs[1,"up"]```とすると1番目のアップスピンの消滅演算子の行列表示が出てきます。```Dict{K,V}```もパラメトリック型でして、```K```が辞書のキーの型、```V```がそのキーに対応する辞書の値の型です。ここでは辞書のキーをタプルにし、値を行列としました。値の型が```Matrix```ではないのは、この後密行列```Matrix```ではない行列を扱う予定のためにより一般的な行列の型である```AbstractMatrix```を指定したからです。
次に、生成演算子も定義してしまいましょう。これはただ```'```をつければOKです。
これで、
```julia 
nsite = 2
Cs = make_Cs(nsite)
display(Cs[1,"up"])
println("")
display(Cs[1,"up"]')
println("")
```
とすると、消滅演算子```Cs[1,"up"]```や生成演算子```Cs[1,"up"]'```がちゃんと計算できていることがわかります。これを使って
```julia 
function make_hamiltonian(nsite,μ,U)
    Cs = make_Cs(nsite)
    H = zero(Cs[1,"up"]) #行列HをCs[1,"up"]と同じサイズのゼロ行列とする
    t = 1
    Eσ = -μ*ones(Float64,2)
    if nsite == 2
        for σ in ["up","down"]
            H += -t*Cs[2,σ]'*Cs[1,σ] - t*Cs[1,σ]'*Cs[2,σ]
        end    
    elseif nsite > 2
        for σ in ["up","down"]
            for isite=1:nsite
                jsite = isite+1
                jsite += ifelse(jsite > nsite,-nsite,0)

                H += -t*Cs[jsite,σ]'*Cs[isite,σ]

                jsite = isite-1
                jsite += ifelse(jsite < 1,nsite,0)

                H += -t*Cs[jsite,σ]'*Cs[isite,σ]
            end
        end 
    end
    for isite =1:nsite
        for (i,σ) in enumerate(["up","down"])
            H += Eσ[i]*Cs[isite,σ]'*Cs[isite,σ]
        end
        H += U*Cs[isite,"up"]'*Cs[isite,"up"]*Cs[isite,"down"]'*Cs[isite,"down"]
    end
    return H
end
```
のようにハミルトニアンを作る関数を作ります。ここで、$E_{\sigma} = - \mu$としました。$\mu$は化学ポテンシャルです。辞書型を使ったことで見た目が数式に近づいたと思います。
1サイトの場合には
```julia
nsite = 1
H = make_hamiltonian(nsite,1,10)
μ = 1
U = 10
display(H)
println("")
```
のようにすると、
```julia 
4×4 Matrix{Float64}:
 0.0   0.0   0.0  0.0
 0.0  -1.0   0.0  0.0
 0.0   0.0  -1.0  0.0
 0.0   0.0   0.0  8.0
```
が出てきますが、最初に示したようにちゃんと対角要素が$0,1=-\mu=,-1=-\mu=,8=U-2\mu$になっていることがわかります。

さて、行列として表示される生成演算子も消滅演算子もそれから作られるハミルトニアンも行列要素のほとんどがゼロになっていることがわかると思います。このような場合には0ではない要素だけを格納できる型「疎行列」を使った方が便利です。疎行列を使うには```SparseArrays```パッケージを読み込む必要があります。つまり、```using SparseArrays```をする必要があります。消滅演算子の行列を作る関数```make_C```の定義の手前で```using SparseArrays```を書いておき、```make_C```の中の```mat_C = zeros(Int8,ndim,ndim)```を```mat_C = spzeros(Int8,ndim,ndim)```に置き換えてください。これによって```H = make_hamiltonian(nsite,1,10)```の出力```display(H)```が
```julia 
4×4 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
  ⋅     ⋅     ⋅    ⋅ 
  ⋅   -1.0    ⋅    ⋅ 
  ⋅     ⋅   -1.0   ⋅ 
  ⋅     ⋅     ⋅   8.0
```
という形に変化します。```3 stored entries```とありますから非ゼロ要素は3しかありません。疎行列は非ゼロ要素だけを格納することメモリを大幅に削減することができます。また、行列の積など演算も非ゼロ要素同士だけの積で表すので高速に実行することができます。消滅演算子の集まりを定義した辞書型の値の型を```AbstractMatrix```型にしたおかげで、普通の行列でも疎行列でも問題なく動くコードになりました。

上で作った1サイトのハミルトニアンの固有状態は電子数が異なる状態の線型結合$|\psi \rangle = a_{00} |00\rangle + a_{01} |01\rangle + a_{10} |10\rangle + a_{11} |11 \rangle$からなると仮定しています。対角化して固有ベクトルを求めると、この係数$a_{ij}$が求まることになります。つまり、系の電子数は変化してもよいと考えてハミルトニアンを作っていることになります。そして、化学ポテンシャル$\mu$の値によって最小の固有エネルギーが変化します。これは統計力学におけるグランドカノニカル分布を考えていることになりますね。最小の固有エネルギーを持つ固有状態は基底状態ですが、その基底状態の電子数が化学ポテンシャル$\mu$に依存してどのように変化するかを調べてみましょう。そのためには、電子数の期待値を計算する必要があります。詳細については述べませんが、「電子数演算子」は
```math
{\cal N} = \sum_{i} \sum_{\sigma} c_{\sigma i}^{\dagger} c_{\sigma i}
```
で定義されます。ですので、絶対零度における電子数の期待値は
```math
N = \langle \psi_0 | \sum_{i} \sum_{\sigma} C_{\sigma i}^{\dagger} C_{\sigma i} | \psi_0 \rangle
```
と書くことができます。ここで$|\psi_0 \rangle$は最小固有エネルギーをもつ固有状態です。生成消滅演算子は行列として表現できていますので、この期待値も行列の演算で簡単に計算ができます。
電子数の期待値の化学ポテンシャル依存性を計算してみましょう。まず、電子数の期待値を計算する関数を
```julia 
function calc_N(Cs,state,nsite)
    N = 0
    for isite =1:nsite
        for (i,σ) in enumerate(["up","down"])
            Cstate =  Cs[isite,σ]*state
            N += Cstate'*Cstate
        end
    end
    return N
end
```
としました。ここで、```Cs```は上で作った消滅演算子の行列が入った辞書型、```state```は波動関数、```nsite```はサイト数です。この関数を使って、サイト数が1の時と2の時の「基底状態エネルギー」と「粒子数」の化学ポテンシャル依存性をまとめてプロットしてみましょう。複数のプロットをまとめるには```@layout```マクロが便利です。これを使うと、
```julia 
function testN()
    nsite = 1
    U = 10
    nsites = [1,2]
    p = [] #プロットを格納するための空の配列
    l = @layout [a b;c d] #行列のようにレイアウトを指定できる

    for nsite in nsites
        nmu = 200
        nstate = 4^nsite
        μs = range(-2U,2U,length=nmu)
        Cs = make_Cs(nsite)
        Ns = []
        es = []
        for (i,μ) in enumerate(μs)
            H = make_hamiltonian(nsite,μ,U)
            e,v = eigen(Matrix(H)) 
            N = calc_N(Cs,v[:,1],nsite)
            push!(Ns,N)
            push!(es,e[1])
        end

        p1 = plot(μs,es,xlabel="μ" ,label="$(nsite)-site ground state energy",legend = :bottomleft)
        push!(p,p1)
        p2 = plot(μs,Ns,xlabel="μ",label="$(nsite)-site electron number",legend = :topleft)
        push!(p,p2)
    end
    plot(p...,layout = l) #pは1次元配列。
    #ここには、p[1],p[2],p[3],p[4]と入れればよいが、p...とやるとタプルとして展開してくれる
    savefig("energy_num.png")

end
testN()
```
というコードができます。ここで```plot```のオプションの```legend = :topleft```は反例を左上に置くと指示したものです。このコードを実行した結果、

![energy_num.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246113/e0ab9de1-ad64-f063-ef04-4ec51a11a651.png)


という図が得られます。電子数が整数で階段状に変化していることからそれぞれの電子数の異なる状態は混ざっていないことがわかります。混ざっていれば中途半端な数になるはずですから。これは、ハミルトニアンが電子数演算子と可換であるためにブロック対角化されているからです。なお、1サイトの時にはハミルトニアンはすでに見たように$4 \times 4$でその固有値は$0$,$-\mu$,$-\mu$,$-2\mu+U$です。つまり、$\mu$を変化させていくと、$\mu < 0$の時はエネルギーはゼロ、$0 < \mu < U$の時エネルギーは$-\mu$、$\mu > U$の時エネルギーは$-2\mu+U$となります。グラフで$\mu > 10$となっているところで傾きが2倍になったのは基底状態となった波動関数が1電子の状態から2電子の状態に変化したからと言えます。
なお、固有値を求めるときに
```julia 
e,v = eigen(Matrix(H)) 
```
を使っていますが、これは```eigen```が密行列用の関数のために疎行列の$H$を```Matrix(H)```とすることで密行列に変換しています。これはせっかくの疎行列の利点を消してしまっています。もし今のように最小エネルギー状態だけが欲しいのであれば、KrylovKitという疎行列用のパッケージを使うとよいです。これは連立方程式や固有値問題を解いてくれます。インストールはいつもと同じようにREPLで]を押してから```add KrylovKit```とすればできます。使う場合は```using KrylovKit```とします。これを使う場合には、コードを
```julia 
e, v, info = eigsolve(H,numeigenvalues,:SR) #:SRはsmallest realを意味する
N = calc_N(Cs,v[1],nsite)
```

と変更すればOKです。密行列の場合は固有ベクトルは2次元配列として格納されていましたが、KrylovKitの場合は1次元配列の1次元配列(今回は求める固有値の数が1なので要素数1の1次元配列の要素が固有ベクトル)となっているので気をつけてください。サイト数が1や2では疎行列の恩恵を感じることができませんが、大きくしていくとわかります。


### 電子数が保存された場合のハミルトニアンの行列表示:Abstract typeの利用
上で考えていたハバード模型のハミルトニアンは、電子数を計算してわかりましたように、電子数ごとにブロック対角化されています。ブロック対角化されているのであれば、全体の対角化をするのではなく自分が考えたいブロックだけを対角化すればより高速に計算ができます。実際の研究では可能な限り計算コストを下げることで大きな計算を行いますので、このブロックだけ取り出して計算することは必須です。
これまでは電子数が異なる状態も基底として用意して生成消滅演算子を作ってきました。せっかく色々コードを作りましたから、それらのコードをなるべく再利用したいですよね。特に、独自に定義した```Electronbasis{N}```型はうまいことやれば電子数が同じ状態だけを取り出すように変更できそうです。そのためには、```Abstract type```を用います。
まず、これまで使ってきた```Electronbasis{N}```型を```Electronbasis_full{N}```型と名前を変えてしまいましょう。名前を変えてしまうと、これまで定義してきた引数を```Electronbasis{N}```とする関数が使えなくなってしまいます。これを使えるように、以下のように```Electronbasis{N}```型を```Abstract type```として定義します。
```julia 
abstract type Electronbasis{N} end

mutable struct Electronbasis_full{N} <: Electronbasis{N}
    state::Array{Bool,1}
    index::Int64

    function Electronbasis_full(N,index)
        state = i2state(index,N)
        return new{N}(state,index)
    end
end

function Electronbasis(N,index;fixed=false,numelectron=div(N,2))
    if fixed
        return Electronbasis_fix(N,index,numelectron)
    else
        return Electronbasis_full(N,index)
    end
end
```
ここで、```Abstract type```である```Electronbasis{N}```がこれまでと同様に呼び出された時の挙動を記述する関数を```Electronbasis(N,index;fixed=false,numelectron=div(N,2))```として定義しています。ここで```fixed=false```というのはパラメータ引数ですので、指定しない場合、例えば前のコードと同じように```Electronbasis(N,index)```とすると、```fixed```には```false```が入り、```Electronbasis_full(N,index)```が呼ばれることになるわけです。そして、```Electronbasis_full{N} <: Electronbasis{N}```の```A <: B```という形は、Aの親の型はBである、ということを意味しています。この指定によって、```Electronbasis{N}```が引数となっている関数を```Electronbasis_full{N}```でも同様に使うことができます。例えば、上で定義した```Base.display(electron::Electronbasis{N}) ```は関数```display(A)```の引数```A```が```Electronbasis{N}```型の時に呼ばれる関数でしたが、```A```が```Electronbasis_full{N}```の時もこの関数が呼ばれることになります。
このように```Electronbasis{N}```と```Electronbasis_full{N}```を定義したことによって、```Electronbasis{N}```型を引数とする関数は```Electronbasis_full{N}```を引数とする関数としても動くようになりました。ですので、この変更を行っても、上で書いたコードの他の箇所を修正する必要はありません。そのまま動くことを確認できると思います。
わざわざこのように```Abstract type```を導入したのは、「電子数が保存する場合の型」である```Electronbasis_fix```を定義するためです。さて、```Electronbasis_fix{N,numelectron}```という型を
```julia 
using Combinatorics
mutable struct Electronbasis_fix{N,numelectron} <: Electronbasis{N}
    state::Array{Bool,1}
    index::Int64

    function Electronbasis_fix(N,index,numelectron)
        allstates = collect(combinations(1:N,numelectron))
        state = i2state_fixed(index,N,numelectron,allstates)
        return new{N,numelectron}(state,index)
    end
end
```
と定義してみます。この形のフィールドの```state```は以前と同じように状態の情報が入っています。```index```も同じように状態を示すインデックスが入ります(前のコードと揃えるために0からスタートとします)。しかし、このインデックスは電子数が保存されている状態だけを数えてつけられるインデックスです。
電子数が保存されている状態について説明します。状態の候補の数(1サイトにアップスピンとダウンスピンがあるので2$\times$サイト数)が```N```で電子数が```numelectron```の時、その状態の総数は、```N```個の箱に```numelectron```個の同じ色のボールを入れる組み合わせの数と同じでして、$_{N}C_{\rm numelectron}$と書けます。Juliaでは順列や組み合わせを計算するパッケージCombinatoricsがありますので、いつもと同じように```add Combinatorics```でパッケージをインストールしてください。今考えたい組み合わせを全て入手するには、```collect(combinations(1:N,numelectron))```とします。```collect(A)```という関数は```A```を配列の形にするものと考えておいてください（正確には、```for ai in A```とした時の```ai```が要素として入る配列）。あとは、この組み合わせの配列の```index+1```番目の状態を作ります。その関数```i2state_fixed(index,N,numelectron,allstates)```は
```julia 
function i2state_fixed(i,N,numelectron,allstates)
    numstates = length(allstates)
    @assert 0 <= i <  numstates "index should be in  0 <= i < $(numstates)!"
    state = zeros(Bool,N)
    for j in allstates[i+1]
        state[j] = true
    end
    return state
end
```
と定義されています。最初に定義した電子数が保存していない時の関数```i2state(i,N)```とほとんど変わっていませんね。
あとは、電子を消す関数
```julia 
function erace!(electron::Electronbasis_fix{N,numelectron},isite,allstates_after) where {N,numelectron}
    @assert check_site(electron,isite) "there is no electron at $(isite)-site!"
    electron.state[isite] = false
    electron.index = state2i_fixed(electron.state,allstates_after)
end
```
と状態からインデックスを取り出す関数
```julia 
function state2i_fixed(state,allstates)
    data = []
    for (i,istrue) in enumerate(state)
        if istrue
            push!(data,i)
        end
    end
    ii = findfirst(x -> x == data,allstates)
    if ii == 0
        error("state $state is not in possible states")
    end
    return ii-1
end
```
の二つを定義さえすれば、消滅演算子を作る関数```make_C```も簡単に書くことができます。上で定義した```make_C(isite,nsite)```に電子数の情報を引数として入れた```make_C(isite,nsite,numelectron)```を多重ディスパッチで定義すれば：
```julia 
function make_C(isite,nsite,numelectron)
    N = 2*nsite
    ndim_before = calc_numstates(N,numelectron)
    ndim_after = calc_numstates(N,numelectron-1)
    allstates_after = collect(combinations(1:N,numelectron-1))

    mat_C = spzeros(Int8,ndim_after ,ndim_before)
    for index=0:ndim_before-1
        electron = Electronbasis(N,index,fixed=true,numelectron=numelectron)
        if check_site(electron,isite) #isiteに電子がいるかどうかチェック
            sign = check_sign(electron,isite) #消滅演算子の符号の計算
            erace!(electron,isite,allstates_after)
            jindex = electron.index
            mat_C[jindex+1,index+1] = sign #行列要素を決める
        end
    end
    return mat_C
end
```
ほとんどコードは完成したようなものです。ここで、電子数が異なる場合は、消滅演算子の行列は正方形ではなく長方形になることに注意してください。つまり、```numelectron```個の電子がある状態に消滅演算子を作用させると電子が一つ消えますので、```numelectron-1```の状態になります。そこで、```ndim_before```と```ndim_after```を定義しています。なお、ここで使った```calc_numstates(N,numelectron)```は組み合わせの総数を返す関数で
```julia 
function calc_numstates(N,numelectron)
    return binomial(N, numelectron)
end
```
と定義しています。```make_C(isite,nsite)```と```make_C(isite,nsite,numelectron)```を見比べてもらえば、ほとんど書き換えていないことがわかると思います。最後に、消滅演算子の行列の組を計算する関数を
```julia {.line-numbers}
function make_Cs(nsite,numelectron)
    Cs = Dict{Tuple,AbstractMatrix{Int64}}() #引数がタプル、値が行列の辞書型
    totalsite = 2*nsite
    for isite=1:totalsite
        spinindex = ifelse(isite <= nsite,"up","down")
        siteindex = ifelse(isite <= nsite,isite,isite -nsite)
        Cs[(siteindex,spinindex)] = make_C(isite,nsite,numelectron)
    end
    return Cs
end
```
と定義しましょう。あとは、```make_hamiltonian(nsite,μ,U)```という関数をコピーして```make_hamiltonian(nsite,μ,U;fixed=false,numelectron=nsite)```を作り、関数の最初に
```julia 
    if fixed
        Cs = make_Cs(nsite,numelectron)
        H = zero(Cs[1,"up"]'*Cs[1,"up"])
    else
        Cs = make_Cs(nsite)
        H = zero(Cs[1,"up"])
    end
```
を追加してください。これで電子数が保存された場合のハミルトニアンの行列表示を作ることができます。

### 多体系の基底状態：Mott絶縁体
さて、上で作ったコードを用いて、基底状態を調べてみましょう。考える系は、6サイトが鎖状に並んだ系とします。また、端同士は繋がっていて、鎖は輪になっている（周期境界条件）としましょう。化学ポテンシャル$\mu$の値は$\mu=U/2$とします。$\mu$と$U$がこの関係を満たすとき、ハミルトニアンは電子正孔対称性を持ち、電子数とサイト数が等しい状態（ハーフフィリング）が基底状態となります。6サイトある系を考えていますので、それぞれにアップスピンとダウンスピンの状態があり、全部で12個の電子が入りうるサイトがあります。$\mu=U/2$の時は、ハミルトニアン行列の次元は、6個の電子がこの12箇所に入りうる状態の数$_{12}C_6=924$となります。この$924 \times 924$の行列を対角化し、固有値を求めてみましょう。ハミルトニアンを作る関数はすでに定義済みですから、固有値を求めるには
```julia 
function testfixed()
    nsite = 6
    numelectron = nsite
    U = 10
    μ = U/2
    println(calc_numstates(nsite*2,numelectron))
    @time H = make_hamiltonian(nsite,μ,U,fixed=true,numelectron=numelectron)
    @time e, v, info = eigsolve(H,10,:SR) #:SRはsmallest realを意味する
    println("energies = $(e[:])")
end
testfixed()
```
とすればよいです。ここでは、固有値の小さい順に10個求めるようにしています。計算を実行すると、```e```に10個の数字が格納されます。一番小さい値```e[1]```が基底状態のエネルギーです。対応する固有ベクトルは```v[1]```に格納されています。この基底状態の固有ベクトルの中身を図示してみましょう。そのためには、
```julia 

function plot_eigfunc(v,N,numelectron;eps=0.5e-2)
    p = []
    ndim = calc_numstates(N,numelectron)
    labels=[]
    vfinite = []
    vsum = 0.0
    vsum2 = 0.0
    for index=0:ndim-1
        electron = Electronbasis(N,index,fixed=true,numelectron=numelectron)
        label = get_label(electron)
        if abs(v[index+1])^2 > eps
            push!(labels,label)
            push!(vfinite,v[index+1])
            vsum2 += abs(v[index+1])^2
        end
        vsum += abs(v[index+1])^2
    end
    plot!(abs.(vfinite).^2,xticks=(1:1:ndim, labels),ylim=(0,0.5),xrotation = 45,label="ground state of $(N ÷ 2)-site Hubbard model",xtickfontsize=6,markershape = :circle, margin = 15mm)
    savefig("groundstate_$N.png")

end
```
という関数を定義します。この関数の内部で呼ばれている関数は全て上で定義済みの関数です。```plot```では、ラベルとして波動関数の2進数表記を使っています。この関数を```testfixed()```で呼び出すと、


![groundstate_12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246113/8db8bd9f-a047-04a7-db11-dd5de6948484.png)

という図になります。これをみると、$|101010010101 \rangle$のような状態の重みが大きいことがわかります。$|101010010101 \rangle$は、左の6つがダウンスピンの状態、右の6つがアップスピンの状態を表しています。ですので、$|101010010101 \rangle$という状態は、ダウンスピンが1,3,5番目のサイトにいてアップスピンが2,4,6番目のサイトにいる状態ということです。スピンの表示で書くと$|\downarrow \uparrow \downarrow \uparrow \downarrow \uparrow \rangle$です。この状態と$|010100101010 \rangle$という状態(スピンの表示で書くと$|\uparrow \downarrow \uparrow \downarrow \uparrow \downarrow \rangle$)が基底状態の重みとして一番大きいようです。スピンの表示からわかりますように、この基底状態は古典的な反強磁性状態の重ね合わせ状態が主成分です。また、他の成分はよく見るとアップスピンの数とダウンスピンの数が等しいものになっています。つまり、全体の磁化はゼロです。また、電子はUが大きいためにアップスピンとダウンスピンが同じ場所にいるとエネルギーを大きく損をします。そのため、基底状態は、各サイトに電子が一つずつ入る状態になっています。この基底状態は絶縁体です。なぜならば、電子を動かそうとすると必ずどこかで一つのサイトに二つの電子が入ることになり、それは非常にエネルギーを損するからです。このように、電子同士の相互作用によって電子がどこへも動けない状態を、Mott絶縁体状態と呼びます。
なお、ここで扱ったハーフフィリングのハバード模型は、Uが大きい極限でスピンが並んだハイゼンベルグ模型というスピンの模型になることが知られています。そして、そのハイゼンベルグ模型のスピン同士の相互作用は反強磁性的であることが摂動論によってわかっています。つまり、ハバード模型のUが大きい時の基底状態は、ハイゼンベルグ模型を解いて得ることもできます。ハイゼンベルグ模型の詳細についてはこの本では述べませんが、Juliaで今回と同じような形で簡単に実装することができます。

今回、生成消滅演算子を行列として表現するという比較的素朴な方法でハバード模型のハミルトニアンを構築しました。この方法でも8サイト位までは計算できますが、それ以上のサイト数ではハミルトニアン行列の作成に時間がかかり過ぎてしまいます。実際に研究でハバード模型を扱う場合には、ハミルトニアン行列を作るのではなく、ハミルトニアン行列とベクトルの演算を定義して扱うことが多いです。Juliaの行列対角化ができるパッケージのKrylovKitであれば、行列を定義しなくても、行列とベクトルの積さえ定義されていれば、固有値を求めることができます。ここで作ったコードを正解だと考えて、より高速なコードのデバッグをするのも良いでしょう。