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

この記事は？

緑Julia本「1週間で学べる！ Julia数値計算プログラミング」の6日目に入れるために書いたのですが、あまりにも分量が多くなってしまって泣く泣くカットした部分です。発売から2年経ちましたし、供養のために公開します。多いので、二つの記事に分けます。この記事はpart2です。part1はこちら。

多体電子系の厳密対角化

part1はこちら。

ハミルトニアンの行列表示

さて、Juliaで多体問題を取り扱う最低限の準備は整いました。ここでは詳細を述べませんが、生成消滅演算子を用いた場合、相互作用のある$N$サイトの1次元強束縛模型の一番単純な模型は

\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サイトの問題の場合は、模型は

-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$を持つ消滅演算子をまとめて

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です。
これで、

nsite = 2
Cs = make_Cs(nsite)
display(Cs[1,"up"])
println("")
display(Cs[1,"up"]')
println("")

とすると、消滅演算子Cs[1,"up"]や生成演算子Cs[1,"up"]'がちゃんと計算できていることがわかります。これを使って

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サイトの場合には

nsite = 1
H = make_hamiltonian(nsite,1,10)
μ = 1
U = 10
display(H)
println("")

のようにすると、

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

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$に依存してどのように変化するかを調べてみましょう。そのためには、電子数の期待値を計算する必要があります。詳細については述べませんが、「電子数演算子」は

{\cal N} = \sum_{i} \sum_{\sigma} c_{\sigma i}^{\dagger} c_{\sigma i}

で定義されます。ですので、絶対零度における電子数の期待値は

N = \langle \psi_0 | \sum_{i} \sum_{\sigma} C_{\sigma i}^{\dagger} C_{\sigma i} | \psi_0 \rangle

と書くことができます。ここで$|\psi_0 \rangle$は最小固有エネルギーをもつ固有状態です。生成消滅演算子は行列として表現できていますので、この期待値も行列の演算で簡単に計算ができます。
電子数の期待値の化学ポテンシャル依存性を計算してみましょう。まず、電子数の期待値を計算する関数を

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マクロが便利です。これを使うと、

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は反例を左上に置くと指示したものです。このコードを実行した結果、

という図が得られます。電子数が整数で階段状に変化していることからそれぞれの電子数の異なる状態は混ざっていないことがわかります。混ざっていれば中途半端な数になるはずですから。これは、ハミルトニアンが電子数演算子と可換であるためにブロック対角化されているからです。なお、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電子の状態に変化したからと言えます。
なお、固有値を求めるときに

e,v = eigen(Matrix(H)) 

を使っていますが、これはeigenが密行列用の関数のために疎行列の$H$をMatrix(H)とすることで密行列に変換しています。これはせっかくの疎行列の利点を消してしまっています。もし今のように最小エネルギー状態だけが欲しいのであれば、KrylovKitという疎行列用のパッケージを使うとよいです。これは連立方程式や固有値問題を解いてくれます。インストールはいつもと同じようにREPLで]を押してからadd KrylovKitとすればできます。使う場合はusing KrylovKitとします。これを使う場合には、コードを

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として定義します。

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}という型を

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

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)とほとんど変わっていませんね。
あとは、電子を消す関数

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

と状態からインデックスを取り出す関数

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)を多重ディスパッチで定義すれば：

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)は組み合わせの総数を返す関数で

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)を作り、関数の最初に

    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$の行列を対角化し、固有値を求めてみましょう。ハミルトニアンを作る関数はすでに定義済みですから、固有値を求めるには

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]に格納されています。この基底状態の固有ベクトルの中身を図示してみましょう。そのためには、

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()で呼び出すと、

という図になります。これをみると、$|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であれば、行列を定義しなくても、行列とベクトルの積さえ定義されていれば、固有値を求めることができます。ここで作ったコードを正解だと考えて、より高速なコードのデバッグをするのも良いでしょう。